基于马尔可夫链-蒙特卡洛法的碳酸盐岩三维数字岩心模型重构

聂昕，胡俊强,万宇，张超谟，张占松

1.油气资源与勘探技术教育部重点实验室(长江大学)，湖北 武汉 430100 2.非常规油气省部共建协同创新中心(长江大学)，湖北 武汉 430100 3.中国石油集团测井有限公司，陕西 西安 710077

数字岩心不仅能够反映储层岩石复杂的孔隙空间，还可以展示岩石的固体骨架特征[1]。岩石各类物理性质(如弹性、导电性等)的数值模拟是以三维数字岩心为基础进行的，因此三维数字岩心模型的准确性直接影响最终岩石物理实验模拟的结果[2]。目前常用的三维数字岩心建模方法主要分为过程法和随机法两大类[3,4]。过程法是模拟颗粒沉积、压实等过程进行建模，适合于普通砂岩储层；随机法是基于已经获得的岩石二维图像，以其信息统计特性为约束，利用数学算法进行三维数字岩心的构建，使其与原始的岩石二维图像的统计特性较为接近，适用范围更广泛。常用的随机法有完全随机法[4]、高斯场法[5-8]、顺序指示模拟法[9-11]、模拟退火算法[12-15]、多点统计法[16，17]以及马尔可夫链-蒙特卡洛法(MCMC法)[18-21]等。MCMC法自2004年被Wu等[18]引进到数字岩心重构以来，因其具有计算速度快、适用范围广及可体现各向异性等优势被广泛用于重构三维数字岩心，在碳酸盐岩、页岩等储层的数字岩心建模方面有所应用。王晨晨等[21]基于MCMC法使用叠加法重构出能同时描述大孔隙和微孔隙性质的碳酸盐岩双孔隙三维数字岩心，并对建立的数字岩心模型的孔喉特征进行了分析，结果表明MCMC法能够基于二维真岩心薄片快速重构出岩石的三维数字岩心，且利用叠加法构建的碳酸盐岩双孔隙数字岩心具有较高的孔隙连通体积比，孔隙之间的连通性相对于单一尺度的建模结果显著提高。聂昕等[22]选用MCMC法并加以改进，分别构建页岩气储层各微观组分(如黏土矿物、孔隙、有机质等)的三维模型，将上述模型叠合嵌套，得到多矿物组分的页岩储层三维数字岩心，对岩心模型分析的结果表明，利用MCMC法和嵌套法重构的三维数字岩心可用于页岩气储层岩石物理性质的数值模拟研究，同时也提出了重构方法中存在的问题和不足，对其他岩石储层数字岩心重构方法的改进提出了建议。目前，针对MCMC法建模的研究仅限于进行岩心的重建，然后选择合适的尺寸进行重构，再研究岩石储层的相关性质。然而，在常用的基于MCMC法的三维数字岩心建模方法中，由于确定马尔科夫键转移概率的邻域模板尺寸有限，仅有几个像素，因此建模结果的尺寸(即边长像素数)的选取直接影响到了建模结果的准确性，不同尺寸的重构结果差异巨大，而前人未对该现象进行过研究和总结。为此，笔者以碳酸盐岩为例，基于二维切片信息，利用MCMC法重构其不同像素数尺寸的三维数字岩心，并对其结果进行孔隙结构特性分析，从而讨论MCMC法重构三维数字岩心模型的准确性。

1 MCMC法原理

MCMC法来自于图像处理过程中马尔可夫随机场(MRF)的使用。马尔可夫链是指一种状态序列，该序列每个位置的状态均取决于该位置前面有限个位置的状态，而存在这个状态的概率称为转移概率[23，24]。利用MCMC法重构图像时，先利用马尔可夫链获得原始图像的转移概率，然后利用转移概率进行新图像赋值重构。

1.1 二维MCMC法原理

二维图像的本质是一个矩阵，图像的高和宽为矩阵的行和列。一般情况下，采用二值化的图像来表示仅包括孔隙和骨架信息的岩心，即岩石孔隙和岩石骨架分别用1和0来表示，因此该图像矩阵的每个位置只有1和0两种可能性。图像的结构特征可以用概率分布函数进行表征。设一个二维图像的像素数为n，x=(x1，x2，…，xn)表示该像素处的值，在二值化的岩心图像中，xi只能为0(骨架)或1(孔隙)。基于原始图像，利用MCMC法重构一个具有相近性质的二维多孔介质的图像，理论上需要获取一个完整的x的概率分布函数，即p(x)。传统马尔可夫链需要从样本空间中进行完美随机抽样来获得概率分布函数，但是由于数字岩心的样本空间非常大，导致该完美随机抽样方案几乎不可能实现。因此，WU等[18]引入邻域模板的概念，并利用其统计结果来代替求取概率分布函，其本质是将模型假设成图像中任意一点的状态由该点附近区域极少数点的状态决定，即对于一个点s，用Λ-s表示该图像中s点以外的其他所有点，则应存在一个满足式(1)的s点的邻近区域集合Ns：

p(xs|x(Λ-s))≈p(xs|x(Ns))

(1)

通常情况下，想要准确获得式(1)中的Ns需要大量的计算，不利于实际应用。WU等[18]在利用MCMC法重构二维岩心图像时，定义了影响某1点状态的邻域为4个点，即该点上面3个点(正上方、左上方和右上方)以及其左边相邻的1个点，而该点右边相邻的点同时受到包括这个点和其邻域4个点的影响，这样就构成了一个5-6点邻域模板系统(见图1(a))，即：

(2)

式中：(i，j)为该图像中位于第i行第j列的体素。

图1 5-6点邻域和2点邻域模板系统示意图Fig.1 Schematic diagram of formwork system of 5-6 point and 2 point neighborhood

在图像最上面边界的点，由于其上方没有点，则简化为该点只受其左边点的影响，即2点邻域模板系统，如图1(b)所示。在重构之前，先进行原始图像的遍历，获得各个邻域模板系统状态情况下的转移概率，再利用该转移概率，采用蒙特卡洛方法对新的图像进行构建。

1.2 三维MCMC法原理

以二维MCMC法原理为基础，将其拓展到三维的情况。设三维数字岩心图像为L行、M列和N层的三维样本空间VLMN，点(i，j，k)为该图像中位于第i行j列k层的体素，x(Vijk)为点(i，j，k)处的状态。三维与二维一样，均采用对原始图像进行遍历扫描的算法来计算转移概率。在重构赋值时，利用转移概率，采用蒙特卡洛方法逐行、逐列、逐层地对体素点进行赋值。第1层采用二维MCMC法，第2层及以上采用三维MCMC法。三维马尔可夫链的模型为：

p(xijk|{xlmn:0

=p(xijk|xi-1,j,k,xi,j-1,k,xi,j,k-1)

(3)

式(3)表示X、Y、Z共3个方向上条件概率的共同作用。

对于任意属于VLMN的点(i，j，k)，联合概率函数为：

(4)

根据式(4)，三维马尔可夫链的转移概率为：

p(xijk|{xlmn:(l,m,n)≠(i,j,k)})

=p(xijk|{xlmn:(l,m,n)∈N(ijk)})

(5)

式中:N(ijk)为点(i，j，k)的邻域，为点(i，j，k)以及其相邻的18个点，即是一个19点邻域模板系统：

(6)

为了避免邻域系统里出现尚未计算到的体素状态，WU等[19]采用了一个简化的只涉及已知体素的邻域系统。与二维情况类似，重构时每次计算为2个体素点(i，j，k)以及其右边的体素点(i，j+1，k)进行赋值。因此，在三维MCMC法中，对二维情况下的6点邻域模板进行延伸，将13个已知点和2个待赋值的点组成邻域模板系统，即15点邻域模板。因此，对于VLMN，就形成了一个固定的矢量马尔可夫链，其总体素为(M+1)(N+1)(L+1)，其维度为3个方向，即(M+1)(N+1)，(N+1)(L+1)，(L+1)(M+1)。

图2 利用3张两两正交的二维图像切面 重构三维数字岩心模型Fig.2 Reconstruction of 3D digital core model by using three orthogonal 2D image sections

将二维MCMC法中对原图进行遍历扫描的算法扩展到三维，用于求取三维情况下的转移概率。在有原始三维图像硬数据的情况下，可直接对原始三维图像进行遍历，获取15点邻域模板的转移概率，并用其进行新图像的重构。在没有三维图像硬数据时，对于各向同性介质，可利用一张二维图像进行三维重构；对于各向异性较强的介质(如碳酸盐岩)，则需要利用能反映3个方向转移概率的3个相互正交的原始二维图像进行图像遍历统计，以获得3个方向的二维邻域模板的转移概率，并用其进行赋值重构(见图2)。基于MCMC法重构三维数字岩心模型的流程如图3所示[2]。在图像边界时，需采用包含较少点的邻域模板。

图3 基于MCMC法重构三维数字岩心模型的流程图Fig.3 Flow chart of reconstruction of 3D digital core model based on MCMC

2 碳酸盐岩三维数字岩心模型重构

2.1 岩心原始二维切面信息分析

碳酸盐岩岩心在各个平面上的原始二维切面是由三维CT扫描图像直接截取获得。通过图像阈值分割，得到岩心的二值化图像(见图4)。

图4 碳酸盐岩岩心原始二维图像切面Fig.4 Original 2D image sections of the carbonate rock core

图4中的白色部分表示孔隙(1)，黑色部分表示岩石骨架(0)。图4中的3张图片是碳酸盐岩岩心两两正交的3个切面，主要用来描述碳酸盐岩各方向上的孔隙特征，各图像尺寸均为400×400像素，3个切面的孔隙度分别为0.212925、0.159969、0.278967，3个切面的平均孔隙度为0.217287。由图4可知，原始图片所展示的碳酸盐岩各向异性较强，孔隙度分布不均匀。

2.2 三维数字岩心模型重构

利用3张两两正交的二维切面图像信息，分别重构图像尺寸为50×50×50体素、100×100×100体素、200×200×200体素、400×400×400体素的三维数字岩心模型(见图5)。

图5 不同图像尺寸的三维数字岩心模型Fig.5 3D digital core models with different image sizes

3 重构结果分析

数值模拟的结果是基于重建的三维数字岩心模型，因此三维数字岩心模型的准确性直接影响最终岩石物理实验模拟的结果。MCMC法具有随机性，其构建的三维数字岩心模型受3个方向的条件概率影响，且其确定马尔科夫链转移概率的邻域模板的长度最多仅有6个点，因此其是否能在较大的图像尺寸的重构中保持原图的孔隙结构值得研究。

笔者利用图像的自相关函数[25]和变差函数对重构的三维数字岩心模型与原始二维切片进行对比和分析。

在重构的碳酸盐岩三维数字岩心模型上选取一个切面，利用自相关函数和变差函数与原始二维图像切面进行对比。在对比前，需对选取的重构三维数字岩心模型的切面进行扩展(见图6)，使其图像尺寸均扩展为400×400×400体素。

图6 重构三维数字岩心模型的切面及扩展Fig.6 The section and extension of the reconstructed 3D digital core model

由原始二维图像切面与重构三维数字岩心模型切面的自相关函数与变差函数对比图(见图7)可以看出，图像尺寸较大时，如200×200×200体素和400×400×400体素，重构三维数字岩心模型切面孔隙度与原始二维图像切面孔隙度差异较大；图像尺寸较小时，如50×50×50体素和100×100×100体素，重构三维数字岩心模型切面孔隙度与原始二维图像切面孔隙度符合度较高，且连通性更好。

图7 原始二维图像切面与重构三维数字岩心模型切面的自相关函数与变差函数对比(部分曲线)Fig.7 Comparison of autocorrelation function and variogram between the original 2D image sections and reconstructed 3D digital core model sections(part of the curves)

4 结语

数字岩心模型是数字岩石物理模拟研究的基础，适用范围较广的MCMC法是建立三维数字岩心模型的常用方法之一。笔者采用三维MCMC法，利用碳酸盐岩储层岩石的二维二值化图像，重构出不同图像尺寸的三维数字岩心模型。通过对重构结果进行分析可知，由于MCMC法中选取用于计算转移概率的邻域模板的长度有限，重构的三维数字岩心模型的准确性受图像尺寸影响较大，建模时选取较小的图像尺寸，建模结果较为理想，如果图像尺寸过大，建模效果不佳。因此，利用MCMC法重构三维数字岩心模型时，应选择合适的图像尺寸，尽可能建立最接近真实岩心的三维数字岩心模型。