低场核磁共振二维谱反演技术及其应用

聂生东，周小龙，王远军

(上海理工大学医疗器械与食品学院 上海 杨浦区 200093)

.生物电子学.

低场核磁共振二维谱反演技术及其应用

聂生东，周小龙，王远军

(上海理工大学医疗器械与食品学院 上海 杨浦区 200093)

针对二维谱技术的特点和优势，将其引入到不同的应用领域中，验证了二维谱技术的可靠性和实用性。首先回顾了一种以将拟合误差控制在与噪声相当的水平为基本目标的二维反演算法，然后通过仿真实验验证了该算法的准确性，最后使用该算法获取的二维谱进行了一系列的应用分析。应用案例包括造影剂浓度选择、工业明胶与食用明胶鉴别以及储层产能评估等，涉及到医学成像、食品药品、能源等多个领域。仿真实验和多个应用案例的结果表明，该方法具有很好的鲁棒性和准确性，能够适用不同信噪比数据的反演，具有实际应用价值。

二维相关谱; 二维反演; 低场核磁共振; 时域核磁共振

低场核磁共振二维谱技术是当前核磁共振领域的一个研究热点，通过使用该技术，不同领域的研究者们进行了大量的应用尝试。文献[1]利用二维谱技术来对现有生物燃料制作工艺进行指导和评价，为寻找新的生物燃料的制作来源提供了一种全新的可靠工具。文献[2]使用时域二维谱技术来评价苹果水心病病变程度。此外，二维谱还用于富铁陶瓷的烤制工艺评价[3]、生物膜通透能力与贵金属的微生物富集技术评价[4]等不同的领域。作为分析测试的新工具，低场核磁共振二维谱技术在一定程度上推动了相关学科的发展。

低场二维谱技术不同于传统的磁共振技术。磁共振成像是Fourier变换的典型应用，图像是通过对相位编码、频率编码之后的数据进行Fourier变换之后得到。一般来说主磁场强度越高，能够获取的图像质量就越高，因此磁共振成像仪的发展趋向于使用更高场强的主磁场。液体核磁共振分析技术中使用的也是Fourier变换，这种技术需要使用较高强度的主磁场，以体现出不同分子结构之间在电子云分布状况上的微小差异对共振频率产生的影响，高分辨率液体核磁共振分析仪一般都是高场设备。高分辨率固体核磁共振技术引入了魔角旋转(magic angle spinning, MAS)等方法来消除固体分子的各向异性，以提高化学位移谱的分辨力，后续处理仍需使用Fourier变换。由此可见，Fourier变换是目前广泛使用的核磁共振技术中基本的数据处理方式，而在低场二维谱技术中，取而代之的是Laplace变换。化学位移谱是典型的频域谱，低场二维谱技术获取的则是与分子动态特性(纵向弛豫时间T1、横向弛豫时间T2、扩散系数D等)相关的时域谱。

高分辨率核磁甚至可以提供与分子结构一一对应的“指纹谱”，这是高场设备高信噪比条件下的优势，但在某些应用条件下，这同时也是劣势。如在食品、农业等领域中，获取自由水、结合水和处于这两种状态之间的水的相关信息对储存、分拣等常规流程具有重要的参考价值，显然不同状态的水的“指纹”是完全一样的。在生物医药、生命科学等领域中，获取分子团簇的大小的信息同样具有很大的实用价值，当组成团簇的物质成分比较单一时，使用化学位移谱极易形成峰的混叠。如果使用低场核磁共振时域谱技术，那么可以直观地通过峰的分布情况来获取这些信息。

通过二维序列进行采样得到的原始信号一般不能直接使用，反演之后得到的二维谱才具有明确、直观的含义。从理论上来说，高场设备上也能获取时域二维谱。但是，高场设备对磁场均匀性的要求更高，还可能在具有铁磁性粒子的样品中诱发强内部梯度，严重影响采样数据的可靠性。然而，使用低场设备的采样信号的信噪比很低，加之Laplace变换对噪声的敏感性[5]，这使得稳定、可靠地进行二维谱反演具有了很大的挑战性。

本文介绍了二维反演的概念以及基于噪声拟合的二维反演算法，通过大量实验验证了算法的可靠性和二维谱技术的广泛应用前景。

1 二维反演简介

以IR-CPMG序列测T1-T2谱为例，采样数据可用具有两个核的第一类Fredholm积分方程表示[6]：

式中，M(τ2,τ1)表示极化时间为τ1、回波时刻为τ2时的信号幅度；M0为初始信号幅度；S(T2,T1)表示横向弛豫时间为T2、纵向弛豫时间为T1的磁性核的含量。显然，在积分的影响下，只能使用数值解来尽可能地进行逼近真实解。通过对T1和T2在某个范围进行布点之后，式(1)可以改为矩阵形式：

其中，有：

vect表示将矩阵按列拼接为列向量。至此，二维反演问题转化为式(3)所示的已知m和K求s的问题。

若实验采集了8条回波个数只有4 096的回波串，两个维度均布128个点，在普通32位机上仅仅将K读入内存就需要消耗2 GiB的空间。由此可见，二维反演处理的数据量非常大。同时，矩阵K一般具有很大的条件数，这造成了反演问题的病态性。

为了解决数据量过大的问题，需要对采样数据进行压缩。非等间隔下采样是一种常用的压缩方式，但是这种方式需要人工调节的参数较多。利用矩阵奇异值在前端聚集的特点，还可以通过对奇异值进行截断来进行压缩。截断位置可以用预设条件数[6]或者直接使用矩阵的秩来确定，能够实现自动压缩，本文也采用这种数据压缩方式。

关于问题的病态性，一般通过正则化处理来解决。目前常用的反演算法大都在正则化方法上有所不同，主要包括子空间(奇异值空间)正则化和使用全空间数据的罚函数正则化[7]。文献[8]实现了基于截断SVD的二维反演方法，截断位置由数据的SNR来确定。文献[9]进一步研究了最佳截断位置与SNR之间的关系并建立了二者的模型。文献[10]进行了进一步改进，引入了LSQR求取初始解。文献[11]将反演谱的熵的相反数作为罚函数来进行正则化，并对传统最大熵法的求解方法进行了改进和完善。

在二维反演中，基于SVD的反演算法极易产生震荡，在二维谱中形成虚假轮廓。罚函数法求解速度慢，且正则化因子的确定方法是基于遍历搜索(L曲线法、广义交叉验证等)的方式，较为费时。因此，二维反演算法还存在一定的发展空间。

2 二维谱反演算法

通过对不同的反演算法进行分析和评价，本文选择了一种高效的基于噪声拟合的二维反演算法，文献[12]对算法原理和算法流程进行了详细的推导，后文将进行简要的回顾。

2.1 基本思想

这种算法的基本思想是将拟合误差控制在与噪声相当的水平。由于噪声的存在，采样数据本身就包含了一些误差，因而在拟合时也通常无法得到误差为0的解。如果不存在其他因素(如采样过程中存在基线漂移等现象)，那么误差将完全由噪声引起，因而将拟合误差控制在与噪声水平相当的量级具有其合理性。

2.2 算法框架

本文采用的二维反演技术框架如图1所示。

噪声提取部分首先采用了小波变换进行滤波，使用sym8小波进行4级分解，然后利用全局阈值对小波系数进行软阈值化和信号重构。通过对比滤波前后信号幅度的变化，可以提取出“噪声”，进而估算出噪声水平。

数据压缩部分使用截断SVD进行压缩，截断位置由矩阵的秩来确定，实现自动压缩。

经过数据压缩之后问题的规模得到大幅缩减，可以提高运算效率，但二维反演往往还是一个高度病态的问题，具有很大的条件数。为了得到稳定且惟一的二维谱，需要进行正则化处理。该算法采用了标准Tikhonov正则化，有：

式中，m、s分别是对矩阵M、S的元素按列拼接之后的向量。带约束的问题式(4)转化为无约束问题[12]，则有：

显然，式(5)是一个凸规划，具有惟一的极值点。

1) 给定α一个偏大的初始值α0,设定最小正则化因子Thresh和最小变化率TOL；

2) 根据给定的α0，使用共轭梯度法求解计算如果α>α0，使用α0=0.5α0进行更新；如果α<α0，α0=α；

3) 当满足以下终止条件中的任何一个时，停机；否则，转至步骤2)。

终止条件如下：

① α≤Thresh ；

② 迭代变化速率

③ 拟合误差与噪声水平相当；

④ α不再变化。

为了进一步提高算法效率，在步骤2)中一般要先对原问题进行条件预优，也可以在迭代过程中结合非精确线搜索来加速收敛。

3 实验与结论

3.1 仿真实验

本文以T1-T2二维谱反演为例进行算法验证，首先构造一个中心在(T2,T1) = (10,100) ms的高斯峰作为预设的T1-T2谱，进行正演得到8条CPMG回波数据串，如图2所示。在此基础之上添加具有不同方差的高斯白噪声，形成不同SNR的数据。

SNR按照CPMG串中的最大幅度与噪声标准差之比进行计算。该实验构造了SNR分别为100、10、5、1等的数据，如图3所示。图中显示了多条CPMG回波串在“τ2-信号量”平面上的投影。从图3可知，当SNR较高时，本文的方法能够得到准确的二维谱，可以进行定量分析。当SNR较低时，本文的方法仍然能够得到很好的反演结果。但是，当SNR=1时，反演谱的分辨率降低(波峰变宽)，而且周围产生了一些虚假峰。实际上，当SNR过小时，任何反演方法得到的反演谱中不确定性都很大(峰越宽，对T1、T2的范围限制也越宽松，不确定性越大)，单纯对反演算法进行改进无法改善反演谱的质量，更好的方法是通过更可靠的采样方式来获取更高SNR的数据。

3.2 应用案例

核磁共振成像造影剂能缩短组织在外磁场作用下的共振时间、增大不同对比信号之间的差异，从而提高成像的对比度、清晰度。核磁共振成像结果是待检组织的T1和T2与成像参数(如回波间隔TE、恢复时间TR等)共同作用的结果，且在理论上具有明确的数学模型。通过在二维谱中观察T1、T2的分布范围，可以为成像参数的选择提供依据，从而得到更清晰的图像。核磁造影一般是通过向目标靶区注入具有一定浓度的金属离子的来改变待检部位的T1、T2值分布(弛豫率)，以增加感兴趣内容与其余内容之间的差异。在一般情况下，需要经过不断尝试使用不同浓度的造影剂，对比成像结果的优劣，然后才能通过主观判断来选择合适的造影剂浓度。如果使用T1-T2二维谱，那么可以直接检测得到不同浓度造影剂的T1、T2分布情况，进而从理论上计算出合适的造影剂浓度。本文的实验以Mn2+造影剂为例来说明通过T1-T2二维谱来帮助选择造影剂浓度的基本流程。图4是将质量体积比分别为10%和15%的Mn2+造影剂进行T1-T2二维谱检测的结果，实验使用的不同浓度造影剂分别密封在不同的色谱瓶中。取出一个色谱瓶后重新获取T1-T2二维谱，根据二维谱中峰的变化情况，可以明确图4所示二维谱中的两个峰所对应的造影剂浓度，得到如下结论：Mn2+造影剂的质量体积比越大，T1、T2弛豫时间越短。在二维谱中，还能直接看出造影剂的T1、T2值分布范围。通过观察造影剂的二维谱，可以直观、准确地对造影作用(对氢质子弛豫率的改变幅度)进行定量分析。

胶囊囊衣的主要成分是明胶，毒胶囊主要是因为在制造囊衣的食用级明胶中少量或者大量添加了廉价工业明胶。由于普通胶囊与毒胶囊组成成分不同，因此二者的T1-T2二维谱也存在差异。本文的实验使用NIUMAG的核磁共振变温分析仪(型号为VT23-010V-T)对某种市购工业明胶和食用明胶进行二维实验。实验在60℃恒温条件下进行(无需溶胀等前处理)，获取了纯样(纯食用明胶、纯工业明胶样品)的T1-T2二维谱，如图5所示。结果表明，相同质量的工业明胶和食用明胶，前者产生的总信号量比后者要大，且工业明胶的二维谱比食用明胶的二维谱多一个峰。对不同混合质量比的样品进行二维谱分析后发现，可以将左侧特征峰作为是否参杂了工业明胶的标志，还能依据特征峰所占比例来进行定量分析。

油和水的弛豫时间分布反演很广，在油水共存的储层岩心中，其T1谱和T2谱通常交叠在一起。但是，油和水在扩散系数上存在显著的差异，因此，可以将扩散系数D作为第二依据来获取D-T2二维谱。由于无法获取真实的岩样，本文的实验在人工岩心中添加了0.1 g某粘度的植物油与0.3 g一定浓度硫酸铜溶液来模拟油水共存的环境，采样仪器为纽迈MR-MD-25，反演算法仍为本文基于噪声拟合的新方法。图6中，水平实线为气线，水平虚线为水线，斜线为油线，这3条线可以通过温度、压强等参数分别计算得到。图中水线穿过的峰为水峰，油峰也刚好被油线穿过。对水、油峰分别积分可以得到二者的信号量之比为71.01%:24.38%。前期的定标实验发现单位质量水和单位质量植物油产生的信号量之比约为1:0.99，进而推断出该岩心中水油比约为3:1，与实际相符。由此可见，使用D-T2二维谱能够在油水共存的岩心环境下准确地计算出含油量。

4 结 束 语

低场核磁共振二维谱反演技术在近年来一直是国内外研究热点，针对数据量大、高度病态两个难题，研究者们进行了不断的尝试和改进。同时，反演过程中涉及大量数值算法，相关理论的完善也能促进二维谱反演技术的发展和壮大。

低场核磁共振二维谱技术拥有着广阔的应用前景，在很多需要分析分子结构和动态特性的领域中都具有应用潜力。该技术目前已经能够满足常规实验的需求，但是对于致密多孔介质(致密页岩等)，由于短弛豫组分与探头死时间相当接近，目前的二维谱技术还不具有能够还原出超短弛豫组分的能力。只有算法(软件)和硬件共同完善才能促使低场核磁共振二维谱技术的进一步发展。

[1] BERMAN P, LESHEM A, ETZIONY O, et al. Novel 1H low field nuclear magnetic resonance applications for the field of biodiesel[J]. Biotechnology for Biofuels, 2013, 6(1): 55.

[2] MELADO-HERREORS A, BARREIRO P, VILLAVALVERDE P, et al. H HR-MAS NMR metabolomic and non-destructive 2D NMR relaxometry to assess internal quality in apples[EB/OL]. [2013-05-11]. http://www. Insidefood.eu/insidefood_web/uk/word/proceedings.

[3] CASIERI C, LUCA F, NODARI L, et al. Effects of time and temperature of firing on fe-rich ceramics studied by Mössbauer spectroscopy and two-dimensional 1H-nuclear magnetic resonance relaxometry[J]. Journal of Applied Physics, 2012, 112: 1-15.

[4] VOGT S. Nuclear magnetic resonance studies of biological and biogeochemical process[D]. Bozeman：Montana State University, 2013.

[5] MITCHELL J, CHANDRASEKERA T, GLADDEN L. Numerical estimation of relaxation and diffusion distributions in two dimensions[J]. Progress in Nuclear Magnetic Resonance Spectroscopy, 2012, 62: 34-50.

[6] SONG Y Q, VENKATARAMANAN L, HURLIMANN M, et al. T1-T2correlation spectra obtained using a fast two-dimensional Laplace inversion[J]. Journal of Magnetic Resonance, 2002(154): 261-268.

[7] 周小龙, 聂生东, 王远军, 等. 核磁共振二维谱反演技术综述[J]. 波谱学杂志, 2013, 30(2): 293-305. ZHOU Xiao-long, NIE Shen-dong, WANG Yuan-jun, et al. 2D NMR inversion method: a review[J]. Chinese Journal of Magnetic Resonance, 2013, 30(2): 293-305.

[8] 顾兆斌, 刘卫. 核磁共振二维谱反演[J]. 波谱学杂志，2007, 24(3): 311-319. GU Zhao-bin, LIU Wei. The inversion of two-dimensional NMR map[J]. Chinese Journal of Magnetic Resonance, 2007, 24(3): 311-319.

[9] LIN F, WANG Z W, LI J Y, et al. Study on algorithms of low SNR inversion of T2spectrum in NMR[J]. Applied Geophysics, 2011, 8(3): 233-238.

[10] TAN M J, ZOU Y L, ZHOU C C. A new inversion method for (T2, D) 2D NMR logging and fluid typing[J]. Computers & Geosciences, 2013, 51: 366-380.

[11] CHOUZENOUX E, MOUSSAOUI S, IDIER J, et al. Efficient maximum entropy reconstruction of nuclear magnetic resonance T1-T2spectra[J]. IEEE Transactions on Signal Processing, 2010, 58(12): 6040-6051.

[12] 周小龙, 聂生东, 王远军, 等. 基于噪声拟合的核磁共振二维谱反演新方法[J]. 仪器仪表学报, 2013, 34(7): 1640-1649. ZHOU Xiao-long, NIE Sheng-dong, WANG Yuan-jun, et al. A new noise-fitting based 2D NMR inversion method[J]. Chinese Journal of Scientific Instrument, 2013, 34(7): 1640-1649.

编 辑 黄 莘

Inversion Technology of 2D NMR Relaxometry in Low Field and Its Applications

NIE Sheng-dong, ZHOU Xiao-long, and WANG Yuan-jun
(School of Medical Instrument and Food Engineering, University of Shanghai for Science and Technology Yangpu Shanghai 200093)

This paper explores the practical ability of 2D nclear magnetic resonance (NMR) relaxometry inversion based on the features and advantages of 2D maps. An inversion algorithm aims to make the misfit comparable to noise level was cast in light firstly. Then some synthetic data was handled with this algorithm to get the accurate 2D relaxometry spectrum. Lastly, this technology was applied to many practical cases in the field of medical imaging, food and medicine, energy and so on. The results of synthetic data and experimental data show that the proposed method has a good robustness and an exact accuracy. The 2D NMR relaxometry including the inversion technology has a great potential value of practical applications.

2D correlations; 2D inversion; low-field NMR; time-domain NMR

TP202+.1

A

10.3969/j.issn.1001-0548.2015.03.025

2013 − 12 − 09；

2014 − 06 − 20

国家重大科学仪器设备开发专项项目(2013YQ17046303)；国家自然科学基金(60972122); 上海市教育委员会科研创新项目(14ZZ135)

聂生东(1962 − )，男，博士，教授，主要从事医学图像处理与信号处理等方面的研究.