程玖兵, 王騰飛, 徐文才

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

0 引言

地球內(nèi)部存在復雜的多尺度的非均勻性.它決定了彈性參數(shù)(如縱、橫波速度與密度等)的空間分布與變化特征，影響地震波的傳播與散射，甚至導致波場衰減與各向異性效應.在主動源地震探測過程中，人們利用地震記錄(包括初至與反射波等震相)去構建地下彈性參數(shù)模型，定位和描述介質(zhì)內(nèi)部的不連續(xù)界面，進而解析地質(zhì)構造與沉積模式，預測巖性、孔隙流體以及應力場的空間分布，服務于油氣及其他礦產(chǎn)資源勘探開發(fā)、區(qū)域地質(zhì)或巖石圈構造演化研究和一些工程與環(huán)境地質(zhì)任務.

以地震波的主波長為參照，彈性非均勻性可以分解到不同尺度，包括長波長、中波長和短波長分量.一般認為，前兩個分量構成光滑的“背景模型”，控制波場運動學(如走時與相位)，而短波長分量構成震蕩的“擾動模型”，主要影響波場動力學(如振幅).因此，在實踐中，背景模型和擾動模型的建立通常是分開進行的(Claerbout,1985).建立背景模型尤其是反演縱、橫波宏觀速度主要依靠走時層析或偏移速度分析.在反射地震勘探領域，基于疊前深度偏移共成像點道集剩余曲率信息的反射走時層析至今仍是速度建模的主要手段(Stork，1992；Woodward et al.,2008).背景模型一旦確定，基于反射或單次散射信號估計彈性參數(shù)擾動可視為一個線性反問題.要么基于線性正演算子的伴隨算子，由射線或波動理論基礎上的地震偏移生成反射界面圖像(Claerbout,1971; Schneider,1978；McMechan,1983)；要么利用正演算子的近似逆，借助線性逆散射(Miller et al.,1987；Beylkin and Burridge,1990; 毛偉建等，2021)或真振幅偏移(Bleistein et al.,2001;張宇,2006)，甚至基于迭代的線性波形反演或最小二乘偏移，即LSM(Clayton and Stolt，1981;Tarantola,1984a；Guy and Rene-Edouard,1999;任浩然等,2013;陳生昌,2018)定量估計阻抗界面的反射率或彈性參數(shù)擾動.

將背景模型與擾動模型反演過程完全分開，不便于考慮兩者的耦合關系，也影響反演的收斂性和分辨率.于是，在高性能計算平臺的支撐下，無需進行模型尺度分解，以模擬與觀測數(shù)據(jù)所有震相擬合誤差最小化為目標的非線性全波形反演，即FWI(Lailly,1983；Tarantola,1984b)在近十年重新引起了關注.它具有“走時層析”和“偏移成像”雙重功能，可捕捉不同角度散射場信號對彈性參數(shù)長、中、短波長分量的敏感性，因此具備高分辨率建模潛力(Mora,1989；Neves and Singh,1996).由于人工源地震數(shù)據(jù)缺少有效的低頻成分，F(xiàn)WI的成功案例多限于透射波、回轉波或折射波等大角度波場信號充足的井間地震、長炮檢距地面地震數(shù)據(jù)(Pratt,1999; Brenders and Pratt,2007).自2005年前后基于雙程波方程的逆時偏移(RTM)投入實際生產(chǎn)以來，F(xiàn)WI作為配套的速度建模技術在一些探區(qū)尤其是針對海洋長纜觀測地震數(shù)據(jù)取得了很好的實用效果(如Ratcliffe et al.,2014).然而，大量常規(guī)地震數(shù)據(jù)的炮檢距與方位角分布范圍有限，大角度波場很難穿透到中、深層，F(xiàn)WI極易受非線性和“周波跳躍”問題的困擾(Virieux and Operto,2009; 董良國等,2013).即便采用從低頻到高頻(Bunks et al.,1995)、層剝離(Wang and Rao,2009)或散射角濾波(Alkhalifah,2016)等多尺度策略，或者采用更凸的目標泛函(如Luo et al.,2016; Yang and Engquist,2018)或者擴展模型(如Huang et al.,2018)，都只能部分克服FWI在處理反射波數(shù)據(jù)時面臨的挑戰(zhàn).

面向反射地震數(shù)據(jù)的非線性反演可以不采用經(jīng)典FWI的理論框架，除了依據(jù)反射波數(shù)據(jù)擬合誤差設計目標函數(shù)，還可以依據(jù)反射波的成像效果(如聚焦性)設計目標函數(shù).數(shù)據(jù)域的非線性反射波形反演，即RWI(Clément et al.,2001;Xu et al.,2012)和成像域的波動方程偏移速度分析,即WEMVA(Sava and Biondi,2004;Shen and Symes,2008)都屬于這一范疇.它們均采用模型尺度分解思想，只不過把重構背景與擾動模型置于一個聯(lián)合的反問題之中，由嵌套在一起的兩個“子反演”交替進行求解.成像域反射走時層析通過疊前深度偏移估計擾動模型或反射率，針對背景模型的目標函數(shù)凸性好，對初始模型要求較低，對應的反演分辨率也低一些(Díaz et al.,2013).數(shù)據(jù)域反射波形擬合目標函數(shù)凸性差，仍受“周波跳躍”影響，故通常在初始迭代階段采用走時擬合標準，當長波長分量得以重構之后再切換到波形擬合標準，進一步彌補一些模型的中波長分量(如Ma and Hale,2013; Wang et al.,2018).

由于反射數(shù)據(jù)對大尺度的密度變化不敏感，而與宏觀速度之間存在很強的非線性關系，曾有學者嘗試用全局優(yōu)化方法去反演背景速度模型(如Snieder et al.,1989).考慮到計算可行性，人們在處理實際問題時幾乎都采用局部優(yōu)化方法.在迭代的局部尋優(yōu)過程中，目標泛函關于模型參數(shù)的一階偏導數(shù)和二階偏導數(shù)(即泛函梯度與海森算子)共同定義了模型的更新方向.為了降低算法復雜度，最速下降與共軛梯度等一階優(yōu)化方法被廣泛用于FWI、LSM和RWI等反問題(Virieux and Operto,2009;黃建平等,2014;姚剛和吳迪,2017).不過，泛函梯度易受觀測孔徑與數(shù)據(jù)頻帶限制引起的模糊化作用的影響，甚至受多參數(shù)串擾腳印與多次散射效應的污染，故導致反演過程收斂較慢，精度也有較大提升空間.二階優(yōu)化思想就是引入海森算子對泛函梯度的預條件處理，從而改善收斂性、提高分辨率、壓制多參數(shù)耦合(Pratt et al.,1998;Métivier et al.,2013;Wang et al.,2016; Pan et al.,2017).根據(jù)采用精確還是近似海森算子，出現(xiàn)了擬牛頓法(如L-BFGS)、高斯-牛頓法以及(全)牛頓法等不同的二階優(yōu)化算法實現(xiàn)(Pratt et al.,1998; Nocedal and Wright,2006).近10年來，已有不少文獻討論針對FWI、LSM等反問題的二階優(yōu)化方法(如高鳳霞等,2013;劉璐等,2013; 任浩然等,2013;Métivier et al.,2013;Wang et al.,2016;Pan et al.,2017;Rao and Wang，2017).但是，針對RWI的二階優(yōu)化方法研究才剛開始，關于海森算子的特征、二階優(yōu)化RWI的實際效果以及所涉及的算法復雜性都值得去調(diào)查(Cheng et al.,2020;Wang et al.,2021;徐文才,2020).

本文面向觀測孔徑有限的常規(guī)反射地震數(shù)據(jù)，在二階優(yōu)化反演理論框架下，推導背景和擾動模型的反射波敏感核、泛函梯度以及海森算子，提出并驗證基于近似海森矩陣的高斯-牛頓反射波形反演方法.論文結構如下：首先采用波動方程的Born近似推導反射波形擬合目標函數(shù)關于背景模型和擾動模型的泛函梯度與海森算子，借助簡單模型分析海森矩陣的特征及對泛函梯度的預條件作用.然后建立基于高斯-牛頓算法的二階優(yōu)化反射波形反演方法.最后利用理論模型合成數(shù)據(jù)和東海實際地震資料，檢驗該方法的有效性與實用性.

1 二階優(yōu)化反射波形反演方法

1.1 正問題

地震波傳播過程通常用波動方程來描述，在頻率域可簡要表示為如下矩陣-矢量形式：

L(x,ω)u(x,ω)=f(x,ω)，

(1)

其中u為波場矢量，f代表震源矢量,x代表空間坐標；L稱為復“阻抗”矩陣，是由頻率ω和地下介質(zhì)彈性參數(shù)m(x)定義的傳播算子，描述波場與參數(shù)之間的非線性關系.在聲學近似意義下，u通常指聲壓場；在彈性介質(zhì)中，u一般包含質(zhì)點的水平、垂直位移(或速度)分量.基于聲壓波場理論的地震成像與反演最關心縱波速度這個模型參數(shù)，有時也考慮介質(zhì)的密度變化和吸收衰減.欲更貼切地描述地震波傳播的物理過程，需采用更復雜的波動理論，引入橫波速度、各向異性參數(shù)或品質(zhì)因子等介質(zhì)參數(shù).

大量觀測與實驗表明，地下介質(zhì)的彈性非均勻性具有復雜的多尺度特征.在地震波成像與反演應用中，彈性參數(shù)模型一般被分解為光滑的背景模型m0和震蕩的擾動模型m1，即

m(x)=m0(x)+m1(x).

(2)

這種尺度分解一般以地震主波長為參照，比如背景模型包含長波長與中波長參數(shù)結構，控制波的傳播，而擾動模型由刻畫短波長特征的參數(shù)分量構成，導致波場散射(Tarantola,1984a;Jannane et al.,1989).假設擾動分量遠小于背景分量，總波場可分解為入射場u0與一階散射場u1之和，即u=u0+u1.于是(1)式演化成一階Born模擬方程：

L(x,ω)u0(x,ω)=f(xs,ω)，

(3a)

L(x,ω)u1(x,ω)=Q(x,ω)u0(x,ω),

(3b)

1.2 反射波形反演

波動方程反射波形反演(RWI)以模擬與觀測反射波數(shù)據(jù)波形殘差最小化為目標，除了重構彈性模型的宏觀(或長波長)結構,還可適當恢復部分中波長分量(Xu et al.,2012;徐文才,2020).由于反射波u1受背景模型m0與擾動模型m1共同影響，同時估計這兩個模型分量面臨極強的非線性與復雜的耦合問題.本文將這個反問題分解成嵌套在一起的兩個子反演，在迭代過程中交替更新模型的背景與擾動分量，結合頻率延拓策略實現(xiàn)寬譜建模目標.

于是，估計背景模型的反問題可定義為

(4)

(5)

其中Δdm0=u1(m0)|m1-dobs代表給定擾動模型條件下因背景模型誤差導致的反射數(shù)據(jù)擬合誤差，“?”代表共軛轉置.由于背景模型與反射數(shù)據(jù)之間存在復雜的非線性關系，由(4)和(5)式定義的反問題的非線性很強.

類似地,估計擾動模型的反問題可定義為

(6)

其中Δdm1=u1(m1)|m0-dobs為基于更新后的背景模型，由擾動模型誤差導致的反射數(shù)據(jù)擬合誤差.由式(3b)可知，擾動模型與反射波響應成線性關系，因此(6)式本質(zhì)上對應一個線性反射波形反演問題,即所謂的最小二乘逆時偏移問題(Tarantola,1984a; Guy and Rene-Edouard, 1999).假定有n次觀測(相當于有n個炮檢對)，有p類不同的物理參數(shù)(比如考慮縱、橫波速度與密度時，p=3)，以q代表每個物理參數(shù)離散結點數(shù)，則模型參數(shù)l=q×p，于是m0和m1是維數(shù)為l的向量，傳播矩陣L維數(shù)為q×q；對單個頻率成分，dobs和Δdm1是維數(shù)為n的向量，u0和u1是維數(shù)為q的向量.

在每個頻段交替求解上述兩個反問題，均可在局部優(yōu)化理論框架下采用如下迭代格式：

(7)

1.2.1 一階優(yōu)化

(8)

(9)

于是，對方程(3a)和(3b)關于背景模型求一階偏導數(shù)，得到相應的Fréchet微商(推導見附錄A)：

(10)

或寫成Jm0=L-1F0，其中F0是維數(shù)為q×l的虛源矩陣，其列向量滿足：

(11)

同理，方程(3a)和(3b)關于擾動模型求一階偏導數(shù)，得到相應的Fréchet微商：

(12)

或寫成Jm1=L-1F1，其中虛源矩陣F1的列向量滿足：

(13)

在特定參數(shù)化方式下，關于擾動模型的Fréchet微商是不依賴于模型m1的.例如，對于常密度聲介質(zhì)，若用慢度平方來表征模型，則(13)式右端偏微分項僅與頻率有關.因此，在Born近似意義下，針對擾動模型的波形反演可視為線性問題.

如果先計算模型空間每個虛源點的Fréchet微商，后按(8)式構建泛函梯度，會涉及太多波場模擬，在實際應用中不可取.人們通常采用伴隨狀態(tài)法，借助正向波場和逆向(或伴隨)波場的零延遲互相關來計算泛函梯度(Tarantola,1984b; Tromp et al.,2005; Plessix,2006).根據(jù)附錄A中的推導，背景模型的泛函梯度滿足：

(14)

根據(jù)附錄B中的推導，擾動模型的泛函梯度可表示成：

(15)

圖1 背景模型與擾動模型泛函梯度計算示意圖

在交替反演過程中，一旦擾動模型獲得更新，就由Born模擬計算正向和伴隨散射場，進而構建關于背景模型的泛函梯度，并按(9)式更新背景模型；接下來重新模擬正向和伴隨入射場，構建擾動模型的泛函梯度并按(9)式更新擾動模型.如此交替進行直至達到終止條件.當初始模型誤差較小時，共軛梯度方法可以穩(wěn)健地逼近真實模型.

1.2.2 二階優(yōu)化

牛頓法認為誤差函數(shù)在初始模型m附近滿足二次型，故可將其做二階泰勒近似：

(16)

其中海森算子H由誤差函數(shù)在初始模型m處的二階偏導數(shù)構成，即：

(17)

假定模型空間包含p類參數(shù)，每類含有q個離散結點(即l=q×p)，則海森算子可表示成維數(shù)為l×l的方陣.若將其排布成p×p個矩陣塊，則對角塊蘊含與觀測孔徑、波場照明以及數(shù)據(jù)帶限相關的信息，非對角塊主要反映非同類參數(shù)之間的耦合效應，也攜帶了帶限數(shù)據(jù)條件下p類參數(shù)空間結點之間的關聯(lián)信息(Pratt et al.,1998;Wang et al.,2016).對于單參數(shù)問題，海森矩陣維度降為q×q，其對角元素與孔徑、照明有關，非對角元素反映帶限數(shù)據(jù)條件下空間結點之間的關聯(lián)性或耦合性.

(18)

引入海森算子可較好地克服因孔徑、照明與頻帶限制對泛函梯度的模糊化作用，使模型更新更均衡、更合理.這將在后文數(shù)值實驗中得到驗證.

一般來講，海森矩陣可拆解成兩項：

(19)

其中第一項只涉及Fréchet微商，定義了所謂的近似海森矩陣：

(20)

第二項涉及波場關于模型參數(shù)的二階偏導數(shù)，即

Rmj=

(21)

當數(shù)據(jù)殘差較大或非線性波場效應(如多次散射)較強時，海森算子第二項的作用不可忽視.對于FWI或LSM問題，考慮這一項有助于處理多次散射效應(Pratt et al.,1998; Métivier et al.,2013).不過，第二項可能會破壞矩陣正定性，在求解牛頓方程時要考慮這一因素.在局部線性假設條件下，可求解如下高斯-牛頓方程：

(22)

獲得模型更新方向.因為近似海森矩陣具有正定性，高斯-牛頓法總能給出合理的下降方向，并避免海森算子非線性項引入的計算復雜性.

1.2.3 海森矩陣數(shù)值分析

雖然本文二階優(yōu)化RWI方法不需要顯式計算海森矩陣及其逆，但為了觀察海森矩陣的特點，本節(jié)以常密度聲介質(zhì)單參數(shù)(即速度)反演為例，先給出反射波敏感核與海森矩陣的計算式，然后基于小模型演示背景與擾動速度海森矩陣的差異.給定震源s和檢波器r的單次觀測，背景慢度的反射波敏感核按(A2)可表示成：

Jm0(x,r,s,ω)=ω2[u0(x,s,ω)G1(x,r,ω)

+u1(x,s,ω)G0(x,r,ω)]

=ω2f(ω)[G0(x,s,ω)G1(x,r,ω)

+G1(x,s,ω)G0(x,r,ω)].

(23)

它涉及入射場與一階散射場的頻率域乘積或時間域褶積，揭示了反射波路徑上震源端、檢波器端的二次散射對背景慢度變化的敏感性.代(23)入(20)式，近似海森矩陣的元素由模型參數(shù)結點x與x′對應的Fréchet微商時間域互相關構成，在頻率域寫成：

×[G0(x′,s,ω)G1(x′,r,ω)]

×[G1(x′,s,ω)G0(x′,r,ω)],

(24)

其中“*”代表復共軛，右邊第一項代表震源端Fréchet微商零延遲互相關，第二項代表檢波器端Fréchet微商零延遲互相關.因為Fréchet微商由震源端和檢波器端兩項組成，所以兩個參數(shù)結點對應Fréchet微商的互相關共包含四項，其中不同端的Fréchet微商在空間上幾乎不重疊，其互相關貢獻很小，故(24)式僅保留了相同端Fréchet微商的互相關，見圖2a和2b.如前文所講，擾動模型的更新是由嵌套在一起的最小二乘偏移來實現(xiàn)的.已有文獻討論這個過程中如何構建和應用海森算子(如Symes,2008;Tang,2009;Chen and Sacchi,2017)，因此僅在附錄B推導擾動模型的反射波敏感核和近似海森矩陣的計算式.

圖2 背景模型與擾動模型近似海森矩陣計算示意圖

在高頻近似意義下，基于射線追蹤的正演方法假定地震波僅沿著所謂的費馬射線傳播.相應地，射線層析利用走時數(shù)據(jù)更新射線路徑上的慢度參數(shù)，其中走時敏感核對應模型每個離散單元內(nèi)射線段的長度(Nolet，2008).波動方程走時層析或全波形反演則依據(jù)地震波能量沿著“波路徑”的傳播，揭示波路徑上模型參數(shù)對波場傳播與散射的作用(Woodward,1992;Tromp et al.,2005)，對應的波形敏感核沿著連接震源、檢波器的費馬射線向四周展開；敏感核的寬度可用來衡量走時層析、偏移成像以及波形反演的分辨率(Woodward,1992).當?shù)卣鹱硬l帶無限寬時，帶限波路徑就蛻變成費馬射線.

本文理論推導采用平方慢度模型參數(shù)化方式，但數(shù)值實驗展示更直觀的速度模型.如圖3，在均勻背景含有一個高斯異常體的模型底部設一反射界面.首先合成反射地震數(shù)據(jù)，然后按前文公式計算背景、擾動模型的泛函梯度和近似海森矩陣.為了方便觀察，圖中僅顯示單次觀測的泛函梯度(或波路徑)，而近似海森矩陣則考慮了所有觀測的貢獻.一方面，背景速度模型的波路徑是由大角度(接近180°)相遇的正向入射場(或散射場)與伴隨散射場(或入射場)零延遲互相關得到的，波路徑的寬度要遠小于其長度(如圖3b).基于反射波形數(shù)據(jù)反演背景速度模型時，沿射線方向分辨率最低，垂直射線方向最高.另一方面，擾動速度模型的波路徑是由小角度(小于90°)相遇的正向與伴隨入射場零延遲互相關得到的，對應于疊前偏移的“微笑型橢圓”，側重更新反射界面或模型短波長分量(圖3c).

圖3 簡單模型RWI泛函梯度與近似海森矩陣

如圖2所示，近似海森矩陣是由反射波敏感核的互相關得到的.上面二維模型參數(shù)結點為100×100，其海森矩陣的維數(shù)為10000×10000.如圖3d，背景速度模型的近似海森矩陣對角元素幅值最大，非對角元素幅值也不小，故而很稠密，這與相應的反射波敏感核(或波路徑)展布有一定寬度有關.相反，在圖3e中，除了左上角之外，擾動速度模型的近似海森矩陣幾乎是對角絕對占優(yōu)的，且沿對角方向幅值逐漸變小，這與幾何擴散效應有關.淺層擾動模型的敏感核涉及到大角度相遇的正向與反向入射場的相互作用(類似FWI中直達波和透射波對應的敏感核)，因此海森矩陣左上角較稠密.在射線走時層析中，海森算子對泛函梯度的預條件相當于對模型網(wǎng)格內(nèi)射線總長度(或照明)進行歸一化校正；而在本文反射波形反演中，背景模型的海森矩陣很大程度可消除孔徑限制、照明不均以及數(shù)據(jù)帶限對泛函梯度的模糊化，使模型更新方向更合理、更均衡，尤其是明顯提升垂向分辨率(見圖4).可見，海森矩陣與泛函梯度共同影響模型更新方向和反演分辨能力.

圖4 背景速度模型更新方向對比

2 數(shù)值實驗

這部分基于理論模型合成數(shù)據(jù)和實際地震資料檢驗本文波動方程反射波形反演二階優(yōu)化方法的有效性.針對反問題的病態(tài)性和零空間問題，我們借助構造傾角約束的正則化處理(Yu et al.,2020)壓制泛函梯度中不合理的高波數(shù)假象，規(guī)范速度模型的更新，提升反演結果的地質(zhì)一致性.

2.1 SEG/EAGE二維推覆體模型

本實驗從SEG/EAGE推覆體模型截取其主體部分，包含逆沖斷層、古河道等典型構造與沉積特征(圖5a).模型空間上有501×187個采樣點，采樣間隔為25×25 m，以主頻為15 Hz的雷克子波以100 m為間隔合成120炮道集作為觀測記錄，最大炮檢距為4 km，時間采樣為2 ms，記錄時長為3 s(見圖5b).為了檢驗RWI方法，將炮記錄初至波切除，僅保留反射波數(shù)據(jù).圖5c和5d為初始速度模型及相應的逆時偏移剖面.盡管初始速度長波長分量大致合理，但其誤差還是導致一些構造成像失真，分辨率較低.接下來從4 Hz到22 Hz分三個頻段(4～8 Hz、10～16 Hz和18～22 Hz)進行多尺度反演建模.

圖5 SEG/EAGE推覆體模型合成數(shù)據(jù)

圖6a顯示了RWI第一次迭代獲得的背景速度更新方向.受觀測孔徑、照明和數(shù)據(jù)帶限效應影響，負梯度方向在模型空間很不均衡，且受一些假象污染.如圖6b所示，高斯-牛頓方向引入了基于海森矩陣的去模糊化處理，提供了比較均衡、合理的更新方向(除了在模型底部兩側照明嚴重不足的區(qū)域).由于每一輪迭代都朝著比較可靠的方向前進，高斯-牛頓法最終重構出了更精確、含更多中波長成分的速度模型(如圖6c與6d).如圖7，共軛梯度法(CG)經(jīng)過100次迭代，誤差函數(shù)僅下降到40%左右，而截斷高斯-牛頓法(TGN)外循環(huán)20次(內(nèi)部循環(huán)最大迭代次數(shù)為5，統(tǒng)計平均將近4次)，誤差函數(shù)下降到了20%以下.把交替反演得到的背景速度和擾動速度疊加起來，重構出寬譜的層速度模型(圖8).連同偽井曲線(圖9)可以看出，高斯-牛頓RWI基本實現(xiàn)高分辨率速度建模目標.

圖6 背景速度模型反演

圖7 反射波形反演殘差下降曲線

圖8 共軛梯度法(左)與高斯-牛頓法(右)反演結果

圖9 偽井速度曲線對比

2.2 東海實際地震資料

接下來基于東海長江坳陷某二維拖纜地震資料檢驗方法實用性.該數(shù)據(jù)共有851炮，炮間距為37.5 m，最大炮檢距為4 km，檢波器間隔為12.5 m，記錄長度為4 s.用于實驗的地震數(shù)據(jù)已完成濾除涌浪噪聲、壓制海面多次波等常規(guī)處理(圖10).為了克服平面外傳播效應的影響，對所有炮記錄采用了三維到二維的振幅和相位校正(Wang and Rao,2009).初始層速度模型由常規(guī)疊前時間偏移處理獲得的均方根速度轉換而成.圖11顯示了基于初始模型的逆時偏移(RTM)剖面及沿測線均勻分布的8個共成像點道集(CIG).淺層CIG同相軸向上彎曲，表明速度偏小，而中深層CIG同相軸明顯向下彎曲，表明速度偏大，導致反射波沒有準確歸位，一些斷層繞射波沒有收斂.

圖10 典型共炮道集

圖11 基于初始速度模型的逆時偏移

作為對比，先用商業(yè)軟件基于疊前深度偏移的反射走時層析技術(Woodward et al.,2008)，從初始模型和相應的共成像點道集出發(fā)，經(jīng)數(shù)輪迭代構建層速度模型，并由逆時偏移得到新的成像結果(圖12).由于層速度模型得到修正，2 km以上平緩地層界面成像連續(xù)性有一定改善，坳陷內(nèi)部復雜構造、基底界面成像也有明顯改觀，CIG同相軸剩余曲率得到大幅度消減.鑒于射線走時層析的分辨率瓶頸，宏觀速度模型還不夠準確，還能觀察到CIG中、深層一些同相軸仍有一定的剩余時差.當初始模型偏離真實情況太遠時，反射波形反演仍然受到周波跳躍問題的困擾.為了從相同的初始模型出發(fā)穩(wěn)健地進行速度更新，我們對波動方程RWI算法略微做了調(diào)整，即在反演前期階段采用凸性較好的走時擬合標準(Ma and Hale,2013；Xu et al.,2019)，在穩(wěn)健地估計出速度模型長波長分量之后，再切換到波形擬合標準，進一步提取一些中波長分量并構建出更準確的層速度模型(圖13a).相應的逆時偏移效果由淺到深都得到了明顯提升，CIG同相軸基本都被拉平(圖13b和13c).由圖12b和13b及其局部放大(圖14)對比可知，本文RWI方法提高了速度建模精度，使逆時偏移除了使2 km以上的水平層狀地層界面更清晰，坳陷內(nèi)部和穿刺基底的一些斷層繞射波收斂也更徹底，從而高分辨率地刻畫了一些復雜的小斷塊，使深部復雜基底的解釋變得更容易.進一步地，圖15對比了基于共軛梯度RWI(CG-RWI)方法和高斯-牛頓RWI(GN-RWI)方法更新速度模型的逆時偏移結果.為了便于區(qū)分二者在模型更新上的差異，圖中偏移剖面疊合了總的速度更新量.可見，按本文二階優(yōu)化思路改善了速度更新，提高了地質(zhì)一致性，使反射波與繞射波歸位、聚焦效果進一步提升，更精細地刻畫了一些中深層的沉積序列和左側基底上覆的斷層結構.

圖12 商業(yè)軟件反射走時層析速度建?；A上的逆時偏移

圖13 本文反射波形反演速度建?；A上的逆時偏移

圖14 逆時偏移局部對比

圖15 不同反射波形反演算法基礎上的逆時偏移(左為速度更新量與偏移剖面疊合圖；右為共成像點道集)

3 討論

本文數(shù)值實驗只涉及聲壓場反射波數(shù)據(jù)縱波速度建模與偏移成像問題.若要逼近實際地下波場傳播的物理過程，有時需要考慮縱-橫波模式轉換與耦合、各向異性乃至吸收衰減等復雜波現(xiàn)象，因此要在建模與成像過程中引入橫波速度、密度、品質(zhì)因子以及各向異性等參數(shù)(如Tarantola,1986;石玉梅等,2010).近年來，隨著聲波方程單參數(shù)FWI、LSM和RWI技術走向工業(yè)應用，基于更復雜波動方程的多參數(shù)反演研究正逐年增多(如Choi et al.,2008; Operto et al.,2013; Alkhalifah and Plessix, 2014; Wang and Cheng,2017; Yang et al.,2016;王騰飛,2017;鄒鵬和程玖兵,2020).為了提高收斂效率、壓制參數(shù)耦合，考慮海森算子的多參數(shù)FWI、LSM方法受到了一些學者的關注(如Wang et al.,2016;Pan et al.,2017;Chen and Sacchi,2017).近來，徐文才(2020)在彈性波RWI理論框架下研究了縱、橫波速度和密度三參數(shù)情況下的海森矩陣，提出了彈性各向同性介質(zhì)反射波形反演高斯-牛頓方法.

本文方法原理和主要理論推導適合聲介質(zhì)也適合(黏)彈性、各向異性介質(zhì).一旦考慮更實際的波場物理和多參數(shù)問題，就要采用更復雜的波動方程正演引擎，泛函梯度、海森-向量積以及模型更新的計算成本會大幅攀升，這是多參數(shù)反射波形反演二階優(yōu)化方法走向實際應用必須應對的挑戰(zhàn).

4 結論

針對常規(guī)地面地震數(shù)據(jù)中、深層速度建模面臨的挑戰(zhàn)，本文在波動方程反射波形反演理論框架下，推導了背景與擾動模型的反射波敏感核、泛函梯度以及海森算子，建立了基于近似海森矩陣的反射波形反演高斯-牛頓方法.基于簡單模型揭示了兩個模型分量近似海森矩陣的特征與差異，證實了海森矩陣對泛函梯度的去模糊化作用可優(yōu)化速度模型更新方向.從SEG/EAGE推覆體模型合成數(shù)據(jù)實驗結果可知，相比于常用的共軛梯度法，利用海森矩陣的高斯-牛頓法提升了反射波形反演的收斂性與寬譜建模能力.東海長江坳陷拖纜地震數(shù)據(jù)處理實例表明，商業(yè)軟件射線走時層析技術的分辨率瓶頸制約了波動方程逆時偏移發(fā)揮其潛力；而本文二階優(yōu)化反射波形反演方法改善了速度長波長與中波長分量的重構，提高了偏移速度模型的可靠性，從而發(fā)揮出逆時偏移的成像分辨率優(yōu)勢，精細刻畫了坳陷內(nèi)部復雜斷裂系統(tǒng)以及深部基底形態(tài).本文方法在處理三維實際地震資料和多參數(shù)問題時，還面臨一些計算層面的挑戰(zhàn)，今后可結合數(shù)據(jù)壓縮、下采樣與隨機優(yōu)化等技術提升它的實用化水平.

致謝本文實例是依托國家科技重大專項子課題“長江坳陷油氣資源潛力評價”研究靶區(qū)的地震數(shù)據(jù)完成的.感謝中海石油有限公司上海分公司提供的支持.

附錄A 背景模型的敏感核函數(shù)和泛函梯度

對方程(3a)和(3b)關于背景模型求一階偏導數(shù)，得到

(A1)

和

(A2)

聯(lián)合(A1)和(A2)式，可得到背景模型對應的Fréchet導數(shù)(矩陣)：

(A3)

該式可寫成：

Jm0=L-1F0,

(A4)

其中F0是維數(shù)為q×l的虛源矩陣，其列向量滿足：

(A5)

因此，就背景模型某個特定網(wǎng)格結點，F(xiàn)réchet導數(shù)可按下式計算：

假定在檢波器r處設一點脈沖，則有背景格林函數(shù)G0=L-1δ(x-r)以及對應的擾動格林函數(shù)LG1=QG0.于是，與檢波器單次觀測對應的Fréchet導數(shù)可表示成：

(A7)

把(A7)代入(8)式，則可構建背景模型對應的泛函梯度：

由于L-1每一列與每個參數(shù)結點處脈沖源的格林函數(shù)波場相對應，且根據(jù)格林函數(shù)的互易性，有(L-1)?=L-1(Virieux and Operto,2009)，于是(A8)式化簡得到

(A9)

(A10a)

(A10b)

附錄B 擾動模型的敏感核函數(shù)、泛函梯度與近似海森矩陣

由(12)與(13)式可知，擾動模型對應的Fréchet導數(shù)滿足：

(B1)

引入與某檢波器r處點脈沖對應的背景格林函數(shù)G0，則有：

(B2)

(B2)代入(8)式，則把擾動模型對應的泛函梯度表示成：

(B3)

同樣，利用格林函數(shù)互易性，(B3)式進一步可寫成：

(B4)

(B5)

給定震源s和檢波器r的單次觀測，擾動速度模型的反射波敏感核按(B2)式可表示成：

Jm1(x,r,s,ω)=ω2f(ω)[G0(x,r,ω)G0(x,s,ω)].

(B6)

它僅涉及入射場的乘積(時間域褶積)，揭示了一次散射對擾動速度變化的敏感性.注意，敏感核是頻率或時間的函數(shù)，只有與數(shù)據(jù)殘差相結合才體現(xiàn)相應數(shù)據(jù)成分對模型更新的作用.

根據(jù)(20)式，對應的近似海森矩陣的元素可按如下公式計算：

×G0(x,r,ω)]*×[G0(x′,s,ω)G0(x′,r,ω)].

(B7)

可見，擾動模型的近似海森矩陣是相應Fréchet微商的互相關構建出來的(圖2c).注意，上式已在前人文獻中給出(如Tang,2009;任浩然等,2013).