基于馬爾可夫鏈-蒙特卡洛法的碳酸鹽巖三維數(shù)字巖心模型重構(gòu)

聶昕，胡俊強(qiáng),萬(wàn)宇，張超謨，張占松

1.油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長(zhǎng)江大學(xué))，湖北 武漢 430100 2.非常規(guī)油氣省部共建協(xié)同創(chuàng)新中心(長(zhǎng)江大學(xué))，湖北 武漢 430100 3.中國(guó)石油集團(tuán)測(cè)井有限公司，陜西 西安 710077

數(shù)字巖心不僅能夠反映儲(chǔ)層巖石復(fù)雜的孔隙空間，還可以展示巖石的固體骨架特征[1]。巖石各類(lèi)物理性質(zhì)(如彈性、導(dǎo)電性等)的數(shù)值模擬是以三維數(shù)字巖心為基礎(chǔ)進(jìn)行的，因此三維數(shù)字巖心模型的準(zhǔn)確性直接影響最終巖石物理實(shí)驗(yàn)?zāi)M的結(jié)果[2]。目前常用的三維數(shù)字巖心建模方法主要分為過(guò)程法和隨機(jī)法兩大類(lèi)[3,4]。過(guò)程法是模擬顆粒沉積、壓實(shí)等過(guò)程進(jìn)行建模，適合于普通砂巖儲(chǔ)層；隨機(jī)法是基于已經(jīng)獲得的巖石二維圖像，以其信息統(tǒng)計(jì)特性為約束，利用數(shù)學(xué)算法進(jìn)行三維數(shù)字巖心的構(gòu)建，使其與原始的巖石二維圖像的統(tǒng)計(jì)特性較為接近，適用范圍更廣泛。常用的隨機(jī)法有完全隨機(jī)法[4]、高斯場(chǎng)法[5-8]、順序指示模擬法[9-11]、模擬退火算法[12-15]、多點(diǎn)統(tǒng)計(jì)法[16，17]以及馬爾可夫鏈-蒙特卡洛法(MCMC法)[18-21]等。MCMC法自2004年被Wu等[18]引進(jìn)到數(shù)字巖心重構(gòu)以來(lái)，因其具有計(jì)算速度快、適用范圍廣及可體現(xiàn)各向異性等優(yōu)勢(shì)被廣泛用于重構(gòu)三維數(shù)字巖心，在碳酸鹽巖、頁(yè)巖等儲(chǔ)層的數(shù)字巖心建模方面有所應(yīng)用。王晨晨等[21]基于MCMC法使用疊加法重構(gòu)出能同時(shí)描述大孔隙和微孔隙性質(zhì)的碳酸鹽巖雙孔隙三維數(shù)字巖心，并對(duì)建立的數(shù)字巖心模型的孔喉特征進(jìn)行了分析，結(jié)果表明MCMC法能夠基于二維真巖心薄片快速重構(gòu)出巖石的三維數(shù)字巖心，且利用疊加法構(gòu)建的碳酸鹽巖雙孔隙數(shù)字巖心具有較高的孔隙連通體積比，孔隙之間的連通性相對(duì)于單一尺度的建模結(jié)果顯著提高。聶昕等[22]選用MCMC法并加以改進(jìn)，分別構(gòu)建頁(yè)巖氣儲(chǔ)層各微觀組分(如黏土礦物、孔隙、有機(jī)質(zhì)等)的三維模型，將上述模型疊合嵌套，得到多礦物組分的頁(yè)巖儲(chǔ)層三維數(shù)字巖心，對(duì)巖心模型分析的結(jié)果表明，利用MCMC法和嵌套法重構(gòu)的三維數(shù)字巖心可用于頁(yè)巖氣儲(chǔ)層巖石物理性質(zhì)的數(shù)值模擬研究，同時(shí)也提出了重構(gòu)方法中存在的問(wèn)題和不足，對(duì)其他巖石儲(chǔ)層數(shù)字巖心重構(gòu)方法的改進(jìn)提出了建議。目前，針對(duì)MCMC法建模的研究?jī)H限于進(jìn)行巖心的重建，然后選擇合適的尺寸進(jìn)行重構(gòu)，再研究巖石儲(chǔ)層的相關(guān)性質(zhì)。然而，在常用的基于MCMC法的三維數(shù)字巖心建模方法中，由于確定馬爾科夫鍵轉(zhuǎn)移概率的鄰域模板尺寸有限，僅有幾個(gè)像素，因此建模結(jié)果的尺寸(即邊長(zhǎng)像素?cái)?shù))的選取直接影響到了建模結(jié)果的準(zhǔn)確性，不同尺寸的重構(gòu)結(jié)果差異巨大，而前人未對(duì)該現(xiàn)象進(jìn)行過(guò)研究和總結(jié)。為此，筆者以碳酸鹽巖為例，基于二維切片信息，利用MCMC法重構(gòu)其不同像素?cái)?shù)尺寸的三維數(shù)字巖心，并對(duì)其結(jié)果進(jìn)行孔隙結(jié)構(gòu)特性分析，從而討論MCMC法重構(gòu)三維數(shù)字巖心模型的準(zhǔn)確性。

1 MCMC法原理

MCMC法來(lái)自于圖像處理過(guò)程中馬爾可夫隨機(jī)場(chǎng)(MRF)的使用。馬爾可夫鏈?zhǔn)侵敢环N狀態(tài)序列，該序列每個(gè)位置的狀態(tài)均取決于該位置前面有限個(gè)位置的狀態(tài)，而存在這個(gè)狀態(tài)的概率稱(chēng)為轉(zhuǎn)移概率[23，24]。利用MCMC法重構(gòu)圖像時(shí)，先利用馬爾可夫鏈獲得原始圖像的轉(zhuǎn)移概率，然后利用轉(zhuǎn)移概率進(jìn)行新圖像賦值重構(gòu)。

1.1 二維MCMC法原理

二維圖像的本質(zhì)是一個(gè)矩陣，圖像的高和寬為矩陣的行和列。一般情況下，采用二值化的圖像來(lái)表示僅包括孔隙和骨架信息的巖心，即巖石孔隙和巖石骨架分別用1和0來(lái)表示，因此該圖像矩陣的每個(gè)位置只有1和0兩種可能性。圖像的結(jié)構(gòu)特征可以用概率分布函數(shù)進(jìn)行表征。設(shè)一個(gè)二維圖像的像素?cái)?shù)為n，x=(x1，x2，…，xn)表示該像素處的值，在二值化的巖心圖像中，xi只能為0(骨架)或1(孔隙)?；谠紙D像，利用MCMC法重構(gòu)一個(gè)具有相近性質(zhì)的二維多孔介質(zhì)的圖像，理論上需要獲取一個(gè)完整的x的概率分布函數(shù)，即p(x)。傳統(tǒng)馬爾可夫鏈需要從樣本空間中進(jìn)行完美隨機(jī)抽樣來(lái)獲得概率分布函數(shù)，但是由于數(shù)字巖心的樣本空間非常大，導(dǎo)致該完美隨機(jī)抽樣方案幾乎不可能實(shí)現(xiàn)。因此，WU等[18]引入鄰域模板的概念，并利用其統(tǒng)計(jì)結(jié)果來(lái)代替求取概率分布函，其本質(zhì)是將模型假設(shè)成圖像中任意一點(diǎn)的狀態(tài)由該點(diǎn)附近區(qū)域極少數(shù)點(diǎn)的狀態(tài)決定，即對(duì)于一個(gè)點(diǎn)s，用Λ-s表示該圖像中s點(diǎn)以外的其他所有點(diǎn)，則應(yīng)存在一個(gè)滿(mǎn)足式(1)的s點(diǎn)的鄰近區(qū)域集合Ns：

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

(1)

通常情況下，想要準(zhǔn)確獲得式(1)中的Ns需要大量的計(jì)算，不利于實(shí)際應(yīng)用。WU等[18]在利用MCMC法重構(gòu)二維巖心圖像時(shí)，定義了影響某1點(diǎn)狀態(tài)的鄰域?yàn)?個(gè)點(diǎn)，即該點(diǎn)上面3個(gè)點(diǎn)(正上方、左上方和右上方)以及其左邊相鄰的1個(gè)點(diǎn)，而該點(diǎn)右邊相鄰的點(diǎn)同時(shí)受到包括這個(gè)點(diǎn)和其鄰域4個(gè)點(diǎn)的影響，這樣就構(gòu)成了一個(gè)5-6點(diǎn)鄰域模板系統(tǒng)(見(jiàn)圖1(a))，即：

(2)

式中：(i，j)為該圖像中位于第i行第j列的體素。

圖1 5-6點(diǎn)鄰域和2點(diǎn)鄰域模板系統(tǒng)示意圖Fig.1 Schematic diagram of formwork system of 5-6 point and 2 point neighborhood

在圖像最上面邊界的點(diǎn)，由于其上方?jīng)]有點(diǎn)，則簡(jiǎn)化為該點(diǎn)只受其左邊點(diǎn)的影響，即2點(diǎn)鄰域模板系統(tǒng)，如圖1(b)所示。在重構(gòu)之前，先進(jìn)行原始圖像的遍歷，獲得各個(gè)鄰域模板系統(tǒng)狀態(tài)情況下的轉(zhuǎn)移概率，再利用該轉(zhuǎn)移概率，采用蒙特卡洛方法對(duì)新的圖像進(jìn)行構(gòu)建。

1.2 三維MCMC法原理

以二維MCMC法原理為基礎(chǔ)，將其拓展到三維的情況。設(shè)三維數(shù)字巖心圖像為L(zhǎng)行、M列和N層的三維樣本空間VLMN，點(diǎn)(i，j，k)為該圖像中位于第i行j列k層的體素，x(Vijk)為點(diǎn)(i，j，k)處的狀態(tài)。三維與二維一樣，均采用對(duì)原始圖像進(jìn)行遍歷掃描的算法來(lái)計(jì)算轉(zhuǎn)移概率。在重構(gòu)賦值時(shí)，利用轉(zhuǎn)移概率，采用蒙特卡洛方法逐行、逐列、逐層地對(duì)體素點(diǎn)進(jìn)行賦值。第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個(gè)方向上條件概率的共同作用。

對(duì)于任意屬于VLMN的點(diǎn)(i，j，k)，聯(lián)合概率函數(shù)為：

(4)

根據(jù)式(4)，三維馬爾可夫鏈的轉(zhuǎn)移概率為：

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

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

(5)

式中:N(ijk)為點(diǎn)(i，j，k)的鄰域，為點(diǎn)(i，j，k)以及其相鄰的18個(gè)點(diǎn)，即是一個(gè)19點(diǎn)鄰域模板系統(tǒng)：

(6)

為了避免鄰域系統(tǒng)里出現(xiàn)尚未計(jì)算到的體素狀態(tài)，WU等[19]采用了一個(gè)簡(jiǎn)化的只涉及已知體素的鄰域系統(tǒng)。與二維情況類(lèi)似，重構(gòu)時(shí)每次計(jì)算為2個(gè)體素點(diǎn)(i，j，k)以及其右邊的體素點(diǎn)(i，j+1，k)進(jìn)行賦值。因此，在三維MCMC法中，對(duì)二維情況下的6點(diǎn)鄰域模板進(jìn)行延伸，將13個(gè)已知點(diǎn)和2個(gè)待賦值的點(diǎn)組成鄰域模板系統(tǒng)，即15點(diǎn)鄰域模板。因此，對(duì)于VLMN，就形成了一個(gè)固定的矢量馬爾可夫鏈，其總體素為(M+1)(N+1)(L+1)，其維度為3個(gè)方向，即(M+1)(N+1)，(N+1)(L+1)，(L+1)(M+1)。

圖2 利用3張兩兩正交的二維圖像切面 重構(gòu)三維數(shù)字巖心模型Fig.2 Reconstruction of 3D digital core model by using three orthogonal 2D image sections

將二維MCMC法中對(duì)原圖進(jìn)行遍歷掃描的算法擴(kuò)展到三維，用于求取三維情況下的轉(zhuǎn)移概率。在有原始三維圖像硬數(shù)據(jù)的情況下，可直接對(duì)原始三維圖像進(jìn)行遍歷，獲取15點(diǎn)鄰域模板的轉(zhuǎn)移概率，并用其進(jìn)行新圖像的重構(gòu)。在沒(méi)有三維圖像硬數(shù)據(jù)時(shí)，對(duì)于各向同性介質(zhì)，可利用一張二維圖像進(jìn)行三維重構(gòu)；對(duì)于各向異性較強(qiáng)的介質(zhì)(如碳酸鹽巖)，則需要利用能反映3個(gè)方向轉(zhuǎn)移概率的3個(gè)相互正交的原始二維圖像進(jìn)行圖像遍歷統(tǒng)計(jì)，以獲得3個(gè)方向的二維鄰域模板的轉(zhuǎn)移概率，并用其進(jìn)行賦值重構(gòu)(見(jiàn)圖2)?；贛CMC法重構(gòu)三維數(shù)字巖心模型的流程如圖3所示[2]。在圖像邊界時(shí)，需采用包含較少點(diǎn)的鄰域模板。

圖3 基于MCMC法重構(gòu)三維數(shù)字巖心模型的流程圖Fig.3 Flow chart of reconstruction of 3D digital core model based on MCMC

2 碳酸鹽巖三維數(shù)字巖心模型重構(gòu)

2.1 巖心原始二維切面信息分析

碳酸鹽巖巖心在各個(gè)平面上的原始二維切面是由三維CT掃描圖像直接截取獲得。通過(guò)圖像閾值分割，得到巖心的二值化圖像(見(jiàn)圖4)。

圖4 碳酸鹽巖巖心原始二維圖像切面Fig.4 Original 2D image sections of the carbonate rock core

圖4中的白色部分表示孔隙(1)，黑色部分表示巖石骨架(0)。圖4中的3張圖片是碳酸鹽巖巖心兩兩正交的3個(gè)切面，主要用來(lái)描述碳酸鹽巖各方向上的孔隙特征，各圖像尺寸均為400×400像素，3個(gè)切面的孔隙度分別為0.212925、0.159969、0.278967，3個(gè)切面的平均孔隙度為0.217287。由圖4可知，原始圖片所展示的碳酸鹽巖各向異性較強(qiáng)，孔隙度分布不均勻。

2.2 三維數(shù)字巖心模型重構(gòu)

利用3張兩兩正交的二維切面圖像信息，分別重構(gòu)圖像尺寸為50×50×50體素、100×100×100體素、200×200×200體素、400×400×400體素的三維數(shù)字巖心模型(見(jiàn)圖5)。

圖5 不同圖像尺寸的三維數(shù)字巖心模型Fig.5 3D digital core models with different image sizes

3 重構(gòu)結(jié)果分析

數(shù)值模擬的結(jié)果是基于重建的三維數(shù)字巖心模型，因此三維數(shù)字巖心模型的準(zhǔn)確性直接影響最終巖石物理實(shí)驗(yàn)?zāi)M的結(jié)果。MCMC法具有隨機(jī)性，其構(gòu)建的三維數(shù)字巖心模型受3個(gè)方向的條件概率影響，且其確定馬爾科夫鏈轉(zhuǎn)移概率的鄰域模板的長(zhǎng)度最多僅有6個(gè)點(diǎn)，因此其是否能在較大的圖像尺寸的重構(gòu)中保持原圖的孔隙結(jié)構(gòu)值得研究。

筆者利用圖像的自相關(guān)函數(shù)[25]和變差函數(shù)對(duì)重構(gòu)的三維數(shù)字巖心模型與原始二維切片進(jìn)行對(duì)比和分析。

在重構(gòu)的碳酸鹽巖三維數(shù)字巖心模型上選取一個(gè)切面，利用自相關(guān)函數(shù)和變差函數(shù)與原始二維圖像切面進(jìn)行對(duì)比。在對(duì)比前，需對(duì)選取的重構(gòu)三維數(shù)字巖心模型的切面進(jìn)行擴(kuò)展(見(jiàn)圖6)，使其圖像尺寸均擴(kuò)展為400×400×400體素。

圖6 重構(gòu)三維數(shù)字巖心模型的切面及擴(kuò)展Fig.6 The section and extension of the reconstructed 3D digital core model

由原始二維圖像切面與重構(gòu)三維數(shù)字巖心模型切面的自相關(guān)函數(shù)與變差函數(shù)對(duì)比圖(見(jiàn)圖7)可以看出，圖像尺寸較大時(shí)，如200×200×200體素和400×400×400體素，重構(gòu)三維數(shù)字巖心模型切面孔隙度與原始二維圖像切面孔隙度差異較大；圖像尺寸較小時(shí)，如50×50×50體素和100×100×100體素，重構(gòu)三維數(shù)字巖心模型切面孔隙度與原始二維圖像切面孔隙度符合度較高，且連通性更好。

圖7 原始二維圖像切面與重構(gòu)三維數(shù)字巖心模型切面的自相關(guān)函數(shù)與變差函數(shù)對(duì)比(部分曲線(xiàn))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 結(jié)語(yǔ)

數(shù)字巖心模型是數(shù)字巖石物理模擬研究的基礎(chǔ)，適用范圍較廣的MCMC法是建立三維數(shù)字巖心模型的常用方法之一。筆者采用三維MCMC法，利用碳酸鹽巖儲(chǔ)層巖石的二維二值化圖像，重構(gòu)出不同圖像尺寸的三維數(shù)字巖心模型。通過(guò)對(duì)重構(gòu)結(jié)果進(jìn)行分析可知，由于MCMC法中選取用于計(jì)算轉(zhuǎn)移概率的鄰域模板的長(zhǎng)度有限，重構(gòu)的三維數(shù)字巖心模型的準(zhǔn)確性受圖像尺寸影響較大，建模時(shí)選取較小的圖像尺寸，建模結(jié)果較為理想，如果圖像尺寸過(guò)大，建模效果不佳。因此，利用MCMC法重構(gòu)三維數(shù)字巖心模型時(shí)，應(yīng)選擇合適的圖像尺寸，盡可能建立最接近真實(shí)巖心的三維數(shù)字巖心模型。