CN113029899A - 一种基于显微图像处理的砂岩渗透率计算方法 - Google Patents

Abstract

本发明公开了一种基于显微图像处理的砂岩渗透率计算方法，包括获取目标砂岩的砂岩显微图像；对砂岩显微图像进行二值化分割处理获取二值化孔隙网络图像、孔喉加强后孔隙网络图像；根据二值化孔隙网络图像确定砂岩显微图像的面孔率；根据孔喉加强后孔隙网络图像确定颗粒直径；根据孔喉加强后孔隙网络图像确定孔隙弯曲度；根据上述确定的面孔率、颗粒直径、孔隙弯曲度计算得到砂岩渗透率。本发明针对难以大量获取柱塞岩心的情况下，仅通过容易获取的由小碎块岩心拍摄的岩石显微照片，通过图像处理提取砂岩显微参数并计算砂岩渗透率，弥补渗透率数据的缺乏，更好地指导砂岩油气藏的勘探与开发。

Description

一种基于显微图像处理的砂岩渗透率计算方法

技术领域

本发明涉及一种基于显微图像处理的砂岩渗透率计算方法，属于石油与天然气勘探与开发技术领域。

背景技术

渗透率指在一定压差下，岩石允许流体通过的能力，它是表征岩石本身传导液体能力的参数，该参数对油气储层的评价至关重要。目前岩石渗透率评价的常用方法主要有两大类：第一类是实验室直接测量，即利用渗透率仪测量柱塞岩心的渗透率，该方法的优点是准确度高，但需要大量的柱塞岩心进行测量。第二类是测井计算获取，根据所采用的测井方法不同，可分为常规测井计算渗透率、核磁共振测井计算渗透率等，该方法的优点是可以连续计算井筒方向的储层渗透率，但也需要大量的岩心气测孔隙度和渗透率数据用于孔渗关系拟合，并且还需要获取测井数据。

当柱塞岩心的获取存在困难(岩心易碎等原因)，岩心气测渗透率数据缺乏的情况下，岩石的渗透率评价会比较困难。已有研究表明，砂岩渗透率的大小与孔隙度和孔隙网络连通性等微观特性密切相关，本发明利用容易获得的铸体薄片、扫描电镜等显微图像(可通过小碎块岩心制备观察样品)，通过图像处理得到岩石微观参数，并用来计算砂岩渗透率。

发明内容

为了克服现有技术中的缺点，本发明提供一种基于显微图像处理的砂岩渗透率计算方法，本发明针对难以大量获取柱塞岩心的情况下，仅通过容易获取的由小碎块岩心拍摄的岩石显微照片，通过图像处理提取砂岩显微参数并计算砂岩渗透率，弥补渗透率数据的缺乏，更好地指导砂岩油气藏的勘探与开发。

本发明解决上述技术问题所提供的技术方案是：一种基于显微图像处理的砂岩渗透率计算方法，包括以下步骤：

获取目标砂岩的砂岩显微图像；

对砂岩显微图像进行二值化分割处理获取二值化孔隙网络图像、孔喉加强后孔隙网络图像；

根据二值化孔隙网络图像确定砂岩显微图像的面孔率；

根据孔喉加强后孔隙网络图像确定颗粒直径；

根据孔喉加强后孔隙网络图像确定孔隙弯曲度；

根据上述确定的面孔率、颗粒直径、孔隙弯曲度计算得到砂岩渗透率。

进一步的技术方案是，所述砂岩显微图像包括铸体薄片照片、扫描电镜照片。

进一步的技术方案是，根据二值化孔隙网络图像确定砂岩显微图像的面孔率中的计算公式为：

φ_a＝m/(L×S)

式中：φ_a为砂岩显微图像的面孔率，小数；L为二值化孔隙网络图像横向上的像素点个数，S为二值化孔隙网络图像纵向上的像素点个数，m为二值化孔隙网络图像中孔隙所占的像素点个数。

进一步的技术方案是，根据孔喉加强后孔隙网络图像确定颗粒直径中的计算公式为：

式中：n为将颗粒按平均直径由大至小排列后，排序前10％的颗粒数量，个；d_i为按平均直径由大至小排列的第i个颗粒的平均直径，μm；d为代表性颗粒直径，μm。

进一步的技术方案是，根据孔喉加强后孔隙网络图像确定孔隙弯曲度包括：

在孔喉加强后孔隙网络图像的基础上，利用倒角距离变换方法提取孔隙网络中轴；

根据孔隙网络中轴确定最短路径；

根据最短路径确定孔隙弯曲度。

进一步的技术方案是，根据最短路径确定孔隙弯曲度中的计算公式为：

τ＝l/L

式中：τ为孔隙弯曲度；L为二值化孔隙网络图像横向上的像素点个数；l为图像中孔隙网络最短路径的像素点个数。

进一步的技术方案是，根据上述确定的面孔率、颗粒直径、孔隙弯曲度计算得到砂岩渗透率中的计算公式为：

式中：φ_a为砂岩显微图像的面孔率，小数；c为地区经验系数；k为砂岩渗透率；τ为孔隙弯曲度；d为代表性颗粒直径，μm。

进一步的技术方案是，上述的地区经验系数可在1-5之间取经验值；或者通过岩心气测渗透率进行标定。

本发明具有以下有益效果：

1、本发明通过砂岩显微图像处理提取砂岩显微参数，并计算砂岩渗透率，其渗透率计算结果与气测渗透率吻合性好；

2、本发明与常用的实验室测量法和测井计算法相比，既不需要大量柱塞岩心，也不需要测井数据，只需要相对容易获取的小碎块砂岩样品(可以是岩心碎块或岩屑等)即可；

3、本发明可很好的弥补砂岩油气藏岩心气测渗透率数据的缺乏，有助于更好地指导砂岩油气藏的勘探与开发；

4、本发明转化为程序后，可集成于各类显微镜系统配套的图像处理软件中，提供多孔介质样品的图像拍摄与渗透率计算的智能一体化解决方案。

附图说明

图1为本发明的流程框图；

图2为砂岩铸体薄片照片及二值化后孔隙网络图像；

图3为孔喉加强后孔隙网络图像编号图；

图4为砂岩颗粒平均直径分布图；

图5为倒角距离变换方法提取二值化孔隙网络图像中轴图；

图6为孔隙网络中轴的最短路径确定图；

图7为P气田6个砂岩铸体薄片显微图像选取与二值化分割结果图；

图8为1-1号砂岩样品砂岩显微图像颗粒分割结果图；

图9为1-2号砂岩样品砂岩显微图像颗粒分割结果图；

图10为2-1号砂岩样品砂岩显微图像颗粒分割结果图；

图11为2-2号砂岩样品砂岩显微图像颗粒分割结果图；

图12为3-1号砂岩样品砂岩显微图像颗粒分割结果图；

图13为3-2号砂岩样品砂岩显微图像颗粒分割结果图；

图14为P气田6个砂岩砂岩样品的颗粒平均直径分布图；

图15为P气田6个砂岩显微图像的孔隙网络最短路径图；

图16为P气田6个砂岩气测渗透率与计算渗透率交会图。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明的一种基于显微图像处理的砂岩渗透率计算方法，包括以下步骤：

(1)砂岩显微图像拍摄与筛选；

砂岩显微图像包括铸体薄片照片、扫描电镜照片等，砂岩显微图像的选择应满足如下两个条件：①孔隙喉道清晰，②有相对较大的视域范围，能容纳砂岩具有代表性的孔隙网络。一般而言，这两个条件是相互制约的，孔隙喉道越清晰，需要放大倍数越大，但同时视域范围就会越小。因此需要选择具有合适视域的砂岩显微图像，一般以放大倍数为50倍左右，视域大小3mm×3mm左右的砂岩显微图像为宜。

(2)砂岩显微图像二值化分割；

砂岩显微图像中最为常见的是铸体薄片照片，其次是扫描电镜照片。铸体薄片将液态染色树脂在针孔加压条件下注入岩石孔隙中，待染色树脂固化后磨成薄片，可在透射偏光显微镜下清楚的观察岩石孔隙结构。将彩色铸体薄片照片转换为灰度图以后，选定合适的阈值，可以将砂岩显微图像分割成骨架和孔隙。图2展示了砂岩铸体薄片照片及二值化后的孔隙网络图像，图2中的(a)为步骤(1)中拍摄与筛选得到的原始铸体薄片照片，蓝色树脂填充孔隙和喉道；图2中的(b)为二值化孔隙网络图像，用于步骤(3)面孔率确定；图2中的(c)为孔喉加强后孔隙网络图像，用于步骤(4)颗粒直径确定和步骤(5)孔隙弯曲度确定。

(3)面孔率确定；

经过图像处理后的二值化孔隙网络图像中，黑色代表孔隙和喉道，二值化孔隙网络图像的本质为二维数组(x，y)，x和y分别为二值化孔隙网络图像横向和纵向上的像素点个数，统计二值化孔隙网络图像中黑色像素点的个数m，则砂岩显微图像的面孔率可表示为(公式1)：

φ_a＝m/(L×S) (1)

式中：φ_a为砂岩显微图像的面孔率，小数；L为二值化孔隙网络图像横向上的像素点个数，S为二值化孔隙网络图像纵向上的像素点个数，m为二值化孔隙网络图像中孔隙所占的像素点个数。

(4)代表性颗粒直径确定；

基于孔喉加强后孔隙网络图像(图2中的(c))，将砂岩颗粒通过孔隙和喉道分割开，可分别对每个砂岩颗粒进行编号(图3)，并提取每个颗粒的直径。由于真实的砂岩颗粒并非圆形，常用的颗粒直径主要有三种：长轴直径、短轴直径和平均直径。长轴直径指过颗粒重心(图3中数字编号所在位置即为该颗粒的重心)的最长直径值，短轴直径指过颗粒重心的最短直径值，平均直径则是指过颗粒重心的多条直径的平均值(图3)。将图像识别得到的颗粒平均直径按照从大到小排序，选取前10％的颗粒(也可选取10％或20％或其他比例的颗粒)，计算其平均直径的平均值作为砂岩代表性颗粒直径(公式2)(图4)。

式中：n为将颗粒按平均直径由大至小排列后，排序前10％的颗粒数量，个；d_i为按平均直径由大至小排列的第i个颗粒的平均直径，μm；d为代表性颗粒直径，μm。

(5)孔隙弯曲度确定；

孔隙弯曲度用来表征砂岩中流体流动路径的弯曲特征，其确定分3步进行：①孔隙网络中轴提取；②最短路径确定；③孔隙弯曲度计算。

①孔隙网络中轴提取；

在孔喉加强后的孔隙网络图像(图2中的(c))的基础上，利用倒角距离变换方法提取孔隙网络中轴(图5)。图5中的(a)为倒角距离变换所用的掩膜。图5中的(b)为铸体薄片图像处理后的孔隙-颗粒二值图某局部区域放大图，图中一个小方块代表一个像素点，白色小方块表示砂岩骨架，黑色小方块表示孔隙，I_i和I_j分别表示x方向和y方向的图像尺寸。通过以下2个步骤的扫描来计算砂岩孔隙网络的中轴图(图5中的(b))：1、从图像左上角到右下角，即坐标(1，1)到坐标(I_i，I_j)，利用正向掩膜进行正向扫描；2、从图像右下角到左上角，即坐标(I_i，I_j)到坐标(1，1)，利用反向掩膜进行反向扫描。正向掩膜和反向掩膜的具体计算方法如下：扫描前先将砂岩颗粒(白色像素点)的值初始化为0，孔隙网络(黑点像素点)的值初始化为某一非常高的值(例如254)。黑色方框所在像素点为当前需要计算并重新赋值的点，与之相邻的4个方框像素点为计算结果的数据来源，重新赋值后的黑色方框所在像素点数值等于相邻方框所在像素点的最小值与4个相邻方框所在像素点数值平均值的和(公式3)。在正向掩模或反向掩模扫描完毕后，均需删除图中低于某一阈值的像素点，并且重新进行二值化，分别得到正向掩模或反向掩模处理后的图像(公式4)。再将这两张图像取交集(公式5)即可得到孔隙网络中轴(图5中的(c))。

a_(i,j)＝min{ω_(i,j-1),ω_(i-1,j-1),ω_(i-1,j),ω_(i-1,j+1)}+(ω_(i,j-1)+ω_(i-1,j-1)+ω_(i-1,j)+ω_(i-1,j+1))/4 (3)

式中：min为求极小值符号；ω_(i,j)为原始图像的m行n列点阵中第i行j列像素点的权重值，1≤i≤m，1≤j≤n；a_(i,j)为掩模处理后图像的m行n列点阵中第i行j列像素点的权重值，1≤i≤m，1≤j≤n。

式中：σ为给定的阈值；β_(i,j)为掩模处理并二值化图像的m行n列点阵中第i行j列像素点的权重值，1≤i≤m，1≤j≤n。

式中：β_1(i,j)为正向掩模处理并二值化图像的m行n列点阵中第i行j列像素点的权重值，1≤i≤m，1≤j≤n；β_2(i,j)为反向掩模处理并二值化图像的m行n列点阵中第i行j列像素点的权重值，1≤i≤m，1≤j≤n；γ_(i,j)为最终结果图像的m行n列点阵中第i行j列像素点的权重值，1≤i≤m，1≤j≤n。

②最短路径确定；

提取得到的孔隙网络中轴，不仅可以用黑色线条表示(图6中的(a))，也可以用灰色线条表示(图6中的(b))。在图6中的(b)中选取某一局部区域得到图6中的(c)，基于图6中的(c)演示如何确定最短路径：将所有线条的交点抽象为节点集合(圆形，代表孔隙)，两节点间使用直线段连接(线条，代表喉道)，该直线段的距离用两节点之间孔隙网络中轴的像素点个数表示，节点集合(V)与带权路径集合(E)就构成一个带权有向连通图G＝<V，E>，图中节点不含数值，路径全为正值(图6中的(d))。

将节点集合V分为两组<S，U>，第一组S为已求出最短路径的节点集合，初始时节点集合S中只有一个源点V，以后每求得一条最短路径，就将其对应的节点加入到节点集合S中，直到全部节点都加入到节点集合S中。第二组U为未确定最短路径的节点集合，按最短路径长度的递增次序依次把第二组U的节点加入第一组节点集合S中。在加入的过程中，总保持从源点V到第一组S中各节点的最短路径长度不大于从源点V到第二组U中任何节点的最短路径长度。此外，每个顶点对应一个距离，第一组S中某节点的距离就是从源点V到此节点的最短路径长度，第二组U中某节点的距离，是从源点V到此节点且仅包括第一组S中的节点为中间节点的最短路径。

对图6中的(d)中的孔隙网络中轴带权有向连通图执行以上最短路径计算方法，选择左上角的节点为起点，右下角的节点为终点。则可以得到从起点到终点的最短路径，即图6中的(e)中节点和线条。将带权有向连通图还原到孔隙网络中轴图中，擦除最短路径外的其他孔隙网络，即可得到孔隙网络中轴的最短路径l(图6中的(f))。对于图6中的(a)或图6中的(b)所示的整个孔隙网络，从图像最左端的节点依次与图像最右端的节点组合成起点和终点，抽象带权连通图并执行最短路径算法，对每一次结果计算路径长度，其中长度最短的路径即为整个孔隙网络的最短路径l。

③孔隙弯曲度计算；

利用公式(6)计算孔隙弯曲度：

τ＝l/L (6)

式中：τ为孔隙弯曲度；L为二值化孔隙网络图像横向上的像素点个数；l为图像中孔隙网络最短路径的像素点个数。

(6)砂岩渗透率计算；

根据前述步骤确定的面孔率、代表性颗粒直径和孔隙弯曲度，利用公式(7)计算砂岩渗透率：

式中：φ_a为砂岩显微图像的面孔率，小数；c为地区经验系数；k为砂岩渗透率；τ为孔隙弯曲度；d为代表性颗粒直径，μm。

公式(7)中地区经验系数c的确定方法有两种：①在1-5之间取经验值；②通过少量岩心气测渗透率进行标定。

方法①在1-5之间取经验值：若有相邻油气田，可以参考其取值，或者直接在1-5之间取经验值，例如取值2.5。

方法②通过少量岩心气测渗透率进行标定：先对每一颗有气测渗透率的岩心利用公式(8)确定单颗岩心的经验系数c_i，再利用公式(9)求所有岩心经验系数的平均值c_i，即得到地区经验系数c。

式中，c_i为第i颗有气测渗透率岩心确定的经验系数，i为1～p之间的整数。

实施例

以P气田6个砂岩样品铸体薄片为例，说明本发明的具体实施方式。

(1)砂岩显微图像拍摄与筛选；

为P气田6个砂岩样品铸体薄片拍摄放大倍数为50倍，视域大小3mm×3mm的显微图像，所拍摄显微图像尽量代表砂岩样品铸体薄片的典型孔隙特征(图7中的(a)、(b)、(c)、(d)、(e)、(f))。

(2)砂岩显微图像二值化分割；

将6个砂岩样品彩色铸体薄片显微图像转换为灰度图以后，选定合适的阈值，将砂岩显微图像分割成骨架和孔隙。图7中的(g)、(h)、(i)、(j)、(k)、(l)为二值化孔隙网络图像，用于步骤(3)面孔率确定；图7中的(m)、(n)、(o)、(p)、(q)、(r)为孔喉加强的二值化孔隙网络图像，用于步骤(4)代表性颗粒直径确定和步骤(5)孔隙弯曲度确定。

(3)面孔率确定；

基于二值化孔隙网络图像(图7中的(g)、(h)、(i)、(j)、(k)、(l))，分别利用公式(1)计算6个砂岩样品的面孔率(表1)：

1-1号样品： φ_a＝202585/(2048×1536)＝0.0644 (10)

1-2号样品： φ_a＝225234/(2048×1536)＝0.0716 (11)

2-1号样品： φ_a＝149422/(2048×1536)＝0.0475 (12)

2-2号样品： φ_a＝167982/(2048×1536)＝0.0475 (13)

3-1号样品： φ_a＝82418/(2048×1536)＝0.0262 (14)

3-2号样品： φ_a＝66060/(2048×1536)＝0.0210 (15)

(4)代表性颗粒直径确定；

基于孔喉加强后二值化孔隙网络图像(图7中的(m)、(n)、(o)、(p)、(q)、(r))，进行颗粒分割(图8-图13)并统计每个颗粒的平均直径d_i；分别将每个砂岩样品颗粒平均直径按照从大到小排序，分别绘制每个砂岩样品的颗粒平均直径分布图(图14中的(g)、(h)、(i)、(j)、(k)、(l))，选取每个砂岩样品前10％的颗粒，利用公式(2)计算每个砂岩样品的代表性颗粒直径d(表1)。

(5)孔隙弯曲度确定；

基于孔喉加强后二值化孔隙网络图像(图7中的(m)、(n)、(o)、(p)、(q)、(r))，利用倒角距离变换方法提取得到孔隙网络中轴(图15中的(a)、(b)、(c)、(d)、(e)、(f))，获得孔隙网络最短路径l(图15中(g)、(h)、(i)、(j)、(k)、(l)中的曲线)，并利用公式(6)分别计算得到6个砂岩铸体薄片的孔隙弯曲度：

1-1号样品： τ＝2744/2048＝1.34 (16)

1-2号样品： τ＝2990/2048＝1.46 (17)

2-1号样品： τ＝3092/2048＝1.51 (18)

2-2号样品： τ＝2990/2048＝1.46 (19)

3-1号样品： τ＝5058/2048＝2.47 (20)

3-2号样品： τ＝5652/2048＝2.76 (21)

(6)砂岩渗透率计算

统计在步骤(1)至步骤(6)中所确定的6个砂岩样品的面孔率、代表性颗粒直径和孔隙弯曲度(表1)，利用公式(7)计算6个砂岩样品的渗透率(C取值3.2)：

1-1号样品：

1-2号样品：

2-1号样品：

2-2号样品：

3-1号样品：

3-2号样品：

6个砂岩样品的计算渗透率与气测渗透率非常接近(表1、图16)，说明本发明提出的这种基于显微图像处理的砂岩渗透率计算方法效果好。

表1P气田6个砂岩显微图像获得的微观参数与渗透率计算结果

以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (8)

1.一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，包括以下步骤：

获取目标砂岩的砂岩显微图像；

对砂岩显微图像进行二值化分割处理获取二值化孔隙网络图像、孔喉加强后孔隙网络图像；

根据二值化孔隙网络图像确定砂岩显微图像的面孔率；

根据孔喉加强后孔隙网络图像确定颗粒直径；

根据孔喉加强后孔隙网络图像确定孔隙弯曲度；

根据上述确定的面孔率、颗粒直径、孔隙弯曲度计算得到砂岩渗透率。

2.根据权利要求1所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，所述砂岩显微图像包括铸体薄片照片、扫描电镜照片。

3.根据权利要求2所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，根据二值化孔隙网络图像确定砂岩显微图像的面孔率中的计算公式为：

φ_a＝m/(L×S)

式中：φ_a为砂岩显微图像的面孔率，小数；L为二值化孔隙网络图像横向上的像素点个数，S为二值化孔隙网络图像纵向上的像素点个数，m为二值化孔隙网络图像中孔隙所占的像素点个数。

4.根据权利要求2所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，根据孔喉加强后孔隙网络图像确定颗粒直径中的计算公式为：

式中：n为将颗粒按平均直径由大至小排列后，排序前10％的颗粒数量，个；d_i为按平均直径由大至小排列的第i个颗粒的平均直径，μm；d为代表性颗粒直径，μm。

5.根据权利要求2所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，根据孔喉加强后孔隙网络图像确定孔隙弯曲度包括：

在孔喉加强后孔隙网络图像的基础上，利用倒角距离变换方法提取孔隙网络中轴；

根据孔隙网络中轴确定最短路径；

根据最短路径确定孔隙弯曲度。

6.根据权利要求5所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，根据最短路径确定孔隙弯曲度中的计算公式为：

τ＝l/L

式中：τ为孔隙弯曲度；L为二值化孔隙网络图像横向上的像素点个数；l为图像中孔隙网络最短路径的像素点个数。

7.根据权利要求6所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，根据上述确定的面孔率、颗粒直径、孔隙弯曲度计算得到砂岩渗透率中的计算公式为：

式中：φ_a为砂岩显微图像的面孔率，小数；c为地区经验系数；k为砂岩渗透率；τ为孔隙弯曲度；d为代表性颗粒直径，μm。

8.根据权利要求7所述的一种基于显微图像处理的砂岩渗透率计算方法，其特征在于，地区经验系数可在1-5之间取经验值；或者通过岩心气测渗透率进行标定。

Images