蒋川东 常星 孙佳 李天威 田宝凤

1)(吉林大学仪器科学与电气工程学院,长春 130026)

2)(地球信息探测仪器教育部重点实验室(吉林大学),长春 130026)

基于L1范数的低场核磁共振T2谱稀疏反演方法∗

蒋川东1)2)常星1)孙佳1)李天威1)田宝凤1)2)†

1)(吉林大学仪器科学与电气工程学院,长春 130026)

2)(地球信息探测仪器教育部重点实验室(吉林大学),长春 130026)

(2016年9月13日收到;2016年11月18日收到修改稿)

低场核磁共振技术(LF-NMR)以其无损、非侵入、原位和绿色等优势被广泛应用在食品、农业、能源和化工等行业,尤其是在食品安全监管领域发挥着越来越重要的作用.在油品品质检测中,常规非负奇异值分解(SVD)弛豫(T2)谱反演方法,只能反映光滑模型的T2谱,对于稀疏模型的反演结果存在较大差异,从而导致T2谱反演分辨率低和品质分析不准确的问题.针对这一问题,本文提出基于L1范数最小化约束的T2谱稀疏反演算法,建立NMR回波曲线的稀疏模型表达式,利用截断牛顿内点法求解L1范数最小化问题,得到稀疏模型的T2谱反演结果.通过构造光滑模型的T2谱、以及不同峰值数和信噪比的稀疏模型的T2谱,对比非负SVD算法和L1稀疏算法的反演效果,得到当信噪比大于20 dB时,L1稀疏算法能精确反演多峰T2谱,峰值幅度和峰位置均优于非负SVD算法结果.最后通过多组煎炸油样品进行低场核磁共振检测实验和不同信噪比数据的反演结果对比,验证了L1范数稀疏反演算法的准确性和优越性.

低场核磁共振,T2谱反演,稀疏模型,L1范数最小化约束

1 引 言

食品安全不仅与广大人民群众的身体健康和生命安全密切相关,而且关系到社会稳定和经济发展.自2010年至今,地沟油、瘦肉精、塑化剂超标、镉大米等事件频发,食品质量和食品安全问题已成为社会反映强烈的热点问题.面对不法食品的泛滥,亟需有效和可靠快速的检测方法监测食品健康,保障食品安全.自1946年发现核磁共振(nuclear magnetic resonance,NMR)现象以来,核磁共振技术已发展成为一种重要的分析和检测技术,并广泛应用于农业食品、能源勘探、高分子材料和生命科学等领域[1-3].以0.5 T的磁场强度为界将核磁共振现象分为高场核磁共振(high field NMR,HF-NMR)和低场核磁共振(low field NMR,LF-NMR),前者主要用来分析物质的化学性质,而后者用来检测物质的物理特性.相对于HF-NMR,LF-NMR的优势在于仪器成本低、体积小,能够实现原位、实时和快速测量[4],同时对检测样品具有非破坏性、非侵入性和无污染等特点,在食品安全监管领域发挥着越来越重要的作用,对于劣质油脂检测、注水肉检测、牛奶掺假、致病菌和重金属离子的检测等具有显著的应用潜力[5].然而,由于LF-NMR磁场强度较低,导致检测分辨率不高,因此通过改进反演算法来提高有害食品成分的检测分辨率和准确性,对于LF-NMR方法能否有效用于食品安全检测方面至关重要.

目前,国内外针对LF-NMR技术进行煎炸油品质测定,主要是基于T2谱分析方法.文献[6]在实验室内对经过60 h煎炸的大豆油进行取样,得到LF-NMR检测的弛豫图谱(简称T2谱)峰值与豆油的酸价、黏度、吸光度和极性组分含量等具有良好的相关性,由此可以评价豆油的油品品质.文献[7,8]对掺兑了煎炸油的油样进行了检测,从LF-NMR的T2谱中可以发现煎炸油的特征峰,且随掺入比例增加而峰值幅度线性增大.然而,由于LF-NMR的磁场强度低,当样品中煎炸油含量少或受噪声影响时,反演得到的T2谱分辨率有限,与煎炸油含量的线性度差,甚至难以区分煎炸油的特征峰.因此,在LF-NMR技术中亟需解决低磁场强度情况下如何提高T2谱反演分辨率的难题.

近20年来,国内外关于NMR反演算法的研究取得了重要进展,提出了包括非负最小二乘法[9]、罚函数法[10]和非负奇异值分解法(singular value decomposition,SVD)等[11]算法.很多学者对这些算法进行了优化,文献[12]等提出了一个基于L2范数最小约束的正则化算法,提高了T2谱分布的连续性,但是由于光滑参数的作用会造成T2谱的畸变.文献[13]通过对系数矩阵添加阻尼项来优化非负SVD分解法,在一定的误差容忍范围内可以接受,混有噪声时得不到理想的结果.文献[14]等通过建立信噪比和奇异值最佳保留个数的关系来改进非负SVD的截断,但需要事先计算信噪比(signal to noise ratio,SNR).可见,非负SVD算法是较为常用的反演算法,但是通过该算法得到的T2谱一般都是连续分布的,对于单一匀质样品或只含有几个离散T2值的样品,非负SVD反演结果过于光滑,不能反映出离散或稀疏的T2谱.文献[15]结合非负最小二乘法和非线性拟合方法进行多指数反演,提高了T2谱的分辨率,然而该方法将T2值也作为反演参数,增加了反演结果的不确定性.文献[16]在此基础上,提出利用线性回归最小二乘方法,改进了非线性优化算法,加快了收敛速度.但是此类算法需已知T2谱线个数,其结果依赖于初值和SNR,当SNR较低时,反演结果存在较大的误差.因此,在未知谱线个数和低信噪比情况下,如何获得稳健的稀疏结果是目前T2谱反演算法急需解决的难点之一.

稀疏表示(sparse representation)方法是信号处理、图像恢复和无线通信等领域研究的热点问题,在信息提取和参数分析等方面具有复杂信号稀疏化、特征差异明显和噪声影响小等优势[17].通过建立NMR回波信号在T2谱上的稀疏表示模型,求解L0范数最小化约束问题获得稀疏系数,即T2谱的幅度.但是L0范数最小化是一个NP难问题,进而转化为L1范数的凸优化问题,因为L1范数最小化存在惟一的全局最优解.因此,随着信号稀疏表示方法的不断发展,L1范数最小化问题的求解算法备受关注,如贪婪算法[18]和迭代收缩算法[19]可以较好地解决这一优化问题,但通常受到数据量和模型维度的限制;同伦算法[20]和内点法[21]可用于处理大型的高维数据,但需要大量的计算时间和迭代次数.截断牛顿内点法[22]是在内点法基础上,采用预处理共轭梯度(preconditioned conjugate gradient,PCG)算法计算搜索步长,减少了计算时间,从而可以快速地处理大型数据,常应用于光学分子影像、无线电通讯、人脸表情识别等方面.文献[23]提出L1范数可应用于NMR回波数据反演,并从数学上证明了迭代软阈值法、最大熵法和最小面积法都等价于最小L1范数反演问题.文献[24]将L1和L2范数同时应用于T2谱反演,改善了峰值失真问题,但该方法需要优化选择正则化参数,且反演结果介于光滑和稀疏模型之间.因此,目前利用L1范数进行低场NMR稀疏反演的研究较少,且T2谱的反演仍存在低信噪比数据反演效果差和分辨率低等问题.

本文以提高低场核磁共振T2谱分辨率为目标,将稀疏理论引入T2谱反演方法,提出了一种基于L1范数最小化约束的T2谱稀疏反演算法(以下简称L1稀疏算法),通过仿真NMR回波数据对比分析了非负SVD算法和L1稀疏算法在光滑和稀疏T2谱模型下的反演结果,最后进行了煎炸油样品的低场核磁共振实验和T2谱反演,验证基于L1范数最小化约束的T2谱稀疏反演算法的准确性和优越性.

2 低场核磁共振原理

低场核磁共振技术通过施加射频脉冲使样品中的氢质子1H发生共振,脉冲停止后1H释放能量,其磁化矢量M不断衰减,其中横向磁化矢量Mxy的衰减过程称为横向弛豫,将Mxy衰减到37%的时间为横向弛豫时间T2,T2及其对应的幅值构成了T2谱,通过对样品T2谱的相关信息进行分析,得到检测样品的品质.

2.1 低场核磁共振测量序列

低场核磁共振实验通常发射两种测量序列：自旋回波(SE)序列和CPMG(Carr-Purcell-Meiboom-Gill)序列.SE序列是先发射90◦脉冲,在半回波间隔τ时刻施加一个180◦脉冲,并在2τ时刻获得自旋回波信号的峰值,通过改变τ值,分别施加90◦和180◦脉冲,即可获得多个自旋回波信号.CPMG序列在SE序列基础上,多次施加了180◦脉冲,从而得到多个衰减的自旋回波信号.CPMG序列及自旋回波信号示意图如图1所示,在t=0时刻施加一个90◦脉冲,经过τ时间后施加第一个180◦脉冲,可在t=2τ时得到第一个回波信号;在t=3τ时继续施加第二个180◦脉冲,在t=4τ时刻得到第二个回波信号.依此类推,每在t=(2n-1)τ时刻施加一个180◦脉冲,在t=nτ时刻即可得到一个回波信号,这些回波信号的峰值是一条按指数衰减的曲线,即NMR回波曲线.

图1 低场核磁共振CPMG序列测量和NMR回波曲线原理图Fig.1.Principle diagram of CPMG sequence and NMR echo curve in LF-NMR measurement.

NMR回波曲线的函数表达式为

其中τ为半回波间隔时间,n为回波个数,y(t)为t时刻所对应的回波信号的峰值幅值,m为弛豫分量个数,xi为各弛豫分量的幅值(i=1,2,···,m);T2i为各弛豫分量的衰减时间,T2i和xi一一对应,构成T2谱.

2.2 T2谱反演方法

T2谱反演是根据NMR回波曲线,通过设定适当的各弛豫分量的衰减时间T2i,计算得到系数矩阵(加其中),从而建立一个y=Ax形式的线性方程组.通过反演算法求解该方程组,即可获得衰减时间T2i所对应的弛豫分量.

目前,T2谱反演常用的方法是非负SVD法[11].通过非负SVD算法反演得到的T2谱是连续分布的,如图2所示.图2中仿真了一个T2谱的光滑模型,布点数为36,按指数分布在区间1—104ms.设回波间隔2τ=12 ms,回波数n=1500.根据(1)式计算系数矩阵A和回波信号y,如图2(a)所示.利用非负SVD算法进行反演,结果如图2(b)所示,其中黑色线为仿真的T2谱模型,蓝色线为非负SVD算法反演得到的T2谱,两者符合较好,仅在T2较小时存在一定的误差.因此,非负SVD算法对于这种光滑模型的T2谱反演效果较好.但是,对于单一匀质样品或只含有几个离散T2值的样品,非负SVD反演结果过于光滑.

图2 (网刊彩色)光滑模型非负SVD反演结果 (a)NMR回波曲线;(b)T2谱Fig.2.(color online)Non-negative SVD inversion result of a smooth model：(a)NMR echo curve;(b)T2spectrum.

3 L1范数稀疏反演算法

基于L1范数最小化的T2谱稀疏反演算法将弛豫分量T2幅值的L1范数进行最小化约束,从而转化为一个凸优化问题,并利用截断牛顿内点法求解,得到T2谱稀疏模型下的最优解.

3.1 稀疏模型的L1范数最小化约束

在低场核磁共振的实际应用中经常需要处理一些稀疏的T2谱,因此在分析和处理NMR回波信号时,建立如下的稀疏模型：

其中‖x‖0是x的L0范数,即x向量中非零元素的个数;x=[x1,x2,···,xm]T是字典A的稀疏系数,表示未知的弛豫分量幅值,L0范数最小化约束保证了x中仅有少量元素不为零;y=[y1,y2,···,yn]T表示NMR回波信号;矩阵A=(a{i,j}){n×m}为稀疏模型的字典.(2)式为n个含有m个未知数x1,x2,···,xm的线性方程组.

由于表达式(2)中的L0范数最小化问题是一个NP难问题,而L1范数是L0范数的最优凸近似,而且更容易优化求解,考虑到实际观测信号中必然会存在噪声,因此根据表达式(2)建立基于L1范数最小化约束的含噪稀疏模型：

3.2 截断牛顿内点法

截断牛顿内点法是在内点法的基础上[22],建立特定的对数障碍函数和中心路径,并通过预处理共轭梯度(PCG)算法,按一种截断机制来计算牛顿系统的近似解,从而减少计算时间.

表达式(3)中的L1范数最小化问题是含有不等式约束的线性规划问题(linear program,LP)[25].根据正则化的思想,并设变量w=Ax-yobs,可将(3)式转化为二次规划问题,作为内点法的原问题(P)：

其中,fp(x)为原目标函数.

根据拉格朗日乘子法可得到相应的拉格朗日函数和对偶问题(D),如表达式(5)和(6)所示：

其中,fd(v)为fp(x)的对偶函数,构造的对偶可行点为τ为比例常数.原目标函数与对偶函数的差值为对偶间隙gap=fp(x)-fd(v),通过最大化对偶函数fd(v),可得到对应的最优解v∗和最小的对偶间隙.

根据截断牛顿法,将表达式(4)中的原问题P转化为一个等价的凸优化问题：

不等式约束-ui≤xi≤ui的对数障碍函数φ(x,u)可表示为

因此,原问题P的无约束形式为

函数ψt(x,u)是一个光滑、有下界的凸函数,有惟一的最小值,在计算其最小值的过程中,产生了一条受参数t影响的曲线,被称为中心路径.参数t的定义为

其中r＞1和smin∈(0,1]是预设的参数,s=βk为回溯线搜索的步长,k为满足回溯搜索条件(11)的最小回溯次数,α,β为预设参数.

利用PCG算法求解如下的牛顿系统(12),便可得到搜索方向：

根据截断牛顿内点法,解决L1范数最小化约束问题(3)的算法步骤如下：

3.3 反演结果

构造一个T2谱的稀疏模型,布点数为36,指数分布于1—104ms,如图3所示.回波间隔2τ=12 ms,回波数n=1500,通过计算得到NMR回波曲线,如图3(a)所示.图3(b)中黑色星线为仿真模型,红色线为利用L1稀疏算法进行反演得到的T2谱.通过对比可以看出,L1稀疏算法得到的T2谱与构造的稀疏模型基本一致,各点的误差均较小.因此,L1稀疏算法对于稀疏模型下的T2谱反演效果较好.

图3 (网刊彩色)稀疏模型和L1稀疏反演结果 (a)NMR回波曲线;(b)T2谱Fig.3.(color online)L1 sparse inversion result of a sparse model：(a)NMR echo curve;(b)T2spectrum.

4 仿真模型反演分析

为了验证L1稀疏算法和传统的非负SVD法对于光滑模型和稀疏模型的反演效果,对比两种方法在处理不同T2峰值数的稀疏模型时反演准确度,以及验证不同信噪比下的稀疏模型反演结果,本文建立了多组仿真模型进行反演实验.

在实验过程中,涉及到信噪比(SNR)的计算公式为

其中,y为根据(1)式得到的模拟理想回波信号,yobs为加入一定比例噪声之后的含噪回波信号.

设构造的T2谱模型中T2峰的峰值数为k,峰值为h0k,对应的弛豫时间为布点数(弛豫分量个数)为m,各弛豫分量的幅值为x0i.采用算法对信号yobs进行反演得到T2谱的峰值为hk,峰值对应的弛豫时间为各弛豫时间点对应的幅值为xi,则峰值误差Δh、峰位置误差ΔT2以及幅值平均误差Δx的计算公式分别为：

4.1 光滑模型与稀疏模型反演对比

首先分别构造光滑模型和稀疏模型的T2谱,布点数均为36,按指数分布在1—104ms,回波间隔2τ=1.2 ms,回波数n=1500,通过计算得到NMR回波曲线,分别加入信噪比为20 dB的高斯白噪声,如图4(a)和图4(b)所示.利用非负SVD算法和L1稀疏算法分别进行T2谱反演,其结果如图4(c)和图4(d)所示.通过对比可以看出：对于光滑模型,非负SVD算法和L1稀疏算法反演结果均与仿真的T2谱基本相符,只在较小的T2处,由于τ值较大,两种方法的反演结果均存在偏差.对于稀疏模型,通过L1稀疏算法得到的T2谱比非负SVD算法更接近仿真模型.这是因为非负SVD算法假设T2谱是连续分布的,对于稀疏模型的反演结果过于光滑.因此L1稀疏算法既适用于光滑模型,也同样适合应用于稀疏模型.

图4 (网刊彩色)光滑模型与稀疏模型下反演结果对比 (a)光滑模型NMR回波曲线;(b)稀疏模型NMR回波曲线;(c)光滑模型T2谱;(d)稀疏模型T2谱Fig.4.(color online)Comparison of the inversion results between the smooth and sparse models：(a)NMR echo curve of the smooth model;(b)NMR echo curve of the sparse model;(c)T2spectrum of the smooth model;(d)T2spectrum of the sparse model.

图5 (网刊彩色)T2峰值数变化模型的反演结果对比 (a)单峰值T2谱;(b)双峰值T2谱;(c)三峰值T2谱;(d)四峰值T2谱Fig.5.(color online)Comparison of inversion results of the multimodal model with different peak numbers：(a)Single-peakT2spectrum;(b)two-peakT2spectrum;(c)three-peakT2spectrum;(d)four-peakT2spectrum.

4.2 T2峰值数变化的反演对比

为了验证两种方法对不同T2峰值数的反演效果,模拟了T2谱峰值数从单峰值变化到四峰值的四个稀疏模型,如图5中黑色星线所示,并分别使用非负SVD算法和L1稀疏算法进行反演,反演结果如图5中蓝色线和红色线所示.

从图5中可以看出,L1稀疏算法在处理不同T2峰值数的稀疏模型时,其反演得到的T2谱都更接近于仿真模型,T2谱的峰值与仿真模型存在较小的误差,对应的T2值与仿真模型基本相等.在图5(d)中,使用L1稀疏算法得到的T2谱四个峰的幅值平均误差Δx仅为8.6%,T2位置与仿真模型相等.而采用非负SVD算法反演的T2谱均为光滑的,且T2谱的峰值误差Δh和峰位置误差ΔT2较大.尤其是四峰值T2谱,非负SVD算法得到的T2谱的幅值平均误差Δx和峰位置误差ΔT2分别为65.7%和19.93 ms,不能反映各峰的峰值大小和位置,而且第四个T2峰已经不明显.因此,可以得出L1稀疏算法反演结果T2峰值及幅度均明显优于非负SVD算法.

4.3 不同信噪比的反演结果

为了验证L1稀疏算法在不同信噪比(SNR)下的反演效果,在稀疏模型的T2谱计算的NMR回波曲线中,分别加入信噪比为50,30,20,15,10及5 dB的高斯白噪声,依次对不同信噪比下的NMR回波曲线利用L1稀疏算法进行反演.为了验证不同噪声数据的统计结果,在不同SNR下随机产生100组噪声数据,将非负SVD算法和L1稀疏算法的反演结果分别描绘在图6中,颜色较深的区域代表反演结果出现的机率较大,颜色较浅的区域代表出现的机率较小,黑色曲线代表构造的T2谱.统计100次反演结果的T2谱平均误差见表1.

图6 (网刊彩色)信噪比不同时的T2谱反演结果对比Fig.6.(color online)Comparison of the inversion resultsT2spectrum with different SNRs.

表1 信噪比不同时的T2谱反演结果误差Table 1.Error analysis of the inversion results ofT2spectrum with different SNRs.

通过对比图6中各SNR下的T2谱反演结果和仿真模型,可以得到：在信噪比较高时L1稀疏算法得到的T2谱更接近于仿真模型,非负SVD算法反演T2谱则一直保持光滑;随着SNR降低,噪声对反演结果的影响逐渐增大,L1稀疏算法结果出现的范围逐渐变大,当SNR小于15 dB时,L1稀疏算法得到的T2谱趋于光滑,与非负SVD算法反演结果类似.从表1中可以看出,当SNR=20 dB时,L1稀疏算法反演T2谱的幅值平均误差Δx为10.45%,而非负SVD算法结果的平均误差Δx为40.16%.SNR越大,峰值误差Δh、峰位置误差ΔT2以及幅值平均误差Δx越小,即L1稀疏算法反演T2谱越准确.在实际应用中应该考虑采取有效预处理手段将信噪比提高至20 dB以上,以确保反演精度.然而非负SVD算法反演T2谱的各项误差始终较大,当SNR=50 dB时,T2谱的幅值平均误差Δx也高达40.10%,可见即使在很小的噪声下,非负SVD算法也不能准确地反演稀疏模型的T2谱.

5 实验反演实例

为了验证L1范数稀疏反演方法在实际应用中的反演效果,通过煎炸红薯制作出煎炸油样品,使用低场核磁共振分析仪器测量,得到各组煎炸油样品的NMR回波曲线,并运用两种反演方法分别对NMR回波曲线进行反演.

5.1 煎炸油低场核磁共振实验机理

低场核磁共振技术可以用来检测煎炸油品质,其原理是通过低场核磁共振仪获得煎炸油样品的NMR回波曲线yobs(t),再反演得到煎炸油的T2谱.若将煎炸油样品视为某单一组分构成的整体,其弛豫时间为T2w,对应的幅值为x∗,则yobs(t)可表示为

通过反演,可得到样品的单组分弛豫时间T2w(ms).

正常油脂由甘油三酯组成,可分为饱和脂肪酸和不饱和脂肪酸,二者的H质子所处的结构不同,不饱和脂肪酸的弛豫时间较长,因此可认为正常油脂的T2谱中第一个峰代表饱和脂肪酸,第二个峰代表不饱和脂肪酸.而经过高温精炼、煎炸等高温氧化的油脂,会产生一系列极性增大、聚合度较高的氧化产物,其弛豫时间较短,即在前两个峰前出现第三个特征峰.因此,反演得到的T2谱具有三个谱峰,分别称为T21峰、T22峰和T23峰[6],其中T21峰的峰面积占总面积的比例S21和T2谱的单组分弛豫时间T2w,与大豆油的酸价、黏度、吸光度和极性组分含量等检测指标之间存在良好的相关性,说明利用低场核磁共振检测的S21和T2w,可以有效反映煎炸油的品质情况[21].

5.2 不同煎炸油品质的T2谱反演实验

首先制作煎炸油.用金龙鱼大豆油煎炸新鲜红薯片,在煎炸炉中倒入3 L大豆油,并将油温控制在190—200◦C范围内,每次放100 g新鲜红薯片于煎炸炉中,5 min后将煎炸好的红薯片捞出,并重新放入新鲜红薯片.每4 h对大豆油取样一次,共煎炸24 h,取样7次,再从各样品中分别取出10 mL煎炸油,根据煎炸时间分为7组(0,4,8,12,16,20,24 h).

通过NMI20台式核磁共振成像分析仪(磁场强度为(0.5±0.08)T,仪器主频率为21.3 MHz),分别对7组实验样品施加硬脉冲CPMG序列,其中90◦脉冲宽度为19µs,重复采样等待时间TR为1500 ms,半回波时间τ为300µs,回波个数根据弛豫过程长度不同设定为1000—4000个.通过循环180◦脉冲,产生多个回波信号,测量得到NMR回波曲线如图7所示.从图7中可以看出,大豆油煎炸的时间越长,得到的NMR回波曲线衰减的速度越快.

图7 (网刊彩色)煎炸油样品测量得到的NMR回波曲线Fig.7.(color online)NMR echo curves from the LFNMR measurements of frying oil samples.

从煎炸油低场核磁共振实验中,得到7组样品的NMR回波曲线,分别对其使用非负SVD算法和L1稀疏算法进行反演,反演得到的T2谱如图8.对比图8中的两组反演结果发现,通过L1稀疏算法反演的T2谱明显存在三个峰值,分别对应煎炸油的T21峰、T22峰和T23峰.而通过非负SVD算法得到的T2谱均只有两个峰值,说明反演结果过于光滑,不能有效区分煎炸油的三个峰值.

同时,L1稀疏算法反演得到的T21峰的峰面积占总面积的比例S21随煎炸时间呈线性增长,单组分弛豫时间T2w随煎炸时间呈线性减小,如图9所示.L1稀疏算法结果经过线性拟合得到：yS21=0.49x+3.11和yT2w=-1.23x+107.76,拟合系数R2分别为0.994和0.995,该数值代表模型对数据的拟合程度较好.因此,根据S21和T2w与大豆油的酸价、黏度、吸光度和极性组分含量等检测指标之间的相关性,可以有效地检测出煎炸油的品质变化[6].而通过非负SVD算法得到的T2谱中,煎炸时间与S21和T2w的线性关系为yS21=0.42x+3.65和yT2w=-6.01x+176.81,拟合系数R2分别为0.971和0.965,说明非负SVD算法对于煎炸油NMR回波曲线的反演结果劣于L1稀疏算法.

图8 (网刊彩色)煎炸油样品的T2谱反演结果对比(a)非负SVD算法;(b)L1稀疏算法Fig.8.(color online)Comparison ofT2spectrum inversion results for frying oil samples：(a)Non-negative SVD algorithm;(b)L1-sparse algorithm.

5.3 不同信噪比的T2谱反演实验

为了验证L1稀疏反演方法在低信噪比数据中的优势,利用煎炸时间为4 h的油样,分别取10,5,2,1和0.5 mL的样品,获得不同信噪比的煎炸油低场核磁共振实验数据,其他实验参数设置与5.2小节中保持一致.图10给出了对应不同样品量的NMR回波曲线,为了对比分析其信噪比,将小样品量的NMR回波曲线等比例放大至与10 mL样品量相当的量级.由图10可以得到,随着样品量的减小,NMR回波曲线的信噪比逐渐降低,当样品量为0.5 mL时,经计算信噪比为17 dB(小于20 dB),NMR回波曲线受噪声影响严重失真.分析原因是低场核磁共振仪器每次测量的噪声水平不变,随样品量减少,NMR回波信号的幅度随之减小,因此信噪比降低.

图9 (网刊彩色)煎炸时间与(a)峰面积比S21和(b)单组分弛豫时间T2w的相关性Fig.9.(color online)The correlation between the frying time and(a)the peak area ratioS21,as well as(b)the single relaxation timeT2w.

图10 (网刊彩色)不同样品量下测得的NMR回波曲线Fig.10.(color online)NMR echo curves from LFNMR measurements of the frying oil samples in different volumes.

利用图10中的NMR回波曲线,分别使用非负SVD算法和L1稀疏算法进行反演计算,得到的T2谱如图11.当样品量较大时(10,5和2 mL),即SNR较高时,L1稀疏算法的T2谱反演结果能够有效区分油样的三个峰值;但是随着样品量减小(1 mL和0.5 mL),T2谱反演结果逐渐趋于光滑,尤其是当样品量为0.5 mL时(信噪比小于20 dB),T2谱反演结果只有两个峰,且峰位置出现较大的偏差.相比之下,通过非负SVD算法得到的T2谱反演结果均是光滑的,无法有效区分油样的三个峰.因此,通过不同样品量的NMR回波曲线反演结果可以证明,L1稀疏算法在SNR较低时,反演结果均优于非负SVD算法,体现了L1稀疏算法的抗干扰能力和稳健性.但是,当SNR小于20 dB时,两种算法的结果均呈现光滑特性,且无法反映煎炸油真实的T2谱,进一步验证仿真结论的正确性.

图11 (网刊彩色)不同信噪比煎炸油样品的T2谱反演结果对比 (a)非负SVD算法;(b)L1稀疏算法Fig.11.(color online)Comparison of inversion results ofT2spectrum for the frying oil samples in different SNRs：(a)Non-negative SVD algorithm;(b)L1-sparse algorithm.

6 结 论

针对常规的非负奇异值分解算法无法准确反演稀疏模型的T2谱的问题,本文提出了基于L1范数最小化约束的T2谱稀疏反演算法,建立NMR回波曲线的稀疏模型表达式,利用截断牛顿内点法求解L1范数最小化问题,实现了稀疏模型下T2谱的高分辨精确反演.

通过模型仿真及低场核磁共振实测数据反演,得到如下结论.

1)光滑模型和稀疏模型的仿真对比实验表明,L1稀疏算法与SVD算法相比,既适用于光滑模型,同时适用于稀疏模型,且具有更高的反演精度.

2)当稀疏模型的T2峰值数由单峰值变化到四峰值时,L1稀疏算法仍能得到准确的反演结果,而SVD算法结果逐渐恶化,甚至不能确定峰值个数.

3)稀疏模型下,当测量NMR曲线的信噪比从5—50 dB逐渐变化时,L1稀疏算法在20 dB以上时可以获得精确的反演结果,峰值误差Δh在10%范围以内,峰位置误差ΔT2和幅度平均误差Δx在5%范围之内,而SVD算法的反演结果在各信噪比下都无法得到精确结果.

4)7组不同煎炸油品质样品的实测数据反演结果表明,L1稀疏算法的结果优于非负SVD算法,获得的T2谱明显呈现出三个峰值,且T21峰面积比S21和单组分弛豫时间T2w随煎炸时间呈线性变化的程度较高,再根据S21和T2w与煎炸油的酸价、黏度和吸光度等参数的关系,可以有效检测煎炸油随时间变化的品质情况.不同信噪比下T2谱的反演结果表明,当信噪比降低时,L1稀疏算法的T2谱反演结果在SNR大于20 dB的条件下仍然优于非负SVD算法.

L1稀疏算法中将求取L0范数弱化为L1范数,简化了计算过程,但是结果仍然不够稀疏.后续可以考虑将L1稀疏算法和非线性拟合相结合进行联合反演,或直接进行L0范数稀疏算法研究,以期获得高分辨率甚至超分辨率的T2谱,对于LF-NMR应用于食品安全领域具有更深远的意义.本文的研究不仅可用于油品品质、鲜乳水分和地沟油辨识等食品快速检测领域,还可以推广应用到核磁测井解释、石油储层的岩心分析、复杂油气藏流体识别等地球物理勘探和石油化工等方面.

[1]Marcone M F,Wang S,Albabish W,Nie S,Somnarain D,Hill A 2013Food Res.Int.51 729

[2]Li X,Xiao L Z,Liu H B,Zhang Z F,Guo B X,Yu H J,Zong F R 2013Acta Phys.Sin.62 147602(in Chinese)[李新,肖立志,刘化冰,张宗富,郭葆鑫,于慧俊,宗芳荣2013物理学报62 147602]

[3]Wang S,Munro R A,Shi L,Kawamura I,Okitsu T,Wada A,Kim S Y,Jung K H,Brown L S,Ladizhansky V 2013Nat.Methods10 1007

[4]Dalitz F,Cudaj M,Maiwald M,Guthausen G 2012Prog Nucl.Mag.Res.Sp.60 52

[5]Shen Y G,Xiao Z Q,Chen S S,Zhang Y L,Jiang W,Lai K Q 2013J.Food Sci.Technol.31 37(in Chinese)[申云刚,肖竹青,陈舜胜,张英力,蒋伟,赖克强 2013食品科学技术学报31 37]

[6]Wang Y W,Wang X,Liu B L,Shi R,Yang P Q 2012Food Sci.33 171(in Chinese)[王永巍,王欣,刘宝林,史然,杨培强2012食品科学33 171]

[7]Zhang Q,Saleh A M,Shen Q 2012Food Bioprocess Technol.6 2562

[8]Zhu W,Wang X,Chen L 2017Food Chem.216 268

[9]Bro R,de Jong S 1997J.Chemom.11 393

[10]Butler J P,Dawson S V 1981Siam J.Numer.Anal.18 381

[11]Wang W M,Li P,Ye C H 2001Sci.China A31 730(in Chinese)[王为民,李培,叶朝辉 2001中国科学A 31 730]

[12]Chen S S,Wang H Z,Yang P Q,Zhang X L 2014J.Biomed.Eng.31 682(in Chinese)[陈珊珊,汪红志,杨培强,张学龙2014生物医学工程学杂志31 682]

[13]Li P J,Ge C,Sun G P,Chen X,Wang Y K 2010Well Logging Technol.34 215(in Chinese)[李鹏举,葛成,孙国平,陈新,王彦凯2010测井技术34 215]

[14]Lin F,Wang Z W,Liu Q H,Ding Y,Li C C 2009J.Jilin Univ.(Earth Sci.Ed)39 1150(in Chinese)[林峰,王祝文,刘菁华,丁阳,李长春 2009吉林大学学报(地球科学版)39 1150]

[15]Wang H,Li G Y 2005Acta Phys.Sin.54 1431(in Chinese)[王鹤,李鲠颖 2005物理学报 54 1431]

[16]Wu L,Chen F,Huang C Y,Ding G H,Ding Y M 2016Acta Phys.Sin.65 107601(in Chinese)[吴量,陈方,黄重阳,丁国辉,丁义明2016物理学报65 107601]

[17]Chen W C,Wang W,Gao J H,Jiang C F,Lei J L 2013Chin.J.Geophys.56 2771(in Chinese)[陈文超,王伟,高静怀,姜呈馥,雷江莉2013地球物理学报56 2771]

[18]Mallat S G,Zhang Z 1993IEEE Trans.Signal Proces.41 3397

[19]Zhou W 2013M.S.Thesis(Guangzhou：South China University of Technology)(in Chinese)[周巍 2013硕士论文(广州：华南理工大学)]

[20]Zhang L Q,Wang J Y 2008Chin.J.Eng.Geophys.5 509(in Chinese)[张丽琴,王家映2008工程地球物理学报5 509]

[21]Wei H,Sasaki H,Kubokawa J,Yokoyama R 1998IEEE Trans.Power Syst.13 870

[22]Koh K,Kim S J,Boyd S 2007J.Mach.Learn.Res.8 1519

[23]Stern A S,Donoho D L,Hoch J C 2007J.Mag.Res.188 295

[24]Berman P,Levi O,Parmet Y,Saunders M,Wiesman Z 2013Concept Mag.Res.A42 72

[25]Karmarkar N 1984Combinatorica4 373

PACS:76.60.Es,82.56.Na DOI:10.7498/aps.66.047601

Sparse inversion method ofT2spectrum based on the L1 norm for low-field nuclear magnetic resonance∗

Jiang Chuan-Dong1)2)Chang Xing1)Sun Jia1)Li Tian-Wei1)Tian Bao-Feng1)2)†
1)(College of Instrumentation and Electrical Engineering,Jilin University,Changchun 130026,China)
2)(Key Laboratory of Geophysical Exploration Equipment,Ministry of Education(Jilin University),Changchun 130026,China)

13 September 2016;revised manuscript

18 November 2016)

The technology of low-field nuclear magnetic resonance(LF-NMR)is commonly used in food,agriculture,energy and chemical sectors due to its non-destructive,non-invasive,in situ,green and other advantages.Recently,this technology played an increasingly large role in the field of food-safety supervision especially.In oil product quality testing,conventionalT2spectrum inversion methods such as the non-negative singular value decomposition(SVD)algorithm can only reflectT2spectrum in a smooth model.However,for a sparse model,the inversion result of non-negative SVD algorithm is expected to be very glossy,leading to low resolution ofT2spectrum and inaccurate analysis of sample property.To solve this problem,we propose a sparseT2spectrum inversion algorithm based on the L1 norm minimization constraint.In this paper,we establish the sparse model expression of NMR echo curve,and obtain theT2sparse spectrum inversion results based on the inner truncated Newton-point method.Furthermore,the effectiveness of L1 sparse inversion algorithm is examined by the synthetic data of the smooth model and the spare model which have different peak numbers and signaltonoise ratios(SNRs).Synthetic results show that compared with the non-negative SVD algorithm,the L1 sparse algorithm is appropriate for both the smooth model and the sparse model with higher inversion accuracy.When the number ofT2peaks in a sparse model changes from a single peak to a quad peak,the L1 sparse algorithm can still obtain accurate inversion results,while the SVD algorithm results in a gradual deterioration,and cannot even determine the peak number.Under the sparse model,when the SNR of the measured NMR curve is gradually changed from 5 dB to 50 dB,the L1 sparse algorithm at 20 dB or more can obtain accurate inversion results which have less than 10%peak error and less than 5%peak position error and amplitude average error.However,the non-negative SVD algorithm cannot obtain accurate results at each SNR.Finally,multiple sets of frying oil samples are utilized to prove the accuracy and robustness of L1 sparse inversion algorithm.Inversion results of seven sets of frying oil samples show that the L1 sparse algorithm prefers the non-negative SVD algorithm.The obtainedT2spectrum by the L1 sparse algorithm shows three peaks obviously,and theT21peak area ratioS21and the single component relaxation timeT2ware higher linear with respect to frying time than the results by non-negative SVD algorithm,which is useful for detecting the frying oil quality change.The inversion results of theT2spectrum at different SNRs are consistent with the synthetic results,i.e.,when the SNR is reduced,theT2spectrum inversion results from the L1 sparse algorithm are better than those from the non-negative SVD algorithm when SNR is greater than 20 dB.

lowfield nuclear magnetic resonance,T2spectrum inversion,sparse representation,L1-norm minimization constraint

:76.60.Es,82.56.Na

10.7498/aps.66.047601

∗国家重大科学仪器设备开发和应用专项项目(批准号：2011YQ030133)、国家自然科学基金青年科学基金(批准号：41604083,41204079,41504086)、吉林省自然科学基金(批准号：20160101281JC)和吉林大学大学生创新创业训练项目(批准号：2015651016)资助的课题.

†通信作者.E-mail：tianbf@jlu.edu.cn

*Project supported by the National Key Foundation for Exploring Scientific Instrument(Grant No.2011YQ030133),the Young Scientists Fund of the National Natural Science Foundation of China(Grant Nos.41604083,41204079,41504086),the Natural Science Foundation of Jilin Province,China(Grant No.20160101281JC),and the College Students’Innovation and Entrepreneurship Training Project of Jilin University,China(Grant No.2015651016).

†Corresponding author.E-mail：tianbf@jlu.edu.cn