含天然氣水合物沉積物三維孔隙結構數值重建

曾建邦　胡高偉　李隆鍵　陳　強　吳能友　王廣君

1. 載運工具與裝備教育部重點實驗室·華東交通大學 2. 青島海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室3. 國土資源部天然氣水合物重點實驗室·青島海洋地質研究所 4. 低品位能源利用技術及系統教育部重點實驗室·重慶大學

含天然氣水合物沉積物三維孔隙結構數值重建

曾建邦1,2,3,4胡高偉2,3李隆鍵4陳強2,3吳能友2,3王廣君1

1. 載運工具與裝備教育部重點實驗室·華東交通大學 2. 青島海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室3. 國土資源部天然氣水合物重點實驗室·青島海洋地質研究所 4. 低品位能源利用技術及系統教育部重點實驗室·重慶大學

曾建邦等.含天然氣水合物沉積物三維孔隙結構數值重建. 天然氣工業，2016,36(5):128-133.

三維孔隙結構重建是揭示沉積物內天然氣水合物（以下簡稱水合物）形成機理及分布規律的基礎和前提。在現有的數值重建方法中，模擬退火法被廣泛應用。為使該方法的重建結構更加接近于實際，對該方法進行了改進，所構建的數值重建模型能夠納入沉積物顆粒形狀（包括圓球狀、橢球狀、扁球狀、長扁球狀、片狀和針狀）、尺寸分布（如正態分布）以及水合物分布模式（如膠結狀、懸浮狀或顆粒狀）等重要結構與統計信息，并利用改進后的模型對含水合物沉積物的三維孔隙結構進行了數值重建。最后，對重建的沉積物（不含水合物）孔隙結構進行特征化分析，可獲得沉積物內孔隙率和孔徑分布等信息；對重建的含水合物沉積物孔隙結構進行特征化分析，可獲得各組分體積分布及其截面平均體積分數沿x軸方向的分布曲線等信息。結論認為：所開發的數值重建模型能夠重建出與實際情況較為接近的含水合物沉積物三維孔隙結構，且能揭示水合物分布模式對沉積物孔隙內水（或冰）和天然氣空間分布的影響。

含天然氣水合物 沉積物 孔隙結構 分布模式 模擬退火法 數值重建 特征化分析

沉積物孔隙結構重建是天然氣水合物介觀孔（或顆粒）尺度模型的基礎和前提，其重建方法有實驗和數值兩種：實驗方法是借助X射線衍射、核磁共振或X-CT等高分辨率儀器對沉積物進行二維或三維成像[1-3]，該方法受實驗設備價格、樣本制備要求、實驗操作人員技能與儀器分辨率等條件的限制。因此，借助計算機對沉積物孔隙結構進行重建是一種行之有效的方法。在現有的數值重建方法中，模擬退火法除了可以考慮孔隙率、組元體積分數以及兩點相關函數等基本信息之外，還可兼顧更多與沉積物孔隙結構相關的重要信息，因而其重建結構與實際更為接近[4]。目前，該方法已在砂巖[5]、鋰離子電池正[6]、負極[7]等復雜孔隙介質的三維重建中得到廣泛地應用。

海底沉積物主要有三大來源[8]：①風化物從內陸巖體被剝蝕后經河水攜帶至出海口，搬運到海中，慢慢沉淀下來的接近球狀的沉積物；②火山灰從火山口噴出后經風攜帶至海上，最后沉積到海底的碎屑巖大致有圓球、橢球、扁球和長扁球狀等4種形狀；③海洋生物死后沉積到海底，被逐漸溶解、分解后剩下的片狀和針狀沉積物。水合物在沉積物內的分布模式大多是根據聲學和熱學探測等手段推測而得[3]，大致有以膠結物的形式存在于沉積物顆粒之間（膠結狀）[2]、游離于沉積物孔隙中（懸浮狀）[9]和與沉積物顆粒一樣成為骨架組成部分（顆粒狀）[1，10]等3種微觀分布狀態。為此，筆者將對模擬退火法實施改進，使其可以考慮沉積物顆粒形狀和尺寸分布以及水合物分布模式等重要信息，并用于重建含水合物沉積物的三維孔隙空間，再對重建孔隙結構進行特征化分析，探討水合物分布模式對沉積物孔隙內水（或冰）和天然氣空間分布的影響。

1　含水合物沉積物孔隙結構數值描述

以Jin等[2]制備的甲烷水合物沉積物為例（圖1），由于受儀器分辨率的限制，圖中僅能分辨出沉積物顆粒、甲烷及水合物和水（或冰）的混合物（文中視為水合物）。在重建過程中，可用相函數I進行描述：

其中r（x, y, z）表示三維空間中任意位置矢量，相i = 0、1、2、3分別代表甲烷、沉積物、水合物和水（或冰）。通過統計平均可得到該結構的基本統計信息，如相i的體積分數εi= ＜Ii（r）＞；兩點相關函數fij（r1, r2） = ＜Ii（r）Ij（r）＞，表示在重構區域內距離為u的兩位置點r1、r2分別屬于i相和j相的概率，代表孔隙結構中各相間的空間分布關系。若假定該孔隙結構各向同性，則相關函數簡化為：其中

基于圖像技術對圖1-a紅色方框內的圖片進行數字化處理得到圖1-b，通過統計平均可得到兩點相關函數（圖2）以及甲烷、沉積物顆粒、水合物和水（或冰）體積分數，分別為0.058 7、0.622 3、0.319 0和0。以10.8 μm作為基本節點尺寸，重建空間區域為Nx= Ny= Nz= 200（即重建區域物理尺寸為Lx= Ly= Lz= 2.16 mm）。重建時，假定重構區域各相同性；沉積物顆粒形狀包括圓球、橢球、扁球、長扁球、片和針狀[8]。這些顆粒形狀均可通過調節橢球3個軸長而得：①當3個軸長相等時為圓球狀；②當3個軸長相近時為橢球狀；③當2個軸長較大，而另一軸長較小/特小時為扁球/片狀；④當1個軸長較長，而另外2個軸長較小/特小時為長扁球/針狀。要對空間橢球進行描述，除了上述3個軸長（各軸長的取值范圍為[10.8 μm，691.2 μm]與圖1中的沉積物顆粒尺寸保持一致，并假定均滿足正態分布）外，還需要考慮3個軸的轉向角取值范圍均為[0°，360°]以及橢球中心坐標，故重建時共需考慮9個參數的變化。

圖1　沉積物內膠結狀水合物分布圖

圖2　含甲烷水合物沉積物內各相兩點相關函數圖

2　模擬退火重建算法及步驟

重建含水合物沉積物三維孔隙空間涉及沉積物骨架、水合物和水（或冰）的重建，具體步驟如下：

1）沉積物顆粒系統初始化。將預設尺寸和形狀的顆粒隨機放置在重建區域內，接收概率為：

式中v表示當前放置顆粒與已有顆粒重疊部分體積的占比；εv表示調節重疊程度的常數。

該步驟直到重建區域內沉積物顆粒體積分數達到預設值時結束。

2）沉積物顆粒的移動。隨機選擇一個顆粒，按隨機產生的距離dr （dx, dy, dz） 移動該顆粒，dx、dy和dz均滿足相同的指數分布，以dx為例：

式中a表示一個控制移動距離的常數。

3）判斷是否接收步驟2）的概率為：

式中T表示系統“退火溫度”；ΔE = Et+1-Et，Et和Et+1分別表示移動前后系統的能量值。

Et可表示為：

式中N = Nx× Ny× Nz；U表示兩位置點r1、r2距離u的最大值；分別表示t時刻含水合物沉積物三維孔隙結構和從圖1-b中提取的作為參照的兩點相關函數（圖2）。

4）重復步驟2）、3），當T = λt+1T0（λ取0.01，為冷卻率；T0取1.0，為初始“退火溫度”）達到目標溫度Tend，則終止計算，得到沉積物顆粒骨架。

5）重建水合物和水（或冰）在沉積物孔隙空間的微觀分布，包括如下3種情況：①對于膠結狀水合物（水合物和水/冰的混合物），只需采用文獻[11]中的模擬退火法對水合物在沉積物內的分布進行重建，但需增加的約束條件有：當水合物移至沉積物顆粒表面時，若不滿足式（5），也接受該次運動；判斷計算是否終止時，在沉積物顆粒表面不允許出現甲烷，除非除沉積物顆粒表面外，孔隙內已沒有水合物。②對于懸浮狀水合物，重建時先后采用文獻[11]中的模擬退火法對水（或冰）和水合物在沉積物內的分布進行重建，但需增加的約束條件有：沉積物顆粒和甲烷的體積分數與圖1-b一致，水合物和水（或冰）的體積分數分別取0.214 5和0.104 5[9]；當水（或冰）移至沉積物顆粒表面時，若不滿足式（5），也接受該次移動；判斷計算是否終止時，在沉積物顆粒表面不允許出現甲烷和水合物，除非除沉積物顆粒表面外，孔隙內已沒有水（或冰）。③對于顆粒狀水合物。先采用步驟1～4對水合物顆粒進行數值重建，初始預設水合物顆粒形狀為球，球半徑的取值范圍為[10.8 μm，162.0 μm]，與文獻[10]保持一致，并假定其滿足正態分布，水合物體積分數仍然取0.214 5；然后采用文獻[11]中構建的模擬退火法對水（或冰）在孔隙內的分布進行數值重建，其體積分數預設值為0.104 5。

3　重建結果分析

3.1重建結果

圖3　數值重建的含膠結狀水合物沉積物圖

基于上節算法分別對含膠結、懸浮和顆粒狀水合物沉積物進行重建。重建過程中先對沉積物顆粒執行退火過程，獲得其在重建區域的分布圖，見圖3-a（結構Ⅰ），圖4-a（結構Ⅱ）和圖5-a（結構Ⅲ）。可見，沉積物骨架均由具有不同大小且形狀各異（圓球、橢球、扁球、長扁球、片和針狀）的沉積物顆粒組成，但三者間的差別較大，說明模型重建結果具有隨機性。然后，對水合物和水（或冰）執行退火過程（其先后順序見前述文字），獲得它們在重建區域的分布（圖3-b、圖4-b和圖5-b），發現相比于含顆粒狀水合物沉積物，在含膠結狀和懸浮狀水合物沉積物內天然氣更為集中；在含懸浮狀水合物沉積物內水（或冰）分布在沉積物顆粒表面，而在含顆粒狀水合物沉積物內其分布規律并不明顯。最后，得到含水合物沉積物孔隙結構，（圖3-c、圖4-c和圖5-c），可清楚地觀察出水合物在沉積物孔隙空間內的分布規律。

為進一步觀察出沉積物顆粒、水合物、天然氣和水（或冰）在含膠結狀、懸浮狀和顆粒狀水合物沉積物內的分布形態，可分別對其進行二維切片（圖3-d、圖4-d和圖5-d），均可發現沉積物顆粒、水合物、天然氣和水（或冰）之間的分界面十分清晰，且分別與實驗拍攝的含膠結狀[2]、懸浮狀[9]和顆粒狀[10]水合物沉積物孔隙結構較為相似。可見，基于筆者所開發的數值重建模型能夠重建出與實際較為接近的含水合物沉積物三維孔隙結構。

圖4　數值重建的含懸浮狀水合物沉積物圖

圖5　數值重建的含顆粒狀水合物沉積物圖

3.2特征化分析

分別對如圖3-a、圖4-a和圖5-a所示的沉積物孔隙結構進行特征化分析：可得到其孔隙率（表1）。

表1　水合物賦存狀態對重建結果的影響表

從表1中可看出，它們與初始設定值均相差不大（最大相對誤差不超過0.32%）。對于沉積物孔徑分析，筆者以歐氏距離轉換法[12]分別計算結構Ⅰ、Ⅱ和Ⅲ內最大內切球直徑[13]，統計分析可得到具有不同尺寸孔的相對體積（具有相同半徑孔的體積之和與總的孔體積的比值）隨孔尺寸大小的分布關系（圖6）。從圖6中可看出：三者的孔徑分布基本一致。另外，筆者還分別計算了它們的平均孔徑（表1）。結合圖6和表1，可以看出不管是最大孔徑還是平均孔徑，三者之間的差別均不大。上述情況與實際相符，即海底沉積物的物理變化過程和現象具有隨機性，但當統計區域達到一定程度時，其內部相關統計信息（如孔隙率、孔徑分布等）卻具有一定的規律性[14]，進一步說明基于本文構建的重建算法能夠較為真實地再現沉積物三維孔隙空間及其內水合物的微觀分布規律。

圖6　沉積物孔隙空間內孔徑分布圖

分別對圖3-c，圖4-c和圖5-c所示的含水合物沉積物微結構進行特征化分析，可得到各組分的體積分數（表1），均與初始預設值相差不大（最大相對誤差不超過1.0%）；還可得到各組分截面平均體積分數沿x軸方向的分布曲線（圖7），發現沉積物、水合物、水（或冰）和天然氣在各自初始額定份額（或平均體積分數）上下波動，其中波動幅度最大的是沉積物顆粒，并按水合物、水（或冰）和天然氣的順序依次減小。然而，波動幅度并不能描述各組分沿x軸方向分布的均勻程度，須引入各組分在截面n上的平均體積分數）與各自在含水合物沉積物內的平均體積分數（）之間的最大相對誤差的絕對值（δ1）和平均相對誤差的絕對值（δ2），即

式中Nx表示重建區域在x軸方向上的節點個數。

基于公式（7），分別對含膠結、懸浮和顆粒狀水合物沉積物中各組分的δ1和δ2進行計算，（表2）。可見，δ1和δ2最小的均為沉積物，且三者之間差別不大；水合物分布規律不同，對其自身、水（或冰）和天然氣的δ1和δ2的影響也不相同，主要體現在：

圖7　各組分體積分數沿x軸方向的分布圖

表2　水合物賦存狀態對δ1和δ2的影響表

①對于含膠結（懸浮或顆粒）狀水合物沉積物，δ1和δ2均按水合物、水（或冰）（膠結狀不予考慮）和天然氣的順序依次增大，即其內部天然氣截面平均體積分數沿x軸方向的變化最為劇烈，并按天然氣、水（或冰）和水合物的順序依次減小；②相比含膠結狀水合物沉積物，在含懸浮狀水合物沉積物中δ1和δ2更大，即其內部各組分截面平均體積分數沿x軸方向的變化更為劇烈；③相比含膠結和懸浮狀水合物沉積物，含顆粒狀水合物沉積物中的水合物截面平均體積分數沿x軸方向的變化最為劇烈，但水（或冰）和天然氣的變化則相對平緩。可見δ1和δ2能較好地衡量水合物分布模式對含水合物沉積物內各組分分布規律的影響。

4　結論

1）基于模擬退火法構建含水合物沉積物孔隙結構數值重建模型，并以實際孔隙率、組元體積分數、沉積物顆粒形狀（圓球、橢球、扁球、長扁球、片狀和針狀）、尺寸（正態分布）分布、水合物分布模式（膠結、懸浮和顆粒狀）和從二維X-CT圖像中提取的兩點相關函數等重要信息作為模型輸入參數，能夠重建出與實際較為相似的含水合物沉積物。

2）根據對重建的沉積物和含水合物沉積物孔隙結構特征化分析，可獲得沉積物內孔隙率、孔徑分布等信息和含水合物沉積物內各組分體積分數及其截面平均體積分數沿x軸方向的分布曲線等信息，可在一定程度上揭示水合物分布模式對沉積物孔隙內水（或冰）和天然氣空間分布的影響規律。

3）由于重建孔隙結構的隨機性及其相關統計信息的規律性，使得所開發的模型能夠重建出與實際更為接近的含水合物沉積物孔隙空間。

[1] Chaouachi M, Falenty A, Sell K, Enzmann F, Kersten M, Haberthur D, et al. Microstructural evolution of gas hydrates in sedimentary matrices observed with synchrotron X-ray computed tomographic microscopy[J]. Geochemistry, Geophysics, Geosystems, 2015, 16(6): 1711-1722.

[2] Jin S, Takeya S, Hayashi J, Nagao J, Kamata Y, Ebinuma T. Structure analyses of artifi cial methane hydrate sediments by microfocus X-ray computed tomography[J]. Japanese Journal of Applied Physics, 2004, 43(8A): 5673-5675.

[3] 胡高偉, 李承峰, 業渝光, 劉昌嶺, 張劍, 刁少波. 沉積物孔隙空間天然氣水合物微觀分布觀測[J]. 地球物理學報, 2014, 57(5): 1675-1682. Hu Gaowei, Li Chengfeng, Ye Yuguang, Liu Changling, Zhang Jian, Diao Shaobo. Observation of gas hydrate distribution in sediment pore space[J]. Chinese Journal of Geophysics, 2014, 57(5): 1675-1682.

[4] Hidajat I, Rastogi A, Singh M, Mohanty KK. Transport properties of porous media from thin-sections[C]//SPE Latin American and Caribbean Petroleum Engineering Conference, 25-28 March 2001, Buenos Aires, Argentina. DOI: http://dx.doi.org/10.2118/69623-MS.

[5] Tang T, Teng Q, He X, Luo D. A pixel selection rule based on the number of different-phase neighbours for the simulated annealing reconstruction of sandstone microstructure[J]. Journal of Microscopy, 2009, 234(3): 262-268.

[6] Zeng JB, Wu W, Jiang FM. Smoothed particle hydrodynamics prediction of effective transport coeffi cients of lithium-ion battery electrodes[J]. Solid State Ionics, 2014, 260(3): 76-85.

[7] 何紹陽, 曾建邦, 蔣方明. 鋰離子電池石墨負極微結構數值重建及特征化分析[J]. 無機材料學報, 2015, 30(9): 906-912. He Shaoyang, Zeng Jianbang, Jiang Fangming. Numerical reconstruction and charactization analysis of microstructure of lithium-ion battery graphite anode[J]. Journal of Inorganic Materials, 2015, 30(9): 906-912.

[8] Lu B, Liu Q, Li GX. Grain and pore factors in acoustic response to seafl oor sediments[J]. Marine Georesources & Geotechnology, 2010, 28(2): 115-129.

[9] Jin S, Nagao J, Takeya S, Jin Y, Hayashi J, Kamata Y, et al. Structural investigation of methane hydrate sediments by microfocus X-ray computed tomography technique under high-pressure conditions[J]. Japanese Journal of Applied Physics, 2006, 45(27): 714-716.

[10] Hyodo M, Yoneda J, Yoshimoto N, Nakata Y. Mechanical and dissociation properties of methane hydrate-bearing sand in deep seabed[J]. Soils and Foundations, 2013, 53(2): 299-314.

[11] Yeong CLY, Torquato S. Reconstructing random media[J]. Physical Review E, Statistical Physics, Plasmas, Fluids & Related Interdisciplinary Topics, 1998, 57(1): 495-506.

[12] Thiele S, Zengerle R, Ziegler C. Nano-morphology of a polymer electrolyte fuel cell catalyst layer-imaging, reconstruction and analysis[J]. Nano Research, 2011, 4(9): 849-860.

[13] Cheolwoong L, Bo Y, Yin LL, Zhu LK. Geometric characteristics of three dimensional reconstructed anode electrodes of lithium-ion batteries[J]. Energies, 2014, 7(4): 2558-2572.

[14] 何幼斌, 王文廣. 沉積巖與沉積相[M]. 北京: 石油工業出版社, 2008. He Youbin, Wang Wenguang. Sedimentary rocks and sedimentary facies[M]. Beijing: Petroleum Industry Press, 2008.

（修改回稿日期 2016-03-08編 輯 羅冬梅）

Numerical reconstruction of 3D pore structure for gas hydrate bearing sediments

Zeng Jianbang1,2,3,4, Hu Gaowei2,3, Li Longjian4, Chen Qiang2,3, Wu Nengyou2,3, Wang Guangjun1
(1. MOE Key Laboratory of Conveyance and Equipment, East China Jiaotong University, Nanchang, Jiangxi 330013, China; 2. Qingdao Institute of Marine Geology, MLR Key Laboratory of Gas Hydrate,Qingdao, Shandong 266071, China; 3. Laboratory of Marine Mineral Resources and Detection Technologies, Qingdao National Laboratory of Marine Science and Technology, Qingdao, Shandong 266071, China; 4. MOE Key Laboratory of Low-grade Energy Utilization Technologies and Systems, Chongqing University, Chongqing 400044, China)

NATUR. GAS IND. VOLUME 36, ISSUE 5, pp.-,5/25/2016. (ISSN 1000-0976; In Chinese)

Reconstruction of 3D pore structures is the basis and prerequisite for revealing the formation mechanisms and distribution laws of gas hydrate in sediments (hereinafter referred to as gas hydrate). Among the existing numerical reconstruction methods, the simulated annealing method is widely used. In this paper, the simulated annealing method was improved so that the actual situations could be presented better by the reconstructed structures. Important structural and statistical information could be introduced into the improved numerical reconstruction model, including the shapes of sediment particles (e.g. spherical, ellipsoidal, oblate spheroid, prolate spheroid, flakey and acicular), size distribution of sediment particles (e.g. normal distribution), and morphology of gas hydrate (e.g. cementing, floating and particle). Then, the 3D pore structures of gas hydrate bearing sediments (hereinafter referred to as sediments) were numerically reconstructed by using the improved model. Conclusions in two aspects were obtained. First, through characterization analysis on the pore structures of sediments, certain key information can be revealed, including the porosity and pore size distribution in sediments, the volume fractions of specific components and the distribution of sectional average volume fraction along x axis. Second, with the proposed numerical reconstruction model, the 3D pore structures of sediments closer to the actual situations can be reconstructed, so that the effect of gas hydrate morphology on the spatial distribution of water (or ice) and gas in the pores of sediments can be revealed.

Gas hydrate; Pore structure of sediments; Hydrate occurrence; Simulated annealing method; Numerical reconstruction; Characterization analysis

10.3787/j.issn.1000-0976.2016.05.019

國家自然科學青年基金項目（編號：51206171）、低品位能源利用技術及系統教育部重點實驗室開放基金項目（編號：LLEUTS-201608）。

曾建邦，1981年生，副研究員，博士；主要從事天然氣水合物成藏機理方面的研究工作。地址：（330013）江西省南昌市經濟技術開發區雙港東大街808號。ORCID:0000-000X-0579-5327。E-mail:jbzeng68@sina.com