張思勤，汪志明，王小秋，李江濤，洪 凱

(中國(guó)石油大學(xué)(北京) 石油工程學(xué)院，北京 102249)

基于MCMC的數(shù)字巖心重建方法

張思勤，汪志明，王小秋，李江濤，洪 凱

(中國(guó)石油大學(xué)(北京) 石油工程學(xué)院，北京 102249)

常規(guī)的多孔介質(zhì)重建方法主要考慮基質(zhì)和孔隙兩部分，而且編程難度大，計(jì)算速度慢。針對(duì)頁(yè)巖氣儲(chǔ)層各向異性且含有有機(jī)質(zhì)的特點(diǎn)，提出了一種基于改進(jìn)后的馬爾科夫鏈-蒙特卡洛(MCMC)方法的數(shù)字巖心重建方法：先從實(shí)際頁(yè)巖氣儲(chǔ)層二維切片圖像中提取孔隙度和有機(jī)質(zhì)，分別重建孔隙和有機(jī)質(zhì)的數(shù)字巖心，再將孔隙和有機(jī)質(zhì)的初始模型組合在一起構(gòu)成最終的頁(yè)巖氣儲(chǔ)層數(shù)字巖心。研究結(jié)果表明：通過(guò)對(duì)比分析二維MCMC重建圖像和原圖自相關(guān)函數(shù)，發(fā)現(xiàn)該方法重建圖像能夠較好反映原圖性質(zhì)；該方法計(jì)算速度快，過(guò)程簡(jiǎn)單，方法所重建的模型能夠較好體現(xiàn)頁(yè)巖儲(chǔ)層各向異性的特點(diǎn)且能反映頁(yè)巖儲(chǔ)層含有有機(jī)質(zhì)的特性。

頁(yè)巖氣儲(chǔ)層； 馬爾科夫鏈-蒙特卡洛；數(shù)字巖心；組合模型

多孔介質(zhì)結(jié)構(gòu)復(fù)雜，普遍存在于自然界中，因此，對(duì)于多孔介質(zhì)的研究具有廣泛的工程意義，如煤層氣的開(kāi)采、頁(yè)巖氣的開(kāi)發(fā)等。在宏觀滲流中，一般應(yīng)用達(dá)西定律來(lái)描述流體運(yùn)動(dòng)特征，由達(dá)西定律得到的宏觀參數(shù)只能用來(lái)表征流體運(yùn)動(dòng)的宏觀特性，而對(duì)于含有納米尺度孔隙的頁(yè)巖氣儲(chǔ)層，無(wú)法應(yīng)用達(dá)西定律進(jìn)行描述，因此，從微觀角度研究多孔介質(zhì)，并建立一套可以描述多孔介質(zhì)特性的模型，從而為進(jìn)一步研究多孔介質(zhì)滲流機(jī)理提供一個(gè)平臺(tái)，具有十分重要的工程價(jià)值和理論指導(dǎo)意義。

多孔介質(zhì)模型主要包括毛管模型[1]、球形顆粒堆積模型、格子模型和數(shù)字巖心網(wǎng)絡(luò)模型,主要分為3類(lèi)：概念模型、統(tǒng)計(jì)模型和圖像重建模型。隨著計(jì)算機(jī)技術(shù)的提高和高分辨率圖像掃描電鏡的發(fā)展，許多學(xué)者提出了不同的方法來(lái)重建多孔介質(zhì)模型，重建的數(shù)字巖心模型也越來(lái)越接近真實(shí)地層。本文在Wu Kejian等學(xué)者[2-3]研究成果的基礎(chǔ)上，結(jié)合MATLAB圖像處理技術(shù)和MCMC方法重建頁(yè)巖氣儲(chǔ)層微觀結(jié)構(gòu)，為進(jìn)一步從微觀角度研究頁(yè)巖氣儲(chǔ)層滲流機(jī)理提供基礎(chǔ)。

1 多孔介質(zhì)模型重建方法

數(shù)字巖心重建方法主要分為2類(lèi)：物理實(shí)驗(yàn)方法和數(shù)值重建方法。物理實(shí)驗(yàn)方法主要有CT掃描法、序列成像組合法和 FIB-SEM法。數(shù)值重建方法發(fā)展較快，主要包括隨機(jī)法和過(guò)程法。

CT掃描技術(shù)首先應(yīng)用于醫(yī)學(xué)領(lǐng)域，Dunsmuir等[4-5]將CT掃描技術(shù)應(yīng)用到石油領(lǐng)域，并通過(guò)CT掃描技術(shù)和計(jì)算機(jī)處理技術(shù)，得到了巖心的三維圖像。但CT掃描儀造價(jià)高，試驗(yàn)成本大，不便于大規(guī)模應(yīng)用。

20世紀(jì)70年代后,隨著計(jì)算機(jī)技術(shù)的發(fā)展，隨機(jī)方法也發(fā)展較快，國(guó)內(nèi)外學(xué)者提出不同的數(shù)值重建方法。迄今為止，發(fā)展的隨機(jī)方法主要包括完全隨機(jī)法、隨機(jī)生長(zhǎng)法、高斯模擬法、模擬退火法、多點(diǎn)統(tǒng)計(jì)法和MCMC法。完全隨機(jī)方法和隨機(jī)生長(zhǎng)方法比較簡(jiǎn)單，一般只能用于理論層面的模擬。Joshi[6]于1974年提出高斯模擬法，該方法使用孔隙度和兩點(diǎn)相關(guān)函數(shù)作為約束條件，Quiblier[7]在Joshi的基礎(chǔ)上，建立了第一個(gè)三維數(shù)字巖心，之后很多學(xué)者在其基礎(chǔ)上進(jìn)行了改進(jìn)。Hazlett[6]于1997年提出模擬退火法。

不同于隨機(jī)建模方法，Bryant、Bakke和Oren等[9-10]提出了通過(guò)模擬真實(shí)地層的沉積、壓實(shí)和成巖3個(gè)過(guò)程來(lái)建立數(shù)字巖心的方法，并運(yùn)用此方法建立了Fontainebleau砂巖的數(shù)字巖心。與其他建模方法進(jìn)行對(duì)比發(fā)現(xiàn)，該模型最大的特點(diǎn)是可以很好地體現(xiàn)真實(shí)地層的傳導(dǎo)性能。此后，許多學(xué)者又進(jìn)一步進(jìn)行了研究，Oren和Bakke[11]于2003年綜合了隨機(jī)法和過(guò)程法的特點(diǎn)，建立了Berea砂巖的孔隙空間數(shù)字巖心模型。

綜上所述，通過(guò)物理實(shí)驗(yàn)方法建立數(shù)字巖心，可以比較準(zhǔn)確地反映巖心的微觀結(jié)構(gòu)，但在實(shí)際應(yīng)用中，由于物理實(shí)驗(yàn)方法成本高，過(guò)程繁瑣，耗時(shí)較長(zhǎng)，不便于大規(guī)模應(yīng)用。而通過(guò)掃描電鏡技術(shù)，巖心的二維切片圖像較易獲得，因而通過(guò)圖像處理技術(shù)得到二維切片圖像信息并配合數(shù)值方法重建數(shù)字巖心更適用于理論研究。由于頁(yè)巖儲(chǔ)層結(jié)構(gòu)復(fù)雜、孔隙類(lèi)型多樣且孔隙小、含有有機(jī)質(zhì)，建模條件比較復(fù)雜，需要考慮的因素較多，綜合考慮建立數(shù)字巖心的時(shí)間、成本以及適用條件，優(yōu)選出MCMC方法來(lái)建立頁(yè)巖氣儲(chǔ)層數(shù)字巖心。

2 MCMC方法圖像重建理論

2.1 二維MCMC方法理論

MCMC法重建圖像的主要思路[2-3,12-13]是利用馬爾可夫鏈思想獲得轉(zhuǎn)移概率，再利用轉(zhuǎn)移概率進(jìn)行賦值重建。首先，將二維圖形看作一個(gè)矩陣，將頁(yè)巖巖心圖像進(jìn)行二值化，每個(gè)點(diǎn)只有0和1兩種狀態(tài)，小于閾值部分代表巖石骨架和大于閾值部分代表巖石孔隙，用公式表示為

(1)

引入領(lǐng)域的思想，認(rèn)為模型中任何點(diǎn)的狀態(tài)只取決于相鄰少數(shù)點(diǎn)的狀態(tài)。具體地說(shuō)，就是對(duì)于一個(gè)特定的點(diǎn)s，用Λ-s表示除s點(diǎn)外的所有點(diǎn)。則存在一個(gè)s點(diǎn)的領(lǐng)域Ns，有

(2)

對(duì)于2點(diǎn)領(lǐng)域系統(tǒng)，定義影響的某點(diǎn)狀態(tài)的領(lǐng)域?yàn)樵擖c(diǎn)左邊1個(gè)點(diǎn)；對(duì)于5點(diǎn)領(lǐng)域系統(tǒng)，定義影響的某點(diǎn)狀態(tài)的領(lǐng)域?yàn)樵擖c(diǎn)上面的3個(gè)點(diǎn)和左邊的1個(gè)點(diǎn)；對(duì)于6點(diǎn)領(lǐng)域系統(tǒng)，定義影響的某點(diǎn)狀態(tài)的領(lǐng)域?yàn)樵擖c(diǎn)上面的3個(gè)點(diǎn)和左邊的2個(gè)點(diǎn)。如圖1所示，其表達(dá)式分別為

(3)

(4)

(5)

圖1 2點(diǎn)、 5點(diǎn)和6點(diǎn)鄰域系統(tǒng)示意圖

獲得條件概率后，就可以利用其為像素賦值，進(jìn)行圖像的重建。重建的步驟為：

(1)利用孔隙度確定(1,1)點(diǎn)的狀態(tài)，然后利用2點(diǎn)領(lǐng)域模板的條件概率，從左向右依次為(1，j)點(diǎn)賦值。

(2)利用步驟(1)中(1,1)點(diǎn)的狀態(tài)，利用2點(diǎn)領(lǐng)域模板的條件概率，從上往下依次為(i,1)點(diǎn)賦值。

(3)利用步驟(2)中(1,1)、(1,2)和(2,1)點(diǎn)的狀態(tài)，采用4點(diǎn)領(lǐng)域模板的條件概率，對(duì)(2,2)賦值。同理，從上往下依次為點(diǎn)(i，j)賦值，其中i>2。

(4)采用6點(diǎn)領(lǐng)域模板，從第2行開(kāi)始從左向右依次為點(diǎn)(2，j)賦值。到第2行結(jié)尾時(shí)，采用同樣方法對(duì)第3行(3，j)賦值，依此類(lèi)推，直至對(duì)最后一行賦值，圖像重建結(jié)束。

(5)對(duì)比重構(gòu)圖像與原圖的孔隙度，如果達(dá)到要求就終止計(jì)算，輸出最終結(jié)果；如果不符合孔隙度要求，則調(diào)整條件概率的加權(quán)因子，再次開(kāi)始重構(gòu)過(guò)程，直到符合要求為止。

流程如圖2所示。

圖2 二維MCMC法重建數(shù)字巖心流程

2.2 MCMC方法驗(yàn)證

為了驗(yàn)證MCMC方法是否能反映頁(yè)巖儲(chǔ)層各項(xiàng)異性特點(diǎn)，首先采用四參數(shù)生長(zhǎng)方法(QSGS)構(gòu)造幾幅各向異性層狀介質(zhì)圖像，四參數(shù)隨機(jī)生長(zhǎng)方法的控制參數(shù)包括孔隙度φ、固相分布概率Pc、方向概率Pi(i=1～8)和概率密度。首先通過(guò)調(diào)整4個(gè)參數(shù)中某個(gè)參數(shù)或某些參數(shù)，得到3種不同形貌的多孔介質(zhì)，如圖3所示；其次，對(duì)應(yīng)于隨機(jī)生長(zhǎng)方法和圖像處理方法產(chǎn)生的多孔介質(zhì)，采用MCMC方法分別重建3種不同多孔介質(zhì)，所得到的重建圖像如圖4所示，對(duì)比重建圖像和原圖， 可以發(fā)現(xiàn)重建圖像能很好代表原圖的特點(diǎn)。為了進(jìn)一步驗(yàn)證二維MCMC方法的可靠性， 在此引入2個(gè)統(tǒng)計(jì)函數(shù)單點(diǎn)

圖3 應(yīng)用隨機(jī)生長(zhǎng)方法和圖像處理方法產(chǎn)生的多孔介質(zhì)(白色為孔隙，黑色為基質(zhì))

圖4 應(yīng)用MCMC方法重建的多孔介質(zhì)(白色為孔隙，黑色為基質(zhì))

概率函數(shù)和自相關(guān)函數(shù)作為約束條件。

(6)

自相關(guān)函數(shù)表示相函數(shù)在不同位置r和r+h的取值之間的相關(guān)程度，對(duì)于只考慮孔隙和巖石骨架的兩相系統(tǒng)，自相關(guān)函數(shù)可以當(dāng)作系統(tǒng)中兩點(diǎn)分布于同一相中的概率，自相關(guān)函數(shù)表達(dá)式為

(7)

式中，r和r+h為系統(tǒng)中相距h的任意兩點(diǎn)。為了計(jì)算方便，令r的單位為像素?cái)?shù)，像素與像素邊長(zhǎng)的乘積即物理長(zhǎng)度。

自相關(guān)函數(shù)是與圖像的結(jié)構(gòu)相關(guān)的函數(shù)，也是評(píng)價(jià)圖像結(jié)構(gòu)性質(zhì)的重要函數(shù)。以圖3中前2種多孔介質(zhì)為例，如圖3(a)和圖3(b)所示，對(duì)比分析2幅多孔介質(zhì)的重建圖和原圖的自相關(guān)函數(shù)，其結(jié)果如圖5所示，從中可以看出，原圖和重建圖像的自相關(guān)函數(shù)吻合性較好，從而說(shuō)明重建圖像可以較好地反映原圖性質(zhì)，即驗(yàn)證了MCMC方法的準(zhǔn)確性。

圖5 2種多孔介質(zhì)的重建圖與原圖自相關(guān)函數(shù)對(duì)比

3 頁(yè)巖儲(chǔ)層數(shù)字巖心重建

3.1 孔隙和有機(jī)質(zhì)二值化

根據(jù)MCMC方法的需要，首先選取具有代表性的頁(yè)巖二維切片圖像(圖6)；根據(jù)建模需要，利用二維頁(yè)巖巖樣切片圖，用MATLAB軟件對(duì)孔隙和有機(jī)質(zhì)分別進(jìn)行二值化，并在二值的過(guò)程中對(duì)圖像進(jìn)行遍歷掃描，基質(zhì)為0(黑色)，孔隙為1(白色)，有機(jī)質(zhì)為2(紅色)，從而可以獲得孔隙相對(duì)于其他組分的分布特征，以便下一步進(jìn)行三維數(shù)字巖心的重建。提取孔隙和有機(jī)質(zhì)時(shí)采用雙閾值法[111，134]，有機(jī)質(zhì)和孔隙空間的圖像如圖7所示，孔隙度為φ0=0.018 1，有機(jī)質(zhì)含量為φc=0.036 24。

圖6 某頁(yè)巖巖樣掃描電鏡圖像 圖7 有機(jī)質(zhì)和孔隙空間二值圖像

3.2 數(shù)字巖心模型的建立

獲取孔隙和有機(jī)質(zhì)的二值圖像后，應(yīng)用MCMC方法分別對(duì)其進(jìn)行遍歷掃描及二維重建，首先利用MCMC方法對(duì)頁(yè)巖氣儲(chǔ)層孔隙空間進(jìn)行重建，重建的孔隙空間模型如圖8(a)所示，重建圖像中孔隙度(白色)φ0=0.020 5，圖中黑色部分代表巖石骨架，白色代表孔隙空間，可以看出，孔隙空間的分布與二維圖像特征吻合，多為小而分散的孔隙；重建的有機(jī)質(zhì)模型如圖8(b)所示，重建圖像中有機(jī)質(zhì)(紅色)含量φc=0.036 38；然后將孔隙二維模型和有機(jī)質(zhì)的二維模型進(jìn)行組合，從而得到完整的頁(yè)巖氣儲(chǔ)層數(shù)字巖心，如圖8(c)所示，該巖心的孔隙度φ0=0.018 5，有機(jī)質(zhì)含量為φc=0.035 79。

3.3 模型的組合

用二值圖像對(duì)孔隙和有機(jī)質(zhì)進(jìn)行重建。獲得初始模型后將孔隙和有機(jī)質(zhì)組合成一個(gè)完整的數(shù)字巖心，不同的值采用不同的方法，其具體思路如下：

(1)將不同的初始模型的組分分別賦于不同的值，將孔隙模型中的孔隙和骨架值分別設(shè)為0和1，有機(jī)質(zhì)模型中有機(jī)質(zhì)和骨架分別設(shè)為2和0；

(2)將不同組分進(jìn)行相加或者相乘的運(yùn)算組合起來(lái)，獲得最終的總模型。

圖8 頁(yè)巖儲(chǔ)層數(shù)字巖心

組合表達(dá)式為

Π=Π1+Π2，

(8)

式中：П為總模型參數(shù)；П1為孔隙模型參數(shù)；П2為有機(jī)質(zhì)模型參數(shù)?？紫赌Ｐ椭邪|(zhì)和孔隙，有機(jī)質(zhì)模型中包含基質(zhì)和有機(jī)質(zhì)2部分，故其組合類(lèi)型有4種：

(9)

在組合過(guò)程中，由于先放進(jìn)模型的組分會(huì)有部分被后放進(jìn)去的所覆蓋，所以在實(shí)際建模過(guò)程中，應(yīng)該結(jié)合各組分含量的特點(diǎn)，根據(jù)各種組分的重要性，設(shè)定優(yōu)先級(jí)。在組合過(guò)程中，設(shè)定有機(jī)質(zhì)的優(yōu)先級(jí)最高，其次為孔隙度，組合完成后，得到的三維數(shù)字巖心如圖8所示，該模型的組分特征為：孔隙度(白色)，有機(jī)質(zhì)(紅色)，與原始圖像吻合性較好。

4 結(jié) 論

(1)對(duì)比分析數(shù)值重建方法，并結(jié)合頁(yè)巖氣儲(chǔ)層孔隙小、結(jié)構(gòu)復(fù)雜、各向異性等特點(diǎn)，優(yōu)選出MCMC方法來(lái)建立數(shù)字巖心，結(jié)果表明：改進(jìn)的MCMC方法計(jì)算速度快，過(guò)程相對(duì)簡(jiǎn)單，并能夠反映頁(yè)巖氣儲(chǔ)層各向異性特點(diǎn)。

(2)利用四參數(shù)隨機(jī)生長(zhǎng)法生成3種不同性質(zhì)的多孔介質(zhì)，然后應(yīng)用MCMC方法進(jìn)行數(shù)值重建，重建的模型能較好體現(xiàn)頁(yè)巖氣儲(chǔ)層各向異性和含有有機(jī)質(zhì)的特點(diǎn)，通過(guò)分析重建圖像和原始圖像的自相關(guān)函數(shù)，發(fā)現(xiàn)重建圖像能較好反映原圖的性質(zhì)。

(3)利用頁(yè)巖氣儲(chǔ)層二維切片圖，對(duì)孔隙和有機(jī)質(zhì)分別進(jìn)行二值化，構(gòu)建各自的數(shù)字巖心，組合到一起，得到頁(yè)巖氣儲(chǔ)層數(shù)字巖心總模型，為頁(yè)巖氣儲(chǔ)層微觀滲流機(jī)理的研究提供了一個(gè)模擬平臺(tái)。

[1]FattI.ThenetworkmodelofporousmediaI,capillarypressurecharacteristics[J].TransAIME,1956,207(a):144-159.

[2]WuKejian,NaoiseNunan,JohnWC,etal.Anefficientmarkovchainmodelforthesimulationofheterogeneoussoilstructure[J].SoilScienceSocietyofAmerican,2004,68:346-351.

[3]WuKejian，MarinusIJ,VanDijke,etal.3Dstochasticofheterogeneousporousmedia-applicationstoreservoirrocks[J].TransportinPorousMedia,2006,65(3):443-467.

[4]DunsmuirJH，F(xiàn)ergusonSR,D'AmicoKL,etal.X-raymicro-tomography:anewtoolforthecharacterizationofporousmedia[C].SPE22860，Proceedingsof66thAnnualTechnicalConferenceandExhibitionoftheSocietyofPetroleumEngineers,Dallas,TX,1991.

[5]CoenenJ,TchouparovaE,JingX.MeasurementparametersandresolutionaspectsofmicroX-raytomographyforadvancedcoreanalysis[C].ProceedingsofInternationalSymposiumoftheSocietyofCoreAnalysts,2004.

[6]JoshiM.AClassofStochasticModelsforPorousMedia[D].Kalnsas:UniversityofKansas,1974.

[7]QuiblierJA.Anewthree-dimensionalmodelingtechniqueforstudyingporousmedia[J].JournalofColloidandInterfaceScience,1984,98(1):84-102.

[8]HazlettRD．Statisticalcharacterizationandstochasticmodelingofporenetworksinrelationtoeluidflow[J].MathematicalGeology,1997,29(6):801-822.

[9]BryantS,BluntM.Predictionofrelativepermeabilityinsimpleporousmedia[J].PhysicalReviewA,1992,46(4):2004-2012.

[10]BakkeS,OrenPE.3-Dpore-scalemodelingofsandstonesandflowsimulationsintheporenetworks[J].SPEJournal,1997,2(2):136-149.

[11]OrenPE,BakkeS.ReconstructionofBereasandstoneandporescalemodelingofwettabilityeffects[J].JournalofPetroleumScienceandEngineering,2003,39(2):177-199.

[12] 施尚明,李照永.蒙特卡洛法在吉祥屯氣藏儲(chǔ)量計(jì)算中的應(yīng)用[J].特種油氣藏,2005,12(2):32-34.SHIShang-ming,LIZhao-yong.ApplicationofMonteCarlomethodinJixiangtungasreservoirreservecalculation[J].SpecialOilandGasReservoirs,2005,12(2):32-34.

[13] 王波，寧正福，姬江，等.多孔介質(zhì)模型的三維重構(gòu)方法[J].西安石油大學(xué)學(xué)報(bào)：自然科學(xué)版,2012,27(4)：54-61.WANGBo,NINGZheng-fu,JIJiang.Studyon3Dreconstructionmethodsofporousmediummodel[J].JournalofXi'anShiyouUniversity:NaturalScienceEdition,2012,27(4):54-61.

[14] 李仁民，劉松玉，方磊,等.采用隨機(jī)生長(zhǎng)四參數(shù)生成法構(gòu)造黏土微觀結(jié)構(gòu)[J].浙江大學(xué)學(xué)報(bào)：工學(xué)版,2010,44(10):1897-1901.LIRen-min,LIUSong-yu,FANGLei,etal.Micro-structureofclaygeneratedbyquartetstructuregenerationset[J].JournalofZhejiangUniversity:EngineeringScience,2010,44(10):1897-1901.

責(zé)任編輯：張新寶

2015-01-16

國(guó)家自然科學(xué)基金創(chuàng)新研究群體“復(fù)雜油氣井鉆井與完井基礎(chǔ)研究”(編號(hào)：51221003)；國(guó)家自然科學(xué)聯(lián)合基金重點(diǎn)支持項(xiàng)目“頁(yè)巖氣鉆探中的井壁穩(wěn)定及高效鉆完井基礎(chǔ)研究”之課題“頁(yè)巖氣儲(chǔ)層壓裂完井優(yōu)化研究”(編號(hào)：U1262201)

張思勤(1988-)，男，碩士研究生，主要從事頁(yè)巖儲(chǔ)層數(shù)字巖心建模及滲流機(jī)理研究。 E-mail:stemenzhang@126.com

1673-064X(2015)05-0069-06

TE319

A