一种通用的岩心核磁共振实验数据分析软件设计与实现

Design and Implementation of a Universal Core NMR Data Analysis Software

**【作者：覃莹瑶，张宫*，何宗斌，张家成】**

摘要

岩心核磁共振实验中，原始采集数据的处理与分析非常关键，但各厂家核磁共振仪器提供的数据处理配套软件的反演方法、参数均不相同，使得同样的岩心测量得到的回波信息无法进行对比分析。研发了一套通用的岩心核磁共振实验数据分析软件，基于统一的数据存储模式、统一的快速高分辨反演算法，允许设置不同反演参数和岩心分析参数，进行岩心核磁数据的反演和参数计算，实现了所有仪器测量回波数据的统一处理和分析。目前软件支持NIUMAG公司、OXFORD公司、Magritek公司和CoreSpec-1000型仪器测量得到的原始数据，处理结果精确，具有操作简单、方便扩展的特点。

0 引言

核磁共振测井是一种有效的测井方法，岩心核磁共振实验分析是NMR测井资料采集、处理、解释及应用的基础^[1]。测井前做少量的实验室测量能节省实际测井时间，和提高测井曲线质量，而且实验室的岩心核磁共振测量实验能够实现岩心与测井之间的标定,建立核磁共振测井解释模型^[2]。利用岩心核磁共振实验可以直接观测到岩石孔隙的流体信号，得到与岩性无关的孔隙度、孔径分布、渗透率、自由流体指数及流体饱和度等参数^[3,4]。岩心核磁共振实验中，原始采集数据的处理与分析非常关键，但各厂家核磁共振仪器提供的数据处理配套软件的反演方法、参数均不相同。例如美国NUMAR公司生产的CoreSpec-1000型岩心核磁共振分析仪由Echofit软件利用非负最小二乘法（NNLS）和基于奇异值分解的MAP-Ⅱ算法进行反演^[5]；英国Resonance Instruments公司的MARAN系列核磁共振岩心分析仪没有提供配套的数据处理分析软件^[6]；NIUMAG公司生产的岩心核磁共振分析仪配套软件使用罚函数法（BRD）和联合重建迭代法（SIRT）进行反演^[7]。仪器性能和反演算法都不相同导致同样的岩心测量得到的回波信息无法进行对比分析，无法确定反演的T₂谱有差异的原因，很大程度上限制了核磁岩心分析仪的作用。王忠东等研发的岩心核磁实验解释分析软件CoreMR可以对MARAN系列谱仪所采集岩心和流体核磁信号进行的解释分析^[6]，但是这款软件只能处理一种谱仪采集的数据，没有实现对不同厂家仪器测量数据的统一处理分析。
研发一种通用的岩心核磁共振实验数据分析软件可以解决这一问题，使用C#语言进行软件开发，支持硬件加速的WPF通用框架进行绘图，采用格式统一、易于交互的XML序列化存储反演参数，易于传输、解析的JSON序列化存储采样数据和处理结果。通过采用统一的模式进行数据存储、统一的快速高分辨反演算法，设置不同参数进行核磁参数计算，本软件可以使所有仪器测量的回波数据得到统一处理和分析，在消除软件处理方法差异的影响上使同样岩心的回波信息进行对比分析。

1 原理

岩心核磁共振实验数据分析软件一般包括加载数据、反演、计算参数和导出报告四个步骤。本文设计的软件力图将不同仪器所采集的数据进行统一处理，业务流程如图1所示，具体步骤为：①数据加载要实现把所有仪器所测的不同格式数据解析成统一的模式；②进行T₂谱反演；③计算参数要根据反演后的T₂谱计算出孔隙度、束缚水饱和度等核磁参数；④数据的处理和分析完成后需要导出指定格式的实验报告。

图1 软件业务流程设计

1.1 数据解析

软件需要加载各类岩心分析仪测量得到不同格式的原始回波串数据，首先要对数据进行解析。以CoreSpec-1000型仪器测量得到的RES数据和NIUMAG公司仪器测量得到的PEA数据为例。首先进行RES文件或PEA、PAR文件的格式分析，如图2所示RES文件中有记录参数信息的头文件和原始回波数据；PEA文件中是原始回波数据，参数信息在PAR文件中。然后对文件进行解析，读取RES文件或PEA、PAR文件的数据后，利用C#面向对象程序设计“CPMGCoreData”类提供的一系列方法将数据都解析成统一的模式。

图2 RES文件数据

1.2 反演方法

对原始回波信号进行反演处理得到分布谱，这个过程叫做反演^[8]。 T₂谱的反演方法较多。其中，最常用的算法有：罚函数法（BRD）、非负最小二乘法（NNLS）、联合迭代重建算法（SIRT）、奇异值分解法(SVD)以及针对SVD法的各种改进型反演算法。
软件的反演模块主要实现了罚函数法（BRD）和联合迭代重建算法（SIRT）改进的快速高分辨反演算法，可以在不抽样压缩回波串的情况下，实时得到T₂谱。
在岩心核磁共振实验中，核磁共振横向弛豫信号y(t)描述为如下的形式：

\[y\left(t\right)=\ \sum_{i}{P_iexp\left(-\frac{t}{T_{2i}}\right)}\]

其矩阵方程的形式为

\[Y=A\bullet\ P (1)\]

Y为测量的横向弛豫信号；P为所要反演计算的T₂谱中各个弛豫时间所对应的谱幅度值。在此利用联合迭代重建算法（SIRT）实现反演^[9]，首先给定谱的初始模型P’，求出预测弛豫信号Y与实测弛豫信号Y’的误差ΔY，即

\[∆Y=Y-\ Y^{‘} (2)\]

\[\sum_{(j=1)}^m\ a_{ij}\ ∆p_j=∆y_i (3)\]

利用ΔY_i求出弛豫谱幅度的修正量Δp_j，这是SIRT算法实现的主要思想。p_j的修正量为

\[{∆p_{j}}^{({q + 1})} = \frac{1}{IA\left( j \right)}{\sum\limits_{i = 1}^{n}\left\lbrack {\lambda a_{ij}{∆y_{j}}^{(q)}} \right\rbrack} (4)\]

式中，IA(j)为矩阵A第j列非零元素的个数；Δy_j^(q)为第q次迭代中第j个弛豫信号的残差。
通过式(4)计算的修正公式为

\[P^{\left(q+1\right)}=P^{\left(q\right)}+∆P (5)\]

SIRT算法简单容易实现，在利用算法处理时不需要用户干预，也不需要设置一些复杂的反演控制参数，减少了人为因素方面的反演结果误差。

1.3 参数计算

（1）核磁孔隙度的确定
用标准样品对饱和岩样测得的T₂谱进行刻度，将核磁信号强度转换成孔隙度^[10]，转换公式如下：

\[\phi_{nmr}=\ \sum_{i}\frac{m_i}{M}\cdot\frac{S}{s}\cdot\frac{G}{g}\cdot\frac{V}{v} (6)\]

式中：Φ_nmr为岩样核磁孔隙度值（%）；M为标准样品T₂谱的总幅度；V为标准样品总含水量（cm³）；S为标准样品在核磁共振数据采集时的累积次数；G为标准样品在核磁共振数据采集时的接受增益；m_i为第i个岩样T₂分量的核磁共振T₂谱幅度；v为岩样的体积（cm³）；s为岩样在核磁共振数据采集时的累积次数；g为样品在核磁共振数据采集时的接受增益。

（2）T₂截止值的确定
采用孔隙度累加法确定T₂截止值^[11]。先将岩样用水100％饱和得到饱水样，进行核磁共振测量得到T₂分布、累加孔隙度曲线和有效孔隙度值（用MPHI表示）。然后对岩样作脱水处理，在一定压力的条件下，使自由水脱出岩样，得到孔隙空间中只剩下束缚水的离心样，再做岩心核磁共振测量，测得T₂分布、累加孔隙度曲线和束缚水孔隙体积值（用MBVI表示）。根据测量结果，以MBVI作一条与横轴平行的平行线，此线与100％饱和水时的累加孔隙度曲线有一交点，以该交点作一条与纵轴平行的平行线，此线与横轴交点对应的T₂值即T_2cutoff值。

（3）T₂平均值的计算
在分析过程中，常用T₂的平均值来表征T₂分布，算术平均值T_2s和几何平均值T_2g按下式计算^[12]：

\[T_{2s}=\frac{\sum T_{2i}\phi_i}{\phi_{nmr}} (7)\]

\[T_{\mathrm{2g}}=\left(\prod T_{2i}^{\phi_i}\right)^\frac{1}{\phi_{nmr}} (8)\]

式中：Φ_nmr为核磁总孔隙度；Φ_i为对应分量T_2i的孔隙度分量。

（4）核磁束缚水饱和度的确定
确定核磁共振束缚水饱和度的方法是以小孔隙束缚水体积模型作为理论基础的T₂截止值法^[13]。在岩心核磁共振实验测得的T₂谱曲线中，核磁束缚水饱和度为T₂谱中小于T₂截止值T_2cutoff的不可动峰下包面积和整个T₂谱曲线下包面积的比值。束缚水体积等于核磁束缚水饱和度和孔隙度的乘积，可动水体积等于岩样的孔隙体积与束缚水体积的差值。

1.4 数据报告导出

利用软件对岩心原始回波串数据设置合适参数进行处理和分析后，可以将文件保存成易于传输和解析的JSON格式文件，还能够将分析结果导出成所选择的固定格式的EXCLE格式报告，如图3所示。岩心测量报告中有具体的实验结果和回波数据，实验结果中包括了岩心的井号、样号、体积、回波间隔等参数信息，孔隙度、饱和度等核磁参数信息，回波衰减曲线和反演的T₂分布谱；回波数据中包括了采集时间、原始回波和T₂分布的时间及幅度。

图3 岩心数据报告

2 软件设计及实现

根据原理和软件业务流程设计，本软件基于.NET Framework 4.6开发框架在Visual Studio2017上使用C#语言进行开发。

2.1 软件架构

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

图4 软件设计结构

软件使用的所有算法统一放在算法库中，独立于交互界面。本文软件实现了罚函数法和联合迭代重建算法改进的反演算法，利用外部算法接口，可以对算法进行持续改进和扩展；应用模块包括数据加载模块、反演模块和参数计算模块；绘图模块和数据管理模块同样相对独立，并专门设置有外部平台接口，可以很容易的将研究成果迁移到其他平台进行使用。其中最核心的三个模块设计了四个类：“CPMGCoreData”类实现了数据的输入与输出；“T2Inversion”类实现反演的功能；“Parameter”类实现参数计算，这三个类都引用了通用算法库的“CommonAnalysis”类。

图5 软件主界面

软件主界面如图5所示，包括菜单栏、绘图区、参数查阅区和状态栏四部分。菜单栏有文件（File）、工具（Tool）、数据分析(Data Analysis)和帮助（Help）。
其中File菜单主要是加载实验数据、输出报告、打开文件和保存文件；Tool菜单中有不同仪器的工具箱，功能主要有批量实验数据处理、导出测井曲线文件；Data Analysis菜单主要实现反演和参数计算功能；帮助菜单中是本软件的一些开发过程及人员信息。

2.2 功能实现

（1）岩心分析数据加载模块
本软件将各类岩心分析仪测量得到的原始回波串数据经过“CPMGCoreData”类提供的一系列方法将数据都解析后，用统一的模式进行存储。软件目前能实现的主要包括：1.CoreSpec-1000型岩心核磁共振分析仪测量得到的RES数据；2.NIUMAG公司生产的各型岩心核磁共振分析仪测量得到的PEA数据；3.OXFORD公司生产的各型岩心核磁共振分析仪测量得到的RIDAT数据；4.Magritek公司生产的仪器测量得到的数据；5.EXCEL通用数据格式。
软件自动解析并加载原始回波数据和对应的参数信息，如图6所示，通过参数表格的形式，显示出了岩心分析仪采集数据的参数及采集得到的具体回波串数据，并在绘图区域绘制原始回波串（图7）。参数表列中主要参数有：等待时间、回波间隔、采集数目等。

图6 数据加载模块界面

图7 原始回波串显示

（2）反演模块
反演模块利用罚函数法（BRD）和联合迭代重建算法（SIRT）改进的快速高分辨反演算法对加载的回波串数据进行反演。核心算法是原理中式(1)—(5)实现的联合迭代重建算法。反演参数设置的参数表区域中，可以直接在里面进行反演参数的修改。主要参数包括：是否需要进行基线漂移校正、反演平滑级别、起始反演回波、终止反演回波、是否进行边界约束、反演布点数目、布点最小值、布点最大值。在反演参数设置区域填写合适的反演参数后，单击右下角反演按钮，即可进行快速高分辨反演，反演结果（图8）会显示在左下角的绘图区域。反演的同时，会对原始回波串进行拟合观察反演结果的好坏。

图8 反演模块及反演结果

（2）参数计算模块
参数计算模块采用了参数计算原理中描述的确定孔隙度、T₂截止值等核磁参数的公式和方法进行计算。图9所示的参数计算设置区域中，主要需要设置的参数包括：样品体积，T₂截止值，渗透率计算模型参数等。反演完成后，在参数设置区域填写正确的岩心分析参数，然后单击分析按钮即可完成分析，其中最重要参数的是岩心的体积。参数计算能够得到T₂几何均值、有效孔隙度、可动孔隙度、渗透率、束缚水饱和度等参数。

图9 参数计算模块及计算结果

3 可靠性验证及应用实例

3.1 可靠性验证

为了验证软件的可靠性，将本文研发的软件和其他软件对回波信息处理的结果进行比对，以CoreSpec-1000型仪器和NIUMAG公司仪器的配套软件为例。
本软件和CoreSpec-1000型仪器配套的软件进行比对。用CoreSpec-1000型仪器测量硫酸铜标样，图10是CoreSpec-1000型仪器的配套软件处理的T₂谱和本软件对原始回波数据进行反演处理得到的T₂谱，可以看出峰值位置非常接近，都在200 ms附近。

图10 CoreSpec-1000配套软件（左）和本软件（右）反演的T2谱比对

本软件和NIUMAG公司生产的仪器配套的软件进行比对。用NIUMAG公司生产的仪器测量一块饱和水砂岩岩心，图11是NIUMAG仪器配套软件反演的T₂谱和用本软件进行反演处理得到的T₂谱。可以观察到两个T₂谱的形态相似，峰值位置也非常接近。
这两个验证说明本文研发的软件和其他软件处理结果一致，验证了本软件的可靠性。

图11 NIUMAG仪器配套软件（左）和本软件（右）反演的T2谱比对

3.2 应用实例

研发的岩心核磁共振实验数据分析软件可以使不同仪器测量同样岩心得到的回波信息进行对比分析。以对一块饱和水的砂岩岩心进行不同仪器的核磁共振比对实验为例，具体步骤为：
第一步：用CoreSpec-1000型仪器和NIUMAG仪器对一块饱和水的砂岩岩心进行核磁共振实验，得到两组岩心数据；
第二步：将两仪器的两组岩心数据导出，用研发的岩心核磁共振实验数据分析软件处理这两组数据，进行反演得到T₂谱；
图12是本文研发的软件对NIUMAG和CoreSpec-1000两种仪器测得同样岩心两组回波数据进行反演得到的T₂谱，在利用本软件对两组数据处理的反演方法和参数相同的情况下，观察到T₂谱的形态有差异，峰值位置不同。左图中主峰峰值位置在100 ms附近，右图主峰峰值位置在200 ms附近，是由于左图的仪器主频较高，右图的仪器主频较低，反映出同样岩心T₂谱的不同是仪器性能不同导致的影响。说明了本软件用统一的反演算法处理岩心数据，能对比同样岩心的回波信息反映出仪器性能的不同。

图12 软件对NIUMAG仪器（左）和CoreSpec-1000型仪器（右）

软件可以一次处理多回波间隔的数据得到长短T_E的核磁T₂谱。以CoreSpec-1000型仪器所测量的两个回波间隔分别为0.1 ms和0.2ms的数据为例，用本软件一次进行处理，得到了如图13所示的回波串和长短T_E的T₂谱。

图13 多回波间隔数据处理结果

4 结语

本文开发的软件处理岩心核磁共振实验数据结果精确，操作简单方便，且具有通用性，实现了对所有仪器测量的回波数据进行统一处理和分析，可以对比分析不同仪器测量同样岩心的回波信息，反映出仪器性能的差异。另外软件预留有外部接口，方便扩展，其中其他系统接口可以对接外部平台，外部算法接口可以迅速实现新的算法。

参考文献

王忠东,汪浩,李能根,等.核磁共振岩心基础实验分析[J].测井技术,2001,25(3):170-174.
Murphy D P. NMR logging and core analysis-simplified[J]. World Oil,1995, 216(4): 65-70.
Morriss C,Rossini D,Straley C,et al. Core analysis by low-field NMR[J]. The Log Analyst, 1997, 38(2): 84-94.
肖立志,石红兵.低场核磁共振岩心分析及其在测井解释中的应用[J].测井技术, 1998, 22(1): 42-49.
肖立志.核磁共振成像测井与岩石核磁共振及其应用[M].北京:科学出版社,1998.
王忠东,魏钢,柯启宇.岩心核磁实验解释分析软件的编制与应用[C]// 第八届核技术与测井技术年会, 2002.
纽迈岩心核磁共振分析测量软件V4.0用户手册[M].
Kevin Munn, Douglas M Smith. A NMR technique for the analysis of pore structure: numerical inversion of relaxation measurements[J]. Journal of Colloid and Interface Science, 1987, 19(1): 117-126.
姚绪刚,王忠东.一种新的核磁共振弛豫谱反演算法[J]. 测井技术,2003,27(5):373-376.
肖秋生,朱巨义.岩样核磁共振分析方法及其在油田勘探中的应用[J].石油实验地质,2009, 31(1): 97-100.
肖亮,肖忠祥.核磁共振测井T2cutoff确定方法及适用性分析[J]. 地球物理学进展, 2008,23(1): 167-172.
白松涛,程道解,万金彬,等.砂岩岩石核磁共振T2谱定量表征[J].石油学报,2016, 37(3):382-391.
Kleinberg R L,Boyd A. Tapered cutoffs for magnetic resonance bound water volume [C]. SPE Annual Technical Conference and Exhibition, 1997,37: 387-388.

版本信息：

编辑：覃莹瑶，2020.03.30