模擬退火法結合兩點簇函數重構與表征頁巖有機質各向異性

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

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

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

0 引言

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

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

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

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

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

1 重構算法統計函數定義

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

1.1 單點相關函數

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

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

1.2 兩點相關函數

在上述多相系統中第i相和第j相的兩點相關函數定義為：

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

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

兩點相關函數只與兩點之間的相對位移矢量有關，即

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

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

1.3 線性路徑函數

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

1.4 兩點簇函數

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

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

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

2 模擬退火法重構頁巖數字巖心

2.1 模擬退火法數值重構原理

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

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

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

2.2 樣品信息及成像處理

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

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

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

2.3 有機質與有機質孔隙重構

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

2.3.1 有機質重構

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

圖3 頁巖樣品有機質FIB-SEM圖像與重構結果Fig. 3 FIB-SEM images and reconstruction results of the organic matters in the shale samples

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

第二步，參數統計。對三維參考模型中的有機質進行參數統計，用于驗證重構模型的準確性。對所選取二維圖像中的有機質的質心坐標﹑面積﹑角度和延伸率等參數進行統計分析，對面積和角度分組后，計算XY﹑XZ和YZ平面的面積分布頻率和角度分布頻率。將有機質簡化為長方體。

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

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

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

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

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

第七步，按照體積值從大到小依次給有機質分配質心點坐標和各參數值，重復執行第五步和第六步，直至生成所有有機質長方體，有機質初始模型優化完成。

第八步，初始模型優化完成后，應用模擬退火法結合兩點簇函數和預判方法進行重構。重構結果三維可視化如圖3c所示。

2.3.2 有機質孔隙重構

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

圖5 頁巖樣品FIB-SEM實驗獲得三個正交方向的有機質二值化圖像示例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 頁巖樣品有機質孔隙的FIB-SEM圖像和重構結果Fig.6 FIB-SEM image and reconstruction result of the organic matter pores in the shale samples

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

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

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

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

第五步，初始化有機質孔隙位置。在模型中隨機給出有機質孔隙的質心點坐標，將上一步得到的最大體積值賦給該有機質孔隙并判斷其是否完全在重構模型的有機質區域內部。若在，則確定該有機質孔隙的完整區域；若不在，則重復此步驟直至完全在有機質區域內部。

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

第七步，重復執行第五步和第六步，直至生成所有有機質孔隙，此時有機質孔隙初始模型優化完成。

第八步，初始模型優化完成后，利用模擬退火法結合兩點簇函數和預判算法采用進行重構，得到能準確反映參考系統的有機質孔隙模型。重構結果三維可視化如圖6b所示。

2.4 頁巖樣品三維重構結果

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

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

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

3 結果分析

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

(1)等體積半徑

統計真實三維數字巖心和重構三維數字巖心中各有機質﹑有機質孔隙和基質孔隙的體積，并計算等體積半徑得到半徑分布曲線，如圖9 所示。通過對比真實巖心和重構巖心的有機質等體積半徑分布曲線﹑有機質孔隙半徑分布曲線和基質孔隙半徑分布曲線，發現吻合程度良好，表明重構所得三維數字巖心表征頁巖有機質﹑有機質孔隙和基質孔隙尺度較準確。

圖8 重構頁巖樣品三維數字巖心，其中深灰色代表有機質，灰色代表有機質孔隙，白色代表基質孔隙，黑色代表基質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)面積

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

(3)有機質延伸率

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

(4)兩點簇函數

圖9 真實頁巖樣品與重構頁巖數字巖心中有機質、有機質孔隙與基質孔隙等體積半徑分布比較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 真實頁巖樣品與重構頁巖數字巖心中有機質面積分布比較Fig. 10 Comparison of the area distribution of organic matter between the real shale sample and the reconstructed shale digital core

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

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

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

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

4 結束語

圖14 真實頁巖樣品與重構頁巖數字巖心中有機質部分兩點簇函數在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

針對頁巖樣品中有機質分布的各向異性特征，本文提出在傳統模擬退火數值重構方法基礎上引進兩點簇函數作為表征有機質分布連續性的統計指示函數，并利用實際頁巖巖心樣品的數值重構對該方法的應用進行說明，通過對頁巖樣品重構結果的有效性評價，證明了兩點簇函數方法表征頁巖有機質各向異性分布特征的有效性。

致謝

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