引言

核磁共振（NMR）技术在工程现场的应用一直是研究者们努力实现的方向.目前，NMR已在石油测井^[1]和电力橡胶绝缘材料老化检测^{[2, 3]}方面显现出较好的实用价值.但便携性NMR仪器的磁场很难达到很高的场强和均匀度，而且受现场复杂工况（比如石油测井中的仪器振动和极端的温度、压力）的影响，NMR信号信噪比相对较低.通常通过增加重复测量次数来提升信噪比，这就意味着需要较多的测量时间，因此对于现场NMR技术来说，提高测量效率显得格外重要.要提高测量效率，测量参数的合理设置是关键.然而，目前参数设置主要依靠操作人员依经验进行多次尝试后，最终设置一个比较合适的值.这种做法费时，且对操作人员素质要求比较高.下面以自扩散系数（D₀）和纵向弛豫时间（T₁）测量为例，进行详细说明.

图 1是低场NMR中，当静态梯度磁场存在时，测量D₀的SGSE序列^[4].其中90˚脉冲与180˚脉冲的时间间隔（T_d）是可变的.在测量扩散特性时，我们逐渐增大T_d，并以T_d为横坐标、每个T_d对应的CPMG回波信号前若干个回波信号的峰值和为纵坐标（带有测量误差），可以拟合出梯度磁场下液体分子的扩散衰减曲线（图 2），从而求出作为曲线参变量的D₀.为了提高回波信号峰值的信噪比，需要重复图 1所示序列若干次并进行叠加，因此相当耗时.

目前实现“增大T_d值”的方式是均匀步进，即设置T_d的“起点值”、“终点值”、“步进数”：（1）T_d的起点值T_d Min是90˚脉冲与第一个180˚脉冲之间的时间间隔最小值，单位为μs；（2）T_d的终点值T_d Max是90˚脉冲与第一个180˚脉冲之间的时间间隔最大值，单位为μs；（3）扩散加权时间T_d的扫描步进数Steps是T_d Min和T_d Max之间细分的扫描测量步进数目.

要想测量准确，就要根据样品的扩散曲线设置相适应的T_d Min、Steps、T_d Max三个参数.实际上面对一个新样品时，测量人员对其扩散曲线是未知的，因此需要不断的改变T_d Min、Steps、T_d Max，去进行多次试测，用这些试测值来大致估计样品的扩散曲线，这种试测无疑会消耗大量时间.同时，大部分测量仪器使用者并不具有测量原理相关的背景知识，很难根据测试值作出针对性的参数调整，甚至根本无法完成测量.因此将测量过程自动化，不仅能节约测量时间、提高测量效率，还能降低测量门槛，有利于NMR测量的推广.

近些年人工智能发展迅速，广泛应用于工业自动化领域，本文介绍了一种基于蒙特卡罗算法的NMR参数智能搜索方法，可以实现T₁、D₀测量中关键参数的自动设置.通过上述介绍我们知道，测量自动化的关键在于通过算法估计真实的扩散曲线，蒙特卡罗模拟则可以有效解决该问题.

1 算法介绍 1.1 蒙特卡罗方法简介

蒙特卡罗方法（Monte Carlo Method），也称统计模拟方法，是20世纪40年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法.蒙特卡罗并不指某一具体算法，而是一类随机方法的统称.这种方法的特点是可以通过对随机采样的计算得到近似结果，随机采样越多，得到的结果越近似最优解，但无论采样多少次它一定能返回一个结果^[5].与之对应的另一类随机算法称为拉斯维加斯算法，它的特点是采样越多越有可能找到最优解，但不保证能返回这个最优解.我们通过蒙特卡罗模拟随机生成大量扩散曲线，再根据已知的测量数据排除掉不合理的扩散曲线，剩下的曲线集合就可以用于估计真实扩散曲线.

1.2 蒙特卡罗算法数学模型

前文提到目前我们一般是通过设定T_d的“起点、终点、步进数”来不断等间距改变T_d值，所以参数设置的自动化实际就是T_d值选择的自动化.本文期望机器能根据已有的测量数据自动选择最优的T_d值，同时不再拘泥于等间距改变T_d.

假设如图 2所示的任意测量曲线为：

$ y = f(t) $

(1)

t是测量数据点的横坐标T_d，单位为ms；y是纵坐标，代表归一化回波幅值. $f(t) $内含的参变量可写为向量形式$\mathit{\boldsymbol{A}} = [{a_0}, {a_1}, {a_2}, \cdots , {a_j}] $，我们要求的D₀、T₁就是其中之一，不妨设为${a_0} $.对每个设置的测量点横坐标$ {\dot t_i}$，机器返回的测量值$ {\dot y_i}$是一个含有噪声的随机变量，上标"·"表示已测得的数据.测量噪声err为白噪声，服从均值为0、方差为${\sigma ^2} $的正态分布^[6]，表示为：

$ err \sim N(0, {\sigma ^2}) $

(2)

对$f(t) $的参数向量A中的各个参量分别按一定的密度函数^[3]进行多次随机抽样，得到新的参数向量${\mathit{\boldsymbol{A}}_1}, {\mathit{\boldsymbol{A}}_2}, \cdots , {\mathit{\boldsymbol{A}}_k} $，其中${\mathit{\boldsymbol{A}}_k} = [{a_0}(k), {a_1}(k), {a_2}(k), \cdots , {a_j}(k)] $.将${\mathit{\boldsymbol{A}}_1}, {\mathit{\boldsymbol{A}}_2}, \cdots , {\mathit{\boldsymbol{A}}_k} $分别带入函数，计算参数向量为${\mathit{\boldsymbol{A}}_k} $的曲线${f_\mathit{\boldsymbol{A}}}(t) $在各测量点$ {\dot t_i}$处与真实测量数据${\dot y_i}$的差值平方和，若该值小于现有测量点个数N与测量误差的方差$ {\sigma ^2}$的乘积，写为：

$ \sum\limits_{i = 1}^N {{{[{{\dot y}_i} - {f_{{\mathit{\boldsymbol{A}}_k}}}({{\dot t}_i})]}^2} < N} {\sigma ^2} $

(3)

则认为该参数$ {\mathit{\boldsymbol{A}}_k}$合理，予以保留^[7].

对一个纵坐标t，参数向量取${\mathit{\boldsymbol{A}}_1}, {\mathit{\boldsymbol{A}}_2}, \cdots , {\mathit{\boldsymbol{A}}_k} $时可分别对应函数值${y_{1, t}}, {y_{2, t}}, \cdots , {y_{k, t}} $，这一系列函数值构成随机变量${Y_t} $.由于取值是等概率的，故${Y_t} $的期望$E{Y_t} $、方差$D{Y_t} $可分别表示如下：

$ \left\{ \begin{array}{l} E{Y_t} = \sum\limits_{i = 1}^N {{y_{i, t}}/N} \\ D{Y_t} = \sum\limits_{i = 1}^N {{{({a_i} - E{Y_t})}^2}/N} \end{array} \right. $

(4)

$E{Y_t} $、$D{Y_t} $都是随t变化的函数.真实扩散曲线在t处的函数值的估计为$\hat y $，

$ \hat y = E{Y_t} $

(5)

同时方差$D{Y_t} $越小，该估计越准确；方差越大，则表明估计值越不准确.很自然的，下一个测量点$ {\dot t_{i + 1}}$应该取在估计最不准确，也就是方差最大的地方，有：

$ {\dot t_{i + 1}} = t|D{Y_t} = {\rm{Max}}(D{Y_t}) $

(6)

$ {\rm{Max}}(D{Y_t})$表示函数$D{Y_t} $的最大值.

一般样品的T₁在1 ms以上，可取测量初始点t₁=1 ms.值得注意的是，初始点的选择对最终结果影响不大，算法会自动根据初始点返回的结果进行后续测量.

每更新一个测量点我们就做一次拟合，求得关键参数$ {a_0}$的估计值$ {\hat a_0}$.明显，随着测量点的增加，$ {\hat a_0}$将趋于稳定.当增加一个测量点对${\hat a_0} $的影响足够小，满足下式时，即可停止测量.

$ |{\hat a_{0(i)}} - {\hat a_{0(i - 1)}}| < \varepsilon $

(7)

${\hat a_{0(i)}} $表示参数$ {\hat a_0}$第i步的估计值，ε是提前设定的容差值.

整个算法流程可用图 3表示.

以下，本文将分析两种最常用的测量情景.

1.2.1 D₀测量

测量D₀曲线时，$f(t) $满足形式^[4]

$ f(t) = B*{{\rm{e}}^{ - {{(\frac{t}{{{D_0}}})}^3}}} $

(8)

此时参数向量A=[D₀, B].考虑到(8)式的形式，若不存在测量误差时，B应等于实测回波峰值和的最大值$ {\rm{Max}}({\dot y_i})$，考虑到测量误差可将参数搜寻范围设为[$ {\rm{Max}}({\dot y_i}) \pm 5\sigma $]；常见样品的D₀在0~10^-8 cm²·s^-1之间，这样不妨设D₀、B服从各自区间内的均匀分布，

$ \left\{ \begin{array}{l} {D_0} \sim U[0, {10^{ - 8}}]\\ B \sim U[{\rm{Max}}({{\dot y}_i}) - 5\sigma , {\rm{Max}}({{\dot y}_i}) + 5\sigma ] \end{array} \right. $

(9)

按上述分布进行多次随机抽样，当满足

$ {\sum\limits_{i = 1}^N {\left[ {{{\dot y}_i} - B*{{\rm{e}}^{ - {{(\frac{t}{{{D_0}}})}^3}}}} \right]} ^2} < N{\sigma ^2} $

(10)

条件的曲线足够多时，可计算$D{Y_t} $，再代入(6)式求得${\dot t_{i + 1}} $.将T_d设为${t_{i + 1}} $的值进行测量，得到新的测量值$ {\dot y_{i + 1}}$，用最小二乘法进行拟合，得到$ {\hat a_{i + 1}}$，满足(7)式时停止迭代.

1.2.2 T₁测量

测量T₁使用的IR序列如图 4所示.与SGSE序列类似，IR序列中的T_d值也是可变的.我们通过设置“T_d Min、T_d Steps、T_d Max”三个参数来不断增大T₁，重复该序列“Scan”次.

测量T₁曲线时，$ f(t)$满足形式

$ f(t) = C*{{\rm{e}}^{ - t/{T_1}}} $

(11)

常见样品的T₁在0~5 000 ms范围内，可设T₁、C服从分布：

$ \left\{ \begin{array}{l} {T_1} \sim U[0, 5\, {\kern 1pt} 000]\\ C \sim U[{\rm{Max}}({{\dot y}_i}) - 5\sigma , {\rm{Max}}({{\dot y}_i}) + 5\sigma ] \end{array} \right. $

(12)

按上述分布进行多次随机抽样，从中筛选出满足下式的曲线计算$D{Y_t} $

$ \sum\limits_{i = 1}^N ( {\dot y_i} - C*{{\rm{e}}^{ - t/{T_1}}}{)^2} < N{\sigma ^2} $

(13)

将计算得到的$D{Y_t} $代入(6)式可得到${\dot t_{i + 1}} $，迭代终止条件同样为(7)式.

2 实验验证

为了验证上述算法的有效性，我们用该算法取点测量纯水的扩散系数D₀和T₁.将算法得出的${t_i} $输入可编程谱仪，实测得到${\dot y_i} $，再将$ {\dot y_i}$返回Matlab计算出$ {t_{i + 1}}$，其流程如图 5所示.

设定误差上限ε= 10^-10，测量12个点后算法停止迭代.随着测量数据增加，预测数据集不断减小，给出已知数据分别为1、5、8、12个时候的预测数据集如图 6所示.实测数据绘出的扩散曲线如图 7所示，最终拟合得到D₀ = 2.447×10^-9 cm²·s^-1.

使用同样的算法测量纯水的T₁，设定误差上限ε= 1，测量13个点后算法停止迭代.已知数据个数分别为1、5、8、13时候的预测数据集如图 8所示.测出T₁曲线如图 9所示，拟合得到T₁=2 754.1 ms.

考虑到温度影响，可以认为该算法下测得的纯水的D₀、T₁参数基本准确.

3 结论

本文通过蒙特卡罗模拟对样品的扩散和弛豫曲线进行估计，算法根据估计值选取最优测量参数，实现了D₀、T₁测量时参数设置的自动化，节省测量时间的同时还大幅降低了测量门槛，利于NMR技术的推广.经验证该算法能够得到较为准确的测量结果.值得注意的是，本文所提的人工智能搜索算法是建立在单一组分样品的前提下，要想适用于多组分样品还需要对上述算法进行扩展，这也是本团队后继的研究方向.

致谢 作者徐征感谢与Schlumberger公司宋一桥博士的讨论，以及宋一桥博士在思路上提供的帮助.