基于馬爾可夫鏈-蒙特卡洛法的碳酸鹽巖三維數字巖心模型重構

聶昕，胡俊強,萬宇，張超謨，張占松

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(孔隙)?；谠紙D像，利用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)?；贛CMC法重構三維數字巖心模型的流程如圖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法重構三維數字巖心模型時，應選擇合適的圖像尺寸，盡可能建立最接近真實巖心的三維數字巖心模型。