胡英，張東，袁建征，黃紹堅(jiān)，姚弟，徐凌，張才，秦前清

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

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

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

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

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

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

0 引言

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

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

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

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

1 Laplace-Fourier域全波形反演

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

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

本次研究Laplace-Fourier域全波形反演在二維聲波方程基礎(chǔ)上，通過對時(shí)間域波場進(jìn)行Laplace-Fourier變換實(shí)現(xiàn)。采用有限差分方法進(jìn)行正演，利用完全匹配吸收層（PML）對模型邊界進(jìn)行處理[15]。對二維標(biāo)量聲波方程進(jìn)行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個(gè)激發(fā)點(diǎn)、第j個(gè)接收點(diǎn)殘差為：

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

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

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

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

反演采用的速度模型越接近真實(shí)模型，估計(jì)的震源就越接近真實(shí)震源。

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

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

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

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

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

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

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

圖1 衰減的時(shí)間域波形記錄

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

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

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

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

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

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

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

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

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

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

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

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

3 模擬結(jié)果

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

3.1 Marmousi模型反演

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

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

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

表1 Marmousi模型的反演策略

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

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

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

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

3.2 Overthrust模型反演

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

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

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

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

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

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

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

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

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

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

表2 Overthrust模型的反演策略

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

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

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

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

圖13 7個(gè)頻率點(diǎn)頻率域反演最終結(jié)果

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

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

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

4 結(jié)論

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

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

符號注釋：

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——計(jì)算區(qū)域深度，m；Hp——近似黑塞矩陣；i，j，k，l——序號；I——單位矩陣；mk——第k次迭代的模型速度參數(shù)；N——波場記錄點(diǎn)總數(shù)；Nr——接收點(diǎn)個(gè)數(shù)；Ns——炮點(diǎn)個(gè)數(shù)；p——Laplace-Fourier域壓強(qiáng)波場；P——壓強(qiáng)波場向量；P*——壓強(qiáng)波場向量的共軛向量；s——復(fù)變量，由實(shí)部和虛部組成；Sm——最大炮檢距，m；t——時(shí)間，s；u——模擬波場記錄；u（t）——時(shí)間域波形記錄；( 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é)進(jìn)展,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ìn)展[J].地球物理學(xué)進(jìn)展,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

國家科技重大專項(xiàng)“大型油氣田及煤層氣開發(fā)”(2011ZX05003-003)；中國石油天然氣股份有限公司配套項(xiàng)目“中西部前陸沖斷帶地震勘探關(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