基于截?cái)嗯ｎD法的VTI介質(zhì)聲波多參數(shù)全波形反演

王義， 董良國(guó)

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

?

基于截?cái)嗯ｎD法的VTI介質(zhì)聲波多參數(shù)全波形反演

王義， 董良國(guó)*

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

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

多參數(shù)； 全波形反演； 各向異性； 耦合； Hessian逆算子； 截?cái)嗯ｎD法

1 引言

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

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

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

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

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

2 截?cái)嗯ｎD法原理

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

(1)

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

mk+1=mk+αkpk，

(2)

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

(3)

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

(4)

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

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

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

算法1 截?cái)嗯ｎD法

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

while反演停止條件未滿(mǎn)足do

whileCG法停止條件未滿(mǎn)足do

計(jì)算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;

計(jì)算步長(zhǎng)αk;

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

k=k+1;

endwhile

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

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

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

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

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

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

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

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

圖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法的對(duì)比Table 1 Comparison of TN method and L-BFGS method for the CST model testing

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

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

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

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

綜合以上對(duì)比可知，TN法能夠較好地平衡Vn和η的同時(shí)更新，得到了比較準(zhǔn)確的反演結(jié)果.而L-BFGS法則過(guò)度集中在一類(lèi)參數(shù)的更新上，反演結(jié)果較差.由第3部分的分析可知導(dǎo)致兩種方法性能差異的原因可能在于此時(shí)Hessian矩陣的性態(tài)比較差.L-BFGS法未能對(duì)Hessian逆矩陣進(jìn)行比較準(zhǔn)確的近似，從而沒(méi)能準(zhǔn)確解釋兩類(lèi)參數(shù)間的耦合效應(yīng)，將Vn引起的數(shù)據(jù)誤差轉(zhuǎn)移到了η上，造成其過(guò)度更新.這也反應(yīng)了多參數(shù)同時(shí)反演的復(fù)雜性.TN法則比L-BFGS法準(zhǔn)確地估計(jì)了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子波，時(shí)間采樣間隔為2.5 ms.對(duì)真實(shí)模型進(jìn)行高斯平滑得到兩類(lèi)參數(shù)的初始模型，如圖11(c，d)所示.試驗(yàn)中也進(jìn)行了L-BFGS法的反演.

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

5 結(jié)論與討論

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

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

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

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

多參數(shù)FWI是一項(xiàng)復(fù)雜的研究工作.開(kāi)發(fā)高效穩(wěn)健的優(yōu)化方法是其中非常重要的一個(gè)方面.同時(shí)，也需要注重優(yōu)化方法與反演策略的有效結(jié)合.例如，將不同類(lèi)型參數(shù)的數(shù)量級(jí)調(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é)果在地表不同距離的垂向曲線(xiàn)對(duì)比黑色和紅色實(shí)線(xiàn)分別為真實(shí)和初始模型，藍(lán)色點(diǎn)線(xiàn)和綠色虛線(xiàn)分別為L(zhǎng)-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)，因而可能會(huì)加速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ù)對(duì)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)

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

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

張美根, 王妙月, 李小凡等. 2003. 時(shí)間域全波場(chǎng)各向異性彈性參數(shù)反演. 地球物理學(xué)報(bào), 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

國(guó)家油氣重大專(zhuān)項(xiàng)(2011ZX05005-005-007HZ)資助.

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

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

10.6038/cjg20150821.

10.6038/cjg20150821

P631

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

王義，董良國(guó). 2015. 基于截?cái)嗯ｎD法的VTI介質(zhì)聲波多參數(shù)全波形反演.地球物理學(xué)報(bào)，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.