利用最小平方基準面延拓壓制自由表面多次波I:海底觀測地震數(shù)據(jù)

朱峰, 程玖兵

同濟大學海洋地質(zhì)國家重點實驗室， 上海 200092

0 引言

主動源海洋地震方法廣泛用于海洋地質(zhì)調(diào)查、油氣與天然氣水合物勘探以及海底殼幔結(jié)構與動力學研究.海洋地震數(shù)據(jù)的采集方式主要分為海面拖纜觀測和海底觀測兩大類，其中后者包括海底電纜(OBC)、海底節(jié)點(OBN)和海底地震儀(OBS)等三種形式.海底檢波器一般由聲壓檢波器(水檢)和三分量質(zhì)點振動速度或加速度檢波器(陸檢)組合而成，所接收的四分量(4C)信號較完整地記錄了地震波三維矢量振動信息，既有利于聯(lián)合縱、橫波改善構造成像、巖性區(qū)分和流體識別，又為從海底觀測記錄分離上-下行波、P/S波創(chuàng)造了有利的數(shù)據(jù)條件.

由于存在強阻抗反差的氣-液自由表面和固-液海底界面，海底多分量地震數(shù)據(jù)波場十分復雜，且發(fā)育多種形式、不同階次的表面多次波，給基于一階散射理論的常規(guī)地震偏移成像或速度反演方法帶來極大的困擾.因此，在海洋地震數(shù)據(jù)處理過程中，壓制表面多次波一直是非常關鍵的環(huán)節(jié).傳統(tǒng)方法要么根據(jù)多次波的周期性或者在某種變換域與一次波能量分布差異進行去除，如預測反褶積法(Peacock and Treitel，1969)、F-K濾波法(Ryu，1982)、Radon變換法(Foster and Mosher，1992)、聚束濾波法(胡天躍和王潤秋，2000)以及擴展成像道集上壓制多次波的方法(Wang et al.，2014)等.要么采用基于波動理論的預測相減法，如模型驅(qū)動的波場外推法(Loewenthal et al.，1974)、數(shù)據(jù)驅(qū)動的反饋迭代法(SRME)(Berkhout and Vershuur，1997)和逆散射級數(shù)法(Weglein et al.，1997).近年來，以WRW傳播模型為正演引擎，在經(jīng)典的SRME方法基礎上發(fā)展出了壓制多次波的稀疏反演法(EPSI)(van Groenestijn and Verschuur，2009a,b；Lin and Herrmann，2013).

自由表面多次波可視為等效面源產(chǎn)生的一次波，比實際點源產(chǎn)生的一次波有更寬的地下照明，可以在地震成像與反演中加以利用(劉伊克等,2018；劉學義等,2021b).如圖1所示，若按照等時原理對多次波進行成像，延拓過程中的上、下行波可能按非物理的波路徑相遇，會在錯誤的位置形成假象.隨著多次波階數(shù)增加，這些非物理串擾會更多、更復雜.若能把多次波分階提取出來，k階多次波作為源端波場進行正傳，k+1階多次波作為檢波器端波場進行反傳，二者的零延遲互相關也可構建地下反射界面或散射點的圖像(劉學建和劉伊克，2016；葉月明等，2019；Liu et al., 2020).此外，相鄰階次多次波之間的互相關可以重構一次波信號(Verschuur and Berkhout，2005；Yu and Schuster，2002；Wapenaar，2006).不過，這些方法的實際應用一直受多次波分階提取精度的制約，而且不同階次多次波成像結(jié)果很難自適應疊加.

圖1 一次反射波與一階多次波聯(lián)合成像示意圖(a) 源端與檢波器端的輸入地震數(shù)據(jù); (b) 互相關成像.Fig.1 Schematic diagram of imaging the primary and first-order multiple reflections(a) Input data at source- and receiver-sides; (b) Imaging via cross-correlation.

迄今為止，處理表面多次波的方法已經(jīng)很多，但都有各自的優(yōu)勢與不足，尤其是對海底觀測地震數(shù)據(jù)的適應性不足(Verschuur and Neumann，1999；劉國昌等，2013；張廣利等，2016).若按產(chǎn)生的位置或時間先后，可把表面多次波分為源端多次波(包括源端鬼波)和檢波器端多次波.在理想情況下，海底聲壓分量(P分量)和質(zhì)點速度垂直分量(Z分量)記錄到的下行波極性相反，上行波極性相同，據(jù)此進行的PZ疊加或PZ相減可經(jīng)上-下行波分離壓制檢波器端的表面多次波.依據(jù)更復雜的彈性波場分離方法，還可在檢波器端上-下行波分離基礎上實現(xiàn)P/S波模式解耦(Amundsen et al., 2000; Muijs et al., 2004).然而，在實際應用中，這些波場分離方法需考慮檢波器在海底的耦合情況和水檢-陸檢儀器響應差異對多分量數(shù)據(jù)進行校正或標定，甚至還需要估計海底介質(zhì)參數(shù)，這就加大了海底觀測數(shù)據(jù)多次波壓制的難度(劉學義等，2021a).在最小二乘理論框架下，把源端正傳的下行波場與檢波器端反傳的上行波場進行迭代反演成像，可以逐步壓制非物理的串擾假象，實現(xiàn)一次波和多次波聯(lián)合成像(Berkhout and Vershuur, 1994; Whitmore et al., 2010；Zhang and Schuster，2014；Wong et al., 2015; Lu et al.，2018).該類方法雖無需分階提取多次波，但需在觀測面上分離上-下行波，或者分離出一次波與多次波.

為了避免多次波分離或分階提取面臨的繁瑣而又困難的工作，本文結(jié)合地震干涉原理(Bakulin and Calvert，2006；Schuster，2009)，借鑒最小平方基準面延拓(LSR)理論框架(Ferguson,2006; Xue and Schuster,2008；Zhu et al., 2020；Zhu and Cheng，2022)，通過迭代反演實現(xiàn)多次波向一次波的轉(zhuǎn)化，從而消除自由表面和海水層對海底地震記錄的影響.本文首先，通過卷積型互易定理推導實際波場狀態(tài)與去除自由表面和上覆介質(zhì)影響的理想波場狀態(tài)之間的數(shù)學物理關系，建立基準面虛擬觀測數(shù)據(jù)同實際觀測數(shù)據(jù)之間的正演方程.然后，結(jié)合海森算子和點擴散函數(shù)揭示常規(guī)干涉法基準面延拓的不足，進而提出基于最小平方反演的基準面延拓方法，并將其用于壓制OBC/OBN數(shù)據(jù)中的表面多次波.最后，通過理論模型合成數(shù)據(jù)實驗證實方法的有效性.

1 考慮自由表面效應的最小二乘基準面延拓方法

1.1 正問題

圖2 兩個獨立的波場狀態(tài)(a) 真實狀態(tài)； (b) 理想狀態(tài).Fig.2 Schematic diagram oftwo wave field states(a) Real state; (b) Ideal state.

根據(jù)附錄A，按照單向波互易定理揭示的傳播不變性(Wapenaar et al., 2004)，可以推導出：

(1)

如圖3，假定在地面或海面激發(fā)震源，在基準面Λ1(如海底或水平井)布放檢波器進行觀測.利用關于R∪的炮檢互易性(Wapenaar，2004)，上式可改寫為

(2)

圖3 受自由表面和上覆介質(zhì)影響的地震觀測示意圖Fig.3 Illustration of seismic observation with the effects of free-surface and overburden medium

1.2 基于干涉原理的基準面延拓

本文基準面延拓的目的是從實際觀測記錄中重構不受自由表面和上覆介質(zhì)影響的、來自目標區(qū)域的反射或背向散射信號R∪(見圖3黑色波路徑).如圖4，方程(2)描述的時間域多維褶積?對應頻率域乘積運算，可寫成如下矩陣形式：

圖4 多維褶積地震干涉正演模擬示意圖Fig.4 Schematic diagram of forward modelingusing multi-dimensional convolution based on seismic interferometry

D-=D+R∪.

(3)

按照地震干涉理論，公式(2)或(3)可利用穩(wěn)相積分原理進行解釋(Schuster，2009).對觀測數(shù)據(jù)施加D+的逆算子，可獲得虛擬觀測數(shù)據(jù)：

R∪=[D+]-1D-.

(4)

這個反問題的直接求解依靠大型矩陣D+的顯式求逆，在實際應用中很難操作.傳統(tǒng)方法一般以伴隨算子D+*(*表示復共軛轉(zhuǎn)置)代替其逆算子，即：

R∪′=D+*D-.

(5)

如圖5，一旦在實際觀測面分離出上、下行波信號，通過(5)式表達的多維互相關，可在一定程度上消除它們在源端的共同波路經(jīng)，通過干涉疊加構建任意兩個檢波器之間的格林函數(shù)場.若把參與互相關的兩個檢波器中的一個視為虛擬震源，對整個實際觀測數(shù)據(jù)如此的干涉疊加最終獲得所謂的虛擬觀測記錄.這種基于干涉原理的基準面延拓方法也被稱為虛源法(Bakulin and Calvert，2006；Schuster，2009).

圖5 多維互相關地震干涉基準面延拓示意圖Fig.5 Schematic diagram of redatuming using multi-dimensional cross-correlation based on seismic interferometry

為了簡潔且不失一般性，分析虛源法輸出信號的具體成分時，暫時僅考慮圖4和圖5中三種有代表性的源端表面多次波路徑.按方程(5)多維互相關得到的基準面延拓數(shù)據(jù)包含九種成分，即

(6)

為了便于理解，將它們以陣列形式展示.如圖6，第一列表示基準面上的上行波信號，第一行表示基準面上的下行波信號的復共軛形式(頻率域的復共軛在時間域相當于對稱翻轉(zhuǎn)，代表從零時刻起順著波路經(jīng)逆時或非因果傳播)，其余主體由3×3矩陣代表(6)式中的數(shù)據(jù)成分.當i=j時，由傳播路徑吻合的上、下行波互相關重構出有效信號，且主對角元素(即R∪′11、R∪′22和R∪′33)走時或相位一致，經(jīng)同相疊加增強有效信號.當i≠j時，傳播路徑不吻合的上、下行波互相關形成串擾噪聲.當?shù)卣鹩涗浐胁豢珊鲆暤母唠A表面多次波時，上述陣列規(guī)模會變大，數(shù)據(jù)成分更多，其中包含階次不對應、路徑不吻合信號互相關形成的噪聲.一般來講，非對角元素的走時或相位差較大，干涉疊加后能量弱于有效信號.

圖6 基準面上-下行波記錄的互相關實線代表正向傳播，虛線代表逆時傳播.Fig.6 Cross-correlation of up- and down-going seismograms at the datumSolid lines denote forward propagation and dashed lines denote reverse-time propagation.

由于自由表面影響以及實際的不完備觀測，使得(5)式中的伴隨算子與下行波逆算子差別較大，導致重構的波場信號包含較多無法干涉相消的串擾噪聲.將(3)式代入(5)式，有

R∪′=HR∪.

(7)

對于完全彈性介質(zhì)，地震波不發(fā)生頻散，基準面延拓正、反問題均可按頻率單獨求解，故有

(8)

注意，等式兩端虛擬觀測數(shù)據(jù)檢波器位置是對應的.方程(7)中H=D+*D+是基準面延拓問題對應的海森算子，可顯式地寫為：

(9)

它攜帶了觀測系統(tǒng)尤其是實際震源的空間分布、時間函數(shù)及其帶限特征，以及自由表面影響下的下行波傳播(甚至散射)效應.也就是說，基于常規(guī)地震干涉原理的基準面延拓數(shù)據(jù)R∪′可以看作是預期的虛擬觀測記錄R∪經(jīng)海森算子濾波(模糊化)的結(jié)果.

1.3 最小平方基準面延拓方法

理論上講，對一次性基準面延拓獲得的波場信號施加海森逆算子H-1，代表反模糊化濾波處理，進而重構高質(zhì)量的虛擬觀測信號，即：

R∪=H-1R∪′，

(10)

其中ε為正則化項的權重因子.如方程(2)和(7)所揭示的，在理想情況下虛擬觀測信號不受震源時間函數(shù)或子波的影響，具有一定的稀疏特性，故采用L1范數(shù)正則化處理.于是，采用梯度類算法迭代估計虛擬觀測信號：

(12)

其中k表示迭代次數(shù)，α表示步長，gk代表數(shù)據(jù)擬合項泛函梯度.注意，正則化處理采用了split Bregman方法(Osher et al.，2005; Yang et al.，2016).詳細計算流程參見算法1.經(jīng)過試算，本文數(shù)值算例中ε取值均為0.01.

參照附錄B，以(2)式為狀態(tài)方程，通過伴隨狀態(tài)法(Plessix，2006)可推導出泛函梯度的表達式：

(13)

它對應觀測與正演信號擬合誤差同下行波信號的多維互相關.第一次迭代的梯度計算與虛源法相當.綜上所述，一旦從實際觀測數(shù)據(jù)中分離出上、下行波信號，就可按這種數(shù)據(jù)驅(qū)動的最小平方基準面延拓(LSR)方法迭代更新基準面上的虛擬觀測數(shù)據(jù)，其中自由表面和上覆介質(zhì)的影響會在這個過程中逐步得到壓制.

算法1: 稀疏約束最小平方基準面延拓

輸入在基準面分離的下行波數(shù)據(jù)D+和上行波數(shù)據(jù)D-;

從k=0到最大迭代次數(shù)N循環(huán)以下計算

計算共軛梯度方向:zk+1=gk+1+βzk;

基準面數(shù)據(jù)更新:qk+1=qk-αzk+1;

2 基于LSR壓制海底地震數(shù)據(jù)表面多次波

海底觀測地震信號可能受到海水層中的導波、海底界面波(Scholte波)以及涌浪等海洋噪聲的干擾.通常會在前處理過程中濾除這些噪聲.就針對反射縱波的數(shù)據(jù)處理而言，利用標定后的陸檢垂直分量與聲壓分量的極性差異，借助PZ疊加或相減可從海底多分量記錄中分離出上行波和下行波信號，從而壓制檢波器端表面多次波對上行波數(shù)據(jù)的影響(如Osen et al., 1999).利用分離出來的上、下行波信號，按前文LSR方法可消除地震波在海水層中的傳播(與散射)效應以及海面相關的源端多次波效應.此時，基準面被選定為實施PZ疊加的海底界面.由于氣槍震源一般位于海表以下，自由表面還會誘發(fā)源端的鬼波效應(圖7a中以序號4標識).當?shù)叵陆橘|(zhì)存在火成巖、膏(鹽)巖以及煤系夾層時，分離出來的下行波數(shù)據(jù)可能包含較強的層間多次散射/反射波.值得強調(diào)的是，本文這種數(shù)據(jù)驅(qū)動的LSR方法可以正確處理源端的鬼波以及層間多次波效應.重構的虛擬觀測信號是源自于基準面之下的背向散射和反射波，其中一階散射/反射波能量最強，是目標區(qū)域地震成像與速度建模最有用的信號.

3 數(shù)值實驗

本文LSR通過數(shù)據(jù)驅(qū)動的迭代反演過程實現(xiàn)在基準面上的虛擬激發(fā)與接收，并剝離上覆介質(zhì)和自由表面的影響.為了不失一般性，首先通過一個簡單的層狀模型演示LSR一次和多次迭代的處理效果，以及海森算子或PSF的作用.為了簡化起見，先僅考慮圖7中序號1、2、3標識的三種一階表面多次波，接著再包含高階表面多次波以驗證它們在稀疏觀測情況下對重構虛擬觀測數(shù)據(jù)的重要貢獻.然后通過一個橫向變速模型合成OBC數(shù)據(jù)，展示LSR方法壓制源端表面多次波(包括鬼波)的應用效果.為了便于控制多次波階次并分析其影響，兩個模型合成疊前炮記錄均采用Berkhout(2014)提出的全波場正演方法，震源采用主頻為20 Hz的Ricker子波.

圖7 海底地震數(shù)據(jù)基準面延拓示意圖(a) 真實狀態(tài); (b) 理想狀態(tài).Fig.7 Schematic diagrams of redatuming for ocean bottom seismograms(a) Real state； (b) Ideal state.

3.1 層狀模型

如圖8，在三層模型的地表等間隔激發(fā)21炮，在深度500 m的水平井中均勻布設101個檢波器，合成炮記錄用于測試LSR算法.圖9顯示了在觀測面分離出的上、下行波信號(暫時只包含一階表面多次波).圖10為經(jīng)LSR重構的、虛擬震源位于基準面中央的反射信號.與圖6呼應，在圖10a顯示的一次迭代獲得的零炮檢距地震道中，通過理論走時計算標識出了各個同相軸對應的波場成分，包括R∪′ii同相疊加構成的強振幅有效信號，以及其他非同相疊加產(chǎn)生的6個串擾噪聲.因為干涉疊加后R∪′21走時趨近于零，波場向零炮檢距聚焦，故其能量也較強(見圖10a與圖10b).隨著迭代次數(shù)增加，串擾噪聲逐漸減弱(圖10c)，殘差快速下降并趨于穩(wěn)定(圖11).經(jīng)過50次迭代，串擾噪聲幾乎被完全壓制，重構反射波信號的垂向分辨率也有明顯提高(見圖10d).

圖8 層狀模型(a) 真實狀態(tài)； (b) 理想狀態(tài). 五角星表示震源位置，三角形表示檢波器位置.Fig.8 A layered model for (a) real state and (b) ideal stateThe stars indicate the source location, and the triangles indicate the geophone location.

圖9 基準面上的聲壓分量(a)和速度垂直分量(b)及分離的下行波(c)和上行波(d)記錄Fig.9 Seismograms of (a) pressure and (b) vertical particle velocity，and the decomposed (c) up-going and (d) down-going signals at the datum

圖10 最小平方基準面延拓實驗(a) 第1次迭代后的零炮檢距地震道； (b) 第1次、(c)第5次和(d)第50次迭代重構的單炮道集.Fig.10 Least-squares redatuming experiment(a) The retrieved zero-offset seismic trace of the first iteration; (b), (c) and (d) represent the retrieved common-shot gathers of the 1st , 5th and 50th iterations, respectively.

圖11 層狀模型實驗殘差收斂曲線Fig.11 Misfit curve of the experiment on the layered model

按前文分析可知，LSR一次迭代的結(jié)果可視為虛擬觀測數(shù)據(jù)被海森算子模糊化濾波得到的，而且海森矩陣某一列與特定炮檢對的PSF相對應.換個角度看，濾波后的點脈沖虛源被改造為以相應PSF為輻射模式的復合型虛源，從而在反射響應中耦合了自由表面、觀測系統(tǒng)以及震源子波等效應.在本例中，下行波信號包含直達波和一階自由表面多次波，沿物理波路徑傳播的上、下行波信號的互相關代表多次波向一次波的干涉轉(zhuǎn)化和同相疊加，其他非對應路徑上、下行波互相關產(chǎn)生無法干涉相消的串擾腳印.LSR迭代反演本質(zhì)上相當于對圖10b一次基準面延拓的炮記錄施加圖12中PSF的逆算子對應的滑動濾波(或虛源輻射模式矯正)，通過去模糊化消除串擾噪聲，重構基準面上虛擬激發(fā)與接收的高分辨率反射波信號.

圖12 在基準面中央激發(fā)和接收的虛擬炮檢對所對應的點擴散函數(shù)Fig.12 A point spread function corresponds to a virtual survey of which both the source and receiver are in the mid of the datum

3.2 高階表面多次波的貢獻

圖13 某合成炮記錄中分離的下行波與上行波信號(a,b) 不包含表面多次波； (c,d)含1階表面多次波; (e,f)含1至3階表面多次波.Fig.13 Separated down- and up-going waves of a synthetic common-shot gatherCommon-shot gather (a,b) without free-surface multiple, (c,d) with 1st-order free-surface multiples, and (e,f) with 1st to 3rd-order free-surface multiples.

為了揭示高階表面多次波在干涉法基準面延拓中的貢獻，下面在模型兩端相隔1 km各激發(fā)一炮數(shù)據(jù)，通過LSR算法重構出虛擬震源等間隔分布在基準面的21炮疊前地震記錄.圖14顯示了經(jīng)50次迭代之后其中代表性的5炮道集.可見，由于高階表面多次波參與基于干涉疊加的基準面延拓，模型中部弱照明問題得以明顯緩解，預期的雙曲型反射信號重構效果得到了明顯改善(箭頭指示).這表明一階乃至高階表面多次波提供的寬廣照明對重構基準面之下的背向散射或反射波信號是非常有益的.否則，僅靠稀疏觀測炮記錄中的直達波和一次反射波信號，難以消除重構炮記錄中目標反射同相軸在零炮檢距附近的運動學誤差(圖中圈定區(qū)域)，甚至在缺少實際震源的模型中央完全無法重構有效信號，出現(xiàn)較強的假信號.

圖14 含不同階次表面多次波炮記錄LSR輸出的虛擬炮道集 (a) 原始炮記錄僅含一次直達波和反射波； (b) 原始炮記錄還包含1階表面多次波； (c) 原始炮記錄還包含1至3階表面多次波.Fig.14 Reconstructed virtual common-shot recording by LSR with shot gathers having different order of free-surface multiplesTheoriginal shot gathers with (a) direct wave and primary reflection; (b) and with 1st-order free-surface multiples; (c) and with 1st- to 3rd -order free-surface multiples.

3.3 復雜模型OBC數(shù)據(jù)

接下來將LSR方法用于處理橫向變速模型合成OBC數(shù)據(jù)，檢驗它能否壓制海水層和自由表面對縱波信號的影響.在圖15a展示的速度模型中，第一層是海水，第三、四層底界面存在較大起伏.假定由氣槍震源從橫向1～9 km范圍內(nèi)在海表以下10 m深處等間隔激發(fā)401炮，震源時間函數(shù)采用主頻為20 Hz的Ricker子波.在海底橫向2～8 km范圍內(nèi)均勻布設301個檢波器，數(shù)值模擬401炮含有鬼波與多次波的原始OBC疊前道集.然后采用PZ疊加/相減，分離出上行波、下行波記錄(如圖16)用于后續(xù)LSR實驗，擬重構基準面(海底)均勻分布151個虛擬震源在301個檢波器接收到的背向散射與反射信號.

圖15 速度模型與OBC觀測系統(tǒng)(a) 真實狀態(tài)； (b) 理想狀態(tài). (五角星代表震源，三角形代表檢波器. Fig.15 Velocity model and OBC survey geometry for (a) real state and (b) ideal stateThe stars denote the sources and the triangles denote the receivers.

圖16 實際震源位于橫向坐標5 km處的OBC炮道集(a)原始聲壓場信號和(b)速度垂直分量，以及分離出的(c)下行波和(d)上行波信號.Fig.16 An OBC common-shot gather for a source at 5 km(a) Recorded pressure and (b) vertical particle velocity seismograms, and the decomposed (c) down-going and (d) up-going wavefield signals.

首先，按照海底的虛擬觀測系統(tǒng)，采用Born模擬合成震源位于橫向坐標5 km處一次反射波單炮記錄(圖17a)，用于對比評估LSR處理結(jié)果.圖17b和17c顯示了LSR第1次和第50迭代輸出的虛擬震源也在5 km處的單炮記錄.而且，LSR最終輸出的反射波波形更窄，頻譜變得更寬(見圖17c與18c).這是因為LSR既壓制了鬼波與表面多次波造成的陷波效應(對比圖18a與18c)，又一定程度去除了子波效應，拓寬了信號頻帶，提高了時間分辨率(見圖17c和圖18c).可見，隨著迭代次數(shù)增加，串擾噪聲逐步得到壓制，殘差穩(wěn)定減弱(圖19)，一次反射波信號得以有效重構.如圖20，對比不同炮道集對應的疊前深度偏移結(jié)果，容易發(fā)現(xiàn)上-下行波分離消除了檢波器端表面多次波對上行波數(shù)據(jù)成像的干擾(白色箭頭)，而多輪迭代LSR處理壓制了源端表面多次波相關的假象(黑色箭頭)，基本消除了對反射界面成像的串擾，提高了成像剖面的分辨率和信噪比.

圖17 基準面重構數(shù)據(jù)(a) Born正演參考數(shù)據(jù)； (b) 一次迭代重構數(shù)據(jù)； (c) 50次迭代重構數(shù)據(jù).Fig.17 Datum reconstructed data(a) Reference data from Born forward modeling; (b) Direct redatuming data and (c) LSR data obtained with Fig.15a as input.

圖18 原始與重構炮道集單道頻譜對比(a) 原始上行波信號(與圖16d對應)； (b) Born模擬反射波(與圖17a對應); (c) LSR重構反射波記錄(對應圖17c).Fig.18 Comparison of the spectra for the up-going wave in the original and the retrieved common-shot gathers(a) Original upgoing waves (corresponding to Fig.16d); (b) Primary reflection synthesized by Born modelling (corresponding to Fig.17a); (c) Retrieved signal by LSR (corresponding to Fig.17c).

圖19 復雜模型實驗的殘差收斂曲線Fig.19 Misfit curve of the experiment on the complex model

圖20 輸入不同炮道集的疊前深度偏移結(jié)果(a) 原始炮道集； (b) 分離后的上行波道集； (c) 經(jīng)1次LSR迭代重構的虛擬炮道集； (d) 經(jīng)50次LSR迭代重構的虛擬炮道集.Fig.20 Comparison of prestack depth migration results for various common-shot gathers(a) Input original gathers without up/down separation; (b) Input up-going waver field signal; Input reconstructed common-shot gathers after the (c) 1st and (d) 50th iteration of LSR.

4 討論

針對OBN/OBC數(shù)據(jù)處理是一個復雜的系統(tǒng)工程，多次波壓制需要一些必要的、合理的前處理步驟進行支撐.在完成檢波器方向與耦合校正以及水檢-陸檢標定之后，聯(lián)合聲壓與垂直速度分量的PZ疊加/相減分離出上、下行縱波信號，這就基本消除檢波器端表面多次波對上行縱波記錄的干擾.在此基礎上，本文LSR方法通過多次迭代，干涉轉(zhuǎn)化或消去的是源端表面多次波.如果數(shù)據(jù)中存在較強的橫波能量，需要在海底觀測面采用彈性波場分離方法(Muijs et al., 2004；劉學義等,2021a)，則可在分離上-下行波之后進一步實現(xiàn)縱波與橫波信號的模式解耦.理論上講，以分離的上、下行P/S波作為輸入，采用本文方法能夠重構出已克服自由表面和海水層影響的PP和PS虛擬反射數(shù)據(jù).這將在今后研究中加以驗證.

值得注意的是，對大多數(shù)OBN/OBC數(shù)據(jù)來說，震源分布較為稠密，滿足干涉法對空間采樣的要求(即不存在圖13那種極端稀疏觀測情況).與OBC采集方式不同，由于OBN節(jié)點或者海底地震儀(OBS)造價比較昂貴，檢波器分布相對稀疏一些.雖然LSR方法比較準確地重構出了任意兩個檢波器之間的格林函數(shù)場，形成海底虛擬觀測記錄，但若檢波器間隔太大，會給后續(xù)偏移成像或參數(shù)反演帶來挑戰(zhàn).對于海洋拖纜采集方式，若采用傳統(tǒng)的單分量聲壓檢波器，從觀測記錄分離上、下行波比較困難.若采用水-陸雙檢(dual sensor)觀測方式，則可用PZ疊加/相減分離上、下行縱波信號.本文姊妹篇將討論如何利用LSR原理消除拖纜地震數(shù)據(jù)中的表面多次波.

本文LSR方法完全是數(shù)據(jù)驅(qū)動的，與地下速度模型無關，因此直接適用于各向異性、黏滯性介質(zhì).它涉及的正演和泛函梯度計算分別通過上-下行波多維褶積和多維互相關實現(xiàn).在CPU型號為Intel Xeon E7-8890的多核工作站上，前文二維簡單模型和復雜模型的每一炮單次迭代平均耗時分別為1.2 s和18.4 s，總的計算量不大，但虛擬觀測數(shù)據(jù)需要較大的內(nèi)存開銷，拓展到三維情況要從計算策略上加以優(yōu)化.對深水或深層目標勘探而言，經(jīng)LSR輸出的虛擬觀測數(shù)據(jù)已經(jīng)剝離了自由表面和上覆介質(zhì)的影響，顯著壓縮了模型空間和數(shù)據(jù)規(guī)模，有利于開展面向局部目標的精細反演成像，預期也可在油藏開發(fā)或二氧化碳填埋的時移地震監(jiān)測中發(fā)揮積極作用.

5 結(jié)論

一次基準面延拓利用多維互相關消除上、下行波場中共同的波路經(jīng)，從而在一定程度上剝離自由表面和上覆海水對基準面重構數(shù)據(jù)的影響，但仍然會產(chǎn)生能量較強的、由非匹配路徑波場信號互相關引起的串擾噪聲.借助最小二乘優(yōu)化反演思想，本文LSR方法不僅能夠很好地消除這些串擾噪音，還可一定程度壓制不完備觀測、鬼波效應和子波帶限效應的影響，提升虛擬觀測信號的分辨率和頻帶寬度.該方法充分利用表面多次波的照明優(yōu)勢和最小平方反演迭代過程，無需對表面多次波進行比較困難的分離或分階提取即可實現(xiàn)向一次波的干涉轉(zhuǎn)化.數(shù)值實驗表明重構數(shù)據(jù)基本消除表面多次波的干擾，改善后續(xù)基于一階散射理論的疊前偏移成像效果.

致謝感謝王騰飛老師、鄒鵬博士、武泗海博士和黃河碩士等在本文研究過程中提供的建議.

附錄A：基于互易定理構建實際觀測和虛擬觀測波場狀態(tài)之間的轉(zhuǎn)化關系

欲利用模型無關的基準面延拓方法克服表面多次波的影響，需建立實際觀測和理想狀態(tài)(無自由表面效應)虛擬觀測之間波場信號的轉(zhuǎn)化關系.這里采用的理論工具是由亥姆霍茲方程和高斯定理推導的單向波互易定理(Wapenaar et al. 2004).它揭示了封閉無源區(qū)域兩個相互獨立狀態(tài)下波場的卷積或相關沿著封閉邊界的積分是守恒的.卷積型互易定理用于描述波場傳播問題，而相關型互易定理主要用于波場反向延拓、逆時聚焦以及干涉成像等(Schuster，2009；Wapenaar，2004).下面利用卷積型互易定理建立從虛擬觀測數(shù)據(jù)到實際觀測數(shù)據(jù)的映射(正演)方程.

根據(jù)文獻(de Hoop，1988；Fokkema and van den Berg，1993)，假定存在兩個相互獨立的波場狀態(tài)A和B，它們在兩個水平無限延伸界面Λi和Λj包圍的封閉區(qū)域內(nèi)介質(zhì)特性保持一致且沒有震源(封閉區(qū)域外介質(zhì)特性可能不同，且可能存在震源)，則單向波卷積型互易定理在頻率-空間域具有如下形式(Wapenaar,1996；Wapenaar et al.,2004)：

(A1)

令Λi=Λ1，Λj=Λ2，將上、下行波格林函數(shù)場代入公式(A1)表征的褶積型互易定理，有：

(A2)

(A3)

它在時間域?qū)亩嗑S褶積給出了地表脈沖震源激發(fā)、基準面檢波器實際記錄的上行波G-同在基準面虛擬激發(fā)、接收的上行波R∪之間的關系.

附錄B：基于伴隨狀態(tài)法推導LSR泛函梯度

(B1)

針對LSR采用的L2范數(shù)目標泛函，形成如下約束優(yōu)化問題：

=0.

(B2)

通過拉格朗日乘子法，將其轉(zhuǎn)化為非約束優(yōu)化問題，即：

(B3)

其中R表示取實部操作，λ為伴隨狀態(tài)變量.由于F(R∪,D-)=0，故有J(R∪)=L(R∪,D-,λ).

方程(B2)對應的最小化問題等價于尋求增廣泛函L(R∪,D-,λ)的鞍點，對應其一階偏導數(shù)為0的解，即：

(B4)

(B5)

(B6)

其中δD-表示模擬數(shù)據(jù)與觀測數(shù)據(jù)的殘差.

由于λ與R∪無關，故有

(B7)

(B8)