Laplace-Fourier域多尺度高效全波形反演方法

胡英，張東，袁建征，黃紹堅，姚弟，徐凌，張才，秦前清

（1.中國石油勘探開發(fā)研究院物探技術(shù)研究所；2.武漢大學(xué)物理科學(xué)與技術(shù)學(xué)院；3.武漢大學(xué)測繪遙感信息工程國家重點實驗室）

Laplace-Fourier域多尺度高效全波形反演方法

胡英1，張東2，袁建征2，黃紹堅2，姚弟2，徐凌1，張才1，秦前清3

（1.中國石油勘探開發(fā)研究院物探技術(shù)研究所；2.武漢大學(xué)物理科學(xué)與技術(shù)學(xué)院；3.武漢大學(xué)測繪遙感信息工程國家重點實驗室）

針對常規(guī)Laplace-Fourier域全波形反演計算量大、耗時長的問題，提出一種多尺度高效Laplace- Fourier域全波形反演方法，并針對Marmousi和Overthrust模型進行數(shù)值模擬。Laplace-Fourier域多尺度高效全波形反演方法對于不同頻率點選擇大小不同的網(wǎng)格，并且根據(jù)Laplace衰減常數(shù)選擇不同的模型計算區(qū)域深度，既不影響反演精度，又顯著提高反演計算效率，同時由于反演網(wǎng)格數(shù)減少，反演穩(wěn)定性也有所提高。Marmousi和Overthrust模型的反演結(jié)果驗證了算法的有效性，同時在缺失低頻信息時，Laplace-Fourier域多尺度高效全波形反演結(jié)果仍較好。圖16表2參27

Laplace-Fourier域；全波形反演；多尺度網(wǎng)格；多尺度計算；反演精度；反演效率

0 引言

通過全波形反演重建地下介質(zhì)速度結(jié)構(gòu)的反演方法突破了傳統(tǒng)射線理論的局限，相比于走時反演，全波形反演結(jié)果的分辨率更高，目前已經(jīng)得到越來越廣泛的研究和關(guān)注。全波形反演通?？稍跁r間域或頻率域中進行[1-3]。由于全局優(yōu)化方法計算量太大，目前只能采用基于梯度的局部優(yōu)化方法求解多維反演問題，但這類方法容易使目標(biāo)函數(shù)陷入局部極小值[4]，從而不能得到正確的反演結(jié)果。多尺度反演方法可以減少局部極小點，降低陷入局部極小點的概率，提高反演精度[5]。從頻率域反演角度來看，相對于高頻分量，低頻信息分量具有較少的局部極小點，但由于實際觀測條件限制，通過地震記錄獲得準確的低頻分量信息存在困難[6]。目前實際應(yīng)用中都是通過其他途徑獲得較準確的初始模型用于全波形反演，從而降低陷入局部極小的風(fēng)險[7]，如通過反射層析、走時反演和Laplace域反演等方法生成初始模型[8-9]。

2008年Shin和Cha第1次提出Laplace域全波形反演方法[10]：利用衰減波形記錄的直流分量進行反演，產(chǎn)生一個對初始模型不敏感的宏觀模型，以此作為初始模型參與頻率域全波形反演。Laplace域全波形反演比頻率域全波形反演更穩(wěn)定，但其反演深度受偏移距和衰減常數(shù)影響較大。隨后Shin又提出了Laplace-Fourier域全波形反演[11]：結(jié)合Laplace域反演和頻率域反演，通過直流分量和低頻衰減波形記錄，能同時產(chǎn)生長波長和中波長數(shù)據(jù)，Laplace-Fourier域反演結(jié)果可作為頻率域全波形反演的初始模型。盡管真實地震記錄中，頻率小于5 Hz的低頻數(shù)據(jù)是不可靠的，但該方法仍然能夠產(chǎn)生一個合理的平滑模型[11]。雖然Shin等人開展了一些關(guān)于Laplace-Fourier域反演的研究[12-13]，但總體上，Laplace-Fourier域全波形反演的理論和方法研究尚不充分。

Laplace-Fourier域全波形反演需選擇不同的衰減常數(shù)與多個頻率點組合進行反演計算，若采用同時計算的并行反演方法，則需占用巨大存儲資源，限制其應(yīng)用范圍。另外也可采用所有頻率點和Laplace衰減常數(shù)點從低到高依次組合的串行計算方法[14]，解決計算存儲量大的問題，但其計算效率低，需要多層循環(huán)，運算時間長。

筆者在Shin等人研究基礎(chǔ)上提出一種多尺度高效反演方法（結(jié)合多尺度網(wǎng)格和多尺度計算區(qū)域兩種方法），既不影響反演精度，同時有效提高Laplace-Fourier域全波形反演串行計算的運行效率。

1 Laplace-Fourier域全波形反演

時間域波形u（t）的Laplace變換u?( s)可以表示成如下形式：

其中s = iω + σ[7]。

本次研究Laplace-Fourier域全波形反演在二維聲波方程基礎(chǔ)上，通過對時間域波場進行Laplace-Fourier變換實現(xiàn)。采用有限差分方法進行正演，利用完全匹配吸收層（PML）對模型邊界進行處理[15]。對二維標(biāo)量聲波方程進行Laplace-Fourier變換，得到Laplace-Fourier域二維聲波方程：

通過有限差分將（2）式離散，寫成矩陣形式：

式中A與復(fù)頻率和介質(zhì)參數(shù)、離散近似格式以及吸收邊界條件有關(guān)[16]。（3）式可用LU分解法或迭代法求解得到波形記錄[17]。假設(shè)d為Laplace-Fourier域的觀測波場記錄，u是通過初始模型正演的模擬波場記錄。使用基于L2范數(shù)（Euclidean范數(shù)）的殘差構(gòu)造目標(biāo)函數(shù)，第i個激發(fā)點、第j個接收點殘差為：

目標(biāo)函數(shù)為：

使用預(yù)處理梯度法實現(xiàn)目標(biāo)函數(shù)最小化，模型參數(shù)m的更新迭代關(guān)系為：

采用拋物線擬合的方法求取αk，?E表示目標(biāo)函數(shù)的梯度方向，Ha為近似黑塞矩陣[18]，用于提高反演的穩(wěn)定性。λ為經(jīng)驗值，本文取為5×10-4。

對實際數(shù)據(jù)進行反演時，真實的震源函數(shù)通常未知。本文采用震源估計方法，利用估計的震源代替原始震源。將（3）式中的震源項乘以一個復(fù)數(shù)比例系數(shù)e作為估計的震源，則（3）式改寫為：

反演采用的速度模型越接近真實模型，估計的震源就越接近真實震源。

2 多尺度高效反演方法原理

常規(guī)Laplace-Fourier域全波形反演迭代時的網(wǎng)格間距和模型計算區(qū)域深度都固定不變，因此計算效率低。本文提出的多尺度高效反演結(jié)合頻率域和Laplace域，通過調(diào)整網(wǎng)格間距大小和模型的計算區(qū)域深度減少反演網(wǎng)格數(shù)，提高反演效率。

2.1 Laplace域多尺度計算區(qū)域算法原理

為了研究衰減波形記錄特點，對單道數(shù)據(jù)的全波形進行分析（見圖1）。σ值足夠大時，波形記錄類似于狄拉克函數(shù)，振幅峰值出現(xiàn)在0.1 s左右，之后振幅趨于0。若σ值小，則存在后續(xù)波形記錄，σ值越小，續(xù)至波越多。續(xù)至波反映模型的深層信息，這是實現(xiàn)Laplace域多尺度計算區(qū)域反演的基礎(chǔ)。

模型最大反演深度取決于最大炮檢距和Laplace衰減常數(shù)σ值[19]。模型最大反演深度一般小于最大炮檢距的一半，若炮檢距小，反演就不能更新模型的深部區(qū)域。另外σ值對反演模型深度影響也很大。從圖1可知，當(dāng)σ值變大時，續(xù)至波基本消失，因此無法更新模型深部區(qū)域。圖2為最大炮檢距為9 km、σ分別為2 s-1和12 s-1、迭代第2次的速度梯度數(shù)值圖。可見，當(dāng)σ值較小時，反演結(jié)果既有淺層信息也有深層信息；當(dāng)σ值較大時，深層區(qū)域速度梯度為0，反演結(jié)果僅有淺層信息，反演只能修正淺層速度。

Laplace-Fourier域全波形反演需要選擇一系列σ值，按照σ值大小依次進行反演。本文按σ值由大到小依次對每個頻率點進行反演，從而實現(xiàn)模型由淺層到深層的修正，并把反演結(jié)果作為下一個頻率點反演的初始模型。常規(guī)Laplace-Fourier域全波形反演方法無論σ值多大，都計算整個速度模型，因此當(dāng)σ值較大時反演效率不高。對于Laplace-Fourier全波形反演，可根據(jù)σ值大小調(diào)整反演計算區(qū)域深度，當(dāng)σ值較大時，反演計算區(qū)域相應(yīng)減小。與固定計算區(qū)域算法相比，多尺度計算區(qū)域算法不僅可顯著提高反演效率，而且由于未知模型參數(shù)減少，反演穩(wěn)定性也相應(yīng)提高。

為了導(dǎo)出σ值與反演區(qū)域深度的關(guān)系，需估計最大炮檢距與最小炮檢距條件下的地震波傳播深度（見圖3）。當(dāng)R1接收到來自S1的初至波時，R接收到的來自S的地震波所能穿透的最大深度為OD（S和R相距很近，為便于分析，標(biāo)識距離較大），則S1與R1之間的波場最大深度范圍為OD1～OD（考慮非初至波的多次反射、折射所造成的能量損失）。

圖1 衰減的時間域波形記錄

圖2 σ值為2 s-1和12 s-1時第2次迭代的目標(biāo)函數(shù)梯度

圖3 理想地下結(jié)構(gòu)中的地震波傳播

若波形記錄振幅衰減為原來的5%，則認為此數(shù)據(jù)已衰減到足夠小，那么對于添加指數(shù)因子e-tσ的衰減波場而言，可忽略時的原波場數(shù)據(jù)[20]。對于最小炮檢距（圖3中S與R的距離，可近似為零），若將地下結(jié)構(gòu)近似為一個均勻介質(zhì)，平均縱波速度為v0，此時地震波能到達的最大深度近似表示為：

對于最大炮檢距，只考慮其中初至波所能傳播的深度，假定觀測數(shù)據(jù)的最大炮檢距為Sm，在此最大炮檢距下，地下介質(zhì)的水平方向縱波平均速度近似為v1，豎直方向縱波平均速度近似為v2，此時地震波傳播的最大深度可近似表示為：

由于實際波場能夠到達的深度范圍為h1～h2，故計算區(qū)域深度公式可表示為：

2.2 頻率域多尺度網(wǎng)格算法原理

1995年Bunks提出一種時間域反演的多尺度方法[21]：通過低通濾波將地震數(shù)據(jù)分成幾個頻率段，從低頻段到高頻段依次進行反演，每個頻率段選擇網(wǎng)格大小不同。時間域反演的多尺度方法首先利用低頻數(shù)據(jù)反演出光滑的初始模型，然后利用高頻數(shù)據(jù)反演出模型的小尺度不均勻信息，最終獲得高分辨率的反演結(jié)果[22]。低頻率點時，目標(biāo)函數(shù)的非線性度低，同時可選的網(wǎng)格間距大，網(wǎng)格點少，局部極小值點相應(yīng)減少，因此不僅可以提高反演效率，還可降低反演陷入局部極小的概率，且容易應(yīng)用于頻域反演[23]。

本文將該方法應(yīng)用于Laplace-Fourier域全波形反演，將頻率點從低到高依次進行反演，低頻率點的反演結(jié)果作為高頻率點反演的初始模型。反演中需要多次正演，正演的穩(wěn)定性直接影響反演的穩(wěn)定性。本文采用Jo等提出的優(yōu)化9點有限差分差分格式進行正演[24]。為保證正演算法穩(wěn)定性，需壓制頻散，因此每個波長至少由4個網(wǎng)格點表示[24]。本文采用正方網(wǎng)格，正演模擬的網(wǎng)格間距滿足以下條件：

由（14）式可得Δd，若vmin不變，fmax越小則λmin越大，Δd越大，對應(yīng)網(wǎng)格點減少，計算效率提高。

常規(guī)Laplace-Fourier域全波形反演方法針對所有反演頻率點采用相同的網(wǎng)格間距，因此反演效率不高。本文提出的多尺度網(wǎng)格反演策略根據(jù)頻率點調(diào)整網(wǎng)格間距：頻率低時選擇粗網(wǎng)格，頻率高時選擇細網(wǎng)格。頻率低時，粗網(wǎng)格可減輕計算負擔(dān)且更容易收斂到全局極小點。從大網(wǎng)格擴展到小網(wǎng)格時，新增加的小網(wǎng)格點若在大網(wǎng)格點連線上，則由最近的2個大網(wǎng)格點速度根據(jù)距離加權(quán)平均計算該點速度；若不在大網(wǎng)格點連線上，則由最近的4個大網(wǎng)格點速度值根據(jù)距離加權(quán)平均計算該點速度。

本文提出的Laplace-Fourier域多尺度高效全波形反演方法流程見圖4。

圖4 Laplace-Fourier多尺度高效全波形反演方法具體流程

3 模擬結(jié)果

為了驗證Laplace-Fourier域多尺度高效全波形反演方法的有效性，選用Marmousi[25]模型和Overthrust模型[26]進行模擬實驗。

3.1 Marmousi模型反演

Marmousi模型的寬度為9 212 m，深度為3 000 m，橫向和縱向網(wǎng)格間距都為12.5 m，共有737×240個網(wǎng)格（見圖5a）。反演使用的合成數(shù)據(jù)集由時間域四階交錯網(wǎng)格有限差分得到，震源為主頻15 Hz的雷克子波，時間采樣間隔為4 ms，記錄時間為3 s，炮點總計360個，均位于地表，炮點間距為25 m，每個炮點對應(yīng)360個接收點，接收點間距為25 m，最大炮檢距為9 000 m。初始模型為簡單均勻遞增模型（見圖5b）。

Laplace-Fourier域全波形反演一般選用小于5 Hz的頻率點。本文采用4個頻率點（0，2.1，2.8，4.2 Hz），σ值為2～14 s-1（間隔為2 s-1）。一般實際數(shù)據(jù)反演中，每個頻率點和Laplace系數(shù)點采用多分量同時反演的群策略，反演結(jié)果會更穩(wěn)健。由于本文討論的重點是如何提升Laplace-Fourier域全波形反演的計算效率，故不考慮該策略。表1為研究采用的多尺度高效反演參數(shù)，頻率低時網(wǎng)格間距大，網(wǎng)格數(shù)少。σ值大時模型計算區(qū)域淺，網(wǎng)格數(shù)減少，當(dāng)σ值為2、4、6、8 s-1時反演計算整個模型；σ值為10、12、14 s-1時反演深度減小一半。常規(guī)Laplace-Fourier域全波形反演方法所有頻率點網(wǎng)格間距一樣，固定為25 m，且對所有σ取值反演計算整個模型。每個Laplace-Fourier點最多迭代10次，故整個反演最多需要4 710次迭代。

圖5 Marmousi模型Laplace-Fourier域全波形反演

表1 Marmousi模型的反演策略

對比常規(guī)Laplace-Fourier域全波形反演和本文提出的多尺度高效全波形反演結(jié)果（見圖5c、5d），以及兩者單道反演結(jié)果（見圖6），發(fā)現(xiàn)兩者反演結(jié)果基本一致（都反演獲得宏觀模型），差別較小，單道反演曲線基本重合，說明兩種方法精度基本相同。常規(guī)反演耗時為10.767 h，而多尺度高效全波形反演耗時3.100 h，耗時縮短至常規(guī)反演耗時的28.8%，反演效率為常規(guī)反演的3.47倍。綜上可知，在不影響反演精度前提下多尺度高效反演計算效率比常規(guī)反演高很多。

圖6 Laplace-Fourier域多尺度高效全波形反演和常規(guī)全波形反演結(jié)果單道對比

σ值不同反演效果也不同：σ值較大，反演深度較淺，反演精度較低；σ值減小，反演深度增大，反演精度提高。圖7為頻率2.1 Hz時σ從14 s-1減小到2 s-1（間隔為2 s-1）的7個點反演結(jié)果，圖8為圖7所對應(yīng)的層間相對標(biāo)準差。從圖8可看出隨著σ值減小，有效反演深度增加，雖然有些衰減系數(shù)條件下反演結(jié)果異常，但總體上深層分辨率逐漸提高。

以多尺度高效反演結(jié)果（見圖5d）作為初始模型進行頻域全波形反演，結(jié)果見圖9。由于實際數(shù)據(jù)中低頻數(shù)據(jù)不可靠，選擇較高的頻率點進行反演，本文選擇的頻率點為5.2，7.3，9.5，12.2，14.9，17.1和20.0 Hz。從圖9c和圖10可看出，除最深層區(qū)域，Marmousi模型的每個細節(jié)均已反演出來，說明本文提出的多尺度高效反演方法產(chǎn)生的初始模型可用于后續(xù)頻域反演。

3.2 Overthrust模型反演

為了進一步驗證本文提出的多尺度高效反演方法的有效性，采用地表速度變化較復(fù)雜的Overthrust模型進行模擬實驗。

Overthrust模型的寬度為20 025 m，深度為4 675 m，網(wǎng)格間距為25 m，共有801×187個網(wǎng)格點（見圖11a）。反演使用的合成數(shù)據(jù)集由時間域四階交錯網(wǎng)格有限差分得到，震源為主頻15 Hz的雷克子波，炮點共計200個，均位于地表，炮點間距為100 m，每個炮點對應(yīng)200個接收點，最大炮檢距為20 km。初始模型是一個簡單的均勻遞增速度模型（見圖11b），速度從地表的2 600 m/s遞增到模型底部的6 000 m/s。

反演采用的頻率點和σ值與Marmousi模型模擬相同，表2為采用的多尺度高效反演參數(shù)。常規(guī)Laplace-Fourier域全波形反演對所有頻率點采用相同網(wǎng)格間距（本文固定為50 m），且對所有σ取值都反演整個模型。

圖7 不同Laplace衰減系數(shù)條件下反演結(jié)果（頻率點為2.1 Hz）

圖8 不同衰減系數(shù)層間相對標(biāo)準差

圖9 不同頻率點的頻率域反演結(jié)果

常規(guī)Laplace-Fourier域全波形反演和多尺度高效全波形反演結(jié)果（見圖11c、11d）顯示，由于Overthrust模型更加復(fù)雜，兩種方法在斷層處的反演效果均較差，其他區(qū)域反演結(jié)果較理想，且兩種方法的反演結(jié)果很接近。對比Laplace-Fourier域多尺度高效全波形反演和常規(guī)全波形反演結(jié)果（見圖12），兩者差別較小，說明兩種方法反演精度基本相同。常規(guī)反演耗時8.533 h，多尺度高效全波形反演耗時2.267 h，耗時縮短為常規(guī)反演耗時的26.6%，反演效率為常規(guī)反演方法的3.76倍。

將多尺度高效反演結(jié)果（見圖11d）作為初始模型進行頻域全波形反演，結(jié)果見圖13。選擇的頻率點分別為5.2，7.3，9.5，12.2，14.9，17.1和20.0 Hz。從圖13、14可見，除最深層區(qū)域，Overthrust模型每個細節(jié)均已反演出來，再次說明本文提出的多尺度高效反演產(chǎn)生的初始模型用于后續(xù)頻域反演的可靠性。

圖10 多尺度高效反演結(jié)果和目標(biāo)模型

圖11 Overthrust模型的Laplace-Fourier域全波形反演

表2 Overthrust模型的反演策略

3.3 低頻數(shù)據(jù)缺失Laplace-Fourier域反演

為了進一步驗證Laplace-Fourier域反演在低頻缺失情況下的反演效果，對主頻15 Hz的雷克子波進行濾波，將完全去除5 Hz以下頻率分量的波形作為震源與Marmousi模型合成波形記錄數(shù)據(jù)。合成數(shù)據(jù)集由時域四階交錯網(wǎng)格有限差分得到，時間采樣間隔為4 ms，記錄時間為3 s。炮點、接收點及初始模型設(shè)置與Marmousi模型反演相同。

圖12 Laplace-Fourier域多尺度高效全波形反演和常規(guī)全波形反演結(jié)果

反演時，震源采用正演中的子波，選擇5 Hz以下的4個頻率點（0，1.5，2.7，4.2 Hz），σ值從12 s-1降低到2 s-1（間隔為2 s-1），反演結(jié)果見圖15a。將此結(jié)果作為頻率域全波形反演的初始模型，選擇頻率點分別為7.3，9.5，12.2，14.9，17.1和20.0 Hz，反演結(jié)果見圖15b。對比缺失低頻震源的反演結(jié)果（見圖15）與未缺失低頻震源的反演結(jié)果（見圖5d、9g），以及缺失低頻震源的反演結(jié)果與目標(biāo)模型（見圖16），說明缺失低頻信息時，Laplace-Fourier域反演結(jié)果仍較好，與理論分析結(jié)果一致[27]。

圖13 7個頻率點頻率域反演最終結(jié)果

圖14 多尺度高效反演結(jié)果和目標(biāo)模型

圖15 缺失低頻震源的Marmousi模型全波形反演

圖16 缺失低頻震源的反演結(jié)果和目標(biāo)模型對比

4 結(jié)論

基于常規(guī)Laplace-Fourier域全波形反演，提出一種多尺度高效反演方法。頻率點的高低決定反演網(wǎng)格間距大小，Laplace衰減常數(shù)大小決定反演模型計算深度，本文將兩者結(jié)合應(yīng)用于Laplace-Fourier域全波形反演，頻率點低，網(wǎng)格間距大；衰減常數(shù)大，反演結(jié)果深度淺。

通過減少網(wǎng)格數(shù)和不同Laplace衰減系數(shù)下反演模型的計算深度，提高了反演效率。Marmousi和Overthrust模型實驗結(jié)果表明，多尺度高效全波形反演與常規(guī)Laplace-Fourier域全波形反演精度基本相同，而效率顯著提高。但多尺度高效全波形反演方法仍需進一步完善和改進：一方面，雖然其反演效率有很大的提升，但對于大型模型的反演耗時仍較長；另一方面，多尺度高效全波形反演與常規(guī)全波形反演結(jié)果仍有些許差異，后期需進一步研究解決。

符號注釋：

A——復(fù)阻抗矩陣；c——聲波在介質(zhì)中的傳播速度（Laplace-Fourier域）；d——觀測波場記錄；diag——取對角元素；Dk——正則化近似黑塞矩陣；e——比例系數(shù)；E——目標(biāo)函數(shù)；f——Laplace-Fourier域震源；fmax——最大頻率，Hz；F——震源向量；gk——目標(biāo)函數(shù)梯度；h1——最小炮檢距的地震波傳播深度，m；h2——最大炮檢距的初至波傳播深度，m；hm——計算區(qū)域深度，m；Hp——近似黑塞矩陣；i，j，k，l——序號；I——單位矩陣；mk——第k次迭代的模型速度參數(shù)；N——波場記錄點總數(shù)；Nr——接收點個數(shù)；Ns——炮點個數(shù)；p——Laplace-Fourier域壓強波場；P——壓強波場向量；P*——壓強波場向量的共軛向量；s——復(fù)變量，由實部和虛部組成；Sm——最大炮檢距，m；t——時間，s；u——模擬波場記錄；u（t）——時間域波形記錄；( s)——Laplace-Fourier域波形記錄；v0——介質(zhì)平均縱波速度，m/s；v1——水平方向縱波平均速度，m/s；v2——垂直方向縱波平均速度，m/s；vmin——模型最小縱波速度，m；x——水平坐標(biāo)，m；z——垂直坐標(biāo)，m；ω——角頻率，rad/s；σ——Laplace衰減常數(shù)，s-1；δr——波場殘差；δr*——波場殘差的共軛向量；αk——搜索步長，無因次；λ——正則化因子；α、β——比例系數(shù)；Δd——網(wǎng)格間距，m；λmin——最小波長，m。

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

[2] Pratt R G.Seismic waveform inversion in the frequency domain,Part 1:Theory and verification in a physical scale model[J].Geophysics,1999,64(3):888-901.

[3] 楊午陽,王西文,雍學(xué)善,等.地震全波形反演方法研究綜述[J].地球物理學(xué)進展,2013,28(2):766-776.Yang Wuyang,Wang Xiwen,Yong Xueshan,et al.The review of seismic full waveform inversion method[J].Progress in Geophysics,2013,28(2):766-776.

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

[5] 李清仁,張向君,易維啟,等.波動方程多尺度反演[J].石油地球物理勘探,2005,40(3):273-276.Li Qingren,Zhang Xiangjun,Yi Weiqi,et al.Multi-scale inversion of wave equation[J].Oil Geophysical Prospecting,2005,40(3):273-276.

[6] 賈承造,鄭民,張永峰.中國非常規(guī)油氣資源與勘探開發(fā)前景[J].石油勘探與開發(fā),2012,39(2):129-136.Jia Chengzao,Zheng Min,Zhang Yongfeng.Unconventional hydrocarbon resources in China and the prospect of exploration and development[J].Petroleum Exploration and Development,2012,39(2):129-136.

[7] Brenders A J,Pratt R G.Full waveform tomography for lithospheric imaging:Results from a blind test in a realistic crustal model[J].Geophysical Journal International,2007,168(1):133-151.

[8] Operto S,Ravaut C,Improta L,et al.Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions:A case study[J].Geophysical Prospecting,2004,52(6):625-651.

[9] 卞愛飛,於文輝,周華偉.頻率域全波形反演方法研究進展[J].地球物理學(xué)進展,2010,25(3):982-993.Bian Aifei,Yu Wenhui,Zhou Huawei.Progress in the frequency-domain full waveform inversion method[J].Progress in Geophysics,2010,25(3):982-993.

[10] Shin C,Cha Y H.Waveform inversion in the Laplace domain[J].Geophysical Journal International,2008,173(3):922-931.

[11] Shin C,Cha Y H.Waveform inversion in the Laplace-Fourier domain[J].Geophysical Journal International,2009,177(3):1067-1079.

[12] Shin C,Ha W.A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains[J].Geophysics,2008,73(5):VE119-VE133.

[13] Bae H S,Pyun S,Shin C,et al.Laplace-domain waveform inversion versus refraction traveltime tomography[J].Geophysical Journal International,2012,190(1):595-606.

[14] Shin C,Koo N H,Cha Y H,et al.Sequentially ordered single-frequency 2-D acoustic waveform inversion in the Laplace-Fourier domain[J].Geophysical Journal International,2010,181(2):935-950.

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

[16] Hustedt B,Operto S,Virieux J.Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling[J].Geophysical Journal International,2004,157(3):1269-1296.

[17] Hustedt B,Operto S,Virieux J.A multi-level direct-iterative solver for seismic wave propagation modelling:Space and wavelet approaches[J].Geophysical Journal International,2003,155(3):953-980.

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

[19] Ha W,Chung W,Park E,et al.2-D acoustic Laplace-domain waveform inversion of marine field data[J].Geophysical Journal International,2012,190(1):421-428.

[20] Agarwal K,Sylvester D,Blaauw D.A simple metric for slew rate of RC circuits based on two circuit moments[J].Computer-Aided Design of Integrated Circuits and Systems IEEE Transactions on,2004,23(9):1346-1354.

[21] Bunks C,Saleck F M,Zaleski S,et al.Multiscale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473.

[22] Boonyasiriwat C,Valasek P,Routh P.An efficient multiscale method for time-domain waveform tomography[J].Geophysics,2009,74(6):wcc59-wcc68.

[23] Ravaut C,Operto S,Improta L,et al.Multiscale imaging of complex structures from multifold wide-aperture seismic data by frequencydomain full-waveform tomography:Application to a thrust belt[J].Geophysical Journal International,2004,159(3):1032-1056.

[24] Jo C H,Shin C,Suh J H.An optimal 9-point finite-difference frequency-space 2-d scalar wave extrapolator[J].Geophysics,1996,61(2):529-537.

[25] Versteeg R.The Marmousi experience:Velocity model determination on a synthetic complex data set[J].The Leading Edge,1994,13(9):927-936.

[26] Aminzadeh F,Brac J,Kunz T.3-D salt and overthrust models:SEG/EAGE 3-D modeling series No.1,SEG[M].Tulsa:Society of Exploration Geophysicsts,1997.

[27] Ha W,Shin C.Proof of the existence of both zero-and low-frequency information in a damped wavefield[J].Journal of Applied Geophysics,2012,83:96-99.

（編輯 林敏捷）

An efficient multi-scale waveform inversion method in Laplace-Fourier domain

Hu Ying1,Zhang Dong2,Yuan Jianzheng2,Huang Shaojian2,Yao Di2,Xu Ling1,Zhang Cai1,Qin Qianqing3
(1.PetroChina Research Institute of Petroleum Exploration &Development,Beijing 100083,China;2.School of Physics and Technology,Wuhan University,Wuhan 430072,China;3.State Key Laboratory of Information Engineering in Surveying,Mapping and Remote Sensing,Wuhan University,Wuhan 430072,China)

Aiming at the problem that large computational resources and long computation time are required for the conventional Laplace-Fourier domain waveform inversion,an efficient multi-scale grid algorithm with variable computed area is proposed,and used in inversion modeling of the Marmousi and Overthrust model.This algorithm can choose a proper grid spacing automatically according to the different frequency,and adjust the depth of computing area according to the Laplace damping constant.This algorithm not only improves inversion efficiency significantly without the loss of inversion precision,but also improves the stability due to the decrease of grid number.Inversion results of the Marmousi and Overthrust model demonstrate the validity of the algorithm.In addition,the inversion results by the algorithm still can be approximate to the real model when low frequency information is missing.

Laplace-Fourier domain;full waveform inversion;multi-scale grid;multi-scale computation;inversion precision;inversion efficiency

國家科技重大專項“大型油氣田及煤層氣開發(fā)”(2011ZX05003-003)；中國石油天然氣股份有限公司配套項目“中西部前陸沖斷帶地震勘探關(guān)鍵技術(shù)應(yīng)用研究”（2011B-0406）

TE122

A

1000-0747(2015)03-0338-09

10.11698/PED.2015.03.10

胡英（1968-），女，四川綿陽人，博士，中國石油勘探開發(fā)研究院物探技術(shù)研究所高級工程師，主要從事地震資料處理技術(shù)應(yīng)用和疊前深度偏移速度建模方法研究工作。地址：北京市海淀區(qū)學(xué)院路20號，中國石油勘探開發(fā)研究院物探技術(shù)研究所，郵政編碼：100083。E-mail：hy@petrochina.com.cn

聯(lián)系作者：張東（1963-），男，湖北武漢人，博士，武漢大學(xué)物理科學(xué)與技術(shù)學(xué)院副教授，主要從事信號與信息處理理論及其在地球物理勘探等領(lǐng)域應(yīng)用研究工作。地址：湖北省武漢市武昌區(qū)八一路299號，武漢大學(xué)物理科學(xué)與技術(shù)學(xué)院，郵政編碼：430072。E-mail：dongz@whu.edu.cn

2014-07-08

2015-04-02