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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 MCMC方法圖像重建理論

2.1 二維MCMC方法理論

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

(1)

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

(2)

對于2點領(lǐng)域系統(tǒng)，定義影響的某點狀態(tài)的領(lǐng)域為該點左邊1個點；對于5點領(lǐng)域系統(tǒng)，定義影響的某點狀態(tài)的領(lǐng)域為該點上面的3個點和左邊的1個點；對于6點領(lǐng)域系統(tǒng)，定義影響的某點狀態(tài)的領(lǐng)域為該點上面的3個點和左邊的2個點。如圖1所示，其表達式分別為

(3)

(4)

(5)

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

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

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

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

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

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

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

流程如圖2所示。

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

2.2 MCMC方法驗證

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

圖3 應(yī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)程度，對于只考慮孔隙和巖石骨架的兩相系統(tǒng)，自相關(guān)函數(shù)可以當作系統(tǒng)中兩點分布于同一相中的概率，自相關(guān)函數(shù)表達式為

(7)

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

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

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

3 頁巖儲層數(shù)字巖心重建

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

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

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

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

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

3.3 模型的組合

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

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

(2)將不同組分進行相加或者相乘的運算組合起來，獲得最終的總模型。

圖8 頁巖儲層數(shù)字巖心

組合表達式為

Π=Π1+Π2，

(8)

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

(9)

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

4 結(jié) 論

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

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

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

[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] 施尚明,李照永.蒙特卡洛法在吉祥屯氣藏儲量計算中的應(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é)報：自然科學(xué)版,2012,27(4)：54-61.WANGBo,NINGZheng-fu,JIJiang.Studyon3Dreconstructionmethodsofporousmediummodel[J].JournalofXi'anShiyouUniversity:NaturalScienceEdition,2012,27(4):54-61.

[14] 李仁民，劉松玉，方磊,等.采用隨機生長四參數(shù)生成法構(gòu)造黏土微觀結(jié)構(gòu)[J].浙江大學(xué)學(xué)報：工學(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

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

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

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

TE319

A