CN110910444B - 一种res尺度页岩等效三维孔隙参数快速提取方法 - Google Patents

Abstract

本发明实施例公开了一种页岩等效三维孔隙参数快速提取方法，包括如下步骤：获取页岩样品扫描电镜成像图像；选取分割阈值，将电镜图像划分为四个组分；对孔隙按照连通区域大小进行划分；对较大区域和较小区域提取的有机孔、无机孔进行拼合；根据图像识别结果进行参数提取；根据升维变换，将获取的二维平面参数转换为等效三维参数；获取RES尺度的样品面积大小。该方法处理速度更快、提取结果更准确、有效信息数据获取更多，有利于准确评估页岩的储层储气能力和输送能力，同时能够快速分辨出有机孔、无机孔，并提取出等效三维参数，此外，该方法操作简单、计算自动化完成、计算效率高、适用性较强，能够对各种页岩样品进行分析研究。

Description

一种RES尺度页岩等效三维孔隙参数快速提取方法

技术领域

本发明实施例涉及石油天然气地质勘探技术领域，具体涉及一种RES尺度页岩等效三维孔隙参数快速提取方法。

背景技术

“页岩气革命”将页岩气大规模勘探开发推向热潮，页岩气作为一种清洁的非常规能源拥有巨大的储量。页岩孔隙结构特征的提取是进行储量精确评估、孔隙网络运移能力计算的基础，也是页岩气勘探开发研究的重点。页岩内成分含量复杂，孔隙结构发育，主要成分为有机质、无机质和黄铁矿。在有机质内发育有有机孔、有机裂缝；在无机质内发育有无机孔、无机裂缝。由于有机孔隙和无机孔隙对气体的吸附、流动的影响有着很大的差别，故通过一定的技术手段将其进行分别提取，是很重要的。

获取多孔介质孔隙结构的方法主要有流体注入法和直接观察法。流体注入法是通过低温氮吸附、低温二氧化碳吸附、压汞等手段，依据Kelvin方程建立起来的孔隙结构表征方法，目前已经得到广泛应用，但无法对有机孔隙和无机孔隙进行区分，也无法对裂缝进行单独提取。

直接观察法是基于以扫描电镜为代表的显微技术发展起来的观察方法，能够得到页岩断面的图像直观信息。基于合理的图像处理技术，能够有效将图像所蕴含的信息进行定量化提取。

但是，现有的用于获取页岩多孔介质孔隙结构的方法存在以下缺陷：

(1)现有的处理方法处理速度较慢、提取结果精准度较低、有效信息数据获取有限，不利于准确评估页岩的储层储气能力和输送能力；

(2)现有的处理方法无法对页岩有机孔隙和无机孔隙进行区分，也无法对裂缝进行单独提取，不利于对页岩样品进行准确的分析研究；

(3)现有的处理方法适用性较低，只能处理部分或成像较好的页岩样品，使得研究结果说服力不足。

发明内容

为此，本发明实施例提供一种RES尺度页岩等效三维孔隙参数快速提取方法，该方法处理速度更快、提取结果更准确、有效信息数据获取更多，有利于准确评估页岩的储层储气能力和输送能力，同时，通过对页岩扫描电镜图像进行孔隙结构参数提取，能够快速分辨出有机孔、无机孔，并提取出等效三维参数，有利于对页岩样品进行准确的分析研究，此外，该方法操作简单、计算自动化完成、计算效率极高、适用性较强，能够对各种页岩样品以及不同成像品质的样品进行分析研究，使得研究结果更具说服力，能有效解决现有技术中存在的问题。

为了实现上述目的，本发明实施例提供如下技术方案：一种RES尺度页岩等效三维孔隙参数快速提取方法，包括如下步骤：

S100、页岩样品处理后通过扫描电镜观察，获取样品电镜图像，并对所述样品电镜图像进行预处理；

S200、选取分割阈值，将所述样品电镜图像划分为孔、有机质、无机质、黄铁矿四个组分并标记，依据黄铁矿的标记结果，并提取其重心位置及长宽比信息，辨别由于扫描电镜的边缘增强效应导致孔隙边缘灰度值异常增大现象，记为“白边”，消除所述的“白边”将所述样品电镜图像还原为真实成分的灰度值；

S300、选定连通区域的面积临界值，对孔隙按照连通区域大小进行划分，对小区域进行二次提取，快速筛选并进行标记，区分出小区域内的有机孔和无机孔，

对大区域进行单孔邻域提取，统计单孔邻域内各组分有效含量，并由此区分大区域内的有机孔和无机孔；

S400、对较大区域和较小区域提取的有机孔、无机孔进行图像拼合，得到总体有机孔、无机孔的提取结果，并用bwlabel方法分别进行标记，进一步提取其孔隙结构参数；

S500、根据图像识别结果进行参数提取，分别统计得到孔径参数、孔隙率参数、各组分含量参数，计算孔频率分布情况；

S600、根据升维变换，将获取的二维平面参数转换为等效三维参数；

S700、进行RES判断，获取RES尺度的样品面积大小。

进一步地，在步骤S100中，具体包括：

S101、选取样品，进行样品切割、表面打磨和抛光处理，得到符合扫描电镜观测所需尺寸的样品薄片；

S102、通过扫描电镜观察，获取页岩样品系列扫描电镜成像图像，并将图像裁剪成为规则小块，单个图像大小为1024pixel×1024pixel；

S103、裁剪除去图像黑边(如果有)，经过均值滤波和中值滤波，除去椒盐噪声，并使图像趋于平滑，具体为：

提取图像灰度值非0的长方形区域四顶角位置，将长方形内的区域保留，长方形外的区域裁除；

对裁除后的图像进行滤波处理，先利用均值滤波消除形貌扫描的非光滑性，再用中值滤波消除图像颗粒噪声。

进一步地，在步骤S200中，将电镜图像划分为四个区域的具体步骤为：

S201、选取所裁切图像中的代表性图像，即含有较多孔、有机质、无机质和黄铁矿的图像；

S202、利用二值化方法分别调节孔划分阈值、有机质与无机质分割阈值、无机质与黄铁矿分割阈值；

S203、根据三个分割阈值，对样品扫描电镜图像进行阈值分割，将每一张图像的要素分为孔、有机质、无机质和黄铁矿四组分，将图像灰度值小于孔阈值的区域划分为孔区域；将图像灰度值大于黄铁矿阈值的区域划分为黄铁矿区域；将图像灰度值小于基本阈值但大于孔阈值的区域划分为有机质区域；将图像灰度值大于基本阈值但小于黄铁矿阈值的区域划分为无机质区域。

进一步地，在步骤S200中，进行“白边”消除的具体步骤为：

将样品图像中各个连通的黄铁矿区域进行次序标记，统计各块面积并提取各个黄铁矿区块的边界坐标、重心坐标、长轴长、短轴长、等效半径等；

计算各个黄铁矿区块的长轴与短轴长之比，并统计重心点附近3×3像素矩阵区域大小内，总计9个像素点的平均灰度值；

若黄铁矿区块的长轴和短轴之比大于10，或者重心附近9个像素点平均灰度值小于黄铁矿阈值，则该黄铁矿区块实质为孔缝周缘的“白边”干扰；

依据边界坐标找出包含该区块的最小长方形，给出长方形四个顶点坐标；

从提取得到的最小长方形往外扩展等效半径尺寸，得到新的长方形作为该孔的邻域；

当某向外扩充的邻域边遇到超出整体图像的边界时，以整体图像边界为准，作为邻域的边界，各边判断过程相互独立进行，这就得到了邻域四边形区域。

对各个连通块的邻域进行统计，分别计算各邻域内有机质、无机质的含量，特别是与黄铁矿区相接触的区域的有机质、无机质含量；

统计邻域内接触区有机组分和内部有机组分的含量，进行加权相加，作为有机组分的整体等效含量，其中接触区有机组分权值是内部区有机组分权值的5倍，对无机质含量的统计也遵循该加权方法处理。进一步根据邻域内有机质组分与无机质组分的相对含量判断该“白边”的真实组成，若有机质含量高，认为该“白边”是有机质，将“白边”区域填充为有机质的灰度值；若无机质含量高，则认为该“白边”的真实成分是无机质，将“白边”区域填充为无机质的灰度值。

进一步地，在步骤S300中，对小区域进行两次提取的具体步骤为：

对有机质区域进行尺寸提取，将连通块面积小于临界面积A(像素值)的区域抽提出来，对该较小域按照孔阈值进行二次阈值分割，分割后如仍存在孔域，则为无机较小孔；

S302、提取孔区域的进行尺寸提取，将连通面积小于临界面积A(像素值)的区域抽提出来，作为较小孔区域，对其扣除S301得到的无机较小孔区域，得到有机较小孔区域；

S303、对较小孔区域取孔域的补集，得到较大孔区域。

进一步地，在步骤S300中，还包括有机质孔、无机质孔的区分，具体为：

对孔域的边界点进行提取，并分别找出边界点坐标的x、y方向各两个最大值和最小值，这两个x坐标和两个y坐标分别进行组合，得到四个点坐标，组成了包含孔域的最小四边形；

从提取得到的最小长方形往外扩展等效半径尺寸，得到新的长方形作为该孔的邻域；

当某向外扩充的邻域边遇到超出整体图像的边界时，以整体图像边界为准，作为邻域的边界，各边判断过程相互独立进行，这就得到了邻域四边形区域。

对各个连通块的邻域进行统计，分别计算各邻域内有效有机质、无机质的含量，特别是与孔区相接触的区域的有机质、无机质含量；

统计邻域内接触区有机组分和内部有机组分的含量，进行加权相加，作为有机组分的整体等效含量，其中接触区有机组分权值是内部区有机组分权值的5倍，对无机质含量的统计也遵循该加权方法处理。进一步根据邻域内有机质组分与无机质组分的相对含量判断该孔的成分属性，若有机质含量高，认为该孔是有机质孔；若无机质含量高，则认为该孔是无机质孔,这样就进行了大孔域有机孔、无机孔的筛分。

进一步地，在步骤S500中，具体包括：

统计孔径参数，取提取的有机孔和无机孔区域进行标记，分别求每一单孔实际等效直径，对有机孔、无机孔和总孔分别求其平均值，并按照5nm的统计间隔做出孔频率分布柱状图；

各组分含量参数，具体包括有机质含量、无机质含量、黄铁矿含量，统计有机质区域的面积，除以样品图像总面积，得到有机质含量；统计无机质区域面积，除以样品图像总面积，得到无机质含量；统计有效的实际黄铁矿区域面积，除以图像样品总面积，得到黄铁矿含量；

统计孔隙率参数，具体包括总孔隙率、有机质孔隙率、无机质孔隙率，将总孔隙面积除以样品总面积，得到总孔隙率；将有机孔面积除以有机质面积得到有机孔隙率；将无机孔面积除以无机质面积得到无机孔隙率。

进一步地，在步骤S600中，将二维平面参数转换为等效三维参数时，

建立三维孔随机截面投影到二维平面时的概率分布函数：

其中r是二维孔径，R是三维孔径，P(r)表示孔径为R的三维孔，经过二维随机截面投影得到孔径为r的孔的概率；

则二维随机截面投影的概率分布的平均值：

将步骤S500中获取的孔径尺寸相关的值，统计除以π/4再进行深入计算求解，即进行升维转换；

对孔径分布进行积分变换，让二维频率分布等效转换为三维频率分布，转换关系为：

式中N(R)是三维孔的孔数量分布函数，n(r)是二维孔的数量分布函数；

通过线性拟合，将根据二维SEM提取得到的n(r)进行多项式拟合，并且假设N(R)与n(r)具有相同的分布形式，通过逐阶积分比对，分别求出N(R)所具有的系数，从而确定N(R)，将二维分布转换到三维。

进一步地，在步骤S700中，具体包括：

S701、对步骤S100至步骤S600的重复过程进行计次，每计算一张小图，为一次重复，当重复次数超过100次时，假定此时重复次数为N；

S702、对前[N/2]个计算数据进行求平均，对后[N/2]个计算数据进行求平均，对中间[N/2]个数据进行求平均，得到三组用于对比的数据；

S703、对获取的三组对比数据进行分析。

进一步地，在步骤S703中：

S704、若三组数据两两之间的偏差均在1％以内，则表示所计算的尺寸有可能达到了RES尺度，按照步骤S705进行验证；若三组数据两两之间的偏差未达到要求，则进行步骤S706；

S705、若此时所取的计算样品数目为Ns个，对中间Ns个数据的求平均，改为对所有数据点中任意位置取Ns个点，对这些点求平均，得到三个平均值进行对比，若两两之间的偏差仍然小于1％，则Ns为取得代表性尺寸的样品数目；若两两之间的偏差大，则选取新图像或者增大Ns，继续重复步骤S701、步骤S702的计算步骤；

S706、若此时N还未达到图像总数，则选取新图像，继续按照步骤S701至步骤S704进行运算；如果N已经达到图像总数，则从i取1开始，取样品计算数目变为[N/2]+i，进行步骤S702至步骤S704的计算过程；

S707、基于上述过程，获取代表性尺寸所满足的样品数目Ns，统计Ns个样品的总像素面积，并将结果乘以分辨率的平方，即可得到RES尺度的样品面积大小。

本发明实施例具有如下优点：

(1)本发明处理速度更快、提取结果更准确、有效信息数据获取更多，有利于准确评估页岩的储层储气能力和输送能力；

(2)本发明通过对页岩扫描电镜图像进行孔隙结构参数提取，能够快速分辨出有机孔、无机孔，并提取出等效三维参数，有利于对页岩样品进行准确的分析研究；

(3)本发明操作简单、计算自动化完成、计算效率极高、适用性较强，能够对各种页岩样品以及不同成像品质的样品进行分析研究，使得研究结果更具说服力。

附图说明

为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。

本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。

图1为本发明的整体流程结构示意图；

图2为本发明的算法流程图；

图3为本发明的页岩样品扫描电镜灰度图；

图4为本发明的页岩样品有机孔提取结果图；

图5为本发明的页岩样品无机孔提取结果图；

图6为本发明的页岩样品综合提取结果图；

图7为本发明的页岩样品的有机孔、无机孔频率分布图；

图8为本发明的页岩样品的有机孔、无机孔孔径分布图；

图9为本发明的页岩样品的孔径分布结果与氮吸附结果对比图。

具体实施方式

以下由特定的具体实施例说明本发明的实施方式，熟悉此技术的人士可由本说明书所揭露的内容轻易地了解本发明的其他优点及功效，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1和图2所示，本发明提供了一种RES尺度页岩等效三维孔隙参数快速提取方法，包括如下步骤：

S100、页岩样品处理后通过扫描电镜观察，获取样品电镜图像，并对所述样品电镜图像进行预处理；

在步骤S100中，具体包括：

步骤S101、选取样品，进行样品切割、表面打磨和抛光处理，得到符合扫描电镜观测所需尺寸的样品薄片。

将所取得的页岩岩心切割成1cm×0.5cm×0.5cm的块体，选尺寸1cm×0.5cm的面作为观测面，进行打磨处理，并放入氩离子抛光设备进行抛光处理。抛光完毕之后，将样品置于扫描电镜设备的观察仓内进行扫描电镜观察，本例的图像大小为398.5μm×398.5μm。

步骤S102、通过扫描电镜观察，获取页岩样品系列扫描电镜成像图像，并将图像裁剪成为规则小块，单个图像大小为1024pixel×1024pixel，图3为本实施例的页岩样品扫描电镜灰度图。

步骤S103、裁剪除去图像黑边(去除图像边界未扫描到样品信息的黑域)，经过均值滤波和中值滤波，除去椒盐噪声，并使图像趋于平滑，具体为：

提取图像灰度值非0的长方形区域四顶角位置，将长方形内的区域保留，长方形外的区域裁除；对裁除后的图像进行滤波处理，先利用均值滤波消除形貌扫描的非光滑性，使得形貌的提取结果更符合真实视觉效果，再用中值滤波消除图像颗粒噪声，除去单像素异常灰度值带来的干扰。

均值滤波具体为：生成一个矩阵大小为3×3(或5×5)的滤波模板，对目标点及周围9(或25)个像素求取平均值，由该平均值代替目标点的像素值，对图像所有像素点依次求均，则得到滤波后的图像。

中值滤波具体为：生成一个矩阵大小为3×3(或5×5)的模板，对目标点及其周围9(或25)个像素求取中值，由该中值代替目标点的像素值，对图像所有像素点依次求中值，则得到滤波后的图像。

步骤S200、选取分割阈值，将电镜图像划分为四个组分，并提取形状因子进行“白边”消除。

找出进行图像分割的三个阈值，初步将页岩的四组分进行划分，这四组分包括：孔、有机质、无机质、黄铁矿。根据黄铁矿的标记结果，提取其重心位置及长宽比，辨别出由边缘增强效应带来的孔隙周边“白边”问题，将其还原为真实成分的灰度值(有机质组分或者无机质组分的灰度值)，其中有机孔、无机孔灰度值最小(图像中呈现黑色)，并且无法通过单一灰度值进行区分，有机质灰度值较大一些(图像中呈现黑灰色)，无机质的灰度值更大一些(图像中呈现灰色)，黄铁矿和“白边”的灰度值最大(图像中呈现灰白色或白色)。

在步骤S200中，将电镜图像样品划分为四个区域的具体步骤为：

步骤S201、选取所裁切后图像中的代表性图像，即含有较多孔、有机质、无机质和黄铁矿的图像；

步骤S202、利用二值化方法分别调节孔划分阈值、有机质与无机质分割阈值、无机质与黄铁矿分割阈值。

首先对孔阈值进行调整提取，取不同大小的灰度临界值作为阈值，按照两个标准判断该灰度值的好坏：一是提取的孔是否与肉眼所见相符合；二是得到的面孔率是否与测得的孔隙度相接近，选取最好的灰度值作为最终划分阈值。

在选择有机质与无机质的分割阈值时，孔阈值调节方法，取大于孔阈值的不同大小的灰度临界值进行提取试验，提取的有机质区域的结果与肉眼所识别结果最接近的灰度值为区分有机质与无机质的基本阈值。

在选择无机质与黄铁矿的分割阈值时，参照有机质与无机质的分割阈值方法，取大于基本阈值的不同大小的灰度临界值进行提取试验，提取的黄铁矿区域的结果与肉眼所识别的结果最接近的灰度值作为区分无机质与黄铁矿的黄铁矿阈值。

步骤S203、根据三个分割阈值，对样品扫描电镜图像进行阈值分割，将每一张图像的要素分为孔、有机质、无机质和黄铁矿四组分，将图像灰度值小于孔阈值的区域划分为孔区域；将图像灰度值大于黄铁矿阈值的区域划分为黄铁矿区域；将图像灰度值小于基本阈值但大于孔阈值的区域划分为有机质区域；将图像灰度值大于基本阈值但小于黄铁矿阈值的区域划分为无机质区域。

在步骤S200中，进行“白边”消除的具体步骤为：

将样品图像中各个连通的黄铁矿区域进行次序标记，统计各块面积并提取各个黄铁矿区块的边界坐标、重心坐标、长轴长、短轴长、等效半径等；

计算各个黄铁矿区块的长轴与短轴长之比，并统计重心点附近3×3区域大小内，总计9个像素点的平均灰度值；

若黄铁矿区块的长轴和短轴之比大于10，或者重心附近9个像素点平均灰度值小于黄铁矿阈值，则该黄铁矿区块实质为孔缝周缘的“白边”干扰；

依据边界坐标找出包含该区块的最小长方形，给出长方形四个顶点坐标；

从提取得到的最小长方形往外扩展等效半径尺寸，得到新的长方形作为该孔的邻域；

当某向外扩充的邻域边遇到超出整体图像的边界时，以整体图像边界为准，作为邻域的边界，各边判断过程相互独立进行，这就得到了邻域四边形区域。

对各个连通块的邻域进行统计，分别计算各邻域内有机质、无机质的含量，特别是与黄铁矿区相接触的区域的有机质、无机质含量；

统计邻域内接触区有机组分和内部有机组分的含量，进行加权相加，作为有机组分的整体等效含量，其中接触区有机组分权值是内部区有机组分权值的5倍，对无机质含量的统计也遵循该加权方法处理。进一步根据邻域内有机质组分与无机质组分的相对含量判断该“白边”的真实组成，若有机质含量高，认为该“白边”是有机质，将“白边”区域填充为有机质的灰度值；若无机质含量高，则认为该“白边”的真实成分是无机质，将“白边”区域填充为无机质的灰度值。

步骤S300、对孔隙按照连通区域大小进行划分，对小区域进行二次提取，快速筛选并进行标记，区分出有机较小孔和无机较小孔，对大区域进行邻域提取，统计邻域内各组分有效含量，并由此区分有机较大孔和无机较大孔。

在步骤S300中，对小区域进行两次提取的具体步骤为：

对有机质区域进行尺寸提取，将连通块面积小于临界面积A(像素值)的区域抽提出来，对该较小域按照孔阈值进行二次阈值分割，分割后如仍存在孔域，则为无机较小孔；

S302、提取孔区域的进行尺寸提取，将连通面积小于临界面积A(像素值)的区域抽提出来，作为较小孔区域，对其扣除S301得到的无机较小孔区域，得到有机较小孔区域；

步骤S303、对较小孔区域取孔域的补集，得到较大孔区域。

在步骤S300中，还包括有机质孔、无机质孔的区分，具体为：

对孔域的边界点进行提取，并分别找出边界点坐标的x、y方向各两个最大值和最小值，这两个x坐标和两个y坐标分别进行组合，得到四个点坐标，组成了包含孔域的最小四边形；

从提取得到的最小长方形往外扩展等效半径尺寸，得到新的长方形作为该孔的邻域；

当某向外扩充的邻域边遇到超出整体图像的边界时，以整体图像边界为准，作为邻域的边界，各边判断过程相互独立建。这就得到了邻域四边形区域。

对各个连通块的邻域进行统计，分别计算各邻域内有效有机质、无机质的含量，特别是与孔区相接触的区域的有机质、无机质含量；

统计邻域内接触区有机组分和内部有机组分的含量，进行加权相加，作为有机组分的整体等效含量，其中接触区有机组分权值是内部区有机组分权值的5倍，对无机质含量的统计也遵循该加权方法处理。进一步根据邻域内有机质组分与无机质组分的相对含量判断该孔的成分属性，若有机质含量高，认为该孔是有机质孔；若无机质含量高，则认为该孔是无机质孔,这样就进行了大孔域有机孔、无机孔的筛分。

步骤S400、对较大区域和较小区域提取的有机孔、无机孔进行拼合，得到总体有机孔、无机孔的提取结果，并用bwlabel方法分别进行标记。

将较大有机孔与较小有机孔的识别结果进行拼合，得到有机孔的整体数据；将较大无机孔与较小无机孔的识别结果进行拼合，得到无机孔的整体数据。经过以上步骤，即可以把有机孔、无机孔、黄铁矿都提取出来。本实施例的页岩样品有机孔提取结果图如图4所示；本实施例的页岩样品无机孔提取结果图如图5所示；本实施例的页岩样品综合提取结果图如图6所示。

步骤S500、根据图像识别结果进行参数提取，分别统计得到孔径参数、孔隙率参数、各组分含量参数，计算孔频率分布情况。

在步骤S500中，具体包括：

统计孔径参数，取提取的有机孔和无机孔区域进行标记，分别求每一单孔实际等效直径，对有机孔、无机孔和总孔分别求其平均值，并按照5nm的统计间隔做出孔频率分布柱状图；

各组分含量参数，具体包括有机质含量、无机质含量、黄铁矿含量，统计有机质区域的面积，除以样品图像总面积，得到有机质含量；统计无机质区域面积，除以样品图像总面积，得到无机质含量；统计有效的实际黄铁矿区域面积，除以图像样品总面积，得到黄铁矿含量；

统计孔隙率参数，具体包括总孔隙率、有机质孔隙率、无机质孔隙率，将总孔隙面积除以样品总面积，得到总孔隙率；将有机孔面积除以有机质面积得到有机孔隙率；将无机孔面积除以无机质面积得到无机孔隙率。

本实施例得到的有机孔隙率为10.7234％，无机孔隙率为0.3303％，总孔隙率为1.7369％，有机质含量为13.552％，无机质含量为85.888％，黄铁矿含量为0.561％。由提取得到的孔径数据，取相应的平均值，计算得到孔隙的平均孔径、有机孔的平均孔径、无机孔的平均孔径，本实施例得到的平均孔径为32.506nm，平均有机孔径为32.235nm，平均无机孔径为36.820nm。对各类孔的提取结果进行孔径统计，得到各类孔的频率分布图谱，如：有机孔、无机孔1-3000nm的孔频率分布图。

步骤S600、根据升维变换，将获取的二维平面参数转换为等效三维参数，得到相应的三维孔隙率、含量参数，三维孔频率分布结果和三维孔径分布结果。图7为本实施例的页岩样品的有机孔、无机孔频率分布图；图8为本实施例的页岩样品的有机孔、无机孔孔径分布图；图9为本发明的页岩样品的孔径分布结果与氮吸附结果对比图。

在步骤S600中，将二维平面参数转换为等效三维参数时，

建立三维孔随机截面投影到二维平面时的概率分布函数：

其中r是二维孔径，R是三维孔径，P(r)表示孔径为R的三维孔，经过二维随机截面投影得到孔径为r的孔的概率；

则二维随机截面投影的概率分布的平均值：

将步骤S500中获取的孔径尺寸相关的值，统计除以π/4再进行深入计算求解，即进行升维转换；

对孔径分布进行积分变换，让二维频率分布等效转换为三维频率分布，转换关系为：

式中N(R)是三维孔的孔数量分布函数，n(r)是二维孔的数量分布函数；

通过线性拟合，将根据二维SEM提取得到的n(r)进行多项式拟合，并且假设N(R)与n(r)具有相同的分布形式，通过逐阶积分比对，分别求出N(R)所具有的系数，从而确定N(R)，将二维分布转换到三维。

步骤S700、进行RES判断，获取RES尺度的样品面积大小。

在步骤S700中，具体包括：

步骤S701、对步骤S100至步骤S600的重复过程进行计次，每计算一张小图，为一次重复，当重复次数超过100次时，假定此时重复次数为N；

步骤S702、对前[N/2[个计算数据进行求平均，对后[N/2[个计算数据进行求平均，对中间[N/2[个数据进行求平均(如果[N/2[是奇数，而N为偶数时，无法对中间取对称点，以方便起见，先以[N/2[-1取对称点，再在数据最前面取一点)，得到三组用于对比的数据；

步骤S703、对获取的三组对比数据进行分析。

在步骤S703中，具体包括：

步骤S704、若三组数据两两之间的偏差均在1％以内，则表示所计算的尺寸有可能达到了RES尺度，按照步骤S705进行验证；若三组数据两两之间的偏差未达到要求，则进行步骤S706；

步骤S705、若此时所取的计算样品数目为Ns个，对中间Ns个数据的求平均，改为对所有数据点中任意位置取Ns个点，对这些点求平均，得到三个平均值进行对比，若两两之间的偏差仍然小于1％，则Ns为取得代表性尺寸的样品数目；若两两之间的偏差大，则选取新图像或者增大Ns，继续重复步骤S701、步骤S702的计算步骤；

步骤S706、若此时N还未达到图像总数，则选取新图像，继续按照步骤S701至步骤S704进行运算；如果N已经达到图像总数，则从i取1开始，取样品计算数目变为[N/2]+i，进行步骤S702至步骤S704的计算过程；

步骤S707、基于上述过程，获取代表性尺寸所满足的样品数目Ns，统计Ns个样品的总像素面积，并将结果乘以分辨率的平方，即可得到RES尺度的样品面积大小。

本发明实施例具有如下优点：

(1)本发明处理速度更快、提取结果更准确、有效信息数据获取更多，有利于准确评估页岩的储层储气能力和输送能力；

(2)本发明通过对页岩扫描电镜图像进行孔隙结构参数提取，能够快速分辨出有机孔、无机孔，并提取出等效三维参数，有利于对页岩样品进行准确的分析研究；

(3)本发明操作简单、计算自动化完成、计算效率极高、适用性较强，能够对各种页岩样品以及不同成像品质的样品进行分析研究，使得研究结果更具说服力。

虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (9)

1.一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，包括如下步骤：

S100、页岩样品处理后通过扫描电镜观察，获取样品电镜图像，并对所述样品电镜图像进行预处理；

S200、选取分割阈值，将所述样品电镜图像划分为孔、有机质、无机质、黄铁矿四个组分并初步标记，依据黄铁矿的标记结果，并提取其重心位置及长宽比信息，辨别由于扫描电镜的边缘增强效应导致孔隙边缘灰度值异常增大现象，记为“白边”，消除所述的“白边”将所述样品电镜图像还原为真实成分的灰度值；

S300、选定连通区域的面积临界值，对孔隙按照连通区域大小进行划分：对小区域进行二次提取，快速筛选并进行标记，区分出有机较小孔和无机较小孔；对大区域进行单孔邻域提取，统计邻域内各组分有效含量，并由此区分有机较大孔和无机较大孔；

在步骤S300中，对小区域进行两次提取的具体步骤为：

S301、对有机质区域进行尺寸提取，将连通块面积小于面积临界值的区域抽提出来，对该较小域按照孔阈值进行二次阈值分割，分割后如仍存在孔域，则为无机较小孔；

S302、提取孔区域的进行尺寸提取，将连通面积小于面积临界值的区域抽提出来，作为较小孔区域，对其扣除S301得到的无机较小孔区域，得到有机较小孔区域；

S303、对较小孔区域取孔域的补集，得到较大孔区域；

S400、对较大区域和较小区域提取的有机孔、无机孔进行图像拼合，得到总体有机孔、无机孔的提取结果，并用bwlabel方法分别进行标记，进一步提取其参数；

S500、根据图像识别结果进行参数提取，分别统计得到孔径参数、孔隙率参数、各组分含量参数，计算孔频率分布情况；

S600、根据升维变换，将获取的二维平面参数转换为等效三维参数；

S700、进行RES判断，获取RES尺度的样品面积大小。

2.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S100中，具体包括：

S101、选取样品，进行样品切割、表面打磨和抛光处理，得到符合扫描电镜观测所需尺寸等要求的页岩样品薄片；

S102、通过扫描电镜观察，获取页岩样品系列扫描电镜成像图像，并将图像裁剪成为规则小块，单个图像大小为1024pixel×1024pixel；

S103、裁剪除去图像黑边，经过均值滤波和中值滤波，除去椒盐噪声，并使图像趋于平滑。

3.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S200中，将电镜图像划分为四个区域的具体步骤为：

S201、选取所裁切图像中的代表性图像，即含有较多孔、有机质、无机质和黄铁矿的图像；

S202、利用二值化方法分别调节孔划分阈值、有机质与无机质分割阈值、无机质与黄铁矿分割阈值；

S203、根据三个分割阈值，对样品扫描电镜图像进行阈值分割，将每一张图像的要素分为孔、有机质、无机质和黄铁矿四组分，将图像灰度值小于孔阈值的区域划分为孔区域；将图像灰度值大于黄铁矿阈值的区域划分为黄铁矿区域；将图像灰度值小于基本阈值但大于孔阈值的区域划分为有机质区域；将图像灰度值大于基本阈值但小于黄铁矿阈值的区域划分为无机质区域。

4.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S200中，进行“白边”消除的具体步骤为：

将样品图像中各个连通的黄铁矿区域进行次序标记，统计各块面积并提取各个黄铁矿区块的边界坐标、重心坐标、长轴长、短轴长、等效半径；

计算各个黄铁矿区块的长轴与短轴长之比，并统计重心点及其附近3×3像素矩阵区域内，总计9个像素点的平均灰度值；

若黄铁矿区块的长轴和短轴之比大于10，或者重心附近9个像素点平均灰度值小于黄铁矿阈值，则判断该黄铁矿区块的实质为孔缝周缘的“白边”干扰，否则为真实黄铁矿；

依据该“黄铁矿区域”的边界坐标找出包含该区块的最小长方形，并给出长方形四个顶点坐标；

从提取得到的最小长方形四条边分别往外扩展“黄铁矿区域”的等效半径尺寸，得到新的长方形作为该孔的邻域；

当某向外扩充的邻域边遇到超出整体图像的边界时，以整体图像边界为准，作为邻域的边界，各孔邻域生成独立进行，这就得到了邻域四边形区域；

对各个连通“黄铁矿”块的邻域进行统计，分别计算各邻域内与黄铁矿区相接触的区域的有机质、无机质含量；

统计邻域内接触区有机组分和内部有机组分的含量，进行加权相加，作为有机组分的整体等效含量，其中接触区有机组分权值是内部区有机组分权值的5倍，对无机质含量的统计也遵循该加权方法处理，进一步根据邻域内有机质组分与无机质组分的相对含量判断该“白边”的真实组成，若有机质含量高，认为该“白边”是有机质，将“白边”区域填充为有机质的灰度值；若无机质含量高，则认为该“白边”的真实成分是无机质，将“白边”区域填充为无机质的灰度值。

5.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S300中，还包括较大孔区域有机质孔、无机质孔的区分，具体为：

对较大孔区域的边界点进行提取，并分别找出边界点坐标的x、y方向各自最大值和最小值，这两个x坐标和两个y坐标分别进行组合，得到四个点坐标，组成了包含孔域的最小四边形；

从提取得到的最小长方形往外扩展孔区域等效半径的尺寸，得到新的长方形作为该孔的邻域；

当某向外扩充的邻域边遇到超出整体图像的边界时，以整体图像边界为准，作为邻域的边界，且各边判断过程互相独立进行，由此得到邻域四边形区域；

对各个连通大孔域的邻域分别进行统计，计算各邻域内有效有机质、无机质的含量，具体包括与孔区相接触的区域的有机质、无机质含量；

统计邻域内接触区有机组分和内部有机组分的含量，进行加权相加，作为有机组分的整体等效含量，其中接触区有机组分权值是内部区有机组分权值的5倍；对无机质含量的统计也遵循该加权方法处理；

进一步根据邻域内有机质组分与无机质组分的相对含量判断该孔的成分属性，若有机质含量高，认为该孔是有机质孔；若无机质含量高，则认为该孔是无机质孔。

6.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S500中，具体包括：

统计孔径参数，取提取的有机孔和无机孔区域进行标记，分别求每一单孔实际等效直径值，按照5nm的统计间隔做出孔频率分布柱状图，并对有机孔、无机孔和总孔分别求其孔径的平均值；

各组分含量参数，具体包括有机质含量、无机质含量、黄铁矿含量，统计有机质区域的面积，除以样品图像总面积，得到有机质含量；

统计无机质区域面积，除以样品图像总面积，得到无机质含量；

统计有效的实际黄铁矿区域面积，除以图像样品总面积，得到黄铁矿含量；

统计孔隙率参数，具体包括总孔隙率、有机质孔隙率、无机质孔隙率，将总孔隙面积除以样品总面积，得到总孔隙率；

将有机孔面积除以有机质面积得到有机孔隙率；

将无机孔面积除以无机质面积得到无机孔隙率。

7.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S600中，将二维平面参数转换为等效三维参数时，

建立三维孔随机截面投影到二维平面时的概率分布函数：

其中r是二维孔径，R是三维孔径，P(r)表示孔径为R的三维孔，经过二维随机截面投影得到孔径为r的孔的概率；

则二维随机截面投影的概率分布的平均值

将步骤S500中获取的孔径尺寸相关的值，统计除以π/4再进行深入计算求解，即进行升维转换；

对孔径分布进行积分变换，让二维频率分布等效转换为三维频率分布，转换关系为：

式中N(R)是三维孔的孔数量分布函数，n(r)是二维孔的数量分布函数；

通过线性拟合，将根据二维SEM提取得到的n(r)进行多项式拟合，并且假设N(R)与n(r)具有相同的分布形式，通过逐阶积分比对，分别求出N(R)所具有的系数，从而确定N(R)，将二维分布转换到三维等效结果。

8.根据权利要求1所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S700中，具体包括：

S701、对步骤S100至步骤S600的重复过程进行计次，每计算一张小图，为一次重复，当重复次数超过100次时，假定此时重复次数为N；

S702、对前[N/2]个计算数据进行求平均，对后[N/2]个计算数据进行求平均，对中间[N/2]个数据进行求平均，得到三组用于对比的数据；

定义[]为对方括号内数值取整；

S703、对获取的三组对比数据进行分析。

9.根据权利要求8所述的一种RES尺度页岩等效三维孔隙参数快速提取方法，其特征在于，在步骤S703中：

S704、若三组数据两两之间的偏差均在5％以内，则表示所计算的尺寸达到了RES尺度，按照步骤S705进行验证；若三组数据两两之间的偏差未达到要求，则进行步骤S706；

S705、若此时所取的计算样品数目为Ns个，对中间Ns个数据的求平均，改为对所有数据点中任意位置取Ns个点，对这些点求平均，得到三个平均值进行对比，若两两之间的偏差仍然小于1％，则Ns为取得代表性尺寸的样品数目；若两两之间的偏差大，则选取新图像或者增大Ns，继续重复步骤S701、步骤S702的计算步骤；

S706、若此时N还未达到图像总数，则选取新图像，继续按照步骤S701至步骤S704进行运算；如果N已经达到图像总数，则从i取1开始，取样品计算数目变为[N/2]+i，进行步骤S702至步骤S704的计算过程；

S707、基于上述过程，获取代表性尺寸所满足的样品数目Ns，统计Ns个样品的总像素面积，并将结果乘以分辨率的平方，即可得到RES尺度的样品面积大小。

Images