隋微波 ，尤園，程思，鄭衣珍

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

2 中國石油大學(xué)(北京) 油氣資源與探測國家重點實驗室，北京 102249

0 引言

隨著數(shù)字巖心技術(shù)的不斷發(fā)展，逐漸成為研究儲層微觀孔隙結(jié)構(gòu)﹑滲流特征參數(shù)和機理的新手段。利用數(shù)字巖心技術(shù)，我們可以準(zhǔn)確觀察目標(biāo)區(qū)域樣品的微觀結(jié)構(gòu)，在此基礎(chǔ)上求取樣品的巖石物理參數(shù)，并以數(shù)字巖心模型為流場研究流體微觀運移規(guī)律[1-2]。近年來隨著微觀成像和材料微觀結(jié)構(gòu)計算表征技術(shù)的共同進(jìn)步，對于非均質(zhì)性較強的材料進(jìn)行定量的微觀結(jié)構(gòu)表征方法也成為了研究熱點。

對于現(xiàn)階段數(shù)字巖心重構(gòu)技術(shù)的兩種主要方法來說，直接成像法利用微米級甚至納米級CT﹑聚焦離子束場發(fā)射掃描電鏡等實驗手段可以獲得超高分辨率巖心三維圖像[3-4]，結(jié)合有限元模擬等方法來預(yù)測其巖石物理特性或進(jìn)行微觀滲流和變形模擬[5-6]。直接成像法實驗一般費用昂貴﹑周期較長，獲得的高分辨率巖石三維結(jié)構(gòu)數(shù)據(jù)體非常龐大[7]，對于后期計算模擬來說包含很多必要性不強的冗余結(jié)構(gòu)信息，且僅對實驗樣品進(jìn)行表征，無法實現(xiàn)快速提取巖心樣品的統(tǒng)計特性并進(jìn)行遷移建模。因此在直接成像技術(shù)不斷進(jìn)步的今天，人們?nèi)匀恍枰l(fā)展高效穩(wěn)定的數(shù)值重構(gòu)方法來建立數(shù)字巖心模型[8]。

數(shù)值重構(gòu)法通常結(jié)合已有的掃描電鏡實驗得到的二維數(shù)字巖心圖像和各種統(tǒng)計學(xué)方法或模擬巖石的形成過程，結(jié)合優(yōu)化算法建立數(shù)字巖心。目前較為典型的數(shù)值重構(gòu)法包括模擬退火法﹑過程模擬法﹑多點地質(zhì)統(tǒng)計法和馬爾可夫鏈蒙特卡洛法。在頁巖數(shù)字巖心重構(gòu)方面，Sui等[9]提出“雙重區(qū)域重構(gòu)”模式并采用非?？焖倌M退火法重構(gòu)頁巖數(shù)字巖心；Chen等[10]基于四川盆地頁巖樣品分辨率為5 nm的場發(fā)射掃描電鏡圖像，利用馬爾可夫鏈蒙特卡羅法建立了頁巖納米級孔隙結(jié)構(gòu)的三維數(shù)字巖心；龐偉[11]以頁巖納米CT數(shù)據(jù)為訓(xùn)練圖像，采用多點地質(zhì)統(tǒng)計法重構(gòu)頁巖數(shù)字巖心；楊永飛等[12]采用模擬退火法構(gòu)建無機孔隙，采用馬爾可夫鏈蒙特卡羅法構(gòu)建有機孔隙，并將二者疊加整合獲得頁巖數(shù)字巖心；劉磊等[8]以多點統(tǒng)計學(xué)為基礎(chǔ)，考慮頁巖儲層微裂縫的影響，利用模式法進(jìn)行多尺度頁巖數(shù)字巖心重構(gòu)。以上關(guān)于頁巖數(shù)字巖心的重構(gòu)方法的研究從不同程度上考慮了頁巖儲層發(fā)育不同孔隙類型和微裂縫，多尺度效應(yīng)明顯的特征，但是對于頁巖儲層中有機質(zhì)分布的各向異性特征在重構(gòu)過程中考慮不足。

有機質(zhì)是頁巖氣儲層的重要組成部分，其類型及含量是影響頁巖氣聚集的重要因素之一[13]。泥頁巖吸附氣含量受有機質(zhì)類型﹑豐度﹑成熟度及埋藏深度等多種因素影響，但在通常情況下與有機質(zhì)豐度呈正相關(guān)關(guān)系[14]。王世謙等[15]認(rèn)為，在相同溫壓條件下，富有機質(zhì)的頁巖較貧有機質(zhì)的頁巖具有更多的微孔隙空間，能吸附更多的天然氣，因此有機質(zhì)豐度成為頁巖氣評價的首要地質(zhì)因素；姚軍等[16]從微觀滲流角度研究了頁巖有機質(zhì)孔隙對頁巖氣流動能力影響，同時提出了有機質(zhì)孔隙分布的三種模式，指出有機質(zhì)孔隙條帶狀分布時氣體流動能力最強；周楓等[17]對四川盆地龍馬溪組頁巖各向異性影響因素進(jìn)行研究，指出有機質(zhì)含量對頁巖各向異性有重要影響；Kwon等[18]通過對美國Wicox頁巖儲層樣品研究，發(fā)現(xiàn)在孔隙尺度下頁巖中的有機質(zhì)孔隙具有很大的長寬比，其結(jié)構(gòu)呈現(xiàn)出明顯的各向異性；Gu等[19]利用中子散射定量研究了頁巖有機質(zhì)孔隙尺度的各向異性；王沫然等[20]考慮頁巖有機質(zhì)各向異性特征對頁巖的微觀滲流機理進(jìn)行了研究。

綜上所述，頁巖儲層有機質(zhì)分布在孔隙尺度存在各向異性特征，且對頁巖宏觀力學(xué)﹑滲流的各向異性存在重要影響，與頁巖油氣高效開發(fā)密切相關(guān)。本文針對目前頁巖數(shù)字巖心重構(gòu)未能充分考慮有機質(zhì)分布各向異性特征的問題，提出在模擬退火隨機重構(gòu)過程中采用兩點簇函數(shù)作為有機質(zhì)三維連續(xù)性的統(tǒng)計指示函數(shù)，從而提高頁巖數(shù)字巖心重構(gòu)的準(zhǔn)確性和可靠性。

1 重構(gòu)算法統(tǒng)計函數(shù)定義

使用模擬退火法進(jìn)行巖心的數(shù)值重構(gòu)時，無法直接用圖像作為約束條件，只能由統(tǒng)計函數(shù)來描述反映真實巖心微觀結(jié)構(gòu)特征的重要信息。而巖心重構(gòu)過程的關(guān)鍵也在于所選取的統(tǒng)計函數(shù)能否有效地反映巖心微觀結(jié)構(gòu)特征。在運用傳統(tǒng)模擬退火法重構(gòu)頁巖數(shù)字巖心時，主要應(yīng)用了3 種統(tǒng)計函數(shù)：單點相關(guān)函數(shù)﹑兩點相關(guān)函數(shù)和線性路徑函數(shù)[21]，本文在以上三種統(tǒng)計函數(shù)的基礎(chǔ)上增加了另外一個微觀結(jié)構(gòu)的描述量即兩點簇函數(shù)。以上四種統(tǒng)計函數(shù)的定義分別如下文，其物理含義的圖形表示見圖1。

1.1 單點相關(guān)函數(shù)

在多相系統(tǒng)中某一i相的指示函數(shù)為Z(i)，則該i相的指示函數(shù)定義為：

其中x表示該系統(tǒng)中的任意一點；i=1, 2, 3, ...,n，其中n表示不同相的數(shù)量。該系統(tǒng)的單點相關(guān)函數(shù)φi表示第i相的體積分?jǐn)?shù)，由上式的集合平均獲得，可表示為

1.2 兩點相關(guān)函數(shù)

在上述多相系統(tǒng)中第i相和第j相的兩點相關(guān)函數(shù)定義為：

該函數(shù)表示系統(tǒng)中隨機選擇的兩點x1和x2分別分布于第i相和第j相的概率。如果系統(tǒng)中具有n個不同相則可寫出n×n個不同的兩點相關(guān)函數(shù)，其中獨立的兩點相關(guān)函數(shù)為n個，其他n×n(- 1)個兩點相關(guān)函數(shù)可由獨立的n個函數(shù)進(jìn)行表示。對于本文中研究的頁巖有機質(zhì)孔隙重構(gòu)的情況n=2，只需考慮有機質(zhì)孔隙的兩點自相關(guān)函數(shù)Sp 和孔隙—基質(zhì)之間的兩點互相關(guān)函數(shù)即可，其中p表示孔隙相，m表示基質(zhì)相。

圖1 兩點相關(guān)函數(shù)S2、線性路徑函數(shù)L2、兩點簇函數(shù)C2Fig. 1 Two-point correlation function S2, linear path function L2, two-point cluster function C2

兩點相關(guān)函數(shù)只與兩點之間的相對位移矢量有關(guān)，即

其中r=x2-x1為相對位移矢量。r=0 時兩點自相關(guān)函數(shù)退化為單點相關(guān)函數(shù)，表示所選點分布于相應(yīng)相中的概率即相應(yīng)相的體積分?jǐn)?shù)；由于某一點不可能同時分布于兩相中，因此r=0 時兩點互相關(guān)函數(shù)值為零。當(dāng)r取值很大時，兩點位于相應(yīng)相中的概率相關(guān)性很小，因此兩點自相關(guān)函數(shù)趨近于，而兩點互相關(guān)函數(shù)趨近于

對于統(tǒng)計意義上的各向同性材料，S2只與兩點間距離相關(guān)，因此無需考慮方向性對S2的影響。對于本文中研究的頁巖有機質(zhì)孔隙來說，實驗數(shù)據(jù)表明其分布具有各向異性特征，因此我們同時采用三個正交方向的S2函數(shù)進(jìn)行重構(gòu)。三個正交方向依次為平行層理面的方向(XY平面)和垂直層理面兩個正交平面(XZ與YZ平面)。具有方向性的兩點相關(guān)函數(shù)可寫為，其中α={p, p-m}，β={XY, XZ, YZ}。

1.3 線性路徑函數(shù)

線性路徑函數(shù)L(i)(r)表示沿向量r方向任意截取長度為的線段完全落入第i相的概率。線性路徑函數(shù)是關(guān)于r的單調(diào)遞減函數(shù)，當(dāng)r=0 時，線性路徑函數(shù)退化為單點相關(guān)函數(shù)。由于在巖心中想找到一條長度非常大的線段完全落入孔隙或基質(zhì)相的概率微乎其微，對于較大的r值，值迅速減小并趨于零即L(i)(∞ )= 0。線性路徑函數(shù)表征系統(tǒng)中各相沿線性路徑的連通特性，是描述多孔介質(zhì)微觀結(jié)構(gòu)的重要函數(shù)。

1.4 兩點簇函數(shù)

如果將系統(tǒng)中第i相中“簇”定義為其中任意兩點都可通過連續(xù)路徑相通的區(qū)域，則兩點簇函數(shù)表示系統(tǒng)中隨機選擇的兩點落入第i相中同一簇的概率[22]。對于統(tǒng)計意義上的均勻介質(zhì)，兩點簇函數(shù)只與兩點間距離有關(guān)，即與線性路徑函數(shù)相比，C2包含某一相中簇的完整聯(lián)通信息，對于材料的系統(tǒng)屬性具有重大影響。在進(jìn)行材料的三維重構(gòu)時，兩點相關(guān)函數(shù)和線性路徑函數(shù)一般均由二維切片的圖像獲取，但是兩點簇函數(shù)需要通過完整的三維圖像獲得，否則很難反映真實的三維系統(tǒng)連通信息。

目前關(guān)于兩點簇函數(shù)的研究已證明，如果將系統(tǒng)中隨機選擇的兩點落入第i相中不同簇的概率記為，則存在以下關(guān)系：

由上式可知，兩點簇函數(shù)實際上是兩點相關(guān)函數(shù)的連通性表征。對于本文中討論的頁巖有機質(zhì)孔隙重構(gòu)的情況，兩點簇函數(shù)寫作，其中α={p, p-m}，β={XY, XZ, YZ}。對于數(shù)字巖心的重構(gòu)來所，孔隙相的C2(r)為短距離函數(shù)，其值隨兩點距離增加迅速下降，當(dāng)r接近連通孔隙中最遠(yuǎn)距離時C2(r)趨于零。兩點簇函數(shù)雖然也屬于兩點概率函數(shù)，但是卻含有高階結(jié)構(gòu)信息，因此兩點簇函數(shù)較之兩點相關(guān)函數(shù)是對結(jié)構(gòu)變化敏感性更強的統(tǒng)計指示函數(shù)。

2 模擬退火法重構(gòu)頁巖數(shù)字巖心

2.1 模擬退火法數(shù)值重構(gòu)原理

本次重構(gòu)頁巖巖心有機質(zhì)孔隙選擇使用模擬退火法，與高斯隨機場法和最近發(fā)展起來的光柵路徑法相比，模擬退火法可以結(jié)合任何類型的關(guān)聯(lián)函數(shù)，靈活性較大。模擬退火法是局部搜索算法的擴展，理論上是一種全局最優(yōu)算法，在本文中簡稱為SA(Simulated Annealing method)。它通過模擬退火過程中原子能量的概率分布進(jìn)行優(yōu)化計算，定義第k+1 次搜索時狀態(tài)的接受概率P為：

其中rand(0,1)為在[0.0,1.0]范圍內(nèi)選取具有等概率的隨機數(shù)。滿足該準(zhǔn)則才是系統(tǒng)第k+1 次可更新的狀態(tài)。當(dāng)ΔE小于預(yù)先設(shè)定的容許誤差時迭代終止，本文中取ΔE= 10-10。

模擬退火法數(shù)值重構(gòu)巖心的主要流程為：基于直接成像實驗所得圖像計算約束函數(shù)值；隨機生成初始模型，設(shè)定初始溫度﹑降溫速率﹑終止溫度和換點失敗率，計算初始模型的相應(yīng)函數(shù)值及能量值；換點并計算新模型的函數(shù)值及能量值，判斷是否滿足更新條件，若滿足則替代原模型，若不滿足判斷是否滿足降溫條件，即同一溫度下模型更新失敗率大于某一臨界值，若滿足則進(jìn)行降溫，若不滿足則繼續(xù)換點更新模型；當(dāng)系統(tǒng)溫度降低到終止溫度，表明此時模型已經(jīng)趨于最優(yōu)狀態(tài)，即可結(jié)束程序。

2.2 樣品信息及成像處理

本次重構(gòu)實驗所用樣品為四川盆地下志留系龍馬溪組頁巖巖心。首先對頁巖樣品進(jìn)行研磨拋光后在表面鍍金，然后分別進(jìn)行場發(fā)射掃描電子顯微鏡(SEM)和聚焦離子束掃描電鏡實驗(FIB-SEM)。場發(fā)射掃描電鏡成像實驗在中國石油大學(xué)(北京)能源材料微結(jié)構(gòu)實驗室進(jìn)行，所使用的設(shè)備是荷蘭FEI公司的Quanta 200F場發(fā)射環(huán)境掃描電子顯微鏡；聚焦離子束掃描電鏡實驗在中國科學(xué)院地質(zhì)與地球物理研究所微納結(jié)構(gòu)成像與數(shù)字巖石物理實驗室進(jìn)行，采用ZEISS Crossbeam 540 型號聚離子束掃描電鏡，該設(shè)備聚焦離子束最優(yōu)分辨率為3 nm/像素，即切片最薄厚度為3 nm；場發(fā)射掃描電鏡最優(yōu)分辨率為0.9 nm/像素，放大倍數(shù)12～2×106倍。

通過觀察大量SEM和FIB-SEM實驗圖像后發(fā)現(xiàn)，有機質(zhì)發(fā)育具有各向異性。大量有機質(zhì)在平行層理面上無定向發(fā)育，而在垂直層理面則多為連續(xù)長條形，如圖2 所示。

基于上述觀察，從頁巖樣品的FIB-SEM實驗結(jié)果中選取400 張連續(xù)的富有機質(zhì)圖像，每張圖像的像幅為1900×3000，分辨率為6.5 nm/像素，實際物理尺寸為12.35 μm×19.5 μm(圖3a)，其中抽提的有機質(zhì)部分三維圖像見圖3b。對圖像進(jìn)行初步分析后使用ImageJ軟件對每張圖像截取400×400 像素大小區(qū)域，然后對圖像進(jìn)行非局部均值濾波和閾值分割并將像幅尺寸粗化至100×100×100，有機質(zhì)占比為10.44%，使用AVIZO軟件進(jìn)行三維可視化，結(jié)果如圖3b所示。從圖中可以看出，有機質(zhì)孔隙發(fā)育具有明顯定向性，與二維圖像觀察結(jié)果相吻合。以此三維數(shù)字巖心作為參考模型，進(jìn)行頁巖有機質(zhì)三維重構(gòu)，并將參考模型用于評價重構(gòu)三維數(shù)字巖心的準(zhǔn)確性，接下來對具體重構(gòu)過程進(jìn)行詳細(xì)說明。

2.3 有機質(zhì)與有機質(zhì)孔隙重構(gòu)

模擬退火法在初始化時會在模型中隨機分布灰度值為0 的像素點，為了在不影響模型效果的情況下加快重構(gòu)速度，提出將初始模型優(yōu)化方法，即將孔隙簡化為規(guī)則形狀取代模型最初的隨機分布，然后結(jié)合兩點簇函數(shù)和預(yù)判方法繼續(xù)進(jìn)行重構(gòu)。下面以頁巖樣品的FIB-SEM實驗圖像為基礎(chǔ)，通過對有機質(zhì)和有機質(zhì)孔隙的重構(gòu)過程對初始模型方法進(jìn)行詳細(xì)說明。

2.3.1 有機質(zhì)重構(gòu)

有機質(zhì)的形狀﹑輪廓復(fù)雜，并且分布散亂，不具有明顯的分布特征規(guī)律。采用面積﹑周長﹑角度﹑延伸率等參數(shù)對其幾何分布特征進(jìn)行描述[24]。在計算特征參數(shù)時，有機質(zhì)面積和周長均為實際有機質(zhì)的面積和周長，同時假設(shè)任何一個有機質(zhì)形狀都可以用一個與其面積相等的橢圓來進(jìn)行等效，從而進(jìn)行有機質(zhì)延伸率的計算。有機質(zhì)延伸率的定義為，該有機質(zhì)團塊與具有相同標(biāo)準(zhǔn)二階中心矩的橢圓的長軸長度與短軸長度之比(圖4)。

圖3 頁巖樣品有機質(zhì)FIB-SEM圖像與重構(gòu)結(jié)果Fig. 3 FIB-SEM images and reconstruction results of the organic matters in the shale samples

第一步，選取二維圖像。選取參考模型在XY﹑XZ和YZ平面上的圖像，如圖5 所示。

第二步，參數(shù)統(tǒng)計。對三維參考模型中的有機質(zhì)進(jìn)行參數(shù)統(tǒng)計，用于驗證重構(gòu)模型的準(zhǔn)確性。對所選取二維圖像中的有機質(zhì)的質(zhì)心坐標(biāo)﹑面積﹑角度和延伸率等參數(shù)進(jìn)行統(tǒng)計分析，對面積和角度分組后，計算XY﹑XZ和YZ平面的面積分布頻率和角度分布頻率。將有機質(zhì)簡化為長方體。

第三步，根據(jù)有機質(zhì)的面積和角度統(tǒng)計結(jié)果，隨機給各有機質(zhì)隨機分配XY﹑XZ和YZ平面的面積值和角度值，并根據(jù)所分配的參數(shù)值計算各有機質(zhì)長方體的三條邊長和體積。逐個隨機分配有機質(zhì)參數(shù)，直至已分配的總有機質(zhì)體積占比與二維圖像中有機質(zhì)平均占比的相對誤差小于0.5%。

圖4 頁巖有機質(zhì)延伸率示意圖Fig. 4 Schematic diagram of the organic matter elongation ratio in the shale sample

第四步，計算已生成有機質(zhì)的面積和角度分布頻率，并與原始二維圖像統(tǒng)計結(jié)果中的面積和角度分布頻率相比較，計算平均相對誤差。若平均相對誤差均小于10%，則將已生成的有機質(zhì)各參數(shù)按體積值從大到小依次排序；否則，重復(fù)第三步直至滿足條件。

第五步，在模型中給有機質(zhì)隨機分配質(zhì)心點坐標(biāo)﹑最大體積值及對應(yīng)的各參數(shù)值，計算該有機質(zhì)長方體的8 個頂點坐標(biāo)和6 個平面方程。判斷各頂點是否在模型范圍內(nèi)，若在范圍內(nèi)，則計算模型中有哪些像素點在這6 個平面所包圍的區(qū)域內(nèi)，確定該有機質(zhì)長方體的完整區(qū)域；若不在范圍內(nèi)，則重復(fù)該步驟直至各頂點均在模型范圍內(nèi)。

第六步，通過計算三維模型中有機質(zhì)相的體素點數(shù)量判斷新生成的有機質(zhì)與之前已生成的有機質(zhì)是否重疊。若重疊，則重復(fù)執(zhí)行第五步；若不重疊，則刪除 該有機質(zhì)對應(yīng)的各參數(shù)值，執(zhí)行下一步。

第七步，按照體積值從大到小依次給有機質(zhì)分配質(zhì)心點坐標(biāo)和各參數(shù)值，重復(fù)執(zhí)行第五步和第六步，直至生成所有有機質(zhì)長方體，有機質(zhì)初始模型優(yōu)化完成。

第八步，初始模型優(yōu)化完成后，應(yīng)用模擬退火法結(jié)合兩點簇函數(shù)和預(yù)判方法進(jìn)行重構(gòu)。重構(gòu)結(jié)果三維可視化如圖3c所示。

2.3.2 有機質(zhì)孔隙重構(gòu)

選擇重構(gòu)有機質(zhì)所截取的400 張400×400 像素大小的原始圖片，對圖像進(jìn)行非局部均值濾波和閾值分割提取其中的有機質(zhì)孔隙，并將像幅尺寸粗化至100×100×100，有機質(zhì)孔隙占比為0.52%，使用AVIZO軟件進(jìn)行三維可視化(圖6a)，以此三維數(shù)字巖心作為有機質(zhì)孔隙參考模型，進(jìn)行頁巖有機質(zhì)孔隙三維重構(gòu)，并將參考模型用于評價有機質(zhì)孔隙重構(gòu)結(jié)果的準(zhǔn)確性。接下來將對具體重構(gòu)過程進(jìn)行詳細(xì)說明，流程框圖如圖7 所示。

圖5 頁巖樣品FIB-SEM實驗獲得三個正交方向的有機質(zhì)二值化圖像示例Fig. 5 Example of the binary images of the organic matters in three orthogonal directions obtained by FIB-SEM experiments of the shale samples

圖6 頁巖樣品有機質(zhì)孔隙的FIB-SEM圖像和重構(gòu)結(jié)果Fig.6 FIB-SEM image and reconstruction result of the organic matter pores in the shale samples

第一步，選取二維圖像。從有機質(zhì)孔隙FIB-SEM三維參考模型(圖6a)的XY﹑XZ和YZ平面上選取圖像，得到二值化的孔隙二維圖像。

第二步，參考模型的參數(shù)統(tǒng)計。巖心結(jié)構(gòu)主要由孔隙的大小和形狀及分布決定的，因此可從孔隙半徑﹑孔隙半徑分布情況等方面來表征巖心結(jié)構(gòu)。對三維參考模型中的有機質(zhì)孔隙進(jìn)行參數(shù)統(tǒng)計，統(tǒng)計結(jié)果用于驗證重構(gòu)模型的準(zhǔn)確性。對所選取二維圖像中的有機質(zhì)孔隙面積進(jìn)行統(tǒng)計，計算對應(yīng)的孔隙半徑，對孔隙半徑分組后計算孔隙半徑分布頻率。

第三步，初始化重構(gòu)模型孔隙半徑。將重構(gòu)三維數(shù)字巖心中所有有機質(zhì)像素點的灰度值暫時設(shè)置為0.5，將有機質(zhì)孔隙簡化為球體。給有機質(zhì)孔隙隨機分配質(zhì)心點坐標(biāo)和孔隙半徑，直至已分配的有機質(zhì)總孔隙度與二維圖像中有機質(zhì)孔隙平均占比的相對誤差小于0.5%。

第四步，計算已生成有機質(zhì)孔隙的孔隙半徑分布頻率，并與原始二維圖像有機質(zhì)孔隙統(tǒng)計結(jié)果中的孔隙半徑分布頻率相比較，計算平均相對誤差。若平均相對誤差均小于10%，則將已生成的有機質(zhì)孔隙體積值從大到小依次排序；否則，重復(fù)第三步。

第五步，初始化有機質(zhì)孔隙位置。在模型中隨機給出有機質(zhì)孔隙的質(zhì)心點坐標(biāo)，將上一步得到的最大體積值賦給該有機質(zhì)孔隙并判斷其是否完全在重構(gòu)模型的有機質(zhì)區(qū)域內(nèi)部。若在，則確定該有機質(zhì)孔隙的完整區(qū)域；若不在，則重復(fù)此步驟直至完全在有機質(zhì)區(qū)域內(nèi)部。

第六步，通過計算三維模型中有機質(zhì)孔隙相的體素點數(shù)量判斷新生成的有機質(zhì)孔隙與之前已生成的有機質(zhì)孔隙是否重疊。若重疊，則重復(fù)執(zhí)行第五步；若不重疊，則刪除最大體積值，執(zhí)行下一步。

第七步，重復(fù)執(zhí)行第五步和第六步，直至生成所有有機質(zhì)孔隙，此時有機質(zhì)孔隙初始模型優(yōu)化完成。

第八步，初始模型優(yōu)化完成后，利用模擬退火法結(jié)合兩點簇函數(shù)和預(yù)判算法采用進(jìn)行重構(gòu)，得到能準(zhǔn)確反映參考系統(tǒng)的有機質(zhì)孔隙模型。重構(gòu)結(jié)果三維可視化如圖6b所示。

2.4 頁巖樣品三維重構(gòu)結(jié)果

頁巖樣品基質(zhì)孔隙可采取與有機質(zhì)孔隙重構(gòu)相同方法進(jìn)行，獲得基質(zhì)孔隙分布模型后采用“雙重區(qū)域”方法進(jìn)行重構(gòu)[9,12]。該方法將頁巖樣品分為2 大區(qū)域：①孔隙發(fā)育程度較好的孔隙發(fā)育區(qū)，主要由有機質(zhì)構(gòu)成；②孔隙發(fā)育程度較低的基質(zhì)區(qū)，主要由粘土礦物和無機質(zhì)礦物及其中的孔隙構(gòu)成。當(dāng)孔隙發(fā)育區(qū)與基質(zhì)區(qū)的內(nèi)部孔隙均重構(gòu)完成后，再根據(jù)各區(qū)域在樣品中的原始位置進(jìn)行耦合疊加，最終形成整體的重構(gòu)巖心模型。該方法可以提高頁巖數(shù)字巖心重構(gòu)的精度并大大降低重構(gòu)時間。

圖7 頁巖樣品有機質(zhì)孔隙重構(gòu)流程框圖Fig.7 Flow chart of the organic matter pores reconstruction in the shale samples

與之前獲得的有機質(zhì)重構(gòu)模型﹑有機質(zhì)孔隙模型相耦合從而建立最終的頁巖樣品數(shù)字巖心模型(圖8)。真實三維數(shù)字巖心的有機質(zhì)﹑有機質(zhì)孔隙和基質(zhì)孔隙體積占比分別為：10.44%﹑0.43%和 0.042%。重構(gòu)三維數(shù)字巖心的有機質(zhì)﹑有機質(zhì)孔隙和基質(zhì)孔隙體積占比分別為：10.40%﹑0.43%和0.042%，與真實三維數(shù)字巖心基本吻合。

3 結(jié)果分析

本節(jié)主要從孔隙度及孔隙半徑﹑面積﹑延伸率和兩點簇函數(shù)等方面對比重構(gòu)數(shù)字巖心與真實頁巖巖心，評價分析重構(gòu)三維數(shù)字巖心的準(zhǔn)確性。

(1)等體積半徑

統(tǒng)計真實三維數(shù)字巖心和重構(gòu)三維數(shù)字巖心中各有機質(zhì)﹑有機質(zhì)孔隙和基質(zhì)孔隙的體積，并計算等體積半徑得到半徑分布曲線，如圖9 所示。通過對比真實巖心和重構(gòu)巖心的有機質(zhì)等體積半徑分布曲線﹑有機質(zhì)孔隙半徑分布曲線和基質(zhì)孔隙半徑分布曲線，發(fā)現(xiàn)吻合程度良好，表明重構(gòu)所得三維數(shù)字巖心表征頁巖有機質(zhì)﹑有機質(zhì)孔隙和基質(zhì)孔隙尺度較準(zhǔn)確。

圖8 重構(gòu)頁巖樣品三維數(shù)字巖心，其中深灰色代表有機質(zhì)，灰色代表有機質(zhì)孔隙，白色代表基質(zhì)孔隙，黑色代表基質(zhì)Fig. 8 The three-dimensional reconstructed digital core of the shale sample, where dark gray stands for organic matter,gray for organic matter pores, white for matrix pores, and black for matrix

(2)面積

統(tǒng)計真實三維數(shù)字巖心和重構(gòu)三維數(shù)字巖心XY﹑XZ和YZ平面上所有有機質(zhì)﹑有機質(zhì)孔隙和基質(zhì)孔隙的面積，對統(tǒng)計結(jié)果分組后得到面積分布曲線，如圖10～12 所示。各面積分布曲線均比較吻合，表明三維重構(gòu)效果較好。

(3)有機質(zhì)延伸率

不管是真實巖心還是重構(gòu)巖心，相比較于平行層理面，兩個垂直層理面上的有機質(zhì)延伸率明顯偏高，如圖13 所示，表明頁巖有機質(zhì)分布具有各向異性。

(4)兩點簇函數(shù)

圖9 真實頁巖樣品與重構(gòu)頁巖數(shù)字巖心中有機質(zhì)、有機質(zhì)孔隙與基質(zhì)孔隙等體積半徑分布比較Fig. 9 Comparison of the equivalent volume radius distribution of organic matter, organic matter pores and matrix pores between the real shale sample and the reconstructed shale digital core

圖10 真實頁巖樣品與重構(gòu)頁巖數(shù)字巖心中有機質(zhì)面積分布比較Fig. 10 Comparison of the area distribution of organic matter between the real shale sample and the reconstructed shale digital core

圖11 真實頁巖樣品與重構(gòu)頁巖數(shù)字巖心中有機質(zhì)孔隙面積分布比較Fig. 11 Comparison of the area distribution of organic matter pores between the real shale sample and the reconstructed shale digital core

圖12 真實頁巖樣品與重構(gòu)頁巖數(shù)字巖心中基質(zhì)孔隙面積分布比較Fig. 12 Comparison of the area distribution of matrix pores between the real shale sample and the reconstructed shale digital core

圖13 真實頁巖樣品與重構(gòu)頁巖數(shù)字巖心中有機質(zhì)延伸率比較Fig. 13 Comparison of the elongation ratio of organic matter between the real shale sample and the reconstructed shale digital core

為了更準(zhǔn)確地重構(gòu)頁巖三維數(shù)字巖心，使用約束函數(shù)時結(jié)合了可表征拓?fù)溥B通性的兩點簇函數(shù)，為研究頁巖各向異性，沿3 個正交方向計算獨立的結(jié)構(gòu)信息，有機質(zhì)部分的兩點簇函數(shù)結(jié)果如圖14 所示。從圖中可以看出，重構(gòu)數(shù)字巖心與真實數(shù)字巖心的在各方向上均擬合較好。有機質(zhì)兩點簇函數(shù)k方向與i﹑j方向的值差距較大，表明有機質(zhì)具有明顯的各向異性。

4 結(jié)束語

圖14 真實頁巖樣品與重構(gòu)頁巖數(shù)字巖心中有機質(zhì)部分兩點簇函數(shù)在3 個正交方向上的分布比較Fig. 14 Comparison of the two-point cluster functions of the organic matters in three orthogonal directions between the real shale sample and the reconstructed shale digital core

針對頁巖樣品中有機質(zhì)分布的各向異性特征，本文提出在傳統(tǒng)模擬退火數(shù)值重構(gòu)方法基礎(chǔ)上引進(jìn)兩點簇函數(shù)作為表征有機質(zhì)分布連續(xù)性的統(tǒng)計指示函數(shù)，并利用實際頁巖巖心樣品的數(shù)值重構(gòu)對該方法的應(yīng)用進(jìn)行說明，通過對頁巖樣品重構(gòu)結(jié)果的有效性評價，證明了兩點簇函數(shù)方法表征頁巖有機質(zhì)各向異性分布特征的有效性。

致謝

感謝美國亞利桑那州立大學(xué)Yang Jiao教授在兩點簇函數(shù)程序編制方面提供的幫助。