王 義，董良國

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

地震全波形反演(full waveform inversion，F(xiàn)WI)是近年來勘探地球物理領(lǐng)域的研究熱點(diǎn)之一。它利用疊前地震記錄的波形信息來重建高精度的地下介質(zhì)的物性參數(shù)(密度、速度和模量等)，可以為油氣勘探中的偏移成像提供更精確的背景速度[1]，可以應(yīng)用到儲層定量研究領(lǐng)域[2]，還可以應(yīng)用到全球局部和區(qū)域地殼和上地幔結(jié)構(gòu)等研究中[3]。因而，在油氣勘探開發(fā)和地球動力學(xué)等領(lǐng)域有著廣闊的應(yīng)用前景。

20世紀(jì)80年代，Lailly[4]和Tarantola[5]利用最小二乘優(yōu)化的思想在時間域建立了FWI的核心框架。FWI希望找到最優(yōu)的地下介質(zhì)模型參數(shù)，使模擬地震記錄與觀測記錄達(dá)到最佳匹配。FWI涉及的數(shù)據(jù)量大、模型參數(shù)多，是典型的大規(guī)模優(yōu)化問題，通常采用局部迭代優(yōu)化方法求解[6]。目標(biāo)函數(shù)對模型參數(shù)的二階導(dǎo)數(shù)，即Hessian矩陣，對加速迭代方法的收斂起著重要作用[7]。但巨大的計算和存儲消耗使顯式獲取Hessian及其逆矩陣變得非常困難。克服這個困難通常有三類方法:預(yù)條件共軛梯度法(preconditioned conjugate gradient method,PCG)、不精確Newton法和擬Newton法[6]。PCG法在FWI中獲得了大量應(yīng)用[8]，它通過使用預(yù)條件矩陣對梯度進(jìn)行校正，來加速迭代的收斂過程[9-13]。而預(yù)條件矩陣通常是先對Hessian矩陣做一些簡單近似(如對角近似或塊對角近似)，然后求逆得到預(yù)條件矩陣。不精確Newton法在每步迭代中通過采用PCG法近似求解Newton方程來獲取本次迭代的更新方向。如果沒有好的預(yù)條件矩陣，每步迭代中的PCG 求解過程收斂很慢，導(dǎo)致不精確Newton法的整體計算效率偏低[6]。擬Newton法是求解無約束優(yōu)化問題的所有利用一階導(dǎo)數(shù)信息的方法中最有效的一類方法[14]，其中的L-BFGS法在FWI中得到了廣泛應(yīng)用[15-17]。

L-BFGS法不需要顯式形成Hessian及其逆矩陣，在非線性優(yōu)化方面表現(xiàn)出比CG法收斂快且穩(wěn)健的優(yōu)勢[18]。該方法在每步迭代中都需要提供Hessian逆矩陣的一個初始近似矩陣。這個矩陣可以隨迭代的進(jìn)行而更新，也可以固定不變。初始矩陣的選取和更新方式對L-BFGS法的收斂性能有重要影響。但初始矩陣的選取并沒有統(tǒng)一標(biāo)準(zhǔn)，需要根據(jù)具體優(yōu)化問題而定[6]。如前文所述，PCG法采用的預(yù)條件矩陣即是對Hessian逆矩陣的近似。因而，選擇初始矩陣時可以借鑒預(yù)條件矩陣的作法。此外，L-BFGS法需要存儲一定數(shù)目的梯度殘差和模型殘差向量對，存儲數(shù)目的多少也會對其收斂性能產(chǎn)生影響。

針對上述問題，我們分析了幾種時間域可行的預(yù)條件矩陣的原理，對采用它們作為初始矩陣的L-BFGS法的反演性能進(jìn)行了對比分析，并與相應(yīng)的PCG法比較。此外，還分析了存儲向量對的數(shù)目對L-BFGS法的影響。最后通過時間域的Marmousi模型FWI反演試驗(yàn)，來說明不同的初始矩陣選擇以及更新方式對L-BFGS法收斂性能的影響。

1 時間域全波形反演理論簡介

時間域波動方程的離散形式為

(1)

式中：u,S分別為波場和震源;m為模型參數(shù)，維數(shù)M;A(m)為正演矩陣?？梢钥闯?，A-1的每一列都對應(yīng)于特定震源位置和特定激發(fā)時刻的Green函數(shù)。

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

(2)

式中：ns為炮數(shù)；ds是單炮觀測記錄；Rs為提取記錄的算子。Newton法的更新方向滿足

(3)

(4)

下面采用虛震源給出Jacobi矩陣的構(gòu)造原理。分別在(1)式兩邊對m求導(dǎo)可得

(5)

定義F為虛震源，它的每一列代表該位置處模型參數(shù)的虛震源。虛震源中除了包含波場的作用外，還涉及到傳播算子對模型參數(shù)的影響?A/?m，描述了數(shù)據(jù)對模型參數(shù)的隨方位角變化的敏感性[19]。以F的每列為震源分別做正演模擬可以得到偏導(dǎo)波場，進(jìn)而構(gòu)造出Jacobi矩陣。因而，Gauss-Newton的近似Hessian矩陣可以寫為

(6)

2 L-BFGS方法原理和初始矩陣的幾種選擇方法

L-BFGS法利用目標(biāo)函數(shù)的一階導(dǎo)數(shù)信息遞歸地計算近似Hessian逆矩陣與梯度的乘積，來得到迭代更新方向，克服了顯式求取Hessian矩陣及其逆的困難，在FWI中應(yīng)用廣泛。

2.1 L-BFGS優(yōu)化方法原理

每步迭代中近似Hessian逆矩陣與梯度向量乘積的具體計算公式為

(7)

(8)

這種初始矩陣在每次迭代中都可以更新[6]。以上2種方式并沒有具體結(jié)合FWI中Hessian矩陣的特點(diǎn)。下面我們討論幾種針對時間域FWI的初始矩陣的選擇方式。

2.2 L-BFGS法中初始矩陣的幾種選擇方式

預(yù)條件矩陣通常是由對Hessian矩陣做簡單的近似，如對角近似或塊對角近似，再求逆后得到。因而，初始矩陣的選取可以借鑒預(yù)條件矩陣的做法。在頻率域也可以對Hessian矩陣做塊對角近似[21]，但在時間域計算塊對角Hessian矩陣的運(yùn)算量很大，所以一般并不可行。下面分析幾種在時間域可行的對角預(yù)條件矩陣方式。

2.2.1 Pseudo Hessian(PH方式)

Shin等人[9]定義FTF為Pseudo Hessian，并用對角矩陣diag(FTF)來近似Hessian矩陣。從虛震源F的定義可知，這種方式只考慮了震源波場的作用和傳播算子對模型參數(shù)的影響。它用單位矩陣替代(6)式的矩陣B，忽略了所有檢波點(diǎn)處的Green函數(shù)對目標(biāo)函數(shù)的影響。

2.2.2 New Pseudo Hessian(newPH方式)

這種方式是對上一種方式的改進(jìn)，由Choi等人在頻率域提出[10]。代替矩陣B的仍是對角矩陣，對角元素的值則是所有震源位置的Green函數(shù)在該空間點(diǎn)處的波場的振幅之和。如果觀測系統(tǒng)中炮點(diǎn)和檢波點(diǎn)有重合的位置，則這種方式能考慮重合部分的檢波點(diǎn)處的Green函數(shù)的作用。在時間域可以采用類似的做法，用所有震源波場在地下同一點(diǎn)的振幅之和，作為代替矩陣B的對角矩陣在該位置的對角元素。

2.2.3 采用震源正傳波場和記錄殘差反傳波場的能量矩陣(RcdIW方式)

該方式采用震源正傳波場與以模擬和觀測記錄的殘差作為震源的反傳波場的能量分布來構(gòu)造梯度的預(yù)條件矩陣，目的是消除正、反傳播過程中的幾何擴(kuò)散效應(yīng)[11]。記u(x;xs)為xs處的震源的正傳波場，λ(x;xr)為記錄殘差反傳的波場，那么，正傳波場的能量矩陣為

(9)

反傳波場的能量矩陣為

(10)

對梯度進(jìn)行預(yù)條件:

(11)

嚴(yán)格地講，這種方式并不是對Hessian逆矩陣的近似，而是從波場能量分布的角度來對梯度進(jìn)行校正。它考慮了記錄殘差反傳波場的作用，但也攜帶了散射點(diǎn)的信息。采用它來預(yù)條件梯度會破壞模型擾動的聚焦。另外，該方式?jīng)]有考慮傳播算子對模型參數(shù)的影響。

2.2.4 采用震源正傳波場和模擬記錄反傳波場的能量矩陣(ResIW方式)

該方法是對上一種方式的改進(jìn)[12]，不同之處在于采用的是模擬記錄反傳波場的能量分布，不會對模型擾動的聚焦造成破壞。但該方式同樣沒有考慮傳播算子對模型參數(shù)的影響。它需要增加一次波場模擬來計算模擬記錄的反傳波場v(x;xr)。正傳波場的能量矩陣仍采用(9)式。反傳波場的能量矩陣以及預(yù)條件后的梯度分別為

2.2.5 用Hessian-vector 乘積構(gòu)造對角矩陣

Korta等[13]在雙參數(shù)FWI中采用這種方式構(gòu)造塊狀對角的Hessian近似。在單參數(shù)FWI中，可以使用有限差分法[6]、二階伴隨狀態(tài)法[22]或基于Born核函數(shù)的矩陣分解方法計算出Hessian-vector乘積Hv，其中v=(1,1,…,1,1)T。以該乘積為對角元素構(gòu)造一個對角矩陣。此時，每個對角元素均為Hessian矩陣該行所有元素的疊加。如

果Hessian矩陣是嚴(yán)格對角占優(yōu)的話，即對角線元素的絕對值大于同一行所有其它元素的絕對值之和，經(jīng)過按行疊加得到的元素會偏離Hessian矩陣的主對角元素的精確值。如果不是嚴(yán)格對角占優(yōu)，計算出的乘積很有可能是零或負(fù)值。即使取絕對值的話，數(shù)值試驗(yàn)表明其反演過程容易出現(xiàn)不穩(wěn)定的現(xiàn)象。所以，在單參數(shù)FWI中，這種方式并不是合適的預(yù)條件矩陣。

3 數(shù)值反演試驗(yàn)

我們通過Marmousi模型合成數(shù)據(jù)的FWI反演數(shù)值試驗(yàn)，對比采用上述前4種方式作為初始矩陣的L-BFGS法的反演性能，同時與對應(yīng)的預(yù)條件共軛梯度法進(jìn)行比較。

從原始Marmousi模型中抽取了5940m×2980m的模型，實(shí)際模型(局部)和初始模型如圖1 所示。網(wǎng)格剖分為298×150，網(wǎng)格間距20m。共使用56炮，間距為100m,檢波器298個，間距為20m。時間采樣間隔2ms，記錄長度3s。震源采用主頻7Hz的Ricker子波。采用空間4階、時間2階有限差分法求解常密度聲波方程，四周均為PML吸收邊界[23]。

圖1 實(shí)際Marmousi模型(a)和初始模型(b)

3.1 殘差向量對的數(shù)目對L-BFGS性能的影響

L-BFGS法構(gòu)造近似Hessian逆矩陣時，需要存儲r對梯度殘差和模型殘差。r越大，構(gòu)造近似Hessian矩陣用到的信息越豐富，但所需內(nèi)存也隨之增大。假設(shè)模型向量的長度為M，則存儲這些殘差向量需要2M×r個浮點(diǎn)類型(float類型)的內(nèi)存。首先，我們測試了L-BFGS法中需要存儲的向量對的數(shù)目r對其性能的影響。

分別取r為 4，8，12，16和20，初始矩陣采用公式(8)的方式，對應(yīng)的目標(biāo)函數(shù)收斂曲線如圖2。全波形反演的主要計算消耗為正演模擬，所以用目標(biāo)函數(shù)隨正演次數(shù)的變化來描述優(yōu)化方法的收斂速度。可以看出，隨著r的增大，L-BFGS法的收斂速度也加快。r=4時收斂最慢，r為8～20時收斂速度逐漸加快，但提升并不明顯。當(dāng)r=20時收斂最快，但存儲殘差向量的內(nèi)存是r=8時的2.5倍。圖3分別給出了r=8時的第324代(972次正演)和r=20時的第291代(873次正演)反演結(jié)果。此時，二者的目標(biāo)函數(shù)均下降到1.01×10-3附近。圖4分別對比了圖3中反演結(jié)果在地表3km和5km處的縱向速度曲線。從圖3 和圖4可以看出，2個反演結(jié)果非常接近。因而，綜合考慮收斂速度和內(nèi)存消耗，參數(shù)r的取值無需過大。在此選取r=8。

圖2 存儲殘差向量對的數(shù)目r對L-BFGS法目標(biāo)函數(shù)收斂的影響

圖3 存儲不同向量對數(shù)目r的L-BFGS法的反演結(jié)果a r=8,第324代; b r=20,第291代

圖4 對應(yīng)圖3中反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對比

3.2 不同初始矩陣情況下的L-BFGS法的反演性能對比

我們分別選擇前4種方式作為L-BFGS法的初始矩陣,即依次為PH,newPH,RcdIW和ResIW。初始矩陣在迭代過程中可以固定不變，也可以每次迭代都更新(Up1)。另外，一并測試單位矩陣(Identity)和公式(8)對應(yīng)的矩陣(Nocedal)的情況。圖5給出了不同初始矩陣、不同更新方式情況下的L-BFGS法的目標(biāo)函數(shù)收斂曲線。

從圖5可以看出，不同初始矩陣、不同更新方式的L-BFGS法的反演性能存在較大差異。按收斂速度快慢大致可以分為3類。

第一類收斂最快，包含以下5種方式。初始矩陣選取PH和newPH且固定不變時(藍(lán)色和粉色虛線)收斂最快，二者非常接近。選RcdIW且固定不變時(綠色虛線)和選PH每次迭代都更新時(藍(lán)色點(diǎn)劃線)也比較接近，相比于前兩種方式則稍微慢一些。選newPH且每代都更新時(粉色點(diǎn)劃線)在開始階段較慢，之后收斂加快。

選RcdIW且每代都更新時(綠色點(diǎn)劃線)與選ResIW固定不變時(黑色虛線)收斂速度比較接近，屬于第二類收斂較快的方式。

選ResIW且每代都更新(黑色點(diǎn)劃線)，Nocedal方式(紅色實(shí)線)和Identity(青色實(shí)線)屬于收斂比較慢的第三類。其中，單位矩陣并未提供有關(guān)Hessian的信息，收斂最慢。

綜合以上對比結(jié)果可知，PH和newPH作初始矩陣且固定不變時，都可以快速收斂。RcdIW方式稍慢于前兩者，可能與其未考慮正演算子對模型參數(shù)的影響有關(guān)。因存在不利于模型擾動聚焦的因素，ResIW在4種方式中收斂最慢。另外一個結(jié)論是，每種方式固定不變時(虛線)比每次迭代都更新(點(diǎn)劃線)的收斂要快。這說明每次迭代更新初始矩陣所增加的計算消耗并未能帶來目標(biāo)函數(shù)的快速下降。

下面對比反演所得的速度模型。表1列出了

目標(biāo)函數(shù)(歸一化后)下降到1.01×10-3左右時，圖5中各種方式的L-BFGS法所需要的迭代次數(shù)和正演次數(shù)的對比，其中單位矩陣(Identity)方式收斂最慢，只下降到1.202×10-3。圖6是PH方式第88代和RcdIW方式第97代的反演結(jié)果，對應(yīng)的目標(biāo)函數(shù)值分別為1.016×10-3和1.011×10-3,正演次數(shù)則為264次和292次。圖7給出了這兩個模型在地表3km和5km處的縱向速度曲線的對比?？梢钥闯觯瑑煞N方式的反演結(jié)果非常接近，與真實(shí)模型也比較吻合。右下方包含高速體的區(qū)域則都存在偏差，對真實(shí)Hessian矩陣逼近的程度不夠可能是原因之一。圖7中還給出了PH方式第34代(102次正演)反演的速度曲線(粉色虛線)，對應(yīng)的目標(biāo)函數(shù)值為1.03×10-3。此時的反演模型與真實(shí)模型有一定偏差，說明仍需迭代更新。從這些對比可以看出，當(dāng)目標(biāo)函數(shù)下降到基本相同的水平時，PH和RcdIW方式的反演結(jié)果也非常接近。而且，PH方式充當(dāng)初始矩陣時在計算效率上占有少量優(yōu)勢。

圖5 不同初始矩陣和不同更新方式的L-BFGS法的目標(biāo)函數(shù)收斂曲線

NocedalPHnewPHRcdIWResIW目標(biāo)函數(shù)值1.017×10-31.016×10-31.039×10-31.011×10-31.018×10-3正演次數(shù)966264250292507迭代次數(shù)322888397169IdentityPH-Up1newPH-Up1RcdIW-Up1ResIW-Up1目標(biāo)函數(shù)值1.202×10-31.015×10-31.007×10-31.019×10-31.011×10-3正演次數(shù)999336336464849迭代次數(shù)33311284116283

圖6 PH方式第88代(a)和RcdIW方式第97代(b)反演的速度模型

圖8對比了PH方式第88代和newPH方式第83代反演結(jié)果的速度曲線。此時它們分別運(yùn)行了264次和250次正演，效率相當(dāng)，反演結(jié)果也非常相似。圖9是PH方式第88代和ResIW方式第169代反演結(jié)果的速度曲線對比，可以看出反演結(jié)果比較接近。但所需正演次數(shù)分別為264次和507次。所以ResIW方式的計算效率偏低。

下面分析初始矩陣每次迭代都更新時的反演結(jié)果。圖10是PH-Up1第112代(336，表示正演次數(shù)，下同)和PH第88代(264)反演結(jié)果的速度曲線對比。圖11是RcdIW-Up1第116代(464)和RcdIW 第97代(292)的對比。可以看出，對于每種方式來說，當(dāng)目標(biāo)函數(shù)下降到基本相同的水平時，初始矩陣固定和每次迭代都更新這兩種做法的反演結(jié)果很接近。但初始矩陣固定時計算效率更高。對于newPH和ResIW兩種方式，結(jié)論相同。

綜上可知，初始矩陣不變時，PH和newPH方式收斂最快，計算效率和反演結(jié)果均相當(dāng);RcdIW方式效率稍慢一些，ResIW方式則收斂較慢。以上4種方式作初始矩陣且每次迭代都更新時，L-BFGS方法的收斂速度要慢于初始矩陣固定不變的情況。因而，選好初始矩陣后，固定不變即可。初始矩陣選Nocedal方式和單位矩陣時收斂很慢。

圖7 對應(yīng)圖6中反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對比

圖8 PH方式和newPH方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對比

圖9 PH方式和ResIW方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對比

圖10 PH-Up1方式和PH方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對比

圖11 RcdIW-Up1方式和RcdIW方式反演結(jié)果在地表3km(a)和5km(b)處的縱向速度曲線對比

此外，L-BFGS法構(gòu)造近似Hessian逆矩陣時使用的是梯度殘差向量。分別嘗試采用PH，newPH，RcdIW和ResIW 4種方式預(yù)條件后的梯度殘差向量來構(gòu)造近似Hessian逆矩陣。如此進(jìn)行的L-BFGS法FWI試驗(yàn)的收斂速度慢于上述作為初始矩陣的情況，而且容易導(dǎo)致反演不穩(wěn)定。所以，不推薦對L-BFGS法中的梯度進(jìn)行額外的預(yù)條件。

3.3 不同預(yù)條件CG法的反演性能對比

分別采用上述4種方式的預(yù)條件共軛梯度(CG)法進(jìn)行同樣的全波形反演試驗(yàn)。圖12給出了4種預(yù)條件CG法與之前不同初始矩陣(固定)L-BFGS法的目標(biāo)函數(shù)收斂曲線的對比?？梢钥闯?，在優(yōu)化過程前期4種預(yù)條件CG法能達(dá)到和L-BFGS法比較接近的收斂速度，但都逐漸變慢。其中，ResIW預(yù)條件CG法收斂最慢。其它3種預(yù)條件CG法也都慢于初始矩陣固定的PH，newPH和RcdIW的L-BFGS法，但優(yōu)于初始矩陣為單位矩陣和Nocedal的方式。所以，在全波形反演中選取適當(dāng)?shù)某跏季仃嚕梢允筁-BFGS法的性能優(yōu)于預(yù)條件CG法。這與Brossier等在頻率域使用PH方式進(jìn)行FWI的結(jié)論一致[16]。

表2列出了不做預(yù)條件的CG法和上述4種預(yù)條件CG法的反演過程中的一組數(shù)據(jù)對比，主要包括目標(biāo)函數(shù)值、所需迭代次數(shù)以及正演次數(shù)。圖13a 給出了PH方式的預(yù)條件CG法第258代反演的速度模型，此時正演次數(shù)為774次，目標(biāo)函數(shù)值約為1.012×10-3。圖13b則給出了RcdIW方式預(yù)條件CG法第203次迭代的結(jié)果，目標(biāo)函數(shù)下降為1.014×10-3，正演次數(shù)為812次。

為了與采用對應(yīng)方式作初始矩陣的L-BFGS法比較(圖6)，我們抽取了相應(yīng)的速度曲線。圖14 顯示了PH方式的預(yù)條件CG法和作為L-BFGS法初始矩陣時反演結(jié)果的速度對比曲線。圖15則是RcdIW方式的對比?？梢钥闯觯瑘D6與圖13對應(yīng)的反演速度模型是非常接近的。從計算消耗上來說，初始矩陣為PH和RcdIW的L-BFGS法分別進(jìn)行了264次和292次正演。可見其計算效率高于相應(yīng)的預(yù)條件CG法。對于newPH方式也可以得到同樣的結(jié)論，這里不再給出具體的圖像和曲線。不做預(yù)條件的CG法和ResIW方式的預(yù)條件CG法的目標(biāo)函數(shù)分別下降到1.604×10-3和1.948×10-3時，已經(jīng)消耗掉約1000次正演，可見其效率偏低。

綜合以上對比分析可以看出，盡管使用預(yù)條件矩陣可以提高CG法的效率，但仍不及采用同樣方式作初始矩陣的L-BFGS法。

圖12 不同預(yù)條件CG法與不同初始矩陣(固定)L-BFGS法的目標(biāo)函數(shù)收斂曲線對比

CGPreCG-PHPreCG-newPHPreCG-RcdIWPreCG-ResIW目標(biāo)函數(shù)值1.604×10-31.012×10-31.038×10-31.014×10-31.948×10-3正演次數(shù)999774698812999迭代次數(shù)333258174203333

圖13 預(yù)條件CG法的反演結(jié)果a PH方式第258代； b RcdIW方式第203代

圖14 PH方式預(yù)條件CG法和充當(dāng)H0且固定不變時的L-BFGS法的反演結(jié)果在地表3km(a)和5km(b)處的vP曲線

圖15 RcdIW方式預(yù)條件CG法和充當(dāng)H0且固定不變的L-BFGS法的反演結(jié)果在地表3km(a)和5km(b)處的vP曲線

4 結(jié)束語

我們討論了全波形反演中L-BFGS法的有效使用問題，主要分析了初始矩陣的選取方式和更新方式對L-BFGS法性能的影響。通過時間域Marmousi模型反演試驗(yàn)的比較，可以得出幾點(diǎn)結(jié)論。

1) 初始矩陣固定不變時，只采用虛震源的PH方式和增加了震源波場振幅的newPH方式的L-BFGS方法收斂速度相當(dāng)且最快。采用震源正傳波場和模擬記錄反傳波場能量的RcdIW方式收斂稍慢。采用震源正傳波場和記錄殘差反傳波場能量的ResIW方式收斂明顯慢于前3種方式，而使用Nocedal方式和單位矩陣的L-BFGS法收斂最慢。另外，在單參數(shù)FWI中，采用Hessian向量乘積構(gòu)造對角矩陣的方式并不適合作為初始矩陣。

2) PH，newPH，RcdIW和ResIW4種方式作為初始矩陣且固定不變的L-BFGS法的收斂速度快于同種方式但初始矩陣每次迭代都更新的情形。

3) 與預(yù)條件CG法相比，采用同種方式作初始矩陣的L-BFGS方法收斂更快。

4) L-BFSG法的收斂速度會隨存儲的殘差向量對的數(shù)目的增大而有所提升，但所需內(nèi)存也會明顯增多。因而綜合計算效率和內(nèi)存消耗的考慮，存儲數(shù)目無需過大(如8左右)。

參 考 文 獻(xiàn)

[1] Boonyasiriwat C,Schuster G T,Paul V,et al.Applications of multiscale waveform inversion to marine data using a flooding technique and dynamic early-arrival windows[J].Geophysics,2010,75(6):R129-R136

[2] 單蕊，卞愛飛，於文輝，等.利用疊前全波形反演進(jìn)行儲層預(yù)測[J].石油地球物理勘探,2011,46(4):629-633

Shan R,Bian A F,Yu W H,et al.Pre-stack full waveform inversion for reservoir prediction[J].Oil Geophysical Prospecting,2011,46(4):629-633

[3] Fichtner A,Kennett B L N,Igel H,et al.Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods[J].Geophysical Journal International,2009,179(3):1703-1725

[4] Lailly P.The seismic inverse problem as a sequence of before stack migrations[C]∥Conference on inverse scattering:theory and application.NewYork:SIAM,1983:206-220

[5] Tarantola A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[6] Nocedal J,Wright S J.Numerical optimization[M].2nd ed.New York:Springer,2006:101-184

[7] Pratt R G,Shin C,Hicks G J.Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

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

[9] Shin C,Jang S,Min D J.Improved amplitude preservation for prestack depth migration by inverse scattering theory[J].Geophysical Prospecting,2001,49(5):592-606

[10] Choi Y,Min D J,Shin C.Frequency-domain elastic full waveform inversion using the new pseudo-Hessian matrix:experience of elastic Marmousi-2 synthetic data[J].Bulletin of the Seismological Society of America,2008,98(5):2402-2415

[11] Zhang Z,Lin Y,Huang L.Full-waveform inversion in the time domain with an energy weighted gradient[J].Expanded Abstracts of 81stAnnual Internet SEG Mtg,2011,2772-2776

[12] Zhang Z,Lin Y,Huang L.A wave-energy-based precondition approach to full-waveform inversion in the time domain[J].Expanded Abstracts of 82ndAnnual Internet SEG Mtg,2012,1-5

[13] Korta N,Fichtner A,Sallaers V.Block-diagonal approximate hessian for Preconditioning in full waveform inversion[J].Expanded Abstracts of 75thEAGE Conference,2013,Tu-P04-16

[14] 袁亞湘.非線性優(yōu)化計算方法[M].北京:科學(xué)出版社，2008:82-88

Yuan Y X.Computation Methods for Nonlinear Optimization[M].Beijing:Science Press,2008:82-88

[15] Guitton A,Symes W W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1319

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

[17] Prieux V,Brossier R,Operto S,et al.Multi-parameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field.Part 1:imaging compressional wavespeed,density and attenuation[J].Geophysical Journal International,2013,194(3):1640-1664

[18] Kelley C T.Iterative methods for optimization[M].New York:SIAM,1999:71-83

[19] Gholami Y,Ribodetti A,Brossier R,et al.Sensitivity analysis of full waveform inversion in VTI media[J].Expanded Abstracts of 72ndEAGE Conference,2010,P371

[20] Plessix R E,Mulder W A.Frequency-domain finite-difference amplitude preserving migration[J].Geophysical Journal International,2004,157(3):975-987

[21] Hu W,Abubabak A,Habashy T M,et al.Preconditioned non-linear conjugate gradient method for frequency domain full waveform seismic inversion[J].Geophysical prospecting,2011,59(3):477-491

[22] Fichtner A,Trampert J.Hessian kernels of seismic data functionals based upon adjoint techniques[J].Geophysical Journal International,2009,185(2):775-798

[23] Berenger J P.A perfectly matched layer for the absorption of electromagnetic waves[J].Journal of Computational Physics，1994,114(2):185-200