本发明涉及磁共振波谱，尤其是涉及一种磁共振波谱的非理想成像因子的估计方法。
背景技术：
1、在实际的磁共振应用中，磁共振波谱通过人体代谢物的波谱来计算人体的代谢物含量，但这些波谱容易受到磁场漂移和不均匀这些非理想成像条件影响，导致代谢物波谱被引入非理想因子(包括整体的零阶相位偏移、每个代谢物的频率漂移、洛伦兹型谱峰的半高宽)，波谱出现失真，进而影响代谢物含量估计准确性。从给定波谱准确预测这些非理想因子，可以去除非理想成像条件对波谱影响，提升代谢物含量的估计准确性。
2、在量化磁共振波谱前，需要对一系列非理想因子做校正的预处理。常用的估计方法是基于均方误差(mean square error，mse)(s.w.provencher,"estimation ofmetabolite concentrations from localized in vivo proton nmr spectra,"magneticresonance in medicine,vol.30,pp.672-679,1993.)的估计。mse多数作为深度学习的损失函数和最优化算法的保真项，表达式如公式(1)。
3、
4、其中，cn为代谢物的浓度因子，mn(t)为代谢物的基集(basis-set)，i为复数的虚部，y为具有t采样采样点的mrs观测信号，为哈达玛乘法，为相位偏移，δωn为每个代谢物的频率漂移，σ为实际的洛伦兹型谱峰半高宽，表示预测值，exp(·)表示指数函数，||·||2表示二范数。
5、文献(y.yang,n.zou,e.lin,f.suo and z.chen,"a neural network method fornonconvex optimization and its application on parameter retrieval,"ieeetransactions on signal processing,vol.69,pp.3383-3398,2021.)研究证实公式(1)梯度的非凸性以及重叠峰或噪声对谱峰干扰导致的多个次优解。因此，公式(1)求解难度大，需额外设计更能多的求解策略或者约束项。在实际应用中，也存在通过峰形估计非理想因子又受限目标峰本身所处的环境影响(如谱峰重叠、噪声、背景信号干扰等)、深度学习估计非理想因子的训练成本高(标签数据难以获得)，泛化性难以验证等一系列问题。
技术实现思路
1、本发明的目的在于针对现有技术存在的上述技术问题，提供非梯度的、结合谱峰先验的启发式搜索的一种磁共振波谱的非理想成像因子的估计方法，本发明快速地、搜索出全局最优解，并且通过自身峰形本身估计初始可行解，保证搜索的解具有合理性与泛化性，减少非理想因子误差对短回波时间(echo time,te)的磁共振波谱的代谢物量化的影响，为提升代谢物量化准确性做好预处理。
2、本发明使用非梯度检索，具体包括以下步骤：
3、1)初始化非理想因子：利用非理想因子的先验信息，随机初始化频率漂移δωn和洛伦兹型谱峰半高宽σ，通过特定目标代谢物频段的最大值定位目标代谢物谱峰，对目标代谢物的峰所在频率估计得到初始化频率漂移δωn；
4、2)非理想因子搜索算法：根据步骤1)得到的初始化频率漂移δωn和洛伦兹型谱峰半高宽σ，搜索出相位设计一种用于评估给定初始频率偏移δωn和半高宽σ时的相位损失函数loss1，确保搜索的解具有合理性与泛化性；
5、3)信号求解算法：根据步骤2)得到的非理想因子参数对磁共振信号进行预估；设计一种迭代求解算法计算mse，设计利用已预估好相位信息的信号损失函数loss2，用于调整理想波谱去适应输入波谱，调整后输出校正谱；以减少非理想因子误差对短回波时间的磁共振波谱的代谢物量化的影响。
6、在步骤1)中，所述对目标代谢物的峰所在频率估计得到初始化频率漂移δωn，其中具有相似官能团如肌酸(cr)与磷酸肌酸(pcr)、磷酸胆碱(pch)与甘油磷酰胆碱(gpc)、n-乙酰天冬氨酸(naa)与n-乙酰天冬氨酰谷氨酸(naag)和谷氨酸(glu)与谷氨酰胺(gln)等代谢物给定相同的频率偏移δωn。
7、在步骤2)中，所述非理想因子搜索算法，具体步骤可为：
8、2.1)循环遍历每个个体i，初始化每个个体的速度vi([-5,5]的均匀分布随机数)以及位置φi([-180,180]的均匀分布随机数)，并初始化每个个体的历史最优位置pbesti为φi，循环结束后初始化种群最优个体最优位置gbest为min(pbesti)，其中min(·)表示最小值；种群个体数g是大于1的整数；
9、2.2)循环遍历每个个体i，根据个体速度更新位置φi；当|φi|>180时，其中|·|表示取绝对值，随机重置位置，重置方法为：
10、φi＝random([-180,180]) (2)
11、其中，random([*min,*max])为从*min到*max的均匀分布随机数；同时计算每个个体的适应度，计算公式为：
12、
13、其中，*real,*imag分别是取复数的实部或虚部，loss1(·)表示非理想因子预估的损失函数，m(·)为17种代谢物基集之和，具体公式如下：
14、
15、其中l是一个正整数；若根据公式(3)所得到的fiti>loss1(m(0),y(0),pbesti)，则更新个体历史最优位置，即将φi赋值给pbesti；若loss1(m(0),y(0),pbesti)>loss1(m(0),y(0),gbest)，则可以得到个体种群最优位置gbest为当前的pbesti；
16、2.3)循环遍历每个个体i，根据得到的pbesti和gbest更新速度，更新公式如下：
17、vi′＝c1×(pbesti-φi)+c2×(gbest-φi)+c3×vi (5)
18、其中，vi′为个体i更新后的速度，个体经验学习因子c1>0，社会经验学习因子c2>0，惯性因子c3>0，当|vi′|＞z时，随机重置速度，重置方法为：
19、vi′＝random([-z+1,z-1]) (6)
20、其中，z是大于0的常数，random([*min,*max])为从*min到*max的均匀分布随机数；
21、2.4)将步骤2.2)与步骤2.3)迭代循环j次。
22、在步骤3)中，所述信号求解算法，具体步骤可为：
23、3.1)循环遍历每个个体i，初始化每个个体的速度vi([vmin,vmax]的均匀分布随机数)以及位置pi([pmin,pmax]的均匀分布随机数)，其中pmax＝[δωn+5,σn+5]，pmin＝[δωn-5,σn-5]，vmax＝(pmax-pmin)/10,vmin＝(pmax-pmin)/20，并初始化初始化每个个体的历史最优位置pbesti为pi，循环结束后初始化种群最优个体最优位置gbest为min(pbesti)，其中min(·)表示最小值，种群个体数g是大于1的整数；
24、3.2)循环遍历每个个体i，根据个体速度更新位置pi；当pi＞pmax或pi＜pmin时，随机重置位置，重置方法为：
25、pi＝random([pmin,pmax]) (7)
26、其中，random([*min,*max])为从*min到*max的均匀分布随机数；计算每个个体的适应度，计算公式为：
27、
28、其中，m＝[m1(t),m2(t),m3(t),...,mn(t)]，y(t)为观测信号，为权利要求3预估的相位，loss2(·)为磁共振校正谱与输入谱之间的损失函数，loss2(·)的具体计算步骤为：
29、(a)对m进行调制得到m'，具体公式为：
30、
31、其中，(·,n)为取矩阵的第n列,为哈达玛乘法,exp(·)表示指数函数，j为虚部；
32、(b)利用线性最小二乘求解其中*h为矩阵的复共轭转置,*-1为矩阵的求逆操作；
33、(c)对m进行调制后求和得到y′(t)，具体公式为：
34、
35、其中cn为向量的n个值，(·,n)为取矩阵的第n列,为哈达玛乘法,exp(·)表示指数函数，j为虚部，y′(t)是通过损失函数loss2校正后的磁共振波谱；
36、(d)计算结构相似指标ssim，计算公式为：
37、
38、其中μ*为*的平均数,为*和·的协方差,为*的方差的平方,c1＝0.0001,c2＝0.0009；
39、若根据公式(8)所得到的则更新个体历史最优位置，即将pi赋值给pbesti；若则可得到个体种群最优位置gbest为当前的pbesti，并更新标志，令flag＝1；
40、3.3)若flag＝0，则更新停滞次数；停滞次数初始值为vmax，更新方法为当前停滞次数减一；若flag不为0，则在这一步中将flag设置为0；
41、3.4)循环遍历每个个体i，根据得到的pbesti和gbest更新速度，更新方法如下：
42、若在前1/3的迭代中，即[1,e/f]，更偏重个体惯性搜索，公式为：
43、v′i＝c1×(pbesti-pi)+(c2-0.7)×(gbest-pi)+(c3+0.7)×vi (12)
44、若在1/3到2/3的迭代中，即(e/f,2×e/f]，更偏重个体经验搜索，公式为：
45、v′i＝(c1+0.7)×(pbesti-pi)+(c2-0.7)×(gbest-pi)+c3×vi (13)
46、若在最后1/3的迭代中，即(2×e/f,3×e/f]，更偏重社会经验搜索，公式为：
47、v′i＝c1×(pbesti-pi)+c2×(gbest-pi)+c3×vi (14)
48、若vi>vmax或vi<vmin，则随机重置速度，公式为：
49、v′i＝random([vmin,vmax]) (15)
50、其中，个体经验学习因子c1>0,社会经验学习因子c2>0,惯性因子c3>0，迭代次数e是一个正整数，飞行习惯系数f>0，vmax＝(pmax-pmin)/10,vmin＝(pmax-pmin)/20，pmax＝[δωn+5,σn+5]，pmin＝[δωn-5,σn-5]，random([*min,*max])为从*min到*max的均匀分布随机数；
51、3.5)将步骤3.2)、步骤3.3)和步骤3.4)的迭代计算进行循环q次；
52、3.6)输出并且利用公式(10)得到经过目标波谱的校准谱。
53、本发明的优点和突出技术效果在于：
54、本发明提出一种磁共振波谱(magnetic resonance spectroscopy,mrs)的非理想成像因子的估计方法；本发明使用非梯度检索，首先利用三种常见的非理想因子(整体的零阶相位偏移每个代谢物的频率漂移δωn、实际的洛伦兹型谱峰半高宽σ)先验信息并进行非逐点遍历，高效地搜索出较好的可行解，然后设计一种包含改进型损失函数loss的磁共振波谱预估算法，结合已经估计出的非理想因子参数来准确估计波谱信号；该算法可用于快速准确地估计磁共振信号中的三种常见的非理想因子，旨在快速地、搜索出全局最优解，并且通过自身峰形本身估计初始可行解，保证搜索的解具有合理性与泛化性，减少非理想因子误差对短回波时间(echo time,te)的磁共振波谱的代谢物量化的影响，为提升代谢物量化准确性做好预处理。
1.一种磁共振波谱的非理想成像因子的估计方法，其特征在于使用非梯度检索，首先利用三种常见的非理想因子先验信息并进行非逐点遍历，高效搜索出较好的可行解，保证搜索的解具有合理性与泛化性；然后设计一种包含改进型损失函数loss的磁共振波谱预估算法，结合已经估计出的非理想因子参数准确估计波谱信号，减少非理想因子误差对短回波时间的磁共振波谱的代谢物量化的影响；所述三种常见的非理想因子为整体的零阶相位偏移每个代谢物的频率漂移δωn、实际的洛伦兹型谱峰半高宽σ。
2.如权利要求1所述一种磁共振波谱的非理想成像因子的估计方法，其特征在于利用非理想因子的先验信息，进行非逐点遍历，利用得到的非理想因子参数高效校正理想波谱，所述方法包括以下步骤：
3.如权利要求2所述一种磁共振波谱的非理想成像因子的估计方法，其特征在于在步骤1)中，所述对目标代谢物的峰所在频率估计得到初始化频率漂移δωn，其中具有相似官能团包括肌酸与磷酸肌酸、磷酸胆碱与甘油磷酰胆碱、n-乙酰天冬氨酸与n-乙酰天冬氨酰谷氨酸、谷氨酸与谷氨酰胺代谢物给定相同的频率偏移δωn。
4.如权利要求2所述一种磁共振波谱的非理想成像因子的估计方法，其特征在于在步骤2)中，所述非理想因子搜索算法，具体步骤为：
5.如权利要求4所述一种磁共振波谱的非理想成像因子的估计方法，其特征在于在所述非理想因子搜索算法写为如下的伪代码形式：
6.如权利要求2所述一种磁共振波谱的非理想成像因子的估计方法，其特征在于在步骤3)中，所述信号求解算法，具体步骤为：
7.如权利要求6所述一种磁共振波谱的非理想成像因子的估计方法，其特征在于所述信号求解算法写为如下的伪代码形式：
