时域核磁共振谱反演的方法
【专利摘要】本发明涉及一种时域核磁共振谱反演的方法，首先，生成反演核并形成L1范数正则化问题：先根据采样序列和采样参数，生成反演核，添加L1范数项进行惩罚，形成L1范数正则化问题；然后，数据拟合：对L1范数正则化问题进行迭代求解，并得到反演谱。大幅提高了反演算法的执行效率与谱的分辨率，同时还具有很好的鲁棒性。
【专利说明】
时域核磁共振谱反演的方法
技术领域
[0001] 本发明涉及一种核磁共振领域中信号处理技术，特别涉及一种时域核磁共振谱反 演的方法。
【背景技术】
[0002] 时域谱是与化学位移谱相对的另一种核磁共振图谱，时域谱的横坐标通常为横向 弛豫时间T2,纵向弛豫时间Tl等。目前时域核磁共振分析技术已经在能源、食品、农业等领 域中得到了广泛应用。
[0003]时域核磁共振设备采集的数据并不能直接使用，需要通过对实验数据进行反演才 能得到时域谱，因此时域谱也可以称为反演谱。目前广泛使用的时域谱为T2谱、Tl谱等一维 谱，但是随着应用的深入，研究者发现一维谱中可能存在峰的交叠问题，二维谱技术应运而 生。二维反演问题具有两个核，通过矩阵重排可以合并为一个核，因而二维反演问题可以转 化为一个等价的一维反演问题。但是，跟一维反演相比，二维反演需要处理的数据量更大， 反演问题的不适定性也越明显。
[0004] 在不考虑数据规模的情况下，二维反演可以转化为一维反演的形式。然而，反演问 题具有不适定性，微小的噪声干扰都可能引起反演谱的巨大变化，因而需要对原始的反演 问题进行正则化。目前广泛使用的正则化类型是L2范数正则化，这种正则化方式要求拟合 误差小的同时反演谱的L2范数也相对较小，从而限制反演谱的复杂性。L2范数正则化的缺 点在于容易引起过度平滑，还有可能擦除小比例的峰。Ll范数正则化能直接限制反演谱的 稀疏性，产生更易于解释的反演谱，还能更容易地还原出真实存在的小比例峰，因而在反演 问题中具有很大的应用潜力。
【发明内容】
[0005] 本发明是针对现在时域核磁共振谱反演存在的问题，提出了一种时域核磁共振谱 反演的方法，可在保证反演谱准确性的前提下，提高反演谱的分辨率、反演算法的效率和鲁 棒性。
[0006] 本发明的技术方案为:一种时域核磁共振谱反演的方法，具体包括如下步骤：
[0007] 1)生成反演核K并形成Ll范数正则化问题:先根据采样序列和采样参数，生成反演 核K，添加 Ll范数项进行惩罚，形成Ll范数正则化问题；
[0008] 2)数据拟合:对Ll范数正则化问题进行迭代求解，并得到反演谱。
[0009] 2、根据权利要求1所述时域核磁共振谱反演的方法，其特征在于，所述步骤1)中生 成反演核K为单个二维矩阵，如采样序列为一维时域核磁共振的实验所得，那就是一维反演 问题，生成的反演核K为单一的二维矩阵，一维反演包括横向弛豫时间T2谱反演、纵向弛豫 时间Tl谱反演、扩散系数D谱反演、内部梯度场G谱反演中任意一项;如采样序列为二维时域 核磁共振的实验所得，存在两个任意一维反演核，对这两个反演核重排到一个二维矩阵中， 作为最终的反演核K。
[0011]
[0010] 所述步骤1)中形成LI范数正则化问题，为了得到反演谱时域谱s，确定采样参数之 后，根据反演核K、采样数据m和正则化因子λ确定一个Ll范数正则化问题：
[0012] I I S I I 丨里，表不向量S中的所有兀素均不小 于0〇
[0013]所述步骤2)数据拟合步骤为：
[0014] I、使用Lipschitz常数对Ll范数正则化问题进行近似表达；
[0015] II、迭代求解的过程为：
[0016]①初始化Lipschitz常数L、正则化因子λ，将中间变量ti初始化为1，即titer=i = l;
[0017] ②设置初始解so，将初始解备份至yi，；
[0018] ③迭代步数记为iter，iter不小于1，每次迭代后的解存入yiter中。PL表示右侧以y 为因变量，L为亚变量的最小值问题，第i ter次迭代的过程为：
[0019；
[0020；
[0021；
[0022] III、迭代的终止条件为以下两条条件中的任意一项：
[0023]第一条:拟合误差达到预设的范围；
[0024]第二条:迭代次数达到预设值。
[0025] 本发明的有益效果在于:本发明时域核磁共振谱反演的方法，大幅提高了反演算 法的执行效率与谱的分辨率，同时还具有很好的鲁棒性。
【附图说明】
[0026] 图1为本发明时域核磁共振谱反演的方法流程图；
[0027]图2为本发明时域核磁共振谱反演的方法中二维采样序列PGSE-CPMG的脉冲时序 图；
[0028] 图3为本发明时域核磁共振谱反演的方法中预设的典型双峰一维T2谱模型图；
[0029] 图4为使用本发明所提出的时域核磁共振谱反演的方法对仿真数据进行一维T2谱 反演的结果图；
[0030] 图5为使用L2范数正则化方法对仿真数据进行一维T2谱反演的结果图；
[0031] 图6为本发明时域核磁共振谱反演的方法中预设的典型T1-T2二维谱；
[0032] 图7为使用本发明方法对仿真数据进行二维T1-T2谱反演的结果图。
【具体实施方式】
[0033]如图1所示时域核磁共振谱反演的方法流程图，包括步骤:第一步，生成反演核K， 形成Ll范数正则化问题:先根据采样序列和采样参数，生成反演核，然后，添加 Ll范数项进 行惩罚，形成Ll范数正则化问题;第二步，数据拟合:对Ll范数正则化问题进行迭代求解，并 得到反演谱。具体步骤展开如下：
[0034] 1、生成反演核K，形成Ll范数正则化问题：
[0035] 确定采样序列和实验参数之后，K就被确定。一维时域核磁共振的实验过程，完成 了 m=Ks中m的采集，而一维反演就是根据采样数据m和反演核K来求解时域谱s的过程，其中 m和s均为向量。一维反演问题的反演核K是一个二维矩阵，不同的采样序列具有不同的反演 核形式。
[0036] 表1列举了常见一维采样序列的反演核形式，其中CPMG序列和IR序列名称相对统 一，而表1提及的PGSE序列在不同厂商的设备中可能具有不同的名称。为了避免混淆，本发 明所述的PGSE脉冲时序如图2中的窗1所示。PGSE序列在SE序列的基础上，往180度射频(RF) 脉冲左右两侧施加了两个相同的梯度，这两个梯度通过控制扩散作用的程度来调节相散程 度，用于对扩散系数D进行编码，故称为相敏梯度。图2中的窗1所示的PGSE序列中，δ为相敏 梯度的持续时间（两个相敏梯度的持续时间相等），A为两个相敏梯度的间隔时间，tgap表 示90度射频脉冲与第一个相敏梯度之间的时间间隔，这个间隔等于第二个相敏梯度与回波 波峰出现时刻之间的时间间隔，因此，回波时刻T〇 = 2(tl+S)+ Δ。
[0037]现以PGSE序列测D为例来介绍一维反演核的合成方法：首先，对D进行布点得到向 量D=%，^ .. .D1,... ,DndID1为合理的扩散系数猜测值，用ND个离散的值去逼近真实 的D谱，因而ND通常设置为一个较大的数;然后，确定实验中改变的实验参数，假设改变的参 数为G，一共变化了NG次，那么将有向量G= [GhGhG3, .. .Gi，...，GNG]T，在这M次采样中，Δ 和δ保持不夺：最后，枏据采样序列所确定的反演核形式，形成NG X ND的反演核
，其中γ表示质子的旋磁比，T表示矩阵转置。
[0038]表1
[0040]二维时域核磁共振实验需要使用两个反演核进行描述，采样完成的足M ::: 中M的采集，二维反演是已知M、K2和Ki这三个二维矩阵求二维矩阵S(二维时域谱）的过程。K2 和心可以根据采样序列类型和采样参数，并通过对未知谱在合理范围之内进行布点之后直 接生成。通过特定的方式，可以将这两个核合并为一个核。二维谱是一维谱的拓展，意义相 同，用相同字母的大小写来区别，小写为一维谱，大写为二维谱。另外本文所有的黑粗体小 写均为向量，黑粗体大写均为二维矩阵。
[0041] 对于核可分离的情况，可以使用Kronecker积将两个核直接合并。以图2所示PGSE- CPMG序列为例，窗1确定的反演核记）
其中
f数意义同前；窗2确定的反演核记，
其 为实验参数，常数；T〇 = 2(tgap+5)+A，r=l，2,3, · · ·， NE，NE为窗2的回波个数;c=l，2,3，.. .，NT2，NT2表示T2的布点个数。由于窗1调节的参数是 G，其余参数保持不变，改变G之后，窗2的反演核不受任何影响（保持不变），因此Kd和1(2为 可分离的独立核。此时，可以将:喊=餐.鼻1<丨重写为111 = 1(8，只需使1( = 1〇'011(1(2，1^)，111 = ￥6(^ (M)，S = Vect(S)，这里的kron表示矩阵的Kronecker积，vect表示将二维矩阵按列拼接为列 向量。
[0042] 对于核不可分离的情况，需要按照特定的顺序进行矩阵元素堆叠与重排。仍以图2所 示的PGSE-CPMG序列为例，假设在窗1中改变的参数是△，一共使用了N个不同的△值，其余参 数保持不变，窗1确定的反演核变仍记为1<?
.~15%3,...5瘦1窗2确定的反演核仍）
其中
但此时，每改变一次△的值，τ〇都将改变，两个核耦合在一起。 对于这种情况，每一个A值对应一个二维矩阵K2(A1),可以直接将Kd的第i行Κ^(?，·表 示第i行所有列）与Κ2( Λ 〇按照kronecker积合并为一个二维矩阵K1，最后将K1拼接到IT1的 下方，保证得到的矩阵的列数维持不变(NDXNT2列）。进行N次重复操作之后，就能得到二维 矩阵K，使_K :SK；·重写为m = Ks，其中m = vect(M)，s = vect(S)。不管是一维实验还是 二维实验，最终都可以转化为相同的表达形式m=Ks。
[0043] 生成的反演核K为单个二维矩阵;对于包括横向弛豫时间T2谱反演、纵向弛豫时间 Tl谱反演、扩散系数D谱反演、内部梯度场G谱反演等在内的一维反演问题，反演核K为单一 的二维矩阵;对于包括T1-T2、D-T2、T2-T2、T2-G等在内的二维反演问题，存在两个反演核， 在这种情况下对这两个反演核重排到一个二维矩阵中，作为最终的反演核K。
[0044] 为了避免产生过于复杂的反演谱，保证反演谱的可解释性，现对反演谱的稀疏性 进行限制，假定真实谱只具有有限的信息，这与目前的公知相符。使用Ll范数作为惩罚项进 行正则化之后，得到正则化问题：
[0045]
(1)
[0046]式（1)中，I |s| I1表示s的LI范数;s>0表示向量s中的所有元素均不小于0;λ为正 则化因子，λ越大，正则化程度越高。
[0047] 2、数据拟合：
[0048] 为了便于描述，记 f(s) = I I Ks-m I |2，g(s)=A| |s| I i，seRn，Rn 表示 η 维向量空间， X、y是该向是奋I、m由的而冬占 -Vf妄5従辟m)所示的最小值问题可以表示为
[0049] (2)
[0050] 显然，f (X)符合Lipschitz条件，那么，对于定义域内的x，yeRn，有如下表达式成
立：
[0051 ] (3)
[0052]式中，L(f)是▽/的Lipschitz 常数。
[0053] 现进一步对F(s) = f (s )+g( s)进行近似表达。对于任意一个Lipschitz常数L>0， 在给定点V处，可以将F( s)做如下近似：
[0054]
(4)
[0055] Ql(s，x)具有单一极小值点pL(y)=arg min{QL(s，y)}，忽略只与X相关的常数项之 后，该极小值点可以进一步简化为
[0056]
(5)
[0057] 可见，可以通过如下基本的格式对式（1)所示问题进行迭代求解：
[0058] Siter = PL(Siter-I) (6)
[0059] 结合前述步骤，本发明数据拟合部分的步骤如下：
[0060] I.使用Lipschitz常数对Ll范数正则化问题进行近似表达；
[0061 ] II.迭代求解的过程为：
[0062]①初始化Lipschitz常数L、正则化因子λ，将中间变量ti初始化为l(titer=i = l);
[0063] ②设置初始解so，将初始解备份至yi;
[0064] ③迭代步数记为iter，iter不小于1，每次迭代后的解存入yiter中。PL表示右侧以y 为因变量，L为亚变量的最小值问题，第i ter次迭代的过程为：
[0068] III.迭代的终止条件为以下条件中的任意一项：
[0065]
[0066]
[0067]
[0069] ①拟合误差达到预设的范围；
[0070] ②迭代次数达到预设值。
[0071] 上述步骤将基本的迭代格式（6)进行了扩展。将f (s)、g(s)的表达式代入（5)和 (6)，就能得到迭代格式(7)。迭代格式(8)和(9)的作用是将前两次的迭代结果进行线性组 合来形成新的结果。迭代过程中s出现的负值将直接置零。
[0072] 本发明使用了 Ll范数作为惩罚项来进行正则化，能够得到具有很好的稀疏性的结 果，更易于对反演谱进行分析与解释。图3为一个典型的双峰T2谱模型，图4和图5处理的数 据均由图3所示T2谱通过计算机仿真得到。图4为使用本发明所提出的方法进行处理得到的 反演T2谱，图5为使用L2范数正则化方法得到的反演T2谱。从图4和图5中可以看出，在组分 比较接近时，使用本发明所提出的方法能够得到更准确的反演谱，而L2范数正则化方法产 生的结果中存在组分交叠。图6为一个典型的T1-T2二维谱模型，图7为使用本发明所提出的 方法进行处理得到的反演T1-T2谱，从图中可以看出，本发明所提出的方法能够很好的还原 真实谱。
[0073] 与现有技术相比，本发明一、使用Ll范数作为惩罚项来直接限制反演谱的稀疏性， 能产生更易于解释的反演谱。同时，本发明产生的结果能有效抑制峰的交叠;二、迭代求解 格式简洁，计算量非常小，能够快速得到反演谱。
【主权项】
1. 一种时域核磁共振谱反演的方法，其特征在于，具体包括如下步骤： 1) 生成反演核K并形成Ll范数正则化问题:先根据采样序列和采样参数，生成反演核K， 添加 Ll范数项进行惩罚，形成Ll范数正则化问题； 2) 数据拟合:对Ll范数正则化问题进行迭代求解，并得到反演谱。2. 根据权利要求1所述时域核磁共振谱反演的方法，其特征在于，所述步骤1)中生成反 演核K为单个二维矩阵，如采样序列为一维时域核磁共振的实验所得，那就是一维反演问 题，生成的反演核K为单一的二维矩阵，一维反演包括横向弛豫时间T2谱反演、纵向弛豫时 间Tl谱反演、扩散系数D谱反演、内部梯度场G谱反演中任意一项;如采样序列为二维时域核 磁共振的实验所得，存在两个任意一维反演核，对这两个反演核重排到一个二维矩阵中，作 为最终的反演核K。3. 根据权利要求1或2所述时域核磁共振谱反演的方法，其特征在于，所述步骤1)中形 成Ll范数正则化问题，为了得到反演谱时域谱s，确定采样参数之后，根据反演核K、采样数 据m和正则化因子λ确定一个Ll范数正则化问题：I I s| Ii表示S的Ll范数，其中m和S均为向量，s>0表示向量S中的所有元素均不小于0。4. 根据权利要求1所述时域核磁共振谱反演的方法，其特征在于，所述步骤2)数据拟合 步骤为： I、 使用Lipschitz常数对Ll范数正则化问题进行近似表达； II、 迭代求解的过程为： ① 初始化Lipschitz常数L、正则化因子λ，将中间变量ti初始化为1，即titer=i = l; ② 设置初始解so，将初始解备份至yi，； ③ 迭代步数记为iter，iter不小于1，每次迭代后的解存入ylter中，PL表示右侧以y为因 变量，L为亚变量的最小值问题，第i ter次迭代的过程为：III、 迭代的终止条件为以下两条条件中的任意一项： 第一条:拟合误差达到预设的范围； 第二条:迭代次数达到预设值。
【文档编号】G01N24/08GK105891249SQ201610187844
【公开日】2016年8月24日
【申请日】2016年3月29日
【发明人】周小龙, 苏冠群, 聂生东, 王丽嘉, 王远军, 张英力, 杨培强
【申请人】上海理工大学