基于图像处理和Kozeny-Carman方程的砂岩储层渗透率预测*
2022-06-23潘卫国
冯 进 石 磊 管 耀 潘卫国 杨 清
(中海石油(中国)有限公司深圳分公司 广东深圳 518054)
渗透率指在一定压差下，岩石允许流体通过的能力，它是表征岩石本身传导液体能力的参数，该参数对油气储层的评价至关重要。目前储层渗透率评价的常用方法主要有两大类：第一类是实验室直接测量[1-3]，即利用渗透率仪测量柱塞岩心的渗透率。该方法的优点是准确度高，但需要大量的柱塞岩心进行测量。第二类是测井计算获取[4-5]，根据所采用的测井方法不同，可分为常规测井计算渗透率、核磁共振测井计算渗透率和偶极声波测井计算渗透率方法。该方法的优点是可以连续计算井筒方向的储层渗透率，但也需要一定数量的岩心气测孔隙度和渗透率数据用于方法建模。
当柱塞岩心不易获取(岩心易碎等原因)，岩心气测渗透率数据缺乏的情况下，储层的渗透率评价会比较困难。已有研究表明，岩石渗透率的大小与孔隙度和孔隙网络连通性等微观特性密切相关。徐艳玲[6]、连会青[7]等对SEM图像进行处理，获取了孔隙面积、孔隙周长、孔隙半径、分形维数等参数，并建立渗透率与孔隙面积(或孔隙周长)和分形维数的拟合关系。考虑到图像处理“盒计数法”获取分形维数的局限性，姜涛 等[8]尝试使用“小岛法”对砂岩SEM图像孔隙的面积和周长进行分析并获取分形维数，并开展了渗透率预测。目前，已有的基于地质图像的渗透率预测方法取得了一定的效果，但也存在一些有待改进的地方，例如粘连孔隙准确分割问题、孔隙连通性与渗流通道弯曲度问题、渗透率计算公式的岩石物理意义等。考虑到以上问题，本文拟利用铸体薄片、扫描电镜等岩石微观孔隙图像(通过小碎块岩心制备观察样品)，获取岩石微观参数，基于Kozeny-Carman方程来计算砂岩渗透率。
1 Kozeny-Carman方程介绍
20世纪20至50年代左右，Kozeny[9]和Carman[10-11]提出了如式(1)所示的半经验、半理论的多孔介质渗透率模型：
(1)
式(1)中：k为岩石渗透率，mD；μ为流体黏度，Pa·s；γ为单位转换因子；φ为岩石有效孔隙度，小数；S为岩石比表面积，1/cm；CK-C为地区经验系数。
在式(1)的基础上，大量学者针对不同特征的多孔介质以及假设条件，衍生推导出了一系列联系渗透率与孔隙度、颗粒半径、比表面、弯曲度等参数的公式，这类结构相似且广泛运用的公式就是著名的Kozeny-Carman方程。以下分别介绍几类具有代表意义的Kozeny-Carman方程。
1.1 基于毛管束模型的Kozeny-Carman方程
根据达西定律，通过岩石的流量Q可表示为式(2)。Kozeny-Carman方程假设岩石可以抽象为一系列相互平行的弯曲毛管束，流体在毛管束中流动(毛管束模型)。假设毛管的截面为圆形，半径为r，毛管中流体的流动可看作层流，单个毛管的流量Q可表示为式(3)。
(2)
式(2)中：Q为流量，m3/s；A为岩石截面积，m2；L为岩石长度，m；Δp为入口端和出口端的压力差，Pa。
(3)
式(3)中：r为毛管半径，m；l为弯曲毛管的实际长度，m。
假设毛管截面形状为圆形，则岩石孔隙度φ和比表面积S可分别表示为式(4)和式(5)：
(4)
(5)
式(4)、(5)中：τ为弯曲度(此处针对毛管)，τ=l/L，无量纲。
联立式(2)、(3)，可得：
(6)
将式(4)、(5)代入式(6)中，可得到基于毛管束模型的Kozeny-Carman方程：
(7)
1.2 基于球状颗粒模型的Kozeny-Carman方程
公式(7)为根据毛管束模型推导的Kozeny-Carman方程，但实际岩石的孔隙结构与毛管束模型还是有一定差距的，相对而言，使用球状颗粒模型来表征岩石孔隙结构会更准确一些。因此假设岩石中球状颗粒的数量为n，颗粒直径为D，则岩石中所有颗粒的体积和表面积分别为nπD3/6和nπD2，颗粒占岩石总体积的百分比为1-φ，岩石总体积可表示为nπD3/[6×(1-φ)]，岩石比表面积可表示为公式(8)：
(8)
式(8)中：D为颗粒直径，m；n为岩石中颗粒数量。
将式(8)代入式(7)中，可进一步得到球状颗粒模型的Kozeny-Carman方程：
(9)
1.3 其他Kozeny-Carman方程
对于黏土，房营光 等[12]采用质量比表面的方式来展示Kozeny-Carman方程(式10)：
(10)
式(10)中：Cs为孔隙形状参数；ρs为黏土密度，g/m3；Ss为黏土质量比表面，m2/g。
Mavko和Nur[10]提出了利用孔隙直径取代颗粒直径的Kozeny-Carman方程(式(11))：
(11)
式(11)中：D为颗粒直径，m。
Mavko等[13]引入孔隙连通性概念，在式(9)的基础上提出了修改后的Kozeny-Carman方程(式(12))。该公式中，新增加了渗流孔隙度下限φp，即低于该孔隙度的岩石渗透率接近0。
(12)
式(12)中：φp为渗流孔隙度下限值，小数。
2 基于图像处理的渗透率模型构建
前述几种Kozeny-Carman方程(式(1)、(7)、(9)～(12))中均使用了孔隙度φ，该参数主要表征渗流空间占岩石总体积的比例，该参数获取相对较容易，并且是渗透率最主要的影响因素之一。另外渗流通道的宽度对岩石渗透率也有非常重要的影响，有些公式采用比表面积S对其进行表征(式(1)、(7)、(10))，有些公式采用颗粒直径对其进行表征(式(9)、(12))，有些公式则采用孔隙直径对其进行表征(式(11))，到底哪一参数最适合？这取决于参数获取的难易程度以及Kozeny-Carman方程所应用岩石的特点。一般而言，获取准确的比表面S相对比较困难，需要开展气体吸附比表面积检测(BET)，且该方法主要适用于微孔和介孔发育的岩石[14-18]，对于宏孔(半径>0.5 μm)发育的岩石效果较差；颗粒直径的获取方法相对成熟且简单，可通过筛析法、沉降法、激光法或铸体薄片观察法等粒度分析方法获得[19-20]；孔隙直径的获取方法较多，常用的有压汞法、铸体薄片观察法和微米CT扫描法[21-22]。另外，为了表征孔隙网络的复杂性，大部分公式还使用了孔隙弯曲度τ(式(7)、(9)、(11)、(12))，一部分公式使用了孔隙形状参数Cs(式(10))。
本文旨在提出一种仅依赖岩石显微图像的渗透率计算模型，因此在综合考虑各种参数对渗透率的影响程度和各种参数获取难易程度的基础上，基于公式(9)提出如下形式的Kozeny-Carman方程用于岩石渗透率计算：
(13)
式(13)中：c为综合了单位转换和地区经验的系数；φa为显微图像面孔率，小数。
3 基于图像处理的岩石微观参数确定
在缺乏岩石柱塞样品的情况下，气体渗透率测试无法进行，但只需要极少量岩石碎块或岩屑即可进行岩石微观孔隙结构观察并拍摄照片。常用的岩石微观孔隙结构观察手段包括铸体薄片、扫描电镜、微(纳)米CT等。本文以铸体薄片图像为例，阐述岩石微观图像处理及式(13)中渗透率计算所需参数的获取方法。
3.1 面孔率确定
铸体薄片将液态染色树脂在真空加压条件下注入岩石孔隙中，待染色树脂固化后磨成薄片，可在透射偏光显微镜下清楚的观察岩石孔隙结构。利用各种图像处理方法，可以基于铸体薄片照片获取微观尺度下的岩石骨架和孔隙喉道信息。图1展示了岩石铸体薄片照片及二值化后孔隙网络图像，图1a中蓝色树脂填充部分代表孔隙和喉道，图1b为经过图像处理后二值化的孔隙网络图像，黑色部分代表孔隙和喉道。图1b本质上为一个二维数组(x,y)，x和y分别为图1b横向和纵向上的像素点个数，统计图1b中黑色像素点的个数n，则图1b的面孔率可表示为式(14)：
(14)
式(14)中：x为图像横向上的像素点个数；y为图像纵向上的像素点个数；n为图像中孔隙所占的像素点个数。
本文的图像处理过程中，会生成两种图像，第一种是直接通过阈值分割并二值化的孔隙网络图像(图1b)，第二种是在该图像基础上对孔喉进行加强(使用粘连颗粒分割算法)后的孔隙度网络图像(图1c)。对于第一种孔隙网络图像，由于铸体薄片图像的分辨率有限，直接采用阈值分割，会导致低于图像分辨率的微孔隙和喉道不计入面孔率中，因此获得的面孔率偏低。对于第二种孔隙网络图像，由于孔喉得到了增强，获得的面孔率大幅度提高，但也存在一个问题，低孔渗岩石的喉道部分会得到额外增强，使得相对误差增大。在对比分析这两种图像结果时，发现如下现象：①第一种面孔率与气测孔隙度呈较好的线性关系，明显优于第二种面孔率；②第一种面孔率基本反映了对渗透率有主要贡献的大孔隙；③第一种面孔率所经历的操作较少，确定性高，有利于后期将算法转化为自动化处理软件。因此综合考虑以后，论文选择使用第一种孔隙网络图像的面孔率。
图1 岩石铸体薄片照片及二值化后孔隙网络图像
3.2 颗粒直径确定
岩石颗粒经孔隙和喉道分割以后(图1c)，可分别对岩石颗粒进行编号(图2a)，并提取其颗粒直径。由于真实的岩石颗粒并非圆形，常用的颗粒直径主要有3种：长轴直径、短轴直径和平均直径。长轴直径指过颗粒重心(图2a中红色数字编号所在位置即为该颗粒的重心)的最长直径值，短轴直径指过颗粒重心的最短直径值，平均直径则是指过颗粒重心的多条直径平均值(图2a)。将图像识别得到的颗粒平均直径按照从大到小排序，选取前10%颗粒，计算其平均直径的平均值作为岩石颗粒直径值(图2b)。
3.3 孔隙弯曲度确定
1937年，Carman[23]最早提出孔隙弯曲度(水力弯曲度)的概念，用来表征多孔介质中流动路径的弯曲特征，此后孔隙弯曲度在多孔介质渗流及导电特性等领域被广泛应用。孔隙弯曲度的公式可简单表示为：
τ=l/L
(15)
图2 岩石颗粒平均直径提取与代表性颗粒直径确定示意图
式(15)中：τ为弯曲度(此处针对孔隙)，无量纲；l为弯曲的孔隙路径长度，m；L为岩石样品长度，m。
对于结构规则的多孔介质，孔隙弯曲度可以很容易通过孔隙几何形态计算得到。但天然岩石孔隙结构复杂，其孔隙弯曲度的获取就要困难得多。本文使用的方法是通过孔隙网络中轴提取并计算最短路径来获得孔隙弯曲度。
3.3.1孔隙网络中轴提取
1967年，Blum[24]提出著名的“火烧草场模型”，形象地阐述了中轴的概念：想象有一块由二维闭合曲线围成的草场，在这条二维闭合曲线上(草场边界)放置火种，同时点燃，火焰以相同的速度向图形内部的各个方向燃烧，火焰相遇时熄灭，火焰熄灭点的集合即为该草场的中轴(也称为骨架)。
常用的中轴提取方法包括Voronoi方法、演化方法、距离场函数方法、拓扑细化方法等几大类。本文采用的是Fouard等[25]在进行血管图像处理时提出的一种倒角距离变换方法。该方法从图像边缘向中心删除血管的像素点，直到无法继续删除像素点为止，最终形成血管的中轴。岩石孔隙网络与血管网络具有一定相似性，因此可以使用该倒角距离变换方法提取岩石孔隙网络的中轴。
图3为基于倒角距离变换方法的岩石孔隙中轴提取示意图。图3a为倒角距离变换所用的掩模。图3b为铸体薄片图像处理后的孔隙-颗粒二值图某局部区域，图中一个小方块代表一个像素点，白色小方块表示岩石骨架，黑色小方块表示孔隙，Ii和Ij分别表示x方向和y方向的图像尺寸。通过以下2个步骤的扫描来计算岩石孔隙网络的距离图：①从图像左上角到右下角，即坐标(1,1)到坐标(Ii,Ij)，利用正向掩模进行正向扫描；②从图像右下角到左上角，即坐标(Ii,Ij)到坐标(1,1)，利用反向掩模进行反向扫描。在正向掩模或反向掩模扫描完毕后，均需删除图中低于某一阈值的像素点，并且重新进行二值化，分别得到正向掩模或反向掩模处理后的图像。再将这两张图像取交集即可得到最终的孔隙网络中轴(图3c)。
图3 基于倒角距离变换的岩石孔隙网络中轴提取示意图
3.3.2最短路径确定
孔隙网络中轴提取以后，下一步是确定流体流动的最短路径(最短优势路径)。最短路径的定义如下：在网络图中任意两点之间的连接线路存在两条及以上的情况下，肯定存在一条总权值最小的线路，即这两点之间的最短路径。本文求取孔隙网络最短路径使用的是Dijkstra算法[26]，具有计算速度快、稳定性好、准确率较高等优势，目前广泛应用于交通、导航、电网、通信等领域。
Dijkstra算法使用类似广度优先搜索的方法解决赋权图的单源最短路径问题。对于已提取中轴的岩石铸体薄片孔隙网络图，可以用黑色线条表示孔隙网络中轴(图4a)，也可以用彩色线条表示中轴，不同颜色表征孔隙或喉道的半径相对值(图4b)。下面在图4b中选取某一局部区域(图4c)演示如何确定最短路径：将所有线条的交点抽象为节点集合(圆形，代表孔隙)，两节点间使用直线段连接(线条，代表喉道)，该直线段的距离用两节点之间孔隙网络中轴的像素点个数表示，节点集合(V)与带权路径集合(E)就构成一个带权有向连通图G=
Dijkstra算法最终可得到孔隙网络中轴图起点和终点之间考虑权重(孔隙喉道大小)的最短路径长度，即流体在孔隙网络中的最短优势通道的像素点个数l。对于固定视域的铸体薄片，其左右边界距离L已知。利用公式(15)即可求得孔隙弯曲度。
图4 孔隙网络中轴的最短路径确定示意图
4 应用实例
以珠江口盆地番禺气田珠江组作为研究对象，其储层包括岩屑长石砂岩、长石岩屑砂岩和长石石英砂岩等，孔隙度分布在8%～24%，渗透率分布在0.1～5 000 mD。选取珠江口盆地番禺气田珠江组12个不同深度位置的铸体薄片样品。
由于薄片图像表征了极小尺度的岩石微观结构，而柱塞气测渗透率则尺度相对较大，因此需要选择合适尺度的薄片图像来进行图像处理。本文在选择薄片图像的视域时，主要考虑了如下2个原则：①以选定视域大小拍摄照片时，不同位置拍摄得到照片的面孔率和颗粒直径差异较小，即该视域下的照片尽量不受岩石非均质性的影响；②以选定视域大小拍摄照片时，要求有足够的分辨率，孔隙、喉道、颗粒等信息要尽量清晰。由于原则①要求视域尽量大，而原则②要求视域尽量小，两者相互制约，因此需要尽量平衡原则①和原则②，以获取最佳效果。
本文对每个样品均拍摄了多张不同放大倍数的照片，最终选择放大倍数为50倍、视域大小为3.072 mm×3.072 mm的照片(该视域下兼顾了原则①和原则②)，用于图像处理和岩石微观参数的确定。12个砂岩样品的最短路径提取结果如图5所示。12个砂岩样品的铸体薄片照片经图像处理后，提取了面孔率、颗粒直径和孔隙弯曲度，具体参数见表1。
本文处理得到的面孔率虽然低于气测孔隙度，但与气测孔隙度与较好的线性关系，满足渗透率计算要求；岩心粒度分析表明样品最大颗粒直径主要分布在0.25～0.71 mm，与本文处理得到的前10%最大颗粒直径平均值(230.22 ～588.50 μm)相吻合。渗透率的计算分两步：首先，选择4个砂岩样品(1-15号样品，K=207.900 mD；1-13号样品，K=83.660 mD；3-7号样品，K=0.994 mD；3-18号样品，K=0.375 mD)，以气测渗透率作为标定，利用公式(13)得到这4个样品的地区经验系数c分别为3.53、2.58、4.22和2.57，取4个样品平均值3.23作为地区经验系数c的最终取值(当研究区暂时无气测孔隙度数据时，一般c可在1～5取值)，因此在珠江口盆地番禺气田珠江组，基于图像处理的Kozeny-Carman方程公式(13)可进一步表示为公式(16)；然后，利用剩下8个样品进行模型的检验(表1)。
图5 珠江口盆地番禺气田珠江组砂岩最短路径提取结果Fig .5 Shortest path extraction results of sandstone of Zhujiang Formation in Panyu gas field, Pearl River Mouth Basin
表1 珠江口盆地番禺气田珠江组砂岩微观参数与渗透率数据Table 1 Microscopic parameters and permeability data of sandstone of Zhujiang Formation in Panyu gas field, Pearl River Mouth Basin
(16)
计算结果表明，利用图像处理获取的微观参数(面孔率、颗粒直径、孔隙弯曲度)与根据式(16)计算得到的珠江口盆地番禺气田珠江组砂岩渗透率与气测渗透率吻合性较好(表1、图6)，表明本文提出的渗透率计算方法效果较好。
图6 珠江口盆地番禺气田珠江组砂岩气测渗透率与 计算渗透率交会图
Pearl River Mouth Basin
5 结论
论文提出了一种基于岩石显微图像(铸体薄片照片、扫描电镜图片等)获取岩石微观结构参数，并利用Kozeny-Carman方程计算岩石渗透率的方法。该方法适用于标准柱塞岩心缺乏的条件下，可以提供一定准确度的渗透率数据，为储层渗透率的评价提供了新思路和新方法。该方法亦可集成于各类显微镜系统配套的图像处理软件中，提供多孔介质样品图像拍摄与渗透率计算的智能一体化解决方案。
