頻率多尺度全波形速度反演

張文生，羅嘉，滕吉文

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

2中國科學(xué)院地質(zhì)與地球物理研究所，北京100029

1 引言

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

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

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

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

2 理論方法

2.1 有限差分正演模擬

在時(shí)間域上，考慮二維聲波方程：

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

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

其中，（ns，ms）為震源坐標(biāo)，這里

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

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

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

相應(yīng)的差分格式是：

其中，

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

2.2 模型修正

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

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

滿足

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

對上式關(guān)于分量pl求偏導(dǎo)，得

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

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

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

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

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

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

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

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

5：計(jì)算新的搜索方向pk；

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

步驟4：輸出結(jié)果v.

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

2.3 梯度計(jì)算

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

則

其中 （xr，zr）為接收點(diǎn)位置坐標(biāo).對J（v）關(guān)于v微分，得：

由鏈?zhǔn)椒▌t，得

成立，所以：

為消除δu，令

即

算出.若L為關(guān)于u的線性算子，則顯然有

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

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

采用下列近似：

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

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

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

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

關(guān)于“初始”時(shí)間的限制條件：

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

3 數(shù)值計(jì)算

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

圖1 Marmousi速度結(jié)構(gòu)模型Fig.1 Marmousi velocity structure model

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

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

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

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

圖11 0～50Hz頻率多尺度反演目標(biāo)函數(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次的全波形反演結(jié)果（紅線）與真實(shí)模型（藍(lán)線）的比較（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 結(jié)論

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

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

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

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

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