藍(lán)坤 朱賀 何濤,? 梁前勇 吳學(xué)敏 董一飛 張毅

天然氣水合物相關(guān)的Slipstream海底滑坡體速度結(jié)構(gòu)模型反演

藍(lán)坤1,2朱賀1,2何濤1,2,?梁前勇3吳學(xué)敏3董一飛3張毅4

1.造山帶與地殼演化教育部重點(diǎn)實(shí)驗(yàn)室, 北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871; 2.北京天然氣水合物國(guó)際研究中心, 北京 100871; 3.中國(guó)地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局, 廣州 510760; 4.中國(guó)地質(zhì)科學(xué)院, 北京 100037; ?通信作者, E-mail: taohe@pku.edu.cn

針對(duì)目前用于建立二維速度結(jié)構(gòu)模型的 RAYINVR 軟件對(duì)四分量海底地震儀(OBS)記錄的轉(zhuǎn)換橫波無(wú)法自動(dòng)反演, 建模效率低的問(wèn)題, 采用 MATLAB 遺傳反演算法, 對(duì) RAYINVR 軟件進(jìn)行改進(jìn), 實(shí)現(xiàn)對(duì)橫波速度結(jié)構(gòu)模型各層和層中區(qū)塊泊松比的自動(dòng)同步反演, 可以為天然氣水合物勘探調(diào)查提供泊松比和楊氏模量等重要彈性參數(shù)。利用采自 Slipstream 海底滑坡的 OBS 數(shù)據(jù), 通過(guò)同步反演方法, 獲得該地區(qū)較為精細(xì)的縱波和橫波速度結(jié)構(gòu)模型, 并可與附近 U1326 鉆孔的測(cè)井?dāng)?shù)據(jù)進(jìn)行對(duì)比, 驗(yàn)證了橫波速度結(jié)構(gòu)同步反演建模方法的有效性。模型揭示了兩個(gè)泊松比差異大的結(jié)構(gòu)面: 代表水合物穩(wěn)定底界的似海底反射界面(BSR)(海底之下230±10m)以及淺部異常高速體(有可能是高飽和度水合物富集的砂體)的底界(海底之下 75～100m), 后者與滑脫面大致重合, 指示水合物與Slipstream海底滑坡的形成有關(guān)。

同步反演; 橫波速度結(jié)構(gòu); Slipstream海底滑坡; 天然氣水合物; 射線追蹤

地震勘探是海域天然氣水合物調(diào)查的重要途徑[1?2], 不但可以對(duì)水合物儲(chǔ)層進(jìn)行成像刻畫(huà), 還可以建立速度結(jié)構(gòu)模型, 從而有效地評(píng)估水合物資源量, 對(duì)后續(xù)開(kāi)采意義重大。常規(guī)的海面拖纜地震勘探僅能產(chǎn)生和接收單一的縱波信號(hào), 對(duì)沉積物孔隙流體中的氣體較為敏感, 少量的游離氣就可能造成大幅度的縱波能量衰減和速度降低[3], 影響水合物地層的成像效果和飽和度估算準(zhǔn)確度。這種現(xiàn)象在作為水合物穩(wěn)定區(qū)底界標(biāo)志的地震似海底反射層(bottom simulating reflector, BSR)上下(表現(xiàn)為雜亂強(qiáng)反射)或攜帶大量游離氣向上運(yùn)輸?shù)睦淙ǖ?氣煙囪)內(nèi)部(表現(xiàn)為條帶狀空白)尤為嚴(yán)重[4]。

橫波主要通過(guò)沉積物骨架傳播, 幾乎不受孔隙流體的影響[5]。因此, 在海域水合物重點(diǎn)勘探區(qū), 通常使用四分量海底地震儀(Ocean Bottom Seismo-meter, OBS)來(lái)記錄由縱波在地層界面上轉(zhuǎn)換而來(lái)的橫波信號(hào)[6?7], 從而建立沉積地層的縱橫波速度結(jié)構(gòu)模型, 以便獲取更加豐富、更加準(zhǔn)確的水合物儲(chǔ)層信息。

加拿大西海岸溫哥華島外的北卡斯卡底(Casca-dia)陸緣增生楔主要由板塊俯沖過(guò)程中刮削下來(lái)的洋盆沉積物堆積而成[8?9], 蘊(yùn)藏著豐富的天然氣水合物資源[10?11], 20 世紀(jì) 80 年代以來(lái)一直是水合物研究的重點(diǎn)區(qū)域, 不但具有豐富的單道地震(single channel seismic, SCS)、多道地震和 OBS 等地球物理勘探數(shù)據(jù)[12?13], 而且有專(zhuān)門(mén)針對(duì)水合物的大洋鉆探計(jì)劃(Ocean Drilling Program, ODP) 146 航次[14]和綜合大洋鉆探計(jì)劃(Integrated Ocean Drilling Program, IODP)X311 航次[15?17]的鉆井?dāng)?shù)據(jù)。特別之處是, 在該增生楔前緣的條狀變形脊上發(fā)育一系列海底滑塌構(gòu)造。Lopez 等[18]和 Yelisetti 等[19]分別運(yùn)用 OBS 數(shù)據(jù)和 RAYINVR 軟件[20], 研究其中規(guī)模較大的 Orca和 Slipstream 海底滑塌的地震速度結(jié)構(gòu), 認(rèn)為其形成與水合物富集層的位置有關(guān)。然而, 受限于轉(zhuǎn)換橫波處理的復(fù)雜性和 RAYINVR 軟件的缺陷, 他們只建立滑塌構(gòu)造的縱波速度模型, 缺乏橫波速度信息, 無(wú)法提供滑塌構(gòu)造力學(xué)分析所需的楊氏模量和泊松比等彈性參數(shù)。

廣泛使用的 RAYINVR 軟件是用射線追蹤方法對(duì)二維速度模型進(jìn)行走時(shí)擬合, 其內(nèi)置的阻尼最小二乘(damped least squares)反演算法能夠方便地建立純縱(橫)波速度結(jié)構(gòu)模型。但是, 對(duì)于 OBS 記錄的轉(zhuǎn)換橫波反射事件, 該軟件只能通過(guò)人工改變模型的泊松比參數(shù), 將縱波速度轉(zhuǎn)換為相應(yīng)的橫波速度, 并且, 一旦存在橫向速度結(jié)構(gòu)不均的現(xiàn)象, 還需要將該結(jié)構(gòu)層進(jìn)一步劃分區(qū)塊, 模型會(huì)變得更加復(fù)雜。因此, 橫波模型需要花費(fèi)大量的時(shí)間進(jìn)行反復(fù)的手動(dòng)測(cè)試, 效率過(guò)低, 嚴(yán)重妨礙 OBS 橫波信號(hào)的處理。

針對(duì)上述問(wèn)題, 本研究在 RAYINVR 軟件的基礎(chǔ)上, 依托 MATLAB 語(yǔ)言建立 OBS 橫波模型的遺傳反演方法, 實(shí)現(xiàn)各結(jié)構(gòu)層和層中區(qū)塊參數(shù)(泊松比)的同步反演。利用該方法建立的較精細(xì)的 Slips-tream 海底滑坡縱、橫波速度結(jié)構(gòu)模型, 指示淺部高速異常體(海底之下 75～100m)和深部 BSR (海底之下 230±10m)兩個(gè)軟弱面為潛在的滑脫面位置。

1 RAYINVR軟件

對(duì) OBS 等廣角地震數(shù)據(jù), 通常采用射線路徑追蹤的算法進(jìn)行速度結(jié)構(gòu)層析反演, 采用正演試錯(cuò)擬合的求解方案對(duì)模型參數(shù)進(jìn)行迭代[21?22], 直至模型計(jì)算的走時(shí)與真實(shí)記錄的地震信號(hào)走時(shí)之間誤差足夠小。對(duì)于二維速度結(jié)構(gòu)模型, 目前普遍使用 Zelt等[20,23]編寫(xiě)的 RAYINVR 軟件, 其模型構(gòu)建采用不規(guī)則四邊形網(wǎng)格, 通過(guò)靈活地設(shè)置速度以及界面節(jié)點(diǎn)的個(gè)數(shù)和位置, 使得模型層界面與地震事件識(shí)別情況一致, 能夠在很大程度上避免過(guò)參數(shù)化引起的反演假象以及欠參數(shù)化引起的細(xì)結(jié)構(gòu)丟失問(wèn)題。模型可以分為若干個(gè)不規(guī)則四邊形, 預(yù)先定義每個(gè)四邊形的頂點(diǎn)速度, 然后采用雙線性內(nèi)插法, 由頂點(diǎn)速度得到四邊形內(nèi)部的速度, 因此在四邊形內(nèi)部的射線追蹤存在解析解。該軟件可以對(duì)折射波、反射波、界面滑行波、轉(zhuǎn)換波和層間反射波等進(jìn)行射線追蹤計(jì)算, 具有運(yùn)算速度快和算法穩(wěn)定等優(yōu)點(diǎn)。

當(dāng)沒(méi)有轉(zhuǎn)換波時(shí), RAYINVR 軟件可以使用自帶的阻尼最小二乘反演工具(DMPLSTSQR 程序), 從上往下逐層對(duì)模型節(jié)點(diǎn)速度和(或)層面深度參數(shù)進(jìn)行自動(dòng)的快速的迭代反演, 每一層的射線走時(shí)擬合流程見(jiàn)圖 1。以建立縱波速度結(jié)構(gòu)模型為例, 首先用 RAYINVR 軟件正演計(jì)算某層的射線路徑和走時(shí), 需要輸入 v.in, tx.in 和 r.in 這 3 個(gè)文件。其中, 文件 v.in 包含模型每一層的深度和速度信息, 文件tx.in 是從地震成像數(shù)據(jù)拾取的真實(shí)地震事件走時(shí), 文件r.in則是關(guān)于射線追蹤的算法控制以及結(jié)果畫(huà)圖的參數(shù)設(shè)置。RAYINVR 正演會(huì)輸出射線追蹤走時(shí)與實(shí)際地震事件走時(shí)殘差的均方根(root mean squares, RMS)和標(biāo)準(zhǔn)卡方值(normalized2), 后者的計(jì)算公式為

其中,為實(shí)際拾取的地震走時(shí)數(shù)量,T為第個(gè)實(shí)際拾取的地震走時(shí),t為相應(yīng)的射線追蹤得到的地震走時(shí),為第個(gè)地震走時(shí)的拾取誤差。在得到標(biāo)準(zhǔn)卡方值后, 通過(guò)人機(jī)交互判斷2是否滿足要求: 當(dāng)2≤1 時(shí), 說(shuō)明走時(shí)殘差達(dá)到地震事件拾取誤差水平, 可以結(jié)束反演; 如果2顯著大于 1, 則將該值輸入 DMPLSTSQR 程序, 同時(shí)輸入反演參數(shù)控制文件d.in (在該文件中選擇需要反演的目標(biāo)界面), 通過(guò)偏微分方法求出模型參數(shù)調(diào)整向量, 據(jù)此更新模型文件 v.in 中的速度和(或)深度參數(shù)。迭代上述正演和反演過(guò)程, 直到獲得模型最優(yōu)解(保存在最后一步的文件 v.in 中)。

可見(jiàn), DMPLSTSQR 程序是通過(guò)改變文件 v.in 中的速度和地層界面深度參數(shù)來(lái)實(shí)現(xiàn)迭代反演的。然而, 當(dāng)對(duì)轉(zhuǎn)換橫波進(jìn)行射線追蹤時(shí), 縱波轉(zhuǎn)換為橫波的反射界面參數(shù)(如轉(zhuǎn)換發(fā)生的層號(hào)和泊松比等)均存儲(chǔ)在文件 r.in 中, 無(wú)法由 DMPLSTSQR 程序自動(dòng)更新, 只能手動(dòng)修改。因此, 橫波速度模型的迭代反演過(guò)程非常耗時(shí), 并且, 當(dāng)研究區(qū)的地質(zhì)構(gòu)造比較復(fù)雜時(shí), 需要對(duì)目標(biāo)層進(jìn)一步細(xì)化分塊, 以至很難手動(dòng)將模型參數(shù)調(diào)節(jié)到令人滿意的程度。

圖1 RAYINVR軟件迭代反演速度模型流程 (修改自文獻(xiàn)[20])

2 轉(zhuǎn)換橫波同步反演方法

根據(jù)震動(dòng)形式, 可將地震波分為體波和面波, 其中體波又分為縱波和橫波。當(dāng)縱波到達(dá)一個(gè)反射界面時(shí), 其反射波中不僅包含縱波, 還包含轉(zhuǎn)換而來(lái)的橫波(圖 2)。由于海面震源激發(fā)的地震波在海水中只能傳播縱波, 因此放置在海底的 OBS 記錄的主要橫波事件(能量最強(qiáng), 可以識(shí)別)是由縱波在其反射界面上轉(zhuǎn)換產(chǎn)生。因此, RAYINVR 對(duì) OBS 橫波的射線追蹤, 不是直接確定橫波速度, 而是在給定的縱波速度結(jié)構(gòu)模型(地層界面深度和縱波速度梯度保持固定)基礎(chǔ)上, 尋找最優(yōu)的地層泊松比。在通過(guò) RAYINVR 軟件得到模型的泊松比之后, 橫波速度S則由對(duì)縱波速度P的變換得到:

如果研究區(qū)域的地質(zhì)構(gòu)造較為簡(jiǎn)單, 可以通過(guò)手動(dòng)調(diào)節(jié)的方法, 給模型各層設(shè)置泊松比參數(shù)。但是, 當(dāng)?shù)刭|(zhì)構(gòu)造比較復(fù)雜時(shí), 地層的泊松比不僅在縱向上可能存在差異, 在橫向上也可能存在較大的差異, 需要在層內(nèi)進(jìn)一步劃分區(qū)塊。此時(shí), 用人工調(diào)節(jié)參數(shù)的方法就會(huì)事倍功半, 因此需要研發(fā)一套能夠讓程序?qū)崿F(xiàn)自動(dòng)同步反演而得到最優(yōu)解的建模方法。

1為P波在反射界面的入射角,2為P波的折射角,1為轉(zhuǎn)換S波的反射角,2為轉(zhuǎn)換S波的折射角

圖2 縱波在彈性界面發(fā)生反射和折射示意圖

Fig. 2 Reflection and refraction diagram of P wave incident on an elastic interface

遺傳算法通過(guò)模擬自然進(jìn)化過(guò)程, 對(duì)值域空間進(jìn)行搜索, 得到全局最優(yōu)解, 是一種可用于復(fù)雜系統(tǒng)優(yōu)化的算法。本文將 RAYINVR 軟件的射線追蹤法與 MATLAB 軟件的遺傳算法相結(jié)合, 用 MATLAB語(yǔ)言編寫(xiě)一套能夠?qū)崿F(xiàn) OBS 橫波速度結(jié)構(gòu)模型在各層以及層間不同區(qū)塊同步反演的方法, 從而節(jié)省人工修改參數(shù)的時(shí)間成本, 實(shí)現(xiàn)復(fù)雜地質(zhì)構(gòu)造環(huán)境下反演橫波速度結(jié)構(gòu)的自動(dòng)化。該同步反演方法的基本流程見(jiàn)圖 3。

1)給定模型各層和層中區(qū)塊的初始泊松比組合(保存在文件 r.in 中)。

2)設(shè)定遺傳算法的最大代數(shù)和每代的種群數(shù)量, 啟動(dòng)自動(dòng)反演。

3)RAYINVR 軟件讀入射線追蹤控制文件 r.in、拾取的地震走時(shí)文件 tx.in 以及速度和層界面深度控制文件 v.in, 通過(guò)正演射線路徑擬合得到2。

圖3 OBS橫波速度模型同步反演流程

4)遺傳算法會(huì)判斷此時(shí)是否達(dá)到預(yù)設(shè)條件: 如果達(dá)到預(yù)設(shè)的迭代次數(shù), 或者2≤1 (意味著此時(shí)的泊松比為全局最優(yōu)解), 反演終止; 如果沒(méi)有達(dá)到預(yù)設(shè)的迭代次數(shù), 并且2>1, 那么程序會(huì)生成一組新的泊松比組合, 并自動(dòng)傳給 RAYINVR 軟件, 進(jìn)行新一輪正演。

5)重復(fù)上述迭代過(guò)程, 直到滿足預(yù)設(shè)條件, 得到的最新泊松比組合將保存在文件 r.in 中。

3 Slipstream海底滑坡區(qū)地質(zhì)背景

在卡斯卡底活動(dòng)陸緣, 胡安德富卡(Juan de Fuca)板塊以約 46mm/a 的速率俯沖到北美板塊之下, 來(lái)自洋盆的大量沉積物在變形前緣位置被刮削下來(lái), 堆積形成增生楔[8]。這些增生體沿著一些巨大的深達(dá)洋殼的逆沖斷層發(fā)生擠壓和縮短, 并相互重疊, 從而在陸坡底部(變形前緣)形成一系列與板塊邊界近平行的長(zhǎng)條狀背斜山脊[19], 本文的研究對(duì)象 Slip-stream 山脊就是其中之一(圖 4)?？ㄋ箍ǖ谆顒?dòng)陸緣一直是天然氣水合物研究的聚焦區(qū)域, 因?yàn)榇罅康卣鹫{(diào)查都觀測(cè)到從變形前緣到陸坡中部的深水沉積物中廣泛存在作為水合物穩(wěn)定區(qū)底界地震標(biāo)志的BSR, 海底之下的深度約為 200～300m, 顯示大量天然氣水合物的富集[5,24?26,28], 也得到 ODP146 和 IODP X311 大洋鉆探結(jié)果的證實(shí)。

近年來(lái), 在高分辨率多波束海底地貌的調(diào)查中發(fā)現(xiàn), 變形前緣的條狀背斜山脊上經(jīng)常發(fā)育滑塌構(gòu)造。例如, IODP X311 U1326 站位所在的山脊發(fā)育一個(gè)規(guī)模較大的 Orca 海底滑塌(圖 4), 其滑脫面深度恰好抵達(dá) BSR, 因此推測(cè)該滑塌構(gòu)造的形成與水合物層有關(guān)[18]。作為對(duì)比, 其相鄰的 Slipstream 滑塌(U1326 站位東南方向約 15km)規(guī)模更大, 推測(cè)其形成也與水合物有關(guān), 但滑脫面的深度卻只有 100m左右[19]。

根據(jù)美國(guó)華盛頓大學(xué) TN175 航次采集的 EM300多波束測(cè)深數(shù)據(jù)[19], Slipstream 海底滑坡體的長(zhǎng)度約為 3.3km, 寬度約為 2.5km, 滑坡角度較大(約為25°), 出露于海底的陡壁高度可達(dá)近 100m。該滑坡體上發(fā)育一系列正斷層, 且部分?jǐn)鄬用娉雎队诤５? 因此滑坡體在地貌上類(lèi)似多個(gè)排列不整齊的梯形箱狀高地。較大的滑坡塊體主要分布在距離山腳 1～2km 的海底平原沉積層之上, 一些粒度較小的滑坡體沉積物則外流到更遠(yuǎn)的地區(qū)(圖 4)。本研究采用的 SCS 和 OBS 數(shù)據(jù)來(lái)自 SeaJade (Seafloor Earthquake Array – Japan Canada Cascadia Experiment)項(xiàng)目。

修改自文獻(xiàn)[27]。紅色虛線框?yàn)?Slipstream 海底滑坡區(qū)域; 紅色三角形為 OBS 測(cè)點(diǎn) 30, 31, 33 和 34 的位置; 4 號(hào)地震測(cè)線平行地穿過(guò)Slipstream 山脊, 其中綠線為 OBS 地震航次, 黃線為 SCS 航次; 紅色五角星為 IODP X311 航次的U1326 站位; 多波束海底地貌圖上疊加的等深線間隔為 25 m

4 縱波速度結(jié)構(gòu)模型

由于 OBS 記錄的橫波是由縱波在地層界面上轉(zhuǎn)換而來(lái), 因此需要建立 Slipstream 海底滑坡區(qū)域的縱波速度結(jié)構(gòu)模型。

傳統(tǒng)的海面拖纜 SCS 數(shù)據(jù)能夠很好地針對(duì)海底沉積層的結(jié)構(gòu)特點(diǎn)進(jìn)行成像, 但受限于炮點(diǎn)與接收器之間的偏移距(source-receiver offset)在深海中太小, 無(wú)法得到準(zhǔn)確的速度和深度信息。布設(shè)于海底的 OBS 采集的是廣角(大偏移距)數(shù)據(jù), 能夠很好地對(duì)海底沉積層的速度信息進(jìn)行分析, 但是對(duì)其結(jié)構(gòu)信息的表征則有所欠缺。因此, 在構(gòu)建海底沉積層縱波速度結(jié)構(gòu)模型時(shí), 我們同時(shí)采用 OBS 廣角數(shù)據(jù)和 SCS 數(shù)據(jù)來(lái)實(shí)現(xiàn)對(duì)地層速度和深度的同步約束。需要注意的是, RAYINVR 軟件最初是用于大陸勘探的, 設(shè)計(jì)的是單震源對(duì)應(yīng)多個(gè)接收器, 而海洋的地震勘探則是單個(gè) OBS 記錄多個(gè)炮點(diǎn)數(shù)據(jù), 因此在建立模型時(shí), 需要將 OBS 作為射線發(fā)射點(diǎn), 將炮點(diǎn)作為接收位置。

前人在卡斯卡底增生楔區(qū)域創(chuàng)建的縱波模型一般為 5～7 層[5,19,25], 用于滑塌區(qū)力學(xué)機(jī)制分析時(shí)精細(xì)度不夠。我們通過(guò)仔細(xì)地辨識(shí) OBS 和 SCS 的反射和界面滑行波事件(圖 5), 在 SCS 數(shù)據(jù)中拾取 10 個(gè)縱波反射事件(line4A_r1～line4A_r10), 在 OBS 數(shù)據(jù)中拾取 7 個(gè)縱波反射事件(RL0～RL6)和 5 個(gè)縱波臨界折射(界面滑行波)事件(RF1～RF5)。

(a)在 4 號(hào)地震測(cè)線的 SCS 數(shù)據(jù)中拾取的共 10 個(gè)縱波反射事件(line4A_r1～line4A_r10), 其中 r1 對(duì)應(yīng)海底, r2 和 r3 分別對(duì)應(yīng)淺部高速層的頂、底界面, r5 對(duì)應(yīng)BSR; (b)在 OBS 測(cè)點(diǎn) 33 的數(shù)據(jù)中拾取的廣角地震反射(RL)和臨界折射事件(RF), 其中 RL0, RL1 和RL2 分別對(duì)應(yīng)海底、高速層頂部和 BSR, RF1 對(duì)應(yīng)高速層頂界面, RF2～RF5 為 BSR 之下的深部臨界折射事件, 下同。折合走時(shí)=走時(shí)?偏移距/折合速度, 本文折合速度為 4000 m/s

綜合 SCS 和 OBS 數(shù)據(jù)的拾取事件, 最終將 Slip-stream 海底滑坡沿 4 號(hào)測(cè)線方向的山脊剖面分為 14層(不包含海水)。反演得到的最終射線路徑追蹤結(jié)果見(jiàn)圖 6(a), 擬合得到的走時(shí)與實(shí)際拾取的走時(shí)相符(圖 6(b)), 模型整體的 RMS 走時(shí)殘差為 15ms, 小于拾取誤差(20ms),2=0.901, 達(dá)到預(yù)期要求。

精細(xì)化的縱波速度結(jié)構(gòu)模型實(shí)現(xiàn)對(duì)最小厚度為25m 左右薄層的識(shí)別(圖 7(a)), 其中 L1 與海水直接接觸, 平均速度約為 1.51km/s, 代表非常松散的海底沉積物; L2 的頂部界面在海底之下約 75m, 厚度約為 25m, 縱波速度突然增加至約 2.0km/s, 推測(cè)可能是有高飽和度天然氣水合物富集的砂體; L3 ～L5, 縱波速度從 1.7km/s 緩慢增加至 2.1km/s, 代表無(wú)(或含少量)水合物沉積物的壓實(shí)趨勢(shì); L5 的底部界面對(duì)應(yīng) BSR, 深度約在海底之下 230m; BSR 之下的L6 縱波速度略有削減, 代表游離氣體的存在; 之后, 縱波速度的增加再次向壓實(shí)趨勢(shì)線回歸。

(a)縱波速度結(jié)構(gòu)模型擬合的射線路徑, 4個(gè)圓圈為OBS測(cè)點(diǎn)位置及編號(hào)(下同); (b)縱波速度結(jié)構(gòu)模型的走時(shí)擬合結(jié)果

5 橫波速度結(jié)構(gòu)模型

當(dāng)沉積物孔隙中含游離氣體時(shí), 不但會(huì)影響SCS 反射數(shù)據(jù)和 OBS 垂直分量的成像數(shù)據(jù), 衰減嚴(yán)重時(shí)還會(huì)出現(xiàn)地震反射空白, 大幅度地降低縱波速度, 影響水合物飽和度的計(jì)算。然而, 橫波傳播不受孔隙流體影響, 可以反映含氣區(qū)域(如 BSR 之下的游離氣層)的內(nèi)部結(jié)構(gòu)。并且, 橫波速度在 BSR之上的大幅度增加可能顯示固體天然氣水合物的存在。此外, 由于 OBS 記錄的橫波是由縱波轉(zhuǎn)換而來(lái), 其頻率與縱波相同, 但速度更慢(橫波在純天然氣水合物中的速度約為 1890m/s, 在不含水合物的淺層海底沉積物中的速度為 100～600m/s)[29], 因此波長(zhǎng)更短, 從而可以解析比縱波更精細(xì)的沉積結(jié)構(gòu)信息。

在上述縱波速度結(jié)構(gòu)模型的基礎(chǔ)上, 保持各層頂、底界面的深度節(jié)點(diǎn)和縱波速度節(jié)點(diǎn)不變, 通過(guò)MATLAB 的遺傳算法對(duì)橫波模型的泊松比進(jìn)行同步反演。模型初始的參考泊松比數(shù)值(圖 8(c))通過(guò)15 km 之外的 IODP U1326D 鉆孔的縱、橫波測(cè)井結(jié)果計(jì)算得到:

需要指出的是, 通過(guò)非線性的遺傳算法得到的滿足條件的泊松比組合是全局最優(yōu)解。作為對(duì)比, RAYINVR 的 DMPLSTSQR 程序是線性優(yōu)化, 在初始參數(shù)與模型最終值偏差較大或地質(zhì)構(gòu)造復(fù)雜時(shí), 容易掉進(jìn)局部解的陷阱, 造成反演結(jié)果不佳, 需要重新手動(dòng)調(diào)節(jié)初始參數(shù)再嘗試。圖 9 展示轉(zhuǎn)換橫波的射線追蹤路徑和走時(shí)擬合結(jié)果, 得到的 RMS 走時(shí)殘差為 27ms, 小于橫波事件的拾取誤差(30ms),2=0.903, 達(dá)到預(yù)期要求。

將橫波速度結(jié)構(gòu)模型(圖 7(b))與縱波速度結(jié)構(gòu)模型(圖 7(a))進(jìn)行對(duì)比, 可見(jiàn)淺部的高速異常體 L2的橫波速度增加幅度(14.3%)明顯小于縱波(24.5%), 推測(cè)該層的高飽和度天然氣水合物很可能以固體填隙物的形式富集, 對(duì)沉積物骨架的力學(xué)性能增加不明顯(泊松比與其上接觸海水的 L1 相近)。BSR 之上的 L5 橫波速度相對(duì)于其下游離氣層的增加幅度(14.0%)明顯大于縱波(5.3%), 說(shuō)明縱波速度受到游離氣的干擾, 沒(méi)有顯示出含水合物層的高速特征。圖 8(a)和(b)以 OBS 測(cè)點(diǎn) 34 所在位置的一維深度?速度剖面形式對(duì)比上述縱、橫波速度的變化情況, 圖8(c)的泊松比則顯示淺部高速層(L1)與 BSR 之下游離氣層相比, 泊松比明顯偏高, 指示這兩個(gè)部位為力學(xué)軟弱層。L2 的泊松比大于游離氣層, 可能是Slipstream滑塌體沿著相對(duì)更軟的高速層底界發(fā)育的原因。作為對(duì)比, Orca 滑塌體旁邊的 U1326D 測(cè)井?dāng)?shù)據(jù)顯示淺部的高速異常體出現(xiàn)在海底之下 60～ 90m 處, 而 Slipstream 海底滑坡區(qū)的高速層是在海底之下 75～100m, 這可能是 U1326D 與 Slipstream滑坡區(qū)之間約 15km 的空間距離引起的差異。此外, U1326D 的淺部高速異常體泊松比為低值, 指示力學(xué)強(qiáng)度相對(duì)較高, 并不是一個(gè)軟弱面, 因此其滑脫面出現(xiàn)在泊松比值對(duì)比更強(qiáng)烈的 BSR 位置(應(yīng)力容易在力學(xué)性質(zhì)相差較大的分界部位集中), 顯示與Slipstream 滑塌體不同的水合物控制作用。

圖8 Slipstream海底滑坡區(qū)縱橫波速度結(jié)構(gòu)模型與U1326D測(cè)井?dāng)?shù)據(jù)[30]對(duì)比

6 結(jié)論

1)將 RAYINVR 軟件射線追蹤算法與 MATLAB遺傳反演相結(jié)合, 實(shí)現(xiàn) OBS 橫波速度結(jié)構(gòu)模型各層和層中區(qū)塊泊松比的同步反演, 極大地提高了地質(zhì)構(gòu)造復(fù)雜區(qū)域的橫波速度結(jié)構(gòu)建模效率。

2)建立 Slipstream 海底滑坡體的精細(xì)縱、橫波速度結(jié)構(gòu)模型, 顯示作為水合物穩(wěn)定區(qū)底界地震標(biāo)志的似海底反射層(BSR)位于海底之下 230～250m, 而淺部約 75～100m 的高速異常體很可能是高飽和度天然氣水合物富集的砂體。模型結(jié)果可與 IODP X311 航次 U1326D 測(cè)井?dāng)?shù)據(jù)進(jìn)行對(duì)比, 驗(yàn)證了同步反演速度結(jié)構(gòu)建模的有效性。

3)淺部高速異常體的深度與滑坡面大致重合, 顯示水合物造成的力學(xué)軟弱面對(duì)滑坡構(gòu)造的控制作用, 為后續(xù)建立三維力學(xué)模型提供了關(guān)鍵證據(jù)。

(a) OBS 橫波速度結(jié)構(gòu)模型擬合的射線路徑, S1～S8代表轉(zhuǎn)換橫波反射事件; (b) OBS 橫波速度結(jié)構(gòu)模型的走時(shí)擬合結(jié)果

致謝 衷心感謝加拿大維多利亞大學(xué) George Spence 教授提供 OBS, SCS 和多波束等數(shù)據(jù), 特別感謝 IHS 全球公司提供地球科學(xué)綜合解釋軟件包“The Kingdom Suite”。

[1] Makogon Y F. Natural gas hydrates — a promising source of energy. Journal of Natural Gas Science and Engineering, 2010, 2(1): 49?59

[2] Helgerud M B, Dvorkin J, Nur A, et al. Elastic-wave velocity in marine sediments with gas hydrates: effective medium modeling. Geophysical Research Letters, 1999, 26(13): 2021?2024

[3] Yang J, Davies R J. Gravity-driven faults: migration pathways for recycling gas after the dissociation of marine gas hydrates. Marine Geology, 2013, 336: 1?9

[4] Davies R J, Clarke A L. Methane recycling between hydrate and critically pressured stratigraphic traps, offshore Mauritania. Geology, 2010, 38(11): 963?966

[5] Dash R, Spence G. P-wave and S-wave velocity struc-ture of northern Cascadia margin gas hydrates. Geo-physical Journal International, 2011, 187(3): 1363? 1377

[6] 張光學(xué), 徐華寧, 劉學(xué)偉, 等. 三維地震與OBS聯(lián)合勘探揭示的神狐海域含水合物地層聲波速度特征. 地球物理學(xué)報(bào), 2014, 57(4): 1169?1176

[7] 劉伊克, 朱偉林, 米立軍, 等. 南海深水多次波成像. 中國(guó)科學(xué): 地球科學(xué), 2015, 45(2): 152?160

[8] Davis E E, Hyndman R D. Accretion and recent defor-mation of sediments along the northern Cascadia sub-duction zone. Geological Society of America Bulle-tin, 1989, 101(11): 1465?1480

[9] Hyndman R D. The Lithoprobe corridor across the Vancouver Island continental margin: the structural and tectonic consequences of subduction. Canadian Journal of Earth Sciences, 1995, 32(10): 1777?1802

[10] Mackay M E, Jarrard R D, Westbrook G K, et al. Origin of bottom-simulating reflectors: geophysical evidence from the Cascadia accretionary prism. Geo-logy, 1994, 22(5): 459

[11] Pohlman J W, Kaneko M, Heuer V B, et al. Methane sources and production in the northern Cascadia mar-gin gas hydrate system. Earth and Planetary Science Letters, 2009, 287(3/4): 504?512

[12] Riedel M, Tréhu A, Spence G. Characterizing the thermal regime of cold vents at the northern Cascadia margin from bottom-simulating reflector distributions, heat-probe measurements and borehole temperature data. Marine Geophysical Researches, 2010, 31(1/2): 1?16

[13] Hyndman R D, Riedel M, Spence G D. Hydrate studies of northern Cascadia margin off Vancouver Island: a reference source // 6th International Con-ference on Gas Hydrates (ICGH 2008). Vancouver, 2008, doi: 10.14288/1.0040999

[14] Westbrook G K, Carson B, Musgrave R J. Proceeding of Ocean Drilling Program, volumn 146 part 1: initial reports [R]. Texas: Ocean Drilling Program, Texas A&M University, 1994

[15] Torres M E, Tréhu A M, Cespedes N, et al. Methane hydrate formation in turbidite sediments of northern Cascadia, IODP Expedition 311. Earth and Planetary Science Letters, 2008, 271: 170?180

[16] Riedel M, Collett T S, Malone M J. Gas hydrate drilling transect across northern Cascadia margin — IODP Expedition 311. Geological Society London Special Publications, 2009, 319(1): 11?19

[17] Expedition 311 Scientists. Cascadia margin gas hyd-rates [R]. Washington, DC: Integrated Ocean Drilling Program Management International, 2005

[18] Lopez C, Spence G, Hyndman R, et al. Frontal ridge slope failure at the northern Cascadia margin: margin-normal fault and gas hydrate control. Geology, 2010, 38(11): 967?970

[19] Yelisetti S, Spence G D, Riedel M. Role of gas hydrates in slope failure on frontal ridge of northern Cascadia margin. Geophysical Journal International, 2014, 199(1): 441?458

[20] Zelt C A, Smith R B. Seismic traveltime inversion for 2-D crustal velocity structure. Geophysical Journal International, 1992, 108(1): 16?34

[21] Catchings R D, Mooney W D. Crustal structure of the Columbia Plateau: evidence for continental rifting. Journal of Geophysical Research: Solid Earth, 1988, 93(B1): 459?474

[22] Boland A V, Ellis R M. Velocity structure of the Kapuskasing Uplift, northern Ontario, from seismic refraction studies. Journal of Geophysical Research: Solid Earth, 1989, 94(B6): 7189?7204

[23] Zelt C A, Ellis R M. Practical and efficient ray tracing in two-dimensional media for rapid traveltime and amplitude forward modeling. Canadian Journal of Exploration Geophysics, 1988, 24(1): 16?31

[24] Hyndman R D, Spence G D. A seismic study of methane hydrate marine bottom simulating reflectors. Journal of Geophysical Research: Solid Earth, 1992, 97(B5): 6683?6698

[25] Chen M. Northern Cascadia marine gas hydrate: con-straints from resistivity, velocity, and AVO [D]. Vic-toria: University of Victoria, 2006

[26] Lopez C. Seismic velocity structure associated with gas hydrate at the frontal ridge of Northern Cascadia Margin [D]. Victoria: University of Victoria, 2007

[27] 隆松伯, 何濤, 梁前勇, 等. 基于ANSYS Work-bench的高精度自動(dòng)化三維地質(zhì)建模方法——以天然氣水合物相關(guān)的Slipstream海底滑坡為例. 北京大學(xué)學(xué)報(bào)(自然科學(xué)版), 2018, 54(5): 994?1002

[28] Hyndman R D, Davis E E. A mechanism for the formation of methane hydrate and seafloor bottom-simulating reflectors by vertical fluid expulsion. Journal of Geophysical Research: Solid Earth, 1992, 97(B5): 7025?7041

[29] Waite W F, Helgerud M B, Nur A, et al. Labora- tory measurements of compressional and shear wave speeds through methane hydrate. Annals of the New York Academy of Sciences, 1999, 912(1): 1003?1010

[30] Riedel M, Collett T S, Malone M J. Expedition 311 summary [R]. Washington, DC: Integrated Ocean Drilling Program Management International, 2006

Model Inversion of Velocity Structure for Slipstream Submarine Slide Related to Gas Hydrate

LAN Kun1,2, ZHU He1,2, HE Tao1,2,?, LIANG Qianyong3, WU Xuemin3, DONG Yifei3, ZHANG Yi4

1. Key Laboratory of Orogenic Belts and Crustal Evolution (MOE), School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Beijing International Center for Gas Hydrate, Peking University, Beijing 100871; 3. Guangzhou Marine Geological Survey, China Geological Survey, Guangzhou 510760; 4. Chinese Academy of Geological Sciences, Beijing 100037; ? Corresponding author, E-mail: taohe@pku.edu.cn

The RAYINVR software, which is widely used in academia for 2-D velocity structure model from four-component OBS (ocean bottom seismometer) data, is unable to invert for the converted S-wave automatically, resulting in the low-efficiency of modeling process. Using MATLAB’s genetic algorithm, the RAYINVR software is improved and able to automatically and synchronously invert for Poisson’s ratios of each layer with all sub-blocks for the S-wave velocity structure model, and thus can provide Young’s modulus, Poisson’s ratio and other important mechanical information for gas hydrate survey. This method is applied to process the OBS data collected at the Slipstream submarine slide, and a fine P- and S-wave velocity structure model is obtained, which is comparable to the logging data of nearby borehole U1326. Therefore, the validity of the auto-synchronous inversion method is verified for the S-wave velocity structure modeling. The optimal velocity model reveals two structural interfaces with large Poisson’s ratio contrast. One is BSR (bottom simulating reflector) at 230±10 mbsf (meter beneath sea floor), which represents the bottom boundary of the gas hydrate stability zone, and the other is the basal boundary of a shallow abnormal high-speed body (possibly a sand body enriched with high saturation gas hydrate) at 75?100 mbsf. The latter agrees roughly with the glide plane of Slipstream submarine slide, indicating that the hydrate is related to the formation of submarine landslide.

synchronous inversion; S-wave velocity structure; Slipstream submarine slide; gas hydrate; ray tracing

10.13209/j.0479-8023.2021.005

國(guó)家自然科學(xué)基金(41676032)和中國(guó)地質(zhì)調(diào)查局國(guó)家天然氣水合物專(zhuān)項(xiàng)基金(DD20190234, DD20190218, DD20189320, DD20160217, HD-JJHT-20)資助

2020–03–10;

2020–05–15