天然氣水合物賦存地層的地震波場模擬及對比

范 金,劉洋廷,張曉波,劉 凱,鄭彥鵬*

(1.自然資源部 第一海洋研究所,山東 青島266061;2.青島海洋科學與技術試點國家實驗室 海洋地質過程與環境功能實驗室,山東 青島266237;3.國家深?；毓芾碇行?山東 青島266237)

在中低緯度地區,天然氣水合物是一種由充足的氣源、水源在高壓(一般大于10 MPa)和低溫(一般小于10℃)條件下絡合而成的籠形結晶化合物[1-2]。因水合物的氣源成分以CH4為主,另有少量CO2和H2S等氣態物質,它又被稱為“甲烷水合物”。天然氣水合物主要分布在陸地的永久凍土帶和大陸邊緣水深超過300 m的陸坡帶[3-4],在我國典型的礦藏區分別是祁連山凍土區和南海北部陸坡區。國際地質調查部門將其認定為地球上尚未大規模開發的可替代性清潔能源,對全球氣溫變化和碳循環過程有至關重要的作用[5]。

由于地質和地球物理條件的不同,天然氣水合物以不同的狀態賦存在地層中。在研究不同區域的水合物時,根據水合物與骨架顆粒的接觸關系、沉積物的粒度等要素,賦存狀態有多種劃分方法[6-9]。胡高偉等[10]發現沉積孔隙中的水合物呈混合分布,在不同的形成階段會以某種分布狀態為主導。近年來,學者發現標志天然氣水合物存在似海底反射(Bottom Simulating Reflections,BSR)、振幅空白、速度異常和極性反轉等地震異常信息[11]。高精度地震勘探技術已成為圈定天然氣水合物成礦區的主導技術,它的基本任務是通過人工地震激發的各種地震波數據反映含水合物沉積層的地質信息,實現對水合物礦藏的精準定位。地震波場的正演模擬是地震勘探的基礎。通過對不同正演模擬方法的改進,Dablain[12]、何兵壽等[13]、Liu和Sen[14]的研究成果不僅保證了計算精度,而且降低了計算成本;基于Biot-Gassmann理論,Lee[15]提出一種可預測未固結的含水合物沉積層物性特征的方法;Cordon等[16]利用正演模擬方法研究水合物儲層中的衰減特征,發現衰減導致的振幅減小可能會影響水合物體積預測的程度。隨著探測方法和技術的不斷提高,天然氣水合物在沉積層中的賦存狀態不斷被重新認識;然而,前人并沒有開展不同正演模擬方法識別不同賦存狀態水合物差異性的研究。

通過研究天然氣水合物在地層中的賦存狀態,本文總結了水合物6種不同的賦存狀態:顆粒間膠結型、顆粒包裹型、顆粒支撐型、孔隙充填型、基質與填充物共存型和結核/裂縫型,并討論了其相互轉化的關系,概述了相應物性參數的計算方法;介紹模擬含水合物沉積層地震波場的聲波模型、彈性波模型和雙相介質模型,根據水合物6種賦存狀態設定適合的孔隙度、飽和度等,基于物性參數計算方法求出地層的彈性參數并分別應用到3種正演模型中;改變沉積層的厚度、子波主頻信息,借助波場快照與地震記錄直觀地反映含水合物沉積層的地震特征,對不同賦存狀態水合物地震響應的差異性進行對比和分析,為識別水合物的勘探研究以及指導地震波反演、成像等提供依據。

1 水合物的賦存狀態概述

天然氣水合物的微觀分布主要指水合物和沉積物顆粒在地層中的接觸關系[17]。微觀分布又稱為賦存狀態。天然氣水合物在地層中的6種賦存狀態及其關系轉化見圖1。狀態1和狀態2屬于膠結狀態,將巖石顆粒視為隨機填充的球體,水合物發生在接觸點(狀態1)或在顆粒周圍生長(狀態2)。當僅存少量的水合物時,預測的彈性性能大幅增加,隨水合物濃度的逐步增加而保持相對接近[7]。在狀態2中,水合物和沉積物膠結良好,這在降低孔隙度的同時,極大改變了沉積物的剛度[18-20]。水合物在狀態3中起著顆粒支撐作用,孔隙一般被水和/或游離氣等非黏滯性流體充填,沉積物的剪切性質不受影響;水合物在固體顆粒邊緣沉積并向孔隙空間內部生長,影響巖石的剛度和熱傳導性質。值得注意的是,隨著飽和度的逐步增加(過程一),水合物可能會由狀態1轉變成狀態2,甚至是狀態3。水合物在狀態4中起填充孔隙的作用,充填物包括水合物、水或游離氣,每種組分的彈性性質各有不同,所以誘發的孔隙壓力也不盡相同;當水合物的飽和度較高時,水合物顆粒與相鄰的固體礦物顆粒相橋連,成為固體基質的一部分,增強固體骨架的力學穩定性,狀態4就可能變為狀態3,導致地層剛度變弱[21]。狀態5是分別將骨架顆粒和水合物作為基質和填充物處理,作為夾雜物嵌入均勻介質中。狀態1～狀態5都被認為天然氣水合物在沉積物中均勻分布,而水合物取芯的證據表明,水合物通常作為淺層泥質沉積物中的結核和裂縫填充物存在,該幾何形態在狀態6中示出[7]。產出于裂縫內的水合物迫使地層巖石張開形成裂縫并充填其中,一般呈網狀、結核狀或脈狀,肉眼可觀測到呈乳白色晶體狀。含水合物層不僅有較大的裂縫傾角,更受區域構造主應力控制[22]。水合物處于狀態3且在局部的飽和度逐漸升高時(過程二),就可能轉變成狀態6。如圖1中藍色箭頭所示,在實際的形成過程中,天然氣水合物可能首先是膠結沉積物顆粒(狀態1),中期以接觸沉積物顆粒或懸浮在孔隙中(狀態3或狀態4)為主,后期又以膠結形式(狀態2)為主,即多種賦存狀態共存。學者認為游離氣的存在和水合物體量增大可能是造成這種變化產生的原因[10,23]。

圖1 天然氣水合物在海底的6種賦存狀態[7,10,21-23]Fig.1 The six main occurrences of natural gas hydrate below the sea floor[7,10,21-23]

不同類型水合物的存在對沉積層的聲學特性有直接的影響,會產生不同的地震波場特征。含不同賦存狀態水合物沉積層的物性參數計算流程見圖2。當沉積層中不含水合物時,按照圖中主干部分(黑色邊框)的計算流程,代入相關的參數值,求出地層整體的彈性模量、密度、速度;當沉積層中儲存水合物時,根據不同的賦存狀態,選擇圖中分支部分(彩色邊框)相應的計算流程,求出含不同賦存狀態水合物地層的彈性模量、密度、速度。然后,根據不同的正演模型(圖2中黑色實心箭頭),確定相應的參數。因篇幅所限,本文對計算公式不做具體介紹,詳見文獻[6,22,24-25]。

圖2 含水合物沉積層的物性參數計算流程Fig.2 Flow chart for calculating physical parameters of hydrate-bearing sedimentary layer

2 正演模擬的理論模型

地震波傳播理論是各類地震學科的基礎,核心內容之一就是地震波傳播和介質模型[26]。彈性波動方程是地震波傳播的動力學核心,解釋了地震波在復雜地層中的傳播機理。不同的模型中介質的速度、相態、密度等參數的計算方法各有不同,主要通過孔隙度和飽和度來界定,進而影響地震波的傳播。本文主要使用二維各向同性的聲波模型、彈性波模型和雙相介質(即孔隙流體介質)模型,其波動方程詳見文獻[27-29]。在地震波場的模擬方法中,有限差分法因實現簡單、計算效率高等優勢被應用于多項研究。董良國等[30]將高階差分法和交錯網格相結合,這一方法目前在地震波場數值模擬研究中被廣泛使用。除此之外,差分格式的相容性、收斂性、穩定性和邊界條件等問題也已被前人研究[27-28,30-33]。

參考Ecker對海底沉積物的參數設定[20],本文模擬的沉積層中各礦相的物性參數值見表1,其中沉積骨架由5%的石英、35%的方解石和60%的黏土共同構成[34]。為對比3種正演模型的差異,設計了一種雙層模型(圖3),物性參數見表2;其中,上部是孔隙度為35%的飽水沉積層,下部是孔隙度為20%的飽水沉積層。計算模擬的波場空間為1 500 m×1 500 m,網格大小為5 m×5 m,時間步長為0.5 ms。采用交錯網格有限差分法,精度為時間2階、空間10階。震源采用主頻為25 Hz的雷克子波。

圖3 雙層模型Fig.3 Two-layer model

表1 沉積物中各礦相的物性參數Table 1 Physical property parameters of various mineral phases in sediments

表2 雙層模型的物性參數表Table 2 Physical property parameters of two-layer model

圖4 是0.25 s時刻地震波在3種模型中傳播的網格化波場快照,展現了清晰的差異性:激發純縱波在聲波介質中波阻抗界面只存在透射和反射;當地震波到達彈性介質的波阻抗界面時,除了產生透射和反射,還會產生轉換波;區別于彈性介質,震源在雙相介質中會產生與快縱波相互伴生的慢縱波。慢縱波在實際地層傳播過程中具有極大的衰減性,在波阻抗界面也會產生更多類型的地震波。觀察波場快照還可知,基于等效彈性模量模型和Biot理論計算的雙相介質參數,模擬的縱波走時與聲波介質、彈性介質中縱波走時存在明顯的差異,主要是由不同的拉梅系數計算方法導致的。

圖4 0.25 s時刻地震波在3種模型中的網格化波場快照Fig.4 Meshed snapshots of seismic waves in the three models at time 0.25 s

3 波場模擬及對比分析

根據4個不同的層狀模型(圖5),設計了合理的地層孔隙度、不同賦存狀態下水合物的飽和度以及游離氣的飽和度,并依據圖2的流程來計算相應層位的物性參數(表3)。采用與前文相同的波場空間、采樣間隔和震源類型。在海水中激發純縱波,炮點位于近海面處,檢波器排列于炮點下方,61道接收,道間距為25 m。采用PML吸收邊界,只考慮地震波在介質中的傳播。對含不同賦存狀態水合物的沉積層進行了大量的正演模擬分析,并挑選出幾種含水合物沉積層的波場快照和地震記錄,展開對比與分析。

層狀模型1從上往下分別是海水層、飽水沉積層、含孔隙充填型水合物沉積層和飽水沉積層(圖5a)。圖6～圖8給出了低頻(子波主頻采用25 Hz)條件下,在聲波介質、彈性介質和雙相介質分別模擬不含游離氣的孔隙充填型水合物沉積層的典型波場快照、去除直達波的單炮地震記錄和含水合物層底界振幅隨偏移距變化的曲線。1)對比波場快照(圖6)可知,地震波在聲波介質中的反射波能量要明顯強于彈性介質和雙相介質,主要是因為縱波在彈性介質中向上部地層反射時,遇到了向下部地層透射的橫波,產生能量耗散;在雙相介質中,地震波不僅發生能量耗散,且孔隙流體的存在導致地震波將能量按孔隙度大小分配給了固、流相,反射振幅進一步減小。2)對比地震記錄(圖7)可知,當含水合物層下伏的沉積層中不含游離氣時,聲波介質和彈性介質模擬含水合物層底界只表現出與海底反射極性相同、與水合物層頂界反射能量相近的反射;雙相介質模擬含水合物層頂界的反射微弱,但底界的負極性反射異常明顯。3)從含水合物層底界振幅隨偏移距變化的曲線(圖8)可知,由于所提取曲線受到了其他波的干擾,造成個別位置產生局部跳動,但不影響整體的變化。隨著偏移距的增大,聲波介質模擬的反射振幅明顯增大,彈性介質模擬的反射振幅整體上減小,而雙相介質模擬的底界振幅是先增大后減小。

圖5 正演模擬與對比所需的層狀模型Fig.5 Forward simulation and layered model for the purpose of comparison

表3 基于等效彈性模量模型和Biot理論計算的層狀模型物性參數Table 3 Physical property parameters of layered model calculated based on equivalent elastic modulus model and Biot theory

圖6 模擬孔隙充填型水合物的網格化波場快照Fig.6 Snapshots of gridded wavefields that simulate pore-filling hydrates

圖7 模擬孔隙充填型水合物的去除了直達波的單炮地震記錄Fig.7 A single shot seismic record that simulates pore-filling hydrates with direct waves removed

圖8 模擬含孔隙充填型水合物層底界振幅隨偏移距變化的曲線Fig.8 The variation curve of the base amplitude of the simulated pore-filling hydrate layer with offset

層狀模型2從上往下分別是海水層、飽水沉積層、含顆粒包裹型水合物沉積層、含游離氣層和飽水沉積層(圖5b)。低頻(子波主頻采用25 Hz)條件下,在3種介質分別模擬含游離氣的顆粒包裹型水合物沉積層的正演結果。圖9～圖11分別展示了正演模擬得到的波場快照、單炮地震記錄和含水合物層底界振幅隨偏移距變化的關系。1)從波場快照(圖9)可知,地震波在3種介質中反射的能量大小關系與圖6相似,且含水合物高速層與含游離氣低速層的界面反射與海底反射極性相反。2)對比地震記錄(圖10)可以看出,3種介質模擬水合物底界都有清晰的BSR響應[35],與波場快照的現象吻合;彈性介質和雙相介質模擬含水合物層底界的振幅在紅色邊框處突然變小,而聲波介質模擬不存在這一現象。究其原因,猜測是由于臨界角的變化或者是極性相反的波形相疊加造成的。3)從振幅隨偏移距變化的曲線(圖11)可知,隨著偏移距的增大,聲波介質模擬水合物層底界的反射振幅整體上先保持穩定,后快速增大;彈性介質模擬水合物底界的振幅先緩慢減小,后產生震蕩式變化;雙相介質模擬水合物底界振幅整體上呈震蕩式減小。與圖6～圖8相比,當含水合物層下伏沉積層中含游離氣時,無論使用哪一種正演模型,水合物層底界的負極性反射都尤為明顯。區別于圖9～圖11,將子波主頻提高為40 Hz后,地震波形變窄。聲波介質模擬水合物底界的反射振幅連續性較好,彈性介質和雙相介質模擬水合物層底界的同相軸出現振幅突然減小的特征。相比于圖9～圖11,提高子波主頻對3種正演模型識別含水合物沉積層的影響不大。

圖9 模擬顆粒包裹型水合物的網格化波場快照Fig.9 Snapshots of gridded wavefields that simulate particle-encapsulated hydrates

圖10 模擬顆粒包裹型水合物的去除了直達波的單炮地震記錄Fig.10 A single shot seismic record that simulates particle-encapsulated hydrates with direct waves removed

圖11 模擬含顆粒包裹型水合物層底界振幅隨偏移距變化的曲線Fig.11 The variation curve of the base amplitude of the simulated particle-encapsulated hydrate layer with offset

層狀模型3從上往下分別是海水層、飽水沉積層、含顆粒間膠結型水合物薄層、含游離氣層和飽水沉積層(圖5c)。低頻(子波主頻采用25 Hz)條件下,在3種介質分別模擬含游離氣的顆粒間膠結型水合物薄層的正演結果。圖12～圖14分別展示了正演模擬得到的波場快照、單炮地震記錄和含水合物層底界振幅隨偏移距變化的關系。1)從波場快照(圖12)可以看出,透射縱波在聲波介質中傳播時的能量最強,彈性介質中次之,雙相介質中最弱。2)由地震記錄(圖13)可知,當含水合物沉積層的厚度較薄時,聲波介質模擬含水合物層頂界的反射較為明顯;在彈性介質中可以清晰地看到薄層頂底面的反射,且縱波走時相差很小,兩種極性相反的同相軸相互干涉;而雙相介質中能在含水合物薄層與含游離氣層之間看到明顯的負極性反射,卻幾乎看不到薄層頂界的反射。3)從振幅隨偏移距變化的曲線(圖14)可以看出,隨著偏移距的增大,聲波介質模擬薄層底界的振幅逐漸增大,而彈性介質和雙相介質模擬的振幅都逐漸減小。

圖12 模擬顆粒間膠結型水合物的網格化波場快照Fig.12 Snapshots of gridded wavefields that simulate interparticle cementation hydrates

圖13 模擬顆粒間膠結型水合物的去除了直達波的單炮地震記錄Fig.13 A single shot seismic record that simulates interparticle cementation hydrates with direct waves removed

圖14 模擬含顆粒間膠結型水合物層底界振幅隨偏移距變化的曲線Fig.14 The variation curve of the base amplitude of the simulated interparticle cementation hydrate layer with offset

將層狀模型3(圖5c)中含水合物薄層中水合物的類型改為顆粒支撐型后進行分析,得到高頻(子波主頻采用40 Hz)條件下,在3種介質中分別模擬含游離氣的顆粒支撐型水合物薄層的正演結果(圖15～圖17)。1)從波場快照(圖15)可知,在提高子波主頻后,地震波形變窄。聲波介質和彈性介質模擬水合物薄層頂底面的反射可以被區分。與之相比,雙相介質模擬水合物薄層頂底界的反射能量要弱很多,不易被識別。2)從地震記錄(圖16)可知,聲波介質和彈性介質模擬的水合物薄層出現了相互疊加的同相軸,但不影響薄層頂底面的識別。3)從振幅隨偏移距變化的曲線(圖17)可知,隨著偏移距的增大,聲波介質和彈性介質模擬薄層底界的振幅整體上逐漸減小,而雙相介質模擬的振幅整體是增加的。與圖12～圖14相比,圖15～圖17中聲波模型對水合物薄層頂底面的識別能力在高頻條件下得以增強,彈性波模型識別水合物薄層的能力差別不大,而雙相介質模型對其識別能力在高頻條件下得以削弱。

圖15 模擬顆粒支撐型水合物的網格化波場快照Fig.15 Snapshots of gridded wavefields that simulate particle-supported hydrates

圖16 模擬顆粒支撐型水合物的去除了直達波的單炮地震記錄Fig.16 A single shot seismic record that simulates particle-supported hydrates with direct waves removed

圖17 模擬含顆粒支撐型水合物層底界振幅隨偏移距變化的曲線Fig.17 The variation curve of the base amplitude of the simulated particle-supported hydrate layer with offset

層狀模型4從上往下分別是海水層、飽水沉積層、含結核型水合物層和飽水沉積層(圖5d)。低頻(子波主頻采用25 Hz)條件下,在3種介質分別模擬含結核型水合物沉積層的正演結果(圖18～圖20)。多個聚集成塊狀體的純水合物(以下稱“結核”)分布在沉積層中,形成多個波阻抗界面。1)從波場快照(圖18)可知,3種介質模擬波場的能量大小關系和前文模擬的相似,在水合物結核處出現明顯的振幅空白特征,主要是因為水合物結核內部阻抗差異較小,反射較弱;各種地震波之間相互干涉,波場異常復雜。2)從地震記錄(圖19)可知,由于水合物下部不存在游離氣,聲波介質模擬水合物底界的反射呈正極性,與圖7a模擬的結果一致;彈性介質模擬水合物底界的反射呈負極性,這與圖7b模擬的結果相反;雙相介質模擬水合物底部的反射與彈性介質模擬的規律一致,與圖7c模擬的結果相似。3)從振幅隨偏移距變化的曲線(圖20)可知,3種介質模擬水合物底界的振幅變化有著相似的規律,且除了雙相介質模擬的在大偏移距處能量較小之外,其他位置基本一樣。

圖18 模擬結核型水合物的網格化波場快照Fig.18 Snapshots of gridded wavefields that simulate nodule hydrates

分析上述試驗,對比結果見表4。其中,3種介質模擬水合物層底界的反射振幅大小整體高于10-3數量級的標記為強振幅,否則,為弱振幅。在利用聲波介質與彈性介質模擬上層含20%水合物且下層不含游離氣的波場時,孫小芳[36]發現含水合物層底界基本不產生負極性的BSR現象且振幅較小;相同的地質模型下,吳志強等[37]也認為很難看到明顯的似海底反射現象。本文圖6～圖8中設定的地質模型和水合物飽和度與之相似,得到的試驗結果也基本一致。當下伏介質含少量游離氣且上層介質含適量的水合物時,含水合物層底界就會產生強振幅的負極性反射[36],與本文圖9～圖11中得到的試驗結果基本一致。吳志強等[37]在研究水合物飽和度的變化對沉積層頂底界振幅的影響時,發現在其他條件相同的情況下,水合物飽和度達到40%以上將使含水合物層的頂界出現明顯的反射。同時,該研究還認為提高子波主頻對分辨薄層水合物十分有效。這兩點在本文圖12～圖14和圖15～圖17的對比中都得到了驗證。

表4 三種介質模擬不同類型水合物地震波場的對比Table 4 Comparison of three media simulating different types of hydrate seismic wave fields

4 結 論

本文總結介紹了天然氣水合物在海洋中的6種主要賦存狀態類型,概述了每種賦存狀態之間的相互轉化關系及其物性參數計算方法,并基于這些方法,求取了層狀模型的各項參數。通過改變沉積層的厚度、子波主頻信息,利用3種正演模擬方法對含不同賦存狀態的水合物沉積層進行波場模擬,對比分析波場快照、地震記錄的差異性,得到以下結論:

1)當地層中存在孔隙充填型水合物且下伏地層不含游離氣時,聲波介質和彈性介質模擬含水合物層底界表現出與海底反射極性相同的弱反射,而雙相介質模擬含水合物層底界表現出負極性特征。當地層中充填結核型水合物且下伏地層不含游離氣時,3種介質模擬水合物底界反射能量相近。聲波介質模擬水合物底界依然表現出與海底反射極性相同的弱反射,而彈性介質和雙相介質模擬水合物底界呈負極性。

2)當地層中充填顆粒包裹型水合物且下伏地層含游離氣時,在低頻條件下,3種介質模擬含水合物層底界的負極性反射明顯;隨著偏移距的增大,聲波介質模擬水合物層底界的振幅先趨于穩定,后快速增大,而彈性介質和雙相介質模擬的振幅呈震蕩式變化。提高子波主頻,地震波形較窄,水合物層底界反射依然明顯。

3)當沉積薄層中充填顆粒間膠結型水合物且下伏地層含游離氣時,不同于聲波介質,彈性介質和雙相介質模擬水合物薄層底界的反射振幅隨偏移距的增大而減小;將水合物類型改為顆粒支撐型并提高子波主頻,不同于雙相介質,聲波介質和彈性介質模擬水合物薄層底界的反射振幅較強,且隨著偏移距的增大而減小。

本文使用的正演模型是二維各向同性的,而實際的地層情況較為復雜,建議將正演模型推廣到三維的、各向異性的情況。同時,還需要不斷優化算法,提高計算效率,設計更多符合實際地層的模型進行測試。文中的研究方法尚未對實際地層中的水合物數據進行應用,建議對實際數據進行應用,以利于后續研究。