随着勘探开发技术的进步及油气需求的增加，非常规油气(如致密气、致密油、煤层气、天然气水合物等)开发受到越来越多人的关注。致密砂岩气作为非常规油气勘探开发的重点，具有极大的资源潜力，是未来相当长时间内天然气增储上产的重要来源^[1-4]。岩石天然裂缝是非常规油气的首选渗流通道，所以具备高渗透性或可改造条件的致密砂岩储层裂缝发育区是非常规油气藏勘探的首要目标^[5]。由于储层孔裂隙结构复杂，进行室内渗流实验难度较大，而随着CT扫描和计算机图像处理技术的发展进步，基于三维数字岩心进行微观渗流模拟成为近年来的研究热点之一^[6-7]。

通过微观渗流模拟能得到渗透率、速度、压力等反映流体在孔裂隙结构中运动及分布状态的重要信息，可充分认识渗流机理的微观本质。目前微观渗流模拟形式包含单相流和两相流两种，前者主要包括流体绝对渗透率模拟^[8-11]及岩石应力敏感性研究^[12-17]，后者主要探究流体的驱替过程及驱替效率的不同影响因素^[18-22]。

在构建三维数字岩心之前，需要对二维CT图像进行分割以提取渗流空间结构。传统分割方法耗时、费力且效率低下，基于深度学习和机器学习自动分割岩心孔裂隙结构成为一种高效的主流方法。其中，U-Net^[23]算法因其独特的结构设计，使得模型能够同时学习到图像的全局信息和局部细节，对于岩心这种具有复杂结构和纹理的对象，能够更好地捕捉到不同尺度的特征；而且能够高效地利用有限的数据进行训练，通过数据增强等技术可以进一步提高模型的泛化能力，从而实现图像的准确分割。相比于一些复杂的深度学习模型，U-Net的训练过程更加稳定，不易出现过拟合等问题，因而成为岩心图像分割任务的理想选择^[24-25]。基于分割完成的三维数字岩心孔裂隙结构模拟计算主要有格子玻尔兹曼方法(Lattice Boltzmann method，LBM)和有限元方法^[26]。LBM方法作为一种高级数值模拟技术，在尝试精准模拟复杂三维岩心中流体渗流行为时，对计算机的计算能力和资源提出了较高的要求。面对真实岩心中错综复杂的孔隙结构和非均质性，该方法在应用中遭遇了多重技术瓶颈和挑战，如计算量大、求解速度慢等^[27-28]。有限元方法中，利用水平集耦合N-S方程能够高精度地跟踪两相流体的界面动态变化。这种方法在解析两相流体的分布特性及驱替流程时，相较于LBM模型，显示出更为显著的优势和有效性^[29]。目前基于数字岩心开展的微观渗流模拟研究大多集中在二维结构模型层面，对于复杂真实三维孔裂隙结构中的流体渗流模拟研究仍有待深入。

本文基于深度学习U-Net算法构建三维数字岩心裂隙结构，首先进行单相流渗透率模拟；其次利用水平集方法耦合N-S方程进行气、水两相流驱替过程模拟，并分析裂隙结构连通性对于水驱气过程的影响。

1.   算法理论

1.1   深度学习U-Net分割算法

U-Net是一种基于FCN的改进模型，其设计灵感来源于字母U，主要用于图像的语义分割任务，其网络结构如图 1所示。该模型由2个主要部分组成，即收缩路径与扩张路径，它们在图像处理过程中相互对称，共同构成U形结构。采用这种方式，即使输入较少的训练图像也能获得良好的分割效果。其关键创新在于将跳跃连接引入解码器中，这种设计帮助解码器更好地利用不同层次的特征信息，使得U-Net在处理图像时能够更好地保留细节信息，提高分割的准确性。

图  1  U-Net网络结构

Figure  1.  U-Net network structure

下载: 全尺寸图片 幻灯片

岩心CT扫描图像展现出多梯度层次、清晰的灰度分级以及丰富的特征细节。针对这一特点，引入U-Net模型来实现数字岩心图像的自动准确分割。该模型的优势在于其训练过程无需大量样本数据，非常适合于岩心CT图像的处理与分析。

1.2   渗透率求解

达西定律是描述流体通过多孔介质流动的基本定律之一，通过在模型的出入口端面对流速积分，可得到流体流经整个模型的体积流量，代入公式(1)即可求得绝对渗透率。它是衡量饱和单相流体通过岩心有效孔隙空间的能力。

式中：Q为流体的体积流量，单位m³/s；S为渗流横截面积，单位m²；K为绝对渗透率，单位m²；Δp为压差，单位Pa；μ为流体动力黏度，单位Pa·s；L为模型的长度，单位m。

1.3   水平集方法与N-S方程

水平集(Level Set)是经典的欧拉界面追踪方法，主要用于描述运动界面的动态变化^[30]。在两相流模拟过程中，引入水平集函数φ来表示两相流体的体积分数，其中0≤φ≤1，φ=0表示流体为气相，φ=1表示流体为水相。

式中：$\frac{\partial \varphi}{\partial t}$为时间累积项；u·$\nabla$φ为对流项；u为流速，单位m/s；γ为重新初始化参数，单位m/s；ε_ls为界面厚度控制参数，单位m。

对于高度复杂孔隙结构模型，选择计算效率较高的N-S方程描述孔隙空间内流体流动。

式中：ρ为流体密度，单位kg/cm³；t为时间，单位s；p为入口压力，单位Pa；I为单位矩阵；K_n为黏性应力张量，n=1时，流体不可压缩，n=2时，流体可压缩；F为表面张力，单位N/m。

2.   数字岩心RVE模型建立

2.1   数据库基础

研究使用的致密砂岩岩心图像样本来自公开数字岩心库(Digital Rock Portal)，岩心为具有天然裂缝(含有胶结物)的Travis Peak砂岩，二维图像尺寸为988像素×988像素，共1 773张，CT扫描分辨率为5.97 μm。切片图像物质组成如图 2所示。

图  2  切片图像物质组成

Figure  2.  Material composition in slice image

下载: 全尺寸图片 幻灯片

2.2   模型建立过程

利用Avizo制作样本数据集，由于相邻像素级切片特征变化不明显，因此隔10个像素点取一张图像，然后经过裁剪、滤波操作，最终模型的输入图像尺寸为128像素×128像素，数量为531张，按照9∶1的比例划分训练集与验证集。训练时采用迁移学习技术，将训练好的VGG16模型移至U-Net编码器部分，使得整个主干特征提取网络具有更强的裂隙结构特征提取能力。图 3为最终训练结果，其中图 3a为精度曲线，说明分割结果的准确性，是预测分割区域与真实分割区域之间的交集与并集的平均比率；图 3b为损失曲线，表示训练过程中模型预测值与真实值间的差异。

图  3  U-Net训练结果

Figure  3.  U-Net training results

下载: 全尺寸图片 幻灯片

由图 3可以看出，模型训练约40轮开始逐渐收敛。为加快训练速度，每训练5轮计算一次精度，最终精度收敛于85%，说明训练效果较好。

代表性体积单元(Representative Volume Element, RVE)是有效表征岩心整体宏观物理性质的最小单元体。参考文献[31]确定砂岩RVE尺寸为400体素×400体素×175体素，整体重构模型如图 4所示。基于上述训练好的深度学习模型，自动分割175张连续切片，以第87张切片为例，分割提取过程如图 5所示。

图  4  砂岩RVE重构模型

Figure  4.  RVE reconstruction model of sandstone

下载: 全尺寸图片 幻灯片

图  5  裂隙分割提取过程

Figure  5.  Fracture segmentation extraction process

下载: 全尺寸图片 幻灯片

图 6a为分割后完整的三维裂隙结构；图 6b则定位出最大的连通裂隙结构。岩石损伤后形成的孔裂隙可以分为孤立孔隙和连通裂隙，其中连通裂隙可视为有效孔隙，对于岩石渗透性有重要作用；孤立孔隙仅对岩石储集性有贡献。因此进行连通性测试，删除孤立、不连通裂隙，提取最大连通裂隙结构(图 7a)。连通性测试采用Avizo中的骨架化模块进行，骨架模型如图 7b所示，其将裂隙区域内的三维结构以骨架的形式表现出来，能够更容易识别和理解材料内部的裂隙分布和连通性。骨架化程度越高，其通道尺寸越大；骨架化程度越低，其通道尺寸越小。计算总裂隙率为6.14%，连通裂隙率为5.72%，连通裂隙占总裂隙的93.16%。

图  6  三维重构结果

Figure  6.  Three-dimensional reconstruction results

下载: 全尺寸图片 幻灯片

图  7  最大连通裂隙结构及其骨架化模型

Figure  7.  Largest connected fracture structure and its skeletonized model

下载: 全尺寸图片 幻灯片

为研究RVE尺度下流体渗流特性，借助AVIZO-COMSOL交互技术，对选取的RVE模型进行流体渗流模拟研究^[32]。几何模型及网格模型如图 8所示。

图  8  COMSOL生成计算模型

Figure  8.  Computational models generated by COMSOL

下载: 全尺寸图片 幻灯片

为确定模型渗透率，模拟分别采用水和甲烷作为流体材料；驱替过程模拟同时采用水、甲烷两种流体材料，流体流动状态均为层流。流体从模型的一端注入，另一端流出，且初始流速为零，如图 8a所示。数字岩心模型中相对立的两面分别为流入及流出边界，其余流动边界及孔壁视为无滑移壁面；流入、流出边界利用压力大小控制。具体模拟参数如表 1所示。

表  1  渗流模拟参数设定

Table  1.  Parameter settings for seepage simulation

流体材料	密度/(g/cm3)	动力黏度/(Pa·s)	压差/Pa
水	1	0.001	1 000
甲烷	6.5×10-4	1.1×10-5	1 000

下载: 导出CSV 

| 显示表格

3.   模拟结果与分析

3.1   渗透率模拟

选择稳态求解器，计算自由度个数为68 164(加上1个内部自由度)，计算得到渗流过程中速度流线和压力分布情况(图 9)。

图  9  模型在不同流体下的速度和压力云图

Figure  9.  Contour maps of model velocity and pressure with different fluids

下载: 全尺寸图片 幻灯片

当流体为水(设为不可压缩流体)时，模型最大流速为1.14 m/s，最大压力为1 070 Pa。文献[31]利用LBM方法计算的平均渗透率为4.95 μm²，而利用有限元方法模拟计算得到的渗透率约为5.1 μm²，二者误差约为3%，验证了模拟结果的准确性。改变流体材料为甲烷(可压缩流体)后，最大压力为1 100 Pa，最大流速为41.9 m/s，模拟计算得到渗透率约为0.14 μm²。由此可见，流体性质的不同影响流动时的压力和速度，进而影响模型整体渗透性。

从图 9可以看出，尽管所有裂隙相互连通，但由于裂隙几何形状的复杂性，不可避免地导致了某些裂隙凝滞(即连通性差)，即流体未流入这些裂隙中，因此速度场并未覆盖整个模型空间，流线曲折交错，在大裂隙连接通道处较为密集，速度也较大；沿渗流方向，裂隙所受压力逐渐减小。

3.2   水、气两相流驱替模拟

致密岩体采气过程是水、气两相流的驱替过程。在确定模型裂隙特征后，进行两相流驱替过程研究(图 10)。

图  10  不同时刻模型速度与压力分布

Figure  10.  Model velocity and pressure distributions at different moments

下载: 全尺寸图片 幻灯片

由图 10可见，驱替过程中流体速度不断增大，最大值为38.7 m/s；压力呈现出先增大后减小的趋势，最大值为1 090 Pa。整体而言，渗流过程中裂隙结构的连通性影响着压力场和速度场的分布特征。连通性较好的裂隙通道先出现高速流线，且越宽、离出口越近的通道流线数量越多，压力分布也更均匀，优先流效应明显^[32]；对于单个裂隙通道，沿出入口通道方向，速度变化不大，且中心位置速度最大、边缘孔壁处速度最小，该现象与方辉煌等^[33]的研究结果一致。连通性越差的裂隙中，“压力堵塞区”范围越大，整体流线越闭塞且数量较少，从而导致渗透率越小。由于应力集中，孔径越小，承受压力越大。

裂隙结构初始状态充满甲烷气体，通过注水加压进行采气。不同时间流体状态分布如图 11所示，红色代表的气相不断减少，蓝色代表的水相逐渐增多，可见气体采收率逐渐增大。图 12为流体饱和度变化曲线，整个过程中，两相饱和度不断变化，通过连通性相对较好的裂隙通道(如图 11中椭圆圈内所示)时，流体饱和度出现明显的改变；通过局部细长、狭窄、连通性较差的裂隙通道(如图 11中圆圈内所示)时，直至渗流结束，其通道内流体饱和度也无明显变化，成为残余气相的主要赋存空间。

图  11  不同时刻流体状态分布

Figure  11.  Fluid state distributions at different moments

下载: 全尺寸图片 幻灯片

图  12  流体饱和度变化曲线

Figure  12.  Fluid saturation variation curves

下载: 全尺寸图片 幻灯片

4.   结论

(1) 基于深度学习U-Net算法，可实现岩心CT扫描图像的自动准确分割，以此为基础建立真实三维数字岩心结构及有限元模型进行微观渗流模拟可行性较高。

(2) 单相流模拟过程中，分别以气、水作为流体材料时，模型的压力差别较小，但速度相差较大，从而导致渗透率有较大变化；连通性差的裂隙，不存在或存在少量低速流线；沿渗流方向，裂隙所受压力逐渐降低。

(3) 基于水平集方法耦合N-S方程进行了气、水两相流模拟，驱替过程中气相不断减少，水相不断增多，水驱气的效率逐渐增大，连通性较好的裂隙成为主渗流通道；渗流结束后仍存在残余气相，其主要分布在模型裂隙连通性较差的通道。

(4) 当流体在压力作用下进入致密储层裂缝时，会优先从宽且笔直的通道填充裂隙空间；流体的流动路径曲折多变，裂隙连通性以及流体的性质都会影响填充的速度和程度，进而影响最终的采收率。