結合Marchenko聚焦的非線性最小平方基準面延拓方法

朱峰, 程玖兵

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

0 引言

地震勘探一般在地表或海面附近激發(fā)、接收地震信號.源自深部的反射或背向散射信號會受到地面/海表、上覆介質(zhì)和觀測系統(tǒng)的復雜影響,具體包括海面或地表相關的自由表面多次波干擾、觀測面起伏或表層橫向變速造成的運動學畸變、不規(guī)則巖體(如火成巖、鹽丘與煤層等)引起的多次散射、障礙物(如鉆井平臺、湖泊等)或禁采區(qū)導致的非規(guī)則采樣等.這些線性或非線性效應嚴重降低了深部有效信號的質(zhì)量,給后續(xù)速度建模和偏移成像帶來挑戰(zhàn).

基準面延拓旨在由地表真實觀測數(shù)據(jù)重構基準面上的虛擬觀測數(shù)據(jù),從而克服上覆介質(zhì)和觀測因素對反射地震數(shù)據(jù)的影響.按照是否依賴地下速度模型,基準面延拓分為模型驅(qū)動和數(shù)據(jù)驅(qū)動兩類方法(Schuster and Zhou, 2006).模型驅(qū)動的基準面延拓根據(jù)地下宏觀速度并利用Kirchhoff積分、單程波或雙程波延拓算子依次把檢波器和炮點校正到基準面上,獲得虛擬觀測的地震數(shù)據(jù)(Berryhill, 1979; Reshef, 1991).它主要被用于消除起伏地表影響(楊鍇等, 2007; 盧回憶等, 2010),數(shù)據(jù)規(guī)則化(Ferguson and Fomel, 2005)、改善異常地質(zhì)體之下的波場照明和構造成像(Dong et al., 2009)等.當層間多次波等非線性波場效應較強時,這類方法需要提供精確的地下速度結構用于模擬全波場響應,以實現(xiàn)嚴格的波場延拓(Mulder, 2005).由于波場傳播物理假設和觀測條件限制,構建精確的速度模型十分困難,因而阻礙了模型驅(qū)動類的基準面延拓方法有效應對非線性波場效應(Planès et al., 2018).

數(shù)據(jù)驅(qū)動的基準面延拓方法不依賴于速度模型.基于井中或海底觀測數(shù)據(jù)攜帶的格林函數(shù)信息,以地震干涉為核心的虛源法(Bakulin and Calvert, 2006; Schuster, 2009)通過多維互相關把實際檢波器轉(zhuǎn)化成虛擬震源,進而消除波場在上覆介質(zhì)中的傳播過程.當上覆介質(zhì)存在強阻抗界面或復雜地質(zhì)體時,波場會發(fā)生多次散射.由于觀測數(shù)據(jù)攜帶的格林函數(shù)包含全波場信息,虛源法可自動考慮上覆介質(zhì)中的非線性波場效應.然而,由于觀測系統(tǒng)不完備(即不在封閉面上充分、規(guī)則地采集信號)和子波帶限效應,干涉法重構的虛擬數(shù)據(jù)僅是目標信號的模糊化近似,在不同程度上受到了同非物理路徑和觀測孔徑限制相關的假信號污染(朱峰和程玖兵, 2022).更重要的是,當檢波器布設面與基準面不重合時,如果沒有宏觀速度模型等必要的先驗信息,很難從實際觀測數(shù)據(jù)獲得干涉法需要的上、下行波格林函數(shù).這很大程度上限制了虛源法的應用領域.幸運的是,隨著解決一維介質(zhì)單邊逆散射問題的Marchenko方程(Marchenko, 1955)在本世紀以來受到地球物理界越來越多的關注,由地表觀測的反射地震數(shù)據(jù)估計觀測面與地下介質(zhì)中任意點之間的格林函數(shù)場形成了比較堅實的理論基礎(Rose, 2001; Broggini and Snieder, 2012; Wapenaar et al., 2014a).利用上覆介質(zhì)宏觀速度模型計算的直達波聚焦函數(shù),Marchenko方法可以從地表觀測的反射數(shù)據(jù)迭代估計地表到基準面的單向波格林函數(shù)(包含上覆介質(zhì)中的多次散射等非線性效應),進而借助上、下行格林函數(shù)的多維反褶積(MDD)預測或衰減層間多次波(Wapenaar et al., 2012; van der Neut, 2012; Behura et al., 2014; Meles et al., 2016; Zhang et al., 2019; 匡偉康等, 2018; 孫紅日等, 2021; 張樂樂等, 2022),或者實現(xiàn)考慮層間多次波的基準面延拓或地震成像(Wapenaar et al., 2014b; Ravasi, 2017; 靳中原等, 2017).不過,MDD涉及顯式計算大量的點擴散函數(shù)和互相關函數(shù),存儲消耗非常大.本質(zhì)上講,Marchenko方法是由模型與數(shù)據(jù)聯(lián)合驅(qū)動的.由于反射地震觀測不完備和數(shù)據(jù)的帶限效應,MDD面臨的超大規(guī)模不適定問題求解起來很不穩(wěn)定,而且對實際觀測面上震源與檢波器的分布有較嚴苛的要求(van der Neut and Herrmann, 2013).

在最小平方反演框架下實現(xiàn)基準面延拓,可在迭代過程中對虛擬觀測數(shù)據(jù)施加恰當?shù)募s束或正則化處理,更好地壓制采集腳印和非物理波場噪聲的干擾.Ferguson(2006)最早將最小平方基準面延拓(LSR)用于數(shù)據(jù)規(guī)則化處理.接著,Xue和Schuster(2008)發(fā)現(xiàn)常規(guī)波動方程基準面延拓中因觀測孔徑限制引起的假信號可由LSR進行壓制.針對模型驅(qū)動的LSR問題,為了改善數(shù)據(jù)擬合條件并盡可能剝離上覆介質(zhì)的影響,Zhu和Cheng(2022)提出了同時反演上覆地層反射率和基準面虛擬觀測數(shù)據(jù)的雙參數(shù)LSR方法,并結合東海油氣田應用實例展示了它在改善深層地震成像方面的潛力.最近,針對典型的海洋地震數(shù)據(jù)采集方式,朱峰和程玖兵(2022),朱峰等(2023)提出了基于數(shù)據(jù)驅(qū)動LSR壓制表面多次波的新方法.

針對檢波器布設面與基準面不一致而且上覆介質(zhì)中發(fā)育較強層間多次波的復雜情況,本文擬把Marchenko與LSR方法相結合,建立可自動處理層間多次波的非線性基準面延拓方法.

1 方法原理

1.1 基準面到地表的波場轉(zhuǎn)化方程

正問題的建立是求解反問題的前提.波場互易定理(Haines, 1988)描述了地下任意兩個水平無限延伸界面上波場之間滿足的積分恒等關系,為建立不同深度界面上波場相互轉(zhuǎn)化方程提供了理論依據(jù).它既被用于構建基準面上虛擬觀測數(shù)據(jù)到單向格林函數(shù)的正演方程,又支撐Marchenko方法以聚焦函數(shù)為橋梁從地面觀測反射信號迭代估計基準面上的單向波格林函數(shù).

假定存在兩個相互獨立的波場狀態(tài)A和B,它們在兩個水平無限延伸邊界Λi和Λj包圍的封閉區(qū)域內(nèi)部介質(zhì)完全相同且無震源分布(封閉區(qū)域外介質(zhì)特性和震源分布無要求),則兩個邊界上的單向波場在頻率域滿足以下卷積型波場互易定理(Wapenaar et al., 2004):

(1)

和相關型波場互易定理:

(2)

其中xi和xj分別是邊界Λi和Λj上的水平坐標,P代表邊界上的通量歸一化單向波場(上標“+”和“-”分別表示下行波和上行波,下標“A”和“B”分別代表上述兩種波場狀態(tài),上標“*”表示對頻率域復數(shù)波場取共軛,對應時間域的波場翻轉(zhuǎn)).注意,相關型互易定理要求封閉區(qū)域內(nèi)無非彈性能量耗散.

表1 下伏介質(zhì)封閉區(qū)域邊界上單向波場的匯總表Table 1 One-way wavefield variables on the boundaries of an enclosed domain beneath the datum

圖1 以基準面下伏介質(zhì)為封閉區(qū)域的波場狀態(tài)示意圖Fig.1 Schematic diagram of the wavefield states in an enclosed domain beneath the datum

將表1中的波場代入(1)式所示的卷積型互易定理,整理后得到:

(3)

該式給出了虛擬觀測數(shù)據(jù)和基準面上的上、下行格林函數(shù)之間的關系.當基準面上恰好布設了足夠多的多分量檢波器,如拖纜雙檢、VSP、海底電纜(OBC)或海底節(jié)點(OBN)觀測,則可通過觀測數(shù)據(jù)的聲壓分量與質(zhì)點垂直振動速度(即P/Z)疊加/相減分離出上、下行波場(Amundsen et al., 2000; 劉學義等, 2021).本文著重考慮檢波器布設面與基準面不重合、上覆介質(zhì)存在層間多次波等非線性波場效應的情況,此時基準面上單向波格林函數(shù)的估計就成為了核心問題.

1.2 重構格林函數(shù)的Marchenko聚焦方法

表2 封閉上覆介質(zhì)邊界上的單向波場匯總表Table 2 One-way wavefields on the boundary of the enclosed overburden media

圖2 以上覆介質(zhì)為封閉區(qū)域的波場狀態(tài)示意圖Fig.2 Schematic diagram of the wavefields in the enclosed overburden media

(4)

方程組(4)含有四個未知變量(即F+、F-、G+和G-),無法直接求解.這里采用Broggini和Snieder(2012)以及Wapenaar等(2014a)提出的迭代求解方法(見附錄A).基本思路如下,將下行波聚焦函數(shù)、格林函數(shù)分別拆分成直達波和尾波兩種成分,其中下行直達波聚焦函數(shù)可依據(jù)已知的上覆介質(zhì)宏觀速度模型正演模擬出來.利用聚焦函數(shù)與格林函數(shù)在時間分布區(qū)間上的差異,可將方程組(4)轉(zhuǎn)化成僅含下行聚焦函數(shù)尾波和上行聚焦函數(shù)兩個待求變量的迭代求解公式.把迭代估計的聚焦函數(shù)代入Marchenko方程組并考慮格林函數(shù)在時間軸上的分布區(qū)間,就可以重構出地表到基準面的上、下行波格林函數(shù).

1.3 基于Marchenko聚焦的最小平方基準面延拓

根據(jù)前文推導,從基準面虛擬觀測數(shù)據(jù)到上行波格林函數(shù)的正演方程可寫成矩陣形式:

G-=G+R∪.

(5)

這表明,基準面上的下行格林函數(shù)G+與虛擬反射數(shù)據(jù)R∪經(jīng)多維褶積運算,能夠合成有更長波路徑的上行格林函數(shù)G-(如圖3).注意,下行格林函數(shù)不僅包含直達透射波(黑色路徑),還含有在上覆介質(zhì)中形成的層間多次波(紅色路徑)以及經(jīng)過下伏介質(zhì)的層間多次波(紫色路徑)等非線性波場效應.盡管(5)式可以考慮高階多次散射,但為了簡潔起見,圖中僅展示了一階層間多次散射.

圖3 多維褶積正演示意圖Fig.3 Schematic diagram of forward modeling by multidimensional convolution

基準面延拓對應(5)式的逆過程,故有:

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

(6)

(6)式表示多維反褶積(MDD)運算.由于直接構建下行格林函數(shù)的稀疏矩陣并求逆非常困難,通常把MDD轉(zhuǎn)化為如下求解公式:

(7)

其中G+*G-和G+*G+分別對應互相關函數(shù)和點擴算函數(shù)的集合(van der Neut,2012).如果按前文Marchenko方法構建單向波格林函數(shù),則(6)或(7)式就是基于MDD的Marchenko基準面延拓方法(Wapenaar et al., 2014b).

Marchenko基準面延拓在實際應用中面臨一系列挑戰(zhàn).首先,對于常規(guī)反射地震觀測系統(tǒng)而言,觀測孔徑有限且炮點、檢波器空間位置不重合,MDD隱含的地震干涉過程不足以消除基準面上覆介質(zhì)中波場效應的影響,估算的虛擬觀測數(shù)據(jù)存在相干的假信號.其次,多維互相關函數(shù)或點擴散函數(shù)集合對應的矩陣規(guī)模為2×Nw×(Ns×Nr)2,其中Nω、Ns和Nr分別是頻率采樣、炮點和檢波點個數(shù).對于典型反射地震觀測系統(tǒng)(取Nω≈50,二維情況下Ns=Nr≈1000,三維情況下Ns=Nr≈1000×1000),以單精度浮點類型進行數(shù)據(jù)存儲時,這些矩陣的理論存儲消耗可達400TB(二維)和4×108TB(三維),對現(xiàn)有的計算機硬件條件構成極大挑戰(zhàn).而且,這種超大規(guī)模MDD運算也很不穩(wěn)定(van der Neut and Herrmann, 2013).

為了避免顯式計算與操作多維互相關函數(shù)和點擴散函數(shù),同時方便對估計的虛擬觀測數(shù)據(jù)施加恰當?shù)恼齽t化處理(以減少Marchenko基準面延拓中的假信號干擾),在最小平方反演理論框架下求解(6)式代表的反問題.于是,構建如下目標函數(shù):

(8)

等式右邊包含基于L2范數(shù)的數(shù)據(jù)擬合項和基于L1范數(shù)的模型(即虛擬觀測數(shù)據(jù))正則化處理或約束優(yōu)化項,其中ε為權重因子.

方程(8)采用梯度類的局部優(yōu)化迭代算法進行求解:

(9)

其中k表示迭代次數(shù),α表示步長,gk代表泛函梯度,其表達式可由伴隨狀態(tài)法(Plessix, 2006; 朱峰和程玖兵, 2022)推導得到:

(10)

(10)式中δG-=G+R∪-G-表示模擬與Marchenko方法估計的上行格林函數(shù)之間的殘差.帶正則化處理的迭代算法可參照朱峰和程玖兵(2022).對于最小平方基準面延拓(LSR)方法,存儲消耗量級約為4×Nt×Ns×Nr(其中Nt是時間采樣點數(shù)),主要包括上、下行格林函數(shù)、泛函梯度以及虛擬觀測數(shù)據(jù).這種最小平方反演算法的存儲規(guī)模遠小于基于MDD的算法.

1.4 算法流程

根據(jù)上述方法原理,基于Marchenko聚焦的非線性LSR方法的實現(xiàn)可歸納為如圖4所示的三步流程:首先估計聚焦函數(shù),即根據(jù)上覆介質(zhì)宏觀速度模型正演模擬透射直達波,以它和反子波(如脈沖反褶積)的地表反射數(shù)據(jù)作為輸入,按附錄A中公式(A9)迭代估算上、下行聚焦函數(shù).然后基于它們和地表反射脈沖響應數(shù)據(jù)按公式(A10)重構基準面上的上、下行格林函數(shù)(場).對于這些單向波格林函數(shù)而言,震源仍然在地表,但檢波器已經(jīng)延拓(沉降)到了基準面.最后將上、下行格林函數(shù)代入約束優(yōu)化LSR反演算法,最終輸出炮點、檢波器均分布在基準面上的虛擬觀測數(shù)據(jù).

圖4 基于Marchenko聚焦的最小平方基準面延拓算法流程Fig.4 Workflow of least-squares redatuming based on Marchenko focusing

2 數(shù)值實驗

本部分通過兩個合成數(shù)據(jù)實驗檢驗非線性LSR方法的有效性.第一,通過含低速夾層的水平層狀模型算例展示處理流程關鍵步驟的結果,結合波路徑分析以及與參考答案的比較揭示方法的工作機制.第二,利用高速透鏡體模型算例考察本文方法對于上覆介質(zhì)強多次散射的適用性.為了避免脈沖反褶積處理,這里在基于全波場模擬方法(Berkhout,2014)正演在地表觀測的反射地震數(shù)據(jù)時,直接采用帶通Ormsby子波(圖5a和5b).

圖5 用于合成脈沖反射數(shù)據(jù)的Ormsby子波(a) 波形; (b) 振幅譜圖.Fig.5 The Ormsby wavelet used to simulate impulsive reflected data(a) Waveform; (b) Amplitude spectrum.

2.1 含低速夾層的水平層狀模型

在如圖6a所示的層狀模型中,101個炮點和檢波器被均勻布設在地表,基準面設定在500 m深處(白色虛線位置).圖6b展示了由真實速度場光滑得到的基準面上覆宏觀速度模型,用于正演直達透射波數(shù)據(jù)以近似下行聚焦函數(shù)的直達波成分(圖7a).由于上覆介質(zhì)含有低速夾層,波場模擬結果中除了含一次反射波,還含有層間多次波(圖7b).為了便于觀察,圖7c顯示了將層間多次波振幅增益10倍后的結果以及其中4個能量最強的一階層間多次波對應的波路徑.

圖6 含低速夾層的水平層狀模型(a) 準確速度場; (b) 上覆宏觀速度場. 圖a中五角星指示了聚焦點位置.Fig.6 The layered model with a low speed interlayer(a) True velocity model; (b) Overlying macro velocity model. In figure a, the pentagram indicates the focusing point.

圖7 用于Marchenko迭代的輸入數(shù)據(jù)(a) 聚焦點到地表的透射直達波; (b) 震源位于地表0.5 km的反射數(shù)據(jù); (c) 層間多次波及主要同相軸波路徑.Fig.7 The input data for the Marchenko method(a) Transmitted primary wave from the focusing point to the surface; (b) Impulsive reflected data for sources at 0.5km on the surface; (c) Internal multiples and wave paths of the main events.

由Marchenko方法構建格林函數(shù)共分為兩步,第一步,計算聚焦函數(shù).將透射直達波和地表反射數(shù)據(jù)代入時窗約束的Marchenko方程((A9)式),迭代10次后得到上、下行聚焦函數(shù)(圖8).在下行聚焦函數(shù)(圖8a)中有兩個同相軸,其中標號1是直達波成分,其波路徑(圖8c)與透射直達波(圖7a)相同但走時相反;標號2是數(shù)值尾波,它和直達波成分極性相反且走時差是低速地層的雙程走時,這樣的時空特征保證了它剛好能夠抵消低速地層頂界面的下行反射作用(圖8c標號5),從而壓制上覆介質(zhì)層間多次波.如圖8b和8c所示,下行聚焦函數(shù)經(jīng)上覆介質(zhì)一次反射后形成上行聚焦函數(shù).第二步通過聚焦函數(shù)實現(xiàn)檢波點端波場聚焦((A10)式),從而得到地表激發(fā)、基準面接收的上、下行格林函數(shù)(圖9).圖9b展示了下行格林函數(shù)及其波路經(jīng),其中能量占主導的是透射直達波(標號1)以及一階層間多次相關的尾波(標號2-4).下行格林函數(shù)經(jīng)過深部地層界面背向反射后形成上行格林函數(shù)(圖9d).與同樣用全波場模擬算法在基準面上以橫向0.5 km的聚焦點為震源另外合成的地表參考數(shù)據(jù)(圖9a和9c)相比,本文Marchenko算法估計的格林函數(shù)和它們相同.這里,上行格林函數(shù)波路徑的前面部分與其對應的下行格林函數(shù)的波路徑一致(圖9b和9d),這部分會在后續(xù)多維互相關運算中被消去.

圖8 橫向0.5 km處的聚焦函數(shù)(a)下行和(b)上行聚焦函數(shù)以及它們的(c)波路徑. 圖c中點線表示對頻率域波場取共軛或時間域波場翻轉(zhuǎn), 在物理上解釋為由負時刻開始傳播的非因果波場.Fig.8 Focusing functions at lateral 0.5 km Down-going (a) and up-going (b) focusing functions and their wave paths (c). The dotted line in Fig. c represents taking conjugate of the wavefield in frequency domain or flipping of the wavefield in time domain. It can be physically explained as non-causal wave which propagates from negative time.

圖9 橫向0.5 km處的格林函數(shù)(a) 參考下行格林函數(shù); (b) 重構的下行格林函數(shù)及其波路徑; (c) 參考上行格林函數(shù); (d) 重構的上行格林函數(shù)及其波路徑.Fig.9 Green′s functions at lateral 0.5 km(a) Synthetic down-going Green′s functions for reference; (b) Constructed down-going Green′s functions and their wavepaths; (c) Synthetic up-going Green′s functions for reference; (d) Constructed up-going Green′s functions and their wavepaths.

利用估算的上、下行格林函數(shù)進行LSR處理.LSR第一輪迭代等同于基于地震干涉的虛源法(Bakulin and Calvert, 2006).在該過程中,多維互相關抵消了上、下行格林函數(shù)波路經(jīng)的重疊部分,隨后同相疊加的是有效的虛擬觀測信號(圖10b),而非同相疊加部分就形成沿非物理路徑傳播的串擾噪聲(圖10b白色箭頭).如圖10c,這些串擾會在LSR迭代30次后基本去除.然而,由于子波旁瓣壓制不徹底,有效信號周圍出現(xiàn)了平行的震蕩噪聲.通過施加稀疏約束正則化處理,震蕩噪聲得到壓制,重構數(shù)據(jù)的分辨率也得到提升(圖10d).

圖10 震源位于基準面上橫向0.5 km的虛擬反射數(shù)據(jù)(a) 參考數(shù)據(jù); (b) 直接干涉數(shù)據(jù); (c) LSR反演數(shù)據(jù); (d) 稀疏約束LSR反演數(shù)據(jù).Fig.10 Virtual reflected data with sources at 0.5 km on the datum(a) Reference data; (b) Direct interferometry data; (c) LSR data; (d) Sparse constrained LSR data.

2.2 高速透鏡體模型

利用近地表強橫向變速與高速透鏡體復合模型合成數(shù)據(jù)測試算法對復雜地質(zhì)情況的適用性.在圖11a展示的速度模型上,采用固定排列觀測系統(tǒng)在地表橫向0～4 km范圍均勻布設401個炮點和檢波器,基準面設置在1.2 km深度(圖中白色虛線位置).在基準面上的虛擬觀測系統(tǒng)與地表實際觀測系統(tǒng)保持一致.由于上覆介質(zhì)存在強反差界面,在地表合成的炮記錄中廣泛發(fā)育層間多次波(圖11b白色箭頭).

圖11 透鏡體模型(a) 速度場; (b) 地表脈沖反射數(shù)據(jù). 圖a中五角星指示了聚焦點位置.Fig.11 The lens model(a) Velocity model; (b) Impulsive reflected data on the surface. In figure a, the pentagram indicates the focusing point.

圖12a和12b展示了利用Marchenko方法經(jīng)20次迭代構建的下行和上行聚焦函數(shù).與簡單模型算例(圖8)相比,因上覆阻抗界面增多,聚焦函數(shù)同相軸也顯著增加.這里的下行聚焦函數(shù)的數(shù)值尾波和直達波極性相同(圖12a),這是由于隨著深度增加,速度(阻抗)逐漸增大,上覆界面的下行反射系數(shù)為負,想要抵消它們產(chǎn)生的極性翻轉(zhuǎn)的下行多次波,數(shù)值尾波的極性需與直達波一致.將聚焦函數(shù)代入(A10)式計算由地表到基準面的上、下行格林函數(shù),其中下行格林函數(shù)(圖12c)由上到下依次為透射直達波(黑色箭頭)和層間多次震蕩形成的尾波.如圖12d,下行直達波和尾波經(jīng)深部地層背向反射回到基準面分別形成上行格林函數(shù)的一次反射波(黑色箭頭)和層間多次波.

圖12 橫向2 km處的聚焦函數(shù)和格林函數(shù)(a) 下行聚焦函數(shù); (b) 上行聚焦函數(shù); (c) 下行格林函數(shù); (b) 上行格林函數(shù).Fig.12 Focusing functions and Green′s functions at lateral 2 kmDown-going (a) and up-going (b) focusing functions; Down-going (c) and up-going (d) Green′s functions.

格林函數(shù)直接干涉重構的虛擬數(shù)據(jù)含有明顯的非物理路徑串擾噪聲(圖13a),經(jīng)稀疏約束LSR迭代30次后進行了有效壓制(圖13b).此外,LSR還抵消了上、下行格林函數(shù)中的子波信息,重構出來自基準面以下的脈沖反射響應.為了方便對比,在虛擬數(shù)據(jù)上褶積了主頻20Hz的雷克子波(圖13c),并抽取了零偏移距剖面(圖14).在地表觀測的反射數(shù)據(jù)中(圖14a),透鏡體頂?shù)捉缑嫦嚓P的層間多次波(白色箭頭)與深部有效信號(黑色箭頭)混疊在一起,不利于深部構造解釋.如圖14b所示,虛源法校正了高速透鏡體對深部反射的扭曲,但也產(chǎn)生了許多串擾噪聲.通過LSR處理,這些噪聲得到有效壓制(圖14c).MDD(Thorbecke et al., 2017)的基準面延拓結果(圖13d和圖14d)同樣具有壓制層間多次波串擾和反子波效果.然而,由于頻率域數(shù)據(jù)矩陣運算難以施加預條件和正則化處理,即便采用了大型稀疏矩陣的穩(wěn)定求解算法(LSQR),反演結果中仍殘留有部分噪聲和干擾.

圖13 基準面反演的虛擬數(shù)據(jù)(a) 直接干涉重構數(shù)據(jù); (b) LSR反演數(shù)據(jù); (c) 褶積雷克子波后的LSR反演數(shù)據(jù); (d) MDD反演數(shù)據(jù).Fig.13 Inverted virtual data on the datum(a) Direct interferometry data; (b) LSR data; (c) LSR data convolved with Ricker wavelet; (d) MDD data.

圖14 零偏移距剖面(a) 地表反射數(shù)據(jù); (b) 直接干涉重構數(shù)據(jù); (c) LSR反演數(shù)據(jù); (d) MDD 反演數(shù)據(jù).Fig.14 Zero offset profiles(a) Surface reflected data; (b) Direct interferometry data; (c) LSR data; (d) MDD data.

3 結論與討論

針對檢波器布設面與基準面不重合且上覆介質(zhì)存在強非線性波場效應的情況,本文構建了基于Marchenko聚焦的非線性最小平方基準面延拓方法.它以聚焦函數(shù)為橋梁,由地表反射數(shù)據(jù)估計格林函數(shù),進而通過約束優(yōu)化的LSR方法反演虛擬數(shù)據(jù).相比于基于MDD的Marchenko基準面延拓方法,LSR一方面避免了龐大的點擴散和互相關函數(shù)集合的顯式計算、存儲,以及大規(guī)模稀疏矩陣的不穩(wěn)定求逆;另一方面,通過對虛擬數(shù)據(jù)施加L1范數(shù)正則化,有效的壓制了非物理路徑串擾噪聲以及觀測孔徑和子波帶限的不利影響,在實現(xiàn)壓制層間多次波的同時,重構了在基準面上虛擬觀測的高分辨率反射數(shù)據(jù),為面向深部目標的成像奠定了良好的基礎.

本文方法在實際應用時還面臨如下問題:首先,上、下行格林函數(shù)的估計依賴于運動學準確的透射直達波場,當上覆宏觀速度模型存在較大誤差時,會導致重構的虛擬數(shù)據(jù)不準確;其次,Marchenko迭代要求輸入反子波的地表單位脈沖響應,這就需要對每一炮數(shù)據(jù)預先估計震源子波,并通過反褶積對其加以消除;最后,經(jīng)典的Marchenko方法要求輸入的觀測數(shù)據(jù)滿足炮檢點耦合,如何擺脫這一限制仍然是當前的研究熱點(Ravasi, 2017).上述挑戰(zhàn)將在今后的研究中加以考慮.

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

附錄A 利用Marchenko聚焦求解單向波格林函數(shù)

由正文可知,基于波場互易定理推導的單邊格林函數(shù)表示定理滿足:

(A1)

參照Broggini和Snieder(2012)以及Wapenaar等(2014a),把下行波聚焦函數(shù)和格林函數(shù)拆分為直達波與尾波,即:

(A2)

(A3)

(A4)

(A5)

(A6)

針對格林函數(shù)及其共軛定義相應的時窗函數(shù):

(A7)

注意Θ(t)與Ψ(t)是互補的.將時窗函數(shù)Ψ(t)作用于(A1)式,可得Marchenko方程:

(A8)

其中I和I-1分別為正、反傅里葉變換.此時,方程組(A8)僅包含兩個未知數(shù),可以聯(lián)立(A5)式構建如下迭代求解公式:

(A9)

將迭代估計的上、下行波聚焦函數(shù)代回(A1)式,并施加時窗函數(shù)Θ(t),可得上、下行波格林函數(shù),即:

(A10)

可見,Marchenko方法是以聚焦函數(shù)為橋梁,從地表觀測反射波脈沖響應估計地表到基準面的上、下行波格林函數(shù).如果基于上覆介質(zhì)宏觀速度正演模擬的直達波含有子波信息,經(jīng)過(A9)式迭代估計的聚焦函數(shù)以及隨后得到的格林函數(shù)均含有相同的子波信息.