结合稀疏和低秩特性的欠采样磁共振扩散谱的重建方法

技术领域

本发明涉及磁共振扩散谱的欠采样重建方法，尤其是涉及一种结合稀疏和低秩特性的欠采样磁共振扩散谱的重建方法。

背景技术

磁共振波谱是化学分析领域中一项十分重要的技术，通过采集的频谱可以测定分子的结构、分析人体组织的代谢产物等。若样品中包含多种化合物，再加上化合物本身J偶合的作用，往往导致一维磁共振谱图出现谱峰拥挤和重叠，这会增大谱峰归属难度，对分析和鉴定样品的化学成分造成不利影响。一个有效的方法是引入额外的谱维度，如引入不同类型的间接时间维可以得到二维J谱(分辨化学位移和自旋偶合)、二维相关谱(表明谱峰之间的相关性)等。引入扩散维则可利用不同物质具有不同的扩散系数从而有效区分不同的化合物。若同时引入间接时间维和扩散维，则可以得到三维扩散谱，其既可以通过扩散系数区分不同的化合物，还可得到化合物本身的结构。三维扩散谱使谱图的拥挤大大减少，谱峰归属更加容易，在化学分析领域有广泛的应用(C.S.Johnson Jr,"Diffusion orderednuclear magnetic resonance spectroscopy: principles and applications,"Progress in Nuclear Magnetic Resonance Spectroscopy,1999，34： 203-256)。

但高维磁共振实验花费时间较长，从几小时到几十天不等，使得实验必须耗费大量的仪器机时，常常令仪器不能满足科研人员的实验需求，采样时间长也增加了不稳定样品的数据采集难度，这些缺点都限制了高维磁共振技术在实际中的应用。将磁共振谱的时间信号建模成指数，利用指数信号构建的汉克尔矩阵具有低秩特性可以大大缩短磁共振实验时间(Xiaobo Qu,Maxim Mayzel,Jian-Feng Cai,Zhong Chen,Vladislav Orekhov,"Accelerated NMR spectroscopy with low-rank reconstruction,"Angewandte ChemieInternational Edition,2015，54(3)：852-854)。各种类型的磁共振谱均可利用低秩汉克尔矩阵方法进行欠采样和信号重建，比如混合时间-频率谱(Hengfa Lu,Xinlin Zhang,Tianyu Qiu,Jian Yang,Jiaxi Ying,Di Guo,Zhong Chen,Xiaobo Qu,VladislavOrekhov,"Low rank enhanced matrix recovery of hybrid time and frequency datain fast magnetic resonance spectroscopy,"IEEE Transactions on BiomedicalEngineering,2018，65(4)： 809-820)和高维张量谱(Jiaxi Ying,Hengfa Lu,QingtaoWei,Jian-Feng Cai,Di Guo,Jihui Wu, Zhong Chen,Xiaobo Qu,"Hankel matrixnuclear norm regularized tensor completion for N-dimensional exponentialsignals,"IEEE Transactions on Signal Processing,2017，65(14)： 3702-3717)等。低秩汉克尔矩阵方法又衍生出了多种改进算法，比如快速汉克尔矩阵重建算法(Di Guo,Hengfa Lu,Xiaobo Qu,"A fast low rank Hankel matrix factorizationreconstruction method for non-uniformly sampled magnetic resonancespectroscopy,"IEEE Access,2017，5：16033-16039)和低强度谱峰重建算法(Di Guo,Xiaobo Qu,"Improved reconstruction of low intensity magnetic resonancespectroscopy with weighted low rank Hankel matrix completion, "IEEE Access,2018，6：4933-4940)等。

在三维磁共振扩散谱中，可以在扩散维和间接时间维进行联合欠采样。由于直接时间维为全采样，当确定直接时间维上一点，则可以得到扩散维和间接时间维联合欠采样平面。可以通过约束此平面对应扩散谱的稀疏性进行重建(Mateusz Urbanczyk,WiktorKozminski,Krzysztof Kazimierczuk,"Accelerating diffusion-ordered NMRspectroscopy by joint sparse sampling of diffusion and time dimensions,"Angewandte Chemie International Edition,2014，126(25)： 6582-6585)。但该方法仅仅利用了扩散谱的稀疏性，未能利用扩散谱的物理模型和特性，易受噪声影响，在重建中产生孤立的噪声点，导致重建谱出现误差。

发明内容

本发明的目的在于提供结合稀疏和低秩特性的欠采样磁共振扩散谱的重建方法。

本发明包括以下步骤：

1)生成拉普拉斯-傅里叶联合变换矩阵；

在步骤1)中，所述生成拉普拉斯-傅里叶联合变换矩阵的具体方法可为：利用已知的扩散谱参数得到拉普拉斯变换矩阵L和傅里叶逆变换矩阵F，拉普拉斯变换矩阵可表示为：

其中，b为常数，与具体实验设置有关，g_m为扩散梯度，D_l为扩散谱的离散扩散系数；

傅里叶逆变换矩阵可表示为：

其中，t_n为扩散实验的离散采样时间，f_n为扩散谱的离散频率；则拉普拉斯-傅里叶联合变换矩阵其中为克罗内克积。

2)建立一种结合稀疏和低秩特性的欠采样重建模型：

其中，为欠采样算子，x为得到的时间维和扩散维的联合欠采样数据，s为待恢复的扩散谱向量，为将向量转为矩阵的算子，表示向量的二范数的平方，||g||₁表示向量的一范数，>*表示矩阵的核范数，λ₁与λ₂是平衡||s||₁和三项重要性的正则化参数；

公式(3)中的算子将向量转为矩阵的具体方法为：

设有向量a＝[a₁>2 L>L×N]^T，表示算子作用在向量a上，作用结果如下：

其中，a_i表示向量a中的第i个元素，的行数为L，列数为N；

3)基于结合稀疏和低秩特性的欠采样重建模型的求解算法；

在步骤3)中，所述基于结合稀疏和低秩特性的欠采样重建模型的求解算法的具体方法可为：利用交替方向乘子法求解(3)中的重建模型，引入变量z,Z和拉格朗日乘子C₁,C₂根据以下公式迭代更新变量：

当达到最大迭代次数K或s在相邻两次迭代中的误差小于设定的正数阈值μ时，迭代结束，其中，s^k+1,z^k+1,Z^k+1,和分别表示变量s，z，Z，D₁和D₂在第k+1次迭代时的值，上标H表示矩阵的复共轭转置，为将矩阵转换为向量的算子，表示奇异值收缩算子，阈值是参数β₁，β₂，τ₁，τ₂，λ₁和λ₂都是正数。

4)由步骤3)得到恢复的扩散谱向量s，经公式(4)中算子作用得到的为最终恢复的扩散谱。

本发明提供了一种可以从欠采样的磁共振扩散谱中重建出高精度和高信噪比的磁共振扩散谱的方法。首先根据实验参数生成拉普拉斯-傅里叶联合变换矩阵；接着建立一种基于结合稀疏和低秩特性的欠采样重建模型；然后通过迭代算法重建扩散谱向量；最后将扩散谱向量转为扩散谱。本发明实现利用少量数据重建完整的磁共振扩散谱，重建精度高，抗噪声能力强。

附图说明

图1是采样模板(白色表示采样到的点，黑色表示未采样到的点)。

图2是本发明重建的扩散谱。

图3是无噪全采样的扩散谱。

具体实施方式

下面通过具体实施例对本发明作进一步的说明，并给出重建结果。本实施例是一个重建时间维和扩散维联合欠采样的扩散谱的模拟实验。全采样的时间维与扩散维的大小分别是N＝128 和M＝64，采样模板(如图1所示)采样20％的数据，噪声标准差为时间信号最大值的0.3％，则本实施例中的核磁共振扩散谱数据点为8192点，采样20％时得到的总采样数据点数为1638 点。具体步骤如下：

1)生成拉普拉斯-傅里叶联合变换矩阵：利用已知的扩散谱参数得到拉普拉斯变换矩阵L 和傅里叶逆变换矩阵F。

拉普拉斯变换矩阵可表示为：

其中，b为常数，g_m为扩散梯度。D_l为扩散谱的离散扩散系数。本实验至为0至3.5×10¹⁰阵列64个点，D_l为5×10^-12至5×10^-9对数等分96个点。

傅里叶逆变换矩阵可表示为：

其中，t_n为扩散实验的离散采样时间，f_n为扩散谱的离散频率，本发明中t_n为0至0.0635s等分128个点，f_n为0Hz至2000Hz等分128个点；则拉普拉斯傅里叶联合变换矩阵其中为克罗内克积。

2)建立一种结合稀疏和低秩特性的欠采样重建模型：

其中，为欠采样算子，x为得到的时间维和扩散维的联合欠采样数据，s为待恢复的扩散谱向量，为将向量转为矩阵的算子，表示向量的二范数的平方，||g||₁表示向量的一范数，||g||_*表示矩阵的核范数。λ₁与λ₂是平衡||s||₁和三项重要性的正则化参数,λ₁与λ₂两者取值均为10^-3。

公式(3)中的算子将向量转为矩阵的具体方法为：

设有向量a＝[a₁>2 L>L×N]^T，表示算子作用在向量a上，作用结果如下：

其中，a_i表示向量a中的第i个元素，的行数为L＝96，列数为N＝128。

3)基于结合稀疏和低秩特性的欠采样重建模型的求解算法：可以利用交替方向乘子法求解(3)中的重建模型。引入变量z,Z和拉格朗日乘子C₁,C₂，其初始值均设为0。根据以下公式迭代更新变量：

当达到最大迭代次数K或s在相邻两次迭代中的误差小于设定的正数阈值μ时，迭代结束,本实验中K为10⁵，μ为10^-6。为其中，s^k+1,z^k+1,Z^k+1,和分别表示变量s，z， Z，D₁和D₂在第k+1次迭代时的值，上标H表示矩阵的复共轭转置，为将矩阵转换为向量的算子，表示奇异值收缩算子，阈值是参数β₁，β₂，τ₁，τ₂均为1。

4)由步骤3)得到恢复的扩散谱向量s，经公式(4)中算子作用得到的为最终恢复的扩散谱(如图2所示)。作为参考，将无噪全采样的信号做拉普拉斯变换得到扩散谱(如图3 所示)(Mateusz Urbanczyk,Diana Bernin,Wiktor Kozminski,KrzysztofKazimierczuk,"Iterative thresholding algorithm for multiexponential decayapplied to PGSE NMR data,"Analytical Chemistry,vol.85,no.3,pp.1828-1833,2013)。可以看出，利用采集到的部分数据和本发明的重建方法，可以重建得到高质量的扩散谱。