殷雅倩 葉沐芊 張靈溪

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

關(guān)鍵詞：全波形反演 最小二乘法 罰函數(shù) 模型

中圖分類號(hào)：P631.4 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1674-098X（2015）08（c）-0036-02

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

20世紀(jì)60年代末，Backus等[3]提出地震波反演理論。受當(dāng)時(shí)較為發(fā)達(dá)的醫(yī)學(xué)層析成像技術(shù)的影響，地質(zhì)學(xué)家從中受到啟發(fā)，利用地震波對(duì)地下不同介質(zhì)的反映，準(zhǔn)確描繪地下地層結(jié)構(gòu)。此后，Claerbout等[4]陸續(xù)提出了建立在波動(dòng)方程基礎(chǔ)上的地球反演理論框架，地球反演技術(shù)得以發(fā)展起來。

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

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

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

（1）

其中，為波場(chǎng)與觀測(cè)數(shù)據(jù)之間的轉(zhuǎn)化算子，為實(shí)驗(yàn)次數(shù)，表示第次實(shí)驗(yàn)，為波場(chǎng)，代表實(shí)際觀測(cè)數(shù)據(jù)，為地層模型，能夠反映地震波在不同介質(zhì)中傳播的差距；為該公式的約束條件，是聯(lián)系波場(chǎng)、波源信息和地層模型的模型方程，其中，是離散之后的Helmhotz算子，是圓頻，是作用在波場(chǎng)上的離散的Laplacian算子，代表波源信息。

求解上述約束最小二乘問題的常用方法有兩種，一種是Lagrange乘子法，其優(yōu)化目標(biāo)函數(shù)更改如下：

。

（2）

其中是Lagrange乘子，表示共軛轉(zhuǎn)置。

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

另外一種方法則是先求解Helmhotz方程，得到，再將求得的代入目標(biāo)函數(shù)得到關(guān)于的單變量函數(shù)：

。 （3）

這種方法稱為約化空間法，將代入到目標(biāo)函數(shù)并對(duì)求導(dǎo)，得到梯度：

。 （4）

滿足共軛方程。

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

2 基于罰函數(shù)的反演優(yōu)化模型提出

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

罰函數(shù)的模型摒棄了伴隨波場(chǎng)，減少了計(jì)算量和存儲(chǔ)量。同時(shí)擴(kuò)大了原有目標(biāo)函數(shù)的搜索空間，有效的限制了出現(xiàn)局部最小點(diǎn)的情況。罰函數(shù)模型將波動(dòng)方程作為懲罰項(xiàng)，令物理?xiàng)l件更加松弛。在對(duì)模型實(shí)際求解時(shí)，可以通過適當(dāng)調(diào)節(jié)懲罰因子，從而平衡約束條件誤差以及目標(biāo)誤差，有利于得到更為便捷有效的求解方法。引入罰函數(shù)之后，可將目標(biāo)函數(shù)改寫為：

。

（5）

首先設(shè)為固定值，將上式展開，目標(biāo)函數(shù)是關(guān)于的函數(shù)，對(duì)求導(dǎo)，最小化，對(duì)每次實(shí)驗(yàn)，均有：

。

（6）

整理之后，波場(chǎng)可由下列公式求解得到：

。 （7）

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

。 （8）

接著將設(shè)為固定值，對(duì)求導(dǎo)最小化，則有：

。（9）

將代入上式，整理可得：

。（10）

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

。 （11）

由此可以寫出相應(yīng)的算法，并利用已知模型進(jìn)行驗(yàn)證。

算法如下：

（1）設(shè)置的初始值；

（2）進(jìn)行下列交替方向迭代；

（3）根據(jù)公式（7）求出；

（4）再根據(jù)公式（11）更新；

（5）不斷迭代直到很小，退出循環(huán)。

3 數(shù)值實(shí)驗(yàn)

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

所有計(jì)算均在一臺(tái)CPU為2.39Hz，內(nèi)存為4GB的計(jì)算機(jī)上運(yùn)行；編程語言為MATLAB R2012b，其中機(jī)器精度為1.1×1016。

圖2是10次迭代之后的反演圖像，已經(jīng)可以比較清楚的反映地層結(jié)構(gòu)。

4 結(jié)語

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

參考文獻(xiàn)

[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ìn)展[J].地球物理學(xué)進(jìn)展，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.