引 言

储层孔隙结构是指岩石所具有的孔隙和喉道的几何形状、大小、分布及其相互连通关系^[1-3]。由于油、气、水在储层岩石连通的孔隙中流动，孔隙结构特征对孔隙中流体的流动具有重要的影响，因此，岩石微观孔隙结构特征是研究流体渗流规律的基础，也是挖掘油气层产能及提高采收率的关键^[4-7]。

然而，涉及气体(如CO₂驱替实验等)以及三相流体的实验难以实施，实验结果在低饱和度区域并不可靠，也难以展示流体全方位的运移过程^[8]。且由于围压条件的限制，目前的实验研究不能准确地模拟现场地层条件。近年来，基于岩石微观图像的孔隙尺度建模被众多学者认可为该领域研究的一个突破口。该领域按照求解思路的不同主要分为两类：(1) 基于LBM方法的等效孔隙网络模型——用球体或圆柱体等规则几何体来表征多孔介质的微观孔隙特征，这种方法无法真实反映复杂的孔隙结构形状^[9-12]；(2) 基于有限元方法的网格模型——该方法有效地反映了真实多孔介质无序复杂的微观孔隙特征，但网格模型多为非结构化网格，网格质量难以控制，收敛性及求解精度较差^[13-16]。

本文通过对岩样三维CT图像进行处理分析，提出一种有效的孔隙度与孔径分布计算方法；针对现阶段孔隙尺度建模的缺陷提出结构化孔隙网格模型建模方案，并将计算结果与实验条件下测定的相关参数进行对比分析，验证了其可行性。

1 岩芯CT图像成像及图像后处理 1.1 微观CT成像过程简介

本文采用西南石油大学油气藏地质及开发工程国家重点实验室的蔡司Xradia MICROXCT-400开展岩芯三维微观成像工作。CT工作时沿z方向成像，每次成像1°并沿顺时针旋转180°。岩芯微观CT图像成像流程如下：

(1) 取样。为了确保获取高分辨率的岩样三维图像，扫描样品必须保持较小的尺寸和匀称的状态。从钻取的标准岩芯上用电钻取出尺寸较小的圆柱形岩样如图 1所示。

(2) 成像。将取出的岩芯样品固定在样品安置台上，关闭样品台外盖，开启成像程序，一般每个样品须扫描15~20 h。本文共选取3个岩样图像作为研究对象，其中岩样S1的三维岩芯孔隙图像三视图如图 2所示。

1.2 图像处理

从获取的图像中可以获得岩样孔隙度、孔径分布等一系列相应数据，并用来构建孔隙尺度的数值计算模型。但在进一步研究之前，需要对原始图像进行一系列前处理工作。

1.2.1 图像提取

如图 2所示，获取的岩样三维图像通常是圆柱形的，且基本上含有1 000³个像素。在开展图像降噪及二值化之前，从原始图像中提取100~500个立方像素的岩样图像作为模型重建数据(如图 2所示的红色框线部分)，这样做有以下优点：(1) 避免图像处理过程中图像的周围黑色区域出现分割错误。(2) 原始岩样为圆柱形，不利于后续孔隙网络的提取以及网格模型表面几何结构的重建，而正方体的岩样图像则不存在这样的问题。(3) 可以有效减少有限元网格数量，降低网格生成难度与计算时间。

1.2.2 降噪

现实中的图像都是带有噪声的，而噪声的存在对后续更高层次的图像处理产生不利影响。在区域数的确定上，图像噪声的存在容易过估区域的数目：在图像的非监督分割中，区域数的准确确定对分割性能产生重要的影响，如果估计图像的区域数过少，在分割中不同的区域不能很好地分离；如果估计图像的区域数过多，具有一致属性的区域可能被分割成许多不同的区域，并对随后的高层次图像处理产生不利影响，因此有必要对图像进行降噪^[17]。

采用中值滤波器进行图像降噪，可有效消除微观岩芯图像中的孤立小孔隙及岩石颗粒。将数字图像中某一点的像素值用该像素点领域中各点值的中值代换，使灰度值相差较大的像素点与其周围的像素值接近，消除了孤立的噪声点，滤除了图像的椒盐噪声^[17]。因此，本文采用中值滤波器实现图像降噪，既能去除图像中的噪声，又可以保护图像的边缘特征，降噪及图像复原效果较好。

1.2.3 图像二值化

图像的二值化处理就是将图像中的像素点灰度值设置为0或1，通过适当的阈值选取获得可以反映图像整体和局部特征(呈现出明显黑白效果)的图像，进而依据识别到的不同像素值实现图像中不同区域的分割。

对数字图像进行二值化处理，使像素点的性质只与其像素值为0或1的位置有关，不再涉及像素的多级值，简化了图像处理与识别过程，大大降低了数据的存储量。为了得到理想的二值图像，一般采用封闭、连通的边界定义不交叠的区域。所有灰度大于或等于阈值的像素被判定为特定物体，其灰度值以1表示，否则这些像素点被排除在物体区域以外，灰度值为0，表示背景或者例外的物体区域^[18]。图像的灰阶直方图及阈值选择见图 3。

图 4展示了原始图像及经降噪后的二值化图像，经中值算法处理后，图像中孤立的小像素数目有效减少，降低了孤立小孔隙及固体颗粒的数目。

2 岩芯孔隙度及孔径分布参数获取

经过图像处理后的岩样图如图 4d所示，图中黑色的部分为岩石孔隙，白色的部分为岩石骨架。从中可以直观地观察到孔隙尺寸及形状的多样性。基于二值化的岩样图像，实现了岩芯孔隙度及孔径分布的计算，得到其孔隙度、孔径分布等数据。

2.1 孔隙度计算

孔隙度指的是岩石中孔隙体积与岩石总体积的比值，它是定量描述岩石储集能力大小的重要参数。表达式如下^[19]

$ \phi = \dfrac{{{V_{\rm{p}}}}}{{{V_{\rm{b}}}}} \times 100\% = \dfrac{{{V_{\rm{p}}}}}{{{V_{\rm{p}}} + {V_{\rm{s}}}}} \times 100\% $

(1)

式中：$\phi$——孔隙度，%；$V_{\rm{b}}$——岩石的总体积，m³；$V_{\rm{p}}$——孔隙体积，m³；$V_{\rm{s}}$——岩石骨架体积，m³。

而在微观岩样图像中，孔隙体积与岩石骨架体积则对应不同颜色区域的像素体积。由于相同分辨率的图像中，每个像素所代表的尺寸是相同的，因此公式(1)可变化为

$ \phi = \dfrac{{{N_{\rm{p}}}}}{{{N_{\rm{b}}}}} \times 100\% = \dfrac{{{N_{\rm{p}}}}}{{{N_{\rm{p}}} + {N_{\rm{s}}}}} \times 100\% $

(2)

式中：$N_{\rm{b}}$——岩样图像总像素数目；$N_{\rm{p}}$——孔隙像素数目；$N_{\rm{s}}$——岩石骨架像素数目。

将得到的二值化图像转换为如图 5b所示的数组矩阵，以图像中每个像素的像素值来代表岩样图像。一个正方体的岩样图像对应一个三维数组矩阵，通过统计该矩阵中“0”和“1”的数目即可分别获得孔隙像素和岩石骨架像素数目。表 1即为岩样的分辨率、像素数目等图像数据及求得的孔隙度以及与原始样品实验测得的孔隙度对照表。由表中可知，图像得到的孔隙度数据与实验数据极为接近。由于图像仅为原始样品的一小部分，该误差在工程允许范围内。

2.2 孔径分布计算

孔径大小和孔径分布是多孔材料的重要结构参数，直接影响多孔材料的性能，因此对孔径及孔径分布的定量表征尤为重要。传统地，在实验室条件下通过气泡法、压汞法、气体透过法、气体吸附法、气体渗透法、液-液法、悬浮液过滤法获取多孔材料的孔径分布数据^[19]。本文基于数值化图像提出了一种基于图像的孔径分布求解方法，具体步骤如下：

(1) 将二值化图像转换为图 5所示数组矩阵。

(2) 将数组矩阵导入Matlab，识别每一幅二维图像中的不同孔隙数目，以及该局部连通的孔隙中所包含的像素数目，以矩阵的形式进行存储；对于三维图像而言，其在空间上等同于N张二维图像在断层扫描方向上的堆叠(N——原始图像所选取的层数)，因而完成N层二维图像的孔径分布统计即可获得三维图像的孔径分布。

(3) 根据等面积原则，某一包含k个分辨率为l的像素的孔隙等效孔径为

$ D = 2l\sqrt k /\sqrt \pi $

(3)

最终，根据得到的孔径数据绘制出孔径分布曲线。通过计算，各个岩芯的三维图像及孔径分布数据分别如图 6所示。

3 基于结构化网格模型的岩芯渗透率预测

采用结构化的网格单元来表征原始微观图像的特定像素(如孔隙或固体颗粒像素)是岩芯孔隙尺度结构化网格模型的构建基础。为此，需要将原始图像的像素数目减少以简化建模过程中的数据处理，降低网格模型中的网格数量。如S1图像中包含了400³个像素，将其分辨率减小1倍，即变为200³个像素，其像素数目减少至原来的1/8，由此构建包含201³个节点和200³个网格单元的节点与网格数据文件。通过查找原图像中的孔隙或固体颗粒的位置参数并提取相应位置的网格与节点参数，即可实现孔隙或固体骨架模型的重构。该结构化网格重构模型不仅能够真实反映原始图像的孔隙结构，而且结构化网格将极大地提高数值模拟的收敛速度与预测的精确性。重构得到的岩样S1结构化有限元模型见图 7 所示。

生成的结构网格模型没有几何实体与之关联，这将导致有限元分析过程中软件无法识别计算域边界，引起计算的不收敛。因此，需要对生成的网格模型进行几何重构并使之与网络模型相关联。同时，由于孔隙网络模型中，孤立的孔隙不参与流体流动，在流体分析中是无效的。因此，在数值分析前删除这些无效孔隙可以有效减少数值计算的工作量。图 8展示了岩样S1的由网格模型重构得到的几何模型与无效孔隙分布。

重构得到的岩样S1，S2和S3的结构化网格模型相关参数见表 2。

数值模拟工作利用Fluent软件完成，流体控制方程采用N-S方程，假设流体为层流，采用速度-压力耦合求解器以及SIMPLE压力修正项，采用默认的松弛因子。边界条件为分别沿x，y，z方向施加压力梯度，四周各面则为不流动边界。设置1×10^-6的绝对收敛因子，各模型均在50个迭代步以内收敛，充分表明了结构化网络模型在收敛速度上的优越性。迭代完成后可以得到压力出口边界的流量，通过达西定律即可求得岩样沿压力梯度方向的渗透率数据。通过数值计算，各岩样的渗透率数据见表 3($K_x$——数值模拟数据；$\overline{K_z}$——实验数据)。岩样S1的速度场图与压力场云图见图 9及图 10。

通过表 3可以看出，岩样结构化单元模型求得的渗透率数据与原始样品的实验渗透率数据误差在±15%以内。由于图像仅为原始样品的一小部分，且多孔介质流动的数学模型未考虑微孔道内强烈的液固作用，该误差在工程允许范围内被认为是合理的。计算结果与实验结果的相互吻合验证了本文建立的结构化孔隙网络模型的正确性与数值模拟的求解精度。

4 结 语

本文基于岩样三维CT图像及图像处理技术，提出了一种利用图像获取岩样孔隙度及孔径分布数据的方法。采用Matlab编程，以3个岩样为例，通过与原始实验测定数据的比对，验证了该方法的合理性与精度。进而，通过编程构建了能够真实反映原始岩样孔隙形态的结构化网格模型。该模型不仅有效克服了传统等效孔隙网络模型无法真实再现多孔介质内部结构的缺陷，更弥补了采用CT图像重构软件建立的非结构化网格模型网格质量差的不足。通过对比分析数值模拟与实验得到的岩样渗透率数据，验证了该建模方法的有效性。虽然本文采用石油工业的相关算例作为例证，但本文所提出的图像分析方法与建模方法以可广泛的应用于医学、材料等一系列与多孔介质相关的领域。