地震波形反演的稀疏約束正則化方法

王 薇，韓 波，唐錦萍

1 哈爾濱工業(yè)大學(xué)數(shù)學(xué)系，哈爾濱 150001

2 復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院，上海 200433

1 引 言

采用彈性動(dòng)力學(xué)中的波動(dòng)方程來(lái)描述地震波在地下介質(zhì)中的傳播過(guò)程，具有能夠利用完全的波場(chǎng)信息的優(yōu)勢(shì)，因此一直都是地震勘探中的研究熱點(diǎn).如在定位地下分布的油氣資源時(shí)，一般的做法是在地面上人工激發(fā)地震波.由于介質(zhì)的非均勻性，當(dāng)?shù)卣鸩ㄔ诘貙咏橘|(zhì)中向各方向傳播時(shí)會(huì)產(chǎn)生反射、衍射、散射和透射等現(xiàn)象，部分地震波返回到地面就構(gòu)成了所接收的地震觀測(cè)數(shù)據(jù).地震勘探的主要工作就是基于波動(dòng)方程，從這些觀測(cè)數(shù)據(jù)中，反演出地層剖面及介質(zhì)的物性參數(shù)，進(jìn)而確定地下構(gòu)造，以此作為油氣勘探或工程物探的基礎(chǔ).從位移響應(yīng)的理論合成數(shù)據(jù)應(yīng)與實(shí)際測(cè)量數(shù)據(jù)相吻合的觀點(diǎn)出發(fā)，應(yīng)用中把波動(dòng)方程速度反演問(wèn)題轉(zhuǎn)化為非線性函數(shù)的極小值問(wèn)題.

近些年，從線性［1－2］到非線性［3－4］，從時(shí)域［5］到頻域［6－7］及 Laplace域［8］，從疊后反演到疊前反演［9－10］，波形反演獲得了廣泛研究.在實(shí)際應(yīng)用中，相比較線性地震波反演，非線性地震波形反演更接近實(shí)際情況，并且Mora也證明了在正常的勘探條件下，只有采用完全非線性地震波形反演方法才能觀測(cè)到地震速度波場(chǎng)的所有波長(zhǎng)分量.Tarantola［3］、Mora［4］等對(duì)完全非線性波形反演進(jìn)行了詳細(xì)研究，力求通過(guò)迭代下降方法，尋求使得目標(biāo)函數(shù)最小的速度模型.隨后，各類優(yōu)化算法，如最速下降法，牛頓法，共軛梯度法等得以應(yīng)用.由于波形反演的非線性性，使得目標(biāo)函數(shù)存在大量的局部極小值，反演結(jié)果對(duì)初值依賴較強(qiáng).為了克服局部極值問(wèn)題，各類方法被討論，其中包括模擬退火法［11］與遺傳算法［12］等統(tǒng)計(jì)學(xué)方法，將多尺度思想、同倫思想與各類優(yōu)化算法相結(jié)合的多種方法［13－15］，以及小波變換的多尺度方法［16］等.

由于波形反演是不適定的，即數(shù)值結(jié)果對(duì)數(shù)據(jù)比較敏感，而觀測(cè)數(shù)據(jù)不可避免地存在噪聲，所以必須采用正則化方法求其近似解，解決這一問(wèn)題首選方法是Tikhonov正則化方法.傳統(tǒng)的Tikhonov正則化方法選取二次罰項(xiàng)，其作用在于減弱原不適定問(wèn)題近似解的震蕩性，使得近似解具有一定的光滑性，從而給出穩(wěn)定的近似解.但同時(shí)也會(huì)導(dǎo)致解的過(guò)度光滑，而偏離實(shí)際.然而，在實(shí)際應(yīng)用中，常常會(huì)遇到解為不連續(xù)函數(shù)或含尖點(diǎn)的函數(shù)，此時(shí)經(jīng)典的正則化方法因其過(guò)度光滑的特點(diǎn)已不再適用，而TV正則化［17］（Total variation regularization）與稀疏正則化方法（Regularization with sparsity constraints）（又稱lp約束正則化方法（1≤p＜2））則能很好地反演出跳躍性較大的參數(shù)部分.為了更好地刻畫(huà)參數(shù)的性態(tài)，本文引入稀疏約束正則化方法.

2 稀疏約束正則化方法

稀疏約束正則化方法可以看作是選取罰項(xiàng)為加權(quán)的lp范數(shù)，適用于求解反演參數(shù)且具有稀疏性的反問(wèn)題.所謂的“稀疏性”是指解序列大部分為零（近似為零）或者解在正交基或框架下具有稀疏表示，即大部分系數(shù)為零（近似為零）.稀疏性成為刻畫(huà)解的一種方式，特別是在壓縮感知中，應(yīng)用l1約束正則化方法，可以有效地減少采樣數(shù)據(jù)，節(jié)省存儲(chǔ)空間，并且通過(guò)少量的信號(hào)實(shí)現(xiàn)信號(hào)的準(zhǔn)確或近似重構(gòu)，因此稀疏正則化方法在信號(hào)處理及圖像去噪中獲得了廣泛應(yīng)用［18－19］，同時(shí)在地震層析反演［20－21］，地 震波阻抗反演［22］中體現(xiàn)了明顯的優(yōu)勢(shì).

2004年Daubechies等［23］給出了稀疏正則化方法與迭代伸縮算法的理論分析，指出其用于求解反演參數(shù)且具有稀疏性的線性反問(wèn)題的有效性.隨后，Ramlau等［24］將稀疏約束正則化方法推廣到非線性問(wèn)題.對(duì)非線性不適定算子方程F（x）＝y(tǒng)，F(xiàn)：X →Y，若X，Y為序列空間 （X＝lp，Y＝l2），相應(yīng)的稀疏約束泛函可表示為

當(dāng)p＝2時(shí)即為經(jīng)典的Tikhonov正則化方法.若X為Hilbert空間，適當(dāng)?shù)匾肟蚣芾碚摚我鈞∈X存在序列 ｛φi｝使x＝，稀疏約束泛函可寫(xiě)為

其中α＞0為正則參數(shù).同時(shí)在數(shù)值算例中指出，隨著p的減少，近似解呈現(xiàn)更加稀疏的性態(tài)；但從優(yōu)化的角度考慮，當(dāng)0＜p≤1時(shí)，罰項(xiàng)部分失去了凸性與零點(diǎn)的可微性，因此本文的討論限定在1＜p＜2上.易知，泛函（1）的極小點(diǎn)滿足

這里p，q＞1為滿足1／p＋1／q＝1的對(duì)偶對(duì)，μk為尺度參數(shù).上式可看作是對(duì)普通的最速下降法（p＝2）（又稱之為L(zhǎng)andweber方法）的推廣.文獻(xiàn)［25］對(duì)這一方法的收斂性進(jìn)行了詳細(xì)分析.本文將應(yīng)用對(duì)偶方法（4）求解稀疏約束泛函（1）的極小點(diǎn)，類似地，為求解稀疏化后的約束泛函（2），迭代公式（3）可寫(xiě)為

3 波形反演模型

本文基于二維聲波方程波形反演問(wèn)題探討稀疏約束正則化方法的可行性和有效性.描述地震波傳播的二維聲波方程的模型是

其中，x，z分別是水平方向和垂直方向，z＝0為地表.u（x，z，t）為質(zhì)點(diǎn)的位移函數(shù).s（x，z，t）為震源函數(shù)，并且s（x，z，t）＝0，t＜0.v（x，z）為介質(zhì)在（x，z）處的速度.

波動(dòng)方程的邊界條件為

初始條件為

這樣，（5）—（8）就構(gòu)成了一個(gè)聲波方程的正演模型.如果加上地表的測(cè)量數(shù)據(jù)

就構(gòu)成了確定速度v的聲波方程反演問(wèn)題.首先對(duì)（5）式進(jìn)行數(shù)值差分離散，即有

相應(yīng)的邊界條件（6），（7）式與初始條件（8）式可寫(xiě)為

其中nx，nz分別為矩形剖分在x，z方向的網(wǎng)格個(gè)數(shù)，nt為時(shí)間間隔個(gè)數(shù)，Δx，Δz，Δt分別為空間、時(shí)間網(wǎng)格上的剖分步長(zhǎng)，記v＝v（i，j），為待反演的速度模型，y＝u（i，1，k）（i＝1，2，…，nx－1，k＝2，3，…，nt）為觀測(cè)數(shù)據(jù)，用F表示模型空間到數(shù)據(jù)空間的映射，即離散正演算子，則通常的地震勘探反演問(wèn)題可描述為有限維非線性算子方程：

具體可參見(jiàn)文獻(xiàn)［14，17］.

4 數(shù)值模擬

本節(jié)將考慮孔洞和分層兩類不同模型，應(yīng)用對(duì)偶方法進(jìn)行數(shù)值模擬，以驗(yàn)證稀疏約束正則化方法的有效性.算例1與算例2記錄了當(dāng)模型本身為稀疏時(shí)稀疏約束正則化方法的反演結(jié)果；當(dāng)模型為不稀疏時(shí)，先應(yīng)用小波變換對(duì)其進(jìn)行稀疏化處理，再應(yīng)用稀疏約束正則化方法反演，如算例3、4與5.

在數(shù)值模擬中，選取空間范圍為800m×800m，空間間隔Δx＝Δz＝40m，時(shí)間采樣間隔為Δt＝0.004s，震源子波選取最常用的雷克（Ricker）子波.其數(shù)學(xué)表達(dá)式為

其中，f＝60Hz為子波主頻.

為了測(cè)試稀疏正則化方法的抗噪能力，分別對(duì)模型1—3添加1%的隨機(jī)噪聲，對(duì)模型4與模型5添加25%的隨機(jī)噪聲.在求解稀疏約束的優(yōu)化問(wèn)題中，應(yīng)用對(duì)偶方法（4）式，選取背景值v0＝2800作為反演的初值，正則參數(shù)為

其中，ξ∈（0，1），α0＝10－3.選取迭代的停止準(zhǔn)則為廣義偏差原則，即

其中，k＊為迭代步數(shù)，τ＞1，具體值可依靠經(jīng)驗(yàn)選取.本文在算例1、2中選取τ＝1.3，算例3～5中選取τ＝3.

為了比較算法精度，定義相對(duì)誤差為這里v為反演結(jié)果，v＊為真實(shí)值，其中范數(shù)皆為l2范數(shù).

4.1 孔洞模型反演結(jié)果

算例1：真實(shí)模型如圖1a所示，在對(duì)偶方法中，選取不同的p值（2，1.6，1.2），得到的反演結(jié)果分別如圖1（b，c，d）所示.結(jié)果表明，隨著p的減少，稀疏約束正則化方法可以獲得更稀疏的解，這與稀疏約束正則化方法的預(yù)期是一致的.表1記錄了不同p相應(yīng)的迭代步數(shù)k＊，CPU運(yùn)行時(shí)間T（s）及近似解與真解的相對(duì)誤差err的計(jì)算結(jié)果；具體的相對(duì)誤差（err）及偏差 （disp＝ ‖F(xiàn)（vδk）－yδ‖）的變化情況可見(jiàn)圖2（a，b），不難看到隨著p的減少，在迭代誤差與計(jì)算效率上，稀疏約束正則化方法都呈現(xiàn)一定的優(yōu)勢(shì).

表1 算例1的數(shù)值結(jié)果Table 1 The numerical results of model 1

算例2：考慮更復(fù)雜的介質(zhì)，真實(shí)模型如圖3a所示，在對(duì)偶方法中，選取不同的p值（2，1.2，1.1），得到的結(jié)果如圖3（b，c，d）所示.結(jié)果表明，對(duì)比經(jīng)典方法，稀疏約束正則化方法反演效果要好的多.

算例3：選取真實(shí)模型如圖4a所示，應(yīng)用尺度為4的Haar小波對(duì)其進(jìn)行稀疏化，再利用對(duì)偶方法，選取不同的p值（2，1.2，1.1）進(jìn)行求解，結(jié)果如圖4（b，c，d）所示.結(jié)果表明，對(duì)比經(jīng)典方法，稀疏約束正則化方法反演效果要清晰得多.但是可以注意到，當(dāng)p減少到一定程度時(shí)，反演結(jié)果會(huì)過(guò)度稀疏化，即出現(xiàn)能量集中的現(xiàn)象，導(dǎo)致某些網(wǎng)格點(diǎn)反演值偏高.因此導(dǎo)致了p＝1.1時(shí)的計(jì)算結(jié)果反不如p＝1.2好.為此，記錄選取不同的p值相應(yīng)的小波系數(shù)重構(gòu)情況，如圖5所示，其中橫坐標(biāo)代表小波系數(shù)指標(biāo)i，縱坐標(biāo)表示小波系數(shù)值.首先可以看到，應(yīng)用小波分解獲得的小波系數(shù)的確是稀疏的；其次對(duì)比真實(shí)模型的小波系數(shù)值（實(shí)線）與近似解的小波系數(shù)值（星號(hào)），隨著p的減少，重構(gòu)系數(shù)確實(shí)呈現(xiàn)更稀疏的性態(tài)；特別注意的是的當(dāng)p＝1.1時(shí)，重構(gòu)的小波系數(shù)確實(shí)更稀疏了，但是真實(shí)模型的小波系數(shù)卻沒(méi)那么稀疏，即稱為過(guò)度稀疏的現(xiàn)象，因此參數(shù)p需適當(dāng)選取；從另一角度，也可以認(rèn)為經(jīng)小波變換后的小波系數(shù)不夠稀疏，因此為了獲得更稀疏系數(shù)，如何選擇適當(dāng)?shù)目蚣芑蚧?，也是我們進(jìn)一步考慮的問(wèn)題.

4.2 分層介質(zhì)反演結(jié)果

考慮分層模型，真實(shí)模型如圖6a，7a所示，此時(shí)模型不稀疏，為了稀疏化本文應(yīng)用小波變換對(duì)其進(jìn)行稀疏化處理，選取尺度為4的Haar小波.

圖3 算例2（a）真實(shí)模型速度；（b）Landweber方法（p＝2）的反演結(jié)果；（c）p＝1.2對(duì)偶方法的反演結(jié)果；（d）p＝1.1的對(duì)偶方法的反演結(jié)果.Fig.3 Model 2（a）Velocity of the real model；（b）Inversion result of Landweber method（p＝2）；（c）Inversion result of dual method with p＝1.2；（d）Inversion result of dual method with p＝1.1.

圖4 算例3（a）真實(shí)模型速度；（b）Landweber方法（p＝2）的反演結(jié)果；（c）p＝1.2對(duì)偶方法的反演結(jié)果；（d）p＝1.1的對(duì)偶方法的反演結(jié)果.Fig.4 Model 3（a）Velocity of the real model；（b）Inversion result of Landweber method（p＝2）；（c）Inversion result of dual method with p＝1.2；（d）Inversion result of dual method with p＝1.1.

算例4：在對(duì)偶方法中，選取不同的p值（2，1.4，1.2），結(jié)果如圖6（b，c，d）所示.

算例5：在對(duì)偶方法中，選取不同的p值（2，1.6，1.2），結(jié)果如圖7（b，c，d）所示.

綜合圖4與圖5的反演結(jié)果，可見(jiàn)對(duì)比經(jīng)典Landweber方法，隨著p的減少，對(duì)偶方法對(duì)不連續(xù)介質(zhì)模型的邊緣具有良好的識(shí)別能力.

5 結(jié) 語(yǔ)

完全的地震波形反演問(wèn)題本質(zhì)上是非線性的，在反演過(guò)程中經(jīng)常會(huì)遇到介質(zhì)不規(guī)則或不連續(xù)的情況，而應(yīng)用傳統(tǒng)的Tikhonov正則化方法又會(huì)出現(xiàn)解過(guò)度光滑現(xiàn)象，從而使反演結(jié)果失效.為此，本文引入了非線性稀疏約束正則化方法，將其應(yīng)用于二維聲波方程波形反演問(wèn)題.數(shù)值結(jié)果表明，在稀疏約束正則化方法中，隨著p的減少，反演結(jié)果會(huì)呈現(xiàn)更稀疏的性態(tài)，且對(duì)非連續(xù)介質(zhì)模型的邊緣具有更好的識(shí)別能力.這項(xiàng)試驗(yàn)無(wú)疑是有實(shí)際意義的，下一步我們將嘗試將這一理論應(yīng)用到更復(fù)雜的地震波形反演問(wèn)題中.同時(shí)，從稀疏約束泛函來(lái)看，當(dāng)p逐漸減少時(shí)，罰項(xiàng)部分凸性變?nèi)?，我們有理由懷疑此時(shí)對(duì)偶方法的收斂域變小，因此如何克服局部極小值問(wèn)題，仍需進(jìn)一步探討.

圖5 算例3的小波系數(shù)結(jié)果（a）p＝2的反演系數(shù)結(jié)果；（b）p＝1.3的反演系數(shù)結(jié)果；（c）p＝1.2的反演系數(shù)結(jié)果；（d）p＝1.1的反演系數(shù)結(jié)果.Fig.5 The wavelet coefficient of Model 3（a）The wavelet coefficient of inversion with p＝2；（b）The wavelet coefficient of inversion with p＝1.3；（c）The wavelet coefficient of inversion with p＝1.2；（d）The wavelet coefficient of inversion with p＝1.1.

（References）

［1］Tarantola A.Linearized inversion of seismic reflection data.Geophysical Prosp.，1984，32（6）：998－1015.

［2］Berkhout A J. Multidimensional linearized inversion and seismic migration.Geophysics，1984，49（11）：1881－1895.

［3］Tarantola A.A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics，1986，51（10）：1983－1903.

［4］Mora P. Nonlinear two－dimensional elastic inversion of multioffset seismic data.Geophysics，1987，52（9）：1211－1228.

［5］Zhou C X，Cai W Y，Luo Y，et al.Acoustic wave equation traveltime and waveform inversion of cross hole seismic data.Geophysics，1995，60（3）：765－773.

［6］Pratt R G，Shin C，Hicks G J.Gauss－Newton and full Newton methods in frequency－space seismic waveform inversion.Geophysical Journal International，1998，133（2）：341－362.

［7］Shin C，Min D J.Waveform inversion using a logarithmic wavefield.Geophysics，2006，71（3）：R31－R42.

［8］Shin C，Cha Y H.Waveform inversion in the Laplace domain.Geophysical Journal International，2008，173（3）：922－931.

［9］Xu T，McMechan G A，Sun R.3－D prestack full－wavefield inversion.Geophysics，1995，60（6）：1805－1818.

［10］苑書(shū)金.疊前地震反演技術(shù)的進(jìn)展及其在巖性油氣藏勘探中的應(yīng)用.地球物理學(xué)進(jìn)展，2007，22（3）：879－886.Yuan S J.Progress of pre－stack inversion and application in exploration of the lithological reservoirs.Progress in Geophysics（in Chinese），2007，22（3）：879－886.

［11］Sen M K，Stoffa P L.Nonlinear one－dimensional seismic waveform inversion using simulated annealing.Geophysics，1991，56（10）：1624－1638.

［12］Stoffa P L，Sen M K.Nonlinear multi－parameter optimization using genetic algorithms inversion of plane－wave seismograms.Geophysics，1991，56（11）：1794－1810.

［13］Bunks C，Saleck F M，Zaleski S，et al.Multiscale seismic waveform inversion.Geophysics，1995，60（5）：1457－1473.

［14］傅紅筍，韓波.二維波動(dòng)方程速度的正則化－同倫－測(cè)井約束反演.地球物理學(xué)報(bào)，2005，48（6）：1441－1448.Fu H S，Han B.A regularization homotopy method for the inverse problem of 2－D wave equation and well log constraint inversion.Chinese J.Geophys.（in Chinese），2005，48（6）：1441－1448.

［15］Han B，F(xiàn)u H S，Liu H.A homotopy method for well－log constraint waveform inversion.Geophysics，2007，72（2）：R1－R7.

［16］孟鴻鷹，劉貴忠.小波變換多尺度地震波形反演.地球物理學(xué)報(bào)，1999，42（2）：241－248.Meng H Y，Liu G Z.Multiscale seismic waveform inversion by wavelet transform.Chinese J.Geophys.（in Chinese），1999，42（2）：241－248.

［17］陳勇，韓波，肖龍等.多尺度全變分法及其在時(shí)移地震中的應(yīng)用.地球物理學(xué)報(bào)，2010，53（8）：1883－1892.Chen Y，Han B，Xiao L，et al.Multiscale total variation method and its application on time－lapse seismic.Chinese J.Geophys.（in Chinese），2010，53（8）：1883－1892.

［18］Donoho D L.Superresolution via sparsity constraints.SIAM J.Math.Anal.，1992，23（5）：1309－1331.

［19］Candes E，Romberg J，Tao T.Stable signal recovery from incomplete and inaccurate measurements.Comm.Pure Appl.Math.，2006，59（8）：1207－1223.

［20］Loris I，Douma H，Nolet G，et al.Nonlinear regularization techniques for seismic tomography.J.Comput.Physics.，2010，229（3）：890－905.

［21］Jafarpour B，Goyal V K，Laughlin D B，et al.Transformdomain sparsity regularization for inverse problems in geosciences.Geophysics，2009，74（5）：R69－R83.

［22］Wang Y F.Seismic impedance inversion using l1norm regularization and gradient descent methods.J.Inv.Ill－Posed Problems，2011，18：823－838.

［23］Daubechies I，Defrise M，Mol C D.An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.Comm.Pure Appl.Math.，2004，57（11）：1413－1457.

［24］Ramlau R，Teschke G.A projection iteration for nonlinear operator equations with sparsity constraints.Numer.Math.，2006，104（2）：177－203.

［25］Kaltenbacher B，Schopfer F，Schuster T.Iterative methods for nonlinear ill－posed problems in Banach spaces：convergence and applications to parameter identification problems.Inverse Problems，2009，25（6）：065003.