二维核磁共振测井正反演模拟软件设计与应用

Development and Application of 2D NMR Logging Simulation Software

【作者：张家成，张宫*， 黄若坤，覃莹瑶】

摘要

   研发了一套二维核磁共振测井正反演模拟软件，能够进行T₂-D、T₂-T₁、T₂-G等多种二维模式的正反演数值模拟。软件包括流体参数设置模块、回波采集模块和二维谱反演模块，可以构建精细的二维分布谱、模拟二维回波采集以及对二维反演进行精细研究。实现了同时模拟四种不同流体的核磁共振弛豫过程，并且可以对各类核磁共振属性进行单参量影响因素分析，通过实时可视化图形显示，可以直观的看到各种参数变化对弛豫过程的影响。研发的软件使用简单、结果准确，可以为二维核磁共振相关工作提供高效的研究手段。

0 引言

   核磁共振（NMR）测井是通过观测地层中流体的氢原子核的NMR信号，来获得用于评价地层流体的信息^[1-2]。从20世纪50年代被引入到石油勘探领域以来，核磁共振技术在饱和度、渗透率计算、流体识别、孔隙度、孔隙结构评价等方面均得到了成功应用^[3-5]。随着2002年Sun、Hurlimann等在石油测井领域引入了核磁共振波谱学中的二维核磁共振概念^[6-7]，二维核磁共振测井技术的出现进一步拓展了核磁共振测井的使用范围，提高了储层测井解释评价的准确率^[8]。核磁反演方法也从对单个回波串数据进行反演发展到同时对多个回波串数据进行反演^[9]，但大量的数据加重了反演过程中的运算负担。因此许多学者为了满足精度和效率上的要求，对二维反演方法进行了大量的研究：

   如顾兆斌等在国内首次利用改进的奇异值分解法对二维核磁共振进行反演^[10]；谢然红等针对二维的多回波串数据处理问题提出了多回波串联合反演方法^[9]；李鹏举等人则针对BRD和TSVD两种反演方法的不足提出了变参量迭代快速反演方法^[11]。但在二维反演方法研究过程中，繁琐的数据构造和结果评测过程增加了二维反演方法研究的时间成本。

   针对这一问题，本文从精度和效率出发，研发了一套二维核磁共振正反演模拟软件，实现了核磁共振二维弛豫及反演的全过程，降低了二维核磁反演方法的研究难度。

1 原理

1.1 核磁共振弛豫测量

   核磁共振测井通过观测氢原子核的NMR信号，反算出地层孔隙中流体的含量及其性质。岩石物理信息，如渗透率、孔径分布、孔隙度和束缚水也都可以通过核磁共振弛豫测量得到。纵向弛豫按照公式（1）进行：

\[M\left(t\right)=\sum_{i}{P_i\left[1-2\exp{\left(\frac{-t}{T_1}\right)}\right]}    (1)\]

   其中，M(t)表示核磁共振测量得到的纵向弛豫信号；P_i 是第i个弛豫信号占总的纵向弛豫信号的百分比；t表示时间；T₁表示纵向弛豫时间，该值反映了自旋系统的磁能从外界吸收或者释放能量的效率^[13]。T₁主要通过反转恢复脉冲序列或饱和恢复法脉冲序列进行测量。   

   横向弛豫则按照公式（2）进行：

\[M\left(t\right)=\sum_{i}{P_i\exp{\left(\frac{-t}{T_2}\right)}}    (2) \]

   其中，M(t) 表示核磁共振测量得到的横向弛豫信号；P_i 是第i个弛豫信号占总的横向弛豫信号的百分比；t表示时间；T₂表示横向弛豫时间^[13]，一般通过CPMG脉冲序列进行测量。在一系列CPMG序列之间，等待一段时间确保所有的自旋完全恢复，再准备下一个激励，这段时间称为 “等待时间”（Wait Time）。实际测量中，不同地层除了横向弛豫时间T₂不同外，纵向弛豫时间T₁、扩散系数D、内部梯度G等均存在差异，会直接影响回波信号的衰减，用一组回波串进行反演往往不能充分地挖掘流体之间的差异。二维核磁共振采用多组回波串组合测量，通过联合反演，得到二维分布谱，使不同流体在处理成果中有更好的区分度。

1.2 二维核磁弛豫构造

   构造二维弛豫谱是采用对数坐标下的二维高斯分布公式进行计算^[14],具体实现按照公式(3)：

\[A_j=F\frac{1}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}exp\left[-\frac{1}{2\left(1-\rho^2\right)}\left(\frac{\left(\ln{T_{2j}}-\ln{T_{2mid}}\right)^2}{2{\sigma_1}^2}+\frac{\left(\ln{D_j}-\ln{D_{mid}}\right)^2}{2{\sigma_2}^2}\ -\ \frac{2\rho\left(\ln{T_{2j}}-\ln{T_{2mid}}\right)(\ln{D_j}-\ln{D_{mid}})}{\sigma_1\sigma_2}\right)\right]    (3) \]

   公式（3）中，j=1,2,…,m,m是构造谱的布点数目；T_2j是第j组分的幅度；F是孔隙度刻度因子；T_2mid是第组分的T₂弛豫时间是T₂弛豫时间峰值的位置；D_j是第j组分的扩散系数；D_mid是扩散系数峰值的位置；σ₁、σ₂可以控制分布的展布宽度； ρ 是T_2j 和D_j 的相关系数。

   为了更加真实的模拟出实际二维谱，二维构造谱必须能够反映多种不同性质流体的叠加效果，所以公式（4）需要进一步改进成如下形式：

\[ A_j=\sum_{g=1}^{4}{F_g{HI}_g\frac{1}{2\pi\sigma_{1g}\sigma_{2g}\sqrt{1-\rho^2}}exp\left[-\frac{1}{2\left(1-\rho^2\right)}{\left(\frac{\left(\ln{T_{2j}}-\ln{T_{2g,mid}}\right)^2}{2{\sigma_{1g}}^2}+\frac{\left(\ln{D_j}-\ln{D_{g,mid}}\right)^2}{2{\sigma_{2g}}^2}-\frac{2\rho\left(\ln{T_{2j}}-\ln{T_{2g,mid}}\right)\left(\ln{D_j}-\ln{D_{g,mid}}\right)}{\sigma_{1g}\sigma_{2g}}\right)}\right]}\ \    (4)\]

   公式(4)中， j=1,2,…,m, m是构造谱的布点数目；设具有4类不同性质的流体，F_g是第g种流体的孔隙度分量；HI_g是第g种流体的含氢指数；T_2g,mid是第g种流体的T₂弛豫中心值；D_g,mid是第g种流体的扩散系数峰值；σ_1g、σ_2g分别是第g种流体T₂峰值的中心展布宽度和扩散系数峰值的展布宽度； ρ 是T_2j 和D_j 的相关系数。

1.3 常见的二维核磁共振脉冲序列

   针对二维核磁共振测井数据采集的特点，国内外学者均提出了不同的CMPG脉冲序列。如Hurlimann 和 Venkataramanan等提出的扩散编辑脉冲序列^[7]。该序列中时序采集分为两个窗口，第一个窗口内仅采集两个回波，通过调整这两个回波的回波间隔来采集孔隙中流体的扩散弛豫信息，而第二个窗口则与改良式CPMG脉冲序列一样，都是采用最小的回波间隔来采集孔隙流体的横向弛豫信息（如图1）。在回波采集过程中，回波的幅度衰减按照公式（5）进行。

\[b_{ik}=\sum_{j=1}^{m}\sum_{p=1}^{n}{f\left(D_p,T_{2j}\right)\exp{\left(-\frac{i \bullet D_p\left(\gamma GT_{Elk}\right)^2}{12}\right)}\exp{\left(-\frac{i \bullet T_{Elk}}{T_{2j}}\right)}}\ \ \ \ \ \ \ i=1,2\cdots\ \ \ \    (5)\]

\[b_{ik}=\sum_{j=1}^{m}\sum_{p=1}^{n}{f\left(D_p,T_{2j}\right)\exp{\left(-\frac{D_p\left(\gamma GT_{Elk}\right)^2}{6}\right)}\exp{\left(-\frac{2T_{Elk}+\left(i-2\right) \bullet T_E}{T_{2j}}\right)}}\ \ \ \ \ \ i=3,4,5\cdots\    (5)\]

   公式（5）中，b_ik表示回波间隔为T_Elk时第一个窗口内的第i回波的幅度；γ 表示氢核的旋磁比；G表示磁场梯度；T_E表示第二个窗口内的回波间隔；当扩散系数为D_p，横向弛豫时间为T_2j 时的孔隙度分量用f(D_P,T_2j) 表示。该脉冲序列通过调整第一个窗口内的回波间隔，采集多条CPMG回波信号，从而实现T₂-D模式下的二维数据采集，利用公式（5）将采集获得的回波数据进行联合反演，便可得到地层孔隙中流体的T₂-D分布。

   谢然红等人则利用回波间隔不同的多条CPMG回波串来完成T₂-D二维数据的采集，相对于前一种方法更加简便，利用现有的一维核磁共振测井仪器便可以完成该类数据的采集。

\[b_{ik}=\ \sum_{j=1}^{m}\sum_{p=1}^{n}{f\left(D_p,T_{2j}\right)\exp{\left(-\frac{i\bullet D_p\left(\gamma G\ T_{Ek}\right)^2}{12}\right)}\exp{\left(-\frac{i\bullet T_{Ek}}{T_{2j}}\right)}}\ \ i=1,2,3\cdots\ \    (6)\]

   公式（6）中，b_ik是回波间隔为T_Ek时第i个回波的幅度^[9]。

   除此之外还有通过调整第一个窗口内的延时来采集多条CPMG回波信号的T₁编辑序列^[4]以及测量多组不同等待时间的多T_W序列等其他常见的脉冲序列^[15]，这里便不再一一赘述。

1.4 二维谱反演方法

   目前常见的二维谱反演方法主要有奇异值分解法（SVD）、非负约束正则化法（BRD）、截断奇异值分解法（TSVD）等。其中SVD法是通过核函数矩阵的特征和反问题的总体特征，来移除能够反映测量误差的最小奇异值，以此实现正则化，求取合理解。顾兆斌等^[10]针对SVD法容易导致谱图不连续的问题，提出一种适用于对信噪比较高的数据进行反演的奇异值分解法。谭茂金等^[15]针对TSVD的分辨率和精度问题提出了一种先用非负最小二乘法（LSQR）进行反演，再将其结果作为初值是使用TSVD算法进行反演的混合反演算法。

   而BRD法是先将求解方程Tikhonov正则化，再 加入罚函数方程的解进行限制，使最终的解在约束条件范围内。但是BRD法反演结果的优劣取决于选用的正则化因子值是否合适^[11]。一般利用L曲线等方法来确定正则化因子，但这些方法在处理低信噪比的数据时常常出现结果没有意义或在反演时不收敛的问题。所以国内外学者提出了很多解决办法，如利用L曲线的斜率来选择正则化参数^[16]或者程晶晶等人提出的一种相对误差较小的基于SVD和BRD的正则化反演算法^[17]等。

2 软件的设计与实现

   本文模拟软件基于. NET4. 6，使用C#语言在Visual Studio2017的开发环境下进行开发。

2.1 结构设计

   软件设计结构如图2所示，整体分为四大部分：算法、应用模块、数据及图形显示。

   软件将所有算法放在统一的算法库中，独立于交互界面，便于算法的持续改进和扩展；应用模块包括二维谱构造模块，模拟回波采集模块，二维谱反演模块；数据管理模块和绘图模块同样相对独立，并专门设置有外部平台接口，用于导入实验室岩心分析数据或实际测井数据。

图3 软件主界面

   主界面分为三部分，分别为菜单栏、绘图区、参数区（如图3）。菜单栏有文件（NMR File）、参数设置（Parameters Set）和关于（About）。其中NMR File菜单主要是新建文件、保存文件、打开文件和退出；Parameters Set菜单用于打开各个模块的参数设置界面，设置参数和导入相应的数据；About菜单是关于本文软件的一些开发过程及人员信息。

2.2 功能实现

(1) 二维谱构造模块

   根据回波串组合序列完成“二维谱构造”模块，该模块是基于二维高斯分布进行设计的。通过选择相应流体类型、填写流体T₂值及其展布、分量孔隙度（各类流体含量）、含氢指数和另一维度参数，以及最后的布点方案：T₂及另一维度参数的布点最小值、最大值和布点个数，来进行二维谱的构造。

   对话框中的各个参数用户可根据实际模拟的流体进行设置（如图4）。并根据用户填写的结果自动计算出总孔隙度和核磁孔隙度。二维构造谱如图5所示。

图4 二维谱构造参数设置对话框

(2) 模拟回波串采集模块

   该模块根据上一步设置的构造二维谱，结合给定的模拟采集指令，来进行回波串的模拟采集。主要参数：回波组数，每组回波串的回波间隔（T_E）、等待时间（T_W）、回波采集个数、采集长度（NE）、信噪比（S/N）。

   运行此模块，会显示出模拟采集得到的回波数据及对应的噪音数据，用户可以根据图形实时对参数进行修订（如图6），回波采集结果如图7。同时，该模块允许用户导入实验室岩心分析数据或实际测井数据进行单点的数据采集与分析。

图6 回波采集参数设置对话框

图7 回波采集结果

(3) 二维谱反演模块

   主要参数包括：T₂谱布点参数（个数，最大值，最小值）、反演使用的回波串选择、第二维参数谱布点的参数（个数，最大值，最小值）、平滑因子以及反演方法的选择(如图8)。二维反演谱如图9所示。

图8 二维谱反演参数设置对话框

图10 二维谱回波串拟合效果

   软件内置了一种基于BRD法的快速高分辨二维谱反演方法，命名为：Default_Method。同时用户可以利用软件中提供的接口来编写自己的反演方法。

   反演完成后，该模块会自动将反演得到的二维谱再次进行正演（拟合），观测拟合结果跟采集得到的回波串是否吻合（如图10），从而判断反演效果好坏。

3 可靠性验证及应用实例

3.1 可靠性验证

   以硫酸铜溶液物理实验为参考，对正反演模拟软件的可靠性进行验证，具体实验步骤如下：

   第一步：配置三种不同浓度的硫酸铜溶液。

   第二步：各取三种不同浓度的硫酸铜溶液3ml放入三只容量较小的试管内，然后同时将三只试管放入到容量较大的试管内进行核磁共振测量。三种浓度的硫酸铜溶液的测量结果中弛豫时间分别为2500ms、80ms、2ms左右（如图11）。

   第三步：利用本文研发的正反演模拟软件采用与物理实验相同的参数进行二维谱构造、回波采集以及二维谱反演，得到数值模拟结果（如图12、图13）。

   将软件模拟硫酸铜实验得到的二维构造谱和二维反演谱与实际测量得到的二维反演谱进行对比。可以观察到数值模拟结果跟物理实验结果无论是分布大小还是分布位置都十分接近，验证了软件的可靠性。

3.2 应用实例

   研发的二维核磁共振正反演模拟软件可以进行多种二维模式的正反演数值模拟。以研究T₂-T₁二维模式下核磁共振方法在含有束缚水和油的地层中识别流体的效果^[14]为例进行模拟。

   首先，根据表1油水的特征参数构造出油水模型的理想分布谱即构造谱（如图14），然后设置不同的等待时间模拟出15组回波串数据。回波间隔（T_E）均为：1ms，等待时间（T_W）分别为10ms、30ms、50ms、70ms、100ms、300ms、500ms、700ms、1000ms、2000ms、4000ms、8000ms、10000ms、12000ms、15000ms。

   同时利用软件的反演模块将模拟数据进行反演处理。通过对比二维反演谱与构造谱，可见油水的分布位置基本重合（如图15）,表明了软件对于T₂-T₁模式下的油水模型的正反演数值模拟具有较高的可信度。

4 结论

   （1）研制的二维核磁共振测井正反演模拟软件可以进行T₂-D、T₂-T₁、T₂-G等多种二维模式的正反演数值模拟，并且其模拟结果与实际岩心核磁共振仪器实验结果非常接近。

   （2）二维核磁共振测井正反演模拟软件可以模拟多达四种不同的核磁共振弛豫过程，并且能够对各类核磁共振属性进行单参量因素分析,包括横向弛豫时间（T₂）、纵向弛豫时间（T₁）、含氢指数（HI）、扩散系数（D）等。

   （3）二维核磁共振测井正反演软件考虑了噪声背景下各类参数模拟，并且可以将实验数据样本和汇总的实际测井数据样本进行处理。

   （4）研发的二维核磁共振测井正反演模拟软件，为二维核磁共振相关研究工作提供了高效的技术手段。

参 考 文 献

1. Coates G R, Xiao L Z, Prammer M G. NMR Logging Principles and Applications[M]. Houston:Gulf Publishing Company,USA,2000.
2. Dunn K-J, Bergman D J, Latorraca G A. Nuclear Magnetic Resonance: Petrophysical and Logging Applications. Handbook of Geophysical Exploration[M].Pergamon, NewYork.2002.
3. 肖立志．核磁共振成像测井与岩石核磁共振及其应用[M]. 科学出版社, 1998.
4. 邹良志, 李伟, 刘清华. 二维核磁共振测井方法[J]. 国外测井技术, 2011(6):18-22.
5. 谢然红, 肖立志, 邓克俊,等. 二维核磁共振测井[J]. 测井技术, 2005, 29(5):430-434.
6. Sun B , Dunn K J . Core analysis with two dimensional NMR[J]. Meteorology & Atmospheric Physics, 2002, 109(1):61-72.
7. Hurlimann, Freedman. Diffusion-Editing: New NMR measurement of Saturation and Pore Geometry[C]// Society of Pertrophysicists & Well Log Analysts Symposium. 2002.
8. 崔海标, 陆凡, 李鑫, 等. D-T₂二维核磁共振测井技术发展现状综述[J]. 石油化工应用, 2016(3):1-5.
9. 谢然红, 肖立志, 刘家军, 等. 核磁共振多回波串联合反演方法[J]. 地球物理学报, 2009, 52(11).
10. 顾兆斌, 刘卫. 核磁共振二维谱反演[J]. 波谱学杂志, 2007, 24(3):311-319.
11. 李鹏举, 谷宇峰. 二维核磁共振变参量迭代快速反演方法[J]. 地球物理学进展, 2015(2):628-635.
12. 王忠东, 肖立志, 刘堂宴. 核磁共振弛豫信号多指数反演新方法及其应用[J]. 中国科学:物理学 力学 天文学, 2003, 33(4):323-332.
13. 彭川,何宗斌,张宫.核磁共振T₂弛豫仿真软件研发及在岩心核磁实验中的应用[J].计算机应用与软件, 2019, 36(03):2 82-286+300.
14. 谢然红, 肖立志, 陆大卫. 识别储层流体的(T₂,T₁)二维核磁共振方法[J]. 测井技术, 2009, 33(1):26-31.
15. 谭茂金,邹友龙.(T₂,D)二维核磁共振测井混合反演方法与参数影响分析[J].地球物理学报,2012,55(02):683-692.
16. 周小龙, 聂生东, 王远军, 等. 核磁共振二维谱反演技术综述[J]. 波谱学杂志, 2013, 30(2):293-305.
17. 程晶晶, 吴磊, 宋公仆. 基于SVD和BRD的二维核磁共振测井正则化反演算法研究[J]. 地球物理学报, 2014(10):3453-3465.