各向异性多孔介质数字岩心模型重构通用算法

宋帅兵, 张通

1 安徽理工大学 深部煤矿采动响应与灾害防控国家重点实验室, 安徽淮南 232001 2 合肥综合性国家科学中心 能源研究院, 合肥 230031 3 中国矿业大学 深部岩土力学与地下工程国家重点实验室, 江苏徐州 221116

0 引言

岩石等多孔介质材料内部的孔隙结构特征直接影响着其自身宏观物理、力学及化学特性,如电导率、渗透率、弹性模量、泊松比、抗压强度及颗粒吸附力等(Ju et al.,2014;Liu et al.,2009;Song et al.,2019a,b;Zhao et al.,2007).为能够对岩石等多孔介质材料内部的孔隙结构特征进行全面表征,常用的特征参数主要包括孔隙率、连通度、迂曲度、孔径分布、孔喉比、朝向分布、配位数及有效路径长度等,上述参数分别从不同角度对孔隙结构的几何分布、拓扑结构及空间发育等特征进行了综合表征(Song et al.,2016,2015;宋帅兵,2020).

近年来,数字岩石物理已经成为一种表征多孔介质微观孔隙结构和预测其宏观物理性质的有效技术手段,并在各类地下工程领域具有十分重要的应用,如地下水流动传输、煤系共伴生油气资源协调开采及二氧化碳封存等(Song et al.,2023;刘向君等,2014;袁亮,2019).数字岩石物理技术主要包括岩石数字岩心的构建和岩石物理特性的数值计算模拟(朱伟和单蕊,2014),作为数字岩石物理工作流程中最为关键的一步,岩石数字岩心的重构精度将直接影响着后续岩石物性数值计算结果的准确性和可靠度.

现阶段,岩石数字岩心主要通过使用如下两种方法来构建(Dong,2007):物理扫描直接重构法和数值模拟等效合成法.前者主要通过借助扫描电子显微镜(SEM)、计算机断层扫描(CT)和聚焦离子束/电子束双束扫描电镜(FIB/SEM)等先进的扫描成像技术来直接获取多孔介质内部孔隙结构的数字图像,其优点在于能够真实呈现所扫描岩样内部的孔隙空间结构,然而该方法存在费用成本高、扫描耗时长等不足,且无法兼顾岩样的扫描视域与分辨率(Ju et al.,2017).为弥补物理扫描直接重构法不足,数值模拟等效合成法现已成为一类极具吸引力和发展潜力的补充替代方法,其主要包括高斯随机场法、模拟退火法及多点统计法等,该类方法只需要基于由有限图像中提取得到的孔隙结构统计信息,并通过使用各类统计数学优化算法便可完成等效数字岩心的合成重构(Zhou and Xiao,2018).在上述方法中,模拟退火算法凭借其较低的费用成本花费和较高的模型重构精度吸引了众多学者的特别关注.国内学者赵秀才等(2007)详细阐述了基于模拟退火算法建立数字岩心的理论方法,并通过实例运算验证了上述理论的适用性;邹孟飞等(2015)针对页岩孔隙分布特征,提出基于非常快速模拟退火法和双重区域的页岩数字岩心重构方法,使重构结果更加接近真实页岩岩心;莫修文等(2016)则针对传统模拟退火法数字岩心重构效率低、且计算结果存在较多孤立点等诸多缺陷,从而设计了一种补充优化方案,使得优化方案结果优于传统模拟退火法的建模结果;杨永飞等(2018)基于模拟退火法和马尔可夫链蒙特卡罗法分别构建了无机孔隙数字岩心和有机孔隙数字岩心,然后按照一定算法叠加两个数字岩心,建立了能同时描述页岩两类孔隙的数字岩心.除此之外,国外学者Kim和Pitsch(2009)、Politis等(2008)及Stenzel等(2013)也分别采用模拟退火法或模拟退火法与其他方法相组合的方式实现了各类多孔介质数字岩心的重构,并通过与物理实验结果的对比验证了重构模型的准确性.

需要特别指出的是,绝大多数学者在应用上述传统模拟退火算法进行数字岩心重构的过程中,只考虑了不同类型孔隙结构统计表征函数对重构结果精度的影响,并未进一步考虑各类孔隙结构特征在不同方向上的差异性,也即假设岩石内部孔隙结构的所有空间特征均满足均匀且各向同性分布.赵秀才和姚军(2007)则通过利用格子Boltzmann方法计算数字岩心在不同方向上的渗透率,证实了传统模拟退火算法所建数字岩心具有良好的各向同性.然而,岩石等绝大多数天然形成的多孔介质材料内部结构事实上并不是各向同性的,其往往呈现出较明显的各向异性特征(Bhandari et al.,2015;Kwon et al.,2004;李滔等,2019;王沫然和王梓岩,2018),且该特征对于岩石的宏观物理力学性能具有十分显著的影响,如岩石的渗透特性(Wang et al.,2016;隋微波等,2021;王超等,2021;王磊等,2021).传统模拟退火算法的各向同性假设诚然可降低数字岩心重构所需的计算量,提高了重构效率,但其同时也导致数字岩心模型重构精度的降低,并对后续岩石物性数值计算结果的准确性产生不利影响.因此,为能够构建更加准确可靠的数字岩心模型,务必需将孔隙结构的各向异性特征融入到数字岩心的重构过程之中,而目前对于各向异性数字岩心重构的研究还相对较少.

为此,本研究拟在现有传统模拟退火算法的基础之上,通过对孔隙结构在不同方向上的分布特征分别进行独立的提取计算表征,同时结合先前已提出的系统增量计算更新策略(Song,2019)和交换像素选择策略(Tang et al.,2009;Zhao et al.,2007)等计算加速技术,提出一种能够用于各向异性多孔介质数字岩心模型重构的改进模拟退火算法,并选用具有显著各向异性特征的人工合成二维切片图像及真实碳酸盐岩和砂岩三维CT扫描图像作为测试样本,对所提出算法的准确性和可靠度进行验证.

1 改进模拟退火算法重构理论

模拟退火算法作为一种解决组合优化问题的近似算法,利用其进行岩石数字岩心重构的常规工作流程如图1所示(宋帅兵,2020).

图1 模拟退火算法数字岩心重构流程Fig.1 Flow chart of digital core reconstruction based on simulated annealing algorithm

1.1 数学统计函数选择

在模拟退火算法中,主要通过使用各类数学统计函数来对参考图像及待重构系统中的孔隙结构进行全面、准确且定量化的表征和描述,故统计函数类型及数量的选择对于数字岩心的重构精度及效率具有重要影响.在众多孔隙结构表征统计函数中,单点概率函数、两点概率函数和线性路径函数是最基本,同时也是应用最为广泛的三类统计函数.因此,本研究也将选用上述三种不同类型的统计函数,来完成对孔隙结构特征的定量化描述和表征,三者的定义分别如下:

φ=〈I(x)〉,

(1)

S(r)=〈I(x)I(x+r)〉,

(2)

(3)

式中,〈〉表示计算其内部值的平均值;I(x)为图像或重构系统内部x位置处像素点的属性值,一般情况下孔隙相取为1,基质相取为0;r为图像内部任意两像素点间的统计距离.

由上述定义可以看出:单点概率函数可以理解为任选一个像素且其属性为孔隙相的概率,该函数直接反映了孔隙相的占比,也即代表了岩样的孔隙度φ,其不会随着像素统计距离的改变而改变;两点概率函数可以理解为任选两个像素且二者属性均为孔隙相的概率,该函数反映了图像内部孔隙相分布的相关性信息,随着像素统计距离的增加,其会逐渐收敛于φ2;线性路径函数则可以理解为任选一条r长度线段且组成该线段的所有像素均为孔隙相的概率,该函数反映了图像内部孔隙相的连通性信息,随着像素统计距离的增加,其会逐渐减小为0(注意:与零值所对应的像素统计距离越大,表明孔隙的连通性越好).

1.2 孔隙结构各向异性特征表征

在利用上述统计函数进行孔隙结构特征信息提取的过程中,为了简化计算工作流程,传统模拟退火算法通常选择沿着正交坐标轴方向来进行统计函数的计算,也即二维数字岩心重构选择水平H和竖直V两正交方向,三维数字岩心重构则选择水平H、竖直V和纵深D三正交方向.每一类统计函数均需在所选择的方向上首先进行相应的计算,然后再进行整合汇总以得到最终的孔隙结构信息,具体计算方式如下:

(4)

(5)

(6)

式中,Ii(x)为图像或重构系统i计算方向上x位置处像素点的属性值;n为所选择的不同计算方向的数目,传统模拟退火算法中二维重构n=2,三维重构n=3.上述处理操作一方面可大幅度降低后续重构系统的优化难度,从而提高了岩心的重构效率,但另一方面其也极大程度的消除了孔隙结构特征不同方向上的差异性,进而无法对真实孔隙空间结构的各向异性特征进行准确表征.

不同于传统模拟退火算法采用式(4)、式(5)和式(6)所引起的孔隙结构各向异性特征消除,本研究将利用式(1)、式(2)和式(3)分别独立完成不同方向上孔隙结构特征的计算提取,从而实现对孔隙结构各向异性特征的表征.除此之外,为进一步提高孔隙结构特征的表征精度,在现有正交坐标轴方向的基础之上,本研究还将增加计算更多其他方向上的孔隙结构特征,以能够适用于更加复杂孔隙结构特征的准确表征,具体内容如下:对于二维数字岩心,将增加计算两正交坐标轴夹角平分线方向上的孔隙结构特征;对于三维数字岩心,将增加计算任意两正交坐标轴夹角平分线以及与三正交坐标轴具有相同夹角方向上的孔隙结构特征,详细计算方向如图2所示(红色指示为传统模拟退火算法所选择的正交坐标轴方向,绿色指示为本研究所补充增加的计算方向).

图2 统计函数计算方向示意图(a) 二维数字岩心; (b) 三维数字岩心.Fig.2 Schematic diagram of calculation direction of statistical function(a) 2D digital core; (b) 3D digital core.

由图2可以看出,二维数字岩心构建共选择了4条不同计算方向,三维数字岩心构建共选择了13条不同计算方向,两者所选定的方向几乎涵盖了孔隙结构特征朝向分布所有可能的情况.通过对上述不同方向上孔隙结构特征的提取,能够做到对孔隙结构各向异性特征的准确表征.

1.3 重构系统能量计算及更新

为与真实物理退火过程中固体自身的能量相对应,岩心重构系统的“能量”取其自身与所参考图像间的所有不同方向上不同类型统计函数值差值平方之和,其具体定义如下:

(7)

式(7)作为模拟退火算法优化的“能量”目标函数,其通过持续不断进行孔隙相与基质相位置的交换操作,并基于Metropolis准则来判断每次位置交换是否接受,进而完成自身内部结构的逐步更新,在此循环操作过程中系统的“能量”值也随之同步降低.当重构系统的能量值低于预先设定的极小临界值时,重构过程结束.此刻,我们认为重构岩心与参考图像在所选用的三个不同类型的统计函数上相互等效.事实上,对于孔隙结构的其他主要特征,如连通度及渗透率等,二者同样具有较高的一致性,也即二者具有统计学意义上较大程度相互等效的孔隙结构特征,下一节中将通过详细对比分析三维数字岩心的连通度及其渗透性等参数来对上述结论进行论证.

2 各向异性数字岩心重构

事实上,任何一种多孔介质材料内部的孔隙结构分布均呈现出一定程度的各向异性,只是其差异程度存在不同,并不存在完全的各向同性分布.为了能够使该项特征更加的显著和突出,从而以更好的验证和展示本文所提出改进模拟退火算法的有效性和准确性,本研究首先采用由人工合成方式得到的具有显著各向异性特征的二维切片图像作为参考图像,分别利用传统算法及上述所提出的改进算法重构其二维各向异性数字岩心模型.此外,为进一步检验改进模拟退火算法的可靠性和扩展应用性,本研究再次以具有各向异性特征的真实碳酸盐岩和砂岩的三维CT扫描图像作为参考图像,按照同样的方式分别重构其三维各向异性数字岩心模型.

2.1 二维各向异性数字岩心

本研究采用由人工合成方式得到的切片图像作为二维各向异性数字岩心重构的参考图像,如图3所示,其像素尺寸为300×300,白色像素代表孔隙相,黑色像素代表基质相,孔隙度为28.07%.可以明显看出,参考图像的孔隙结构形态多为狭长形,且其朝向以水平方向和负对角线斜向分布为主,呈现出了较显著的各向异性特征.

图3 二维各向异性数字岩心重构二值切片参考图像Fig.3 Binary slice reference image for 2D anisotropic digital core reconstruction

2.1.1 初始二维重构系统生成

设定初始重构系统的尺寸与参考图像保持一致,即同样为300×300,其孔隙度也与参考图像相同,由此可计算得到系统内部孔隙相的像素点数目.随后,在初始重构系统内部采用随机算法选择相同数目且不重复的位置,并设定当前所选择位置处像素点的属性为孔隙,而剩余未被选择位置处的像素点的属性则全部设定为基质,由此最终完成初始重构系统的生成,具体生成结果如图4所示,白色像素代表孔隙相,黑色像素代

图4 初始二维重构系统Fig.4 Initial 2D reconstruction system

表基质相.可以十分清楚的看出,初始重构系统内部的孔隙相分布处于一种相当杂乱无章的混沌状态,其所包含的各类孔隙结构特征无法得到准确的识别和提取.

为了能够更好的对比和展示后续两种不同算法重构结果的差异,本研究令两者均基于同一个起点,也即使用相同的初始重构系统来分别完成其对应类型岩心的构建.

2.1.2 传统模拟退火算法二维岩心重构

基于图1所示的模拟退火算法数字岩心重构流程,首先沿着预先选定的4条不同方向,利用式(5)和式(6)分别对参考图像和初始重构系统完成其两点概率函数和线性路径函数值的计算,得到其各自与方向无关的单一统计函数值,其中像素点的最大统计距离R=50.参考图像和初始重构系统两者统计函数值随像素点统计距离的变化曲线分别如图5a和图5b所示.

图5 二维重构统计函数曲线(a) 参考图像; (b) 初始重构系统.Fig.5 Statistical function curves of 2D reconstruction(a) Reference image; (b) Initial reconstruction system.

可以明显看出,参考图像与初始重构系统的统计函数曲线形态具有十分显著的差异,这与先前对比图3和图4所示的二者图像所得到的直观视觉感受相一致:前者两点概率函数值在像素统计距离为30时趋于稳定,线性路径函数值则在像素统计距离为40时降为0,表明本研究像素点最大统计距离R=50的选择是合理的;而后者两点概率函数和线性路径函数则分别在极小的像素统计距离时(r=2和r=5)便已经迅速达到稳定和降为0,再次表明初始重构系统内部孔隙结构的杂乱无序性.

随后,通过持续不断进行孔隙相与基质相位置的交换来实现重构系统的更新.由于传统模拟退火算法岩心重构所需的计算量相对较少,故本研究对其设定一个较低的能量临界值1×10-6,即当重构系统能量小于上述临界值时,重构过程结束.传统模拟退火算法所获得的孔隙结构最终重构结果如图6a所示,此外图6b还给出了其相应的统计函数曲线,可以看出其与图5a所示的参考图像统计函数曲线近乎一致,表明二者具有统计学意义上相互等效的孔隙结构特征.

图6 传统模拟退火算法二维数字岩心重构(a) 孔隙结构最终重构结果; (b) 最终统计函数曲线.Fig.6 Reconstruction of 2D digital core based on conventional simulated annealing algorithm(a) Final reconstruction result of pore structure; (b) Final statistical function curve.

2.1.3 改进模拟退火算法二维岩心重构

改进模拟退火算法岩心重构采用与传统模拟退火算法相类似的操作流程,即同样也首先分别对参考图像和初始重构系统完成其在4条不同方向上两点概率函数和线性路径函数值的计算,其中像素点的最大统计距离保持不变R=50.不同之处在于,改进模拟退火算法不再使用式(5)和式(6)将不同方向上的统计函数值进行整合汇总平均,而是采用式(2)和式(3)直接计算每一方向上的统计函数值,并将其作为一个独立的孔隙结构特征进行考虑.参考图像和初始重构系统两者不同方向上统计函数值随像素点统计距离的变化曲线分别如图7和图8所示.

图8 初始重构系统不同方向统计函数曲线(a) 两点概率函数; (b) 线性路径函数.Fig.8 Statistical function curves in different directions of the reference image(a) Two-point probability function; (b) Lineal-path function.

可以清晰看出:参考图像不同方向上孔隙结构的均质性和连通性均具有明显的差异,也即呈现出了较显著的各向异性特征,而传统模拟退火算法提取得到的如图5a所示的统计函数计算结果则未能对此差异进行准确的捕获和表征;对于初始重构系统而言,不同方向上孔隙结构特征近乎不存在任何差异,这主要是由于该系统是基于随机算法生成的均质系统.

随后,同样通过持续不断地进行孔隙相与基质相位置的交换来实现重构系统的更新,而由于此时重构所使用的统计函数数量较多(共计为8个),其导致每次系统更新所需的计算量相对较大,故本研究对其设定一个相对较高的能量临界值1×10-4,当重构系统能量小于该临界值时,重构过程结束.改进模拟退火算法所获得的孔隙结构最终重构结果如图9a所示,此外图9和图9c还分别给出了其不同方向上两点概率函数和线性路径函数曲线,可以看出其与图7a和图7b所示的参考图像统计函数曲线近乎一致,同样表明重构图像与参考图像具有统计学意义上相互等效的孔隙结构特征.

图9 改进模拟退火算法二维数字岩心重构(a) 孔隙结构最终重构结果; (b) 最终两点概率函数曲线; (c) 最终线性路径函数曲线.Fig.9 Reconstruction of 2D digital core based on improved simulated annealing algorithm(a) Final reconstruction result of pore structure; (b) Final two-point probability function curve; (c) Final lineal-path function curve.

2.2 三维各向异性数字岩心

对于三维各向异性数字岩心重构,本研究选用具有各向异性特征的真实碳酸盐岩和砂岩的三维CT扫描图像作为参考图像,其三维孔隙结构的渲染结果如图10所示.

需要特别说明的是,上述两者的三维CT扫描图像均取自于帝国理工学院孔隙尺度模拟及成像研究团队的公开出版数据(其可通过如下网站https:∥www.imperial.ac.uk/earth-science/research/research-groups/pore-scale-modelling/micro-ct-images-and-networks/获取),本研究使用该数据目的在于进一步增强和提高模型重构结果的说服力和可靠度,详细信息请参见表1.可以明显看出,两者不同方向上的渗透率具有明显的不同,也即表明其内部孔隙结构具有典型的各向异性特征.

2.2.1 三维重构计算加速策略

三维各向异性数字岩心模型的主要重构流程与上述二维模型大体相同,但由于重构维度和计算方向数量(由4条变为13条)的双重增长,从而导致三维模型重构所需的计算量急剧增加.此时,若仍然按照原有重构框架进行三维模型的重构,其重构效率将变得极为低下,也即必须要花费更多的时间成本才能顺利完成三维各向异性数字岩心模型的构建.

为解决上述问题,通常需在原有重构框架的基础上进行相应的调整和改进,具体为通过采用各种计算加速策略来提高三维数字岩心模型的重构效率,如图像缩放策略(Song et al.,2017)、系统增量计算更新策略(Song,2019)、交换像素选择策略(Tang et al.,2009;Zhao et al.,2007)和多线程并行计算策略(Ju et al.,2017)等,详细信息请参见对应文献.三维数字岩心模型的实际重构结果表明:上述加速策略的使用能够有效消除计算量增加对重构效率造成的不利影响,重构系统的能量也能够快速且有效地收敛至预先所设定的临界值,计算的时间代价成本处于一种可接受区间.

因此,在后续进行碳酸盐岩和砂岩三维各向异性数字岩心重构的过程中,本研究选择采用了图像缩放、交换像素选择和系统增量计算更新三者相结合的加速策略,以实现高效快速重构三维数字岩心模型的目标.

2.2.2 传统模拟退火算法三维岩心重构

得益于上述重构计算加速策略的加持,以及传统模拟退火算法自身相对较少的计算量,碳酸盐岩和砂岩重构系统的能量值可以快速地降至一个极低的水准,故本研究对其临界值分别设定为2.4×10-10和1×10-7.传统模拟退火算法所获得的碳酸盐岩和砂岩三维数字岩心模型的最终重构结果如图11(三维孔隙结构渲染)和图12(不同方向切面二维孔隙结构切片)所示.

图11 传统模拟退火算法三维数字岩心重构(从左至右分别为:孤立孔隙、连通孔隙和总孔隙)(a) 碳酸盐岩; (b) 砂岩.Fig.11 Reconstruction of 3D digital core based on conventional simulated annealing algorithm (From left to right: Isolated pore, Connected pore and Total pore)(a) Carbonate rock; (b) Sandstone.

图12 传统模拟退火算法三维数字岩心孔隙结构切片图像(从左至右分别为:X-Y、X-Z、Y-Z和X-Y-Z正交切片)(a) 碳酸盐岩; (b) 砂岩.Fig.12 Pore structure slice image of 3D digital core based on conventional simulated annealing algorithm (From left to right: X-Y, X-Z, Y-Z and X-Y-Z orthogonal slice)(a) Carbonate rock; (b) Sandstone.

2.2.3 改进模拟退火算法三维岩心重构

同样的,考虑到改进模拟退火算法重构所使用的统计函数数量较多(共计为26个),导致其计算量相对较大.在确保重构精度的前提下,本研究对碳酸盐岩和砂岩的能量临界值进行适当提高,具体分别为7×10-8和1.5×10-7.

改进模拟退火算法所获得的碳酸盐岩和砂岩三维数字岩心模型的最终重构结果如图13(三维孔隙结构渲染)和图14(不同方向切面二维孔隙结构切片)所示.

图14 改进模拟退火算法三维数字岩心孔隙结构切片图像(从左至右分别为:X-Y、X-Z、Y-Z和X-Y-Z正交切片)(a) 碳酸盐岩; (b) 砂岩.Fig.14 Pore structure slice image of 3D digital core based on improved simulated annealing algorithm(From left to right: X-Y, X-Z, Y-Z and X-Y-Z orthogonal slice)(a) Carbonate rock; (b) Sandstone.

需要特别说明的是,不同于对二维重构岩心各向异性特征的直观评价衡量,三维重构岩心由于孔隙结构的复杂性,无法仅通过对其切片孔隙结构形态的定性直观观察,进而来对其各向异性特征进行准确的评估.因此,必须要采用客观定量的方法来对其进行表征计算,在下一节中将对此进行详细说明.

3 结果对比分析及讨论

3.1 二维孔隙结构形态及统计特征

对比图6a和图9a所示的二维重构岩心图像可以发现,由于传统和改进模拟退火算法重构所选择孔隙结构表征统计函数数量和计算方向的不同,两者所重构岩心孔隙结构的形态及分布也存在较大的差异:直观来看,前者的孔隙结构形态整体较为圆润,狭长形孔隙数量较少,且孔隙的空间朝向分布也较为均匀,并不存在较明显的倾向性;而后者则呈现出了另外一种截然不同的特征,其内部孔隙结构多为狭长形,且相应朝向分布呈现出了十分明显的倾向性,多沿水平方向和负对角线方向分布,这可由其相应方向上的统计函数曲线得到验证.此外,以图3所示的参考图像内部孔隙结构形态作为基准,通过与参考图像进行直观对比,可以看出采用改进模拟退火算法所构建的数字岩心较传统模拟退火算法能够对孔隙结构进行更加准确的复现和表征.

在上述直观定性对比的基础之上,为能够以一种定量的方式进一步对传统模拟退火算法所构建的如图6a所示的数字岩心进行评估,此时不再对其不同方向上的孔隙结构特征进行汇总整合,而是采用改进模拟退火算法孔隙结构特征的提取模式,分别对其不同方向上的孔隙结构特征进行独立计算表征.不同方向上两点概率函数和线性路径函数的具体结果分别如图15a和图15b所示,可以看出,尽管传统模拟退火算法在岩心重构过程中并未考虑孔隙结构不同方向上的差异性,但其最终重构岩心依然表现出较轻微程度的各向异性特征,并非是完全意义上的各向同性分布,具体表现为:两类孔隙结构特征统计函数在水平(竖直)方向与正(负)对角线方向上的差异性较大,而处于两正交方向上(如水平-竖直方向或正-负对角线方向)的统计函数值则差异较小,这也在一定程度上对图6a所示的最终孔隙结构形态进行了解释说明.

图15 传统模拟退火算法重构岩心不同方向统计函数曲线(a) 两点概率函数; (b) 线性路径函数.Fig.15 Statistical function curves in different directions of digital core based on conventional simulated annealing algorithm(a) Two-point probability function; (b) Lineal-path function.

然而,通过进一步与图7所示参考图像的统计函数曲线进行对比,可以十分清晰的发现二者主要存在如下两点不同:① 参考图像内部4条不同方向上的孔隙结构特征彼此之间均具有较显著的差异,且其差异程度远大于图15所示传统模拟退火算法所构建的数字岩心,也即表明两者之间的孔隙结构分布特征存在根本性的不同;② 传统模拟退火算法所构建岩心不同方向上的孔隙结构特征未能够与参考岩心完全准确的对应,以能够表征孔隙长程连通性的线性路径函数为例,可以看出参考岩心不同方向上的孔隙连通性由好到差依次分别为:水平>负对角线>竖直>正对角线,而重构岩心则近似为:水平≈竖直>正对角线>负对角线,也即表明重构岩心未能够将参考岩心不同方向上的孔隙结构特征进行准确复原.上述统计函数曲线的差异也再次定量验证了先前直接对比两者图像所得到的直观视觉感受.

3.2 三维孔隙结构连通性

孔隙空间结构的连通性是评价重构岩心精度的一个关键指标,同时也是开展后续相关研究的一个重要必备前提条件,例如岩心的渗透特性研究,其值等于岩心内部连通孔隙与总孔隙的体积之比.在本论文中,我们主要通过采用式(3)所示的线性路径统计函数,来确保三维重构岩心具有与参考岩心相同的孔隙连通性.

为验证改进模拟退火重构算法的准确性,根据上述孔隙空间结构连通性的定义,基于图11和图13所示由不同算法重构得到的碳酸盐岩和砂岩三维数字岩心,分别对其孔隙空间结构的连通性进行定量化的计算提取,具体结果见表2.

表2 重构三维数字岩心的孔隙结构连通性Table 2 Pore structure connectivity of reconstructed 3D digital core

由表2对比可以看出,两种算法所重构岩心均能够对参考岩心的孔隙连通性进行极大程度的复原和重建:尤其对于碳酸盐岩而言,其重构岩心的孔隙结构连通性优于参考岩心;而由于砂岩参考岩心本身具有极高的孔隙连通性,导致重构岩心的孔隙结构连通性略低,事实上这是由统计随机重构算法的本质属性决定的,其只能使重构岩心的连通性尽可能逼近参考岩心.此外,进一步对比传统算法和改进算法重构岩心的孔隙连通性,可以发现,后者的孔隙结构连通性略小于前者,这可能是由于后者需要考虑更多不同方向上的孔隙像素分布所引起的:对于碳酸盐岩岩心而言,改进模拟退火算法重构岩心的孔隙连通性与其参考岩心更为接近,也即表明改进算法重构岩心具有相对更高的准确性.

3.3 三维孔隙结构渗透特性

众所周知,岩石等多孔介质材料的渗透特性与其内部孔隙结构紧密相关,微观孔隙结构特征直接决定了其宏观的渗透特性.为了进一步检验改进模拟退火算法重构岩心模型的精确度和可靠度,现对碳酸盐岩和砂岩不同方向上的渗透特性进行数值求解预测.

3.3.1 流体流动控制方程

在进行数值模拟计算时,流体通常采用连续介质模型假设,也即假定流体微元是连续紧密排列的,且彼此间不存在任何的间隙.孔隙中流体的流动传输需满足两类物理控制方程:连续性方程(质量守恒方程)和运动方程(动量守恒方程).

在本研究中,我们选用不可压缩的牛顿流体-水作为传输介质,使其在多孔介质内部孔隙结构中进行低速稳态流动.此时,连续性控制方程和运动控制方程的具体形式分别如下所示:

(8)

(9)

式中,vx、vy和vz分别为流体在X、Y和Z方向上的速度分量;p为流体所受压强;μ为流体动力黏度,μ水=0.001 Pa·s.

3.3.2 绝对渗透率预测

在已构建的碳酸盐岩和砂岩三维数字岩心模型的基础上,分别沿着X、Y和Z方向在岩心模型的两端面施加一恒定的压力差值(本研究设定入口端1 Pa,出口端0 Pa),以使流体能够在孔隙结构中进行稳定流动传输.此外,孔隙空间区域内流体的初始压力和初始速度全部设置为零;基质与孔隙边界处流体的速度设置为零,也即采用无滑动边界条件;岩心模型与流体传输方向相平行的四个端面设定为速度对称边界条件,也即流体在与该端面垂直方向上的速度分量为零.基于上述初始和边界条件,通过对式(8)和式(9)进行有限元数值求解,分别完成其内部孔隙结构中不同方向上流体流动传输过程的仿真,进而求解得到整个孔隙空间中流体的压力场和速度场分布.

基于上述求解结果,任取一与流体流动方向相垂直的横截面,通过对该截面上的流体流速进行平面积分,计算得到流体通过该横截面的总流量,最后由达西线性渗流定律计算得到岩心各个不同方向上的绝对渗透率,具体结果见表3.

表3 重构三维数字岩心不同方向上的绝对渗透率Table 3 Absolute permeability in different directions of reconstructed 3D digital core

对比分析表3的计算结果,可以清晰且直观的发现传统算法和改进算法重构岩心不同方向上绝对渗透率的差异:前者不同方向上绝对渗透率数值差异程度较小,可近似认为相同,这与先前赵秀才和姚军(2007)的研究结果相一致,但其与参考岩心的绝对渗透率数值差异较大,表明其未能对参考岩心不同方向上的孔隙结构特征准确表征;而后者不同方向上的绝对渗透率数值则呈现出了较大程度的差异,且同一方向上的绝对渗透率也与参考岩心表现了较高程度的一致性,进而表明改进模拟退火算法能够准确复现参考岩心的孔隙结构特征,也再次验证了其所重构岩心的精确度和可靠性.

3.4 分析及讨论

综合上述对比结果,可以看出本研究所提出的改进模拟退火算法的数字岩心模型重构精度远高于传统模拟退火算法,其准确表征和复现了不同方向上的孔隙结构特征,能够应用于各向异性多孔介质材料数字岩心模型的构建.两种不同算法的重构流程近似相同,差异主要在于孔隙结构特征统计函数的选择和使用.事实上,孔隙结构特征统计函数作为岩心重构过程中的约束限制条件,其对于孔隙结构的全面准确表征起着至关重要的作用,直接影响着模型的最终重构精度.本研究仅从统计函数数量角度对原有算法进行了改进,以使其能够捕获更多不同方向上的孔隙结构特征,后续将考虑采用更多不同类型及数量的高阶孔隙结构表征统计函数,以完成具有更加复杂孔隙空间结构的高精度数字岩心模型构建.

此外,对于岩石内部多组分均呈现各向异性特征的岩心重构(如页岩岩心重构,其内部包含孔隙、层理、有机质及黄铁矿等众多各向异性组分),本研究所提出的重构算法同样适用,具体可通过如下三步操作:首先,对页岩等复杂岩石的原始扫描图像进行多相分割,从将需要进行重构的各组分准确的分割出来(Song et al.,2020);然后,将上述分割后的各组分单独提取出来,并分别作为主要待重构对象,直接使用本文所提出的方法完成其相应各向异性岩心的重构;最后,将上述重构所得到的各组分的各向异性岩心进行组合和叠加(Ji et al.,2019;Song et al.,2021),从而最终完成页岩等较复杂且多组分均呈现各向异性岩心的重构.

4 结论

(1) 基于采用传统模拟退火算法进行数字岩心重构的理论,通过独立引入表征不同方向上孔隙结构特征的统计函数,实现了对真实孔隙结构特征的全方位准确提取表征,进而提出了一种能够适用于各向异性数字岩心模型重构的改进模拟退火算法.

(2) 天然多孔介质材料内部的孔隙结构分布均存在一定程度的各向异性特征,传统模拟退火算法所构建的数字岩心模型理论上应为完全各向同性,实际处理中发现其仍存在较轻微的各向异性特征.然而,上述差异对于其所参考的具有显著各向异性特征的孔隙结构而言,可以近乎忽略不记,也即差异性不够明显和突出,且其不同方向上的孔隙结构特征也未能够与参考岩心完全准确的对应,综合表明传统模拟退火算法重构岩心的精度较差.

(3) 以人工合成二维切片图像及真实碳酸盐岩和砂岩三维CT扫描图像作为参考样本,通过将传统模拟退火算法和改进模拟退火算法所构建的数字岩心模型分别与参考图像中的孔隙结构进行综合的定性和定量对比分析,结果表明:无论是二维孔隙结构的形态和朝向分布,还是三维孔隙结构的连通性和不同方向上的渗透特性,改进模拟退火算法重构岩心均与参考图像表现出了较高的一致性,从而证实了本研究所提出的各向异性多孔介质数字岩心模型重构方法的准确性、有效性和可靠性.

致谢作者感谢帝国理工学院孔隙尺度模拟及成像研究团队分享的碳酸盐岩和砂岩三维CT扫描图像,以及本文所引用文献作者的前期研究工作.此外,作者还要特别感谢责任编辑老师和三位匿名评审专家对本论文提出的建设性修改意见和建议.