近年来, CT扫描成像技术在地球科学、土壤科学、石油工程和许多其他领域的研究中发挥着越来越重要的作用^[1-6]。基于CT扫描成像的数字岩芯分析技术为岩芯物理性质的模拟分析提供了一种新的方法, 岩芯三维可视性是其主要特点之一。三维数字岩芯图像可以较准确和直观地反映岩石内部的结构特征和孔隙的空间展布, 基于数字岩芯图像的岩石物理性质模拟为岩芯的物性分析提供了一种新的方法。与传统的岩芯实验分析方法相比, 由于 CT扫描成像技术对样本的无损性以及数字岩芯图像的重复利用性特点^[7], 可以对岩芯进行数次模拟而不破坏岩芯的结构和性质, 从而避免岩芯经多次实验后结果的偏差, 也可减少对岩芯的浪费^[8]。

孔‒渗关系是岩芯分析中最重要的问题之一。传统的孔隙度和渗透率的测量计算建立在实验的基础上, 基于数字岩芯图像的分析为岩芯的孔‒渗关系分析提供了一种新的方法^[3,8-12]。基于分相的数字岩芯图像, 可以较准确地获取岩芯的孔隙度。随着岩芯图像分辨率的提升及岩芯处理方法的优化, 准确度还会进一步提高^[4,7]。在渗透率分析方面, 近十年发展起来的格子玻尔兹曼方法(lattice Boltz-mann method, LBM)作为探究微观尺度下复杂孔隙结构中流体渗流和热传导的有力工具, 具有简单高效、适合大规模并行运算等优势, 自诞生之日起就受到广泛关注。除在一般简单流动问题中得到成功应用外, LBM方法还可以应用于多相流、化学反扩散微尺度粒子悬浮流、多孔介质及磁体力学等领 域^[13-16]。LBM方法与数字岩芯分析的有机结合为岩芯的渗透率模拟提供了一种新的方法。目前, 国内大部分学者主要采用数学方法进行岩芯的建模, 并采用LBM方式进行模拟^[10-12], 但存在数字岩芯尺寸过小等问题。国外常见基于物理方法的岩芯建模^[9], 但存在建模后岩芯体数据过大, 无法整体模拟等问题^[17-19]。

本研究首先基于砂岩数字岩芯样本, 应用闵可夫斯基泛函和 K-means 聚类法, 对岩石进行精细岩石类别的划分; 然后基于分层后的岩芯图像, 在层次系统的基础上, 应用驱替模拟的方法, 模拟在不同驱替压力下湿相与非湿相(水油两相)的分布情况; 最后基于不同驱替压力下湿相与非湿相的分布情况, 结合格子玻尔兹曼方法, 对各个小层的渗透率进行模拟, 采用 3 种方法对岩芯渗透率进行粗化, 从而获得岩芯渗透率和孔‒渗关系方程。

1 基于三维数字岩芯的精细分层

1.1 岩芯样本简介

本研究使用的样本为某砂岩岩芯, 如图 1(a)所示。真实尺寸为30 mm×7.5 mm×7.5 mm, 垂向为30 mm, 横截面为7.5 mm×7.5 mm。经物理法成像后获得数字岩芯图像, 取其中一段作为分析对象, 尺寸为11.5 mm×7.5mm×7.5 mm。去除边界无效体元后, 体元尺寸为1760×1193×1514, 分辨率为5.6 μm (每个体元的边长), 该图像某一XZ平面上的图像切片如图 1(b)所示。分别将孔隙、泥质和骨架颗粒的体元孔隙度赋值为1, 0.5和0, 求得岩芯的总孔隙度为16.08%,岩芯真实照片和数字图像反映了岩芯垂向上较强的非均质性及分层性。

1.2 基于K-means聚类法的岩芯精细分层

闵可夫斯基泛函是对不同维度样本几何形态的表征。对特定的样本, 无论是三维还是二维^[17], 闵可夫斯基泛函可以很好地反映其空间上的拓扑结构和形态特征^[20]。对一个三维样本, 闵可夫斯基泛函包括4个函数变量^[21], 分别为体积函数(V)、表面积函数(S)、平均积分曲率函数(M)和高斯积分曲率函数(X), 如式(1)~(5)所示。

\[{{V}_{i}}=f[V(Y),\text{ }S(Y),\text{ }M(Y),\text{ }X(Y)] \ \ (1)\]

,

\[V(Y)=V, (2)\]

\[S(Y)=\int\limits_{\partial y}{\text{d}s}\text{ ,}\ \ (3)\]

\[M(Y)=\int\limits_{\partial y}{\frac{1}{2}}\left[\frac{1}{{{r}_{1}}(s)}+\frac{1}{{{r}_{2}}(s)} \right]\text{d}s\text{ },\ \ (4)\]

\[X(Y)=\int\limits_{\partial y}{\frac{1}{{{r}_{1}}(s)\cdot {{r}_{2}}(s)}}\ \text{d}s\text{ ,}\ \ (5)\]

其中, V_i为闵可夫斯基泛函, V(Y)为体积函数, S(Y)为表面积函数, M(Y)为平均积分曲率函数, X(Y)为高斯积分曲率函数, r₁(s)和r₂(s)分别为表面 ∂y上的最小和最大曲率。

岩芯的闵可夫斯基泛函响应可以综合表征岩芯的物性。通常情况下, 体积函数表示岩芯孔隙度的大小, 其值越大, 代表颗粒所占的体积越大, 从而孔隙体积越小。表面积函数与岩石的比表面积相关, 体积函数值越大, 表面积函数值越小, 岩石的比表面积越小。平均积分曲率函数值与岩石颗粒的圆度、球度、均一度等特征有关, 一般情况下, 对于凸集均为负值, 且其值越大, 表示岩石的磨圆度和球度越好。高斯积分曲率函数与孔隙结构的拓扑结构以及连通性有关, 一般情况下, 负值表示连通情况良好, 正值表示连通情况较差^[22], 且其值越小表示连通性越好。将岩芯的闵可夫斯基泛函进行线性标准化(式(6))^[21,23]后, 应用 K-means 聚类方法对岩芯进行分层处理。

\[V_{i}^{(N)}=\frac{{{V}_{i}}-\mu }{\sigma }\ \ (i=1,\ 2,\ ...),\ \ (6)\]

其中, V_i^(N)是线性标准化的闵可夫斯基泛函; i 为样 本数量; V_i是原始的闵可夫斯基泛函; μ为闵可夫斯基泛函中变量函数的平均值; σ为闵可夫斯基泛函中变量函数的标准差或协方差矩阵。

显示原图| 下载原图ZIP| 生成PPT

图1   岩芯照片(a)与数字岩芯的某个XZ平面图像(b)

Fig. 1   Photo of sandstone core sample (a) and XZ slice of part of that sample (b)

考虑到岩芯尺寸和局部效应等因素, 在保留其主要特征的基础上, 最小化岩芯图像边界区域和局部效应的影响。应用一个单位测量区域, 使尺寸为900×900×48 的体元沿着岩芯 Z 方向自上而下移动, 移动速率为Δ Z=1, 获得岩芯的局部闵可夫斯基泛函, 如图 2 所示, 曲线的波动幅度可以明显地反映岩芯在 Z 方向上较强的非均质性, 说明在 Z 方向上进行岩芯的精细分层十分必要。

K-means 聚类又称 Lloyd 算法, 是一种典型的基于质心距离的无监督或半监督型迭代聚类方法。K-means聚类是运用数据分类算法, 将数据样本(标量、向量或矩阵)划分成不同类别, 每个数据样本属于唯一确定的类, 属于硬聚类算法^[24]。K-means算法通常输入聚类样本数和聚类数据样本, 输出聚类样本的类别。

应用K-means聚类算法, 将数据样本分为k类, 对岩石类型进行划分。首先, 根据数据样本的平均值 μ₁⁽⁰⁾, μ₂⁽⁰⁾, …, μ_k⁽⁰⁾, 任意给定一系列的质心, 将各个数据样本依据式(7)归类:

\[V_{i}^{(t)}=\left\{ {{v}_{p}}\text{:}{{\left\| {{v}_{p}}-\mu _{i}^{(t)} \right\|}^{2}}\le {{\left\| {{v}_{p}}-\mu _{j}^{(t)} \right\|}^{2}}\forall 1\le j\le k \right\}\ \ , (7)\]

其中, V_i^(t)为第i类中所有的闵可夫斯基泛函, k为分类数量, v_p为闵可夫斯基泛函, μ_i^(t)和μ_j^(t)分别为第$i$类和第$j$类闵可夫斯基泛函第t次的平均值。然后, 重新计算各组质心并更新各组质心, 直至质心不再发生改变, 如式(8)^[25]所示:

\[\mu _{i}^{(t+1)}={{\left| V_{i}^{(t)} \right|}^{-1}}\sum\limits_{{{v}_{p}}\in V_{i}^{(t)}}{{{v}_{p}}},\ \ (8)\]

其中, μ_i^(t+1)为第$i$类闵可夫斯基泛函第 t+1 次的平均值, |V_i^(t)|为第i类在第t次迭代中的所有的闵可夫斯基泛函的数量。

基于岩芯垂直方向上的闵可夫斯基泛函响应值(图 2), 结合数字岩芯图像(图 1)和基于岩芯图像的碎屑岩颗粒结构特性图(图 3), 可将岩石划分为粗粒区、细粒区和过渡区3种类型, 获得岩芯的精细分层, 结果如表1和图4所示。从图4可以明显地看出, 岩芯的分层结果与岩芯图像之间呈现较好的匹配性, 说明基于 K-means 聚类算法完成的岩石精细分层具有较好的准确性。

显示原图| 下载原图ZIP| 生成PPT

图2   沿垂直地层面方向砂岩岩芯局部闵可夫斯基泛函

Fig. 2   Local Minkowski functionals along Z-direction of sandstone core sample

2 分层岩芯的渗透率模拟及粗化

2.1 岩芯渗透率模拟

渗透率描述多孔介质中流体的流动能力, 其值越高, 代表流体在多孔介质之越容易流动^[26]。渗透率的常用工程单位为达西(1 D≈10^-12 m²), 一般采用达西渗流定律来定量表征渗透率^[27]。

本文采用格子玻尔兹曼方法对岩芯的渗透率进行模拟, 将宏观的纳维斯托克斯方程和微观的离散波尔兹曼函数相结合, 从而模拟流体的流动, 计算岩层渗透率。本文采用的格子模型为 D3Q19^[13-14], 将自由边界作为边界处理条件(流体自由流动区域)。考虑到边界效应和局部效应对计算结果的影响, 采用的自由边界厚度为 15。在该计算条件下, 经变换后的小层 Z 坐标值和岩层厚度体元值如表2所示。各岩层的计算起始坐标上移 15 体元, 厚度增加 30 体元(顶底部各 15 体元)。由于岩芯尺寸和坐标限制, 第一层和最后一层变换的结果与其他层不同。

显示原图| 下载原图ZIP| 生成PPT

图3   碎屑颗粒的结构特性(改编自文献[25])

Fig. 3   Various textural properties of siliciclastic grains (after Ref. [25])

采用驱替模拟方法, 对在驱替过程中数字岩芯湿相与非湿相的两相分布进行模拟^[9,28-29]。赋予每一体元非湿相的最大的球形半径来表征其达到不同驱替压力时湿相与非湿相的分布情况, 同时假设所有非湿相所能到达的区域均是通过沿某一方向渗流而到达的, 因此通过设置不同的非湿相体元半径来模拟驱替过程中不同驱替压力下湿相与非湿相的分布情况, 岩芯某一XZ平面的切片如图5所示。

采用驱替模拟的方法, 结合基于 K-means 聚类法的岩石类型精细分层结果, 对每一小层在 3 个方向(X, Y, Z)上进行渗透率的数值模拟运算, 并绘制渗透率模拟结果与该层孔隙度的散点图。如图6所示, 不同的半径(r 值)对应不同的驱替压力, 渗透率与孔隙度之间较明显地正相关。同时, 也明显地呈现不同的分布状况: 粗粒区的孔渗性最好, 过渡区其次, 细粒区相对较差。模拟结果很好地印证了岩石类型划分结果的准确性和必要性。

显示原图| 下载原图ZIP| 生成PPT

图4   基于K-means聚类法的岩石类型划分

Fig. 4   K-means cluster result by overlapping the segmentation image

2.2 基于模拟结果的粗化

采用以下3种方法对岩芯渗透率模拟结果进行粗化, 从而建立基于岩芯分层的数据模拟结果与岩芯间的相关关系。1) 基于简单数学平均法的岩石渗透率粗化; 2) 基于各岩石类型孔‒渗关系曲线的渗透率粗化; 3) 去除薄层后, 基于各岩石类型孔‒ 渗关系曲线的渗透率粗化。

2.2.1 基于简单数学平均法的岩石渗透率粗化

采用加权平均的方式, 对各个小层的渗透率数据直接进行粗化。对水平方向采用算术平均法(式(9)), 对垂直方向采用调和平均法(式(10))。

\[{{k}_{\text{horizontal}\left| x,y \right.}}=\frac{\sum\limits_{i=1}^{N}{{{k}_{i}}{{T}_{i}}}}{{{T}_{t}}}\ , (9) \ \ \]

\[{{k}_{vertical\left| z \right.}}=\frac{{{T}_{t}}}{\sum\limits_{i=1}^{N}{\frac{{{T}_{i}}}{{{k}_{i}}}}}\ ,\ \ (10)\]

其中, N 为小层的数目, k_i为各小层的渗透率, T_i为小层的体元厚度或者真实厚度, T_t为整个岩芯的体元厚度或真实厚度。应用上述算法, 对岩芯渗透率进行粗化, 结果如表 3 示。可以看出, 水平方向上渗透率明显优于垂直方向, 并且随着驱替压力的增大, 湿相渗透率呈现明显的下降趋势。

显示原图| 下载原图ZIP| 生成PPT

图5   岩芯驱替模拟下的基于相划分图像对应不同半径下沿X方向非湿相驱替湿相过程中流体的分布状况(同一XZ平面切片)

Fig. 5   Simulation process with capillary drainage transform map based on segmented-images (corresponding XZ slice) with in-creasing the radius value of X-direction invading

2.2.2 基于各岩石类型孔‒渗关系曲线的渗透率粗化

由于\(k\propto \phi\)(k为渗透率, \(\phi\)为孔隙度), 在一定的 孔隙度约束条件下, 对各类别岩石的渗透率与孔隙度关系曲线进行回归, 得到回归方程k = a · exp(b ·$\phi$)的参数a和b, 结果如表4和图7所示。

从表4和图7可看出, 除垂向上渗透率偏差较大外, 水平方向上的孔‒渗关系回归结果呈现较好的相关性; 过渡区岩石类型的曲线回归结果明显比粗粒区和细粒区差。

显示原图| 下载原图ZIP| 生成PPT

图6   岩芯驱替模拟不同半径下岩芯各小层渗透率与孔隙度的相关性

Fig. 6   Correlation of permeability vs. porosity of increasing the radius simulated by capillary drainage transform

显示原图| 下载原图ZIP| 生成PPT

图7   所有小层孔‒渗关系回归曲线

Fig. 7   Cross-relationship of permeability vs. porosity of all layers

2.2.3 去除薄层后基于各岩石类型孔‒渗关系曲线的渗透率粗化

由于薄层局部效应过大, 其数据对曲线回归质量有较大影响, 因此我们除去薄层(体元厚度小于100的1, 3, 4, 13小层)数据后再进行曲线回归, 结果如表5和图8所示。

从表5和图8可以看出, 去除薄层后, 回归质量有较大的提升, 与水平方向相比, 垂向上的回归质量仍存在一定差距, 且粗粒区和细粒区岩石类型的曲线回归结果仍明显优于过渡区。

显示原图| 下载原图ZIP| 生成PPT

图8   去除薄层后的孔‒渗关系回归曲线

Fig. 8   Cross-relationship of permeability vs. porosity without the thin layers

对比2.2.2节和2.2.3节两种方法的曲线回归结果可以看出, 去除薄层后, 回归结果质量有较大幅度的提高, 尤其在垂向上。因此, 如果岩层划分得过细, 对岩石物性的模拟会产生一定程度的影响。

3 结论

本文基于砂岩数字岩芯图像, 采用闵可夫斯基泛函与 K-means 聚类算法, 对岩石类型进行精细划分, 并将岩芯精细层次系统与 LBM 方法相结合, 对岩芯渗透率进行模拟和粗化, 获得粗化后的渗透率参数和孔‒渗关系曲线方程, 得到以下结论。

1) 闵可夫斯基泛函可以较全面地反映岩芯的物理性质, 其函数变量的数值与岩芯孔隙度、比表面积、颗粒的磨圆程度、孔隙的联通情况等具有较强的相关性。应用K-means方法可以较准确地将岩芯样本划分成性质迥异的粗粒、细粒和过渡 3 种岩石类型。

2) 利用基于数字岩芯图像的驱替模拟变换, 可以较准确和直观地模拟不同驱替压力条件下湿相与非湿相在岩芯中的分布状况, 可以模拟岩芯水驱油的过程以及不同油水饱和度条件下的相分布, 为孔‒渗和相‒渗曲线的分析提供一种新的方式。

3) 基于数字岩芯图像, 结合 LBM 方法, 可以较准确地模拟岩芯渗透率。利用基于模拟结果的直接加权平均以及基于岩芯分层的孔‒渗关系回归等方法, 可以较好地对岩芯渗透率模拟结果进行粗化。获得的岩芯渗透率粗化结果和基于岩石类型的孔‒渗关系曲线, 可以为实验分析和现场生产提供新的依据。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子