基于等效介質理論的天然氣水合物巖石物理建模方法

楊文強 宗兆云*② 姜 曼 劉欣欣

(①中國石油大學(華東)地球科學與技術學院，山東青島266580；②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島266580；③中國地質調查局青島海洋地質研究所，山東青島266071)

1 問題的提出

天然氣水合物又稱可燃冰，是在低溫高壓下，由水的冰晶格架及其吸附的天然氣分子組成的一種籠形化合物，廣泛分布于大陸邊緣海底、深水盆地和永久凍土帶中，儲量豐富，具備成為21世紀新型能源的潛力[1]。目前在天然氣水合物勘探中，地震方法起到了主要作用。如何從地震響應中預測天然氣水合物的空間展布特征，是一個需要不斷探索的問題。

在地震資料上，天然氣水合物通常表現出一些與常規油氣儲層不同的響應特征，如似海底反射層(BSR)、振幅空白帶和極性反轉等[2]，這些獨特的地震反射特征已經被廣泛地應用于天然氣水合物的識別，是海域天然氣水合物存在的重要地球物理標志。但是，這些反射特征與水合物地層的關聯程度還存在爭議，已發現的水合物與這些響應并不具嚴格的對應關系[3]。

巖石物理模型可用于預測含天然氣水合物地層的內部結構和天然氣水合物的含量。將沉積層的彈性波速度與密度、孔隙度、各種微觀礦物組分及水合物、水和氣體的物性參數等聯系起來[4]，是水合物地層的地震響應分析、高分辨率疊前地震反演及水合物飽和度的定量估算的有效工具[5-6]。因此，眾多學者提出了一些水合物巖石物理經驗公式和理論模型。Lee等[7]提出基于雙相和三相介質時間平均方程和的Wood方程，早期被廣泛使用，這類方法簡單有效但缺乏實際物理意義，模型參數的選擇需要大量的測井資料作為先驗信息進行約束。Helgerud等[8]基于等效介質模型將水合物作為基質礦物的一部分構成巖石骨架，臨界孔隙度通常取經驗值，但實際水合物地層臨界孔隙度很難準確確定。Dvokin等[9-10]考慮水合物膠結作用及其對彈性參數的影響，當水合物含量較低時，預測的地層速度會發生急劇增大的現象，計算得到的彈性參數值較實測值高。Zimmerman等[11]使用Kuster-Toks?z模型將水和水合物的混合物作為巖石骨架，黏土作為填充物，該方法適用于低孔隙度的地層。Lee等[12]、Qian等[13]應用Gassmann方程，結合接觸理論，計算了水合物儲層彈性參數，該方法適用于未固結的海洋沉積物，對于固結沉積物的地層應用效果無法令人滿意。劉欣欣等[4]在考慮水合物剪切性質和孔隙微觀結構情況下，建立了孔隙填充和顆粒支撐兩種模式的水合物模型，在實際資料測試中取到了較好的效果。顆粒支撐模式[14-15]和孔隙填充模式[16-18]是物理意義上的兩種極端分布情況，當天然氣水合物在固體顆粒邊緣沉積作為孔隙空間充填物的一部分時，如圖1a所示，為孔隙填充模式；當水合物飽和度較高時，孔隙填充模式可能變為顆粒支撐模式，增加固體骨架的穩定性，成為固體基質的一部分，如圖1b所示。

圖1 水合物微觀模式示意圖(a)孔隙填充模式； (b)顆粒支撐模式 白色表示孔隙，黃色表示水合物，藍色表示固體顆粒

孔隙填充和顆粒支撐兩種模式的巖石物理模型是對水合物微觀賦存情況的極端假設，但實際上水合物同時分布于孔隙和骨架基質中。因此，本文考慮水合物實際的微觀分布模式，基于等效介質理論[19-20]和Gassmann固體替換方程[21]，將實測縱波速度作為約束，得到更加合理的水合物微觀分布情況，構建天然氣水合物地層的巖石物理模型。通過測試數據和實際資料驗證該模型，效果較好。該方法為后續利用地震數據進行天然氣水合物地層物性參數的定量解釋提供了技術支持。

2 方法原理

天然氣水合物具有類固體特性，因此巖石物理建模時需要考慮水合物的剪切性質及微觀分布情況，具體流程如圖2所示。

圖2 本文巖石物理建模流程 Vh為水合物體積含量，p為比例參數，VP為縱波速度

(1)基于等效介質理論，固體基質由石英、方解石、黏土等礦物組分均勻混合構成，向基質中加入孔隙結構可構成干燥巖石骨架。假設初始狀態水合物作為孔隙填充物的一部分，并構成部分孔隙流體形成飽和沉積物，此時，水合物的體積含量可通過孔隙度和含水飽和度計算得到。

(2)引入比例參數p(0≤p≤1)，對水合物體積含量進行分配，將一部分水合物分配到基質中，成為基質礦物的一部分。部分天然氣水合物成為基質礦物后，孔隙度降低，因此認為孔隙度下降的程度即水合物在孔隙中減少的程度。水合物分配到基質中后將影響各基質礦物組分的體積含量，重新計算各個組分的模量，可以求得飽和巖石的體積模量、剪切模量，進而計算得到縱、橫波速度。

(3)受p的影響，各組分的彈性模量處于動態變化中。最后以測井實測的縱波速度約束巖石物理模型的速度，并比較計算值與實測縱波速度的差異；如果這種差異不可接受，則通過迭代更新p，調整天然氣水合物微觀分布，直到模型計算的縱波速度與實測的誤差在允許接受范圍內，進而可以計算得到合理的縱波以及橫波速度。

在建模過程中，關鍵步驟包括固體基質、干燥骨架、孔隙流體等效模量及飽和巖石彈性模量的計算。

2.1 固體基質等效模量計算

水合物地層涵蓋了泥巖、砂巖或礫巖等多種類型，主要由陸源碎屑礦物、黏土礦物和碳酸鹽礦物組成[22]。其中，碎屑礦物以石英、白云母、長石等輕礦物為主；黏土礦物主要為伊利石、蒙脫石等；碳酸鹽礦物主要為方解石。在計算水合物地層固體基質的等效模量時，如果基質礦物各組分的體積含量和彈性模量已知，那么基質的有效彈性模量必定為Voigt上限和Reuss下限之內[15]。對于各向同性彈性基質，各個固體組分之間幾何細節未知時，Hashin-Shtrikman為最窄的上下限[23]。當兩種以上礦物組分進行混合時，基質的彈性模量必將位于廣義的Hashin-Shtrikman-Walpole界限之內[4]，即

(1)

其中

式中：K±表示各固體組分中最大和最小體積模量；μ±表示各固體組分中最大和最小剪切模量；Ki、μi和fi分別是第i種礦物組分的體積模量、剪切模量、體積含量；N為基質礦物組分所包含的種類。

(2)

式中Kma和μma分別為固體基質的體積模量和剪切模量。

2.2 干燥骨架等效模量計算

對水合物地層彈性性質的影響因素主要包括孔隙系統的微觀結構、礦物組分(影響孔隙縱橫比)和孔隙形狀等。因此，本文重點考慮孔隙的形狀和連通性對水合物地層彈性性質的影響。首先參照前人對孔隙系統的分類，將孔隙分為連通性差的硬幣狀孔隙和連通性較好的橢球型孔隙[15]；然后通過DEM(微分等效介質)模型向固體基質中加入孔隙構成干燥骨架，求取干燥骨架的彈性模量。DEM理論是通過向固體礦物相中逐漸加入包含物相模擬雙相混合物，等效體積K*和剪切模量μ*的耦合微分方程組為

(3)

式中：y為包含物相的含量；K2、μ2分別為包含物相的體積模量、剪切模量，初始條件K*(0)=K1，μ*(0)=μ1，其中K1和μ1為固體基質相的體積模量和剪切模量；P和Q為幾何因子，與包含物的形狀、體積含量有關；上標“*2”表示包含物相。由于P和Q二者與骨架的體積模量、剪切模量有關，彈性模量與微分方程是耦合的，因此需要將微分等效方程解耦為常微分方程，再進行求解。

2.3 孔隙流體等效模量計算

當孔隙充填物為水和(或)游離氣等非黏滯性流體時，流體的非黏滯性不會影響沉積物的剪切性質，即孔隙流體的剪切模量為0。但孔隙空間中的等效流體包含天然氣水合物時，由于水合物的類固體特性，導致孔隙流體的等效剪切模量不為0。此時，孔隙流體的等效彈性模量可以由體積平均的方法計算，即

(4)

式中：Kf、μf分別為孔隙流體等效體積模量、剪切模量；Sw為含水飽和度；Kw、Kh分別為水、天然氣水合物的體積模量；μh為天然氣水合物的剪切模量。

2.4 飽和巖石彈性模量計算

對于充填物為具有黏滯性特征的水合物孔隙，劉欣欣等[4]將Brown-Korringa各向異性Gassmann方程進行推廣，得到了適用于孔隙空間中包含固體的替換公式，并將該公式進行各向同性假設，得到廣義的Gassmann方程，即

(5)

式中：Ksat、μsat分別為飽和巖石的體積模量、剪切模量；Kdry、μdry分別為巖石骨架的體積模量、剪切模量；Kgr、μgr分別為基質的體積模量、剪切模量；Kφ和μφ是與孔隙有關的體積模量和剪切模量；Kif和μif是與孔隙內部填充物有關的體積模量和剪切模量。

利用廣義Gassmann方程可以得到固體(如黏滯性流體)充填于巖石孔隙情況下的飽和巖石的彈性模量。同時，廣義Gassmann方程也可以解釋孔隙流體在具有不可忽視的剪切模量時，根據等效介質理論計算得到的飽和巖石的剪切模量與干燥骨架的剪切模量差異較大。

當巖石為各向同性時，以上的模量關系為：Kgr=Kφ；μgr=μφ。當孔隙填充的是自由流體時，即μif=0。根據等效介質理論，在均勻各向同性彈性介質中，縱、橫波速度可由飽和巖石的彈性模量計算得到，即

(6)

式中ρsat為飽和巖石的密度。

3 模型測試

為了驗證本文所構建模型的合理性，定量分析泥質含量、孔隙度、水合物飽和度對含天然氣水合物地層彈性參數的影響。

3.1 泥質含量的影響

假設水合物地層的礦物組分為石英、黏土和方解石，水合物存在于孔隙空間中，流體為水，各組分的彈性模量如表1所示。控制變量參數為泥質含量，且泥質含量變化范圍為5%～70%，石英與方解石的比例為4∶1，孔隙度為45%，含水飽和度為40%。

表1 各組分彈性參數表

當參數p不斷發生變化時，模型計算的縱波速度、橫波速度、縱橫波速度比的變化如圖3所示。由圖可見，隨著泥質含量的增加，縱波速度、橫波速度減小，縱橫波速度比變化范圍反映縱、橫波速度變化合理。隨著p的增加，縱、橫波速度呈下降趨勢，這是因為p的增加使基質中的水合物含量增加、孔隙度降低，孔隙度降低的部分即是水合物被分配到基質中的體積含量；同時，水合物被分配到基質中，導致基質中原礦物的體積含量降低，因此使礦物基質的彈性模量降低。

圖3 不同比例參數情況下縱波速度(a)、橫波速度(b)、縱橫波速度比(c)隨泥質含量的變化

3.2 孔隙度的影響

假定沉積物中石英體積含量、方解石體積含量、泥質含量分別占基質礦物的50%、10%、40%，水合物的飽和度為40%，孔隙度的變化范圍為15%～60%。縱波速度、橫波速度、縱橫波速度比隨孔隙度的變化如圖4所示。由圖可見：孔隙度的增加使縱、橫波速度減小；當孔隙度大于50%時，縱、橫波速度降低的主控因素為孔隙度，而參數p的影響變小，縱橫波速度比變化范圍反映縱、橫波速度變化合理。

圖4 不同比例參數情況下縱波速度(a)、橫波速度(b)、縱橫波速度比(c)隨孔隙度的變化

3.3 水合物飽和度的影響

保持沉積物各礦物組分的體積含量不變，孔隙度為45%，水合物飽和度的變化范圍為20%～70%。縱波速度、橫波速度、縱橫波速度比隨水合物飽和度的變化如圖5所示。由圖可見，當p小于0.55時，縱波速度隨著水合物飽和度的增加而增加；當p大于0.55時，縱波速度隨著水合物飽和度的增加緩慢降低(圖5a)。這是由于水合物飽和度即使不發生改變，p的變化也會引起孔隙度、基質礦物各組分體積含量發生變化。同理，當p保持不變時，水合物飽和度的增加也會引起孔隙度、基質礦物各組分體積含量發生變化。因此從測試結果看，當p大約在0.55時，縱波速度不隨水合物飽和度發生變化。同樣，當p大約為0.34時，橫波速度不隨水合物飽和度發生變化(圖5b)。同時，縱橫波速度比變化范圍也表明了縱、橫波速度變化的合理性。

圖5 不同比例參數情況下縱波速度(a)、橫波速度(b)、縱橫波速度比(c)隨水合物飽和度的變化

4 實際數據應用

利用某海域SH2井的實測資料驗證巖石物理模型。SH2井揭示，水合物位于海底面之下195～220m，厚度約為25m，平均飽和度為25%。沉積物組分以碎屑、黏土、碳酸鹽三類礦物為主，其中，碎屑礦物以石英、斜長石、正長石、白云母為主，陸源碎屑礦物含量較高，為40.00%左右；黏土礦物以伊利石、高嶺石和蒙脫石為主，平均含量為19.64%；碳酸鹽礦物以方解石為主，平均含量為16.46%。

模型所使用的彈性模量及密度參數如表2所示，水合物壓力為5MPa，溫度為273K；礦物混合物的具體成分石英、長石、云母、黏土、方解石體積含量分別為28%、12%、26%、20%、14%。

表2 SH2井水合物地層的彈性模量和密度參數

利用現有測井資料，可計算地層孔隙度和水合物的飽和度。利用密度測井曲線計算孔隙度[24-25]，即

(7)

式中ρb、ρw和ρm分別為地層、水和骨架密度，水和骨架的密度分別取1.03g/cm3、2.65g/cm3。水合物的飽和度可以應用Archie公式[22，26]計算，可利用電阻率數據估算水合物飽和度，即

(8)

式中：Rw和Rt分別為水和地層電阻率；a、m和n均為常數，本文取a=1.3，m=n=2.0。

圖6為應用本文模型計算得到的縱、橫波速度與實際測井實測速度的對比。結果表明，本文模型計算得到的縱波速度與實測值吻合較好(圖6a)；縱波速度誤差分布符合正態分布，且在4%以內(圖7)；p的變化范圍集中在0～0.4(圖8)，均值為0.157，說明該地區水合物的微觀分布更集中于孔隙中。本文模型計算得到的橫波速度與實測橫波速度變化趨勢大體相同，吻合較好(圖6b)，進一步驗證了本文構建模型的合理性。

圖7 SH2井縱波速度相對誤差分布

圖8 SH2井p變化范圍

5 結束語

本文建立了一種基于等效介質理論的天然氣水合物地層巖石物理模型，并分析了礦物組分、孔隙度、水合物飽和度對模型的影響。通過實測縱波速度作為模型計算縱波速度的約束，可以達到不斷調整參數p的大小的目的，進而更新水合物在孔隙與基質中的分配比例。模型考慮了水合物的微觀賦存狀態，使計算得到的縱、橫波速度更加合理。將該模型應用于某海域SH2井的實際數據，結果顯示，應用效果較好，驗證了模型的適用性。

在前人方法的基礎上，本文將水合物的微觀賦存模式表達成動態形式，并考慮這種變化對孔隙度、含水飽和度、礦物組分的體積含量的影響和孔隙流體的剪切性質，構建了更為合理的巖石物理模型，為天然氣水合物地層的速度等物性參數的預測提供了方法支持。