基于生成对抗网络的页岩三维数字岩芯构建

杨永飞，刘夫贵，姚 军，宋华军，王 民

1.深层油气重点实验室·中国石油大学（华东），山东 青岛 266580；2.中国石油大学（华东）石油工程学院，山东 青岛 266580；3.中国石油大学（华东）海洋与空间信息学院，山东 青岛 266580；4.中国石油大学（华东）地球科学与技术学院，山东 青岛266580

引言

随着中国对常规油气资源的不断开发、利用，常规油气藏普遍进入开发中后期，油气勘探与开发逐渐转向非常规油气藏。其中，页岩气分布范围广、资源量大、稳产周期长，成为了油气资源的有力补充[1]。页岩孔隙结构复杂，具有强烈的非均质性，页岩孔隙空间包括纳米级有机粒内孔隙、纳米—微米级无机粒间空隙以及微米—毫米级天然裂缝[2]，常规岩芯实验很难对页岩储层孔隙空间特征进行定量描述。而页岩储层宏观性质由其微观孔隙空间性质及渗流规律决定，因此，准确表征页岩储层孔隙结构、研究页岩纳微孔隙中的流体渗流规律是页岩油气藏高效开发的关键。数字岩芯作为孔隙级微观渗流理论研究的基础平台，可以使岩芯样品可视化、定量化，为精细描述页岩储层孔隙空间、研究页岩储层内流体运移机制提供了有效途径。目前，数字岩芯的构建方法主要分为两大类：物理实验法和数值重建法[3]。通过聚焦离子束扫描电子显微镜（Focused Ion Beam-Scanning Electron Microscope，FIB-SEM）、纳米CT 可以直接获取页岩的三维数字岩芯[4]。这样获得的数字岩芯分辨率高、结构准确，但是这些方法在获取岩芯图像时对实验设备要求高，价格昂贵，且完成一块小尺寸的岩芯的扫描往往需要很长的时间。借助扫描电子显微镜（SEM）等实验方法可以获取页岩岩芯的二维切片数据，然后通过数值重建法重构三维数字岩芯。数值重建法主要有高斯模拟法、模拟退火法、马尔科夫链蒙特卡洛法、多点地质统计学方法、过程模拟法以及混合法[5-8]。数值重建法重构的三维数字岩芯具有良好的孔隙结构性质。重建结果受不同的约束条件控制，约束条件越多，重建结果越好，但计算成本也就越高。近年来，机器学习的发展给图像处理、图像重建研究提供了新的思路。支持向量机是最早用于数字岩芯重建的机器学习方法[9]。迁移学习[10]、卷积神经网络[11]等深度学习方法也被很多学者应用到多孔介质重建和微观孔隙结构性质预测等方面。自2014 年Goodfellow 等[12]提出生成对抗网络（Generative Adversarial Network，GAN）之后，因其具有的强大的图像生成能力，很多学者都开始着眼于应用GAN 模型进行多孔介质微观结构的重建工作研究。Mosser 等[13-14]将深度卷积生成对抗网络（Deep Convolutional Generative Adversarial Network，DCGAN）应用到了均质性较好的多孔介质的三维重建。Shams 等[15]结合了自动编码器（Auto Encoder，AE）和生成对抗网络，重建了同时包含粒间孔隙信息和粒内孔隙信息的多尺度数字岩芯。此外，条件生成对抗网络（Conditional Generative Adversarial Network，CGAN）[16]、循环一致性生成对抗网络（Cycle-Generative Adversarial Networks，Cycle-GAN）[17]、基于样式的生成对抗网络（Style-Generative Adversarial Networks，Style-GAN）[18]等不同形式的生成对抗网络也被用来重建数字岩芯，包括由岩芯二维切片重建三维数字岩芯、由三维岩芯图像重建三维数字岩芯，但构建的对象往往是具有较好均质性的岩芯。Zha等[19-20]利用DCGAN、WGAN（Wasserstein Generative Adversarial Network）进行了页岩二维图像的重建，证明了可通过生成对抗网络重建非均质性强的岩芯。考虑页岩岩芯不易获取，并且孔隙结构复杂、非均质强等问题，本文基于真实页岩岩芯三维聚焦离子束扫描电子显微镜图像，对原始生成对抗网络模型的结构重新设计。为保证重建结果充分反映页岩岩芯的孔隙结构信息，增大了训练样本的尺寸。通过对原始数字岩芯和重建数字岩芯的孔隙度以及孔隙空间的几何结构参数和拓扑结构参数进行对比分析，验证重建岩芯的准确性以及构建方法的可靠性。并通过训练好的生成模型构建了多个数字岩芯，计算了孔隙度、孔隙结构参数均值以及变化区间，验证了生成模型的稳定性。

1 生成模型与评价方法

1.1 生成对抗网络模型

本文对生成对抗网络模型的结构进行重新设计，以此构建三维页岩数字岩芯。生成对抗网络是非监督式的深度生成模型，由两个完全卷积的神经网络组成，即生成器和判别器。生成器的作用是实现从一个服从高斯分布的低维随机分布到高维图像的映射，生成器的输入是随机分布，输出是一定尺寸的图像。判别器的作用是区分图像是来自真实的训练图像集还是来自生成器生成的“假”图像集，判别器的输入是图像，输出是一个0～1 的标量，代表这个图像属于真实图像集的概率。判别器的目标是尽可能准确地区分一个图像是来自真实图像集还是由生成器生成的图像，而生成器的目标是尽量生成判别器无法区分真假的图像[12]，如图1 所示。

图1 生成对抗网络原理图Fig.1 Schematic of GAN

判别器是一个二分类器，输入真实图像x或者生成器生成的“假”图像，输出一个标量表示输入图像是真实图像的置信度，即输入图像来自于满足preal_data(x)分布的真实图像集的概率，该过程用D(x,θd)表示，其中，θd是决定判别器网络结构的权重参数。

判别器输出的标量在0～1，输出为1 表示输入的图像为真实图像，输出为0 表示输入的图像为生成器生成的“假”图像。输出的值越接近于1 表示输入的图像越符合真实图像，判别器以此判断生成器生成的图像是否符合真实图像的数据分布。

生成器是一个深度生成模型，将一个低维的随机分布pz(z) 映射成一个高维的满足数据分布pfake_data(x) 的“假”图像，该过程用G(z,θg) 表示，其中，θg是决定生成器神经网络结构的权重参数。

本文使用的随机分布pz(z)是（0，1）正态分布。生成对抗网络的训练过程是不断地交替优化判别器和生成器。真实图像的标签为1，“假”样本的标签为0，那么可以用二元交叉熵定义生成对抗网络的损失函数为

式中：

D—判别器函数；

G—生成器函数；

x—真实图像，无因次；

preal_data(x)—真实图像集的数据分布，无因次；

pfake_data(x)—“假”图像集的数据分布，无因次；

z—随机噪声数组，无因次；

pz(z)—随机噪声分布，无因次。

对于判别器的优化，目的是使判别器的判别能力不断提高，即判别器的输入为x时，D(x)应该尽可能接近于1；判别器的输入为时，D()应该尽可能接近于0。于是优化过程是使判别器的损失函数最大化，即

对于生成器的优化，目的是使生成器能够生成以假乱真的图像，使判别器无法识别“假”图像，即判别器的输入为时，D() 应该尽可能接近与1。于是优化过程是使生成器的损失函数最小化，即

在实际训练过程中，当生成器的生成能力很差时，1−D(G(z))趋近于1，此时生成器损失函数的梯度接近饱和，不利于梯度的更新，从而用式（4）作为损失函数进行生成器的优化。

1.2 孔隙结构评价参数

为了评价重建三维数字岩芯的准确性，本文选取了能表征岩芯孔隙度的单点概率函数作为重建数字岩芯的初步评价参数。同时，针对三维数字岩芯可提取孔隙网络模型[8]，以孔隙网络模型为研究平台进一步分析岩芯的孔隙结构特征，包括几何结构特征参数和拓扑结构特征参数。

依据上述几类结构评价参数，对原始岩芯和重建岩芯的孔隙结构参数进行对比分析。

1.2.1 孔隙度

基于三维岩芯图像，可利用单点概率函数对岩芯孔隙度进行表征。假设在三维数字岩芯系统中孔隙相所占据的总区域为v，其相对于整个系统的体积分数，那么可定义孔隙相的相函数为

式中：P—孔隙相的相函数；

r—三维数字岩芯中的一点；

v—三维数字岩芯中孔隙相占据的总区域。

依据孔隙相的相函数定义，可以得到三维数字岩芯的孔隙度，即

式中：φ—孔隙度，%。

1.2.2 几何结构参数

孔隙空间的几何特征参数是反映重建岩芯结构准确性的重要评价参数。孔隙空间的几何结构参数包括孔隙半径、喉道半径、喉道长度、孔隙空间总体半径、形状因子及迂曲度等参数。

其中，孔隙半径R为孔隙内切球的半径，喉道定义为两个孔隙相连接通道，喉道长度L的定义为

式中：d—喉道所连接的两个孔隙中心点的实际距离，nm；

R1，R2—喉道所连接的两个孔隙的半径，nm。

形状因子可以表征孔喉截面形状、定量研究孔喉截面形状的复杂性，定义为孔喉截面面积与孔喉截面周长平方的比值，即

式中：g—形状因子，无因次；

A—孔喉截面面积，nm2；

p—孔喉截面周长，nm。

迂曲度可以描述两个连通孔隙间喉道的曲折程度，定义为两个连通孔隙之间的实际距离与两个连通孔隙之间的最短距离之间的比值，可表达为

式中：τ—迂曲度，无因次；

d′—两个连通孔隙之间的最短距离，nm。

1.2.3 拓扑结构参数

孔隙拓扑结构参数可以进一步评价重建岩芯在孔喉连通特征方面的准确性。

采用配位数、网络连通性函数作为评价重建岩芯拓扑结构特征的参数。配位数是指与孔隙连通的喉道的数量，描述孔隙空间的连通程度。

网络连通性函数可以用比欧拉示性数来表征，可以描述孔隙网络模型中孔隙喉道单元的拓扑结构信息[21-22]。

2 三维页岩数字岩芯构建

2.1 数据准备

使用的页岩岩芯图像样本来源于公开的数字岩芯库（Digital Rock Portal），该岩芯图像是由FIB-SEM 扫描获取的三维真实页岩岩芯图像，Andrew[23]基于该数据开展了相关的研究工作。

经Andrew 处理，原始岩芯图像被分割成5 相，即孔隙相和其他4 种矿物组分的骨架相。岩芯图像体素尺寸为1 024×861×606，分辨率为5 nm×5 nm×5 nm，样品实际尺寸为5.12µm×4.30µm×3.03µm。

考虑到页岩油气藏中流体的渗流发生在孔隙空间中，因此，本文的研究着眼于页岩孔隙空间的重建。对数据体进行预处理，将4 种组分的骨架相重新划分为一相，和孔隙相重组为一个新的三维数据体，即包含骨架相和孔隙相两相的三维数字岩芯，如图2 所示。

图2 原始页岩数字岩芯Fig.2 Initial shale digital rock

为了确保训练生成对抗网络模型有足够的训练图像，从预处理后的数字岩芯中提取一定尺寸的子数据体作为训练图像。页岩岩芯本身具有很强的非均质性，在提取子数据体作为训练图像时需要确保单个训练图像可以反映岩芯足够完整的信息。从分割得到的1 024×861×606 三维数据体中每隔64 个体素步长提取一个体素尺寸为128×128×128 的子数据体，共得到1 873 个数据体作为训练样本图像。

2.2 模型构建及训练

2.2.1 模型构建

本文中生成器和判别器均采用卷积神经网络对岩芯图像特征进行表示，该深度卷积生成对抗网络（DCGAN）结构是Radford 和Luke[24]针对原始生成对抗网络训练不稳定提出的改进，在DCGAN 的基础上结合所用的数据进行相关神经网络层结构的重新设计。其中，卷积是将一个卷积核大小的图像像素加权平均，实现特征提取的功能；卷积核是一个由权值组成的数组；卷积核的数量表示用几个卷积核对一张图像进行处理，得到一定数量的特征图，得到的特征图的数量和卷积核的数量一致。图3 给出了本文中的生成器的网络结构，判别器的结构与之类似。

图3 生成器结构示意图Fig.3 General structure diagram of generator

判别器的输入为三维岩芯图像，该图像的体素尺寸为128×128×128，该岩芯图像有两个来源：训练图像集和生成器生成的图像。通过5 层的体积卷积层（Conv3D）进行特征提取与表征，得到一定数量的体素尺寸为4×4×4 的特征图。将得到的所有特征图传递给最后一个只包含一个卷积核的卷积层。将最后一层卷积层的输出经Sigmoid 激活函数作用，获得一个0～1 的标量作为输入图像来自于训练图像集的概率。除了最后一层卷积层之外，其他卷积层均采用批量归一化和LeakyReLU 激活函数。生成器从输入的随机噪声数组中进行上采样，生成三维岩芯图像，岩芯图像的体素尺寸与训练图像一致。上采样过程是通过体积转置卷积进行的。生成器对随机噪声数组进行重组，将随机噪声数组转换成体素尺寸为4×4×4 的特征图。每一层特征图数量如图3 所示，其中，ng为生成器初始特征图数量，数值为64。经过5 层的转置卷积（ConvT3D）、批量归一化以及ReLU 激活函数的作用后，得到体素尺寸为128×128×128 的岩芯图像。生成器和判别器网络结构参数设置见表1。

表1 生成器和判别器的网络结构参数Tab.1 Network structure parameters of generator and discriminator

2.2.2 模型训练过程

生成对抗网络的训练分为两步：第一步，固定生成器，训练判别器。保证生成器的参数不变，调整判别器的参数，使得式最大化，提高判别器的判别能力，使判别器更加准确的区分输入图像的来源。第二步，固定判别器，训练生成器，提高生成器的生成能力，使生成器能够生成更加符合训练图像数据分布preal_data(x)的图像。在此过程中，保证判别器的参数不变，生成器从正态分布的潜在空间中进行上采样，生成与训练图像尺寸相同的图像，并传递给判别器，通过最大化式调整生成器的参数，提高生成器的生成能力。

在判别器和生成器不断交替优化的过程中，判别器区分真实图像和生成图像的能力及生成器“欺骗”判别器的能力均不断提升。

最终收敛时，preal_data(x)=pfake_data(x)，生成器和判别器达到纳什均衡。此时，判别器的返回值为0.5，表明判别器认为它接收到的图像满足preal_data(x)的概率为0.5，无法区分真实图像和生成图像。训练过程如图4 所示。

图4 生成对抗网络训练过程Fig.4 The training procedure of the generative adversarial network

3 重建数字岩芯结果与讨论

3.1 三维页岩数字岩芯重建

针对预处理的1 873 张训练图像，完成了生成对抗网络模型的训练。

为了提高模型训练的稳定性，采用标签平滑的原理，设置真实标签为0.9～1.2 的随机数，“假”标签为0～0.3 的随机数。模型训练好之后，在计算资源允许的条件下，理论上可以生成任意尺寸的数字岩芯。本文利用训练好的生成模型，生成了16 个体素尺寸为400×400×400 的三维数字岩芯，这里的16是训练过程中用到的小批量样本的数量。

图5 给出了部分三维数字岩芯和二维切片的展示，从图中可以看出，生成器生成的数字岩芯很好的描述了原始岩芯的孔隙空间。

图5 生成器生成的三维数字岩芯及其二维切片Fig.5 The generated 3D digital rocks and 2D slices by generator

生成器不仅生成了原始岩芯的不同尺寸的孔隙，同时捕捉到了孔隙的空间分布信息，比如非均质性：没有孔隙相，只有被骨架相占据的大片区域。这初步验证了训练好的生成对抗网络模型具有很好的生成能力，可以生成满足给定孔隙结构信息的页岩数字岩芯。

为了进一步验证生成器对页岩数字岩芯的重建效果，对重建数字岩芯进行定量评价，并与原始岩芯的结构参数进行对比分析。

3.2 孔隙结构参数对比分析

3.2.1 孔隙度

利用生成的16 个体素尺寸为400×400×400的三维页岩数字岩芯，基于数字岩芯研究平台，借助单点统计函数定义，对数字岩芯系统孔隙相进行了统计平均，计算出三维数字岩芯的孔隙度，并对原始数字岩芯和重建数字岩芯的两个孔隙度进行简单的对比分析。

两个重建数字岩芯分别记为1#、2#。表2 给出了这3 块数字岩芯的基本孔隙结构参数。

表2 不同数字岩芯基本孔隙结构参数Tab.2 Basic pore structure parameters of different digital rock

经统计计算，原始页岩数字岩芯的孔隙度为0.054，重建岩芯1#、2# 的孔隙度为0.053、0.055。可以看出重建数字岩芯的平均孔隙度与原始岩芯的孔隙度只相差0.001，这说明训练好的模型可以生成满足孔隙度要求的数字岩芯。

同时，这16 个重建数字岩芯的最大孔隙度为0.055，最小孔隙度为0.051，重建数字岩芯的孔隙度变化不大，表明模型重建岩芯的稳定性较好。

为进一步对重建数字岩芯孔隙结构特征进行评价，从重建的数字岩芯中任选两个岩芯进行孔隙网络模型的提取，进而计算其几何结构参数和拓扑结构参数，进一步从几何结构参数和拓扑结构参数角度进行对比分析，研究重建数字岩芯的孔隙结构的准确性。

3.2.2 几何结构参数对比分析

重建岩芯1#、重建岩芯2#及原始岩芯的孔隙空间几何结构参数对比见图6。

图6 不同数字岩芯孔隙空间几何结构参数对比Fig.6 Comparison of geometrical pore structure parameters of different digital rocks

可以看出，重建岩芯1#、重建岩芯2#和原始岩芯的孔隙半径、喉道半径、喉道长度、总体半径、形状因子和迂曲度等孔隙几何结构参数的概率分布具有很高的一致性。

此外，从表2 给出的不同数字岩芯基本孔隙结构参数可以得到，原始岩芯的平均孔隙半径为18.77 nm，而重建岩芯1#、2#的平均孔隙半径分别为18.22 nm、18.07 nm，与原始岩芯平均孔隙半径差异不足1.00 nm。同样，重建岩芯1#、2#的平均喉道半径与原始岩芯仅相差不足0.5 nm。这表明，通过训练好的生成器生成的数字岩芯可以很好地再现原始岩芯的孔隙结构特征，包括孔隙半径、喉道半径、喉道长度及孔隙空间的整体尺寸分布。重建岩芯1#、2#的平均形状因子分别为0.042 61、0.042 48，原始岩芯的平均形状因子为0.042 68。结合形状因子和迂曲度对比结果，重建岩芯可以很好地表征岩芯孔隙喉道截面的形状，描述岩芯喉道的弯曲程度。

3.2.3 拓扑结构参数对比分析

不同数字岩芯的孔隙拓扑结构参数对比结果如图7 所示。图7a 给出了3 块数字岩芯的配位数概率分布对比，原始岩芯的孔隙配位数以2 和3 为主，平均配位数为2.54。重建岩芯#1、#2 的平均配位数为2.61、2.30，3 块岩芯的配位数概率分布趋势一致，峰值呈现的区间一致。通过生成对抗网络模型生成的数字岩芯再现了原始岩芯中孔喉配比关系，很好的表征了岩芯的连通程度。

图7b 是不同数字岩芯的比欧拉示性数对比情况，原始岩芯的比欧拉示性数为0 时对应的临界最小半径为8.61 nm，重建岩芯#1、#2 的临界最小半径分别为8.11 nm、7.39 nm，与原始岩芯相差不足1.00 nm，表明重建岩芯可以很好地描述原始页岩岩芯的孔隙连通性特征。

图7 不同数字岩芯孔隙空间拓扑结构参数对比Fig.7 Comparison of topological pore structure parameters of different digital rocks

通过对重建岩芯与原始岩芯的孔隙度、孔隙几何结构参数、拓扑结构参数的对比分析，可以发现，基于生成对抗网络模型可以学习到原始页岩岩芯的孔隙空间特征，利用训练好的生成网络可以对页岩岩芯进行重建。重建的数字岩芯不仅可以保留原始岩芯的孔隙度、孔隙和喉道半径分布、喉道长度、喉道弯曲程度等几何结构特征，还可以表征岩芯孔隙喉道配比关系、孔隙喉道连通程度等拓扑结构特征，充分证明了重建的数字岩芯具有很高的准确性。

4 生成模型稳定性分析

为了验证训练好的生成模型可以稳定地生成数字岩芯，即训练好的模型可以任意生成孔隙结构特征与实际相符的数字岩芯。基于16 个数字岩芯提取孔隙网络模型，并分析了其孔隙结构特征参数。

表3 给出了原始岩芯孔隙结构参数和生成器生成的16 个数字岩芯的平均孔隙结构参数。可以看出，这一组重建岩芯的平均孔隙半径的平均值为19.44 nm，与原始岩芯的平均孔隙半径仅相差0.68 nm。重建岩芯的平均喉道半径、平均喉道长度、平均形状因子以及平均配位数等孔隙结构参数的平均值与原始岩芯的差异也很小，表明重建的数字岩芯符合期望。

表3 不同数字岩芯孔隙结构参数Tab.3 Pore structure parameters of different digital rock

此外，这一组重建岩芯的平均孔隙半径最小为18.07 nm，最大值为20.81 nm，变化幅度不超过7.05%。重建岩芯的平均喉道半径变化幅度不超过14.18%，平均喉道长度的变化幅度不超过6.61%，平均形状因为的变化幅度不超过1.97%，平均配位数的变化幅度不超过6.70%。

通过上述孔隙结构特征参数的分析可以发现，重建数字岩芯的孔隙结构特征稳定在一个范围。

通过这组重建页岩数字岩芯的孔隙结构参数的变化范围可以发现，训练好的生成器可以稳定地生成与训练所用数字岩芯结构参数特征一致的数字岩芯。重建的数字岩芯具有原始岩芯孔隙特征参数分布特征，但具体的孔隙结构又不完全相同，证明了重建数字岩芯具有多样性与可靠性，可以满足不同性质的研究需要。

5 结论

（1）基于真实页岩岩芯的三维FIB-SEM 图像获取大尺寸训练图像，分别用6 层的深度卷积神经网络作为生成对抗网络的生成模型和判别模型，完成了模型的训练及页岩数字岩芯的构建。

（2）对重建岩芯的孔隙度、几何结构参数和拓扑结构参数进行了对比分析，证明了重建的页岩数字岩芯可以再现页岩岩芯的孔隙结构性质，描述页岩岩芯的孔隙连通关系以及孔隙喉道配位关系。

（3）在计算资源允许的条件下，训练好的生成模型可以快速生成任意数量、任意尺寸的数字岩芯。利用训练好的模型生成了多个数字岩芯，统计计算了其结构参数的平均值，证明了生成的数字岩芯具有稳定的孔隙空间特征，该生成模型具有良好的稳定性。