CN107993261B - 一种基于三维岩心扫描图像的孔隙与孔喉识别方法 - Google Patents

Abstract

本发明涉及石油天然气地质勘探技术，涉及一种基于三维岩心扫描图像的孔隙和孔喉识别方法，所述孔隙和孔喉识别方法为从三维岩心扫描图像中分割出孔隙空间，从所述孔隙空间对孔隙和孔喉进行分割实现对孔隙和孔喉的识别；所述孔喉为连通两个所述孔隙的通道。这一发明可以克服现有技术中无法精确且高效地分割孔隙和孔喉，能够对孔隙和孔喉特征进行定量分析，可以更好地评价储层储油储气能力、孔隙内的油和气的赋存状态、运移能力。

Description

一种基于三维岩心扫描图像的孔隙与孔喉识别方法

技术领域

本发明涉及石油天然气地质勘探技术，特别是一种基于三维岩心扫描图像的孔隙与孔喉识别的方法。

背景技术

储集层岩石的微观孔隙结构特征(孔隙空间的大小、几何形状、分布及相互连通情况)直接反应储集层的渗流能力，可为储集层评价提供重要依据。储集层岩石的孔隙、孔喉特征对油气的勘探开发至关重要，直接关系到储集层岩石的储油和储气能力、油和气的赋存状态、运移机理和能力以及储集层岩石的沉积演化过程，最终影响油气开采的难易程度和勘探开发价值。因此精确厘定储集层岩石的孔隙和孔喉特征意义重大。只有准确了解了储集层岩石孔隙和孔喉特征才能正确反应储集层岩石的渗流能力，并制定出正确的勘探开发方案。

储集层岩石中未被矿物颗粒、胶结物或者其它固体物质填集的空间称为岩石的孔隙空间。孔隙空间可以分为孔隙和孔喉。由岩石颗粒包围着的较大的空间称为孔隙，而连通两个孔隙的狭窄的通道称为孔喉，如图1所示。孔隙、孔喉的储油和储气能力不同，油和气的赋存状态、运移机理和能力也不相同，因此需要分别分析。

目前对于孔隙和孔喉的分割已开展了一些研究，其中主要方法有多向搜索算法、基于中轴的算法和最大球算法。多向搜索算法是向有限的几个面进行搜索，所得结果并不精确，部分孔喉无法搜到；基于中轴的算法的分割效果依赖中轴提取的结果，若提取中轴效果不好，孔隙和孔喉也无法准确分割；基于最大球算法对于噪声非常敏感，会出现孔喉过度分割的情况。鉴于此，分割的方法必须创新。本文采用基于三维图像骨架化的方法分割孔隙和孔喉，简言之，采用细化算法(thinner)将孔隙空间骨架化，根据孔隙空间骨架上点的厚度信息将厚度相对较大的区域定为孔隙，将连接两个孔隙且厚度相对较小的区域定为孔喉，由此，将孔隙和孔喉区分开。

发明内容

针对上述问题，本发明提供一种基于三维岩心扫描图像的孔隙和孔喉识别的方法。本发明可以克服现有技术中无法精确且高效地分割孔隙和孔喉，能够对孔隙和孔喉特征进行定量分析，可以更好地评价储层储油储气能力、孔隙内的油和气的赋存状态、运移能力。

本发明是通过以下技术方案实现的：

一种基于三维岩心扫描图像的孔隙和孔喉识别方法，所述孔隙和孔喉识别方法为从三维岩心扫描图像中分割出孔隙空间，从所述孔隙空间对孔隙和孔喉进行分割实现对孔隙和孔喉的识别；所述孔喉为连通两个所述孔隙的通道。

进一步地，从所述孔隙空间对孔隙和孔喉进行分割具体包括以下步骤：

骨架线提取：在所述孔隙空间提取出相互连通区域的骨架线；所述骨架线上的每个点是严格的一个体素宽度，每一个所述点包含空间坐标信息和厚度；

所述空间坐标信息为X、Y、Z值；所述厚度为该点距离最近的非孔隙边界的距离；

分段：将所述骨架线的端点和交叉点定义为节点，根据节点断开，将所述骨架线断开成若干线段，每一所述线段都包括多个所述点；单线段分割：对其中一个所述线段进行孔隙和孔喉分割，所述单线段分割包括疑似孔喉点边界判断和孔喉确定；

所述疑似孔喉点判断为，当所述点的厚度与相邻点的最大厚度的比值小于或者等于a，即被认为是相对狭小区域的边界点，所述相对狭小区域的边界点为疑似孔喉边界点；

所述孔喉确定为，

在所述线段上，存在依次排列的点o、p、q、r，其中p点和q点重合或不重合，p点和q点为所述疑似孔喉边界点；

o、p、q、r的厚度满足：d_p/d₀小于或者等于a；且d_q/d_r小于或者等于a；且p、q之间不存在点s使得d_p/d_s或d_q/d_s小于或者等于a，则，认定o、q两点之间为孔喉；所述0<a<1；

全部分割：针对其他所述线段，采用所述单线段分割的方法对线段内的孔隙和孔喉进行分割。

进一步地，所述单线段分割具体包括以下步骤：

第一步，提取某一所述线段内所有所述点的厚度，按照点的顺序将所有厚度存于一个数列中得到厚度数列d＝{d₁,…,d_i,…d_n}，d_i为第i个点的厚度，i∈(1,n)，n为所述线段内点的个数；

若所述厚度数列中存在相邻的厚度相同时只保留一个厚度组成简化厚度数列d′＝{d′₁,…,d′_i,…d′_m}，1<m≤n；

第二步，将所述简化厚度数列进行极值化，只保留所述简化厚度数列中的极大值和极小值得到极值序列d″＝{d″₁…,d″_i,…d″_f}，1<f<m；

所述极值序列只包含极大值和极小值，且极大值和极小值交替出现；

所述极大值是指厚度比与其数列中相邻的两个点都大的点；

所述极小值是指厚度比与其数列中相邻的两个点都小的点；

其中端点只与相邻的一个点进行比较；

第三步，孔喉点判断，分析所述极值序列，引入两个虚拟值：V_min和V_max，

将所述极值序列中的第一个极小值点赋值给V_min，将所述极值序列中的第一个极大值点赋值给V_max，计算V_min/V_max；

针对d″_i的判断为：若d″_i为极大值点，则比较d″_i与V_max，若d″_i>V_max，则重新赋值V_max＝d″_i；若d″_i为极小值点，则比较d″_i与V_min，若d″_i<V_min，则重新赋值V_min＝d″_i；重新赋值后计算V_min/V_max；i∈[1，f]

当V_min/V_max大于a值时，进行d″_i+1的判断，d″_i+1的判断方法与d″_i的判断方法相同；

当V_min/V_max小于等于a值时，判定d″_i点为疑似孔喉点，同时相应的将V_min或V_max清空后进行d″_i+1的判断；

其中所述将V_min或V_max清空选择为：当d″_i为极大值时，将V_min清空；当d″_i为极小值时，将V_max清空；

依次对所述极值序列中所有点进行判断，将所述单线段中疑似孔喉边界点全部找出。

进一步地，所述孔喉确定中，

p点和q点仅选择在所述极值序列中的极小值点；

o、s、r点选择在所述极值序列中的极大值点。

进一步地，从三维岩心扫描图像中分割出孔隙空间具体为采用图像阈值分割法对三维岩心扫描图像的矿物区和孔隙空间进行分割。

进一步地，所述三维岩心扫描图像利用岩心扫描设备对岩心进行扫描得到；

进一步地，在采用图像阈值分割法前对所述三维岩心扫描图像进行图像滤波处理。

进一步地，完成孔隙和孔喉分割后进一步分析可得到孔隙参数和孔喉参数。

进一步地，所述孔隙参数包括孔隙体积、孔隙等效直径、孔隙平均直径以及孔隙包含的具体点的信息。

进一步地，所述孔喉参数包括孔喉体积、孔喉等效直径、孔喉平均直径以及孔喉包含的具体点的信息。

本发明的有益技术效果：本发明可以克服现有技术中无法精确且高效地分割孔隙和孔喉，能够对孔隙和孔喉特征进行定量分析，可以更好地评价储层储油储气能力、孔隙内的油和气的赋存状态、运移能力。

附图说明

图1为本发明孔隙孔喉示意图；

图2为本发明基本流程图；

图3为本发明实施例砂岩的CT扫描图(一个断面)；

图4为本发明实施例砂岩的滤波平滑效果图(一个断面)；

图5为本发明实施例整体滤波后的效果图；

图6为本发明实施例提取孔隙效果图；

图7为本发明实施例骨架线示意图；

图8为本发明实施例节点示意图；

图9为本发明实施例节点局部放大图；

图10为本发明实施例骨架线上含有厚度信息的点(黑色点)；

图11为本发明实施例分割孔隙和孔喉的示意图；

图12为本发明实施例砂岩孔隙直径——数量示意图；

图13为本发明实施例砂岩孔隙直径——体积示意图；

图14为本发明实施例砂岩孔喉直径——数量示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细描述。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。

相反，本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。进一步，为了使公众对本发明有更好的了解，在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。

实施例1

本实施例提供一种基于三维岩心扫描图像的孔隙和孔喉识别的方法，孔隙和孔喉的示意图如图1所示，孔隙和孔喉识别的方法流程图如图2所示，主要包括以下步骤：三维图像获取、三维图像预处理、三维图像分割、孔隙孔喉的分割，通过连通区域内的孔隙空间结构特征，将孔隙和孔喉分割开，分别得出孔隙与孔喉的空间特性。

本实施例以CT扫描砂岩数据为例，包括如下步骤：

1.利用微米CT获得数据图像：

通过微米CT获取的图像大小为：550×550×800，每个像素代表3um。即扫描区域大小长为1.65mm，宽为1.65mm，高为2.4mm。图3为本发明实施例砂岩的CT扫描图。所用数据图像并不只是针对CT图像，也可以是通过聚焦粒子束‐扫描电镜获得的电镜图像等。图像大小不受限制。每个像素大小也不受限制。

2.图像预处理：

包括图像滤波，在保证原始图像所有重要细节信息特征情况下，将扫描成像设备得到的带有一定噪点的灰度图像进行滤波。我们采用非局部均值滤波(Non‐Local‐Means)，对图像进行平滑处理。图4为本发明实施例砂岩的滤波平滑图，可以看出通过该方法，图像在去除噪声的同时很好地保留图像的细节特征，整体效果示于图5。

图像滤波也可采用有均值滤波、非局部均值滤波、中值滤波、维纳滤波、高斯滤波等。

3.图像分割：

采用图像阈值分割将孔隙空间和矿物区分开来。阈值分割是一种传统的最常用的图像分割方法，因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。原理是：通过设定不同的特征阈值，把图像的象素值分成若干类。这种方法不仅可以极大的压缩数据量，也可以大大简化了分析和处理步骤，因此是进行图像分析之前必要的图像预处理过程。图像阈值化的目的是要按照灰度级，对像素集合进行划分，得到的每个子集形成一个与现实景物相对应的区域，各个区域内部具有一致的属性，而相邻区域不具有这种一致属性。

由于阈值分割不考虑孔隙的真实存在的可能性，因此，需要将异常点剔除，例如：分割出单个像素的点是否为孔隙，在此我们认为此情况下不是孔隙，需要将其剔除。

对图5采用阈值分割，分割得到孔隙、矿物的分布图，并将异常点1个像素的孔隙和矿物剔除。孔隙空间效果示于图6。

4.孔隙和孔喉的分割：

我们对提取出的孔隙空间图像进行距离变换，通过图像细化算法提取出相互连通区域的骨架，提取的骨架线示于图7。骨架线上的每个点是严格的一个体素宽度，且每一个点除了包含空间坐标信息：X、Y、Z值，还包括距离最近的非孔隙边界的距离(即厚度)。骨架线中的端点为孔隙空间开始和结束的点，骨架线中的交点为孔隙空间交叉的点。骨架线中的端点和交点称为节点。节点图示于图8，节点局部放大图示于图9。根据节点断开，将相互连通孔隙区域的骨架线断开成若干部分，一部分都包括若干含有厚度信息的点，如图10所示，将每一部分单独分析，寻找每一部分相对狭小的区域作为孔喉，相对宽大的部分为孔隙。

以其中相邻2个节点间的数据序列为例：

第一步，提取相邻节点间所有点的厚度，按照顺序将所有厚度存于一个数列中，并将相邻的厚度相同的点合并成一个点。

第二步，将厚度的数列进行极值化，数列中只保留极大值(比周围2个数都大的点)和极小值点(比周围2个数据都小的点)，其中端点只与相邻的一个点进行比较，所得的数据序列只包含极大值和极小值，且极大值和极小值交替出现且呈锯齿状；极值化后的图像示于图11，点A…K为筛选出的极值点，其厚度分别为V_a…V_k。

第三步，我们设定一个孔喉厚度与相邻孔隙最大厚度的比值a(0<a<1)，一旦计算的比值小于或者等于a值，即被认为是相对狭小区域边界点即疑似孔喉边界点。这个比值在分割孔隙和孔喉过程中起着非常重要的作用，各领域可以根据各领域应用标注设定a值。本文以a＝0.4为例，即，一旦出现计算的比值不高于a＝0.4，即该极小值点被认为是疑似孔喉边界点。

接下来分析极大值与极小值序列，引入两个虚拟值：V_min和V_max，将开始的极小值和极大值这两个值分别赋值给V_min和V_max，将V_min和V_max进行比较，看其比值是否大于a。①如大于a，再往下一个点进行查找，下一个点为极小值点，首先将其厚度与V_min值进行比较，如果该厚度比V_min值小，则V_min值等于更小的那个的厚度，然后再将V_min和V_max的比值与a值进行比对；如果该厚度不小于V_min值，则V_min值不变，再进行搜寻下一个极大值点，首先将搜寻的极大值点厚度与V_max值进行比较，如果该值大于V_max值，则V_max值等于更大的极大值点的厚度，之后再将V_min和V_max比值与a值进行比对；如果该厚度不大于V_max值，则V_max值不变，再进行搜寻下一个极小值点，以此类推。②如小于等于a值，说明该极小值点可能为孔喉边界点；将V_max值清空，进行搜寻下一个极大值点，以此类推。

所述孔喉确定具体为：针对在所述极值序列为极小值的所述疑似孔喉边界点d″_p，当存在d″_p/d″_o≤a且d″_p/d″_q≤a时则认定为孔喉边界点；其中o、p、q∈(1，f)，且o<p<q；

则，在所述厚度数列中，相应的在d″_o和d″_q之间的区域为孔喉。

近一步，为了方便直观的看出各部分的孔喉区域，对孔喉进行标注，具体实施步骤如下：根据11所示，将V_a赋给V_max，V_b赋给V_min。计算V_min和V_max的比值，发现其比值小于a＝0.4，则从B点至A点标以“→”，箭头的方向指向极大值A点。下一个搜寻的点为极大值C点，此时将V_max值清空，V_min值保留，将C点厚度赋于V_max，这时V_min＝V_b，V_max＝V_c，再进行计算V_min和V_max的比值，其比值也小于a＝0.4，则从B点至C点标以“→”，箭头的方向指向极大值C点。下一个搜寻的点为极小值D点，此时将V_min值清空，V_max值保留，将搜寻的下一个极小值D点厚度赋于V_min，计算V_min和V_max的比值，其比值大于a＝0.4。继续往下搜寻极大值E点，由于V_e大于V_c，所以V_max由V_c变为V_e，此时，V_min等于V_d，计算V_min和V_max的比值，其比值大于a＝0.4，继续搜寻下一个极小值F点，由于V_f小于V_d，V_min由V_d变为V_f，此时，V_max等于V_e，计算V_min和V_max的比值，其比值小于a＝0.4，则从F点至E点标以“→”，箭头的方向指向极大值E点。这时将V_max值清空，V_min值保留，将搜寻的下一个极大值点厚度V_g赋于V_max，再计算V_min和V_max的比值，发现比值大于a＝0.4，继续往下搜寻。以此类推，计算完全部厚度极值序列。最终出现的箭头为A←B、B→C、E←F、H→I和J→K。具体情况示于图11。通过分析相邻的箭头只会出现3种情况：①相交(“→”和“←”)，②背离(“←”和“→”)和③追逐(“→”和“→”或者“←”和“←”)，只有②背离情况，两个极小值点之间的部分为孔喉，①相交和③追逐情况直接合并为孔隙。由此可知B点和F-H部分为孔喉，A-B、B-F和H-K为3个孔隙，然后再将之前去掉的非极值点和重复的点按顺序加到极值序列中，按照箭头指示的极大值和极小值点进行标记，最终完成该序列中孔隙和孔喉分割。

将所有的相邻的节点中的数据进行统计，重复上面的步骤即得到最终的分割结果，最终得到孔隙和孔喉参数。

5.分析：

得到该数字岩心孔隙孔喉参数为：孔隙方面包括：孔隙体积、孔隙等效直径、孔隙平均直径以及孔隙包含的具体点的信息；孔喉方面包括：孔喉体积、孔喉等效直径、孔喉平均直径以及孔喉包含的具体点的信息。

本实施例分割出孔隙43959个，孔隙体积最大为3.6×10⁵μm³，最小为78.4μm³。孔隙数目呈现两头小中间大的趋势，如图12所示，且绝大多数孔隙等效直径分布于12μm至60μm间。直径对应的孔隙体积百分比随直径的增大而升高，在直径为60μm处孔隙体积百分比达到最大值，随后降低，降至直径为72μm处又转为升高，如图13所示。

本实施例分割出的孔喉个数为9638个，其中最小孔喉平均直径为3μm，最大孔喉平均直径为130μm，最小孔喉等效直径为3μm，最大孔喉等效直径为45μm。绝大多数孔喉等效直径在3至12μm之间，且孔喉直径在8μm左右时，个数最多，为2643个。如图14所示。

Claims (9)

1.一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，所述孔隙和孔喉识别方法为从三维岩心扫描图像中分割出孔隙空间，从所述孔隙空间对孔隙和孔喉进行分割实现对孔隙和孔喉的识别；所述孔喉为连通两个所述孔隙的通道；

从所述孔隙空间对孔隙和孔喉进行分割具体包括以下步骤：

骨架线提取：在所述孔隙空间提取出相互连通区域的骨架线；所述骨架线上的每个点是严格的一个体素宽度，每一个所述点包含空间坐标信息和厚度；

分段：将所述骨架线的端点和交叉点定义为节点，根据节点断开，将所述骨架线断开成若干线段，每一所述线段都包括多个所述点；

单线段分割：对其中一个所述线段进行孔隙和孔喉分割，所述单线段分割包括疑似孔喉点边界判断和孔喉确定；

全部分割：针对其他所述线段，采用所述单线段分割的方法对线段内的孔隙和孔喉进行分割；

所述空间坐标信息为X、Y、Z值；所述厚度为该点距离最近的非孔隙边界的距离；

所述疑似孔喉点判断为，当所述点的厚度与相邻区域的最大厚度的比值小于或者等于a，即被认为是相对狭小区域的边界点，所述相对狭小区域的边界点为疑似孔喉边界点；

所述孔喉确定为，在所述线段上，存在依次排列的点o、p、q、r，其中p点和q点重合或不重合，p点和q点为所述疑似孔喉边界点；o、p、q、r的厚度满足：d_p/d₀小于或者等于a；且d_q/d_r小于或者等于a；且p、q之间不存在点s使得d_p/d_s或d_q/d_s小于或者等于a；则，认定o、q两点之间为孔喉；所述0<a<1。

2.如权利要求1所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，所述疑似孔喉边界点判断具体包括以下步骤：

第一步，提取某一所述线段内所有所述点的厚度，按照点的顺序将所有厚度存于一个数列中得到厚度数列d＝{d₁,…,d_i,…d_n}，d_i为第i个点的厚度，i∈(1,n)，n为所述线段内点的个数；

若所述厚度数列中存在相邻的厚度相同时只保留一个厚度组成简化厚度数列d′＝{d′₁,…,d′_i,…d′_m}，1＜m≤n；

第二步，将所述简化厚度数列进行极值化，只保留所述简化厚度数列中的极大值和极小值得到极值序列d″＝{d″₁…,d″_i,…d″_f}，1＜f≤m；

所述极值序列只包含极大值和极小值，且极大值和极小值交替出现；所述极大值是指厚度比与其数列中相邻的两个点都大的点；

所述极小值是指厚度比与其数列中相邻的两个点都小的点；

其中端点只与相邻的一个点进行比较；

第三步，判断，分析所述极值序列，引入两个虚拟值：V_min和V_max，将所述极值序列中的第一个极小值点赋值给V_min，将所述极值序列中的第一个极大值点赋值给V_max，计算V_min/V_max；i∈[1，f]；

针对d″_i的判断为：若d″_i为极大值点，则比较d″_i与V_max，若d″_i>V_max，则重新赋值V_max＝d″_i；

若d″_i为极小值点，则比较d″_i与V_min，若d″_i<V_min，则重新赋值V_min＝d″_i；重新赋值后计算V_min/V_max；

当V_min/V_max大于a值时，进行d″_i+1的判断，d″_i+1的判断方法与d″_i的判断方法相同；

当V_min/V_max小于等于a值时，判定d″_i点为疑似孔喉点，同时相应的将V_min或V_max清空后进行d″_i+1的判断；

其中所述将V_min或V_max清空选择为：当d″_i为极大值时，将V_min清空；当d″_i为极小值时，将V_max清空；

依次对所述极值序列中所有点进行判断，将所述单线段中疑似孔喉边界点全部找出。

3.如权利要求2所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，所述孔喉确定中，

p点和q点仅选择在所述极值序列中的极小值点；

o、s、r点选择在所述极值序列中的极大值点。

4.如权利要求1所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，从三维岩心扫描图像中分割出孔隙空间具体为采用图像阈值分割法对三维岩心扫描图像的矿物区和孔隙空间进行分割。

5.如权利要求1所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，所述三维岩心扫描图像利用岩心扫描设备对岩心进行扫描得到。

6.如权利要求4所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，在采用图像阈值分割法前对所述三维岩心扫描图像进行图像滤波处理。

7.如权利要求1-6任一所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，完成孔隙和孔喉分割后进一步分析可得到孔隙参数和孔喉参数。

8.如权利要求7所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，所述孔隙参数包括孔隙体积、孔隙等效直径、孔隙平均直径以及孔隙包含的具体点的信息。

9.如权利要求7所述一种基于三维岩心扫描图像的孔隙和孔喉识别方法，其特征在于，所述孔喉参数包括孔喉体积、孔喉等效直径、孔喉平均直径以及孔喉包含的具体点的信息。