基于罰函數的全波形反演優化模型及其算法研究

殷雅倩 葉沐芊 張靈溪

摘 要：地震波反演技術在地質勘探中具有重要意義，該文在總結歸納拉格朗日乘子法和約化空間法的基礎上，提出了基于罰函數的求解地震波反演方法的優化模型及其算法，利用共軛梯度法求取地震波反演中最小二乘問題的最優解。在實際模型中的數值實驗表明，該模型和算法是行之有效的。

關鍵詞：全波形反演 最小二乘法 罰函數 模型

中圖分類號：P631.4 文獻標識碼：A 文章編號：1674-098X（2015）08（c）-0036-02

全波形反演技術在近年來不斷受到很多學者的關注，其研究在石油、工程等諸多勘探問題中具有重要意義。該文對全波形反演問題求解中的拉格朗日乘子法[1]和約化空間法[2]做了簡要介紹，建立了基于罰函數的反演優化模型，并且提出利用交替算法求解全波形反演的最小二乘問題的方法，最后，本文用Matlab語言編寫出算法，在實際模型中驗證了模型和算法的有效性。

20世紀60年代末，Backus等[3]提出地震波反演理論。受當時較為發達的醫學層析成像技術的影響，地質學家從中受到啟發，利用地震波對地下不同介質的反映，準確描繪地下地層結構。此后，Claerbout等[4]陸續提出了建立在波動方程基礎上的地球反演理論框架，地球反演技術得以發展起來。

全波形反演是利用完整波場信息進行地震波反演的方法[5]，其優勢在于可以精確刻畫模型細節，使得反演效果更加出色。Tarantola等[6]首先在1982年提出了基于最小二乘法進行反演的時間域全波形反演，隨后Pratt[7]又將其擴展至頻率域。頻率域全波形反演較之時間域全波形反演，擁有更快的計算速度和更高的計算精度。在求解頻率域全波形反演時，該文在研究拉格朗日乘子法和約化空間法的基礎上，建立了基于罰函數的反演優化模型，利用交替方向算法對其最小二乘問題進行求解。在數值實驗中得到體現。

1 全波形反演模型及主要方法

頻率域的全波形反演主要是基于Helmhotz波動方程的模型求解，通過對反射地震波的數據收集，進行反演的模擬，在這個過程中，要求使得模擬值與真實值之間的差值最小。也可簡要描述為如下約束最小二乘問題：

（1）

其中，為波場與觀測數據之間的轉化算子，為實驗次數，表示第次實驗，為波場，代表實際觀測數據，為地層模型，能夠反映地震波在不同介質中傳播的差距；為該公式的約束條件，是聯系波場、波源信息和地層模型的模型方程，其中，是離散之后的Helmhotz算子，是圓頻，是作用在波場上的離散的Laplacian算子，代表波源信息。

求解上述約束最小二乘問題的常用方法有兩種，一種是Lagrange乘子法，其優化目標函數更改如下：

。

（2）

其中是Lagrange乘子，表示共軛轉置。

該方法的優越性在于以Lagrange乘子的形式將約束條件增加到目標函數中，使得搜索區域擴大，減少局部最小值對優化問題的影響。同時，Lagrange乘子法能夠有效避免每次迭代都要對波動方程進行精確求解的弊端，減少計算量，提高了運算效率。但是在實際問題中，往往對較大范圍的區域進行反演，如果每次迭代更新都對的波場信息進行儲存，可能會造成存儲量太大。因此，Lagarange乘子法對于全波形反演并不完全適用。

另外一種方法則是先求解Helmhotz方程，得到，再將求得的代入目標函數得到關于的單變量函數：

。 （3）

這種方法稱為約化空間法，將代入到目標函數并對求導，得到梯度：

。 （4）

滿足共軛方程。

約化空間法與Lagrange乘子法相比，在每次更新時無需對波場信息進行存儲，即存儲量較小，但是在隨后的梯度計算時，則需要對波動方程以及共軛方程進行精確求解。由于求解這兩步的難度較大，當在處理大型問題時，計算代價就會較高。

2 基于罰函數的反演優化模型提出

由于Lagrange乘子法和約化空間法在頻率域的全波形反演上均有優勢，但也有各自的不足，本文將兩者的優勢綜合，將約束條件引入到罰函數中[8]，建立了基于罰函數的頻率域全波形反演優化模型。

罰函數的模型摒棄了伴隨波場，減少了計算量和存儲量。同時擴大了原有目標函數的搜索空間，有效的限制了出現局部最小點的情況。罰函數模型將波動方程作為懲罰項，令物理條件更加松弛。在對模型實際求解時，可以通過適當調節懲罰因子，從而平衡約束條件誤差以及目標誤差，有利于得到更為便捷有效的求解方法。引入罰函數之后，可將目標函數改寫為：

。

（5）

首先設為固定值，將上式展開，目標函數是關于的函數，對求導，最小化，對每次實驗，均有：

。

（6）

整理之后，波場可由下列公式求解得到：

。 （7）

其中。將上式改寫成矩陣形式，則有：

。 （8）

接著將設為固定值，對求導最小化，則有：

。（9）

將代入上式，整理可得：

。（10）

由于上式中的的系數矩陣是一個對角陣，我們可以直接求出。同時，令，可用牛頓法求解的迭代格式：

。 （11）

由此可以寫出相應的算法，并利用已知模型進行驗證。

算法如下：

（1）設置的初始值；

（2）進行下列交替方向迭代；

（3）根據公式（7）求出；

（4）再根據公式（11）更新；

（5）不斷迭代直到很小，退出循環。

3 數值實驗

通過數值實驗將算法應用到求解地震波反演問題中，從而驗證算法的有效性。如圖1所示，建立一個均勻速度的初始地層模型，并在其中心放置一個正方形塊體。在區域左端放置51個震源，在同層右端的相應位置放置51個接收器，多個結果疊加可以使模擬結果更加接近于真實值。將深度y方向和區域長度x方向分別劃分出51條網格線，形成51*51的網格區域，設定震源頻率為10 Hz，網格間距為20 m，網格邊界條件取吸收邊界條件。

所有計算均在一臺CPU為2.39Hz，內存為4GB的計算機上運行；編程語言為MATLAB R2012b，其中機器精度為1.1×1016。

圖2是10次迭代之后的反演圖像，已經可以比較清楚的反映地層結構。

4 結語

該文在拉格朗日乘子法和約化空間法的基礎上，提出了基于罰函數的反演優化模型，并且利用交替方向算法進行迭代；數值實驗中，該模型和算法能夠反演出較好的結果。但是如何改善算法的速度，以及如何提高算法的實用性，還有待進一步研究。

參考文獻

[1] Haber E， Ascher U M， Oldenburg D. On optimization techniques for solving nonlinear inverse problems[J].Inverse problems， 2000， 16（5）：1263-1280.

[2] Pratt R G. Frequency-domain elastic wave modeling by finite differences： A tool for crosshole seismic imaging[J].Geophysics， 2012，55（5）：626-632.

[3] Backus G E，Gilbert F.The Resolving Power of Gross Earth Data[J].Geophysical Journal of the Royal Astronomical Society，1968，16（2）：169-205.

[4] Claerbout，J.F.Toward a unified theory of reflector mapping[J].Geophysics，1971（3）：467-481.

[5] 卞愛飛，於文輝，周華偉.頻率域全波形反演方法研究進展[J].地球物理學進展，2010（3）：982-993.

[6] Tarantola A，Valette B.Generalized non-linear inverse problems solved using the least squares criterion[J]. Review of geophysics and space physics，1982，20（2）：219-232.

[7] Pratt G R， Shin C S， Hicks G J. Gaussnewton And Full Newton Methods In Frequencyspace Seismic Waveform. Inversion[J]. Geophysical Journal International，1998，133（2）：341-362.

[8] Tristan Van Leeuwen， Herrmann F J. Mitigating local minima in full-waveform inversion by expanding the search space[J]. Geophysical Journal International，2013，195（1）：661-667.