波動(dòng)方程初至波多信息聯(lián)合反演方法

張建明， 董良國， 王建華

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

0 引言

在山地地震勘探中，建立一個(gè)相對精確的近地表速度模型是后續(xù)地震資料處理和解釋的關(guān)鍵步驟之一.包含有直達(dá)波、折射波和回轉(zhuǎn)波的地震初至波，最先被檢波器記錄到，并且通常能量較強(qiáng)，容易識別，因此，在構(gòu)建近地表速度模型時(shí)應(yīng)用最為廣泛.

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

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

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

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

1 方法原理

JTEW方法同時(shí)匹配模擬和觀測的地震初至波的走時(shí)、包絡(luò)和波形，通過權(quán)重因子調(diào)節(jié)三種信息的比例.定義JTEW的目標(biāo)函數(shù)為

(1)

其中，ω1和ω2分別是走時(shí)和包絡(luò)的權(quán)重因子，1-ω1-ω2為波形的權(quán)重因子，且0≤ω1≤1，0≤ω2≤1，JTEW的權(quán)重因子選擇是一個(gè)重要的課題，將在下一節(jié)中詳細(xì)討論.Jt、Je和Jw分別為走時(shí)、包絡(luò)和波形的目標(biāo)函數(shù)值，因?yàn)槿咄ǔ２皇且粋€(gè)數(shù)量級，所以我們對它們進(jìn)行了歸一化處理，其中Jt0、Je0和Jw0為歸一化因子，它們分別表示第一次迭代的走時(shí)、包絡(luò)和波形的目標(biāo)函數(shù)值.

在JTEW的目標(biāo)函數(shù)(1)中，走時(shí)、包絡(luò)和波形各自的目標(biāo)函數(shù)分別定義為

(2)

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

(3)

針對上述JTEW的目標(biāo)函數(shù)(1)，基于伴隨狀態(tài)法(Plessix, 2006)，通過正傳波場與伴隨波場的零延遲互相關(guān)就可以得到目標(biāo)函數(shù)(1)的梯度(具體推導(dǎo)過程見附錄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)分別為走時(shí)、包絡(luò)和波形的伴隨波場，與之對應(yīng)的伴隨震源分別為

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)

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

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

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

(9)

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

(10)

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

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

2 聯(lián)合反演權(quán)重因子的選擇策略

在JTEW方法中，權(quán)重因子的選擇對反演結(jié)果有很大的影響.下面通過對目標(biāo)函數(shù)性態(tài)分析，提出了一種合理有效的可變權(quán)重因子選擇策略.

為了檢驗(yàn)不同初至波信息的目標(biāo)函數(shù)的性態(tài)，本文設(shè)計(jì)了如圖1所示的模型和觀測方式，在地表的中間放炮，模型的底部接收，共501個(gè)檢波器，相鄰檢波器間距10 m.真實(shí)速度v=2500 m·s-1，速度擾動(dòng)范圍為1500～3500 m·s-1，計(jì)算出不同速度攝動(dòng)時(shí)的目標(biāo)函數(shù)如圖2所示.

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

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

圖2 目標(biāo)函數(shù)性態(tài)(a) 不同信息目標(biāo)函數(shù)的性態(tài); (b) 不同權(quán)重系數(shù)時(shí)的聯(lián)合反演目標(biāo)函數(shù).圖a中，綠線：波形反演的目標(biāo)函數(shù); 藍(lán)線：包絡(luò)反演的目標(biāo)函數(shù); 紅線：走時(shí)反演的目標(biāo)函數(shù); 黑線：聯(lián)合反演的目標(biāo)函數(shù)(w1=0.3，w2=0.4).圖b中，洋紅線：w1=0.7，w2=0.2; 藍(lán)線：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目標(biāo)函數(shù)可以兼具三者的優(yōu)點(diǎn)，在反演中可以通過調(diào)整不同信息的權(quán)重因子得到凸性不同、和對速度敏感性不同的目標(biāo)函數(shù).總體來說，走時(shí)權(quán)重從大變小，波形權(quán)重從小變大，聯(lián)合反演目標(biāo)函數(shù)凸性變差而對速度的敏感性變強(qiáng).根據(jù)這個(gè)特征，在反演初期階段設(shè)置較大的走時(shí)權(quán)重，主要更新背景速度模型；隨著迭代的進(jìn)行，逐漸減小走時(shí)權(quán)重而先增大包絡(luò)的權(quán)重，更新模型的中波數(shù)成分；最后減小走時(shí)和包絡(luò)的權(quán)重，增大波形權(quán)重更新模型的高波數(shù)成分，從而自然地實(shí)現(xiàn)分階段、多尺度反演，最終得到高分辨率的反演速度模型.這也為本文聯(lián)合反演權(quán)重因子的選擇策略的設(shè)計(jì)提供了依據(jù).

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

圖3 聯(lián)合反演權(quán)重因子系數(shù)Fig.3 Weighting factors in JTEW

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

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

表1 不同反演策略權(quán)重因子表Table 1 Weighting factors of different inversion strategies

3 理論模型試驗(yàn)

利用橫向速度周期變化模型和Foothill模型兩個(gè)數(shù)值實(shí)例，說明本文所提出的可變權(quán)重聯(lián)合反演方法相較于單信息反演、串聯(lián)反演以及固定權(quán)重聯(lián)合反演的優(yōu)勢.

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

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

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

圖4 橫向速度周期模型及不同方法反演結(jié)果(a) 真實(shí)速度模型; (b) 初至波波形反演結(jié)果; (c) 走時(shí)反演結(jié)果; (d) 走時(shí)+波形變權(quán)重聯(lián)合反演結(jié)果; (e) 走時(shí)+包絡(luò)+波形固定權(quán)重因子聯(lián)合反演結(jié)果; (f) 走時(shí)+包絡(luò)+波形變權(quán)重因子聯(lián)合反演結(jié)果.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反演結(jié)果(紅線)總是比波形反演(綠線)和固定權(quán)重聯(lián)合反演結(jié)果(藍(lán)線)更加接近真實(shí)模型(黑實(shí)線).這表明，JTEW方法在反演初期階段充分利用了走時(shí)信息，為后續(xù)反演階段提供了更加穩(wěn)定的初始模型，而在聯(lián)合反演的后期階段更加充分地利用了波形信息，有效提高了速度模型的分辨率.可見，與單信息反演和固定權(quán)重聯(lián)合反演相比，JTEW方法在反演精度方面具有明顯的優(yōu)勢.

圖5 不同深度處的反演速度模型的水平速度剖面(a) 150 m深度; (b) 200 m深度. 黑色實(shí)線和虛線分別表示真實(shí)模型和初始模型；藍(lán)色和綠色實(shí)線分別表示固定權(quán)重因子聯(lián)合反演和波形反演結(jié)果；紅線表示變權(quán)重因子聯(lián)合反演結(jié)果.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顯示了不同反演速度模型與真實(shí)模型第101炮的初至波波形擬合情況.可以看到，對應(yīng)于初始模型的理論波形(藍(lán)色)與“觀測”波形(紅色)匹配較差(圖6a)，在圖6a中的局部區(qū)域A和B，“觀測”波形和理論波形相位差異超過了1/2周期.對應(yīng)于初至波波形反演結(jié)果(圖4b)的初至波，小偏移距的初至波波形與“觀測”初至波波形擬合較好，但是大偏移距的初至波波形仍然擬合較差(圖6b).這是由于初始模型的估計(jì)偏離了真實(shí)模型，導(dǎo)致反演陷入了局部極值.

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

圖6 不同模型模擬的第101炮初至波形與觀測初至波形的擬合情況(a) 初始模型； (b) 波形反演模型； (c) 固定權(quán)重聯(lián)合反演模型; (d) 變權(quán)重聯(lián)合反演模型. 紅色為觀測初至波形；藍(lán)色為不同反演結(jié)果對應(yīng)的正演初至波形.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模型試驗(yàn)

復(fù)雜Foothill模型橫縱向網(wǎng)格數(shù)為1255×100，網(wǎng)格間距20 m.地面均勻激發(fā)126炮，每炮最多401個(gè)檢波器均勻分布在地表，最大偏移距為4 km.正演模擬采用基于縱向坐標(biāo)變換的起伏地表地震波模擬方法，其中第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)，進(jìn)行了四種初至波反演方法的對比試驗(yàn)，包括波動(dòng)方程初至波包絡(luò)反演、初至波波形反演、走時(shí)+包絡(luò)+波形串聯(lián)反演、以及走時(shí)+包絡(luò)+波形變權(quán)重因子聯(lián)合反演.

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

圖9 地表下不同深度處的反演速度模型的速度剖面(a) 地表下200 m； (b) 地表下400 m. 黑色實(shí)線和虛線分別表示真實(shí)模型和初始模型；藍(lán)色和綠色實(shí)線分別表示波形反演和包絡(luò)反演結(jié)果；紅線表示聯(lián)合反演結(jié)果.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.

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

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

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

圖11顯示了不同反演模型與真實(shí)模型之間初至波互相關(guān)提取的走時(shí)殘差.通過對比我們可以發(fā)現(xiàn)三個(gè)現(xiàn)象：

圖11 在不同反演結(jié)果上模擬的第66炮初至波和“觀測”初至波的互相關(guān)函數(shù)(a) 初始模型; (b) 走時(shí)反演; (c) 走時(shí)+包絡(luò)+波形串聯(lián)反演; (d) 走時(shí)+包絡(luò)+波形聯(lián)合反演.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發(fā)現(xiàn)，經(jīng)過初至波波動(dòng)方程走時(shí)反演，互相關(guān)函數(shù)的極大值更加靠近零值，說明經(jīng)過波動(dòng)方程走時(shí)反演后初至波走時(shí)信息已經(jīng)得到較好的匹配.

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

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

4 實(shí)際資料測試

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

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

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

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

圖13 速度模型(a) 初始模型; (b) 伴隨狀態(tài)法走時(shí)層析迭代65次; (c) 聯(lián)合反演迭代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)聯(lián)合反演模型的計(jì)算走時(shí)(藍(lán)點(diǎn))與拾取走時(shí)(紅線)Fig.14 The calculated traveltime (blue dot) and picked traveltime (red line) on (a) initial and (b) JTEW model.

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

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

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

5 結(jié)論

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

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

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

(3)模型復(fù)雜時(shí)，在串聯(lián)反演過程中，走時(shí)的目標(biāo)函數(shù)會(huì)發(fā)生漂移現(xiàn)象(本質(zhì)上是周波跳躍問題).本文提出的JTEW方法可以一定程度上緩解這個(gè)問題.

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

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

附錄A 聯(lián)合反演梯度推導(dǎo)

時(shí)間域聲波方程為

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

(A1)

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

(A2)

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

傳統(tǒng)全波形反演的目標(biāo)函數(shù)一般定義為模擬數(shù)據(jù)和觀測數(shù)據(jù)殘差的L2泛函：

(A3)

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

(A4)

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

在地震初至波包絡(luò)反演中，目標(biāo)函數(shù)定義為包絡(luò)對數(shù)殘差的L2泛函：

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

(A5)

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

fe(xr,t;xs)=

在波動(dòng)方程走時(shí)反演中，目標(biāo)函數(shù)定義為相對走時(shí)差的L2泛函(Luo and Schuster, 1991)：

(A7)

Δτ(xs,xr)為對應(yīng)于炮點(diǎn)位置xs和檢波點(diǎn)位置xr的觀測初至波和模擬初至波的走時(shí)殘差.該走時(shí)殘差是通過互相關(guān)函數(shù)的最大值自動(dòng)拾取的，即：

(A8)

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

(A9)

在JTEW方法中，目標(biāo)函數(shù)中同時(shí)匹配地震初至波的走時(shí)、包絡(luò)和波形信息.三種信息的伴隨源結(jié)合聯(lián)合目標(biāo)函數(shù)中的權(quán)重加在一起，構(gòu)成一個(gè)多信息聯(lián)合的復(fù)合伴隨源，通過反傳伴隨源得到的伴隨波場與正傳波場做零延遲互相關(guān)即為(4)式的聯(lián)合反演的梯度表達(dá)式.