一种基于指数分解约束的磁共振谱去噪方法

1.本发明涉及磁共振谱去噪方法，尤其是涉及一种基于指数分解约束的求解磁共振谱去噪问题方法。


背景技术：

2.众所周知，磁共振谱是生物、化学、医学等领域的重要分析工具之一。然而，磁共振谱的采集过程往往会受到噪声的影响，存在信噪比低的问题。在实际应用中，往往会采用重复采集取平均的方式来降低噪声，然而，多次重复采集会增加采样时间，造成时间的浪费。利用信号处理的方式进行磁共振谱去噪是从带有噪声的原始信号中提取出无噪的信号，通常以高斯噪声建模，从而提高核磁共振信号的信噪比。
3.在过去的几十年里，研究学者们提出很多磁共振谱的去噪方法。比如，cadzow(yung,ya lin,lian pin,hwang,"nmr signal enhancement based on matrix property mappings,"journal of magnetic resonance,series a,103,109-114,1993)是一种经典的去噪方法，这种方法通过规定截断矩阵的汉克尔矩阵结构进行去噪，但是得到的解不一定为全局最优解。在磁共振谱去噪中，许多研究人员利用磁共振信号的数学特性进行磁共振谱去噪重建，比如通过磁共振谱时域信号的汉克尔矩阵的低秩特性(xiaobo qu,maxim,mayzel,jian-feng cai,zhong chen,vladislav orekhov,"accelerated nmr spectroscopy with low-rank reconstruction,"angewandte chemie international edition,54,852-854,2015)。qu等人的方法通过构建基于汉克尔矩阵核范数的拉格朗日凸优化求解模型，对信号进行降噪，并应用到更高维谱重建(jiaxi ying,hengfa lu,qingtao wei,jian-feng cai,di guo,jihui wu,zhong chen,xiaobo qu,"hankel matrix nuclear norm regularized tensor completion for n-dimensional exponential signals,"ieee transactions on signal processing,65(14),3702-3717,2017)和时域频域信号混合信号重建(hengfa lu,xinlin zhang,tianyu qiu,jian yang,jiaxi ying,di guo,zhong chen,xiaobo qu,"low rank enhanced matrix recovery of hybrid time and frequency data in fast magnetic resonance spectroscopy,"ieee transactions on biomedical engineering,65(4),809-820,2018)。最近，ying等人提出了结合范德蒙分解的方法对原有算法进行优化(jiaxi ying,jian-feng cai,di guo,gongguo tang,zhong chen,xiaobo qu,"vandermonde factorization of hankel matrix for complex exponential signal recovery—application in fast nmr spectroscopy,"ieee transactions on signal processing,66(21),5520-5533,2018)，使低秩约束更为合理，但该方法在采样点数较多的情况下去噪时间较长。(tianyu qiu,wenjing liao,yihui huang,jinyu wu,di guo,dongbao liu,xin wang,jian-feng cai,bingwen hu,xiaobo qu,"an automatic denoising method for nmr spectroscopy based on low-rank hankel model,"ieee transactions on instrumentation and measurement,70,1-12,2021)提出一种自动选取参数的去噪方法，自主得到最优解。(xiaobo qu,yihui huang,
hengfa lu,tianyu qiu,di guo,vladislav orekhov,zhong chen,"accelerated nuclear magnetic resonance spectroscopy with deep learning,"angewandte chemie international edition,59(26),10297-10300,2020)，使用指数模型构建仿真信号来训练神经网络，可以快速地得到高保真的磁共振谱重建结果。(yihui huang#,jinkui zhao#,zi wang,di guo,xiaobo qu,"exponential signal recovery with deep hankel matrix factorization,"ieee transactions on neural networks and learning systems,1-13,2021)提出一种采用低秩汉克尔矩阵分解的深度学习方法，其神经网络结构是通过模仿指数信号重建方法的迭代过程构建的。该方法可以获得更低的重建误差，并且能够有效的保留弱峰。
4.在磁共振谱去噪中，现有方法的去噪质量仍有待提高，并且基于汉克尔矩阵奇异值的性质并未得到充分应用。


技术实现要素：

5.本发明的目的在于提供一种可同时利用磁共振谱时域信号的汉克尔矩阵低秩特性以及指数分解约束，以实现磁共振谱高质量去噪的一种基于指数分解约束的磁共振谱去噪方法。
6.本发明包括以下步骤：
7.1)提出一个同时利用磁共振谱时域信号的汉克尔矩阵低秩特性和指数分解约束的磁共振谱去噪模型；
8.2)通过迭代算法对磁共振谱进行去噪。
9.在步骤1)中，所述提出一个同时利用磁共振谱时域信号的汉克尔矩阵低秩特性和指数分解约束的磁共振谱去噪模型的具体方法为：
10.用表示不含噪声的磁共时域信号；
11.磁共振的自由衰减信号(free induction decay，fid)可以表示为
12.其中，n表示信号的长度；可以建模为衰减指数信号的叠加，将xn表示为：
[0013][0014]
其中，n＝0,
…
,n-1.cr、fr和τr分别为第r个峰的幅度、共振频率以及衰减系数，r＝1,2,
…
,r，其中r表示信号谱峰的总个数，δt表示采样间隔。
[0015]
定义算子p＝n-q 1，表示将fid信号转换为汉克尔矩阵的过程，汉克尔矩阵可以表示为：
[0016][0017]
当n为奇数时，q＝(n 1)/2；当n为偶数时，q＝n/2。
[0018]
该fid信号的汉克尔矩阵的秩与指数信号个数(也即谱峰个数)相同；当指数个数r
小于汉克尔矩阵的维度时，具有低秩特性；
[0019]
设所有谱峰的索引集为|
·
|表示索引集中的元素个数。可以将公式(1)中的指数的幅度进一步表示为并将cr排列为向量信号中的所有zr可以排列成范德蒙矩阵其中将表示为
[0020]
设为待去噪信号，由x转化成的汉克尔矩阵满足低秩特性，同时拥有x＝zc的分解形式，即具有指数分解约束。
[0021]
提出的去噪模型为：
[0022][0023]
其中，为带噪声的测量信号，为待去噪信号,是一个范德蒙矩阵。‖
·
‖
*
是核范数，是低秩约束的凸松弛，‖
·
‖2表示向量的2范数，η表示数据保真的范围。
[0024]
这个模型可以进一步被表示为：
[0025][0026]
其中，λ和μ为正则化参数，都是正实数。
[0027]
步骤2)中，所述通过迭代算法对磁共振谱进行去噪的具体方法可为：
[0028]
公式(4)可以通过交替乘子法进行计算；首先，式(4)的增广拉格朗日形式为:
[0029][0030]
其中d为对偶变量，《
·
,
·
》表示内积，表示矩阵f范数的平方，λ＞0,β＞0以及μ＞0都是正则化参数。
[0031]
模型(5)的子问题可以表述为：
[0032][0033]
公式(6)的各个子问题通过以下方法进行求解：
[0034]
(1)x
i 1
可以通过求解下式获得：
[0035][0036]
对上式求导，可以的到其解为
[0037][0038]
其中，为的伴随算子。
[0039]
(2)z
i 1
可以通过求解下式获得：
[0040][0041]
因此，可以得到z
i 1
的求解步骤为：
[0042]
第一步：计算汉克尔矩阵的奇异值分解其中(
·
)h表示矩阵的共轭转置然后对分解矩阵u∑vh截断为其中下标r表示取原矩阵的前r列；
[0043]
第二步：计算u
rl
mr≈u
rf
的最小二乘解mr，上标l和f分别表示删除矩阵ur的第一行和删除矩阵ur的最后一行。计算mr的特征值，并将其视作r＝1,2,
…
,r。
[0044]
第三步：将r＝1,2,
…
,r排列成范德蒙矩阵
[0045]
(1)c
i 1
可以通过对下式求导获得
[0046][0047]
其解为
[0048]ci 1
＝((z
i 1
)hz
i 1
)-1
(z
i 1
)hx
i 1
ꢀꢀꢀ
(11)
[0049]
其中，(
·
)-1
表示求伪逆。
[0050]
(2)u
i 1
可以通过求解下式获得：
[0051][0052]
得到解为：
[0053][0054]
(3)d
i 1
的解为：
[0055][0056]
其中t表示更新的步长。
[0057]
整个迭代过程为：
[0058]
初始化：输入y，以及谱峰总数r，设置β＝t＝1，μ＝0.1，数据保真的范围η＝10-8
。对偶变量d0＝1，变量x0＝y,z0＝0,c0＝0,u0＝0，δx＝10-10
，其中1和0分别表示与相应矩阵维度相同的元素全为1和元素全为0的矩阵或向量。
[0059]
当δx＞η时，重复计算步骤a)-e)，并计算出相邻两次的去噪磁共振谱的fid的误差当相邻两次的去噪磁共振谱的fid的误差缩小到指定范围时，即δx＜η时，停止迭代运算，将最后一步迭代的去噪结果x作为最终的去噪结果。
[0060]
本方法设计磁共振谱的去噪方法，尤其同时利用汉克尔矩阵奇异值分解的平移不变性以及其低秩特性。所述方法包括：1)提出一个同时利用磁共振谱时域信号的汉克尔矩阵低秩特性以及指数分解约束的磁共振谱去噪模型；2)通过迭代运算对磁共振谱进行去噪。本发明用于降噪磁共振谱，具有降噪质量高的特点。
附图说明
[0061]
图1为本发明实施例中的无噪波谱。
[0062]
图2为本发明实施例中输入的带噪波谱。
[0063]
图3为本发明实施例中的带噪波谱经过降噪后的波谱。
具体实施方式
[0064]
为使本发明的技术方案更清楚，以下实施例将结合附图对本发明作进一步的说明。应当理解，优选实施例仅为了说明本发明，而不是为了限制本发明的保护范围。
[0065]
本发明实施例用指数函数生成不含噪声的磁共振时域信号，该信号长度为n＝512，谱峰总个数r＝5。对其进行傅里叶变换，即可得到不含噪声的波谱，如图1所示，并在该时域信号上加上了标准差为0.02的高斯白噪声，对其进行傅里叶变换，即可得到加上噪声的波谱，如图2所示。
[0066]
本发明实施例包括以下步骤：
[0067]
1)提出一个同时利用磁共振谱时域信号的汉克尔矩阵低秩特性和指数分解约束的磁共振谱去噪模型：
[0068]
用表示不含噪声的磁共时域信号。
[0069]
磁共振的自由衰减信号(free induction decay，fid)可以表示为其中n表示信号的长度。可以建模为衰减指数信号的叠加，将xn表示为
[0070][0071]
其中，n＝0,
…
,n-1.cr、fr和τr分别为第r个峰的幅度、共振频率以及衰减系数，r＝1,2,
…
,r，其中r表示信号谱峰的总个数，δt表示采样间隔。
[0072]
定义算子p＝n-q 1，表示将fid信号转换为汉克尔矩阵的过程，汉克尔矩阵可以表示为
[0073][0074]
当n为奇数时，q＝(n 1)/2；当n为偶数时，q＝n/2。
[0075]
该fid信号的汉克尔矩阵的秩与指数信号个数(也即谱峰个数)相同。当指数个数r小于汉克尔矩阵的维度时，具有低秩特性。
[0076]
设所有谱峰的索引集为|
·
|表示索引集中的元素个数。可以将公式(1)中的指数的幅度进一步表示为并将cr排列为向量信号中
的所有zr可以排列成范德蒙矩阵其中将表示为
[0077]
设为待去噪信号，由x转化成的汉克尔矩阵满足低秩特性，同时拥有x＝zc的分解形式，即具有指数分解约束。
[0078]
提出的去噪模型为：
[0079][0080]
其中，为带噪声的测量信号，为待去噪信号,是一个范德蒙矩阵。‖
·
‖
*
是核范数，是低秩约束的凸松弛，‖
·
‖2表示向量的2范数，η表示数据保真的范围。
[0081]
这个模型可以进一步被表示为：
[0082][0083]
其中，λ和μ为正则化参数，都是正实数。
[0084]
2)通过迭代算法对磁共振谱进行去噪：
[0085]
模型(4)可以通过交替乘子法进行计算。首先，(4)式的增广拉格朗日形式为:
[0086][0087]
其中d为对偶变量，《
·
,
·
》表示内积，表示矩阵f范数的平方，λ＞0,β＞0以及＞0都是正则化参数。
[0088]
模型(5)的子问题可以表述为：
[0089][0090]
公式(6)的各个子问题通过以下方法进行求解：
[0091]
a)x
i 1
可以通过求解下式获得
[0092][0093]
对上式求导，可以的到其解为
[0094][0095]
其中，为的伴随算子。
[0096]
b)z
i 1
可以通过求解下式获得
[0097]
[0098]
因此，可以得到z
i 1
的求解步骤为
[0099]
第一步：计算汉克尔矩阵的奇异值分解其中(
·
)h表示矩阵的共轭转置然后对分解矩阵u∑vh截断为其中下标r表示取原矩阵的前r列；
[0100]
第二步：计算u
rl
mr≈u
rf
的最小二乘解mr，上标l和f分别表示删除矩阵ur的第一行和删除矩阵ur的最后一行。计算mr的特征值，并将其视作r＝1,2,
…
,r。
[0101]
第三步：将r＝1,2,
…
,r排列成范德蒙矩阵
[0102]
c)c
i 1
可以通过对下式求导获得：
[0103][0104]
其解为：
[0105]ci 1
＝((z
i 1
)hz
i 1
)-1
(z
i 1
)hx
i 1
ꢀꢀꢀ
(11)
[0106]
其中，(
·
)-1
表示求伪逆。
[0107]
d)u
i 1
可以通过求解下式获得：
[0108][0109]
得到解为：
[0110][0111]
e)d
i 1
的解为
[0112][0113]
其中t表示更新的步长。
[0114]
整个迭代过程为：
[0115]
初始化：输入y，以及谱峰总数r，设置β＝t＝1，μ＝0.1，数据保真的范围η＝10-8
。对偶变量d0＝1，变量x0＝y,z0＝0,c0＝0,u0＝0，δx＝10-10
，其中1和0分别表示与相应矩阵维度相同的元素全为1和元素全为0的矩阵或向量。
[0116]
当δx＞η时，重复计算步骤a)-e)，并计算出相邻两次的去噪磁共振谱的fid的误差当相邻两次的去噪磁共振谱的fid的误差缩小到指定范围时，即δx＜η时，停止迭代运算，将最后一步迭代的去噪结果x作为最终的去噪结果。
[0117]
对降噪后的时域信号x
out
进行傅里叶变换，即可得到降噪后的波谱，如图3所示。