王義， 董良國

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

?

基于截斷牛頓法的VTI介質(zhì)聲波多參數(shù)全波形反演

王義， 董良國*

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

不同類別參數(shù)間的相互耦合使多參數(shù)地震全波形反演的非線性程度顯著增加，地震波速度與各向異性參數(shù)取值數(shù)量級的巨大差異也會使反演問題的性態(tài)變差.合理使用Hessian逆算子可以減弱這兩類問題對反演的影響，提高多參數(shù)反演的精度，而截斷牛頓法是一種可以比較準(zhǔn)確地估計Hessian逆算子的優(yōu)化方法.本文采用截斷牛頓法在時間域進(jìn)行了VTI介質(zhì)的聲波雙參數(shù)同時反演的研究.不同模型的反演試驗(yàn)表明，在VTI介質(zhì)聲波雙參數(shù)同時反演中，截斷牛頓法比有限內(nèi)存BFGS(Limited-memory Broyden-Fletcher-Goldfarb-Shanno, L-BFGS)法能更準(zhǔn)確地估計Hessian逆算子，進(jìn)而較好地平衡兩類不同參數(shù)的同時更新，得到了比較精確的反演結(jié)果.

多參數(shù)； 全波形反演； 各向異性； 耦合； Hessian逆算子； 截斷牛頓法

1 引言

地震全波形反演(Full waveform inversion，F(xiàn)WI)綜合利用觀測地震記錄的波形信息來估計地下介質(zhì)的物性參數(shù)(如密度、速度和各向異性參數(shù)等)，在油氣勘探、開發(fā)和地球動力學(xué)等領(lǐng)域都有廣闊的應(yīng)用前景(劉國峰等，2012；石玉梅等，2014；Fichtner et al.，2009).隨著對復(fù)雜地下介質(zhì)研究需求的增長和大規(guī)?？茖W(xué)計算能力的提升，F(xiàn)WI逐漸由單一縱波速度的反演向多個參數(shù)同時反演的方向發(fā)展(Virieux and Operto，2009).

地下介質(zhì)的各向異性使地震波沿不同方位角以不同的速度傳播，在地下存在各向異性情況時，使用各向同性FWI會得到錯誤的速度結(jié)構(gòu)(Plessix and Cao, 2011).各向異性參數(shù)在巖性識別和裂縫特征描述等方面均有重要作用.利用各向異性FWI能夠更為精確地反演地下介質(zhì)的速度和各向異性參數(shù)(張美根，2003；Gholami，2012).但是，不同類別參數(shù)間存在的相互耦合效應(yīng)增加了多參數(shù)FWI的非線性程度.同時，速度與各向異性參數(shù)取值數(shù)量級的巨大差異也會使反演問題的性態(tài)變差(Operto et al., 2013).

這兩個問題都可以通過合理使用Hessian逆算子來部分解決.Hessian逆算子在FWI中起著重要作用.在單參數(shù)反演中，該算子可以對照明不足位置處的參數(shù)重新聚焦(Pratt et al.,1998).在多參數(shù)反演中，Hessian算子的非對角塊矩陣中的元素描述了不同類別參數(shù)間的耦合效應(yīng)(Operto et al., 2013)，使用Hessian逆算子可以減弱這種耦合效應(yīng)對FWI的影響(Metivier et al., 2014a).但由于巨大的計算和存儲消耗，顯式計算和存儲Hessian及其逆算子是不現(xiàn)實(shí)的，這使得牛頓法和高斯牛頓法不能得到廣泛應(yīng)用(Virieux and Operto，2009).常規(guī)的梯度類方法如梯度法和非線性共軛梯度(Conjugate gradient, CG)法通過使用對角近似的Hessian逆算子對梯度進(jìn)行預(yù)條件來加速反演的收斂過程(Shin et al., 2001；Lee et al.，2010).L-BFGS法則更進(jìn)一步，利用一定數(shù)目的模型和梯度信息構(gòu)造正定矩陣來近似Hessian逆算子(Nocedal and Wright, 2006).這種近似已不再是簡單的對角近似，使得L-BFGS法的穩(wěn)健性和收斂速度都優(yōu)于常規(guī)的梯度類方法，在FWI中得到了廣泛應(yīng)用(Guitton and Symes, 2003； Brossier et al., 2009).但是在很多病態(tài)反演問題中，Hessian矩陣的非正定性使得L-BFGS法的性能變差，截斷牛頓(Truncated Newton, TN)法則可能更加有效(Santosa and Symes, 1988；Nash and Nocedal, 1991；Metivier et al.,2013).它在每步迭代中通過線性CG法求解牛頓方程來獲取模型增量, 以此來更加準(zhǔn)確地利用Hessian逆算子的信息.

TN法在地震FWI領(lǐng)域逐漸得到研究和應(yīng)用(Santosa and Symes, 1988； Akcelik, 2002；Epanomeritakis et al.,2008).最近兩年，Metivier等(2014b) 在頻率域進(jìn)行了縱波速度的TN法全波形反演的研究，并用試驗(yàn)說明了存在強(qiáng)多次反射波情況下，TN法比截斷高斯牛頓法、L-BFGS法以及非線性CG法能夠更加準(zhǔn)確地估計Hessian逆算子.在多參數(shù)方面，Bae等(2012)采用截斷高斯牛頓法研究了各向同性聲彈耦合介質(zhì)的多參數(shù)反演，得到了縱橫波速度和密度，但并未進(jìn)一步研究TN法的反演性能.

在VTI介質(zhì)FWI方面，Lee 等(2010)用預(yù)條件梯度法反演了4個彈性系數(shù)，他們更注重采用分步反演的策略來克服多參數(shù)反演中的病態(tài)性.Plessix和Cao(2011)在聲波VTI介質(zhì)中分別對NMO速度和水平速度使用L-BFGS法進(jìn)行單參數(shù)反演，前者的反演結(jié)果偏高，而用后者的結(jié)果轉(zhuǎn)換出的等價各向異性參數(shù)的結(jié)果與真實(shí)模型相差較多.Gholami(2012)運(yùn)用L-BFGS法進(jìn)行了聲波和彈性波的VTI介質(zhì)多參數(shù)同時反演.他使用參數(shù)標(biāo)準(zhǔn)化(Normalization)將不同類別的參數(shù)調(diào)整到相同的數(shù)量級上，使反演過程更加穩(wěn)健，但得到的往往是過度估計的速度參數(shù)和低估的各向異性參數(shù).Metivier等(2014b) 運(yùn)用TN法在頻率域進(jìn)行了聲波VTI介質(zhì)的縱波速度的單參數(shù)全波形反演.

本文將截斷牛頓法應(yīng)用到VTI介質(zhì)聲波的雙參數(shù)同時反演中，并分析了它在反演中的性能.在介紹了TN法的基本原理以及分析了Hessian矩陣的性態(tài)之后，通過時間域不同模型的合成數(shù)據(jù)數(shù)值反演試驗(yàn)，對比研究了TN法與L-BFGS法的性能.

2 截斷牛頓法原理

設(shè)最小二乘目標(biāo)函數(shù)為

(1)

其中，m為模型參數(shù)，dobs和f(m)分別為觀測和模擬記錄，f是聯(lián)系模型參數(shù)和模擬記錄的算子，在時間域*代表轉(zhuǎn)置.求解(1)式常采用局部迭代優(yōu)化方法，其模型更新公式為

mk+1=mk+αkpk，

(2)

其中，pk和αk分別為第k次迭代的更新方向和步長.對目標(biāo)函數(shù)的二階Taylor展開進(jìn)行最小化，可以得到牛頓法的更新方向所滿足的方程：

(3)

由于顯式Hessian及其逆矩陣的計算和存儲消耗巨大，在大規(guī)模FWI中求取精確的牛頓更新方向是不現(xiàn)實(shí)的，使得具有局部快速收斂優(yōu)勢的牛頓法無法獲得廣泛應(yīng)用.TN法則可以克服這個困難，它在每步迭代中使用線性CG法求解(3)式來獲取比較準(zhǔn)確的牛頓更新方向，而線性CG法只需整體使用Hessian矩陣和已知向量的乘積.每個Hessian向量乘積都可以利用有限差分法(NocedalandWright, 2006)或者二階伴隨狀態(tài)法(FichtnerandTrampert,2011；Metivieretal.,2013)通過四次正演計算得到.本文使用有限差分法計算Hessian向量乘積，計算公式為

(4)

其中，v為已知向量，h為一個很小的正數(shù).相比于牛頓法，TN法的計算量已經(jīng)可以承受，在大規(guī)模反演問題中獲得了應(yīng)用并得到了較好的反演結(jié)果(NashandNocedal, 1991；Metivieretal.,2014b).

下面的算法1給出了TN法的流程，內(nèi)層while循環(huán)代表了線性CG迭代過程.為避免計算量過大，CG迭代需要設(shè)定合理的停止條件(Wang等, 2013).當(dāng)CG迭代中遇到負(fù)曲率方向(即算法1中β1<0)，表示當(dāng)前的Hessian矩陣為非正定矩陣，繼續(xù)迭代會得到不下降的更新方向，因而將CG迭代停止.通過這種CG過程，TN法主要利用了Hessian逆算子

的大特征值的作用，因而對模型更新的計算有內(nèi)在的正則化作用(Kaltenbacher et al., 2008).

算法1 截斷牛頓法

給定初始模型m0，令k=0

while反演停止條件未滿足do

whileCG法停止條件未滿足do

計算H(mk)x;

β1=(H(mk)x,x);β2=‖r‖2;

d=d+(β2/β1)x;

r=r+(β2/β1)H(mk)x;

x=-r+(‖r‖2/β2)x;

l=l+1;

endwhile

更新方向pk=x;

計算步長αk;

更新模型mk+1=mk+αkpk;

k=k+1;

endwhile

3 Hessian矩陣性態(tài)分析

Hessian矩陣的性態(tài)會影響到局部迭代優(yōu)化方法的性能.當(dāng)Hessian矩陣的條件數(shù)很大時(病態(tài)矩陣)，會導(dǎo)致局部迭代優(yōu)化方法收斂緩慢甚至不收斂(NocedalandWright，2006).下面我們采用一個簡單的模型來分析參數(shù)Vn和η的Hessian矩陣的結(jié)構(gòu)和性態(tài).

考慮圓形異常模型，如圖1所示.模型網(wǎng)格剖分為51×51，間距均為10m.η真實(shí)模型的四周是寬度為100m的橢圓各向異性層，震源和檢波器置于此層內(nèi)可以壓制S波的影響(PlessixandCao,2011).在模型的每條邊上分別均勻布置12個炮點(diǎn)，炮間距為40m,每炮都采用四周觀測方式接收數(shù)據(jù).震源使用主頻為8Hz的Ricker子波，時間采樣間隔為1ms.初始模型為沒有異常體的背景模型，我們使用有限差分法近似計算出初始模型中心部分(31×31)的Hessian矩陣(NocedalandWright, 2006)，如圖2a所示.

由于考慮了兩類參數(shù)，Hessian矩陣由4個961×961大小的分塊矩陣組成，在圖2a中用黃色虛線進(jìn)行了劃分.左上角塊矩陣中的非對角元素代表兩個不同位置處的Vn對目標(biāo)函數(shù)的影響，對角線上的元素則代表該位置的Vn對目標(biāo)函數(shù)的二階影響，右下角塊矩陣類似地對應(yīng)著參數(shù)η的作用.而非對角塊矩陣(右上和左下)中的元素則描述了不同類別參數(shù)間的相互耦合效應(yīng).圖2b單獨(dú)給出了Hessian矩陣的右上角塊矩陣，它的主對角元素比較突出，說明耦合效應(yīng)主要集中在相同位置的Vn和η上.從圖2a可以看出，右下角塊非常突出.經(jīng)計算Hessian矩陣的條件數(shù)為3.58×1010，病態(tài)比較嚴(yán)重.圖2c給出了該矩陣的特征值分布，該矩陣沒有0特征值(絕對值均大于1×10-5)，是非奇異矩陣.可以看出有負(fù)特征值存在且部分幅度較大，因而該Hessian矩陣不是正定矩陣.若此時采用L-BFGS法進(jìn)行參數(shù)Vn和η的同時反演，由于它構(gòu)造的近似Hessian逆算子為正定矩陣(NocedalandWright, 2006)，可能會收斂緩慢甚至失敗.如果使用TN法則可能較好地處理這類病態(tài)問題.下面我們通過數(shù)值試驗(yàn)進(jìn)行檢驗(yàn).

4 數(shù)值試驗(yàn)

本節(jié)采用兩個不同的二維模型測試TN法的反演性能.每個試驗(yàn)中都將TN法與L-BFGS法(使用最近8代的模型和梯度向量)的收斂速度和反演結(jié)果進(jìn)行了對比分析.

4.1 簡單模型試驗(yàn)

真實(shí)模型包含了圓形(半徑400 m)、正方形(邊長800 m)和三角形(底邊長為800 m的等腰三角形)三個異常體，如圖3所示.η真實(shí)模型的四周各有厚度為100 m的橢圓各向異性層，反演時固定不變.網(wǎng)格剖分為401×401，間距均為10 m.為便于敘述，記該模型為CST模型.在模型的每條邊上分別均勻布置16個炮點(diǎn)，炮間距為240 m,每炮都采用四周觀測方式接收記錄.震源采用主頻為8 Hz的Ricker子波，時間采樣間隔為1 ms.兩類參數(shù)的初始模型均為沒有異常體的背景模型.

圖1 真實(shí)模型 (a) Vn (b) ηFig.1 True model (a) Vn (b) η

圖2 (a) 初始模型的Hessian矩陣及其(c) 特征值分布，(b)為單獨(dú)顯示的Hessian矩陣的右上角塊矩陣Fig.2 (a) Hessian matrix of initial model and its (c) Eigenvalue distribution, (b) Display for the right-upper block matrix of the Hessian

圖3 CST真實(shí)模型 (a) Vn (b) ηFig.3 True CST model (a) Vn (b) η

表1 CST模型試驗(yàn)中TN法和L-BFGS法的對比Table 1 Comparison of TN method and L-BFGS method for the CST model testing

兩種方法的歸一化的目標(biāo)函數(shù)的下降曲線如圖4所示，其橫坐標(biāo)為迭代次數(shù).可以看出，TN法目標(biāo)函數(shù)(紅色實(shí)線)隨迭代次數(shù)的下降速度明顯快于L-BFGS法(黑色實(shí)線)，說明了TN法在單次迭代中能利用更多的Hessian矩陣信息.但從第2部分的原理介紹可知，下降速度的加快需要以增加模擬計算次數(shù)為代價.表1給出了兩種方法在反演過程中的計算數(shù)據(jù)對比.表中第2行數(shù)據(jù)表明在第25次迭代時，TN法的目標(biāo)函數(shù)下降到0.0039，遠(yuǎn)低于L-BFGS法的0.1762，但使用的正演次數(shù)也增加了很多.這反應(yīng)出TN法的總體計算消耗較大.同時表中第1行給出了TN法第6代和L-BFGS法第12代的對比，它們都使用了48次正演.TN法的目標(biāo)函數(shù)已經(jīng)下降到0.052，但L-BFGS法只下降到0.2332，因而在同等計算消耗下，TN法下降更快.

圖4 CST模型的歸一化的目標(biāo)函數(shù)的下降曲線Fig.4 The decrease curve of normalized misfit function for CST model

反演結(jié)果的好壞是評判優(yōu)化方法的關(guān)鍵因素.圖5給出了兩種方法各自對應(yīng)于表1第1行的反演結(jié)果.可以看出，兩組反演結(jié)果差異較大：在3個異常體位置，L-BFGS法對模型η的更新過度，對Vn的更新則很小，有可能收斂不到全局最優(yōu)模型；而TN法對兩類參數(shù)都有較好更新，反演結(jié)果在異常體結(jié)構(gòu)和幅度上都已經(jīng)比較接近真實(shí)模型.接下來進(jìn)一步對比兩種方法的第25代反演結(jié)果，如圖6所示.在TN法的結(jié)果中，Vn模型(圖6c)的3個異常體的結(jié)構(gòu)比第6代時更加清晰，幅度更逼近真實(shí)模型，η模型(圖6d)的幅度比第6代也有細(xì)微改善.對于L-BFGS法，模型Vn的更新仍舊非常緩慢，而模型η偏離真實(shí)模型的程度加重.

為了更清晰地觀察兩種方法的差別，我們對比了反演結(jié)果在不同位置(圖6c中白色點(diǎn)線所示)的模型參數(shù)曲線.圖7給出了對應(yīng)于圖5的反演結(jié)果在深度2km的橫向?qū)Ρ惹€，圖8則是相應(yīng)反演結(jié)果在地表距離1km，2km和3km處的垂向?qū)Ρ惹€.這些曲線也表明TN法第6代的反演結(jié)果，尤其是η模型，已經(jīng)比較接近真實(shí)模型.而L-BFGS法第12代所反演的η模型已較大幅度偏離真實(shí)模型.兩種方法第25代反演結(jié)果(即圖6中模型)的對比曲線分別由圖9和圖10提供.可以看出，TN法的Vn模型曲線已經(jīng)非常接近真實(shí)模型，η模型曲線雖然在異常體的區(qū)分上仍有不足，但在異常體處的幅度與真實(shí)模型基本一致.對于L-BFGS法，迭代次數(shù)的增多并未使反演結(jié)果好轉(zhuǎn)：η模型與真實(shí)模型的偏差增大，Vn模型幾乎沒有更新.這說明L-BFGS法很可能收斂到遠(yuǎn)離真實(shí)模型的極小值點(diǎn)，甚至迭代發(fā)散.

綜合以上對比可知，TN法能夠較好地平衡Vn和η的同時更新，得到了比較準(zhǔn)確的反演結(jié)果.而L-BFGS法則過度集中在一類參數(shù)的更新上，反演結(jié)果較差.由第3部分的分析可知導(dǎo)致兩種方法性能差異的原因可能在于此時Hessian矩陣的性態(tài)比較差.L-BFGS法未能對Hessian逆矩陣進(jìn)行比較準(zhǔn)確的近似，從而沒能準(zhǔn)確解釋兩類參數(shù)間的耦合效應(yīng)，將Vn引起的數(shù)據(jù)誤差轉(zhuǎn)移到了η上，造成其過度更新.這也反應(yīng)了多參數(shù)同時反演的復(fù)雜性.TN法則比L-BFGS法準(zhǔn)確地估計了Hessian逆矩陣，展現(xiàn)出良好性能.

4.2 復(fù)雜模型試驗(yàn)

真實(shí)模型如圖11(a，b)所示，網(wǎng)格剖分為301×71，間距均為20 m.在地表激發(fā)56炮，炮間距為100 m，在地表所有網(wǎng)格點(diǎn)處接收地震記錄.采用主頻為9 Hz的Ricker子波，時間采樣間隔為2.5 ms.對真實(shí)模型進(jìn)行高斯平滑得到兩類參數(shù)的初始模型，如圖11(c，d)所示.試驗(yàn)中也進(jìn)行了L-BFGS法的反演.

兩種方法的歸一化的目標(biāo)函數(shù)隨迭代次數(shù)的下降曲線由圖12給出.類似于CST模型試驗(yàn)，TN法的目標(biāo)函數(shù)(紅色)隨迭代次數(shù)的下降同樣明顯快于L-BFGS法(黑色).圖13分別給出了TN法的最終反演結(jié)果(第46代，正演次數(shù)超過1000，目標(biāo)函數(shù)值為0.0112)與L-BFGS法在相同迭代次數(shù)時(第46代，正演次數(shù)為184，目標(biāo)函數(shù)值為0.35)的反演結(jié)果.可以看出，TN法反演的Vn模型(圖13c)已經(jīng)與真實(shí)模型比較接近，η模型中(圖13d)也得到了比較準(zhǔn)確的界面位置.L-BFSG法則主要對η模型進(jìn)行了過度更新，對Vn模型改善很小.圖14分別給出了兩種方法的反演結(jié)果在地表1.5km，3km和4.5km處(具體位置如圖13c中的白色點(diǎn)線所示)的垂向曲線對比.從中可以更清晰地看出，TN法真實(shí)地反演出了Vn模型的結(jié)構(gòu)和幅度.雖然η模型的幅度與真實(shí)模型仍有差距，但能夠得到比較準(zhǔn)確的界面結(jié)構(gòu).L-BFSG法的模型曲線表明，它對Vn模型的更新非常有限，而主要對η模型的淺部進(jìn)行了過度更新，因此該方法未能收斂到全局最優(yōu)模型.

圖5 L-BFGS法第12代的 (a) Vn (b) η 與TN法第6代的 (c) Vn (d) ηFig.5 Inversion results at the 12th iteration of L-BFGS method (a) Vn (b) η and those at the 6th iteration of TN method (c) Vn (d) η

圖6 第25代反演結(jié)果. L-BFGS法的 (a) Vn (b) η, TN法的 (c) Vn (d) ηFig.6 Inversion results at the 25th iteration . L-BFGS method (a) Vn (b) η and TN method (c) Vn (d) η

圖7 L-BFGS法第12代和TN法第6代反演結(jié)果在深度2 km的橫向曲線對比黑色和紅色實(shí)線分別為真實(shí)和初始模型，藍(lán)色和綠色虛線分別為L-BFGS法和TN法的結(jié)果.Fig.7 Horizontal profile comparisons at depth 2 km between inversion results at the 12th iteration of L-BFGS method (blue dash line) and those at the 6th iteration of TN method (green dash line)The black and red solid lines represent true and initial models, respectively.

圖8 L-BFGS法第12代和TN法第6代反演結(jié)果在地表不同距離的垂向曲線對比黑色和紅色實(shí)線分別為真實(shí)和初始模型，藍(lán)色和綠色虛線分別為L-BFGS和TN法的結(jié)果.Fig.8 Vertical profile comparisons at different surface distances between inversion results at the 12th iteration of L-BFGS method (blue dash line) and those at the 6th iteration of TN method(green dash line)The black and red solid lines represent true and initial models, respectively.

圖9 兩種方法第25代反演結(jié)果在深度2 km處的橫向曲線對比黑色和紅色實(shí)線分別為真實(shí)和初始模型，藍(lán)色和綠色虛線分別為L-BFGS法和TN法的結(jié)果.Fig.9 Horizontal profile comparisons at depth 2 km between inversion results at the 25th iteration of both L-BFGS (blue dash line) and TN method (green dash line)The black and red solid lines represent true and initial models, respectively.

圖10 兩種方法第25代反演結(jié)果在地表不同距離的垂向曲線對比黑色和紅色實(shí)線分別為真實(shí)和初始模型，藍(lán)色和綠色虛線分別為L-BFGS和TN法的結(jié)果.Fig.10 Vertical profile comparisons at different surface distances between inversion results at the 25th iteration of both L-BFGS (blue dash line) and TN method(green dash line) The black and red solid lines represent true and initial models, respectively.

圖11 真實(shí)模型(a) Vn (b) η和初始模型(c) Vn (d) ηFig.11 True models (a) Vn (b) η and initial modes (c) Vn (d) η

以上復(fù)雜模型試驗(yàn)同樣說明了TN法能夠比L-BFGS法更為準(zhǔn)確地估計Hessian逆矩陣，較好地平衡了兩類不同參數(shù)的同時更新.

5 結(jié)論與討論

本文采用TN法開展了時間域聲波VTI介質(zhì)的雙參數(shù)同時反演的研究.地震波速度與各向異性參數(shù)間存在著耦合效應(yīng)，增加了多參數(shù)FWI的復(fù)雜程度.同時，地震波速度和各向異性參數(shù)間取值數(shù)量級的巨大差異也使VTI介質(zhì)多參數(shù)反演問題的性態(tài)變差.合理使用Hessian逆算子能夠減弱這兩類問題對FWI的影響.不同模型的反演試驗(yàn)表明，在地震波速度與各向異性參數(shù)的同時反演中，L-BFGS法對各向異性參數(shù)進(jìn)行了過度更新，對地震波速度的更新量很小，反演效果較差，這也驗(yàn)證了多參數(shù)同

圖12 復(fù)雜模型的歸一化的目標(biāo)函數(shù)下降曲線Fig.12 The decrease curve of normalized misfit function for complex model

時反演問題的復(fù)雜性.而TN法則較好地平衡了這兩類不同參數(shù)的同時更新，得到了比較準(zhǔn)確的反演結(jié)果，這說明了TN法比L-BFGS法能更準(zhǔn)確地估計Hessian逆算子.因而，TN法在多參數(shù)反演中的應(yīng)用潛力較大.

從文中的反演試驗(yàn)可以看出，雖然TN法得到了比較好的反演結(jié)果，對多參數(shù)FWI有重要的理論價值，但由于每步迭代中的CG循環(huán)都需要增加額外的模擬過程，也使方法本身付出了較大的計算代價，這也是限制TN法大規(guī)模實(shí)際應(yīng)用的主要原因所以，采用適當(dāng)?shù)氖侄?如對CG循環(huán)進(jìn)行預(yù)條件，或與其他方法混合迭代等)來提高TN法的計算效率需要進(jìn)一步研究.

多參數(shù)FWI是一項復(fù)雜的研究工作.開發(fā)高效穩(wěn)健的優(yōu)化方法是其中非常重要的一個方面.同時，也需要注重優(yōu)化方法與反演策略的有效結(jié)合.例如，將不同類型參數(shù)的數(shù)量級調(diào)整到相同取值范圍的參

圖13 第46代反演結(jié)果L-BFGS法的(a) Vn (b) η, TN法的 (c) Vn (d) η.Fig.13 Inversion results at the 46th iteration L-BFGS method (a) Vn (b) η and TN method (c) Vn (d) η.

圖14 兩種方法第46代反演結(jié)果在地表不同距離的垂向曲線對比黑色和紅色實(shí)線分別為真實(shí)和初始模型，藍(lán)色點(diǎn)線和綠色虛線分別為L-BFGS和TN法的結(jié)果.Fig.14 Vertical profile comparisons at different surface distances between inversion results at the 46th iteration of both L-BFGS (blue dotted line) and TN method(green dash line) The black and red solid lines represent true and initial models, respectively.

數(shù)標(biāo)準(zhǔn)化策略可以改善Hessian矩陣的性態(tài)，因而可能會加速TN法的收斂.此外，VTI介質(zhì)聲波方程有多種不同的參數(shù)化方式，每種方式的Hessian算子的性態(tài)可能不同.因此，TN法在不同參數(shù)化方式下的性能表現(xiàn)也需要具體研究.

附錄A

pn(t=0)=ph(t=0)=0

?tpn(t=0)=?tph(t=0)=0

(A1)

目標(biāo)函數(shù)對Vn和η的梯度分別為(PlessixandCao，2011)：

(A2)

λn(t=T)=λh(t=T)=0

?tλn(t=T)=?tλh(t=T)=0

(A3)

Akcelik V. 2002. Multiscale Newton-Krylov methods for inverse acoustic wave propagation [Ph. D. thesis]. Pittsburgh: Civil and Environmental Engineering, Carnegie Mellon University.

Alkhalifah T, Plessix R é. 2014. A recipe for practical full-waveform inversion in anisotropic media: An analytical parameter resolution study.Geophysics, 79(3): R91-R101.

Bae H S, Pyun S, Chung W, et al. 2012. Frequency-domain acoustic-elastic coupled waveform inversion using the Gauss-Newton conjugate gradient method.GeophysicalProspecting, 60(3): 413-432.

Brossier R, Operto S, Virieux J. 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion.Geophysics, 74(6): WCC105-WCC118.

Epanomeritakis I, Ak?elik V, Ghattas O, et al. 2008. A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion.InverseProblems, 24(3): 034015.Fichtner A, Kennett B L N, Igel H, et al. 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods.GeophysicalJournalInternational, 179(3): 1703-1725.Fichtner A, Trampert J. 2011. Hessian kernels of seismic data functionals based upon adjoint techniques.GeophysicalJournalInternational, 185(2): 775-798.

Gholami Y, Brossier R, Operto S, et al. 2011. Acoustic VTI full waveform inversion: sensitivity analysis and realistic synthetic examples. ∥81th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2465-2470.

Gholami Y. 2012. Two-dimensional seismic imaging of anisotropic media by full waveform inversion [Ph. D. thesis]. Nice: Université de Nice-Sophia Antipolis.

Guitton A, Symes W W. 2003. Robust inversion of seismic data using the Huber norm.Geophysics, 68(4): 1310-1319. Kaltenbacher B, Neubauer A, Scherzer O. 2008. Iterative Regularization Methods for Nonlinear Ill-posed Problems. New York: Walter de Gruyter.

Lee H -Y, Koo J M, Min D -J, et al. 2010. Frequency-domain elastic full waveform inversion for VTI media.GeophysicalJournalInternational, 183(2): 884-904.

Liu G F, Liu H, Meng X H, et al. 2012. Frequency-related factors analysis in frequency domain waveform inversion.ChineseJ.Geophys. (in Chinese), 55(4): 1345-1353, doi: 10.6038/j.issn.0001-5733.2012.04.030.

McGarry R, Moghaddam P. 2009. NPML boundary conditions for second-order wave equations. ∥79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3590-3594.

Métivier L, Brossier R, Virieux J, et al. 2013. Full waveform inversion and the truncated Newton method.J.Sci.Comput., 35(2): B401-B437.

Métivier L, Brossier R, Operto S, et al. 2014a. Multi-parameter FWI- an illustration of the Hessian operator role for mitigating trade-off between parameter classes. 6th EAGE Saint Petersburg International Conference & Exhibition. Saint Petersburg, Russia.Métivier L, Bretaudeau F, Brossier R, et al. 2014b. Full waveform inversion and the truncated newton method: quantitative imaging of complex subsurface structures.GeophysicalProspecting, 62(6): 1353-1375.Nash S G, Nocedal J. 1991. A numerical study of the limited memory BFGS method and the truncated Newton method for large-scale optimization.SIAMJ.Optim., 1(3): 358-372.

Nocedal J, Wright S J. 2006. Numerical optimization. New York: Springer.

Operto S, Brossier R, Gholami Y, et al. 2013. A guided tour of multiparameter full waveform inversion for multicomponent data: from theory to practice. The Leading Edge,32(9):1040-1054.

Plessix R é. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications.GeophysicalJournalInternational, 167(2): 495-503.

Plessix R é, Cao Q. 2011. A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium.GeophysicalJournalInternational, 185(1): 539-556.Pratt R G, Shin C, Hicks. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362.

Santosa F, Symes W W. 1988. Computation of the Hessian for least-squares solutions of inverse problems of reflection seismology.InverseProblems, 4(1): 211-233.Shi Y M, Zhang Y, Yao F C, et al. 2014. Methodology of seismic imaging for hydrocarbon reservoirs based on acoustic full waveform inversion.ChineseJ.Geophys. (in Chinese), 57(2): 607-617, doi: 10.6038/cjg20140224.Shin C, Jang S, Min D -J. 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory.GeophysicalProspecting, 49(5): 592-606.

Thomsen L. 1986. Weak elastic anisotropy.Geophysics, 51(10):1954-1966.

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6): WCC1-WCC26.

Wang Y, Dong L G, Liu Y Z. 2013. Improved hybrid iterative optimization method for seismic full waveform inversion.AppliedGeophysics, 10(3): 265-277.

Zhang M G, Wang M Y, Li X F, et al. 2003. Full wavefield inversion of anisotropic elastic parameters in the time domain.ChineseJ.Geophys. (in Chinese), 46(1): 94-100.

附中文參考文獻(xiàn)

劉國峰, 劉洪, 孟小紅等. 2012. 頻率域波形反演中與頻率相關(guān)的影響因素分析. 地球物理學(xué)報, 55(4): 1345-1353, doi: 10.6038/j.issn.0001-5733.2012.04.030.

石玉梅, 張研, 姚逢昌等. 2014. 基于聲學(xué)全波形反演的油氣藏地震成像方法. 地球物理學(xué)報, 57(2): 607-617, doi: 10.6038/cjg20140224.

張美根, 王妙月, 李小凡等. 2003. 時間域全波場各向異性彈性參數(shù)反演. 地球物理學(xué)報, 46(1): 94-100.

(本文編輯 汪海英)

Multi-Parameter full waveform inversion for acoustic VTI media using the truncated Newton method

WANG Yi, DONG Liang-Guo*

StateKeyLaboratoryofMarineGeology,TongjiUniversity,Shanghai200092,China

It is important to reconstruct both seismic velocities and anisotropy parameters by anisotropic full waveform inversion (FWI). But the trade-off between different parameter classes increases its nonlinearity severely. The significant difference of magnitude orders between seismic velocities and anisotropy parameters makes the inversion problem poorly conditioned. Judiciously incorporating the inverse Hessian operator within the FWI scheme can help to alleviate the negative impacts of these two problems on FWI. The truncated Newton (TN) method can better estimate the inverse Hessian operator. Therefore, we use the TN method to tackle these two problems and simultaneously reconstruct both seismic velocities and anisotropy parameters for acoustic VTI media.The TN method has two advantages. On one hand, it does not need explicit Hessian matrix, which is prohibitive to build due to high computational cost and large memory requirements. To get the updated direction, the TN method solves Newton equation at each time of iteration with the conjugate gradient method (CG), which only uses Hessian-vector products. We compute these products by the finite difference method. On the other hand, the TN method can adjust flexibly the updated direction when the Hessian matrix is non-positive definite. The CG method at each step of the TN iteration is terminated as soon as a direction of negative curvature is generated. With this strategy the TN method may perform better than Limited-memory Broyden-Fletcher-Goldfarb-Shanno method (L-BFGS) for ill-conditioned problems.We test the inversion capability of the TN method with two numerical experiments. In each experiment, the NMO velocity and anisotropy parameter eta are inverted simultaneously. The L-BFGS method is also tested for comparison. The first experiment is on a simple 2D VTI homogenous model with three different anomalies inside. The inversion is performed in the time domain with a full acquisition system. Inversion results show that the TN method can balance the updates of the two different parameters, and reconstruct more accurate models. In contrast, the L-BFGS method overestimates parameter eta but underestimates NMO velocity. The second experiment is on a more complex 2D model with a surface acquisition system. The two methods perform the same as they do in the first experiment.The reason for their different performances may be found from the Hessian analysis part that we have finished before the two experiments. We design a small size (51×51) inclusion model and compute explicitly the Hessian matrix of the background initial model. The analyses show that the Hessian matrix has several negative eigenvalues and a large condition number (ill-conditioned matrix), which may result in the failure of the L-BFGS method. The TN method better overcomes these difficulties. Hence, in the context of simultaneous inversion of two parameters for acoustic VTI media, the TN method can account for more accurate inverse Hessian operator than the L-BFGS method and get more accurate inversion results. It is a promising method for multi-parameter FWI.

Multi-Parameter；Full waveform inversion；Anisotropy；Trade-off；Inverse Hessian operator；Truncated Newton method

國家油氣重大專項(2011ZX05005-005-007HZ)資助.

王義，男，1985年生，博士，研究方向?yàn)榈卣鹑ㄐ畏囱莸睦碚摵蛻?yīng)用.

*通訊作者董良國，男，1966年生，同濟(jì)大學(xué)教授，博士生導(dǎo)師，主要從事地震波傳播與數(shù)值模擬、地震波波形反演與層析成像等方面的研究.E-mail:dlg@#edu.cn

10.6038/cjg20150821.

10.6038/cjg20150821

P631

2014-08-17，2015-06-16收修定稿

王義，董良國. 2015. 基于截斷牛頓法的VTI介質(zhì)聲波多參數(shù)全波形反演.地球物理學(xué)報，58(8):2873-2885,

Wang Y, Dong L G. 2015. Multi-Parameter full waveform inversion for acoustic VTI media using the truncated Newton method.ChineseJ.Geophys. (in Chinese),58(8):2873-2885,doi:10.6038/cjg20150821.