波動方程初至波多信息聯合反演方法

張建明， 董良國， 王建華

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

0 引言

在山地地震勘探中，建立一個相對精確的近地表速度模型是后續地震資料處理和解釋的關鍵步驟之一.包含有直達波、折射波和回轉波的地震初至波，最先被檢波器記錄到，并且通常能量較強，容易識別，因此，在構建近地表速度模型時應用最為廣泛.

初至波中的不同信息(如走時、相位、振幅、包絡、波形等)對速度參數的敏感性和反演能力不同，在速度反演時對初始速度模型的依賴程度不同，反演的收斂性也不同.其中，地震波走時隨速度變化具有更好的線性特征，它是近地表速度建模中最常用的一種信息.初至波走時層析方法可以基于不同的正演引擎，如射線追蹤(Zhu et al., 1992)、Eikonal方程求解(Sei and Symes, 1994, 1995；Taillandier et al., 2009；Waheed et al., 2014; Waheed et al., 2016)，也可以基于波動方程進行初至波走時反演(WTI)(Luo and Schuster, 1991)等.

然而，由于地震波走時主要受到大尺度速度變化的影響，走時層析方法往往只能得到分辨率較低的宏觀速度模型，有時無法滿足實際地震勘探對地下介質精細刻畫的要求.為了充分利用初至波中的動力學信息，進一步提高近地表速度建模的精度，有學者提出了初至波波形層析成像(EWT)方法(Sheng et al., 2006)，這種方法是全波形反演(Lailly, 1883；Tarantola, 1984；Pratt et al., 1998)在近地表建模方面的一種應用.EWT可以利用包含在地震初至波中的全部信息，從理論上講，EWT比走時層析具有更高的反演精度，但由于地震波傳播的復雜性，使用波形殘差構造目標函數，必然導致反演的強烈非線性，在初始模型遠離真實模型、觀測和模擬的初至波走時差大于1/2周期時會產生周波跳躍而致使EWT反演陷入局部極值(Virieux and Operto, 2009)，子波未知以及噪聲干擾等一系列實際問題也嚴重阻礙了波形反演的實用化進程.因此，在現階段的實際速度建模中，波形反演還無法完全取代傳統的走時層析技術，而是應該根據勘探的不同階段，采取分階段、分尺度逐漸提高反演精度的策略(Bunks et al,. 1995; 董良國等, 2013, 2015).在此指導思想下，在反演的初期階段，應充分利用初至波的走時信息(Luo and Schuster, 1991； Zhang and Toks?z，1998)或相位信息(Choi and Alkhalifah, 2013)來更新近地表速度模型的長波長分量，得到宏觀的背景速度場；然后逐漸利用地震波的低中頻信息進一步提高反演精度，例如利用初至波包絡反演近地表速度的方法(EEI)(敖瑞德等, 2015)；最后使用波形信息，得到高精度的速度模型.

這種常規的走時、包絡和波形串聯反演策略雖然也可以得到較高精度的速度模型，但是，在后期利用初至波的包絡或波形的反演階段，由于缺少走時信息的約束，在初至波的包絡或波形匹配時會產生一定程度的周波跳躍(Liu and Zhang, 2017).因此，Liu和Zhang(2017)提出了一種近地表速度結構建模的走時、包絡和波形的聯合反演方法，在不同的反演階段，每一次迭代過程中走時、包絡和波形信息都同時被利用.然而，在Liu和Zhang(2017)所提出的方法中，走時匹配部分采用的是基于射線追蹤的層析方法；包絡和波形匹配部分基于波動方程正演模擬.這導致他們所提出的聯合反演方法有3個不足之處：(1)在每一次迭代中，射線追蹤和波動方程正演需要分別計算；(2)走時匹配基于射線理論，使得反演方法的精度受到射線理論的高頻假設制約，致使他們所提出的聯合反演方法難以適用于復雜速度變化的情況；(3)走時匹配部分依賴于精確的走時拾取，在數據質量不佳的情況下，拾取走時必然耗費巨大的人力和時間.

為此，本文提出了一種統一基于波動方程正演引擎的初至波走時、包絡和波形聯合反演方法(Joint traveltime, envelope and waveform inversion, JTEW)來建立近地表速度模型.該方法統一使用波動方程正演引擎計算正傳波場和伴隨波場，正傳波場可以同時應用于走時、包絡和波形的匹配，伴隨波場應用于梯度的計算.使用初至波形的互相關可以自動提取模擬和觀測初至波的時差，無需再進行額外的射線追蹤.通過權重因子可以分階段自動調節不同信息的比例，自然地實現分階段、多尺度反演，為初至波波形反演的實用化提供了一種有效途徑.

1 方法原理

JTEW方法同時匹配模擬和觀測的地震初至波的走時、包絡和波形，通過權重因子調節三種信息的比例.定義JTEW的目標函數為

(1)

其中，ω1和ω2分別是走時和包絡的權重因子，1-ω1-ω2為波形的權重因子，且0≤ω1≤1，0≤ω2≤1，JTEW的權重因子選擇是一個重要的課題，將在下一節中詳細討論.Jt、Je和Jw分別為走時、包絡和波形的目標函數值，因為三者通常不是一個數量級，所以我們對它們進行了歸一化處理，其中Jt0、Je0和Jw0為歸一化因子，它們分別表示第一次迭代的走時、包絡和波形的目標函數值.

在JTEW的目標函數(1)中，走時、包絡和波形各自的目標函數分別定義為

(2)

其中，Δτ(xs,xr)為對應于炮點xs和檢波點xr的觀測和模擬初至波的走時差；E和E0分別表示觀測與模擬的初至波形的包絡, ln表示取對數；Δu(xr,t;xs)=u(xr,t;xs)-u0(xr,t;xs)為觀測與模擬初至波的波形殘差.其中理論模擬波場u(x,t;xs)滿足下列點震源子波為f(t)的波動方程：

(3)

針對上述JTEW的目標函數(1)，基于伴隨狀態法(Plessix, 2006)，通過正傳波場與伴隨波場的零延遲互相關就可以得到目標函數(1)的梯度(具體推導過程見附錄A):

+u′e(x,t;xs)+u′w(x,t;xs)]}dt，

(4)

這里，c(x)為空間位置x處的速度.u′t(x,t;xs)、u′e(x,t;xs)和u′w(x,t;xs)分別為走時、包絡和波形的伴隨波場，與之對應的伴隨震源分別為

f′t(xr,t;xs)=

(5)

(6)

f′w(xr,t;xs)=

(7)

令

u′all(x,t;xs)=u′t(x,t;xs)+u′e(x,t;xs)

+u′w(x,t;xs),

(8)

為總的伴隨波場，它是由三種信息聯合的復合伴隨震源f′(xr,t;xs)以波動方程正演引擎反向傳播得到，其中

f′(xr,t;xs)=f′t(xr,t;xs)+f′e(xr,t;xs)

+f′w(xr,t;xs).

(9)

這樣，在每次迭代中只需要使用波動方程(3)對上述(9)式的復合伴隨震源f′(xr,t;xs)進行一次反向傳播，將得到的伴隨波場u′all(x,t;xs)與震源端的正傳波場u(x,t;xs)進行零延遲互相關(公式4)即可得到JTEW目標函數的梯度.然后通過線性搜索算法可以很方便地獲得一個優化的步長αn，這樣就可以利用下面(10)式的局部優化算法進行速度迭代，得到最終的速度反演結果.

(10)

另外，由于(5)、(6)、(7)式的計算量與波場模擬相比幾乎可以忽略不計，因此，本文提出的JTEW方法與單獨采用波形反演的計算量基本相當.

從以上方法可以看出，本文的JTEW方法彌補了Liu和Zhang(2017)聯合反演方法的不足，主要優勢體現在以下三個方面：(1)在每一次迭代過程中，統一使用波動方程正演引擎計算正傳波場和伴隨波場進行梯度計算，無需再進行額外的射線追蹤正演計算；(2)走時、包絡和波形匹配都是基于波動方程正演模擬結果，提高了方法在復雜介質情況下的適用性；(3)通過初至波形互相關自動提取模擬和觀測初至波的時差，無需人工拾取走時，省時省力.

2 聯合反演權重因子的選擇策略

在JTEW方法中，權重因子的選擇對反演結果有很大的影響.下面通過對目標函數性態分析，提出了一種合理有效的可變權重因子選擇策略.

為了檢驗不同初至波信息的目標函數的性態，本文設計了如圖1所示的模型和觀測方式，在地表的中間放炮，模型的底部接收，共501個檢波器，相鄰檢波器間距10 m.真實速度v=2500 m·s-1，速度擾動范圍為1500～3500 m·s-1，計算出不同速度攝動時的目標函數如圖2所示.

圖1 模型和觀測方式示意圖Fig.1 Schematic diagram of the model and survey geometry

圖2a表明，波形目標函數(綠線)凸性最差，說明波形反演是一個強烈的非線性問題，只有初始模型足夠接近真實模型時，目標函數才能正確收斂到全局極值.包絡目標函數(藍線)凸性較好，但是當初始模型距離真實模型較遠時，特別是當速度偏低時，目標函數對速度變化不敏感，速度很難向正確的方向更新.走時目標函數(紅線)具有很好的凸性，但是同時它對速度變化的敏感性最低，因此，走時反演通常只能得到低分辨率的速度模型.

圖2 目標函數性態(a) 不同信息目標函數的性態; (b) 不同權重系數時的聯合反演目標函數.圖a中，綠線：波形反演的目標函數; 藍線：包絡反演的目標函數; 紅線：走時反演的目標函數; 黑線：聯合反演的目標函數(w1=0.3，w2=0.4).圖b中，洋紅線：w1=0.7，w2=0.2; 藍線：w1=0.3，w2=0.5; 紅線：w1=0.1，w 2=0.6; 綠線：w1=0.1，w 2=0.3; 黑線：w1=0.05，w 2=0.1.Fig.2 Objective function behavior(a) The behavior of different information objective functions; (b) The JTEW objective functions with different weighting factors. Green line: EWT objective function; blue line: EEI objective function; red line: WTI objective function; black line: JTEW objective function, where w1=0.3, w2=0.4. Magenta line: w1=0.7, w 2=0.2; Blue line: w1=0.3, w 2=0.5; Red line: w1=0.1,w 2=0.6; Green line:w1=0.1, w 2=0.3; Black line: w1=0.05, w 2=0.1.

圖2b表明,初至波JTEW目標函數可以兼具三者的優點，在反演中可以通過調整不同信息的權重因子得到凸性不同、和對速度敏感性不同的目標函數.總體來說，走時權重從大變小，波形權重從小變大，聯合反演目標函數凸性變差而對速度的敏感性變強.根據這個特征，在反演初期階段設置較大的走時權重，主要更新背景速度模型；隨著迭代的進行，逐漸減小走時權重而先增大包絡的權重，更新模型的中波數成分；最后減小走時和包絡的權重，增大波形權重更新模型的高波數成分，從而自然地實現分階段、多尺度反演，最終得到高分辨率的反演速度模型.這也為本文聯合反演權重因子的選擇策略的設計提供了依據.

依據上述分析，本文設計了一種全自動的可變權重因子的聯合反演策略(圖3)，總共15個反演階段，三個固定權重因子對應一個反演階段.在第一個反演階段，走時、包絡和波形的權重分別為1、0、0，當該階段的目標函數不再下降時，自動跳到第二個反演階段，新階段的權重因子按圖3自動調節，以此類推，直到最后一個反演階段結束.整個反演過程中走時的權重因子從1指數衰減到0，波形權重因子與走時權重因子的變化是反對稱的，即從0指數增大到1，其余部分是包絡的權重因子.

圖3 聯合反演權重因子系數Fig.3 Weighting factors in JTEW

這樣設計權重因子，充分考慮到了走時與速度之間的弱非線性關系，使得JTEW在反演的初期階段能夠獲得一個可靠的宏觀速度模型.在后續的反演階段，走時與速度之間的這種弱非線性不會被丟棄，而是作為包絡和波形的約束加載在目標函數中，并且逐漸減小其比重而加大波形信息的權重，以便在聯合反演的后期階段，波形信息被充分利用，從而自動實現分階段、多尺度反演，最終得到一個高精度的近地表速度模型.

在本文的JTEW方法中，也可以根據實際情況靈活地設計不同的權重策略而將反演分成多個階段，表1展示的是幾種不同反演策略的權重設置及其實際意義.可以看出，無論是單信息反演、固定權重聯合反演或者是常規串聯反演，均是可變權重聯合反演的一個特例.下文中的JTEW若沒有特別說明，均是指可變權重的聯合反演.聯合反演最佳權重策略該如何設計(如設置多少個反演階段，走時、包絡和波形的權重因子是何種變化形式，第一個和最后一個階段的權重值分別設置多大，等等)，還有待進一步優化.

表1 不同反演策略權重因子表Table 1 Weighting factors of different inversion strategies

3 理論模型試驗

利用橫向速度周期變化模型和Foothill模型兩個數值實例，說明本文所提出的可變權重聯合反演方法相較于單信息反演、串聯反演以及固定權重聯合反演的優勢.

(1)橫向速度周期變化模型試驗

模型的橫縱向網格數為501×51，水平和垂直網格間距均為10 m.在試驗過程中，采用有限差分法正演模擬地震波，震源為主頻15 Hz的雷克子波，101炮均勻分布于地表，炮間距50 m，每炮501個檢波器在地表均勻分布，相鄰檢波器間距10 m.真實模型如圖4a所示，初始模型是速度值為2500 m·s-1的常速初始模型.使用相同的常速初始模型，進行了五種不同反演方法的數值試驗，包括波動方程初至波走時反演(Luo and Schuster, 1991)、初至波波形反演(Sheng et al., 2006)、固定權重因子(ω1=0.7，ω2=0.2)聯合反演、走時+波形變權重因子聯合反演、以及JTEW反演.

由于常速初始模型與真實模型相差較大，初至波波形反演結果明顯地陷入了局部極值(圖4b)，反演結果完全失真.初至波走時反演只能得到分辨率較低的背景速度模型(圖4c).與固定權重因子聯合反演結果(圖4e)相比，走時+波形變權重因子聯合反演結果(圖4d)和JTEW反演結果(圖4f)，不僅有效更新深度更大，而且在近地表分辨率更高.這是因為在固定權重因子聯合反演的整個迭代過程中，各部分的權重保持不變(ω1=0.7，ω2=0.2)，導致走時、包絡和波形三部分信息只有權重最大的走時信息被充分利用，波形信息利用不足，若減小走時的權重，又會造成反演的不穩定.而JTEW反演方法則可以通過改變不同信息的權重，在反演的不同階段強調利用不同的信息，自然地實現多尺度反演，最終得到分辨率更高、有效更新深度更大的反演速度模型.值得說明的是，盡管只利用走時和波形兩種信息進行變權重因子聯合反演，也得到了較好的反演效果(圖4d)，但包絡信息的加入，使得JTEW方法得到了精度更高的反演結果(圖4f).

圖4 橫向速度周期模型及不同方法反演結果(a) 真實速度模型; (b) 初至波波形反演結果; (c) 走時反演結果; (d) 走時+波形變權重聯合反演結果; (e) 走時+包絡+波形固定權重因子聯合反演結果; (f) 走時+包絡+波形變權重因子聯合反演結果.Fig.4 Lateral periodic velocity model and inversion results with different inversion methods(a) True velocity model; (b) EWT inversion result; (c) WTI inversion result; (d) WTI+EWT Joint inversion result with variable weighting factors; (e) JTEW inversion result with fixed weighting factors; (f) JTEW with variable weighting factors.

圖5顯示的是150 m和200 m深度處的水平速度曲線，可以看出，JTEW反演結果(紅線)總是比波形反演(綠線)和固定權重聯合反演結果(藍線)更加接近真實模型(黑實線).這表明，JTEW方法在反演初期階段充分利用了走時信息，為后續反演階段提供了更加穩定的初始模型，而在聯合反演的后期階段更加充分地利用了波形信息，有效提高了速度模型的分辨率.可見，與單信息反演和固定權重聯合反演相比，JTEW方法在反演精度方面具有明顯的優勢.

圖5 不同深度處的反演速度模型的水平速度剖面(a) 150 m深度; (b) 200 m深度. 黑色實線和虛線分別表示真實模型和初始模型；藍色和綠色實線分別表示固定權重因子聯合反演和波形反演結果；紅線表示變權重因子聯合反演結果.Fig.5 Horizontal velocity profiles of the inversion velocity model at (a) y=150 m and (b) y=200 mThe black solid line and the dashed line represent the true model and the initial model, respectively; The blue and the green solid line represent the JTEW result with fixed weighting factors and EWT result, respectively; The red line represents the JTEW results with variable weighting factors.

圖6顯示了不同反演速度模型與真實模型第101炮的初至波波形擬合情況.可以看到，對應于初始模型的理論波形(藍色)與“觀測”波形(紅色)匹配較差(圖6a)，在圖6a中的局部區域A和B，“觀測”波形和理論波形相位差異超過了1/2周期.對應于初至波波形反演結果(圖4b)的初至波，小偏移距的初至波波形與“觀測”初至波波形擬合較好，但是大偏移距的初至波波形仍然擬合較差(圖6b).這是由于初始模型的估計偏離了真實模型，導致反演陷入了局部極值.

固定權重的聯合反演結果較初至波波形反演結果的波形擬合有了大幅提升，但仍然存在一定殘差(圖6c)，而JTEW得到了比較完美的波形匹配效果(圖6d).這是因為固定權重的聯合反演考慮了走時的信息，非線性程度降低.然而，固定權重的聯合反演只強調利用三種信息中權重因子比較大的走時信息，導致權重值小的波形信息不能被充分利用.JTEW則能夠在反演初期階段充分利用走時信息反演得到較光滑的宏觀速度模型，在反演后期期階段充分利用波形信息提高模型的分辨率.從上面的波形匹配結果也可以看出，JTEW方法可以最大程度擬合“觀測”數據，最終提高了速度反演的精度.

圖6 不同模型模擬的第101炮初至波形與觀測初至波形的擬合情況(a) 初始模型； (b) 波形反演模型； (c) 固定權重聯合反演模型; (d) 變權重聯合反演模型. 紅色為觀測初至波形；藍色為不同反演結果對應的正演初至波形.Fig.6 The 101th shot waveform of (a) initial model; (b) EWT inversion model; (c) JTEW inversion model with fixed weighting factors and (d) JTEW inversion model with variable weighting factorsThe red line is the observed waveform; The blue line is the forward waveform corresponding to different inversion results.

(2)Foothill模型試驗

復雜Foothill模型橫縱向網格數為1255×100，網格間距20 m.地面均勻激發126炮，每炮最多401個檢波器均勻分布在地表，最大偏移距為4 km.正演模擬采用基于縱向坐標變換的起伏地表地震波模擬方法，其中第65炮的完整波形及反演所用的初至波波形如圖7所示.

圖7 Foothill模型第65炮“觀測”記錄以及反演所用的初至波(a) 完整波形; (b) 反演使用的初至波波形.Fig.7 Waveform of the 65th shot for Foothill model(a) Full waveform; (b) First arrival waveform used in inversion.

采用相同的常梯度初始速度模型(圖8b)，進行了四種初至波反演方法的對比試驗，包括波動方程初至波包絡反演、初至波波形反演、走時+包絡+波形串聯反演、以及走時+包絡+波形變權重因子聯合反演.

初至波波形反演結果(圖8c)和初至波包絡反演結果(圖8d)都比較差，在圖9的地下不同深度的速度曲線上，初至波波形反演結果(藍線)和初至波包絡反演結果(綠線)都是在初始模型(水平黑虛線)附近抖動，這是由于初始模型與真實模型相差較大，致使波形和包絡反演均陷入了局部極值，模型沒有得到有效更新.采用相同的常梯度初始模型(圖8b)，盡管走時+包絡+波形串聯反演也得到了較理想的反演結果(圖8e)，但本文提出的JTEW方法得到的反演結果精度更高(圖8f和圖9中的紅色曲線)更接近于實際模型，說明了多信息、分階段JTEW反演相對于單信息反演具有明顯的優越性.

圖9 地表下不同深度處的反演速度模型的速度剖面(a) 地表下200 m； (b) 地表下400 m. 黑色實線和虛線分別表示真實模型和初始模型；藍色和綠色實線分別表示波形反演和包絡反演結果；紅線表示聯合反演結果.Fig.9 The velocity profiles of the inversion velocity model at (a) 200 m and (b) 400 m below the surfaceThe black solid line and the dashed line represent the true model and the initial model, respectively；The blue solid line and the green solid line represent the single EWT and EEI results, respectively；The red line represents the JTEW result.

另外，目標函數下降曲線(圖10a)顯示，在走時+包絡+波形串聯反演的前120次迭代過程中，在走時目標函數值(紅線)下降的同時，包絡(黑線)和波形(藍線)的目標函數值并不是持續下降，而是在一些局部出現了增大的現象(第30次迭代附近).第120次走時反演迭代結束后，包絡和波形目標函數值下降的同時，走時目標函數值卻又出現了升高的現象(圖10a中第120～320代)，我們稱之為走時目標函數的漂移現象.這說明在走時+包絡+波形串聯反演過程中，走時反演階段結束后，走時信息本來得到了較好的匹配，然而在后期的包絡和波形反演過程中則放棄了走時信息，致使包絡和波形反演過程中缺少了走時信息的約束，波形匹配變好的同時走時匹配變差.圖10a中的這些現象，均說明多信息串聯反演無法同時較好地匹配走時、包絡和波形信息.

圖10 目標函數下降曲線(a) 串聯反演； (b) 聯合反演.Fig.10 Decrease of objective function with iterations for (a) sequential inversion and (b) JTEW

而在聯合反演過程中，走時、包絡和波形目標函數值均持續下降(圖10b).這是由于在聯合反演后期的迭代過程中都使用走時信息作為包絡和波形的約束，使得波形匹配變好的同時走時匹配也變好.因此，變權重因子聯合反演可以有效緩解走時目標函數漂移問題，這也是聯合反演較串聯反演的優勢所在.

圖11顯示了不同反演模型與真實模型之間初至波互相關提取的走時殘差.通過對比我們可以發現三個現象：

圖11 在不同反演結果上模擬的第66炮初至波和“觀測”初至波的互相關函數(a) 初始模型; (b) 走時反演; (c) 走時+包絡+波形串聯反演; (d) 走時+包絡+波形聯合反演.Fig.11 The cross-correlation function of the first arrival wave and the observed first arrival wave simulated by the different inversion results of the 66th shot(a) Initial model; (b) Travel time inversion; (c) Travel time, envelope, and waveform sequential inversion; (d) JTEW inversion.

(1)對比圖11a和圖11b發現，經過初至波波動方程走時反演，互相關函數的極大值更加靠近零值，說明經過波動方程走時反演后初至波走時信息已經得到較好的匹配.

(2)對比圖11b和圖11c發現，基于走時反演模型，再進行包絡和波形反演，互相關函數極大值在某些位置(特別是遠偏移距)反而偏離零值，說明從走時反演模型出發，串聯反演中的包絡和波形反演速度模型的某些位置(例如圖8e紅色箭頭所指位置)錯誤更新，造成速度反演結果對走時的預測變差.

(3)對比圖11c和圖11d發現，基于JTEW反演結果模擬的初至波和“觀測”初至波的互相關函數的極值基本在零值附近，說明在聯合反演過程中，由于走時信息的約束，波形匹配變好的同時走時匹配也變得更好.走時殘差對比說明了聯合反演可以同時使走時、包絡和波形匹配變好，充分體現了聯合反演較串聯反演的優勢.

4 實際資料測試

選擇某工區的一條二維測線試驗聯合反演的實際效果.該測線共有712炮，相鄰炮間距約50 m.每炮最多384道，道間距20～30 m，最大偏移距為4800 m.其中，第480炮原始地震數據和拾取的初至走時(藍線)如圖12所示.

圖12 實際地震資料第480炮原始記錄Fig.12 The original field data record of the 480th shot

首先對數據進行了必要的預處理，包括異常道的剔除、地震數據三維至二維的轉換、4～30 Hz的帶通濾波等.根據初至斜率，初步估算一個平均的地表速度值，并據此設計了如圖13a所示的常梯度初始速度模型，其中地表速度為1800 m·s-1，模型的最底部速度為2280 m·s-1.然后，將不依賴子波的JTEW方法應用于該測線，并將近地表反演結果與伴隨狀態法初至波走時層析(AST)反演結果進行對比.

在JTEW反演過程中，我們使用包絡互相關提取觀測初至波與模擬初至波的時差，以便盡可能地克服噪聲對初至波時差的影響.整個模型剖分為1574×101個網格，網格間距25 m×10 m.AST與JTEW反演分別迭代了65次和44次，反演結果分別見如圖13b和圖13c.相比于AST反演得到的宏觀模型(圖13b),在靠近地表以及某些的深部JTEW反演結果的分辨率更高一些.由于初始速度模型不夠準確，根據初始模型計算的初至走時(圖14a藍線)在遠偏移距處與拾取的走時差別比較大(圖14a紅線)，而根據JTEW最終反演的速度模型計算的理論走時與拾取走時的誤差明顯減小(圖14b).同時，利用JTEW最終反演的速度模型模擬的初至波與實際觀測的初至波擬合良好(圖15).

圖13 速度模型(a) 初始模型; (b) 伴隨狀態法走時層析迭代65次; (c) 聯合反演迭代44次的反演速度模型.Fig.13 Velocity model(a) Initial model; (b) The 65th Invertion velocity model of AST; (c) The 44th Invertion velocity model of JTEW.

圖14 (a)初始模型和(b)聯合反演模型的計算走時(藍點)與拾取走時(紅線)Fig.14 The calculated traveltime (blue dot) and picked traveltime (red line) on (a) initial and (b) JTEW model.

圖15 第700炮實際初至波與聯合反演反演模型模擬的初至波擬合情況Fig.15 Comparison between the observed and the calculated early waveform on JTEW model of the 700st shot

盡管JTEW結果的疊加剖面(圖16b)與走時反演結果的疊加剖面(圖16a)相比，剖面質量并沒有本質的改善(深度域剖面存在明顯的不同，構造形態也會有差別，因為兩種方法反演的速度模型圖13b和圖13c存在較大的差異)，但同相軸連續性還是得到了一定程度的提高(見黑色框部分)，JTEW反演結果的分辨率也更高一些，這都說明了JTEW方法提高了該測線近地表速度建模的精度.

圖16 (a)走時層析與 (b)聯合反演結果疊加剖面對比Fig.16 Comparison of stacked profiles between (a) traveltime inversion and (b) JTEW

5 結論

本文提出的基于波動方程的初至波多信息聯合反演方法(JTEW)在統一波動理論的框架下，分階段、分尺度地考慮利用地震初至波的不同信息，有效地提高了近地表速度建模的精度.綜合理論分析和兩個數值算例以及一條二維測線實際地震資料處理結果，證明了提出的JETW方法的有效性.

(1)基于統一的波動方程正演引擎，一次正演得到的波場可同時用于走時、包絡、波形反演梯度的計算，無需額外的射線追蹤，提高了反演效率.

(2)走時、包絡和波形匹配部分都是基于波動方程正演引擎，克服了射線理論中高頻假設的限制，提高了方法在復雜地質情況下的適用性.

(3)模型復雜時，在串聯反演過程中，走時的目標函數會發生漂移現象(本質上是周波跳躍問題).本文提出的JTEW方法可以一定程度上緩解這個問題.

(4)在初始速度模型較差的情況下，單獨使用包絡或波形反演容易陷入局部極值.本文提出的初至波JTEW方法，可分階段綜合利用初至波的走時、包絡和波形信息，有效地降低反演對初始模型的依賴，提高近地表速度模型的反演精度.

(5)初至波波動方程單信息走時反演、包絡反演、波形反演、以及固定權重聯合反演和串聯反演方法，均是本文所提出的JTEW方法框架下的一種特殊情況，通過選擇不同的權重策略，即可靈活地實現不同的反演方法.

附錄A 聯合反演梯度推導

時間域聲波方程為

L[u(x,t;xs)]=δ(x-xs)f(t)，

(A1)

其中，xs是震源位置，u(x,t;xs)空間位置x處的波場，L為聲波算子

(A2)

其中，c(x)為地下空間位置x處的速度.

傳統全波形反演的目標函數一般定義為模擬數據和觀測數據殘差的L2泛函：

(A3)

利用伴隨狀態法(Tarantola，1984；Plessix，2006)可以求出梯度表達式.

(A4)

其中u′(x,t;xs)為伴隨源fw(xr,t;xs)=u(xr,t;xs)-u0(xr,t;xs)的反傳波場.

在地震初至波包絡反演中，目標函數定義為包絡對數殘差的L2泛函：

-ln(E0(xr,t;xs))]2dt，

(A5)

基于伴隨狀態法可以推導出包絡反演的伴隨源為(Chi et al., 2013，2014；Wu et al.,2014;敖瑞德等，2015)：

fe(xr,t;xs)=

在波動方程走時反演中，目標函數定義為相對走時差的L2泛函(Luo and Schuster, 1991)：

(A7)

Δτ(xs,xr)為對應于炮點位置xs和檢波點位置xr的觀測初至波和模擬初至波的走時殘差.該走時殘差是通過互相關函數的最大值自動拾取的，即：

(A8)

其中，T為觀測初至波和模擬初至波的最大時差.借助連接函數的概念(Luo and Schuster,1991)，可以推導出波動方程走時反演的伴隨源：

(A9)

在JTEW方法中，目標函數中同時匹配地震初至波的走時、包絡和波形信息.三種信息的伴隨源結合聯合目標函數中的權重加在一起，構成一個多信息聯合的復合伴隨源，通過反傳伴隨源得到的伴隨波場與正傳波場做零延遲互相關即為(4)式的聯合反演的梯度表達式.