一种基于核磁共振成像的岩石非均质性定量评价方法

摘要

本发明公开了一种基于核磁共振成像的岩石非均质性定量评价方法，通过选层脉冲、相位编码脉冲和频率编码脉冲实现岩石三维空间定位，应用自旋回波序列得到成像信号，通过实验刻度进行成像实验参数的优化选择。在此基础上，对实验测量所得核磁成像信号进行数字图像处理生成伪彩图，通过标样的孔隙度与核磁成像信号强度的关系即可得到单层面的总孔隙度和孔隙度分布谱，将多层成像结果进行对比并定义孔隙度非均质性系数，即可得到岩石纵向的孔隙度分布特性及非均质性。此外，应用一阶球状变差函数模型及网格化搜索方法得到变差函数的特征参数，定义非均质性系数和相对非均质性系数实现岩石纵、横向的非均质性定量表征。

说明书

技术领域

本发明属于岩石物理实验及石油勘探领域，具体地，涉及一种利用核磁共振成像数据进 行岩石纵、横向非均质性定量评价的方法。

背景技术

随着国际能源需求的进一步进展，油气勘探正由常规向低孔低渗、致密储层发展。低孔 低渗、致密储层物性差、孔隙结构复杂、非均质性强导致测井响应特征复杂，油气评价难度 大，极大制约了油气钻探的成功率。

非均质性评价是低孔低渗、致密储层岩石物理实验和测井评价的重要内容。目前国内外 对岩石非均质的评价方法主要有岩心观测法、测井识别法、X-CT扫描法、铸体薄片法、扫描 电镜法等，分别对应不同的尺度。X-CT扫描、铸体薄片、扫描电镜是实验室研究岩石微观尺 度非均质性的常用方法，但由于它们只能对小尺度的样品进行研究分析，难以实现与地层岩 石的配套分析；岩心观察法对研究人员的经验性要求较强，且大都以定性描述为主，难以进 行定量的评价；测井识别法主要通过成像测井、偶极声波测井等新技术方法显示井下岩石的 结构特征，并通过图像处理技术进行岩石非均质性定量表征，受探测深度、仪器响应和测井 环境的影响。此外，测井非均质性定量表征方法只能实现岩石某一层面的非均质性特征。

发明内容

针对以上问题，本发明提供一种基于核磁共振成像的的岩石非均质性定量评价方法。

其技术解决方案是：

一种基于核磁共振成像的岩石非均质性定量评价方法，其特征在于包括以下步骤：

1.1岩心核磁共振成像测量方法和参数优化

1.1.1测量序列的优化

采用自旋回波序列作为基本序列，包括选层脉冲、相位编码脉冲和频率编码脉冲三部分；

1.1.2测量参数的优化

（1）标样选择：介质为水，孔隙度分别为1%、5%、8%、10%、15%和20%；

（2）对标样进行不同回波间隔TE下的一维核磁共振测试得到核磁共振T2谱；

（3）计算每块标样的核磁共振T2调和平均值；

（4）建立不同回波间隔TE下的核磁共振T2调和平均值与标样孔隙度的关系，选择相关 系数最高的组所对应回波间隔TE为最佳回波间隔TE；

（5）选择一系列恢复时间TR进行核磁共振成像测试，并得到不同恢复时间TR下的核磁 共振成像总信号与标样孔隙度的关系；

（6）选择核磁共振成像总信号与标样孔隙度信号线性最强，相关性最好所对应组的TR 作为核磁共振成像实验的恢复时间；

（7）针对待测样品，同样先开展一维核磁共振实验测试得到T2谱，重复步骤（3）～（4） 得到样品最佳回波间隔TE；

（8）根据步骤（6）所确定的恢复时间TR进行岩心的核磁共振成像实验；

1.2核磁共振成像的图像处理和非均质性定量表征

1.2.1标准样品的测试

选择一组标样：介质为水，孔隙度分别为1%、5%、8%、10%、15%、20%；对所选标样进 行核磁共振成像测试，采集信号的实部和虚部，设(i,j)点的实信号为Real(i,j)、虚信号为 Imaginary(i,j)，则该点的信号强度$\mathrm{Amplitude}(i,j)=\sqrt{{\mathrm{Real}}^2(i,j)+{\mathrm{Imaginary}}^2(i,j)});$通过线 性回归得到像点信号强度与孔隙度的关系；

1.2.2核磁共振图像的生成

核磁共振成像同时采集信号的实部和虚部；设(i,j)点的实信号为Real(i,j)、虚信号为 Imaginary(i,j)，则该点信号可表示为：Datacomplex(i,j)=Real(i,j)+iImaginary(i,j)，对 Datacomplex(i,j)进行二维傅里叶变换，可得核磁共振成像图；

1.2.3孔隙度分布

由于信号强度直接与孔隙度相关，根据标样测试结果的信号强度与孔隙度的刻度关系， 可得每个像素点所表征的孔隙度，进而得到该层面的孔隙度分布，而总信号所对应的孔隙度， 则为该成像面的孔隙度；

1.2.4非均质性定量表征

应用选层梯度控制层面，沿着岩心轴向进行多层采集，得到岩石轴向上的孔隙度分布及 总孔隙变化；定义孔隙度非均质性系数为φ_{heterogeneity}(i)=φ(i)/min(φ)，min(φ)为轴向最小孔隙 度；非均质性系数越大，则岩心非均质性越强；

或基于一阶球状变差函数拟合的非均质性定量表征：应用一阶球状模型进行拟合分别得 到变差函数的特征参数，将归一化的实验变差函数写为：

r(h)=[S₂(0)-S₂(h)]/S₂(0)

其中S₂(0)为两点相关函数矩阵的方差；S₂(h)为两点相关函数；

将一阶球状模型函数写为：

$r\left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(\frac{3h}{2a}-\frac{h^3}{2a^3})& 0<h\le a\\ C_0+C& h>a\end{array}\right)$

其中a为变程，表征变量在其邻域内的最大影响距离；C₀为块金常数，是基于滞后尺度下的 变异性大小的标志；C为拱高；C₀+C称为基台值，是表征变异性大小的极限值；首先将a进 行网格化，然后进行不同变程下的最小二乘拟合并对比不同变程下的拟合效果，确定最优参 数；

变程与非均质性成反比，基台值与非均质性成正比，因此定义非均质性系数I为：

$I=\frac{C_0+C}{a}$

从上式可知，非均质性系数越大，孔隙的非均质性越强；如果轴向最小非均质性系数为 min(I)，则可定义相对非均质性系数为I_{heterogeneity}(i)=I(i)/min(I)，相对均质性系数越大，则 岩心非均质性越强。

本发明具有以下效果：

1、本发明提供的基于核磁共振成像的非均质性定量表征方法，可以得到岩心尺度的孔隙 结构特征和非均质性特征，有助于低孔低渗、致密等非常规储层的孔隙结构分析及非均质性 定量研究，比铸体薄片、扫描电镜、X-CT扫描等方法的尺度要大，可以克服非均质岩心取样 无代表性的困难。样品预处理及成像分析等过程简单，计算方便，实用性强。

2、利用本发明得到的岩石非均质性系数，将有助于定量研究岩石非均质性对声波、电阻 率、毛管压力等属性的影响，结合“岩心刻度测井”，能实现地层条件下的岩石非均质性定量 表征，为低孔低渗、致密等非常规储层的孔隙结构评价和非均质性研究带来新的思路。

附图说明

图1是本发明提供的一种基于核磁共振成像的岩石非均质性定量评价方法的技术流程 图；

图2是核磁共振成像实验中选层梯度对岩心进行切片分析的示意图；

图3是某砂岩进行10层切片得到的核磁共振T2加权成像伪彩图；

图4是图3所述砂岩第一层切面的两点相关函数图；

图5是图3所述砂岩第一层切面的变差函数及一阶球状模型拟合效果图。

图6是图3所述砂岩岩心第一层切面的孔隙度分布图。

图7是图3所述砂岩岩心10层切面的孔隙度非均质系数分布图。

图8是图3所述砂岩岩心10层切面的相对非均质性系数分布图。

具体实施方式

一种基于核磁共振成像的岩石非均质性定量评价方法，其包括以下步骤：

1.1岩心核磁共振成像测量方法和参数优化

1.1.1测量序列的优化

核磁共振测井和一维核磁共振一般测量的是横向驰豫信号（T2），而核磁共振成像采集的 是与质子密度、横向驰豫时间（T2）、纵向驰豫时间（T1）有关的信号。为了与常规核磁共振 进行对比，需要通过相关序列，突出横向驰豫信号。因此本发明采用自旋回波序列作为基本 序列，其包括选层脉冲、相位编码脉冲和频率编码脉冲三部分。

1.1.2测量参数的优化

通过改变测量序列的关键参数如恢复时间、回波间隔等即可实现不同条件下的成像信号。 当恢复时间T_R很长，回波间隔T_E约等于T₂时，像点信号强度近似等于此时像中增加 了横向弛豫时间的权重，称为T₂加权成像。

参数优化方案步骤如下：

（1）标样选择：介质为水，孔隙度分别为1%、5%、8%、10%、15%和20%；

（2）对标样进行不同回波间隔TE下的一维核磁共振测试得到核磁共振T2谱；

（3）计算每块标样的核磁共振T2调和平均值；

（4）建立不同回波间隔TE下的核磁共振T2调和平均值与标样孔隙度的关系，选择相关 系数最高的组所对应回波间隔TE为最佳回波间隔TE；

（5）选择一系列恢复时间TR进行核磁共振成像测试，并得到不同恢复时间TR下的核磁 共振成像总信号与标样孔隙度的关系；

（6）选择核磁共振成像总信号与标样孔隙度信号线性最强，相关性最好所对应组的TR 作为核磁共振成像实验的恢复时间；

（7）针对待测样品，同样先开展一维核磁共振实验测试得到T2谱，重复步骤（3）～（4） 得到样品最佳回波间隔TE；

（8）根据步骤（6）所确定的恢复时间TR进行岩心的核磁共振成像实验。

1.2核磁共振成像的图像处理和非均质性定量表征

1.2.1标准样品的测试

选择一组标样：介质为水，孔隙度分别为1%、5%、8%、10%、15%、20%；对所选标样进 行核磁共振成像测试，采集信号的实部和虚部，设(i,j)点的实信号为Real(i,j)、虚信号为 Imaginary(i,j)，则该点的信号强度$\mathrm{Amplitude}(i,j)=\sqrt{{\mathrm{Real}}^2(i,j)+{\mathrm{Imaginary}}^2(i,j)});$通过线 性回归得到像点信号强度与孔隙度的关系。

1.2.2核磁共振图像的生成

核磁共振成像同时采集信号的实部和虚部；设(i,j)点的实信号为Real(i,j)、虚信号为 Imaginary(i,j)，则该点信号可表示为：Datacomplex(i,j)=Real(i,j)+iImaginary(i,j)，对 Datacomplex(i,j)进行二维傅里叶变换，可得核磁共振成像图。

1.2.3孔隙度分布

由于信号强度直接与孔隙度相关，根据标样测试结果的信号强度与孔隙度的刻度关系， 可得每个像素点所表征的孔隙度，进而得到该层面的孔隙度分布，而总信号所对应的孔隙度， 则为该成像面的孔隙度。

1.2.4非均质性定量表征

方法一：应用选层梯度控制层面，沿着岩心轴向进行多层采集，得到岩石轴向上的孔隙 度分布及总孔隙变化。如果轴向最小孔隙度为min(φ)，则可定义孔隙度非均质性系数为 φ_{heterogeneity}(i)=φ(i)/min(φ)，非均质性系数越大，则岩心非均质性越强。

方法二：基于一阶球状变差函数拟合的非均质性定量表征。在地质统计学中，通常用变 差函数进行非均质性定量表征。本发明在研究变差函数性质的基础上，通过相关函数理论得 到图像的实验变差函数，应用一阶球状模型进行拟合分别得到变差函数的特征参数。将归一 化的实验变差函数写为：

r(h)=[S₂(0)-S₂(h)]/S₂(0)

其中S₂(0)为两点相关函数矩阵的方差；S₂(h)为两点相关函数；

对于离散的二维数据，其两点相关函数可写为：

$S_2(x,y)=\underset{i=1}{\overset{M-x}{\Sigma }}\underset{j=1}{\overset{N-y}{\Sigma }}\frac{f(i,j)f(i+x,j+y)}{(M-x)(N-y)}$

其中：M、N为图像在x、y方向的像素点数。为提高运算速度，应用傅立叶变换将图像 从空域转换成频域后进行处理。图像f(x,y)的二维傅立叶变换为：

Z_f(δ_x,δ_y)=FFT₂[f(x,y)]

式中：δ_x、δ_y为空域（x,y）的频域表示。

图像在频域的相关函数矩阵可写为：

$C({\delta }_x,{\delta }_y)=Z_f({\delta }_x,{\delta }_y)Z_f^{*}({\delta }_x,{\delta }_y)$

式中：为Z_f(δ_x,δ_y)的复共轭。

空域的相关函数矩阵是频域相关函数矩阵的反傅立叶变换：

$S(x,y)=\frac{{\mathrm{IFFT}}_2\left[C({\delta }_x,{\delta }_y)\right]}{(M-x)(N-y)}$

将一阶球状模型函数写为：

$r\left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(\frac{3h}{2a}-\frac{h^3}{2a^3})& 0<h\le a\\ C_0+C& h>a\end{array}\right)$

其中a为变程，表征变量在其邻域内的最大影响距离；C₀为块金常数，是基于滞后尺度 下的变异性大小的标志；C为拱高；C₀+C称为基台值，是表征变异性大小的极限值；首先 将a进行网格化，然后进行不同变程下的最小二乘拟合并对比不同变程下的拟合效果，确定 最优参数；

变程与非均质性成反比，基台值与非均质性成正比，因此定义非均质性系数I为：

$I=\frac{C_0+C}{a}$

从上式可知，非均质性系数越大，孔隙的非均质性越强；如果轴向最小非均质性系数为 min(I)，则可定义相对非均质性系数为I_{heterogeneity}(i)=I(i)/min(I)，相对均质性系数越大，则 岩心非均质性越强。

通过核磁共振成像测量参数优化和数字图像处理、非均质性定量表征，可得到岩石在纵、 横向的孔隙变化特征及非均质性特征，实现岩石的非均质性定量表征。

下面结合附图对本发明的具体应用实例进行说明。

一种基于核磁共振成像的岩石非均质性定量评价方法，通过选层脉冲、梯度编码和相位 编码脉冲实现岩石空间定位，应用自旋回波序列得到成像信号，通过实验刻度进行成像实验 参数的优化选择。在此基础上，对实验测量所得核磁成像信号进行数字图像处理生成伪彩图， 通过标样的孔隙度与核磁成像信号强度的关系即可得到单层面的总孔隙度和孔隙度分布谱， 将多层成像结果进行对比并定义孔隙度非均质性系数，即可得到岩石纵向的孔隙度分布特性 及非均质性。此外，应用一阶球状变差函数模型及网格化搜索方法得到变差函数的特征参数， 定义非均质性系数和相对非均质性系数实现岩石纵、横向的非均质性定量表征。一般而言， 孔隙度非均质性系数越大，岩石横向非均质性程度越高；相对孔隙度非均质系数越大，岩石 纵向非均质性程度越高。

图1是基于核磁共振成像的岩石非均质性定量评价方法流程图，主要包括核磁共振成像 信号产生、核磁共振信号检测与编码、核磁共振信号采集与存储、核磁成像图像显示与处理、 核磁成像图像非均质性分析五部分，这五部分缺一不可，且顺序不可颠倒。

图2是核磁共振成像实验中选层梯度对岩心进行切片分析的示意图。在选层脉冲的作用 下进行岩石切面和层面选择。选择好层面后，通过相位编码脉冲和频率编码脉冲实现单层面 的空间定位，将定位点称为像素（voxel)，像素点坐标与岩石层面的物理地址呈一一对应关 系。如图2所示，我们沿着岩石纵向（层理方向）进行选层，本示例中共选择10个层面。

图3是某砂岩进行10层切片得到的核磁共振T2加权成像伪彩图。从(a)到（j）分别表 示10个切面。T2加权成像伪彩图的成像采用傅里叶变换得到。为了实现海量数据的快速分 析，我们在这里采用快速傅里叶变换方法。值得注意的是，由于核磁共振成像分别采集信号 的实部和虚部，在快速傅里叶变换时我们首先需要将信号进行合成。设(i,j)点的实信号为 Real(i,j)、虚信号为Imaginary(i,j)，则该点信号可表示为： Datacomplex(i,j)=Real(i,j)+iImaginary(i,j)。

图4是核磁共振T2加权成像图的两点相关函数，为了实现实验变差函数的拟合，两点相 关函数的计算是至关重要的。对于离散的二维数据，其两点相关函数可写为：

$S_2(x,y)=\underset{i=1}{\overset{M-x}{\Sigma }}\underset{j=1}{\overset{N-y}{\Sigma }}\frac{f(i,j)f(i+x,j+y)}{(M-x)(N-y)}$

其中：M、N为图像在x、y方向的像素点数。为提高运算速度，应用傅立叶变换将图像从空 域转换成频域后进行处理。图像f(x,y)的二维傅立叶变换为：

Z_f(δ_x,δ_y)=FFT₂[f(x,y)]

式中：δ_x、δ_y为空域（x,y）的频域表示。

图像在频域的相关函数矩阵可写为：

$C({\delta }_x,{\delta }_y)=Z_f({\delta }_x,{\delta }_y)Z_f^{*}({\delta }_x,{\delta }_y)$

式中：为Z_f(δ_x,δ_y)的复共轭。

空域的相关函数矩阵是频域相关函数矩阵的反傅立叶变换：

$S(x,y)=\frac{{\mathrm{IFFT}}_2\left[C({\delta }_x,{\delta }_y)\right]}{(M-x)(N-y)}$

图5是图3所述砂岩第一层切面的变差函数及一阶球状模型拟合效果图。从图中分析可 知，本发明采用的拟合算法的计算结果精确度高，拟合值与实际值误差较小。

图6是图3所述砂岩岩心第一层切面的孔隙度分布图。从图中可知，该岩石在第一层切 面上的核磁共振孔隙度呈单峰分布特征，主峰分布在0.0004%之间。

图7是图3所述砂岩岩心10层切面的孔隙度非均质系数分布图。10个切面的计算结果 表明：该岩心孔隙度分布在6.42%～7.86%之间，平均为7.11%。相对孔隙度在纵向上存在较 小幅度的变化。

图8是图3所述砂岩岩心10层切面的相对非均质性系数分布图。该岩心的非均质性系数 分布在0.0024～0.0054之间，平均为0.0041，分析可知，岩心横向非均质程度较弱。但相 对非均质系数在纵向上存在较大幅度的变化。

通过图3～图8的实例分析可知，该样品的整体非均质性较弱，但在纵向上仍存在一定 的非均质性。

需要说明的是，上述实施例仅仅是本发明一部分实施例，而不是全部实施例。基于本发 明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例都 属于本发明的保护范围。