朱峰, 程玖兵

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

0 引言

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

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

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

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

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

1 方法原理

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

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

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

(1)

和相關(guān)型波場互易定理:

(2)

其中xi和xj分別是邊界Λi和Λj上的水平坐標(biāo),P代表邊界上的通量歸一化單向波場(上標(biāo)“+”和“-”分別表示下行波和上行波,下標(biāo)“A”和“B”分別代表上述兩種波場狀態(tài),上標(biāo)“*”表示對(duì)頻率域復(fù)數(shù)波場取共軛,對(duì)應(yīng)時(shí)間域的波場翻轉(zhuǎn)).注意,相關(guā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ǔn)面下伏介質(zhì)為封閉區(qū)域的波場狀態(tài)示意圖Fig.1 Schematic diagram of the wavefield states in an enclosed domain beneath the datum

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

(3)

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

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

1.3 基于Marchenko聚焦的最小平方基準(zhǔn)面延拓

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

G-=G+R∪.

(5)

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

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

基準(zhǔn)面延拓對(duì)應(yīng)(5)式的逆過程,故有:

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

(6)

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

(7)

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

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

為了避免顯式計(jì)算與操作多維互相關(guān)函數(shù)和點(diǎn)擴(kuò)散函數(shù),同時(shí)方便對(duì)估計(jì)的虛擬觀測數(shù)據(jù)施加恰當(dāng)?shù)恼齽t化處理(以減少M(fèi)archenko基準(zhǔn)面延拓中的假信號(hào)干擾),在最小平方反演理論框架下求解(6)式代表的反問題.于是,構(gòu)建如下目標(biāo)函數(shù):

(8)

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

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

(9)

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

(10)

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

1.4 算法流程

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

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

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

本部分通過兩個(gè)合成數(shù)據(jù)實(shí)驗(yàn)檢驗(yàn)非線性LSR方法的有效性.第一,通過含低速夾層的水平層狀模型算例展示處理流程關(guān)鍵步驟的結(jié)果,結(jié)合波路徑分析以及與參考答案的比較揭示方法的工作機(jī)制.第二,利用高速透鏡體模型算例考察本文方法對(duì)于上覆介質(zhì)強(qiáng)多次散射的適用性.為了避免脈沖反褶積處理,這里在基于全波場模擬方法(Berkhout,2014)正演在地表觀測的反射地震數(shù)據(jù)時(shí),直接采用帶通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個(gè)炮點(diǎn)和檢波器被均勻布設(shè)在地表,基準(zhǔn)面設(shè)定在500 m深處(白色虛線位置).圖6b展示了由真實(shí)速度場光滑得到的基準(zhǔn)面上覆宏觀速度模型,用于正演直達(dá)透射波數(shù)據(jù)以近似下行聚焦函數(shù)的直達(dá)波成分(圖7a).由于上覆介質(zhì)含有低速夾層,波場模擬結(jié)果中除了含一次反射波,還含有層間多次波(圖7b).為了便于觀察,圖7c顯示了將層間多次波振幅增益10倍后的結(jié)果以及其中4個(gè)能量最強(qiáng)的一階層間多次波對(duì)應(yīng)的波路徑.

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

圖8 橫向0.5 km處的聚焦函數(shù)(a)下行和(b)上行聚焦函數(shù)以及它們的(c)波路徑. 圖c中點(diǎn)線表示對(duì)頻率域波場取共軛或時(shí)間域波場翻轉(zhuǎn), 在物理上解釋為由負(fù)時(shí)刻開始傳播的非因果波場.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) 重構(gòu)的下行格林函數(shù)及其波路徑; (c) 參考上行格林函數(shù); (d) 重構(gòu)的上行格林函數(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ù)進(jìn)行LSR處理.LSR第一輪迭代等同于基于地震干涉的虛源法(Bakulin and Calvert, 2006).在該過程中,多維互相關(guān)抵消了上、下行格林函數(shù)波路經(jīng)的重疊部分,隨后同相疊加的是有效的虛擬觀測信號(hào)(圖10b),而非同相疊加部分就形成沿非物理路徑傳播的串?dāng)_噪聲(圖10b白色箭頭).如圖10c,這些串?dāng)_會(huì)在LSR迭代30次后基本去除.然而,由于子波旁瓣壓制不徹底,有效信號(hào)周圍出現(xiàn)了平行的震蕩噪聲.通過施加稀疏約束正則化處理,震蕩噪聲得到壓制,重構(gòu)數(shù)據(jù)的分辨率也得到提升(圖10d).

圖10 震源位于基準(zhǔn)面上橫向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 高速透鏡體模型

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

圖11 透鏡體模型(a) 速度場; (b) 地表脈沖反射數(shù)據(jù). 圖a中五角星指示了聚焦點(diǎn)位置.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次迭代構(gòu)建的下行和上行聚焦函數(shù).與簡單模型算例(圖8)相比,因上覆阻抗界面增多,聚焦函數(shù)同相軸也顯著增加.這里的下行聚焦函數(shù)的數(shù)值尾波和直達(dá)波極性相同(圖12a),這是由于隨著深度增加,速度(阻抗)逐漸增大,上覆界面的下行反射系數(shù)為負(fù),想要抵消它們產(chǎn)生的極性翻轉(zhuǎn)的下行多次波,數(shù)值尾波的極性需與直達(dá)波一致.將聚焦函數(shù)代入(A10)式計(jì)算由地表到基準(zhǔn)面的上、下行格林函數(shù),其中下行格林函數(shù)(圖12c)由上到下依次為透射直達(dá)波(黑色箭頭)和層間多次震蕩形成的尾波.如圖12d,下行直達(dá)波和尾波經(jīng)深部地層背向反射回到基準(zhǔn)面分別形成上行格林函數(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ù)直接干涉重構(gòu)的虛擬數(shù)據(jù)含有明顯的非物理路徑串?dāng)_噪聲(圖13a),經(jīng)稀疏約束LSR迭代30次后進(jìn)行了有效壓制(圖13b).此外,LSR還抵消了上、下行格林函數(shù)中的子波信息,重構(gòu)出來自基準(zhǔn)面以下的脈沖反射響應(yīng).為了方便對(duì)比,在虛擬數(shù)據(jù)上褶積了主頻20Hz的雷克子波(圖13c),并抽取了零偏移距剖面(圖14).在地表觀測的反射數(shù)據(jù)中(圖14a),透鏡體頂?shù)捉缑嫦嚓P(guān)的層間多次波(白色箭頭)與深部有效信號(hào)(黑色箭頭)混疊在一起,不利于深部構(gòu)造解釋.如圖14b所示,虛源法校正了高速透鏡體對(duì)深部反射的扭曲,但也產(chǎn)生了許多串?dāng)_噪聲.通過LSR處理,這些噪聲得到有效壓制(圖14c).MDD(Thorbecke et al., 2017)的基準(zhǔn)面延拓結(jié)果(圖13d和圖14d)同樣具有壓制層間多次波串?dāng)_和反子波效果.然而,由于頻率域數(shù)據(jù)矩陣運(yùn)算難以施加預(yù)條件和正則化處理,即便采用了大型稀疏矩陣的穩(wěn)定求解算法(LSQR),反演結(jié)果中仍殘留有部分噪聲和干擾.

圖13 基準(zhǔn)面反演的虛擬數(shù)據(jù)(a) 直接干涉重構(gòu)數(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) 直接干涉重構(gòu)數(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 結(jié)論與討論

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

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

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

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

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

(A1)

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

(A2)

(A3)

(A4)

(A5)

(A6)

針對(duì)格林函數(shù)及其共軛定義相應(yīng)的時(shí)窗函數(shù):

(A7)

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

(A8)

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

(A9)

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

(A10)

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