考慮水層鳴震的海底多分量數(shù)據(jù)轉(zhuǎn)換波走時(shí)反演方法

劉學(xué)義, 程玖兵, 王騰飛

同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200092

0 引言

在現(xiàn)代地震觀測(cè)和高性能計(jì)算技術(shù)的推動(dòng)下，多分量地震方法在油氣資源勘探領(lǐng)域越來(lái)越受到重視.它更全面地記錄了地震縱波(P)與橫波(S)的矢量振動(dòng)信息，用于改善“氣云區(qū)”和復(fù)雜地質(zhì)體(如鹽丘、火成巖等)地震成像，也有助于區(qū)分巖性，刻畫含油氣儲(chǔ)層及其內(nèi)部非均質(zhì)性，檢測(cè)定向發(fā)育的裂縫系統(tǒng)以及分析應(yīng)力分布特征等(Stewart et al.,2002; 趙邦六,2008).

在多分量數(shù)據(jù)處理流程中，PS轉(zhuǎn)換波偏移成像處于核心地位，其成敗很大程度取決于地震記錄P/S分離的可靠性，以及P波和S波偏移速度模型的合理性(蘆俊等,2018).通常，在分離出PP波和PS波地震記錄之后，采用傳統(tǒng)的偏移速度分析或走時(shí)層析技術(shù)就能得到P波速度模型，而估計(jì)S波速度則要遵循PS波入射-反射路徑非對(duì)稱性，采用相適應(yīng)的速度分析或走時(shí)反演方法(Broto et al.,2003;杜啟振等,2011;Wang et al.,2019).隨著基于波動(dòng)方程的逆時(shí)偏移(RTM)和全波形反演(FWI)等高精度成像與建模技術(shù)的推廣應(yīng)用，面向多分量地震數(shù)據(jù)的RTM與FWI方法也逐漸受到關(guān)注(Sears et al.,2008;Ren and Liu,2016).當(dāng)?shù)卣鹩涗浫鄙倏煽康牡皖l成分和大偏移距信號(hào)時(shí)，F(xiàn)WI很難合理地重構(gòu)中深層的縱、橫波宏觀速度模型，影響反射PP波與PS波成像.

這些年來(lái)，強(qiáng)化FWI的反射層析功能或發(fā)展針對(duì)反射波信號(hào)的反演方法成為了勘探地球物理界的研究熱點(diǎn)，其中有代表性且與本文相關(guān)的工作包括Xu 等(2012)提出的反射波形反演(RWI)以及由它退化而成的波動(dòng)方程反射走時(shí)反演(RTI)方法 (Ma and Hale,2013).雖然理論上這些方法經(jīng)過(guò)拓展也能處理多分量地震數(shù)據(jù)，但都會(huì)面臨P-S波串?dāng)_、多參數(shù)(如縱-橫波速度)耦合造成的強(qiáng)烈非線性以及由此引起的收斂性問題.為此，除了基于海森矩陣的二階優(yōu)化方法(Epanomeritakis et al.,2008)或子空間方法(Baumstein,2014)，近年一些學(xué)者提出了P/S模式解耦預(yù)條件的彈性波FWI(Wang and Cheng,2017)、彈性波RTI/RWI(Wang et al.,2018)以及彈性波動(dòng)方程偏移速度分析方法(Wang et al.,2019).在此基礎(chǔ)上，Xu 等(2019)基于P/S模式解耦的數(shù)據(jù)或梯度預(yù)條件，建立了串聯(lián)標(biāo)量聲波方程反射走時(shí)反演、彈性波方程RWI和FWI的多級(jí)速度建模方法，有效地利用了多分量地震數(shù)據(jù)的走時(shí)與波形信息，構(gòu)建出高精度的縱、橫波速度模型.值得注意的是，上述方法處理實(shí)際海底多分量地震資料時(shí)要克服一系列挑戰(zhàn).例如，海面的自由界面效應(yīng)和海底的上-下行波耦合會(huì)極大地增加海底多分量地震信號(hào)的復(fù)雜性.

海底上-下行P/S波分離對(duì)多分量地震數(shù)據(jù)偏移速度分析是非常必要的.當(dāng)前，海底多分量記錄的波場(chǎng)分離方法分為兩大類：一類是基于聲波方程的分離方法，如PZ疊加(Dragoset and Barr,1994;Osen et al.,1999)，它可以很好地壓制鬼波對(duì)水聽器P波記錄的干擾；另一類是基于彈性波方程的分離方法(Schalkwijk et al.,2003;Muijs et al.,2007)，它同時(shí)完成上-下行波分離和P/S解耦，壓制檢波器端的自由表面多次波，但分離的上、下行P波和上行S波記錄仍包含源端自由表面多次波.比如，針對(duì)由繩系海底地震儀(OBS)在我國(guó)東海QY探區(qū)采集的多分量地震資料，劉學(xué)義等(2021)根據(jù)噪聲特點(diǎn)對(duì)傳統(tǒng)海底數(shù)據(jù)P/S模式分離方法與流程進(jìn)行了調(diào)整，發(fā)現(xiàn)分離后的P波和S波記錄仍受源端多次波干擾，給后續(xù)處理尤其是速度分析和偏移成像帶來(lái)極大的困擾.

所以，在完成上-下行P/S分離之后，海底多分量地震數(shù)據(jù)處理還需克服源端自由表面多次波的影響(Brown and Yan,1999).對(duì)于P波記錄，通行的做法是利用預(yù)測(cè)濾波(Wiggins,1988)或基于波動(dòng)方程的自由表面多次波預(yù)測(cè)-相減法(SRME: Verschuur et al.,1992)進(jìn)一步濾除源端多次波，然后進(jìn)行速度分析與偏移成像.近年來(lái)，人們也提出了一些視多次波為有效信號(hào)的偏移成像方法(劉伊克等,2018;Liu et al.,2020).不過(guò)，所有這些方法都僅關(guān)注反射PP波.在PS波速度分析與偏移成像中考慮源端P波水層鳴震，需克服海底介質(zhì)的復(fù)雜影響和PS波傳播路徑的非對(duì)稱性，挑戰(zhàn)很大，但相關(guān)研究比較少(Dawson et al.,2010).

本文基于海洋聲學(xué)與地球物理文獻(xiàn)資料，分析地震波在典型海底界面的能量分配規(guī)律，調(diào)查軟海底情況下源端多次波的類型，明確源端水層鳴震對(duì)P/S分離之后獲得的上行S波數(shù)據(jù)(即本文所說(shuō)的轉(zhuǎn)換波數(shù)據(jù))的主導(dǎo)性干擾(見圖1至圖3).若已知縱波速度模型和PP波成像結(jié)果，針對(duì)橫波速度建模問題，我們提出一種考慮源端水層鳴震的波動(dòng)方程轉(zhuǎn)換波反射走時(shí)反演方法,以擺脫常規(guī)方法因鳴震干擾帶來(lái)的數(shù)據(jù)擬合困難和收斂性問題.

1 上行S波數(shù)據(jù)多次波干擾分析

海底多分量地震記錄受到海面和海底這兩個(gè)特殊界面的雙重影響.大量海洋聲學(xué)調(diào)查和地質(zhì)、地球物理資料表明，除了基巖或老地層直接出露的局部區(qū)域，海底大都覆蓋一層未固結(jié)軟性沉積物，具有速度隨深度逐漸增加、泊松比或縱-橫波速度比較大的特性(Carlson et al.,1984; Berge et al.,1991).海底介質(zhì)的縱波速度一般低于兩千米每秒,橫波速度在幾十到幾百米每秒之間(Hamilton,1979;Sidler and Holliger,2010).海底表層沉積物孔隙度較大，隨著深度增加介質(zhì)孔隙度逐漸變小，密度逐步增大，以孔隙度為35%的砂質(zhì)沉積為例，密度在2100 kg·m-3左右(Buckingham,2005).針對(duì)較普遍的軟質(zhì)海底環(huán)境，下面分析地震波在海底界面的能量分配規(guī)律，推斷在海底觀測(cè)的上行S波中的源端多次波的主導(dǎo)成分.在給定海水P波速度為1.5 km·s-1、密度為1.02 g·cm-3的情況下，依據(jù)海底介質(zhì)縱、橫波速度和密度取值范圍，各四次均勻采樣構(gòu)成64種參數(shù)組合，然后用Zoeppritz方程的線性近似(Aki and Richards,2002)計(jì)算海底界面處P波、轉(zhuǎn)換S波的透射與反射系數(shù)，最后統(tǒng)計(jì)并繪制其臨界角之前的平均值(如圖1所示，四個(gè)大圖板代表四種橫波速度，每個(gè)圖板中的四列代表四種縱波速度，每列四組代表四種密度取值).圖2展示了其中兩種典型軟海底界面的透射/反射系數(shù).不難發(fā)現(xiàn)，盡管大部分P波能量經(jīng)海底透射到地層中，但其中的透射轉(zhuǎn)換S波并不強(qiáng)；即便對(duì)于橫波速度較大、泊松比較小的海底地層，海底界面的PP波反射也強(qiáng)于PS透射波.所以在軟海底情況下，在海底界面發(fā)生的模式轉(zhuǎn)換對(duì)最終分離出的S波記錄貢獻(xiàn)可以忽略(如圖2).這種能量分配特點(diǎn)表明，軟海底多分量地震記錄的源端多次波干擾主要是一些P波水層鳴震或自由表面多次波.這與前人一些研究相吻合(Caldwell,1999; Suarez,2000).因此，在海底界面實(shí)施上-下行P/S分離之后，S波數(shù)據(jù)中除了源自于海底之下一些阻抗界面的反射PS波，也包含一些源端多次波，其中由海面與海底各反射一次的PPS波(圖3粗實(shí)線)路徑最短、能量最強(qiáng).而當(dāng)海底以下不存在鹽丘等強(qiáng)波阻抗界面時(shí)，地下界面反射波場(chǎng)相關(guān)的PPS波能量比較弱(圖3細(xì)虛線),本文不予考慮.后文數(shù)值試驗(yàn)部分基于一階速度-應(yīng)力方程模擬的海底多分量地震記錄也與上述分析相吻合.

圖1 海底界面反射/透射系數(shù)平均值統(tǒng)計(jì)圖板 縱軸代表透射/反射系數(shù)平均值；橫軸為縱/橫波速度四次采樣；每一組含4次密度采樣， 取值依次為1.5 g·cm-3、1.7 g·cm-3、1.9 g·cm-3和2.1 g·cm-3.Fig.1 Average statistics of transmission/reflection coefficient at seabed interface The vertical axis represents the average of transmission/reflection coefficient; the horizontal axis represents four samples of P-/S- wave velocities, and each group has 4 samples (1.5 g·cm-3、1.7 g·cm-3、1.9 g·cm-3 and 2.1 g·cm-3) of density.

圖2 兩種典型軟海底界面的反射/透射系數(shù)曲線 (a) vP=1.7 km·s-1,vS=0.1 km·s-1,ρ=1.9 g·cm-3; (b) vP=1.7 km·s-1,vS=0.5 km·s-1,ρ=1.9 g·cm-3.Fig.2 Transmission/reflection coefficient curves of two typical soft seabed interfaces

圖3 海底上行S波傳播路徑示意圖 黑線表示P波傳播路徑，紅線表示S波傳播路徑.Fig.3 Schematic diagram of propagation paths of up-going S-waves The black lines are the paths of P-waves, and the red lines are the paths of S-waves.

2 考慮源端水層鳴震的轉(zhuǎn)換波走時(shí)反演方法

國(guó)內(nèi)外迄今采集了很多海底多分量地震數(shù)據(jù).為了利用縱波數(shù)據(jù)，最常用的前處理手段是聯(lián)合水檢和陸檢信號(hào)壓制檢波器端海水鳴震(李維新等,2018).為了充分利用記錄到的縱波與橫波信號(hào)，通常要在海底界面上實(shí)施上-下行P/S 波分離，壓制檢波器端的水層鳴震或海面多次波，殘余的源端多次波還需要其他多次波消去法進(jìn)行壓制(Edme and Sing,2009).近年來(lái)，陸續(xù)出現(xiàn)了一些可同時(shí)處理一次反射波、表面多次波(甚至層間多次波)的偏移成像和速度分析方法(Lu et al.,2015;李志娜和李振春,2016;葉月明等,2019;Liu et al.,2020),但這些研究仍局限于縱波數(shù)據(jù).

根據(jù)前文分析，在上-下行P/S 波分離得到的上行S波數(shù)據(jù)中，除了人們需要的PS波，還包括PPS波等源端水層鳴震和其他噪聲.我們?cè)谇捌谘芯恐邪l(fā)現(xiàn)，這些源端多次波會(huì)給基于轉(zhuǎn)換波數(shù)據(jù)(上行S波數(shù)據(jù))的橫波速度分析帶來(lái)極大的困擾.為此，下文在波動(dòng)方程反射走時(shí)反演理論框架下，首先給出反射PS波及含源端水層鳴震路徑的S波(即PPS波)的Born模擬方程，然后建立轉(zhuǎn)換反射波走時(shí)擬合目標(biāo)泛函，分析PS波和PPS波這兩種成分對(duì)泛函梯度或反射波敏感核的影響，進(jìn)而提出合適的梯度預(yù)條件方法.這里假設(shè)前期處理已獲得可靠的縱波速度模型、海底界面位置以及PP波成像結(jié)果，聚焦宏觀橫波速度模型構(gòu)建問題.

2.1 正問題

考慮到走時(shí)反演對(duì)彈性效應(yīng)不敏感，故采用聲波方程來(lái)相對(duì)經(jīng)濟(jì)地描述縱波和轉(zhuǎn)換橫波的傳播過(guò)程.常密度標(biāo)量聲波方程在頻率域可寫成：

L(mp(x),ω)u(x,xs,ω)=-F(ω)δ(x-xs),

(1)

L(mp0(x),ω)up0(x,xs,ω)=-F(ω)δ(x-xs),

(2)

L(mp0(x),ω)up1(x,xs,ω)=Kp(x,ω)up0(x,xs,ω),

(3)

式中Kp(x,ω)=ω2mp1(x)是高波數(shù)縱波慢度擾動(dòng)形成的二次虛源，一階散射場(chǎng)up1通常稱為PP波.對(duì)于一次轉(zhuǎn)換/反射PS波，替換(3)式中的模型參數(shù)，也可用一階Born模擬合成運(yùn)動(dòng)學(xué)可靠的數(shù)據(jù)，表示為

L(ms0(x),ω)us0(x,xs,ω)=Ks(x,ω)up0(x,xs,ω),

(4)

式中Ks(x,ω)=ω2ms1(x)，ms0和ms1分別代表橫波慢度平方模型的背景和擾動(dòng)分量，us0表示PS波，其傳播路徑如圖3所示.

對(duì)于源端水層鳴震中占主導(dǎo)的PPS波，可認(rèn)為是經(jīng)海底反射到達(dá)海面的一次反射波作為新的震源，經(jīng)水層傳播到地層中在一些阻抗界面發(fā)生P-S模式轉(zhuǎn)換與反射，最后上傳到海底被檢波器接收.假定已通過(guò)海洋聲學(xué)調(diào)查或地球物理方法獲得海底界面信息(Yang et al.,2015;Perdomo et al.,2018)，則可由公式(2)和(3)模擬僅僅由海底界面反射到達(dá)海面的PP波；然后借助下面兩個(gè)方程模擬出運(yùn)動(dòng)學(xué)可靠的PPS波，即

L(mp0(x),ω)up2(x,xs,ω)=-up1(x,xs,ω)δ(x-x0),

(5)

L(ms0(x),ω)us1(x,xs,ω)=Ks(x,ω)up2(x,xs,ω),

(6)

(7)

其中加權(quán)系數(shù)w0和w1用于調(diào)節(jié)兩種波場(chǎng)成分的相對(duì)幅值，其大小主要與兩種轉(zhuǎn)換波傳播路程引起的幾何擴(kuò)散效應(yīng)差異(在淺水區(qū)可忽略不計(jì))以及PPS波發(fā)生在海底、海面的兩次反射有關(guān).加權(quán)系數(shù)的設(shè)置方法參見附錄A.

2.2 反射波走時(shí)擬合目標(biāo)泛函及其梯度計(jì)算

波動(dòng)方程反射走時(shí)反演的關(guān)鍵在于定量估計(jì)模擬反射波與觀測(cè)反射波的走時(shí)擬合誤差.主流方法包括動(dòng)態(tài)圖像匹配(DIW)(Hale,2013)和局部互相關(guān)(Van Leeuwen & Mulder,2010).Xu 等(2019)將Ma 和Hale(2013)提出的基于DIW的反射走時(shí)反演方法拓展到了多分量地震數(shù)據(jù)，但他們沒有考慮水層鳴震等多次波干擾.在將分離出來(lái)的上行S波記錄作為觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)進(jìn)行走時(shí)匹配時(shí)，必然受到一些較強(qiáng)的源端水層鳴震(如PPS波)的影響(可見數(shù)值試驗(yàn)部分的例證).為此，我們認(rèn)為既要在模擬轉(zhuǎn)換反射波時(shí)同時(shí)考慮PS波和PPS波，又要在反演過(guò)程中避免二者相互串?dāng)_對(duì)模型更新方向的畸變.下面介紹本文的方法和處理步驟.

τ=argminlD(l(xr,xs,t)),

(8)

式中l(wèi)(xr,xs,t)為時(shí)空變化的時(shí)移量，D是l的函數(shù)，用來(lái)描述觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的L2范數(shù)殘差，即：

(9)

于是將轉(zhuǎn)換波走時(shí)反演的目標(biāo)函數(shù)定義為

(10)

在局部?jī)?yōu)化算法基礎(chǔ)上，逐步更新橫波速度模型使目標(biāo)函數(shù)最小化.

借鑒Ma 和Hale(2013)的PP反射波走時(shí)反演方法，可推導(dǎo)出轉(zhuǎn)換波(含PS波和PPS波)反射走時(shí)反演的泛函梯度表達(dá)式：

(11)

(12)

注意，公式(11)表示泛函梯度是由正向擾動(dòng)場(chǎng)與伴隨入射場(chǎng)零延遲互相關(guān)獲得，代表地下反射界面與海底檢波器之間上行波路徑對(duì)橫波速度的敏感性，這與PS、PPS轉(zhuǎn)換波運(yùn)動(dòng)學(xué)特征一致.

假設(shè)背景橫波速度擾動(dòng)較小，真實(shí)反射數(shù)據(jù)可由時(shí)移后的模擬反射數(shù)據(jù)替代(Luo and Schuster,1991)，即：

(14)

則伴隨源可簡(jiǎn)化為DIW估計(jì)的走時(shí)殘差τ(xr,xs,t)對(duì)模擬記錄偏導(dǎo)數(shù)波場(chǎng)的加權(quán)，進(jìn)而把公式(12)改寫為

2.3 轉(zhuǎn)換反射波走時(shí)敏感核分析

這部分基于水平層狀介質(zhì)模型(圖4)合成OBC數(shù)據(jù)進(jìn)行泛函梯度或敏感核分析.震源位于(360 m,15 m)處，采用主頻為20 Hz的雷克子波，在海底(深度200 m)以5 m的間隔接收模擬的多分量地震信號(hào)(圖4a).圖5顯示按本文Born模擬方法合成的含PS波和PPS波的轉(zhuǎn)換波記錄.

圖4 水平層狀介質(zhì)模型 (a) 縱波速度； (b) 橫波速度.五角星表示炮點(diǎn)位置，倒三角表示檢波器位置.Fig.4 Horizontal layered model (a) P-wave velocity； (b) S-wave velocity. The star is for source, and the triangles are for receivers.

圖5 共炮點(diǎn)轉(zhuǎn)換波數(shù)據(jù) 右側(cè)子圖中黑線表示P波路徑，紅線表示S波路徑.Fig.5 Converted waves of common-shot gather In the rightsubgraphs, black lines represent P-waves path, and red lines represent S-waves path.

圖6顯示按前文公式計(jì)算的單炮檢對(duì)(炮點(diǎn)、檢波器位置已在圖4b中標(biāo)識(shí))對(duì)應(yīng)的轉(zhuǎn)換波走時(shí)關(guān)于橫波速度的敏感核，包含低波數(shù)的轉(zhuǎn)換反射波路徑和高波數(shù)的偏移脈沖響應(yīng).圖6a最為復(fù)雜，這是由PS波、PPS波的偏移脈沖響應(yīng)以及轉(zhuǎn)換反射波路徑相互混疊造成.對(duì)于模擬數(shù)據(jù)，可在梯度計(jì)算過(guò)程中對(duì)參與互相關(guān)的伴隨入射波場(chǎng)和正向擾動(dòng)波場(chǎng)進(jìn)行PS波與PPS波分離，進(jìn)而獲得四種分解后的敏感核(圖6b-e).可見，通過(guò)對(duì)轉(zhuǎn)換波記錄進(jìn)行波模式分離，可以摒棄串模式互相關(guān)產(chǎn)生的高波數(shù)干擾(圖6d和6e)，有助于壓制對(duì)低波數(shù)轉(zhuǎn)換反射波路徑的污染(見圖6b和6c).注意，圖6b和6c中的低波數(shù)敏感核形狀雖然相似，但傾斜度有差異(后者稍陡一些)，這與該炮檢對(duì)PS和PPS波在檢波器端反射路徑的傾角差異相吻合.圖中箭頭指示的殘余高波數(shù)干擾是單一模式轉(zhuǎn)換波記錄的偏移脈沖響應(yīng)，對(duì)更新宏觀速度模型沒有幫助，經(jīng)多炮檢對(duì)敏感核疊加后明顯減弱，也可在波場(chǎng)互相關(guān)時(shí)濾除小角度相遇的正向與伴隨波場(chǎng)進(jìn)一步加以壓制.

2.4 梯度預(yù)條件處理

前文敏感核數(shù)值分析表明，對(duì)參與泛函梯度計(jì)算的正向和伴隨波場(chǎng)進(jìn)行PS波與PPS波模式分離，可以明顯改善泛函梯度的合理性.基于聲波方程Born模擬算子可以單獨(dú)模擬正向傳播的這兩種震相，但由于直接在觀測(cè)記錄中區(qū)分PS與PPS很困難，因此獲得相應(yīng)的伴隨波場(chǎng)并不容易.

由于公式(14)引入了對(duì)觀測(cè)記錄的運(yùn)動(dòng)學(xué)替代處理，因此只要在Born模擬時(shí)區(qū)分這兩種震相，就可以按下面方程分別實(shí)現(xiàn)PS和PPS伴隨波場(chǎng)的逆時(shí)傳播，即

(16)

(17)

其中αi取值范圍為0～1，且滿足α0+α1=1.當(dāng)α0=1且α1=0時(shí)，表示只用PS波走時(shí)信息；當(dāng)α0=0且α1=1時(shí)，表示只用PPS波走時(shí)信息.考慮到波場(chǎng)中PS波能量占主導(dǎo)，α0取值通常不小于α1.上式表示將類似圖6b和圖6c中的單震相轉(zhuǎn)換反射波敏感核加權(quán)疊加，獲得合理的橫波速度更新方向.

圖6 單炮檢對(duì)情況下橫波速度對(duì)應(yīng)的泛函梯度及其分解 (a) 原始泛函梯度； (b) PS波泛函梯度； (c) PPS波泛函梯度； (d) 正向擾動(dòng)PPS波與伴隨入射PS波互相關(guān)產(chǎn)生的梯度假象； (e) 正向擾動(dòng)PS波與伴隨入射PPS波互相關(guān)產(chǎn)生的梯度假象.Fig.6 The gradient with respect to S-wave velocity and its decomposition for a single pair of shot and receiver (a) Original functional gradient； (b) Functional gradient of PS waves； (c) Functional gradient of PPS waves； (d) Gradient artifacts generated by the cross-correlation between forward perturbed PPS waves and adjoint incident PS waves； (e) Gradient artifacts generated by the cross-correlation between forward perturbed PS waves and adjoint incident PPS waves.

一旦確定更新方向，則按局部尋優(yōu)方法進(jìn)行背景橫波速度模型的更新：

(18)

3 數(shù)值實(shí)驗(yàn)

轉(zhuǎn)換波走時(shí)反演流程圖如圖7所示.本文聚焦于利用轉(zhuǎn)換波走時(shí)信息進(jìn)行宏觀橫波速度建模，所以假設(shè)前期處理已獲得可靠的縱波速度模型、海底界面位置以及PP波成像結(jié)果，并且已經(jīng)完成海底多分量數(shù)據(jù)的P/S分離工作.

圖7 轉(zhuǎn)換波走時(shí)反演流程圖Fig.7 Flowchart of converted wave travel time inversion

3.1 水平層狀模型

首先基于水平層狀介質(zhì)模型檢驗(yàn)方法有效性.該模型水深為200 m，海底界面以下縱、橫波速度隨深度線性增加，在海底附近泊松比較大(圖8).海水表面采用自由表面邊界條件，基于彈性波速度-應(yīng)力方程有限差分算法合成多分量觀測(cè)數(shù)據(jù)，震源為主頻25 Hz的雷克子波，炮間距為50 m，共40炮均勻分布于水深5 m處，檢波器置于海底(如圖8a中三角形所示)，間隔為5 m，最大炮檢距為700 m，記錄時(shí)長(zhǎng)1.15 s.圖9a至9c顯示了炮點(diǎn)位于(1250 m，5 m)處的三分量記錄.基于彈性波理論進(jìn)行上-下行P/S波分離(劉學(xué)義等,2021)，獲得的海底上行S波數(shù)據(jù)如圖9d所示.可見，它仍然受到源端自由表面多次波的強(qiáng)烈干擾，也存在一些較弱的層間多次波噪聲.總體看來(lái)，經(jīng)海底、海面各反射一次的水層鳴震形成了在源端多次波中占主導(dǎo)的PPS波信號(hào).

圖8 水平層狀介質(zhì)模型 (a) 縱波速度； (b) 橫波速度； (c) 密度.Fig.8 Horizontal layered medium model (a) P-wave velocity； (b) S-wave velocity； (c) Density.

圖9 數(shù)值模擬的三分量記錄和分離后的轉(zhuǎn)換波數(shù)據(jù) (a) 聲壓分量； (b) vz分量； (c) vx分量； (d) 海底分離出的上行S波.Fig.9 Synthetic 3C common shot gathers and the decomposed up-going S-waves (a) Pressure； (b) Vertical particle velocity； (c) Horizontal particle velocity； (d) Decomposed up-going S-waves.

為了反演橫波速度，假定已知準(zhǔn)確的縱波背景速度、海底界面模型和反射PP波成像結(jié)果.在這里，縱波背景速度由真實(shí)速度垂向平滑獲得(圖10a).考慮走時(shí)反演對(duì)波場(chǎng)模擬振幅精度要求不高，而且可以由權(quán)系數(shù)調(diào)節(jié)合成轉(zhuǎn)換波記錄中的PS波和PPS波相對(duì)振幅關(guān)系，因此要求海底界面模型深度準(zhǔn)確、能大致反映從水層到海底的速度擾動(dòng)即可，文中不再展示.PP成像結(jié)果由反射PP波記錄逆時(shí)偏移獲得(圖10b).初始橫波速度模型是真實(shí)層速度垂向大尺度平滑后并減去一個(gè)常量得到(圖11a)，密度設(shè)為常數(shù).公式(17)中兩種震相權(quán)重相等，即α0=α1=0.5.接下來(lái)進(jìn)行兩組反演試驗(yàn)：第一組的正演引擎只能模擬PS波，則出現(xiàn)模擬記錄與含有較強(qiáng)PPS波信號(hào)的觀測(cè)記錄很難匹配的問題，獲取不到比較合理的走時(shí)殘差，導(dǎo)致反演結(jié)果不合理(圖11b和圖12藍(lán)色實(shí)線).第二組的正演引擎同時(shí)模擬出PS波和PPS波信號(hào)，轉(zhuǎn)換波走時(shí)殘差估計(jì)變得相對(duì)容易，最終得到了合理的反演結(jié)果(圖11c和圖12紅色實(shí)線).注意，模型底部速度無(wú)法更新是因?yàn)闆]有反射波場(chǎng)照明這一層.

圖10 先驗(yàn)?zāi)Ｐ?(a) 背景縱波速度； (b) 反射PP波逆時(shí)偏移結(jié)果.Fig.10 The prior model (a) Background P-wave velocity； (b) Reverse time migration (RTM) of the reflected PP waves.

圖11 含源端水層鳴震的轉(zhuǎn)換波數(shù)據(jù)橫波速度反演結(jié)果 (a) 初始背景橫波速度； (b) 正演引擎無(wú)法模擬PPS波時(shí)的反演結(jié)果； (c) 本文方法走時(shí)反演結(jié)果.Fig.11 The inversion S-wave velocity with converted wave recordings including water-layer source-side reverberations (a) The initial background S-wave velocity； (b) The result of travel-time inversion with a forward modeling engine that can not simulates PS waves； (c) The proposed method.

圖12 偽井橫波速度曲線Fig.12 Pseudo-well curves of the S-wave velocity models

圖13把兩組實(shí)驗(yàn)由初始和反演的橫波速度模型模擬的轉(zhuǎn)換波記錄與原始的觀測(cè)記錄進(jìn)行對(duì)比.當(dāng)正演引擎不能模擬PPS波時(shí)，會(huì)導(dǎo)致走時(shí)匹配紊亂，第二個(gè)PS波同相軸錯(cuò)誤地匹配了觀測(cè)記錄的PPS波同相軸.一旦正演引擎可以同時(shí)模擬PS波和PPS波，即便初始模型對(duì)應(yīng)的模擬數(shù)據(jù)與觀測(cè)數(shù)據(jù)存在較大運(yùn)動(dòng)學(xué)誤差，通過(guò)多輪反演速度模型也能朝著正確方向進(jìn)行更新，最終使模擬和觀測(cè)記錄基本匹配.注意，由于觀測(cè)記錄由彈性波方程合成，而反演中采用的正演引擎是聲波方程Born模擬算子，模擬和觀測(cè)數(shù)據(jù)在振幅、波形方面存在一定偏差是正常的.

圖13 第20炮模擬與觀測(cè)數(shù)據(jù)匹配情況 (a) 正演引擎無(wú)法模擬PPS波的方法； (b) 本文方法.Fig.13 Data matching of 20th common-shot gathers (a) The forward modeling engine can not model PPS waves； (b) The proposed method.

3.2 斷塊模型

通過(guò)斷塊模型合成數(shù)據(jù)檢驗(yàn)方法對(duì)復(fù)雜介質(zhì)情況的適用性.圖14顯示了相應(yīng)的縱、橫波速度與密度模型，其中水深為245 m.這里基于彈性波速度-應(yīng)力方程Born模擬算子合成含PS波與PPS波的轉(zhuǎn)換波數(shù)據(jù)，視為觀測(cè)記錄.震源為主頻25 Hz的雷克子波，共100炮，炮點(diǎn)以40 m的間隔均勻分布于海面，檢波器以5 m的均勻間隔布置于海底，最大炮檢距為1.0 km.第50炮轉(zhuǎn)換波記錄如圖15所示，觀測(cè)數(shù)據(jù)中PS波受到PPS波的嚴(yán)重干擾(如15c箭頭指示).

圖14 斷塊模型 (a) 縱波速度； (b) 橫波速度； (c) 密度.Fig.14 The fault model (a) P-wave velocity; (b) S-wave velocity; (c) Density.

圖15 第50炮轉(zhuǎn)換波記錄 (a) PS波； (b) 源端水層鳴震相關(guān)的PPS波； (c) 合成的轉(zhuǎn)換波觀測(cè)記錄(含PS和PPS波).Fig.15 Converted waves of No.50 common-shot gathers (a) PS waves； (b) PPS waves related to source-side water layer reverberations； (c) The synthetic recordings containing PS and PPS waves.

基于可靠的背景縱波速度(圖16a)、海底界面模型以及反射PP波逆時(shí)偏移剖面(圖16b)進(jìn)行橫波速度反演建模.如圖17a所示，初始橫波速度模型是真實(shí)模型大尺度光滑而成，密度則設(shè)為常數(shù).同樣進(jìn)行兩組試驗(yàn)，第一組反演算法采用的正演引擎僅模擬PS波，同樣發(fā)現(xiàn)因?yàn)闊o(wú)法與含有PPS波的觀測(cè)記錄匹配，導(dǎo)致反演的橫波速度很不合理，尤其在中間斷塊附近存在明顯異常(圖17b).本文方法由于同時(shí)模擬PS與PPS波，使得轉(zhuǎn)換反射波走時(shí)殘差估計(jì)比較可靠，最終獲得了合理的反演結(jié)果(圖17c).圖18顯示的目標(biāo)函數(shù)下降曲線表明，本文方法收斂性更好，最終走時(shí)殘差更小.如圖19所示，在水平位置分別為1 km、2 km和3 km處抽取的偽井速度曲線也表明，本文方法合理地考慮了源端水層鳴震對(duì)轉(zhuǎn)換反射波走時(shí)反演的影響，使橫波速度反演過(guò)程更穩(wěn)健，更準(zhǔn)確.相反，如果模擬的轉(zhuǎn)換波記錄不包含PPS波，反演過(guò)程中走時(shí)殘差最小化過(guò)程發(fā)生模擬PS波與觀測(cè)PPS波串相匹配，導(dǎo)致錯(cuò)誤的模型更新方向，使得橫波速度沿深度交替出現(xiàn)非常異常的劇烈更新(見圖19藍(lán)色曲線).圖20展示了不同橫波速度模型情況下轉(zhuǎn)換PS波逆時(shí)偏移結(jié)果.由于初始模型誤差大，反射界面形態(tài)和深度存在明顯誤差，深層幾個(gè)散射點(diǎn)對(duì)應(yīng)的繞射波沒有收斂.在反演算法中僅模擬PS波導(dǎo)致不合理的走時(shí)殘差估計(jì)進(jìn)而影響速度更新，使得斷塊相關(guān)的界面以及深層散射點(diǎn)成像存在明顯失真.本文方法反演的橫波速度模型基本保證所有界面形態(tài)和深度的準(zhǔn)確性，取得與真實(shí)速度模型偏移結(jié)果相當(dāng)?shù)某上裥Ч?

圖16 (a) 背景縱波速度； (b) PP波逆時(shí)偏移結(jié)果Fig.16 (a) Background P-wave velocity； (b) RTM of reflected PP waves

圖17 (a) 初始橫波速度； (b) 正演引擎無(wú)法模擬PPS波時(shí)的反演結(jié)果； (c) 本文方法反演結(jié)果Fig.17 (a) The initial background S-wave velocity； (b) The result of travel-time inversion with a forward modeling engine that can not simulates PS waves； (c) The proposed method

圖18 目標(biāo)泛函下降曲線Fig.18 Normalized descent curve of objective function

圖19 偽井橫波速度曲線 (a) x=1 km; (b) x=2 km; (c) x=3 km.Fig.19 Pseudo-well curves of the S-wave velocity models at the lateral locations of (a) x=1 km, (b) x=2 km, (c) x=3 km

圖20 反射PS波逆時(shí)偏移結(jié)果 (a) 準(zhǔn)確速度模型； (b) 初始橫波速度模型； (c) 基于圖17b速度模型； (d) 基于圖17c速度模型.Fig.20 Reverse time migration of the reflected PS waves The PS image built by (a) true velocity, (b) the initial background S-wave velocity, (c) the S-wave velocity shown in Fig.17b, (d) the S-wave velocity shown in Fig.17c.

4 討論

本文方法建立在已獲取準(zhǔn)確海底位置、較好完成多分量數(shù)據(jù)P/S分離以及獲得較準(zhǔn)確縱波背景速度和PP波成像結(jié)果基礎(chǔ)上，這些數(shù)據(jù)處理結(jié)果的準(zhǔn)確性直接影響最終橫波速度建模精度.海底位置影響源端水層鳴震P波旅行時(shí)，在縱波背景速度固定不變情況下，由海底位置不準(zhǔn)確引起的鳴震P波走時(shí)擾動(dòng)在反演迭代過(guò)程中會(huì)逐漸投影到橫波背景速度中.此外，縱波背景速度誤差同樣會(huì)迭代投影到橫波速度反演結(jié)果中，導(dǎo)致橫波速度反演結(jié)果誤差增大.在震源子波頻帶不變的情況下，PP波成像剖面分辨率直接影響轉(zhuǎn)換波反偏移數(shù)據(jù)的頻帶.當(dāng)PP波成像結(jié)果分辨率較低時(shí)，相應(yīng)反偏移轉(zhuǎn)換波數(shù)據(jù)整體頻帶會(huì)偏低(表現(xiàn)為記錄中同相軸變粗)，從而影響DIW走時(shí)殘差拾取精度，最終影響反演結(jié)果.數(shù)值實(shí)驗(yàn)表明，當(dāng)橫波初始模型離真實(shí)值偏差較大時(shí)，會(huì)影響一次波與多次波的相對(duì)走時(shí)關(guān)系，給基于DIW算法的走時(shí)反演帶來(lái)挑戰(zhàn)，導(dǎo)致最終無(wú)法獲得準(zhǔn)確的反演結(jié)果.

從正問題描述中發(fā)現(xiàn)，我們利用公式(2)和(4)模擬PS波數(shù)據(jù)，利用公式(2)、(3)、(5)和(6)模擬PPS波數(shù)據(jù).其中公式(3)用來(lái)描述水層散射場(chǎng)(PP波)，僅需包含海底和水層的速度模型，而不需要全空間模型，所以本文方法數(shù)值實(shí)驗(yàn)中正演數(shù)據(jù)模擬階段計(jì)算成本大約為只考慮PS波時(shí)的兩倍.本文提出的梯度預(yù)處理方法需要單獨(dú)求取PS波、PPS波的梯度.相比于只利用PS波，本文方法求取泛函梯度時(shí)計(jì)算成本加倍.

從圖3中我們發(fā)現(xiàn)，當(dāng)海底較淺時(shí)PPS波與PS波路徑基本一致，兩者照明范圍基本相同.隨著水深增加PPS波照明范圍越來(lái)越大，對(duì)橫波背景速度反演的幫助可能會(huì)逐步增大，這一表現(xiàn)后續(xù)需做更詳細(xì)研究.另外，不同水深和海底介質(zhì)條件源端多次波對(duì)橫波速度分析的影響不盡相同.結(jié)合本文工作，今后還應(yīng)當(dāng)從如下幾方面開展深入的研究：首先，將本文方法應(yīng)用到一些淺水探區(qū)的實(shí)際多分量地震資料，在實(shí)踐中進(jìn)行檢驗(yàn)與完善，且與對(duì)轉(zhuǎn)換波記錄進(jìn)行源端多次波壓制，再實(shí)施反射走時(shí)反演的方法流程進(jìn)行對(duì)比，明確不同技術(shù)路線的適用性.其次，雖然不對(duì)轉(zhuǎn)換波數(shù)據(jù)做源端多次波壓制時(shí)，利用本文方法能反演獲得較精確的橫波速度模型.但基于現(xiàn)有轉(zhuǎn)換波成像方法，數(shù)據(jù)中如果含有PPS波等多次波會(huì)造成偏移剖面假象嚴(yán)重.所以后續(xù)研究面向轉(zhuǎn)換波數(shù)據(jù)中包含多次波的相關(guān)成像方法具有重要實(shí)用意義.再次，在硬海底環(huán)境下，在海底界面直接透射轉(zhuǎn)換形成的PSS波會(huì)使轉(zhuǎn)換波走時(shí)反演問題更復(fù)雜.盡管可參照本文思路引入針對(duì)PSS波的Born模擬并探究相應(yīng)的梯度預(yù)條件技術(shù)，但實(shí)際可行性有待驗(yàn)證.此外，在深水(水深超過(guò)500 m)或超深水(水深超過(guò)1500 m)情況下，源端自由表面多次波在水層中傳播時(shí)間長(zhǎng)，主要干擾一些較深層的一次反射與轉(zhuǎn)換波信號(hào).針對(duì)深水區(qū)海底多分量地震橫波速度建模與轉(zhuǎn)換波成像問題開展相關(guān)研究也是很有意義的.

5 結(jié)論

海底多分量地震數(shù)據(jù)經(jīng)過(guò)上-下行P/S波分離之后，獲得的P波和S波記錄仍然受到源端多次波尤其是水層鳴震的干擾，給基于反射PP波和PS波的速度分析和偏移成像造成極大的困擾.近年來(lái)，在全波形反演或反射波形反演理論框架下，出現(xiàn)了同時(shí)考慮一次反射波與自由表面多次波的反演成像方法，但主要還是面向P波數(shù)據(jù)，缺乏針對(duì)S波數(shù)據(jù)的相關(guān)研究.針對(duì)數(shù)據(jù)域P/S分離之后的轉(zhuǎn)換反射波走時(shí)反演問題，本文聚焦淺水、軟海底情況，由海底界面地震波能量分配規(guī)律推斷海底波場(chǎng)分離后的上行S波數(shù)據(jù)主要受源端水層鳴震(尤其是PPS波)干擾.在反演過(guò)程中，如果正演引擎僅考慮PS波，很難與同時(shí)包含PS波、PPS波的觀測(cè)記錄進(jìn)行走時(shí)擬合與殘差估算，導(dǎo)致不合理的模型更新.為此，本文構(gòu)建了描述PS波和PPS波的Born模擬方程，利用伴隨狀態(tài)法提出了同時(shí)考慮這兩種轉(zhuǎn)換波震相的反射走時(shí)反演方法，并引入震相分離梯度預(yù)條件方法，使橫波速度模型的更新方向合理化.水平層狀模型和復(fù)雜斷塊模型合成數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明，該方法為淺水、軟海底環(huán)境基于轉(zhuǎn)換波數(shù)據(jù)的橫波速度建模提供了一種有潛力的技術(shù)手段.

致謝感謝多位審稿人對(duì)本文提出的寶貴建議.

附錄A 基于PS波和PPS波合成轉(zhuǎn)換S波的加權(quán)系數(shù)

源端不同階次水層鳴震P波產(chǎn)生的下行P波場(chǎng)如附圖1 所示.圖中紅線表示一次透射P波，在S波速度擾動(dòng)處產(chǎn)生PS反射波.假設(shè)源端直達(dá)波場(chǎng)能量為1，則其能量為T；藍(lán)線表示源端一階水層鳴震P波，在S波速度擾動(dòng)處產(chǎn)生PPS反射波，其能量為-TR；綠線表示源端二階水層鳴震P波，在S波速度擾動(dòng)處產(chǎn)生PPPS反射波，其能量為TR2.

附圖1 反射S波能量關(guān)系示意圖 T為下行P波在海底處的透射系數(shù)；R為反射系數(shù)； -1為自由表面反射系數(shù).Fig.A1 Energy relationship of the reflected S- waves T represents the transmission coefficient of the down-going P- waves at the seabed；R is reflection coefficient；-1 is the reflection coefficient of free surface.

S波速度擾動(dòng)不變情況下，如果不考慮幾何擴(kuò)散效應(yīng)，反射波能量大小取決于入射波場(chǎng)能量，所以轉(zhuǎn)換波中任意階自由表面多次波與一次波能量關(guān)系表示為

(A1)

式中，n=1,2,3,…，E表示波場(chǎng)能量求取算子，R表示下行P波在海底界面處的反射系數(shù)，可由多分量數(shù)據(jù)P/S分離(劉學(xué)義等,2021)過(guò)程中估計(jì)出的海底介質(zhì)縱橫波速度和密度計(jì)算獲得，或利用其他海底參數(shù)估計(jì)方法(Amundsen and Reitan,1995;Muijs et al.,2004)計(jì)算獲得.公式(7)中的不同波場(chǎng)成分之間相對(duì)幅值需滿足公式(A1)，所以將其分別代入公式(A1)求得的wn即為待求加權(quán)系數(shù).由于走時(shí)反演對(duì)模擬數(shù)據(jù)的振幅精度的要求并不高，根據(jù)介質(zhì)參數(shù)估計(jì)出海底反射系數(shù)之后，結(jié)合實(shí)際數(shù)據(jù)情況給出相對(duì)合理的加權(quán)系數(shù)即可.