有机页岩高分辨率数字岩心构建与分析方法与流程

本发明涉及石油测井岩石物理研究领域，具体地说，涉及一种有机页岩高分辨率数字岩心构建与分析方法。

背景技术：

页岩油气是一种重要的非常规油气类型，针对有机页岩岩石物理特性研究是储层测井评价与地震预测的基础。岩石物理研究主要包括理论方法研究，岩石物理实验以及近年来发展起来的数字岩石物理分析技术。由于页岩孔隙小，矿物成分复杂，数字岩心分析及其数值模拟技术是一个有效的分析手段。因此，如何结合有机页岩的特点构建有效的接近实际的数字岩心模型则是微观数值模拟的基础与关键。

目前，数字岩心建模方法有两大类：物理实验方法和数值重建方法。物理实验方法借助扫描电镜和ct成像仪等高精度仪器获取岩心二维图像，经三维重建即可得到三维数字岩心(hazlettetal.,1995)。近年来，更高分辨率的fib-sem(focusedionbeamcombinedwithsem)图像扫描技术，其分辨率高达3nm，能够满足页岩的岩石物理实验。数值重建方法主要是基于岩心切片和重建算法建立数字岩心，如模拟退火法(王晨晨等，2012)、过程法(jin等，2011)、多点地质统计学以及一些混合方法，如过程法和模拟退火法(刘学峰等，2010)等。对于有机页岩大多采用fib-sem方法。

对于有机页岩，目前采用fib-sem成像技术比较昂贵，限制了该项技术的使用。而且，该技术不只针对页岩，实验的结果是针对具体的岩心，不能代表形形色色的岩石，不能反映不同润湿性的情况，也不能反映不同饱和度的情形，更不能模拟不同含气饱和度的情形。因此，目前的技术进展限制了后续微观数值模拟中影响因素和微观响应规律的研究。

技术实现要素：

本发明正是为了解决上述技术问题而设计的一种针对页岩油气的数字岩心模型的制作方法，能够用于模拟不同饱和度，不同润湿性的含气或含油页岩的微观数值模拟，能够为后续的电学、声学、核磁共振等物理特性的微观数值模拟提供基础模型。

本发明解决其技术问题所采用的技术方案是：

一种有机页岩高分辨率数字岩心构建与分析方法，有机页岩高分辨率数字岩心构建与分析方法按如下步骤实现：

a.页岩样品制备与同步辐射扫描

a1.借助高分辨率电子显微镜选取页岩样品，放入ct仪器载物台进行ct扫描；

a2.扫描图像重构。从扫描图像中剔除背景图像，得到x-y平面的单帧16位页岩岩心灰度图像；

a3.用图像处理软件把多帧16位图像构建为三维图像，观察岩心样本的整体结构特征和孔隙特征；

b.数字图像滤波。对步骤a2所得页岩岩心灰度图像进行图像滤波处理，使岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于ct扫描重建过程中产生的噪声点；

c.ct图像组分分割。对步骤b所得ct图像进行阈值分割，划分岩石的骨架和孔隙：

c1.阈值计算及骨架与孔隙的准确分割

采用阈值分割法，其原理是：岩心ct图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值，通过阈值选取技术选取合理的阈值，将ct图像中的骨架和孔隙准确地分割开来；

c2.骨架矿物的精细分割

对于一般岩石ct图像，只要准确地分割ct图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩ct图像来说，由于其成分复杂，还需要对页岩ct图像的矿物成分进行进一步分割，其主要步骤为：

第一、导入ct图像数据，设置数据格式，将三维ct图像裁剪为需要的尺寸；

第二、依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；

第三、使用灰度阈值分割程序，精细地划分其矿物成分；

第四、导出数据并构建数字岩心，即岩心三维成像图。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤a中同步辐射扫描是指对页岩岩心的高分辨率扫描，其高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm，大视场模式成像视场大小为60μm×60μm。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤b中阈值滤波是指采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理：

二维中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑；图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:

进行滤波处理时，常常选用不同窗口大小[用n表示]，窗口大小的数值为奇数，优先采用窗口大小为3×3的正方形窗口作为中值滤波的模板，其模板定义公式如下：

式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值；具体计算步骤如下:

第一、将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；

第二、读取出二维窗口中各个像素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值，即中值；

第三、用计算得到的中值替代中值滤波窗口中心点的灰度值。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤c中的阈值选取技术包括最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值技术。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤c中的阈值选取技术采用最大类间方差阈值技术，采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过选取不同的阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值；

设imax，imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为n1(k)、n2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1²(k)、σ2²(k)，全部像素的平均灰度值为μ1，则组内方差σw²(k)：

σw²(k)＝n1σ1²+n2σ2²[3]

组件方差为：

σb²(k)＝n1(μ1-μt)²+n2(μ2+μt)²[4]

则满足的下式的k^*即为所求的阈值，

f(k^*)＝max{f(k)＝σb²(k)/σw²(k)}[5]

式中，

本发明的有益效果是：本发明完整介绍了有机页岩高分辨率成像与数字岩心构建的方法与流程。在骨架与孔隙的阈值求取中，采用两组间性质差异与组内差异之比达到最大值确定阈值，准确实现了孔隙组分、矿物、干酪根等组分的分辨。该方法得到的数字岩心已经在核磁共振、声学孔隙尺度数值模拟与微观特性得到了成功应用。为孔隙结构研究以及后续的微观电学、声学、核磁共振、渗流等岩石物理特性提供了条件。

附图说明

图1为页岩岩心x射线ct扫描单帧灰度图像。

图2为页岩岩心x射线ct扫描3d图像。

图3为页岩岩心x射线ct扫描图未滤波的图像。

图4为页岩岩心x射线ct扫描图滤波后的图像。

图5为页岩ct图像阈值计算界面图。

图6为页岩岩心矿物成分、孔隙流体的分割图像。

图7为重构的页岩三维数字岩心图像。

图8为另一块页岩三维数字岩心图像。

图9为页岩孔隙空间。

图10为页岩干酪根。

图11为页岩岩石骨架。

图12为页岩黄铁矿。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

如图1-12所示，本发明一种有机页岩高分辨率数字岩心构建与分析方法，有机页岩高分辨率数字岩心构建与分析方法按如下步骤实现：

a.页岩样品制备与同步辐射扫描

a1.借助高分辨率电子显微镜选取页岩样品，放入ct仪器载物台进行ct扫描；

a2.扫描图像重构。从扫描图像中剔除背景图像，得到x-y平面的单帧16位页岩岩心灰度图像；

a3.用图像处理软件把多帧16位图像构建为三维图像，观察岩心样本的整体结构特征和孔隙特征；

b.数字图像滤波。对步骤a2所得页岩岩心灰度图像进行图像滤波处理，使岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于ct扫描重建过程中产生的噪声点；

c.ct图像组分分割。对步骤b所得ct图像进行阈值分割，划分岩石的骨架和孔隙：

c1.阈值计算及骨架与孔隙的准确分割

采用阈值分割法，其原理是：岩心ct图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值，通过阈值选取技术选取合理的阈值，将ct图像中的骨架和孔隙准确地分割开来；

c2.骨架矿物的精细分割

对于一般岩石ct图像，只要准确地分割ct图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩ct图像来说，由于其成分复杂，还需要对页岩ct图像的矿物成分进行进一步分割，其主要步骤为：

第一、导入ct图像数据，设置数据格式，将三维ct图像裁剪为需要的尺寸；

第二、依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；

第三、使用灰度阈值分割程序，精细地划分其矿物成分；

第四、导出数据并构建数字岩心，即岩心三维成像图。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤a中同步辐射扫描是指对页岩岩心的高分辨率扫描，其高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm，大视场模式成像视场大小为60μm×60μm。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤b中阈值滤波是指采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理：

二维中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑；图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:

进行滤波处理时，常常选用不同窗口大小[用n表示]，窗口大小的数值为奇数，优先采用窗口大小为3×3的正方形窗口作为中值滤波的模板，其模板定义公式如下：

式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值；具体计算步骤如下:

第一、将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；

第二、读取出二维窗口中各个像素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值，即中值；

第三、用计算得到的中值替代中值滤波窗口中心点的灰度值。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤c中的阈值选取技术包括最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值技术。

所述有机页岩高分辨率数字岩心构建与分析方法，其步骤c中的阈值选取技术采用最大类间方差阈值技术，采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过选取不同的阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值；

设imax，imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为n1(k)、n2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1²(k)、σ2²(k)，全部像素的平均灰度值为μ1，则组内方差σw²(k)：

σw²(k)＝n1σ1²+n2σ2²[3]

组件方差为：

σb²(k)＝n1(μ1-μt)²+n2(μ2+μr)²[4]

则满足的下式的k^*即为所求的阈值，

f(k^*)＝max{f(k)＝σb²(k)/σw²(k)}[5]

式中，

实施例：

1.页岩样品制备与同步辐射扫描

页岩岩心高分辨率扫描是在中国科学院高能物理研究所北京同步辐射装置试验站下的4w1a-x射线成像实验站进行。实验站所用光源是4w1a光束线是北京同步辐射装置最长的一条光束线，可以提供准直性较好的x射线光束，其中纳米成像设备可获得的空间分辨率最高可达30纳米。x射线纳米分辨三维成像实验可以提供高分辨模式和大视场模式这两种实验模式：高分辨模式成像视场大小为15μm×15μm，可达到的分辨率为50nm(标准样品)，大视场模式成像视场大小为60μm×60μm。因为页岩岩样比较致密，同时页岩对于x射线的吸收较强，采用高分辨模式进行ct图像扫描，主要步骤为：

(1)借助高分辨率电子显微镜选作页岩样品，放入ct仪器载物台进行ct扫描；

(2)扫描图像重构。从扫描图像中剔除背景图像，得到x-y平面的单帧16位灰度图像，如图1；

(3)用txm3dviewer软件把多帧16位图像构建为三维图像，以展示其整体特征及形态，如图2。

可以看到，页岩样本的ct扫描图像中小孔隙占整个样本的大多数，多为独立的单独孔隙。

2.阈值滤波

三维页岩岩样切片图像中含有一定的噪声(图3)，需要对图像进行滤波处理。采用二维中值滤波算法对页岩岩心灰度图像进行滤波处理。中值滤波器是一种非线性的滤波技术，它是基于对于信号数值的排序统计理论，对信号进行重新恢复的技术，它可以很好地抑制噪声，同时保留信号中真实的信息和图像的边缘特征。

图像中值滤波是对图像中的灰度点进行运算，将改点灰度值变为该点周围附近窗口内的所有点的中值，实现图像边界的平滑。图像中值的定义如下：假设一组灰度值分别为x1,x2,...,xn，将各个灰度值重新按从大到小的方向排列，其值为:x1',x2',...,xn'，则中值的计算式为:

进行滤波处理时，常常选用不同窗口大小(用n表示)，窗口大小的数值为奇数，采用窗口大小为3*3的正方形窗口作为中值滤波的模板，其模板定义公式如下：

式中为中值滤波后该点的数值，am，n为滤波前窗口内的灰度值。具体计算步骤如下:

(1)将二维窗口在立体图像中按一定方向移动，使二维窗口中心与立体图像中的某个体素重合，该点为中值滤波的目标；

(2)读取出二维窗口中各个体素点所对应的灰度值，将这些灰度值由小到大进行排序，找出排在中间的一个值(中值)；

(3)用计算得到的中值替代中值滤波窗口中心点的灰度值。

通过编写二维中值滤波程序，对岩心ct图像进行滤波处理，其对比效果如图4所示。可以看出，经过中值滤波处理过后，岩石孔隙与骨架之间的边界变得更加明显，更加平滑，也剔除了由于ct扫描重建过程中产生的噪声点，使得孔隙和骨架的灰度值分布更加明显。

3.ct图像组分分割

岩石ct图像的矿物组分分割是生成三维数字岩心的非常重要的一步，对ct图像进行阈值分割，划分岩石的骨架和孔隙。

1)骨架与孔隙的分割与阈值计算

最常用的方法就是阈值分割法。其原理是，岩心ct图像都是由许多的灰度点构成，骨架和孔隙有着明显的灰度差别，通过数学统计的相关方法，求取一个特定的灰度数值，用于划分骨架和孔隙，该值称为阈值。

如何选取合理的阈值就是三维数字岩心能否符合实际岩心的关键。常见的阈值选取技术主要有：最大熵阈值、直方图阈值、模糊阈值和最大类间方差阈值。采用最大类间方差法来划定孔隙和骨架的阈值，将灰度图像中所有像素点的灰度值为一个集合，通过不断选取阈值，将灰度值集合划分成两组，当两组灰度值的组间方差和组内方差之比最大的阈值即为合理值。

设imax，imin二分别代表图像的最小灰度值和最大灰度值，p(i)为灰度值为i的像素个数，k为阈值，它将像素点分为两组，两组的像素数数分别为n1(k)、n2(k)，平均灰度值分别为μ1(k)、μ1(k)，方差分别为σ1²(k)、σ2²(k)，全部像素的平均灰度值为μ1，则组内方差σw²(k)：

σw²(k)＝n1σ1²+n2σ2²(3)

组件方差为：

σb²(k)＝n1(μ1-μt)²+n2(μ2+μt)²(4)

则满足的下式的k^*即为所求的阈值，

f(k^*)＝max{f(k)＝σb²(k)/σw²(k)}(5)式中，

通过所求的阈值，就可以将ct图像中的骨架和孔隙准确地分割开来。

根据上述原理，编写了ct图像阈值求取程序，通过该程序可以实现ct图像不同大小的裁剪，对ct图像再用最大类间方差法求取阈值，并将该代码编写为gui界面，如图5所示。

2)骨架矿物的精细分割

对于一般岩石ct图像，只要准确的分割ct图像的岩石骨架和孔隙空间就已经足够了，但是对于页岩ct图像来说，由于其成分复杂，还需要对页岩ct图像的矿物成分进行进一步分割。进行ct图像处理的主要步骤为：

(1)导入ct图像数据，设置数据格式，将三维ct图像裁剪为需要的尺寸；

(2)依据计算的阈值，对三维数字岩心进行骨架和孔隙的划分，完成初步骨架和孔隙的识别；

(3)使用灰度阈值分割程序，精细的划分其矿物成分，骨架为红色，孔隙为蓝色，黄色为干酪根，浅绿色为黄铁矿，如图6所示；

(4)导出数据并构建数字岩心,绘制岩心三维成像图，如图7所示。

由x射线ct扫描原理可知，如果岩石组分密度较高，则在灰度图像上显示为高亮度，其灰度值也低，如黄铁矿；而密度较低的组分，在灰度图像上显示为低亮度，其灰度值也高，如孔隙内的流体；密度介于两者之间的组分，则在灰度图像上表现为灰白色，例如脆性矿物和黏土矿物。

图8为另一块有机页岩的数字岩心构建实例。从同步辐射图像扫描图像，到纳米成像、数字岩心构建以及组分分析。图中最终处理结果为：黏土矿物为92.2％，石英为3.1％，两者为骨架如图9所示；孔隙约2.4％(图10)，干酪根为2.0％(图11)，黄铁矿为0.3％(图12)。

本发明不局限于上述最佳实施方式，任何人在本发明的启示下得出的其他任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。