基于垂直構造梯度優化的全波形反演方法

張靈莉， 任玉曉， 劉斌,3*， 王凱， 許新驥， 陳磊

1 山東大學巖土與結構工程中心， 濟南 250061 2 山東大學齊魯交通學院， 濟南 250002 3 山東大學數據科學研究院， 濟南 250100 4 山東高速集團有限公司， 濟南 250102

0 引言

自Lailly(1983)和Tarantola(1984,1986)等學者于20世紀80年代提出基于最小二乘法的全波形反演(Full Waveform Inversion, FWI)方法以來，全波形反演技術以其速度建模精度高、地下介質重構分辨率高和復雜構造適應性強等優勢成為勘探地球物理領域的重要研究方向之一.全波形反演技術受計算能力影響較大(Pratt,1990; Bunks et al.,1995)，近年來，得益于計算機技術的飛速提升，全波形反演技術獲得了迅速的發展(Prieux et al.,2011,2013; Sirgue et al.,2010; 卞愛飛等,2010; 董良國等,2013; 楊積忠等,2014; 孫史磊等,2020; 潘冬雪等,2021).

全波形反演基于波動方程反演框架，能夠充分利用地震波的動力學和運動學信息.以聲波為例，全波形反演問題可以簡化為實際觀測數據和模擬數據之間的誤差極小值求解問題，即以反演目標函數(Guitton and Symes,2003; Liu et al.,2017)為收斂判斷依據，從初始速度模型開始，計算梯度方向并應用優化理論優化梯度，隨后計算相應步長，不斷更新速度模型最終獲得誤差極小值，該極小值對應的速度模型即為最終的反演結果(Tarantola et al.,1987).梯度計算及優化是影響全波形反演收斂效果的主要因素之一.在全波形反演梯度計算方面，伴隨狀態法是主要方法之一(Plessix,2006; Virieux and Operto,2009).伴隨狀態法求解目標函數梯度方向類似于逆時偏移成像過程，首先正演模擬獲得震源正傳波場的時間二階導項和觀測地震記錄，然后以觀測地震記錄與實際記錄的殘差為震源逆傳獲得殘差逆傳波場，震源正傳波場的時間二階導項和殘差逆傳波場的互相關即為近似梯度方向.在全波形反演梯度優化方面，相關研究主要集中在包括共軛梯度法、牛頓法、擬牛頓法等方法在內的局部優化方法(Shin et al.,2006; Dagnino et al.,2014; 王義和董良國,2015).

高陡構造普遍存在于油氣藏勘探和金屬礦勘察中，具有探測難度大、定位精度低的特點，制約了高分辨率地震勘探的發展(Hale et al.,1992; 曲英銘等,2020).為解決高陡構造逆時偏移成像效果不佳的問題，Singh和Curtis(2019)提出了一種新的檢波波場自相關成像條件，改善了以垂直構造為代表的高陡構造的成像效果.研究過程中，本文發現全波形反演方法盡管可以實現水平構造的較準確反演，但是對垂直構造反演效果較差.以單層階梯模型為例(圖1a)，初始速度模型(圖1b)由真實速度模型經100的平滑因子平滑得到(100的平滑因子本質上是真實模型以100×100的模板進行均值濾波，本文采取Matlab中的fspecial函數具體實現平滑過程).反演結果(圖1c)與逆時偏移成像結果類似，除獲得了階梯構造上下界面信息以外，反演結果未體現垂直構造信息.

圖1 (a) 單層階梯模型； (b) 初始速度模型； (c) 應用伴隨狀態法的全波形反演結果

受Singh和Curtis(2019)啟發，針對垂直構造全波形反演效果不佳的問題，本文提出了一種基于垂直構造梯度優化的全波形反演方法.該方法通過構建模擬數據與觀測數據之間殘差逆傳波場的自相關項，實現了垂直構造位置處的梯度補償；并將其以自相關梯度項能量最大值為歸一準則與伴隨狀態法所得梯度相結合，在改善垂直構造反演效果的同時提升了水平構造的反演效果.本文首先從波場傳播角度分析了伴隨狀態法在垂直構造梯度計算方面的局限性.隨后，介紹了基于垂直構造的梯度優化方法，給出了基于殘差逆傳波場自相關項補償梯度的改進梯度計算公式、權重因子及改進全波形反演方法流程.通過單層階梯模型、雙層階梯模型和局部Marmousi模型測試了本文方法的有效性和實用性.最后給出了討論和結論.

1 改進全波形反演梯度計算方法

1.1 傳統伴隨狀態法局限性分析

全波形反演方法通過最小化模擬數據與觀測數據的殘差反演獲得地下模型參數.全波形反演一般采取伴隨狀態法(Tarantola,1984)計算梯度，以避免計算Fréchet導數的巨大代價.本文研究基于二維常密度聲波方程，相應的伴隨狀態法梯度計算公式可表示為：

(1)

其中，t、U、Ub、g、xs分別表示時間、正傳波場、殘差逆傳波場、梯度和震源；殘差逆傳波場Ub由模擬數據和觀測數據的殘差逆時傳播獲得，偽保守形式下，Ub滿足：

(2)

其中x、z分別表示水平方向變量、垂直方向變量；v表示速度模型；Serr為模擬數據與觀測數據的殘差，可以表示為：

Serr=RU-dobs,

(3)

其中R代表檢波器序列，RU代表檢波器接收到的地震記錄，dobs為觀測到的地震記錄.式(1)同時包含正傳波場對時間的二階導數項和逆傳波場項，為統一格式以便簡化下文關于伴隨狀態法在計算垂直構造方面局限性的理論闡述，本文對式(1)定義的梯度方向進行了一定的近似(胡光輝等,2014)，省去了正傳波場對時間的二階導數(實際伴隨狀態法梯度計算時仍采用公式(1))，此時，可得到一個形似逆時偏移的梯度解：

g=∑xs∑tUUb.

(4)

由式(4)可知，模型空間內任意(x,z)位置處梯度不為零的充分必要條件是存在某一時刻，該位置處的震源正傳波場場值與殘差逆傳波場場值均不為零.當正演模型與真實模型參數完全一致時，模擬數據與觀測數據之差方為零，此時，所對應的殘差逆傳波場值也相應為零.正演模型與真實模型稍有偏差，殘差逆傳波場均不為零.受限于反演本身，反演結果無法與真實模型完全一致，因此，只要反射波可以被檢波器接收到，模擬數據與觀測數據之差便不為零，即殘差逆傳波場不為零.綜上所述，轉換到波場幾何空間中，基于伴隨狀態法，模型空間中某處梯度值不為零的充分必要條件是該位置處產生的反射波可以被檢波器接收到且與震源激發入射波在該位置幾何相交.

本文從波場傳播角度分析了伴隨狀態法的局限性，為具體分析正傳及逆傳過程中的波場傳播情況，本文再次引入了圖1a所示單層階梯模型.為了便于分析殘差逆傳波場，初始模型設置為均一模型，即模擬數據中不包含反射波，因此，模擬數據與觀測數據的殘差可認為僅由單一觀測記錄組成，殘差可以表示為：

Serr=-dobs.

(5)

此外，為了簡化波傳狀態以便于分析波場(使波通過介質分界面時依然按照直線路徑傳播)，本文認為單層階梯上下介質僅存在密度變化，波速保持相同.本文方法同樣適用于變速度模型.

根據震源與垂直構造的相對位置關系，波場傳播情況可分為以下三種情況.(1)當震源位于階梯拐角正上方位置時，如圖2a所示；當激發地震波沿sp1向左傳播時，遇上層水平界面產生可以被檢波器R1接收到的一次反射波rp1，R1逆傳時rp1路徑與正傳sp1路徑在上層界面存在交點，即對水平構造反演有貢獻、對垂直構造無貢獻，本文將該情況命名為一次反射波梯度計算效應.當地震波沿sp2向右傳播時滿足一次反射波梯度計算效應條件.當地震波沿sp3到達階梯拐角位置時，將產生發散狀繞射波(rp3-rpn所示)，部分可以被檢波器接收到的上行繞射波逆傳時與激發地震波sp2在拐角位置相交，所以拐角處梯度不為零，且由于疊加效應，此位置梯度值將相對較大，本文將此種現象命名為繞射波梯度計算效應，此效應等同于對垂直構造反演無貢獻.(2)當震源位于垂直構造左側時(圖2b)，激發地震波除傳播至階梯拐角位置和上、下層水平界面以外(對垂直構造無貢獻)，還可能沿sp1傳播至垂直界面產生下行反射波rp1(無法被檢波器接收)和繼續傳播的下行透射波sp1，此時，逆傳波場為零，仍然對垂直構造反演無貢獻.(3)當震源位于垂直構造右側時(圖2c)，激發地震波除傳播至拐角和上層水平界面(對垂直構造無貢獻)，還可能沿sp1直接傳播至垂直界面產生一次反射波rp1，rp1遇下層水平界面產生二次反射波rp1′，sp1與rp1′平行無交點(反射定律)，此時無法反演垂直構造.或沿sp2傳播至下層水平界面產生一次反射波rp2直至被檢波器R2接收，對下層水平構造產生貢獻.同時rp2遇垂直界面將在五角星處產生二次反射波rp2′被檢波器R3接收，rp2′和sp2平行無交點，對反演無貢獻.綜上所述，當震源位于階梯拐角正上方、垂直構造左側和垂直構造右側時，應用基于伴隨狀態法的全波形反演方法均無法有效反演垂直構造.

圖2 震源在單層階梯模型階梯構造(a)拐角正上方、 (b)垂直構造左側、(c) 垂直構造右側的波傳示意圖

1.2 改進梯度計算方法原理

Singh和Curtis(2019)針對垂直(近垂直)構造，提出了基于檢波波場自相關的新型逆時偏移成像條件，可以實現垂直(近垂直)構造位置處的有效成像.由本文1.1節內容可知，傳統伴隨狀態法在計算垂直(近垂直)構造時存在一定的局限性，基于逆時偏移成像條件和全波形反演梯度更新項的相似性，檢波波場自相關項為全波形反演垂直(近垂直)位置處梯度的修正提供了一種新思路.本節將重點討論解決伴隨狀態法在垂直(近垂直)構造位置處梯度計算的局限性問題.

由圖2c可知，盡管入射波sp1與遇垂直構造產生的反射波rp1無法在垂直構造處幾何相交，但是可以被檢波器R2接收到的一次反射波rp2和可以被檢波器R3接收到的二次反射波rp2′在垂直構造五角星標注位置相交，因此，本文考慮利用一次反射波與二次反射波的這種時空一致性特征來反演垂直構造.據此提出了通過殘差逆傳波場自相關來反映垂直構造位置的數值梯度，表達式為：

(6)

其中，gauto為通過殘差逆傳波場自相關構造的數值梯度，然而，由于是殘差逆傳波場的自相關，gauto將恒為正值.為解決這一問題，考慮在公式(6)的右側增加梯度方向控制因子γ，梯度gdir進一步表示為：

(7)

胡光輝等(2014)指出在保證全波形反演梯度方向大致正確(不至于出現符號反轉)的前提下，可以對梯度作相應的改變即形成方式上作一些近似.因此，本文研究基于以下基本假設：傳統伴隨狀態法計算所得梯度方向是正確的，僅需在數值上做一定的修正以解決伴隨狀態法在垂直(近垂直)構造位置處梯度計算的局限性.基于以上假設，方向控制因子應與傳統伴隨狀態法所得梯度的方向相一致，方向控制因子可以表示為：

(8)

其中，symol運算符為本文定義的取符號算子.

圖1a所示的單層階梯模型再次被引入以測試上文所述梯度修正思路的正確性.單層階梯模型上層介質和下層介質波速分別設置為2500 m·s-1和3000 m·s-1，中心頻率、網格間距、時間步長、時窗分別設置為40 Hz、5 m、0.01 ms以及0.35 s；震源設置在X=250 m，Z=5 m處，100個檢波器均布在Z=5 m處.為與圖2所示情景相對應，本節采用均一速度模型為反演初始模型.以公式(4)、公式(6)、公式(7)為梯度計算準則的第一次迭代梯度結果分別如圖3a、b、c所示.讀圖可知，傳統伴隨狀態法所得梯度(圖3a)在單層階梯上、下界面梯度值較大、垂直構造位置處的梯度值基本為零；殘差逆傳波場自相關所得梯度(圖3b)可以有效反映垂直構造位置處存在波速變化，但恒為正值；通過引入方向控制因子，梯度值(圖3c)得到了修正，獲得了同時可以反映水平構造和垂直構造的梯度值，這說明以傳統伴隨狀態法所得互相關梯度為方向主導、以殘差逆傳波場自相關為數值主導的梯度計算思路是正確的.

觀察圖3c可得，式(7)所得梯度顯著補償了單層階梯模型垂直構造位置處的梯度，但同時也犧牲了部分水平構造處的梯度值.為此，考慮結合互相關梯度數值與自相關梯度數值，以在保證水平構造反演效果的同時，提升垂直構造反演效果.梯度可更新為：

+β*γ*Ub(x,z,t)*Ub(x,z,t)),

(9)

其中，gimp為本文最終優化所得全波形反演梯度，α和β分別為互相關項加權因子和自相關項加權因子.由于上文證實殘差逆傳波場自相關項為數值主導的梯度計算思路是正確的，因此傳統伴隨狀態法所得互相關梯度數值應以自相關梯度項能量最大值為歸一準則經能量級匹配之后與之相加，即α、β按如下方式取值：

圖3d為經公式(10)計算所得梯度示意圖，圖3c中垂直構造梯度得以保留的同時，水平構造位置處的梯度得到了有效增強.需要說明的是，圖3d中左下方及右上方的能量弧可經多炮疊加后得到有效壓制.

綜上所述，以激發地震波遇垂直界面產生的一次反射波和二次反射波的波傳時空一致性為突破口，以實現垂直構造有效反演兼具水平構造反演效果為目標，基于多次波信息，本文提出了一種新的基于殘差逆傳波場自相關項補償垂直構造位置梯度的改進全波形反演梯度計算方法，進而形成一套基于垂直構造梯度優化的全波形反演方法，該方法將伴隨狀態法所求梯度(正傳波場的時間二階導項和殘差逆傳波場互相關)和殘差逆傳波場自相關項以自相關梯度項能量最大值為準則進行歸一加權.

1.3 基于垂直構造梯度優化的全波形反演方法實現流程

本文提出的基于垂直構造梯度優化的全波形反演方法具體實現流程如下(圖4):

圖4 基于垂直構造梯度優化的全波形反演方法流程

(1)反演初始參數設置，初始模型輸入；

(2)計算震源正傳波場，獲得模擬地震記錄；

(3)依據目標函數，計算模擬地震記錄與實際地震記錄間的的殘差；

(4)將步驟(3)所得殘差作為震源逆時傳播獲得殘差逆傳波場；

(5)根據公式(9)計算梯度方向；

(6)利用局部優化理論優化梯度方向，獲得最終模型更新方向；

(7)基于步驟(6)所得模型更新方向，依據步長計算準則，計算反演收斂步長；

(8)更新模型；

(9)根據目標函數或最大迭代次數判斷是否繼續迭代，若停止迭代，步驟(8)所得模型即是最終反演模型；否則繼續迭代，重復步驟(2)～(9)，直至迭代結束，獲得最終反演模型.

2 數值算例

為測試本文提出的全波形反演方法在垂直構造、水平構造以及傾斜構造反演方面的能力，本文引入了單層、雙層階梯模型和局部Marmousi模型進行了算例試算.本文相關算例測試均基于聲波全波形反演、L2范數目標函數、共軛梯度法優化和拋物線步長計算方法等應用相對廣泛的全波形反演理論.同時，考慮到低頻信息對全波形反演的重要性，為保證反演效果，本文相關研究均基于利用維納濾波的多尺度全波形反演方法(張文生等, 2015).

2.1 單層階梯模型

為了初步測試本文提出的全波形反演方法應對垂直構造的反演能力，本文再次引入了圖1a所示的單層階梯模型.相關參數與1.2節所述相同，此外，為增加反演的低頻信息以保證反演效果，本文選擇了包括5 Hz、15 Hz、30 Hz、45 Hz、60 Hz五個頻率在內的多尺度策略，每個頻率段反演迭代十次.共設置50個炮點.反演初始模型如圖1b所示.

基于伴隨狀態法的單層階梯模型全波形反演結果如圖1c所示，應用本文方法(公式(9)) 的反演結果如圖5b 所示.除在水平界面位置出現類似于成像界面的能量弧以外，圖1c基本保留了初始模型特征，多次迭代后仍與真實模型偏差較大.圖5a經反演之后基本恢復了單層階梯模型的結構特征，上下層介質波速均已趨近于真實模型波速.就垂直構造反演效果而言，圖1c中未見明顯垂直構造特征；僅在垂直構造頂端點和底端點位置出現輕微波速擾動現象，這是由于1.1節所述“繞射波梯度計算效應”造成的.相比圖1c，圖5a中垂直構造反演的效果已經有效提升，垂直構造已經初現輪廓；垂直構造左側區域波速已與實際模型波速較為接近；就水平構造而言，圖1c僅在上下界面位置處出現兩處水平狀能量弧，能量弧上下方介質波速與初始模型波速基本保持一致,且能量弧僅出現在模型中間位置.相比較圖1c僅出現能量弧且能量弧不連續的情況，圖5a中階梯構造上下界面已完整呈現，上下介質波速均已趨近于真實模型；即使是反射信息不能被完全接收的水平構造兩端，也已獲得了較好的反演效果.單層階梯模型的反演效果對比初步證明了本文所提全波形反演方法的有效性.為了更加定量的驗證本文提出方法的有效性，本節引入了應用傳統伴隨狀態法和本文方法的目標函數對比圖(圖6，目標函數為常規最小二乘目標函數.為使二者區別更為明顯，應用傳統伴隨狀態法和本文方法的目標函數均以前者的最大值為基準做歸一化).讀圖可知，本文方法的應用顯著提升了全波形反演的收斂速度，且擬合誤差在反演全過程中均低于傳統方法.

圖6 目標函數值隨迭代次數變化圖(單層階梯模型)

同時，本節還給出了依據公式(7)計算梯度所得的反演結果(圖5b)，以期支撐上文1.2中所得結論.讀圖可知，盡管圖5b相較圖1c的反演質量得到了明顯提升，但是較圖5a的反演結果仍存在一定的差距：首先，對于上下水平界面來說，圖5a更為清晰、界限感更強；其次，相較于圖5a，圖5b反演所得模型速度較不均勻，這一點尤其體現在下層介質；另外，在垂直構造位置處，圖5a的收斂效果要明顯優于圖5b，尤其體現在垂直構造的中下部.本算例再次證實了1.2節所述最終梯度計算方法(公式(9))是對公式(7)的進一步完善.

圖5 單層階梯模型應用(a)公式(9)和(b)公式(7)計算梯度所得反演結果

2.2 雙層階梯模型

為了測試本文提出的方法應對稍復雜模型的反演能力，本文引入了圖7a所示的雙層階梯模型，該模型共有三種介質組成，從上至下波速分別設置為：2500 m·s-1、3000 m·s-1和3500 m·s-1.初始速度模型(圖7b)由真實速度模型以100的平滑因子平滑得到，其余反演參數設置與單層階梯模型相同.

基于伴隨狀態法的全波形反演結果和應用本文方法的反演結果分別如圖7c、d所示.盡管圖7c中的上層介質已出現向真實模型收斂的趨勢，但仍與真實模型差別較大；除水平構造位置處出現水平能量弧外，中層介質和下層介質均基本保留了初始模型的速度特征.圖7d中上層介質基本趨同于真實模型波速，僅出現波速輕微反演偏差；中層介質和下層介質均已跳脫出初始速度模型，出現向真實模型收斂的趨勢.就垂直構造反演效果而言，相較圖7c，圖7d中垂直構造的反演效果已經明顯提升.對于上層階梯的垂直構造來說，圖7c僅出現了由于繞射波造成的“偽垂直構造”現象，這是由于上層階梯尺度較小，文中1.1節中所述“繞射波梯度計算效應”在小尺度構造上對比明顯造成的；相比之下，圖7d中上層階梯垂直構造輪廓已以較高分辨率完整呈現.對于下層階梯垂直構造來說，圖7c未出現垂直構造響應特征，其兩側區域均保留了初始模型波速特征，可以認為下層階梯垂直構造未收斂；而圖7d中，下層垂直構造盡管沒有完全收斂，但其相較初始速度模型已發生明顯變化，出現向真實速度模型收斂的趨勢；對比圖7d上、下層階梯垂直構造反演效果可知，相較上層階梯垂直構造的完整呈現，下層階梯垂直構造反演僅出現收斂趨勢未能完全收斂，這是因為在偏移距固定的情況下，可以接收到的反射波信息隨著垂直構造深度的增加而逐漸減少.從水平構造來看，相比圖7c，本文所提全波形反演方法的應用極大提升了水平界面(包含上下階梯構造共計四個水平界面)的反演效果，主要體現在界面感增強、分辨率提高、邊緣反演效果提升三方面.雙層階梯模型的反演效果對比證明了本文方法具有應對稍復雜模型的能力.

圖7 (a) 雙層階梯模型； (b) 反演初始模型； (c) 伴隨狀態法反演結果； (d)本文方法反演結果

圖7c和圖7d的對比可以直觀看出本文方法相較傳統方法可以明顯提升反演效果，為定量對比本文方法與傳統方法效果，這里同樣給出了利用傳統方法與本文方法反演雙層階梯模型時目標函數隨迭代次數的變化趨勢圖(圖8，歸一化方式與單層階梯模型一致).與單層階梯模型類似，應用本文方法的目標函數值在反演全過程均小于傳統方法，且收斂速度更快.通過對目標函數值(反演誤差)的定量分析，進一步驗證了本文方法的有效性.

圖8 目標函數值隨迭代次數變化圖(雙層階梯模型)

2.3 局部Marmousi模型

為了測試本文方法對復雜傾斜模型的適用性，本文引入了如圖9a所示的局部Marmousi模型(該模型經原始Marmousi模型抽道、截取所得)，采取真實模型波速最小值到最大值為波速上下限的漸變速度模型為初始模型(圖9c).在反演過程中，網格間距設置為4 m，時間步長設置為0.25 ms，時窗設置為0.35 s，中心頻率為30 Hz；為保證反演效果，本文依然采取多尺度反演策略，頻率分別設置為5 Hz、15 Hz、30 Hz、45 Hz、60 Hz，每個頻率的反演迭代次數設置為10.在地面位置均布70個炮點和140個檢波器；本算例中未額外拓展檢波器和炮點，涉及所有計算均在本文所述空間內，因此，反演結果未實現較深部區域的反演，后續可以通過拓寬模型(X軸所示方向)以實現較深部的反演，本文僅涉及理論初步探討，在此不再贅述.

圖9 (a) 局部Marmousi模型； (b) 圖(a)矩形框選區域放大圖； (c) 初始速度模型； (d) 圖(c)矩形框選區域； (e) 伴隨狀態法反演結果； (f) 圖(c)矩形框選區域放大圖； (g) 本文方法反演結果； (h) 圖(e)矩形框選區域放大圖. i、ii、iii標注了模型空間內的三個主要斷層; a,b標注了重點分析的兩個界面, A和B分別為a、b兩個斷層在實線框區域內的起始點

圖9e為基于伴隨狀態法全波形反演方法的反演結果，圖9g為應用本文方法的反演結果，圖9b、d、f、h分別為圖9a、c、e、g中白色實線框所圈區域的局部放大圖.相比較圖9e，圖9g的反演結果整體均有所提升.本文重點分析圖9a中由“i”、“ii”和“iii”標注的三條主斷層的反演效果.讀圖可知，不論是圖9e還是圖9g均實現了斷層“i”的較好反演，這是因為該斷層左側和右側均由小傾角層狀界面堆疊而成，斷層“i”的有效反演等同于其兩側緩傾角構造的較準確反演.而基于伴隨狀態法的全波形反演方法和本文方法均可以實現水平構造(小傾角構造)的較精確反演，因此兩種方法對斷層“i”的反演效果不分伯仲.斷層”ii”左側部分(圖9b、d、f、h中由實線“b”標記)存在相對低速-相對高速(記為高速1)-相對低速-相對高速(記為高速2)-相對低速的大致速度變化；由于在b方向上發育較窄，圖9f中基本未體現高速1，應用本文方法后，高速1在圖9h中得以顯現；同時，對于高速2，相比圖9f，圖9h的反演結果也更趨近于真實高速2的分布范圍；得益于高速1和高速2的較準確反演，本文方法對斷層2的反演效果優于傳統方法.受限于檢波器分布，斷層“iii”(圖9b、d、f、h中由“a”大致標記)上半部構造未得到有效反演，但是，相較圖9f接近勻質的反演結果，圖9h已出現向斷層“iii”收斂的趨勢；就斷層“iii”下部高速體反演效果而言，圖9h的反演結果也更接近于真實模型.局部Marmousi模型的反演結果證明了本文所提方法不僅適用于垂直構造，同樣適用于傾斜構造.

為更加清晰的展示本文方法對局部Marmousi模型斷層構造反演的提升作用，沿圖9中a、b線抽取了其速度值，并以沿斷層距起始點的距離和速度作為橫、縱坐標，分別作圖如圖10a、b所示.讀圖10可知，不管是在a線位置還是在b線位置，相比較由點劃線標記的傳統方法速度反演結果，由短劃線所示的速度反演結果(應用本文方法)在大多數區域更接近于相應位置處的真實速度分布情況.

圖10 (a)沿圖9中 a線抽道和(b)沿b線抽道的速度擬合圖

相比較2.1節單層階梯模型、2.2節雙層階梯模型應用本文方法反演效果的顯著提升而言，本節局部Marmousi模型的效果提升相對較弱，這主要是由于Marmousi本身構造特性造成的.如上文所述，Marmousi模型斷層由兩側緩傾角構造密集堆疊而成，因此，斷層兩側的緩傾角構造可以有效反演，即等同于Marmousi模型的斷層構造可以有效反演，本文稱之為斷層構造的“被動反演”.當采用全Marmousi模型時，由于偏移距足夠大，斷層兩側的緩傾角構造反射信息可以被分布范圍較廣的檢波器獲得，在利用傳統伴隨狀態法時即可以得到有效反演斷層兩側的緩傾角構造，此時，斷層實現“被動反演”.因此，本節未采用全Marmousi模型，以避免斷層構造的“被動反演”掩蓋了本文方法應對垂直(近垂直)構造的反演能力.

單層階梯模型和雙層階梯模型的算例對比證明了本文提出的全波形反演方法在垂直構造反演、水平構造反演和整體速度反演方面均具有較大優勢.相較于基于伴隨狀態法的全波形反演方法垂直構造無法反演、水平構造邊界收斂不佳、速度反演不勻質的特點，本方法可以較準確刻畫垂直構造、較完整呈現水平構造、較準確重現速度；證明了本文方法的有效性.局部Marmousi算例結果表明了本方法在應對傾斜構造時同樣可以提升反演效果，進一步拓寬了本文應用范圍，證明了本方法的實用性和普適性.

3 討論

本文所提方法利用殘差逆傳波場自相關以檢波波場為準則歸一化補償伴隨狀態法所得梯度，在一定程度改善垂直構造反演效果的同時，提升了水平構造的反演效果.局部Marmousi模型的反演結果證明本文所提方法同樣適用于傾斜構造.本方法可以提升全波形反演效果，但仍存在不足和深入研究空間，在此作如下討論：

(1) 本文所提方法中的殘差逆傳波場自相關實際上是垂直構造處一次反射波及二次反射波的互相關，推廣可至不同階次多次波的互相關.上述多次波信息即為垂直構造反演的有效信息，因此垂直構造的反演效果將隨多次波信息的增加而有所提升.由單層階梯模型反演結果(見圖1c、圖5a)和雙層階梯模型反演結果對比(見圖7c、d)可知：雙層階梯模型上部垂直構造反演效果明顯優于下層垂直構造和單層階梯模型垂直構造反演效果.這是因為雙層階梯模型上層階梯垂直構造處除存在由上層階梯下界面產生的多次波還存在由下層階梯下界面產生的多次波，該界面有效多次波信息相比單層階梯模型垂直構造和雙層階梯模型下階梯垂直構造明顯增多，因此雙層階梯模型上層階梯垂直構造反演效果較好.此外，多次反射波的能量級可能隨反射次數的增加而出現不在相同量級的情況，而本文僅實現了多次反射波的直接自相關，尚未考慮多次反射波的能量級是否匹配問題.因此，在后續全波形反演研究中，可以從不同階數多次波與反演效果關系研究、不同階多次波能量級匹配及相應梯度優化、及含多次波信息的目標函數建立三方面進一步完善本方法.

(2) 由圖1c和圖7c可知，基于伴隨狀態法的反演結果在水平界面邊緣位置效果較差，這是由于偏移距不足造成的.對比圖5a、圖7d可知，應用本文所提方法，除在垂直構造位置處反演效果得到明顯提升以外，水平構造的邊緣位置反演結果界面更為清晰、分辨率也相對增高.綜上所述，本文認為基于垂直構造梯度優化的全波形反演方法對偏移距要求更低、對小偏移距適應性更強.基于單層階梯模型，本節就偏移距對反演效果的影響進行了初步測試.圖11a、b、c，圖1a、b、c、圖5a、b，圖11d、e、f分別代表模型橫向尺寸為400 m、500 m、600 m的相關結果.震源和檢波器的間隔分別設置為10 m和5 m，其余參數與2.1節保持一致.通過對比不同偏移距的反演結果，本節可以得出如下結論：

圖11 小偏移距時 (a) 單層階梯模型反演初始模型; (b) 伴隨狀態法反演結果; (c) 本文方法反演結果及大偏移距時; (d) 單層階梯模型反演初始模型; (e) 伴隨狀態法反演結果; (f) 本文方法反演結果

①不管偏移距如何變化，相同偏移距情況，應用本文方法均對垂直構造有響應，且反演結果均優于傳統方法.

②隨著偏移距的增加，應用本文方法的反演效果得以進一步提升.提升效果可以體現在垂直構造收斂效果更好、介質內部波速反演更為均勻、更為準確等多個方面.

③隨著偏移距的增加，傳統方法的反演效果也得以提升.對于水平構造反演來說，偏移距的增加會使得有效反射信息變多，從而增強反演效果；對于垂直構造反演來說，水平構造的有效反演一定程度上會增強垂直構造的反演視覺效果.此外，如果偏移距足夠長的話，入射波遇垂直構造產生的反射淺水波有可能被檢波器接收到，從而增強垂直構造位置處的反演效果.

本文僅從偏移距與反演效果方面進行了初步的測試和分析，未來，可以從偏移距與反演效果評價方面繼續改進伴隨狀態所得梯度補償項及加權因子，以進一步提升方法對小偏移距數據的反演效果.在此基礎上，可以進一步將本文方法應用于現場實際數據中(在實際探測中，由于現場條件所限，實際數據往往偏移距較小)，以進一步測試本文所提方法小偏移距反演的能力.

(3) 式(9)中α和β的取值會影響最終的反演結果，這一點由圖1c、圖5a和圖5b可見一斑，上述三圖對應的α和β取值可以分別表示為式(11)、式(10)和式(12).式(11)和式(12)為：

α=1,β=0，

(11)

α=0,β=1.

(12)

為探討加權因子對反演效果的影響，本文另外測試了常規標準歸一化加權因子和以互相關梯度項能量最大值為歸一準則的加權因子的反演效果.常規標準歸一化加權因子的具體實現方法為，將互相關項除以互相關項絕對值的最大值以實現互相關項的標準歸一化，將自相關項除以自相關項絕對值的最大值以實現自相關項的標準歸一化，此時α和β的取值可以表示為公式(13):

(13)

當以互相關梯度項為歸一化準則進行非常規歸一化加權因子計算時, 互相關梯度項的數量級為最終計算梯度數量級；因此，互相關梯度項加權因子為1，自相關梯度項歸一至以互相關梯度項最大值為最大值，表達式為：

圖1a所示的單層階梯模型和圖1b所示的初始速度模型再次被引入以測試上述加權因子的反演效果，相關反演參數設置與上文引言部分及2.1節所述參數一致，常規標準歸一化加權因子和以互相關梯度項能量最大值為歸一準則的加權因子反演效果分別如圖12a和圖12b所示，圖12a所示反演結果盡管已經顯示了階梯構造的輪廓，但是卻在上部出現了較強的干擾噪聲；相比圖1c所示的傳統梯度計算反演結果，圖12b反演效果得到明顯提升，但是效果卻不敵圖5a和圖5b，對此我們做分析如下：

(1)圖1c和圖12b均為互相關能量占優，圖1c僅包含互相關能量，圖12b除包含互相關能量外，還包含少部分的自相關能量.相較圖1c僅在水平位置出現兩條水平能量弧、其上下方介質波速未見明顯收斂的情況，圖12b水平界面上下方介質的速度得到了明顯收斂.這表明在互相關能量占優的情況下，自相關能量的加入可以明顯提高反演效果.

(2)圖12a為互相關和自相關能量平均化的反演結果，強干擾噪聲的出現說明反演已接近不穩定狀態，這說明使互相關能量和自相關能量占比相同的加權因子(公式(13))是不恰當的；圖5a和圖5b為自相關能量占優的反演結果，圖5b對應僅包含自相關能量，這是自相關能量占優的極限最大情況；圖5a為本文方法對應的加權因子，其以自相關值能量值的最大值為能量匹配準則.相比較圖12a的不穩定反演狀態，圖5a反演穩定且效果較好，而圖5b的質量卻低于圖5a，具體表現可以參考本文2.1節；這表明，自相關能量占比不僅會影響反演穩定性也會影響反演質量.

圖12 單層階梯模型應用 (a) 公式(13)和 (b) 公式(14)加權因子所得反演結果

綜上所述，在合理的互相關項和自相關項加權因子的作用下，自相關項的加入可以提升全波形反演效果，進一步地，加權因子的選擇會影響反演穩定性和反演質量.本節僅對加權因子的選擇與反演穩定性以及反演質量的關系進行了粗淺的初步探討，僅給出初步結論，明確自相關能量正作用反演效果的范圍和更加科學的加權因子選擇將是本論文下一步研究的重點.

4 結論

本文提出的基于垂直構造梯度優化的全波形反演方法，是對基于伴隨狀態法的全波形反演方法做出的垂直構造適應性改進.

(1)該方法通過構建模擬數據與觀測數據之間殘差逆傳波場的自相關項，實現了垂直構造位置處的梯度補償；同時，利用殘差逆傳波場能量最大值為準則使之與伴隨狀態法所得梯度結合，在顯著改善垂直構造反演效果的同時提升了水平構造的反演能力.

(2)多次正演和波場逆傳是全波形反演的主要計算量來源，對比基于伴隨狀態法的全波形反演方法流程，本文提出的方法不需要額外的波場正傳、反傳，僅需要額外計算殘差逆傳波場的自相關，僅增加十分有限的計算量，并不影響全波形反演計算效率.

(3)單層階梯模型、雙層階梯模型和局部Marmousi模型反演結果表明本文方法不僅適用于垂直構造，同樣適用于更具實際意義的傾斜構造，證明了本文方法的實用性和普適性.

(4)在后續研究中，可以從不同階數多次波信息能量級匹配與反演效果評價、含多次波信息目標函數建立、偏移距與反演效果評價、實際小偏移距數據測試以及梯度加權因子選擇等方面展開深入研究，以期進一步提升本方法反演效果.