CN105891249A - 时域核磁共振谱反演的方法 - Google Patents

Abstract

本发明涉及一种时域核磁共振谱反演的方法，首先，生成反演核并形成L1范数正则化问题：先根据采样序列和采样参数，生成反演核，添加L1范数项进行惩罚，形成L1范数正则化问题；然后，数据拟合：对L1范数正则化问题进行迭代求解，并得到反演谱。大幅提高了反演算法的执行效率与谱的分辨率，同时还具有很好的鲁棒性。

Description

时域核磁共振谱反演的方法

技术领域

本发明涉及一种核磁共振领域中信号处理技术，特别涉及一种时域核磁共振谱反演的方法。

背景技术

时域谱是与化学位移谱相对的另一种核磁共振图谱，时域谱的横坐标通常为横向弛豫时间T2，纵向弛豫时间T1等。目前时域核磁共振分析技术已经在能源、食品、农业等领域中得到了广泛应用。

时域核磁共振设备采集的数据并不能直接使用，需要通过对实验数据进行反演才能得到时域谱，因此时域谱也可以称为反演谱。目前广泛使用的时域谱为T2谱、T1谱等一维谱，但是随着应用的深入，研究者发现一维谱中可能存在峰的交叠问题，二维谱技术应运而生。二维反演问题具有两个核，通过矩阵重排可以合并为一个核，因而二维反演问题可以转化为一个等价的一维反演问题。但是，跟一维反演相比，二维反演需要处理的数据量更大，反演问题的不适定性也越明显。

在不考虑数据规模的情况下，二维反演可以转化为一维反演的形式。然而，反演问题具有不适定性，微小的噪声干扰都可能引起反演谱的巨大变化，因而需要对原始的反演问题进行正则化。目前广泛使用的正则化类型是L2范数正则化，这种正则化方式要求拟合误差小的同时反演谱的L2范数也相对较小，从而限制反演谱的复杂性。L2范数正则化的缺点在于容易引起过度平滑，还有可能擦除小比例的峰。L1范数正则化能直接限制反演谱的稀疏性，产生更易于解释的反演谱，还能更容易地还原出真实存在的小比例峰，因而在反演问题中具有很大的应用潜力。

发明内容

本发明是针对现在时域核磁共振谱反演存在的问题，提出了一种时域核磁共振谱反演的方法，可在保证反演谱准确性的前提下，提高反演谱的分辨率、反演算法的效率和鲁棒性。

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

1)生成反演核K并形成L1范数正则化问题：先根据采样序列和采样参数，生成反演核K，添加L1范数项进行惩罚，形成L1范数正则化问题；

2)数据拟合：对L1范数正则化问题进行迭代求解，并得到反演谱。

2、根据权利要求1所述时域核磁共振谱反演的方法，其特征在于，所述步骤1)中生成反演核K为单个二维矩阵，如采样序列为一维时域核磁共振的实验所得，那就是一维反演问题，生成的反演核K为单一的二维矩阵，一维反演包括横向弛豫时间T2谱反演、纵向弛豫时间T1谱反演、扩散系数D谱反演、内部梯度场G谱反演中任意一项；如采样序列为二维时域核磁共振的实验所得，存在两个任意一维反演核，对这两个反演核重排到一个二维矩阵中，作为最终的反演核K。

所述步骤1)中形成L1范数正则化问题，为了得到反演谱时域谱s，确定采样参数之后，根据反演核K、采样数据m和正则化因子λ确定一个L1范数正则化问题：

$\underset{s>0}{\arg }\min \{F\left(s\right)\≡||Ks-m|{|}^2+\mathrm{\λ}\left|\right|s\left|{|}_1\right\}$

||s||₁表示s的L1范数，其中m和s均为向量，s＞0表示向量s中的所有元素均不小于0。

所述步骤2)数据拟合步骤为：

I、使用Lipschitz常数对L1范数正则化问题进行近似表达；

II、迭代求解的过程为：

①初始化Lipschitz常数L、正则化因子λ，将中间变量t₁初始化为1，即t_iter＝1＝1；

②设置初始解s₀，将初始解备份至y₁，；

③迭代步数记为iter，iter不小于1，每次迭代后的解存入y_iter中。p_L表示右侧以y为因变量，L为亚变量的最小值问题，第iter次迭代的过程为：

$s_{iter}=p_L\left(y_{iter}\right)=\underset{s}{\arg \min }\left\{\mathrm{\λ}\right|\left|s\right|{|}_1+\frac{L}{2}\left|\right|s-(y_{iter}-\frac{2}{L}{K}^T({\mathrm{Ky}}_{iter}-m\left)\right)\left|{|}^2\right\}$

$t_{iter+1}=\frac{1+\sqrt{1+4t_{iter}^2}}{2}$

$y_{iter+1}=s_{iter}+\left(\frac{t_{iter}-1}{t_{iter+1}}\right)(s_{iter}-s_{iter-1});$

III、迭代的终止条件为以下两条条件中的任意一项：

第一条：拟合误差达到预设的范围；

第二条：迭代次数达到预设值。

本发明的有益效果在于：本发明时域核磁共振谱反演的方法，大幅提高了反演算法的执行效率与谱的分辨率，同时还具有很好的鲁棒性。

附图说明

图1为本发明时域核磁共振谱反演的方法流程图；

图2为本发明时域核磁共振谱反演的方法中二维采样序列PGSE-CPMG的脉冲时序图；

图3为本发明时域核磁共振谱反演的方法中预设的典型双峰一维T2谱模型图；

图4为使用本发明所提出的时域核磁共振谱反演的方法对仿真数据进行一维T2谱反演的结果图；

图5为使用L2范数正则化方法对仿真数据进行一维T2谱反演的结果图；

图6为本发明时域核磁共振谱反演的方法中预设的典型T1-T2二维谱；

图7为使用本发明方法对仿真数据进行二维T1-T2谱反演的结果图。

具体实施方式

如图1所示时域核磁共振谱反演的方法流程图，包括步骤：第一步，生成反演核K，形成L1范数正则化问题：先根据采样序列和采样参数，生成反演核，然后，添加L1范数项进行惩罚，形成L1范数正则化问题；第二步，数据拟合：对L1范数正则化问题进行迭代求解，并得到反演谱。具体步骤展开如下：

1、生成反演核K，形成L1范数正则化问题：

确定采样序列和实验参数之后，K就被确定。一维时域核磁共振的实验过程，完成了m＝Ks中m的采集，而一维反演就是根据采样数据m和反演核K来求解时域谱s的过程，其中m和s均为向量。一维反演问题的反演核K是一个二维矩阵，不同的采样序列具有不同的反演核形式。

表1列举了常见一维采样序列的反演核形式，其中CPMG序列和IR序列名称相对统一，而表1提及的PGSE序列在不同厂商的设备中可能具有不同的名称。为了避免混淆，本发明所述的PGSE脉冲时序如图2中的窗1所示。PGSE序列在SE序列的基础上，往180度射频(RF)脉冲左右两侧施加了两个相同的梯度，这两个梯度通过控制扩散作用的程度来调节相散程度，用于对扩散系数D进行编码，故称为相敏梯度。图2中的窗1所示的PGSE序列中，δ为相敏梯度的持续时间(两个相敏梯度的持续时间相等)，Δ为两个相敏梯度的间隔时间，tgap表示90度射频脉冲与第一个相敏梯度之间的时间间隔，这个间隔等于第二个相敏梯度与回波波峰出现时刻之间的时间间隔，因此，回波时刻τ₀＝2(t1+δ)+Δ。

现以PGSE序列测D为例来介绍一维反演核的合成方法：首先，对D进行布点得到向量D＝[D₁，D₂，D₃，...D_i，...，D_ND]^T，D_i为合理的扩散系数猜测值，用ND个离散的值去逼近真实的D谱，因而ND通常设置为一个较大的数；然后，确定实验中改变的实验参数，假设改变的参数为G，一共变化了NG次，那么将有向量G＝[G₁，G₂，G₃，...G_i，...，G_NG]^T，在这M次采样中，Δ和δ保持不变；最后，根据采样序列所确定的反演核形式，形成NG×ND的反演核其中γ表示质子的旋磁比，T表示矩阵转置。

表1

二维时域核磁共振实验需要使用两个反演核进行描述，采样完成的是中M的采集，二维反演是已知M、K₂和K₁这三个二维矩阵求二维矩阵S(二维时域谱)的过程。K₂和K₁可以根据采样序列类型和采样参数，并通过对未知谱在合理范围之内进行布点之后直接生成。通过特定的方式，可以将这两个核合并为一个核。二维谱是一维谱的拓展，意义相同，用相同字母的大小写来区别，小写为一维谱，大写为二维谱。另外本文所有的黑粗体小写均为向量，黑粗体大写均为二维矩阵。

对于核可分离的情况，可以使用Kronecker积将两个核直接合并。以图2所示PGSE-CPMG序列为例，窗1确定的反演核记为其中 参数意义同前；窗2确定的反演核记为其中TE为实验参数，常数；τ₀＝2(tgap+δ)+Δ，r＝1，2，3，...，NE，NE为窗2的回波个数；c＝1，2，3，...，NT₂，NT₂表示T₂的布点个数。由于窗1调节的参数是G，其余参数保持不变，改变G之后，窗2的反演核不受任何影响(τ₀保持不变)，因此K_D和K₂为可分离的独立核。此时，可以将重写为m＝Ks，只需使K＝kron(K₂，K_D)，m＝vect(M)，s＝vect(S)，这里的kron表示矩阵的Kronecker积，vect表示将二维矩阵按列拼接为列向量。

对于核不可分离的情况，需要按照特定的顺序进行矩阵元素堆叠与重排。仍以图2所示的PGSE-CPMG序列为例，假设在窗1中改变的参数是Δ，一共使用了N个不同的Δ值，其余参数保持不变，窗1确定的反演核变仍记为其中 窗2确定的反演核仍为其中但此时，每改变一次Δ的值，τ₀都将改变，两个核耦合在一起。对于这种情况，每一个Δ值对应一个二维矩阵K₂(Δ_i)，可以直接将K_D的第i行(i，·表示第i行所有列)与K₂(Δ_i)按照kronecker积合并为一个二维矩阵Kⁱ，最后将Kⁱ拼接到K^i-1的下方，保证得到的矩阵的列数维持不变(ND×NT₂列)。进行N次重复操作之后，就能得到二维矩阵K，使重写为m＝Ks，其中m＝vect(M)，s＝vect(S)。不管是一维实验还是二维实验，最终都可以转化为相同的表达形式m＝Ks。

生成的反演核K为单个二维矩阵；对于包括横向弛豫时间T2谱反演、纵向弛豫时间T1谱反演、扩散系数D谱反演、内部梯度场G谱反演等在内的一维反演问题，反演核K为单一的二维矩阵；对于包括T1-T2、D-T2、T2-T2、T2-G等在内的二维反演问题，存在两个反演核，在这种情况下对这两个反演核重排到一个二维矩阵中，作为最终的反演核K。

为了避免产生过于复杂的反演谱，保证反演谱的可解释性，现对反演谱的稀疏性进行限制，假定真实谱只具有有限的信息，这与目前的公知相符。使用L1范数作为惩罚项进行正则化之后，得到正则化问题：

$\underset{s>0}{\arg }min\{F\left(s\right)\≡||Ks-m|{|}^2+\mathrm{\λ}\left|\right|s\left|{|}_1\right\}---\left(1\right)$

式(1)中，||s||₁表示s的L1范数；s＞0表示向量s中的所有元素均不小于0；λ为正则化因子，λ越大，正则化程度越高。

2、数据拟合：

为了便于描述，记f(s)＝||Ks-m||²，g(s)＝λ||s||₁，s∈Rⁿ，Rⁿ表示n维向量空间，x、y是该向量空间中的两个点，表示梯度，式(1)所示的最小值问题可以表示为

$\underset{s>0}{\arg }min\{F\left(s\right)\≡f\left(s\right)+g\left(s\right)\}---\left(2\right)$

显然，f(x)符合Lipschitz条件，那么，对于定义域内的x，y∈Rⁿ，有如下表达式成立：

$\left|\right|\▿f\left(x\right)-\▿f\left(y\right)\left|\right|\≤L\left(f\right)\left|\right|x-y\left|\right|---\left(3\right)$

式中，L(f)是的Lipschitz常数。

现进一步对F(s)≡f(s)+g(s)进行近似表达。对于任意一个Lipschitz常数L＞0，在给定点y处，可以将F(s)做如下近似：

$Q_L(s,y)\&equiv;f\left(y\right)+<s-y,\&dtri;f\left(y\right)>+\frac{L}{2}\left|\right|s-y|{|}^2+g\left(s\right)---\left(4\right)$

Q_L(s，x)具有单一极小值点p_L(y)＝arg min{Q_L(s，y)}，忽略只与x相关的常数项之后，该极小值点可以进一步简化为

$p_L\left(y\right)=\mathrm{argmin}\{g\left(s\right)+\frac{L}{2}||s-(y-\frac{1}{L}\&dtri;f(y\left)\right)|{|}^2\}---\left(5\right)$

可见，可以通过如下基本的格式对式(1)所示问题进行迭代求解：

s_iter＝p_L(s_iter-1) (6)

结合前述步骤，本发明数据拟合部分的步骤如下：

I.使用Lipschitz常数对L1范数正则化问题进行近似表达；

II.迭代求解的过程为：

①初始化Lipschitz常数L、正则化因子λ，将中间变量t₁初始化为1(t_iter＝1＝1)；

②设置初始解s₀，将初始解备份至y₁；

③迭代步数记为iter，iter不小于1，每次迭代后的解存入y_iter中。p_L表示右侧以y为因变量，L为亚变量的最小值问题，第iter次迭代的过程为：

$s_{iter}=p_L\left(y_{iter}\right)=\underset{s}{\arg \min }\left\{\mathrm{\&lambda;}\right|\left|s\right|{|}_1+\frac{L}{2}\left|\right|s-(y_{iter}-\frac{2}{L}{K}^T({\mathrm{Ky}}_{iter}-m\left)\right)\left|{|}^2\right\}---\left(7\right)$

$t_{iter+1}=\frac{1+\sqrt{1+4t_{iter}^2}}{2}---\left(8\right)$

$y_{iter+1}=s_{iter}+\left(\frac{t_{iter}-1}{t_{iter+1}}\right)(s_{iter}-s_{iter-1})---\left(9\right)$

III.迭代的终止条件为以下条件中的任意一项：

①拟合误差达到预设的范围；

②迭代次数达到预设值。

上述步骤将基本的迭代格式(6)进行了扩展。将f(s)、g(s)的表达式代入(5)和(6)，就能得到迭代格式(7)。迭代格式(8)和(9)的作用是将前两次的迭代结果进行线性组合来形成新的结果。迭代过程中s出现的负值将直接置零。

本发明使用了L1范数作为惩罚项来进行正则化，能够得到具有很好的稀疏性的结果，更易于对反演谱进行分析与解释。图3为一个典型的双峰T2谱模型，图4和图5处理的数据均由图3所示T2谱通过计算机仿真得到。图4为使用本发明所提出的方法进行处理得到的反演T2谱，图5为使用L2范数正则化方法得到的反演T2谱。从图4和图5中可以看出，在组分比较接近时，使用本发明所提出的方法能够得到更准确的反演谱，而L2范数正则化方法产生的结果中存在组分交叠。图6为一个典型的T1-T2二维谱模型，图7为使用本发明所提出的方法进行处理得到的反演T1-T2谱，从图中可以看出，本发明所提出的方法能够很好的还原真实谱。

与现有技术相比，本发明一、使用L1范数作为惩罚项来直接限制反演谱的稀疏性，能产生更易于解释的反演谱。同时，本发明产生的结果能有效抑制峰的交叠；二、迭代求解格式简洁，计算量非常小，能够快速得到反演谱。

Claims (4)

1.一种时域核磁共振谱反演的方法，其特征在于，具体包括如下步骤：

1)生成反演核K并形成L1范数正则化问题：先根据采样序列和采样参数，生成反演核K，添加L1范数项进行惩罚，形成L1范数正则化问题；

2)数据拟合：对L1范数正则化问题进行迭代求解，并得到反演谱。

2.根据权利要求1所述时域核磁共振谱反演的方法，其特征在于，所述步骤1)中生成反演核K为单个二维矩阵，如采样序列为一维时域核磁共振的实验所得，那就是一维反演问题，生成的反演核K为单一的二维矩阵，一维反演包括横向弛豫时间T2谱反演、纵向弛豫时间T1谱反演、扩散系数D谱反演、内部梯度场G谱反演中任意一项；如采样序列为二维时域核磁共振的实验所得，存在两个任意一维反演核，对这两个反演核重排到一个二维矩阵中，作为最终的反演核K。

3.根据权利要求1或2所述时域核磁共振谱反演的方法，其特征在于，所述步骤1)中形成L1范数正则化问题，为了得到反演谱时域谱s，确定采样参数之后，根据反演核K、采样数据m和正则化因子λ确定一个L1范数正则化问题：

$\underset{s>0}{\arg }min\{F\left(s\right)\&equiv;||Ks-m|{|}^2+\mathrm{\&lambda;}\left|\right|s\left|{|}_1\right\}$

||s||₁表示s的L1范数，其中m和s均为向量，s＞0表示向量s中的所有元素均不小于0。

4.根据权利要求1所述时域核磁共振谱反演的方法，其特征在于，所述步骤2)数据拟合步骤为：

I、使用Lipschitz常数对L1范数正则化问题进行近似表达；

II、迭代求解的过程为：

①初始化Lipschitz常数L、正则化因子λ，将中间变量t₁初始化为1，即t_iter＝1＝1；

②设置初始解s₀，将初始解备份至y₁，；

③迭代步数记为iter，iter不小于1，每次迭代后的解存入y_iter中，p_L表示右侧以y为因变量，L为亚变量的最小值问题，第iter次迭代的过程为：

$s_{iter}=p_L\left(y_{iter}\right)=\underset{s}{\mathrm{argmin}}\left\{\mathrm{\&lambda;}\right|\left|s\right|{|}_1+\frac{L}{2}\left|\right|s-(y_{iter}-\frac{2}{L}{K}^T({\mathrm{Ky}}_{iter}-m\left)\right)\left|{|}^2\right\}$

$t_{iter+1}=\frac{1+\sqrt{1+4t_{iter}^2}}{2}$

$y_{iter+1}=s_{iter}+\left(\frac{t_{iter}-1}{t_{iter+1}}\right)(s_{iter}-s_{iter-1});$

III、迭代的终止条件为以下两条条件中的任意一项：

第一条：拟合误差达到预设的范围；

第二条：迭代次数达到预设值。