張文生，羅嘉，滕吉文

1中國科學院數(shù)學與系統(tǒng)科學研究院，計算數(shù)學與科學工程計算研究所，LSEC，北京100190

2中國科學院地質與地球物理研究所，北京100029

1 引言

地震全波形反演是利用地表或鉆孔中觀測到的疊前波場記錄信息，推測地球內部介質的物性參數(shù)，以確定含油氣構造.速度是重要的物性參數(shù)之一，速度反演或建模是地球物理偏移成像和地震資料解釋的重要基礎.走時反演或成析成像（例如，F(xiàn)awcett et al.，1985；Dahlen et al.，2000；Hung et al.，2004；Gautier et al.，2008；Tong et al.，2011）是速度反演的一種重要方法，但該方法只利用了波場的到時信息，本質上是一種高頻近似方法，適于對區(qū)域背景速度場成像，對復雜速度構造模型成像精度不足.全波形反演（例如，Tarantola，1984；Mora，1987；Yang，1993；Pratt et al.，1998）利用波場的走時、振幅和相位信息來反演，具有更好的精度和更高的分辨率，但計算量大，難度大，特別是二維和三維情形.從20世紀80年代以來，人們一直致力于這方面的研究，發(fā)展了多種地震全波形反演方法，這些方法相互交叉或綜合，并無嚴格界限，大致可如下分類.基于求解波動方程計算域的不同，可分時間域反演（例如，Tarantola，1984；Bleistein，1987；Zhang et al.，2013）、頻率域反 演 （Pratt et al.，1988；Song，1995；Pratt et al.，1998；Prat，1999；Sirgue et al.，2004）和Laplace域反演（Shin et al.，2008，2009；Ha et al.，2012）.基于求解方式的不同，有局部線性化梯度類迭代反演（Tarantola，1984，1987；Gauthier et al.，1986；Bleistein，1987；Mora，1987）、混沌反演（Yang，1993）和全局搜索優(yōu)化反演等.

線性反演不能恢復速度模型的長波場分量（Devancy，1984），波動方程反問題的主要方法是非線性優(yōu)化迭代反演方法和Born反演（Bleistein，1987）.Tarantola首先 （Tarantola，1984，1986）對完全非線性全波形反演進行了詳細的研究，數(shù)學上歸結為一個局部最優(yōu)化問題.求解這類問題時，擾動模型沿目標函數(shù)的梯度搜索，其中梯度可由震源激發(fā)的入射波場和殘差反傳播后的波場計算得到，這種方法可以避免直接通過模型參數(shù)擾動來計算Frechet導數(shù)，大大減少了計算量，使得二維時間域波動方程全波形反演成為可能，該思想已用到聲波和彈性波方程的全波形反演中（Gauthier et al.，1986；Mora，1987）；此外，Sheen等（Sheen et al.，2006）對彈性介質在時間域用高斯牛頓法進行了全波形反演，Choi等（Choi et al.，2008）對聲學彈性復合介質中的多分量數(shù)據(jù)用共軛梯度法進行了二維全波形反演，Epanomeritakis等（Epanomeritakis et al.，2008）結合高斯牛頓法和共軛梯度法對三維彈性波方程進行了全波形反演.將時間域的聲波方程和彈性波方程變換到頻率域，然后對給定的頻率再進行正演模擬，反演也在頻率域中求解（Pratt et al.，1988；Song，1995；Pratt et al.，1998；Prat，1999；Sirgue et al.，2004），這就是頻率域全波形反演.Pratt等（Pratt et al.，1988）首先將全波形反演思想用于頻率域中，頻率域反演和時間域反演的思想是完全一致的，只是頻率域的波場是互相解耦的，可以對每個頻率單獨進行反演，頻率域反演對初始模型也存在較大的依賴性.類似地，全波形反演還可在Laplace域中進行（Shin et al.，2008；Shin et al.，2009；Ha，2012）.

全波形反演數(shù)值上常用迭代法來求解，總體上是一個梯度類的優(yōu)化迭代求解過程，有精度高的優(yōu)點，但對初始模型有嚴重依賴性，這是梯度類反演方法的共性.Gauthier等（1986）曾證明當初始速度模型與真實模型有10%的誤差時，全波形反演就失效.在這種情況下，梯度類方法趨向找到局部極值點.為了提高全波形反演對初值的穩(wěn)健性，Bunks等（1995）提出了多尺度全波形反演方法，通過將問題分解為不同尺度求解以保證反演過程穩(wěn)定收斂.此外，人們還研究了用全局優(yōu)化搜索的方法來求解.全局優(yōu)化搜索基于隨機過程，不需目標函數(shù)的導數(shù)信息，也不考慮目標函數(shù)的形狀，不需要一個好的初值模型，而且不會陷入局部極值點，例如，Jin等（1994）提出了蒙特卡洛法反演，Smabridge等（1992）提出了遺傳算法反演，Sen和Stoffa（Stoffa et al.，1991；Sen et al.，1995）提出了模擬退火法.全局優(yōu)化搜索由于為求得目標函數(shù)的全局極值所導致的計算量極大，特別是二維和三維反演.

依據(jù)線性化反演理論，當一個初始模型在觀測數(shù)據(jù)和模擬數(shù)據(jù)之間引起大的動力學殘差時，速度模型的擾動對殘差的平方就沒有影響，從而數(shù)據(jù)對模型參數(shù)的梯度為零，但這時求得的極值點顯然不是全局極值點，這就造成速度反演的大量極值點，不但難以求得好的收斂結果或甚至發(fā)散，而且反復迭代計算也導致巨大的計算量.為此，本文研究了時間域頻率多尺度全波形反演方法，提出了基于不同尺度的頻率數(shù)據(jù)的逐級反演策略，取得了明顯效果.

2 理論方法

2.1 有限差分正演模擬

在時間域上，考慮二維聲波方程：

其中x為地面橫向坐標，z為深度縱向坐標，u（x，z，t）表示壓力，v（x，z）為介質速度，f（t）是震源函數(shù)，（xs，zs）為震源位置.若介質在激發(fā)前靜止，則初始條件為

有限差分法是數(shù)值模擬波場傳播的重要方法之一，具有計算效率高的優(yōu)點.本文用有限差分法來求解波動方程.設Δx和Δz分別為x和z方向的離散步長，Δt為時間步長，用表示lΔt時刻及空間位置為 （nΔx，mΔz）處的波場值，vn，m表示空間位置（nΔx，mΔz）處的速度離散值.用二階中心差分格式對方程（1）進行離散，得

其中，（ns，ms）為震源坐標，這里

記M＝NxNz為空間離散點數(shù)，在反演波速時也即優(yōu)化問題的未知量個數(shù).初始條件（2）的離散格式為

應用離散Fourier級數(shù)法分析易得，該二階差分格式的數(shù)值穩(wěn)定條件為

由于計算都是在有限區(qū)域上進行的，為了消除邊界反射模擬無界區(qū)域中的波場傳播，需要應用吸收邊界條件（Clayton et al.，1977；Berenger，1994；Yang et al.，2002，2003）.通常計算區(qū)域都是規(guī)則的矩形域，不妨設為Ω＝ （0，X）×（0，Z），我們采用局部化的二階邊界吸收條件（Clayton et al.，1977）：

相應的差分格式是：

其中，

顯然，邊界條件的差分格式（11）—（14）也是二階精度.

2.2 模型修正

全波形反演是一個迭代極小化觀測數(shù)據(jù)與擬合數(shù)據(jù)之間的殘量即目標函數(shù)來反演速度v的過程.殘量使用l2范數(shù)來度量，目標函數(shù)可寫成：

其中，ucal表示數(shù)值模擬的波場值，uobs表示已知的地震記錄數(shù)據(jù).模型參數(shù)v表示整個離散區(qū)域的速度，是一個M維向量.目標函數(shù)J（v）的極小化是從一個給定的初始值v0出發(fā)進行搜索，這是一個局部優(yōu)化問題.由于數(shù)據(jù)與模型參數(shù)間的非線性關系，需要多次迭代才能收斂到目標函數(shù)在v0附近的局部極小點.每次迭代都尋找一個下降方向p和搜索步長α，以確保新的迭代步

滿足

即滿足目標函數(shù)值下降，注意這里p也是一個M維向量.步長α可用一維線搜索方法（袁亞湘等，2005）得到.計算搜索方向的方法有很多，一般都要基于求解目標函數(shù)的一階導數(shù)信息，即梯度；但單純的梯度法在解附近收斂緩慢，需要用目標函數(shù)的二階導數(shù)信息，即Hessian矩陣來修正搜索方向，這類通稱為牛頓類方法.為此，將目標函數(shù)J（vk＋p）在迭代點vk處作泰勒近似展開，忽略高階項，得

對上式關于分量pl求偏導，得

實際計算中用擬牛頓法來計算，如DFP法（Fletcher et al.，1963；Davidon，1991）和BFGS法（Broyden，1970；Fletcher，1970；Goldfarb，1970；Shanno，1970），即不直接計算式（23）中的二階導數(shù)，而是在每個迭代步用某個對稱非奇異且易求逆的矩陣Bk代替該二階導數(shù)矩陣.實際計算表明，BFGS方法是擬牛頓法中最有效的一個方法，Bk和其逆均有遞推的修正公式.對傳統(tǒng)的BFGS方法，由于Hessian陣初值選取不當或者當前迭代步Hessian陣條件數(shù)很壞，會出現(xiàn)收斂慢的情況.為了克服這種情況，可在更新Hessian陣之前對矩陣進行一個尺度變換，即用τkBk代替Bk，τk為變尺度因子，這種變尺度的方法能加速BFGS法在當前迭代步的下降率（Nocedal et al.1993）.在具體計算中，為節(jié)省計算內存，常采用有限內存的BFGS方法即L－BFGS方法來計算.下面給出全波形反演的算法流程.

步驟1：給定初始模型v0，迭代中止精度ε，最大迭代次數(shù)kmax，令迭代步k＝0；

步驟2：對v0，計算目標函數(shù)值J和梯度g；設定初始搜索方向p0＝－g；令α0＝1／‖p‖2；

步驟3：如果k＜kmax，迭代循環(huán)：

1：用一維線搜索算法，得到搜索步長αk；

2：速度模型修正vk＋1＝vk＋αkpk；

3：對修正的vk＋1，計算新的目標函數(shù)值J和梯度g：

4：若梯度或目標函數(shù)值滿足終止條件，則算法終止，否則轉步驟4；

5：計算新的搜索方向pk；

6：令k＝k＋1，返回步驟3繼續(xù).

步驟4：輸出結果v.

在該算法中，需要計算目標函數(shù)的梯度g，對大規(guī)模模型的計算，我們基于下面的共軛方法來計算.

2.3 梯度計算

在BFGS方法中需要梯度的信息，這是一個M×（Nt×Nr）矩陣，這里Nr為接收點數(shù).若用差分法求梯度需要求解M次的正問題，計算量太大.差分法計算梯度適合求解規(guī)模較小的問題（Zhang et al.，2013）.共軛法即求解原問題的共軛問題，是大規(guī)模問題中最常用的梯度計算方法，它只需要兩次正問題的計算量就可以得到一個較為準確的梯度.下面給出這種方法的一般形式的推導.記

則

其中 （xr，zr）為接收點位置坐標.對J（v）關于v微分，得：

由鏈式法則，得

成立，所以：

為消除δu，令

即

算出.若L為關于u的線性算子，則顯然有

對于波動方程（1），算子L關于u是線性的，且顯而易見有

則v＋δv和u（v＋δv）同樣滿足聲波方程：

采用下列近似：

將式（38）減去式（1），再將式（39）代入，且采用弱散射近似，可以得出方程：

這里Δ是Laplace算子.考慮式（41）最后一個等號右端第1項，由Green公式可知：（42）式等號右端第2項中，已知且可令則第2項為0.等號右端第3項中，已知且可令則第3項為0.考慮式（41）最后一個等號右端第2項，由Green公式可知：

注意到（43）式右端的邊界項，在計算中希望邊界項為0，因為δu的邊界已知，所以根據(jù)情況選擇使得邊界項為0.要求上邊界為零邊界條件，其他三邊使用吸收邊界條件，顯然零邊界條件使得邊界項為0，對于吸收邊界條件，只要吸收得比較徹底，其實可以認為區(qū)域是無界的，這樣就不需要考慮邊界項的影響，可以看出好的吸收邊界對于梯度的計算是很重要的，吸收不完全會對梯度計算的準確性帶來影響.

綜合式（41）、（42）和（43），可以得到：

關于“初始”時間的限制條件：

邊界條件與u相同，上邊界為0，其他三邊為吸收邊界（8）—（10）.對于上述方程和定解條件，我們用時間逆推的方法進行求解，過程和計算量與求解正問題基本相同.最后給出梯度公式：

3 數(shù)值計算

本節(jié)用Marmousi復雜構造模型來進行反演數(shù)值計算.圖1是Marmousi速度模型，常用來做偏移成像，用于驗證各種成像方法的能力和效果（張文生，2009）.數(shù)值算例實驗設計如下：計算的實際規(guī)模是地表寬度2952m，地下深度為1488m，記錄時間為1.8s.空間步長Δx＝6m，Δz＝6m，由穩(wěn)定性條件取Δt＝0.6ms.震源的最大頻率為60Hz，中心頻率為30Hz.離散后區(qū)域的空間點數(shù)為Nx＝493和Nz＝249，時間取樣點數(shù)為Nt＝3001.實驗中共設80個炮點，每個炮點有40個接收點，雙邊接收.接收點之間的距離均為24m，炮檢距離為24m.炮點和接收點置放在同一水平線上，深度為6m.

圖1 Marmousi速度結構模型Fig.1 Marmousi velocity structure model

觀測數(shù)據(jù)是由Marmousi模型正演得到的人工數(shù)據(jù).圖2是某炮的炮記錄，圖2a無邊界吸收，圖2b有邊界吸收，比較可以看到，邊界反射得到了明顯的消除，這有助于提高梯度計算的精度.圖3是用于反演的初始速度模型圖，速度范圍是1500m·s－1至4300m·s－1，中間每層線性遞增得到的，每層速度相同；可以看到初始模型與真實模型的差異很大，已經(jīng)完全沒有原來真實模型的構造信息了.

頻率多尺度方法就是根據(jù)問題對不同尺度的頻率進行分解.為簡單起見，我們采用截去某個頻率以上的所有頻率，只保留這個頻率以下的頻率進行反演，然后增加頻率上界，依次進行反演.在本文計算中，選取的頻率多尺度實驗分為5個尺度：0～5Hz，0～10Hz，0～25Hz，0～35Hz，0～60Hz；每個尺度有逐級包含關系，這樣可以更有效地恢復模型波數(shù).為作為比較，在不采用本文頻率多尺度方法的情況下也直接對原始數(shù)據(jù)反演，圖4是迭代250步的結果，可以看到結果較差，再繼續(xù)迭代也是如此.下面我們對5個頻率尺度的波場進行逐級反演，每個尺度分別迭代50次.首先用0～5Hz頻率尺度波場進行反演，圖5a、5b和5c分別是迭代5步、20步和50步的反演結果，初始模型仍為圖3.可以明顯看到，反演結果隨迭代次數(shù)增加逐步向好的方向收斂.為了節(jié)省篇幅，下面直接給出其他尺度迭代50次的反演結果.圖6用是0～10Hz頻率尺度的波場反演迭代50次的結果，初始模型是圖5c；圖7是用0～25Hz數(shù)據(jù)反演迭代50次的結果，初始模型是圖6；圖8是用0～35Hz數(shù)據(jù)反演迭代50次的結果，初始模型是圖7；圖9是用0～60Hz數(shù)據(jù)反演迭代50次的結果，初始模型是圖7.可以看到，反演結果從圖6到圖9逐步變好；圖8和圖9的結果比較接近，不過比較可以看出在細節(jié)上還是有些差別（比較下面的圖10e和圖10f更清楚）.為了進一步與真實模型作細致的比較，我們取CDP號為180處的位置，來與真實模型作比較；圖10是5個不同頻率尺度迭代50次的全波形反演結果（紅線）與真實模型（藍線）的比較.比較可知，在模型淺部，有非常好的收斂，在深部，反演的精度雖不如淺部，但仍能反映出良好的構造信息.原因有兩方面，一是深部波場信息相對與淺層較弱，在計算中我們未對波場作任何處理；二是從數(shù)值計算的角度看，目標函數(shù)的未知量是一個向量，分量是每個網(wǎng)格點的速度值，目標函數(shù)梯度的分量對應每個網(wǎng)格點的梯度，從而也表示了不同深度的梯度，負梯度方向是未知量的修正方向，由于淺部接收信號強，淺部的變化量即梯度就會比較大，修正的時候著重于修正淺部.目標函數(shù)的二階導數(shù)反映每個分量變化對梯度的影響，對深淺有平衡作用，但二階導數(shù)是近似計算的，因此不可避免會導致深部反演結果的誤差大于淺部.圖11給出了0～50Hz頻率尺度反演的目標函數(shù)值隨迭代步的變化情況，可以看到前30次迭代，目標函數(shù)值下降很快，在80次之后，下降很小，表明繼續(xù)迭代反演結果改善也不大，其他尺度的情況大致類似，考慮到計算量，每級反演設置最大迭代次數(shù)為50.此外，我們還對加噪聲的情況也進行了數(shù)值計算，結果表明對20%的隨機噪聲也有校好的反演結果，為節(jié)省篇幅，不再詳述.

圖2 某炮的炮記錄（a）無邊界吸收；（b）有邊界吸收.Fig.2 A shot gather record with（a）no boundary absorption，and（b）boundary absorption

圖3 初始速度模型Fig.3 Initial velocity model

圖4 不用頻率多尺度方法，用0～60Hz的波場數(shù)據(jù)迭代250次的反演結果（初始模型為圖3）Fig.4 The inversion result after 250iterations with 0～60Hz wavefield data.The frequency multiscale method is not adopted（The initial model is Fig.3）

圖5 0～5Hz的波場數(shù)據(jù)的全波形反演結果（初始模型為圖3）（a）迭代5次；（b）迭代20次；（c）迭代50次.Fig.5 The full－waveform inversion results with 0～5Hz wavefield data（The initial model is Fig.3）（a）5iterations；（b）20iterations；（c）50iterations.

圖6 0～10Hz的波場數(shù)據(jù)迭代50次的全波形反演結果（初始模型為圖5c）Fig.6 The full－waveform inversion result with 0～10Hz wavefield data after 50iterations（The initial model is Fig.5c）

圖7 0～25Hz的波場數(shù)據(jù)迭代50次的全波形反演結果（初始模型為圖6）Fig.7 The full－waveform inversion result with 0～25Hz wavefield data after 50iterations（The initial model is Fig.6）

圖8 0～35Hz的波場數(shù)據(jù)迭代50次的全波形反演結果（初始模型為圖7）Fig.8 The full－waveform inversion result with 0～35Hz wavefield data after 50iterations（The initial model is Fig.7）

最后指出，全波形反演是一個典型的大規(guī)?？茖W計算問題，計算量巨大，我們采用了MPI并行算法來實現(xiàn).在計算中，對計算區(qū)域進行劃分，分別賦予不同進程計算，這樣，問題就分布存儲在各個進程中，從而滿足大規(guī)模計算對內存的要求.正問題是整個區(qū)域上關于時間遞推的計算，每個離散點的計算都有它附近點進行參與，所以在進程管理區(qū)域邊界與相鄰的進程之間有信息交換.顯然，邊界點數(shù)與子區(qū)域點數(shù)比值越小，通信時間與計算時間的比值就越小計算效率就越高.通信過多會導致并行效率變差.由BFGS算法可知，搜索方向的計算只有向量參與，沒有矩陣，向量點積只需在各進程中進行計算，再將所有進程所求的值進行相加即可，注意重疊區(qū)域不要進行多次運算.下面分析一下并行效率.我們知道，串行程序的執(zhí)行時間近似等于程序指令執(zhí)行花費的CPU時間，但并行程序相對復雜，其執(zhí)行時間等于從并行程序開始執(zhí)行，到所有進程執(zhí)行完畢，墻上時鐘走過的時間，也稱之為墻上時間（Walltime）.對各個進程，墻上時間包括計算CPU時間、通信CPU時間、同步開銷時間、同步導致的進程空閑時間.表1列出了4進程、8進程、16進程、32進程迭代反演10次的墻上計算時間的比較，由此可算出并行可擴展性.可以看到，平均的并行可擴展性大于0.9，是比較高的.計算在“科學與工程計算國家重點實驗室”三號機群系統(tǒng)上完成，該機群有282個計算刀片，每個刀片包含兩顆Intel X5550四核處理和24GB內存，單核雙精度浮點峰值性能為10.68Gflops.用32個進程進行計算，總迭代250步，計算時間約47h.

圖9 0～60Hz的波場數(shù)據(jù)迭代50次的全波形反演結果（初始模型為圖8）Fig.9 The full－waveform inversion result with 0～60Hz wavefield data after 50iterations（The initial model is Fig.8）

圖11 0～50Hz頻率多尺度反演目標函數(shù)隨迭代步k的變化情況Fig.11 The variation of objective function vale of 0～50 Hz frequency multiscale inversion with the iterative number k

表1 并行算法效率分析Table 1 Analysis of the efficiency of parallel algorithm

圖10 在CDP180位置處，不同頻率尺度波場迭代50次的全波形反演結果（紅線）與真實模型（藍線）的比較（a）初始模型；（b）0～5Hz；（c）0～10Hz；（d）0～25Hz；（e）0～35Hz；（f）0～60Hz.Fig.10 A comparison between the full－waveform inversion results after 50iterations on five different frequency scales（red line）and the exact model（blue line）at the position of CDP 180（a）Initial model；（b）0～5Hz；（c）0～10Hz；（d）0～25Hz；（e）0～35Hz；（f）0～60Hz.

4 結論

全波形反演是一個極小化模擬數(shù)據(jù)與已知記錄之間殘差的迭代過程.反演在時間域中進行，針對反演易陷入局部極小值的問題，我們發(fā)展了基于不同尺度頻率數(shù)據(jù)的逐級反演方法，即用前一尺度的反演結果作為下一尺度反演的初值模型，有效地解決了初值離真解較遠時，反演不收斂的問題.文中詳細闡述和推導了理論公式和相應的算法，并基于Marmousi復雜構造模型進行了大規(guī)模反演計算，得到了較好的反演結果，說明了該方法的有效性和對初始模型具有較高的穩(wěn)健性.波動方程全波形反演是一個典型的大規(guī)模科學計算問題，計算量巨大，文中用MPI并行方法來實現(xiàn)，提高了計算效率，這為進一步應用提供了基礎.

致謝本文的計算在科學與工程計算國家重點實驗室三號機群系統(tǒng)上完成，在此表示感謝，此外還要感謝實驗室很多老師的大力幫助和支持，在此不一一指出.最后，還要感謝審稿者提出的寶貴建議，使得本文能以目前的形式呈現(xiàn).

Berenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.J.Comput.Phys.，114（2）：185－200.Bleistein N，Cohen J K，Hagin F G.1987.Two and one－h(huán)alf dimensional Born inversion with an arbitrary reference.Geophysics，52（1）：26－36.

Broyden C G.1970.The convergence of a class of double rand minimization algorithm：2.The new algorithm.IMAJ.Appl.Math.，6（3）：222－231.

Bunks C，Saleck F，Zaleski S，et al.1995.Multiscale seismic waveform inversion.Geophysics，60（5）：1457－1473.

Bursteddle C，Ghattas O.2009.Algorithmic strategies for full waveform inversion：1Dexperiments.Geophysics，74（6）：WCC37－WCC46.

Choi Y，Min D J，Shin C.2008.Two－dimensional waveform inversion of multi－component data in acoustic－elastic coupled media.GeophysicalProspecting，56（6）：863－881.

Clayton R，Engquist B.1977.Absorbing boundary conditions for acoustic and elastic wave equations.BulletinoftheSeismological SocietyofAmerica，67（6）：1529－1540.

Dahlen F A，Hung S H，Nolet G.2000.Frechet kernels for finite frequency traveltimes—Ⅰ.Theory.Geophys.J.Int.，141（1）：157－174.

Davidon W C.1991.Variable metric method for minimization.SIAMJ.Opt.，1（1）：1－17.

Devancy A J.1984.Geophysical diffraction tomography.IEEE TransactionsonGeoscienceandRemoteSensing，22（1）：3－13.Epanomeritakis I，Akcelik V，Ghattas O，et al.2008.A Newton－CG method for large－scale three－dimensional elastic full waveform seismic inversion.InverseProblems，24（3）：1－26.

Fawcett J，Keller H B.1985.Three－dimension ray tracing and geophysical inversion in layered media.SIAMJ.Appl.Math.，45（3）：491－501.

Fletcher R，Powell M J D.1963.A rapidly convergent descent method for minimization.TheComputerJournal，6（2）：163－168.

Fletcher R.1970.A new approach to variable metric algorithms.TheComputerJournal，13（3）：317－322.

Gauthier O，Virieux J，Tarantola A.1986.Two－dimensional nonlinear inversion of seismic waveforms：Numerical results.Geophysics，51（7）：1387－1403.

Gautier S，Nolet G，Virieux J.2008.Finite－frequency tomography in a crustal environment：application to the western part of the Gulf of Corinth.Geophys.Prospect.，56（4）：493－503.

Goldfarb D.1970.A family of variable－metric methods derived by variational means.MathematicsofComputation，24（109）：23－26.Ha W，Chung W，Park E，et al.2012.2－D acoustic Laplace－domain waveform inversion of Marine field data.Geophys.J.Int.，190（1）：421－428.

Hung S H，Shen Y，Chiao L Y.2004.Imaging seismic velocity structure beneath the Iceland hot spot：a finite frequency approach.J.Geophys.Res.，109（B8）：doi：10.1029／2003JB002889.

Jin S，Madariaga R.1994.Nonlinear velocity inversion by a twostep Monte－Carlo method.Geophysics，59（4）：577－590.

Mora P.1987. Nonlinear two－dimensional elastic inversion of multioffset seismic data.Geophysics，52（9）：1211－1228.

Nocedal J，Yuan Y X.1993.Analysis of a self－scaling quasi－Newton method.MathematicalProgramming，61（1－3）：19－37.

Pratt R G，Worthington M H.1988.The application of diffraction tomography to cross－h(huán)ole seismic data.Geophysics，53（10）：1284－1294.

Pratt R G，Shin C，Hick G J.1998.Gauss－Newton and full Newton methods in frequency－space seismic waveform inversion.Geophys.J.Int.，133（2）：341－362.

Pratt R G.1999.Seismic waveform inversion in the frequency domain PartⅠ：Theory and verification in a physical scale model.Geophysics，64（3）：888－901.

Sambridge M，Drijkoningen G.1992.Genetic algorithms in seismic waveform inversion.Geophys.J.Int.，109（2）：323－342.

Sen M K，Stoffa P L.1991.Nonlinear one－dimensional seismic waveform inversion using simulated annealing.Geophysics，56（10）：1624－1638.

Sen M K，Stoffa P L.1995.Global Optimization Methods in Geophysical Inversion.The Netherlands：Elsevier Science B.V.

Shanno D F.1970.Conditioning of quasi－Newton methods for function minimization.Math.Comput.，24（111）：647－656.

Sheen D H，Tuncy K，Baag C E，et al.2006.Time domain Gauss－Newton seismic waveform inversion in elastic media.Geophys.J.Int.，167（3）：1373－1384.

Shin C，Cha Y H.2008.Waveform inversion in the Laplace domain.Geophys.J，Int.，173（3）：922－931.

Shin C，Ha W.2008.A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains.Geophysics，73（5）：VE119－VE133.

Shin C，Cha Y H.2009.Waveform inversion in the Laplace－Fourier domain.Geophys.J.Int.，177（3）：1067－1079.

Sirgue L，Pratt R G.2004.Efficient waveform inversion and imaging： A strategy for selecting temporal frequencies.Geophysics，69（1）：231－248.

Song Z M，Williamson P R，Pratt R G.1995.Frequency－domain acoustic－wave modeling and inversion of crosshole data：Part 2—Inversion method，synthetic experiments and real－data results.Geophysics，60（3）：796－809.

Stoffa P L，Sen M K.1991.Nonlinear multiparameter optimization using genetic algorithms：Inversion of plane－wave seismograms.Geophysics，56（11）：1974－1810.

Tarantola A.1984.Inversion of seismic reflection data in the acoustic approximation.Geophysics，49（8）：1259－1266.

Tarantola A.1986.A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics，51（10）：1893－1903.

Tong P，Zhao D P，Yang D H.2011.Tomography of the 1995 Kobe earthquake area：Comparison of ray and finite－frequency approaches.Geophys.J.Int.，187（1）：278－302.

Yang D H，Liu E，Zhang Z J，et al.2002.Finite－difference modeling in two－dimensional anisotropic media using a fluxcorrected transport technique.GeophysicalJournalInternational，148（2）：320－328.

Yang D H，Wang S Q，Zhang Z J，et al.2003.N－times absorbing boundary conditions for compact finite difference modeling of acoustic and elastic wave propagation in the 2－D TI medium.BulletinoftheSeismologicalSocietyofAmerica，93（6）：2389－2401.

Yang W C.1993.Nonlinear Chaotic inversion of seismic traces，Part 1：Theory and numerical experiments.ChineseJournalofGeophysics，36（2）：222－232.

Yang W C.1993.Nonlinear Chaotic inversion of seismic traces，Part 2：Lyapunov experiments and attractors.ChineseJournal ofGeophysics，36（3）：376－386.

Yuan Y X，Sun W Y.1997.Optimization Theory and Methods（in Chiese）.Beijing：Science Press.

Zhang W S.2009.Wave Equation Imaging Methods and Computations（in Chinese）.Beijing：Science Press.

Zhang W S，Luo J.2013.Full－waveform velocity inversion based on the acoustic wave equation.AmericanJournalofComputational Mathematics，3：13－20.

附中文參考文獻

袁亞湘，孫文瑜.1997.最優(yōu)化理論與方法.北京：科學出版社.

張文生.2009.波動方程成像方法及其計算.北京：科學出版社.