殷雅倩 葉沐芊 張靈溪
摘 要:地震波反演技術在地質勘探中具有重要意義,該文在總結歸納拉格朗日乘子法和約化空間法的基礎上,提出了基于罰函數(shù)的求解地震波反演方法的優(yōu)化模型及其算法,利用共軛梯度法求取地震波反演中最小二乘問題的最優(yōu)解。在實際模型中的數(shù)值實驗表明,該模型和算法是行之有效的。
關鍵詞:全波形反演 最小二乘法 罰函數(shù) 模型
中圖分類號:P631.4 文獻標識碼:A 文章編號:1674-098X(2015)08(c)-0036-02
全波形反演技術在近年來不斷受到很多學者的關注,其研究在石油、工程等諸多勘探問題中具有重要意義。該文對全波形反演問題求解中的拉格朗日乘子法[1]和約化空間法[2]做了簡要介紹,建立了基于罰函數(shù)的反演優(yōu)化模型,并且提出利用交替算法求解全波形反演的最小二乘問題的方法,最后,本文用Matlab語言編寫出算法,在實際模型中驗證了模型和算法的有效性。
20世紀60年代末,Backus等[3]提出地震波反演理論。受當時較為發(fā)達的醫(yī)學層析成像技術的影響,地質學家從中受到啟發(fā),利用地震波對地下不同介質的反映,準確描繪地下地層結構。此后,Claerbout等[4]陸續(xù)提出了建立在波動方程基礎上的地球反演理論框架,地球反演技術得以發(fā)展起來。
全波形反演是利用完整波場信息進行地震波反演的方法[5],其優(yōu)勢在于可以精確刻畫模型細節(jié),使得反演效果更加出色。Tarantola等[6]首先在1982年提出了基于最小二乘法進行反演的時間域全波形反演,隨后Pratt[7]又將其擴展至頻率域。頻率域全波形反演較之時間域全波形反演,擁有更快的計算速度和更高的計算精度。在求解頻率域全波形反演時,該文在研究拉格朗日乘子法和約化空間法的基礎上,建立了基于罰函數(shù)的反演優(yōu)化模型,利用交替方向算法對其最小二乘問題進行求解。在數(shù)值實驗中得到體現(xiàn)。
1 全波形反演模型及主要方法
頻率域的全波形反演主要是基于Helmhotz波動方程的模型求解,通過對反射地震波的數(shù)據(jù)收集,進行反演的模擬,在這個過程中,要求使得模擬值與真實值之間的差值最小。也可簡要描述為如下約束最小二乘問題:
(1)
其中,為波場與觀測數(shù)據(jù)之間的轉化算子,為實驗次數(shù),表示第次實驗,為波場,代表實際觀測數(shù)據(jù),為地層模型,能夠反映地震波在不同介質中傳播的差距;為該公式的約束條件,是聯(lián)系波場、波源信息和地層模型的模型方程,其中,是離散之后的Helmhotz算子,是圓頻,是作用在波場上的離散的Laplacian算子,代表波源信息。
求解上述約束最小二乘問題的常用方法有兩種,一種是Lagrange乘子法,其優(yōu)化目標函數(shù)更改如下:
。
(2)
其中是Lagrange乘子,表示共軛轉置。
該方法的優(yōu)越性在于以Lagrange乘子的形式將約束條件增加到目標函數(shù)中,使得搜索區(qū)域擴大,減少局部最小值對優(yōu)化問題的影響。同時,Lagrange乘子法能夠有效避免每次迭代都要對波動方程進行精確求解的弊端,減少計算量,提高了運算效率。但是在實際問題中,往往對較大范圍的區(qū)域進行反演,如果每次迭代更新都對的波場信息進行儲存,可能會造成存儲量太大。因此,Lagarange乘子法對于全波形反演并不完全適用。
另外一種方法則是先求解Helmhotz方程,得到,再將求得的代入目標函數(shù)得到關于的單變量函數(shù):
。 (3)
這種方法稱為約化空間法,將代入到目標函數(shù)并對求導,得到梯度:
。 (4)
滿足共軛方程。
約化空間法與Lagrange乘子法相比,在每次更新時無需對波場信息進行存儲,即存儲量較小,但是在隨后的梯度計算時,則需要對波動方程以及共軛方程進行精確求解。由于求解這兩步的難度較大,當在處理大型問題時,計算代價就會較高。
2 基于罰函數(shù)的反演優(yōu)化模型提出
由于Lagrange乘子法和約化空間法在頻率域的全波形反演上均有優(yōu)勢,但也有各自的不足,本文將兩者的優(yōu)勢綜合,將約束條件引入到罰函數(shù)中[8],建立了基于罰函數(shù)的頻率域全波形反演優(yōu)化模型。
罰函數(shù)的模型摒棄了伴隨波場,減少了計算量和存儲量。同時擴大了原有目標函數(shù)的搜索空間,有效的限制了出現(xiàn)局部最小點的情況。罰函數(shù)模型將波動方程作為懲罰項,令物理條件更加松弛。在對模型實際求解時,可以通過適當調節(jié)懲罰因子,從而平衡約束條件誤差以及目標誤差,有利于得到更為便捷有效的求解方法。引入罰函數(shù)之后,可將目標函數(shù)改寫為:
。
(5)
首先設為固定值,將上式展開,目標函數(shù)是關于的函數(shù),對求導,最小化,對每次實驗,均有:
。
(6)
整理之后,波場可由下列公式求解得到:
。 (7)
其中。將上式改寫成矩陣形式,則有:
。 (8)
接著將設為固定值,對求導最小化,則有:
。(9)
將代入上式,整理可得:
。(10)
由于上式中的的系數(shù)矩陣是一個對角陣,我們可以直接求出。同時,令,可用牛頓法求解的迭代格式:
。 (11)
由此可以寫出相應的算法,并利用已知模型進行驗證。
算法如下:
(1)設置的初始值;
(2)進行下列交替方向迭代;
(3)根據(jù)公式(7)求出;
(4)再根據(jù)公式(11)更新;
(5)不斷迭代直到很小,退出循環(huán)。
3 數(shù)值實驗
通過數(shù)值實驗將算法應用到求解地震波反演問題中,從而驗證算法的有效性。如圖1所示,建立一個均勻速度的初始地層模型,并在其中心放置一個正方形塊體。在區(qū)域左端放置51個震源,在同層右端的相應位置放置51個接收器,多個結果疊加可以使模擬結果更加接近于真實值。將深度y方向和區(qū)域長度x方向分別劃分出51條網(wǎng)格線,形成51*51的網(wǎng)格區(qū)域,設定震源頻率為10 Hz,網(wǎng)格間距為20 m,網(wǎng)格邊界條件取吸收邊界條件。
所有計算均在一臺CPU為2.39Hz,內存為4GB的計算機上運行;編程語言為MATLAB R2012b,其中機器精度為1.1×1016。
圖2是10次迭代之后的反演圖像,已經(jīng)可以比較清楚的反映地層結構。
4 結語
該文在拉格朗日乘子法和約化空間法的基礎上,提出了基于罰函數(shù)的反演優(yōu)化模型,并且利用交替方向算法進行迭代;數(shù)值實驗中,該模型和算法能夠反演出較好的結果。但是如何改善算法的速度,以及如何提高算法的實用性,還有待進一步研究。
參考文獻
[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.