0 引言

数字岩心技术是基于扫描电镜图像或CT扫描图像，结合计算机技术对图像进行处理，重新构建三维数字化图像，研究岩心物理性质的一种新兴技术^[1]。近年来，随着科学技术的不断发展，X射线CT扫描与电子显微镜等技术在石油领域方面的应用越来越广泛。同时，基于CT扫描图像或扫描电镜图像的数字岩心技术也得到了极大的发展^[2-5]。目前，数字岩心的构建主要有物理实验与模拟重构两种方法。物理实验方法主要通过扫描获取的岩心二维图像构建岩心三维图像，如X射线CT扫描法^[6]、聚焦离子束扫描电镜等；模拟重构方法是在拥有较少的岩心图像时，利用数学方法对岩心图像进行处理，随机生成新的图像，进而完成岩心三维图像的构建，如高斯模拟法^[7]、模拟退火法^[8]和过程模拟法^[9]等。

由于数字岩心能够完整、真实地展现岩心的结构特征并用于渗流模拟，因此许多学者利用数字岩心技术对岩心结构特征等进行了研究。姚军等^[10]利用整合法建立了碳酸盐岩双孔隙网络模型，对碳酸盐岩不同尺度孔隙的渗流规律进行了研究。赵明等^[11]基于多孔介质的分形理论，对9块数字岩心孔隙度和渗透率等进行了预测，证明了分形理论在利用数字岩心预测岩心物性方面的准确性。Sun等^[12]利用二维扫描电镜图像重构数字岩心，构建了基于数字岩心技术的页岩气藏与致密气藏表观渗透率计算模型，并对表观渗透率的影响因素进行了分析。林承焰等^[13]利用CT扫描图像进行数字岩心三维建模，对岩心孔喉特征进行分析计算，并利用模型进行了微观渗流模拟。Golparvar等^[14]评价了基于CT扫描图像提取的孔隙网络模型(PNM)在研究多孔介质中流固耦合作用方面的优缺点，并分析了多相流建模面临的挑战。

渗透率作为岩心的一个重要物理参数，其与岩心的微观孔隙结构密切相关，可以直观地表现岩心的渗流能力。合理确定岩心的渗透率对岩心孔隙结构的研究至关重要。本文基于数字岩心技术，结合分形理论，通过求取岩心分形特征参数，构建等效分形介质模型，对岩心渗透率进行预测。

1 数字岩心图像处理及孔隙结构分析

CT扫描图像的本质是衰减系数成像。不同物质对X射线的吸收能力不同，锥形X射线从不同方向穿过岩心后得到不同的投影值，利用所得到的投影值求取各个体素的衰减系数，并将其转化为灰度值即可得到岩心CT扫描图像。本次研究选取两块砂岩岩心进行微米CT扫描测试，岩心参数如表 1所示。

① mD(毫达西)为非法定计量单位，1 mD=0.987×10^-3 μm²，下同。

在岩心CT扫描图像孔隙发育较好的区域截取600像素×600像素×600像素的立方体进行分析。对图像进行中值滤波、阈值分割处理从而得到二值化图像，如图 1所示。二值化图像中黑色部分为孔隙相，其值为1，白色部分为基质相，其值为0。从二值化图像中可以看出，两块岩心孔隙分布具有明显区别：岩心B1的孔隙分布较为分散，孔隙大小较为均匀；岩心B2的孔隙分布较为集中，且存在大孔隙。

岩心孔隙的连通性是评价岩心的重要指标，为更加直观地观察岩心的孔喉分布及其连通性，利用二维扫描图像构建三维数字岩心，利用最大球法提取岩心的孔隙网络球棍模型^[15]，如图 2所示。从图 2中可以发现，虽然岩心B1孔隙数目较多，但孔隙、喉道半径小，孔隙连通性差，因此其岩心渗透率较小；而岩心B2存在着较多尺寸较大的孔隙与喉道，连通性较好，因此渗透率较大。

图 3为两块岩心孔隙网络模型的喉道半径分布曲线、孔隙半径分布曲线以及孔喉配位数分布曲线。从图 3中可以看出：岩心B2的孔隙及喉道半径分布偏右，喉道半径分布峰值在3.2 μm左右，且存在较大半径的喉道，喉道半径可高达10 μm，而岩心B1的喉道半径分布峰值在2 μm左右，喉道半径均值小于岩心B2(图 3a)；岩心B2的平均孔隙半径大于岩心B1，且两块岩心的孔隙半径分布相差较大，岩心B2的孔隙半径集中分布在4~6 μm处，多数孔隙半径大于岩心B1最大孔隙半径(图 3b)；两块岩心的孔喉配位数分布也存在着一定的差异，岩心B1的孔喉配位数主要集中分布在2~6之间，而岩心B2的孔喉配位数分布较岩心B1宽，大于6的孔喉配位数仍占有很大比例(图 3c)。以上分析可以很好地解释岩心B1的渗透率小于岩心B2的渗透率的原因。孔喉半径越大，孔喉配位数越高，孔喉连通性越强，渗流阻力越小，岩心渗透率越大。

2 分形渗透率模型

假设砂岩孔隙由一系列迂曲毛管束组成，根据分形理论可知，毛管束的数目与毛管直径呈幂指数关系^[16]：

式中：λ为孔隙直径；N(≥λ)为孔隙直径大于等于λ的孔隙数目；λ_max为最大孔隙直径；D_f为分形维数。

假设孔隙分布为连续的，则孔隙直径介于λ与λ+dλ之间的孔隙数目为^[17]

假设毛管束的长度也具有分形特征，毛管的迂曲长度L_t与岩心特征长度L₀符合下列关系^[18]：

式中，D_t为迂曲度分形维数。D_t、L₀分别由以下公式求得^[18-19]：

式中：λ为平均孔隙直径；τ为平均迂曲度；φ为孔隙度。

Wang等^[20]结合分形理论与哈根泊肃叶方程，构建了毛管束模型的绝对渗透率计算模型：

式中：K_L为岩心绝对渗透率；A为岩心的横截面积，且A=L₀²。该模型表明岩心的渗透率受分形维数、迂曲度分形维数以及岩心最大孔隙直径的影响。模型中各岩心的分形维数与迂曲度分形维数是二维的，其值介于1~2之间。迂曲度分形维数作为指数，其对渗透率有较大的影响，随着迂曲度分形维数增加，岩心渗透率降低^[18]。由于岩心渗透率与最大孔隙直径呈幂指数关系，最大孔隙直径对渗透率影响显著。下面将介绍各模型参数的具体计算方法。

3 分形渗透率模型参数计算方法 3.1 分形维数

岩心分形维数可利用不同方法进行求取，本次研究采用盒计数法求取^[21-23]。盒计数法是利用不同尺寸的方形盒子对岩心图片进行覆盖(图 4)，统计获得不同尺寸盒子下含有孔隙相盒子的数目，在双对数图上所求斜率为分形维数：

式中，N(δ)为尺寸为δ的盒子含有孔隙相盒子的数目。

本文利用MATLAB对岩心的分形维数进行求取，流程如图 5所示。图 6为岩心B2单张图片分形维数的求取过程及所有图片分形维数的汇总图。如图 6b所示，每张图片计算得到的分形维数不同，取其平均值作为模型的分形维数。

3.2 迂曲度分形维数

迂曲度对岩心的渗透率具有一定影响，在分形渗透率模型中其以迂曲度分形维数的形式出现。由式(4)可知，为求取迂曲度分形维数，需计算岩心的平均迂曲度。目前求取岩心迂曲度的模型有积分法和孔隙度函数法等，每种方法在计算孔隙迂曲度时具有一定的适应性。而基于流线法取得的迂曲度能更为精确地反映岩心中连通孔径的迂曲程度。Al-Raoush等^[24]基于CT扫描图像，利用MATLAB编程，提出了一种迂曲度计算方法。该方法通过对孔隙相周围的像素点进行扫描，标记岩心各个孔隙最小直径(即喉道直径)的位置，如图 7所示。随后选择起始点，搜索连通路径，通过喉道的中点进行连接，将其视为流线，流线长度与岩心长度的比值即为该连通孔径的迂曲度。对岩心中所有连通孔径进行扫描，获取岩心所有连通孔径的迂曲度，求取平均值，将其视为岩心的平均迂曲度。本文中所选取的迂曲度是沿z方向进行搜索所确定的。

获得岩心的平均迂曲度后，将其与岩心的平均孔隙直径、岩心特征长度带入式(4)，即可得到迂曲度分形维数。

3.3 最大孔隙直径

在分形渗透率模型中，最大孔隙直径的指数项放大了岩心中最大孔隙直径对岩心渗透率的影响，因而需要慎重地选择岩心最大孔隙直径。利用Image J软件对岩心CT扫描图像进行处理，可以得到岩心的最大孔隙直径。图 8为确定岩心最大孔隙直径的流程图。首先，确定岩心扫描图片的像素值，根据岩心像素值的不同，划定岩心内孔隙的边界，并对每个孔隙进行编号，部分孔隙标记如图 9所示；然后，根据定义不同，Image J可计算各孔隙面积、周长、最大直径、最小直径等参数，选择需要计算的孔隙参数运行；最后，将计算好的数据导出到Excel表中。利用Image J计算的部分岩心孔隙参数如表 2所示。

由于岩心中孔隙的不规则性，在求取岩心孔隙直径时，可求得几种不同的等效直径。通过比较各种直径的计算结果，选取岩心孔隙的内切圆直径(孔隙等效最小直径)为岩心孔隙的等效直径预测渗透率效果较好。统计每张岩心CT图中各孔隙的直径，获取各扫描图片中的最大孔隙直径。考虑岩心孔隙的连通性质以及岩心孔隙分布的偶然性，求取所有图片的最大孔隙直径，求取其平均值并将其视为岩心的最大孔隙直径。

4 渗透率预测

为评估基于数字岩心技术的分形渗透率模型计算方法，对上述两块岩心进行渗透率计算，并将计算得到的结果与岩心渗透率、孔隙网络模型计算得到的渗透率进行对比，如表 3所示。从表 3可以看出，分形渗透率模型计算的岩心渗透率与孔隙网络模型计算的渗透率相差不大，表明通过构建数字岩心的等效分形介质模型继而预测岩心渗透率是可行的。分形渗透率模型计算的渗透率与岩心气测渗透率存在一定的差异，这是由于用于计算的数字岩心只是岩心一个区域，岩心内部存在一定的非均质性，因此数字岩心计算得的渗透率与岩心气测渗透率具有一定的差异。

为更加准确地验证分形渗透率模型在计算数字岩心渗透率方面的可行性，对帝国理工学院网站公开的9块砂岩数字岩心图像^[25]进行处理，求取分形渗透率模型参数并进行渗透率预测。9块数字岩心参数如表 4所示，其中渗透率为孔隙网络模型计算得到的渗透率。基于CT扫描图像计算得到的等效分形介质模型参数如表 5所示。

2块微米CT扫描构建的数字岩心与9块帝国理工学院网站公开的数字岩心的等效分形渗透率模型计算结果与孔隙网络计算渗透率对比结果如图 10所示。利用分形渗透率预测模型所计算的渗透率与岩心孔隙网络计算渗透率有着良好的相关性，相关系数高达0.971 5。但结果也表明，随着岩心渗透率的增加，分形渗透率模型计算的渗透率误差增大。表 5中岩心S8的计算渗透率远高于岩心的渗透率，对比其模型参数可以发现，Image J软件计算得到的岩心S8的最大孔隙直径为343.66 μm，最大孔隙直径计算值偏大是造成其预测渗透率值过高的原因。

5 结论

基于微米CT扫描测试得到的数字岩心分析孔隙结构特征，结合分形理论构建数字岩心的等效分形介质模型，并预测岩心渗透率，研究结果表明：

1) 利用CT扫描图像提取的孔隙网络模型能够直观地反映岩心内部的孔喉分布及孔隙的连通性。孔隙与喉道尺寸分布、孔喉配位数分布影响岩心渗透率。孔喉尺寸越大、分布越宽，且孔喉配位数越大、连通性越强，岩心渗透率越高。

2) 通过构建数字岩心的等效分形介质模型预测岩心渗透率是可行的。等效分形介质渗透率计算结果受分形维数、迂曲度分形维数与最大孔隙直径的影响。分形渗透率模型计算得到的渗透率与孔隙网络模型计算得到渗透率具有很好的相关性，可以用于岩心渗透率预测。

3) 最大孔隙直径对分形渗透率模型计算结果影响较大。随着岩心渗透率的增大，利用分形渗透率模型计算的数字岩心渗透率与孔隙网络模型计算的渗透率偏差增大，这可能是由于最大孔隙直径计算值偏大造成的。