CN106199474B - 一种低场核磁共振二维谱反演算法 - Google Patents

Abstract

本发明涉及一种低场核磁共振二维谱反演算法，在PDCO（一种求解凸规划的原始对偶内点法）基础上推广到二维空间，用于求解二维反演问题。与目前国内外文献报道的其他低场核磁共振反演算法相比，本发明其使用L1正则化和L2正则化而不是通常的只使用L2正则化，更贴近原始信号的稀疏性，得到的解更能体现真实的谱分布；原始对偶内点法在计算时就保证了解的非负性，不需要进行额外的非负约束计算，简化了反演流程；对紧密相连的临近峰也能够区分，计算精度高；鲁棒性好，在不同信噪比数据情况下，能够得到稳定的反演结果。

Description

一种低场核磁共振二维谱反演算法

技术领域

本发明涉及一种核磁共振信号处理技术，特别涉及一种基于PDCO的低场核磁共振二维谱反演算法。

背景技术

在世界范围内，核磁共振技术发展迅猛，在很多领域得到了广泛的发展，比如高场强核磁共振技术应用于对人体的临床诊断等。随着核磁共振技术的不断发展，人们发现核磁共振技术不仅可以应用于临床，还可以在其他领域(如食品科学、农业、石油能源、材料科学、纺织化工等)发挥其他科学仪器所不能发挥的作用。很多科学研究已经证明，在这些领域中应用核磁共振分析技术，可以解决现有的其他科学仪器所不能解决的问题，对这些领域的科学进步起到了举足轻重的推动作用。与用于临床的高场核磁共振技术不同的是，在食品科学、农业、石油能源、材料科学、纺织化工等领域采用低场核磁共振技术就可以很好的解决相关问题。此外，低场核磁共振分析仪器由于采用低场强磁体，可以使所研发的分析仪器体积较小，安装、调试、维护、操作都很方便。因此，低场核磁共振分析仪器得到了科学界的关注。

CPMG(Carr-Purcell-Meiboom-Gill)序列速度快，是低场核磁共振中最常用的序列之一。研究者常常利用CPMG等序列的原始数据和样品横向弛豫时间、纵向弛豫时间的分布特点，进行时域谱反演的相关研究。但是核磁共振采集到的原始信号并不能得到直观的样品结构信息，需要通过反演技术才能得到易于理解的时域谱。一般，一维时域谱就可以为研究者分析样品的组成、性质等提供重要的依据。但是，随着低场核磁共振应用的不断深入，研究者们发现实验得到的一维谱可能存在峰的交叠。解决这一问题通常需要进行多次额外实验或者引入其它辅助手段，二维时域谱技术应运而生。二维谱不仅简化了原有流程，还可以为分析样品提供更多有价值的信息。

国内二维反演技术由于起步晚，国外相关企业的技术垄断，加上二维谱反演需要处理的数据量巨大等问题，想要发展具有极大的挑战性。目前二维反演算法，国外研究者以显式增加罚函数的正则化算法为主，对解的某些性质进行限制，以得到合理的反演结果；国内研究者最早使用基于截断奇异值分解的二维反演算法，并取得了很多成果，在信噪比(Signal to Noise Ratio,SNR)较高的情况下能够得到可靠的反演结果。现在一般采用的是标准Tikhonov正则化形式，惩罚项用于限制解的范数，这种方式的问题求解策略也相对成熟，此外还有对斜率和曲率进行平滑的正则化形式。这些正则化算法的关键在于如何选择一个合适的正则化因子，现有正则化因子的选择算法需要一定的人为干预，还待完善。

发明内容

本发明是针对核磁共振信号处理中反演技术存在的问题，提出了一种低场核磁共振二维谱反演算法，在PDCO(a Primal-Dual interior method for Convex Objectives一种求解凸规划的原始对偶内点法)基础上推广到二维空间，用于求解二维反演问题。PDCO迭代使用L1范数与L2范数，使得反演结果更贴近于真实的谱分布，并且保证了解的非负性，不再需要额外的非负约束计算。该算法能够提高反演运算的计算精度和鲁棒性，可以得到清晰、准确的二维谱。

本发明的技术方案为：一种低场核磁共振二维谱反演算法，具体包括如下步骤：

1)读取低场核磁共振设备采集得到的原始数据文件，提取出数据文件中包含的采样时间和对应时刻的采样数据M；

2)对原始数据进行横向弛豫时间、纵向弛豫时间布点，矩阵的拼接、张量积预处理操作得到反演核心矩阵K和信号幅值m；

3)计算采集数据的信噪比，根据信噪比计算反演参数，采用求解凸规划的原始对偶内点法PDCO解决反演问题表达式如下：

s.t.Ks+r＝m

s≥0

其中s＝vect(S)，vect表示将矩阵按列拼接，形成一个列向量，S(T₂，T₁)表示横向弛豫时间为T₂、纵向弛豫时间为T₁的物质的含量，反演参数λ₁＝α₁||m||₁/SNR,反演参数λ₂＝α₂/SNR，r表示拟合残差，||*||₁表示取某向量的L1范数，||*||₂表示取某向量的L2范数，信噪比SNR使用的计算方法为：采样数据信号最大值除以衰减后数据的方差，α₁与α₂参数值由大量仿真实验测得：α₁取值为0.0001，α₂为20；

4)进行迭代计算，每次迭代并保证解的非负性，得到最优解s'；

5)生成网格，绘出二维谱，用T₂、T₁这两组数据分别表示Y轴X轴信息，将这两组信息生成网格，根据T₂、T₁中元素的个数并对最优解s'进行重新排序成s”，s”表示二维谱反演结果信号的幅值，利用T₂、T₁、s”这3组数据绘出二维谱。

本发明的有益效果在于：本发明一种低场核磁共振二维谱反演算法，与目前国内外文献报道的其他低场核磁共振反演算法相比，其使用L1正则化和L2正则化而不是通常的只使用L2正则化，更贴近原始信号的稀疏性，得到的解更能体现真实的谱分布；原始对偶内点法在计算时就保证了解的非负性，不需要进行额外的非负约束计算，简化了反演流程；对紧密相连的临近峰也能够区分，计算精度高；鲁棒性好，在不同信噪比数据情况下，能够得到稳定的反演结果。

附图说明

图1为本发明低场核磁共振二维谱反演算法步骤示意图；

图2-1为本发明信噪比为1000的构造高斯峰二维谱仿真实验结果图；

图2-2为本发明信噪比为1000的CPMG数据串仿真实验结果图；

图2-3为本发明信噪比为1000的反演仿真实验结果图；

图2-4为本发明信噪比为1000的标准Tikhonov正则化算法的反演仿真实验结果图；

图3-1为本发明信噪比为100的构造高斯峰二维谱仿真实验结果图；

图3-2为本发明信噪比为100的CPMG数据串仿真实验结果图；

图3-3为本发明信噪比为100的反演仿真实验结果图；

图3-4为本发明信噪比为100的标准Tikhonov正则化算法的反演仿真实验结果图；

图4-1为本发明实验案例低浓度溶液反演结果图；

图4-2为本发明实验案例加入高浓度溶液反演结果图。

具体实施方式

下面结合附图对本发明提供的具体实施方式作详细说明。

如附图1所示，一种基于PDCO的低场核磁共振二维谱反演算法，包括如下步骤：

a.读取低场核磁共振设备采集得到的原始数据文件。

读取低场核磁共振设备采集得到的原始数据文件，提取出数据文件中包含的采样时间和对应时刻的采样数据M等信息。

b.对原始数据进行横向弛豫时间、纵向弛豫时间布点，矩阵的拼接、张量积等预处理操作得到反演核心矩阵K和信号幅值m。

二维反演问题就是求解如式(1)所示的具有两个核的Fredholm积分方程，τ₂表示各回波波峰的回波时刻，τ₁表示等待时间，S(T₂，T₁)表示横向弛豫时间为T₂、纵向弛豫时间为T₁的物质的含量，M表示特定时刻的采样数据。

将公式(1)中的弛豫时间、离散化后可以用矩阵的形式表示如式(2)所示，其中(上标i,j表示矩阵中元素位置的索引，TE表示CPMG序列中的回波时间，TW为等待时间)。

上标T表示矩阵转置，对于公式(2)，可以利用矩阵的张量积将两个核K₁、K₂合成一个核，合并后公式如(3)式所示，其中m＝vect(M)，s＝vect(S)，vect表示将矩阵按列拼接，形成一个列向量，表示两个矩阵的张量积。

m＝Ks (3)

这样转换之后，二维反演问题就变成一个已知m和K求s的问题，其中矩阵K称为反演的核心矩阵。

为了得到符合真实情况的谱分布，要对T₂、T₁进行布点，布点方法是在以10为底的对数空间内均匀布点，布点数目根据具体实验情况而定，实验精度要求越高，布点数目越多，一般为了节约计算时间将布点数目设为64。

该步骤中所述的数据预处理主要包括：T₂、T₁布点；按照公式(2)中K₁、K₂的定义计算得到K₁、K₂；利用矩阵的张量积转化K₁、K₂得到K；使用矩阵的拼接方法得到信号幅值m等操作。

c.计算采集数据的信噪比，根据信噪比计算反演参数。PDCO解决上述反演问题表达式如下：

s.t.Ks+r＝m

s≥0

其中反演参数λ₁＝α₁||m||₁/SNR,反演参数λ₂＝α₂/SNR，r表示拟合残差，||*||₁表示取某向量的L1范数，||*||₂表示取某向量的L2范数。信噪比SNR使用的计算方法为：采样数据信号最大值除以衰减后数据的方差。α₁与α₂参数值由大量仿真实验测得：α₁取值为0.0001，α₂为20(信噪比较高时，可将λ₂直接设为0.5)。

d.进行迭代计算，每次迭代并保证解的非负性，得到最优解。

迭代计算时主要求解如下问题:

s.t.Ks+r＝m:y

s-s₁＝0:z₁

其中μ,s₁≥0(μ表示惩罚因子，通常取一个很小的正参数，s₁表示松弛变量，通常第一次迭代时所有元素初始化为1)，l表示s₁中元素的位置，n表示s₁中元素总个数，由此实现累加求和。引入惩罚因子与松弛变量，保证了PDCO算法能够拓展到二维空间，并且增加了对数项保证了迭代过程中解的非负性。问题转换之后，在求解时将上式与Karush-Kuhn-Tuucker(KKT)条件联立，使用牛顿法进行搜索，即可进行迭代计算。阈值误差设置为0.001，达到终止条件时，停止迭代，得到最优解s'。

e.生成网格，绘出二维谱。用T₂、T₁这两组数据分别表示Y轴X轴信息，将这两组信息生成网格。根据T₂、T₁中元素的个数并对最优解s'进行重新排序成s”，s”表示二维谱反演结果信号的幅值，利用T₂、T₁、s”这3组数据绘出二维谱。

本发明的效果能够通过以下实验进一步说明。

1.仿真实验：

实验首先构造一个中心位于(T₁，T₂)＝(100,10)ms处的高斯峰作为理想的T1-T2谱，然后向正演结果中加入一定程度的高斯白噪声，得到不同信噪比的仿真数据。

2.仿真实验结果及结果分析

图2-1～2-4和图3-1～3-4分别是信噪比为1000(仿真低场下的高信噪比数据)和信噪比为100(仿真低场下的低信噪比数据)的仿真结果，图2-1/图3-1显示的是构造的高斯峰二维谱；图2-2/图3-2是原始数据正演得到的8条CPMG数据串；图2-3/图3-3是本发明提出的基于PDCO的二维谱反演算法反演结果；图2-4/图3-4是使用标准Tikhonov正则化算法的反演结果。

将这两组实验进行对比分析，可以看出在高信噪比环境下，2种算法都能得到反演结果，这两种算法都有一定的干扰峰出现，但明显可以看出基于PDCO的二维反演算法结果优于传统的标准Tikhonov正则化算法。在低信噪比环境下，传统的Tikhonov正则化算法，反演结果比较稳定，谱峰仍然明显被拓宽，有虚假轮廓出现；基于PDCO的混合反演算法得到的反演结果中也产生了虚假轮廓，反演结果相对于另一种算法明显更贴近仿真数据。

为了进一步验证本发明算法的抗噪性能，表1列举了信噪比分别为3000、1000、300、200、100、50下的两种算法对拟合误差比较结果。由表1可知，本发明的算法与标准Tikhonov正则化算法能够在不同信噪比数据的情况下得到稳定的结果；本发明的算法与标准Tikhonov正则化算法相比，拟合误差更小，反演结果更贴近真实的仿真数据。

表1

3.实验案例：

本实验使用NIUMAG公司的NMI-20低场核磁共振分析仪进行采样，采样序列为IR-CPMG，TW设置为0～2500ms内对数均匀分布的8个点(对应界面参数DL1)，TE设为0.235ms，累加次数NS为4。实验样品为预先制备的两种不同浓度的CuSO4溶液，分别密封在两个无核磁信号的色谱瓶中。

4.实验案例结果及结果分析：

实验首先在无磁试管中放入装有低浓度溶液的色谱瓶进行采样，使用本发明算法反演得到了图4-1所示的T1-T2谱；然后再向试管中放入装有高浓度溶液的色谱瓶(色谱瓶足够小，保证样品均处于设备的匀场区域)，得到了图4-2所示的T1-T2谱。

通常，在水中添加无磁性重金属离子会影响质子的弛豫时间，添加的离子浓度越高，采集到的CPMG信号衰减越快，质子的弛豫时间越短，在二维谱中的位置就越靠近原点。该实验中，加入高浓度的CuSO4溶液后，在T1-T2谱中更靠近原点的位置又出现了一个谱峰，说明样品中出现了更短弛豫时间的成分，进而证明了无磁性重金属离子对质子影响的论断。

Claims (1)

1.一种低场核磁共振二维谱反演算法，其特征在于，具体包括如下步骤：

1)读取低场核磁共振设备采集得到的原始数据文件，提取出数据文件中包含的采样时间和对应时刻的采样数据M；

2)对原始数据进行横向弛豫时间、纵向弛豫时间布点，矩阵的拼接、张量积预处理操作得到反演核心矩阵K和信号幅值m；

3)计算采集数据的信噪比，根据信噪比计算反演参数，采用求解凸规划的原始对偶内点法PDCO解决反演问题表达式如下：

s.t.Ks+r＝m

s≥0

其中s＝vect(S)，vect表示将矩阵按列拼接，形成一个列向量，S(T₂，T₁)表示横向弛豫时间为T₂、纵向弛豫时间为T₁的物质的含量，反演参数λ₁＝α₁||m||₁/SNR,反演参数λ₂＝α₂/SNR，r表示拟合残差，||*||₁表示取某向量的L1范数，||*||₂表示取某向量的L2范数，信噪比SNR使用的计算方法为：采样数据信号最大值除以衰减后数据的方差，α₁与α₂参数值由大量仿真实验测得：α₁取值为0.0001，α₂为20；

4)进行迭代计算，每次迭代并保证解的非负性，得到最优解s'；

5)生成网格，绘出二维谱，用T₂、T₁这两组数据分别表示Y轴X轴信息，将这两组信息生成网格，根据T₂、T₁中元素的个数并对最优解s'进行重新排序成s”，s”表示二维谱反演结果信号的幅值，利用T₂、T₁、s”这3组数据绘出二维谱。