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

张家成 张 宫* 黄若坤 覃莹瑶

1(长江大学油气资源与勘探技术教育部重点实验室 湖北 武汉 430100)2(塔里木油田分公司勘探开发研究院 新疆 库尔勒 841000)

0 引 言

核磁共振(NMR)测井通过观测储层中含氢核流体的NMR信号来获得用于评价储层参数的信息[1-2]。自20世纪50年代被引入到石油测井领域以来，核磁共振技术在饱和度、渗透率计算、流体识别、孔隙度、孔隙结构评价等方面均得到了成功应用[3-5]。2002年Sun等[6]和Hurlimann等[7]在石油测井领域引入了核磁共振波谱学中的二维核磁共振概念，进一步拓展了核磁共振测井的使用范围，提高了储层测井解释评价的准确率[8]。核磁反演方法也从对单个回波串数据进行反演发展到同时对多个回波串数据进行反演[9]，但大量的数据加重了反演过程中的运算负担。一直以来，国内外学者为满足反演精度和反演速度上的要求，对二维核磁反演方法进行了不断的探索。顾兆斌等[10]利用改进的奇异值分解法对二维核磁共振进行反演；谢然红等[9]针对二维的多回波串数据处理问题提出了多回波串联合反演方法；李鹏举等[11]针对BRD和TSVD两种反演方法的不足提出了变参量迭代的快速核磁反演方法。在二维核磁反演方法研究过程中，繁琐的数据构造和结果评测过程增加了研究的时间成本。

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

1 原 理

1.1 核磁共振弛豫测量

核磁共振测井通过观测流体中氢原子核的NMR信号，处理得到地层孔隙中流体的含量及其性质。岩石物理信息，如渗透率、孔径分布、孔隙度和束缚水，也可以通过核磁共振弛豫测量得到。核磁共振弛豫过程有两种，其中纵向弛豫计算如下：

(1)

式中：Pi是第i个弛豫信号占总的纵向弛豫信号的百分比；t表示时间；T1表示纵向弛豫时间，反映了自旋系统从外界吸收或者释放能量的效率，主要通过反转恢复脉冲序列或饱和恢复法脉冲序列进行测量[13]。

横向弛豫计算如下：

(2)

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

1.2 二维核磁弛豫构造

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

(3)

式中：j=1,2,…,m，m是构造谱的布点数目；Aj是第j组分的幅度；F是孔隙度刻度因子；T2j是第j组分的T2弛豫时间；T2mid是T2弛豫时间峰值的位置；Dj是第j组分的扩散系数；Dmid是扩散系数峰值的位置；σ1、σ2可以控制分布的展布宽度；ρ是T2j和Dj的相关系数。

为了更加真实地模拟出实际二维谱，二维构造谱必须能够反映多种不同性质流体的叠加效果，所以式(3)需要进一步改进为：

(4)

式中：j=1,2,…,m，m是构造谱的布点数目；设具有4类不同性质的流体，Fg是第g种流体的孔隙度分量；HIg是第g种流体的含氢指数；T2g,mid是第g种流体的T2弛豫中心值；Dg,mid是第g种流体的扩散系数峰值；σ1g、σ2g分别是第g种流体T2峰值的中心展布宽度和扩散系数峰值的展布宽度；ρ是T2j和Dj的相关系数。

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

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

图1 扩散编辑T2-D测量序列

在回波采集过程中，回波的幅度衰减计算如下：

(5)

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

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

(6)

式中：bik是回波间隔为TEk时第i个回波的幅度[9]。

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

1.4 二维谱反演方法

目前常见的二维谱反演方法主要有奇异值分解法(SVD)、非负约束正则化法(BRD)和截断奇异值分解法(TSVD)等。

SVD法是通过反演测量误差的核函数矩阵中的最小奇异值，以实现正则化，求取合理解。顾兆斌等[10]针对SVD法容易导致谱图不连续的问题，提出一种适用于对信噪比较高的数据进行反演的奇异值分解法。谭茂金等[15]提出了同时考虑TSVD的分辨率和精度的混合反演算法。

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

2 软件设计与实现

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

2.1 结构设计

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

图2 软件设计结构

软件设计时采用模块化编程方式，将所有数学算法函数集中打包放在统一的算法库中，从而独立于交互界面，便于算法的持续改进和扩展；应用模块包括二维谱构造模块(用于根据用户设置的流体类型及参数，构造出满足二维高斯分布的构造谱)，模拟回波采集模块(用于根据用户设置的采集模式信息，模拟仪器进行回波采集的过程，同时保存采集的回波串数据)，二维谱反演模块(实现了SVD和BRD二维谱反演算法，用于将采集得到的回波串数据反演成二维谱)；数据管理模块和绘图模块同样相对独立，并设置有外部平台接口，用于导入实验室岩心分析数据或实际测井数据。

主界面分为菜单栏、绘图区、参数区三部分，如图3所示。菜单栏有文件(NMR File)、参数设置(Parameters Set)和关于(About)。其中：NMR File菜单主要功能包括新建模拟文件、保存模拟文件、打开模拟文件和退出；Parameters Set菜单主要功能包括加载各个应用模块的参数设置界面、设置参数和导入其他平台的数据；About菜单用于显示本文软件的开发过程及人员信息。

图3 软件主界面

2.2 功能实现

(1) 二维谱构造模块。根据回波串组合序列完成“二维谱构造”模块，该模块是基于二维高斯分布进行设计的。通过选择相应流体类型、填写流体T2值及其展布、分量孔隙度(各类流体含量)、含氢指数和另一维度参数(纵向弛豫时间T1、扩散系数D、或内部梯度G)，以及最后的布点方案：T2及另一维度参数的布点的起始值、终止值和布点数目，来进行二维谱的构造。

对话框中的各个参数用户可根据实际模拟的流体进行设置，如图4所示，完成后软件会根据用户填写的结果自动计算出总孔隙度和核磁孔隙度。参数设置分为流体参数(各类流体的核磁属性信息及孔隙度分量)和全局参数(模拟模式构造谱分布参数)，所有参数设置完成后，可以实时地显示出二维构造谱，如图5所示。

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

图5 构造的T2-T1二维分布谱

(2) 模拟回波串采集模块。该模块根据上一步设置构造的二维谱数据，结合给定的模拟采集指令模拟采集回波串，并存储回波数据。参数包括：需要采集的回波组数，各组回波串对应的回波间隔(TE)值、等待时间(TW)值、回波采集个数(NE)、磁场梯度，以及设置的信噪比(S/N)。

参数设置完成后，软件会自动显示出模拟采集得到的回波数据及对应的噪音数据，用户可以根据图形实时对参数进行修订，如图6所示，模拟多TE二维采集得到的回波串数据如图7所示。另外，该模块允许用户导入实验室岩心分析得到的二维回波数据，以及实际测井得到的二维回波数据，并对数据进行后续的反演处理和分析。

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

图7 回波采集结果

(3) 二维谱反演模块。该模块的主要功能是将上一步模拟采集得到的回波串数据进行二维反演，得到最终的二维反演谱。需要设置的参数包括：T2维度布点参数(布点的数量，起始值和终止值)、第二维参数谱布点的方式(布点的数量，起始值和终止值)、反演使用的回波串选择、平滑因子以及反演方法的选择，参数设置对话框如图8所示。图9为反演得到的T2-T1二维谱。

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

除了软件内置的二维反演方法(SVD及BRD)外，还可以利用软件中提供的C#接口(该接口命名为：INMR_InversionMethod)来编写自己的反演算法，从而研究不同二维谱反演算法的适用性及影响因素。反演完成后，软件会再次将反演得到的二维谱进行正演(拟合)，观测拟合结果与采集得到的回波串是否吻合，从而判断反演效果好坏。二维回波串拟合效果如图10所示。

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

3 可靠性验证及应用实例

3.1 可靠性验证

软件编写完成后，以硫酸铜溶液的核磁物理实验为例，对正反演模拟软件的可靠性进行验证，具体实施步骤如下：

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

第二步分别取三种不同浓度的硫酸铜溶液3 ml放入玻璃质的集气瓶中，然后将三只集气瓶封口后放入试管内进行二维核磁共振(T2-T1)测量。实验结果显示：三种浓度的硫酸铜溶液的弛豫时间分别约为2 500、80、2 ms，如图11所示。

图11 二维反演谱(物理实验)

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

图12 二维构造谱(软件模拟)

图13 二维反演谱(软件模拟)

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

3.2 应用实例

研发的二维核磁共振正反演模拟软件可以进行多种二维模式的正反演数值模拟。现以T2-T1二维模式下核磁共振方法在含有束缚水的油层中识别流体的效果为例进行模拟，验证软件的可靠性。

油水的特征参数如表1所示，根据该表参数构造出油水模型的理想分布谱即构造谱，如图14所示，然后设置不同的等待时间模拟出15组回波串数据：TE均设为1 ms，TW分别为10、30、50、70、100、300、500、700、1 000、2 000、4 000、8 000、10 000、12 000、15 000 ms。

表1 油水模型参数表[14]

图14 二维构造谱

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

图15 二维反演谱

4 结 语

本文研发的二维核磁共振测井正反演模拟软件可以进行多种二维模式的正反演数值模拟，其模拟结果与实际岩心核磁共振仪器实验结果非常接近。该软件可以同时模拟4种不同的核磁共振弛豫过程，并且能够对各类核磁共振属性进行单参量因素分析，为二维核磁共振相关研究工作提供了高效的技术手段。