基于L1范數(shù)的多次波自適應(yīng)減方法及應(yīng)用分析

熊繁升， 黃新武， 高孝巧， 蔡雙霜， 雷海波， 何 江

(中國地質(zhì)大學(xué) 工程技術(shù)學(xué)院，北京 100083)

0 引言

多次波問題在地震勘探中普遍存在，如果不能對多次波進行有效地壓制，則將影響地震成像的真實性和可靠性，并干擾對有效波的識別，從而導(dǎo)致速度分析、疊前及疊后偏移的極大困難。因此，對多次波壓制方法及應(yīng)用進行研究不僅具有重要理論意義，而且有一定實際意義和應(yīng)用價值。針對多次波壓制問題，不同學(xué)者提出了不同的方法，可分為兩個主要類別：①基于一次波和多次波在空間上特性差異的去噪方法；②基于多次波的周期性和可預(yù)測性的去噪方法[1]。

由Verschuur等[2-4]提出的自由表面相關(guān)多次波壓制(SRME)方法是多次波壓制技術(shù)的突破，亦是目前壓制多次波方法和技術(shù)的主流。SRME方法使用一個匹配算子將通過反饋迭代法[1]預(yù)測得到的地震記錄的多次波模型從輸入數(shù)據(jù)中減去，其要求預(yù)測的多次波與實際數(shù)據(jù)中的多次波匹配較好，否則將直接影響多次波的壓制效果甚至損害有效波，因此做好多次波自適應(yīng)減對其壓制效果十分重要。

目前常用的多次波自適應(yīng)減方法有L2范數(shù)多次波自適應(yīng)減方法。該方法基于剩余能量最小原則，即有效波與多次波正交這一假設(shè)，而實際的地震數(shù)據(jù)并不是都能滿足該假設(shè)。同時，L2范數(shù)對于誤差較大的情況即出現(xiàn)異常值較為敏感[5]，因此應(yīng)用此方法有可能會導(dǎo)致大的多次波剩余，或者衰減有效波的能量。L1范數(shù)對于出現(xiàn)較大異常值時也能給出穩(wěn)定解[5-6]，且無須滿足有效波與多次波正交這一條件,從而可以利用L1范數(shù)來代替L2范數(shù)建立多次波自適應(yīng)匹配濾波器。

作者采用基于迭代重加權(quán)最小二乘法(IRLS算法)的混合L1/L2范數(shù)來近似L1范數(shù)解，該方法能有效地逼近Lp范數(shù)[6]。作者先就該方法的基本原理進行闡述，然后通過對復(fù)雜海底模型、Pluto 1.5模型數(shù)據(jù)和實際地震數(shù)據(jù)進行多次波壓制處理，并將處理結(jié)果與L2范數(shù)多次波自適應(yīng)減方法進行比較。此外還對影響多次波壓制效果的若干因素進行了簡要分析。

1 方法原理

1.1 基于L2范數(shù)的多次波自適應(yīng)減方法

目前常用的多次波自適應(yīng)減方法基于最小二乘準(zhǔn)則，即L2范數(shù)多次波自適應(yīng)減方法，其目標(biāo)是使輸入數(shù)據(jù)與濾波后的多次波之差的平方和即輸出能量E(s)達(dá)到最小[1]。L2范數(shù)多次波自適應(yīng)減方法的目標(biāo)函數(shù)為式(1)。

(1)

式中d為原始含多次波地震記錄；M為由反饋迭代法得到的預(yù)測多次波模型的矩陣形式；s為自適應(yīng)濾波算子。由于實際的地震數(shù)據(jù)并不是都能滿足輸出能量E(s)最小，即有效波與多次波正交，因此應(yīng)用此方法有可能會導(dǎo)致大的多次波剩余，或者衰減有效波的能量。

1.2 基于L1范數(shù)的多次波自適應(yīng)減方法

對于L1范數(shù)解，則不需要滿足有效波與多次波正交。基于L1范數(shù)的多次波自適應(yīng)減方法的目標(biāo)函數(shù)為式(2)[7]：

E(s)=|d-Ms|1

(2)

式中 各符號意義同式(1)。

秦鐵崖痛飲一杯，繼續(xù)道：“秦某就不同了。緝拿真兇、抓捕惡徒是我的職責(zé)，只要能抓到真兇，自己落敗受傷也在所不惜，對手是死是活更無足輕重，只要能不辱使命，就算完勝。再說，我是公門中人，受傷可以休養(yǎng)，名正言順。年年辦案，月月追兇，我也會厭煩，也會煩躁，暗地里也會發(fā)牢騷：這苦差，什么時候才是盡頭？這營生，什么時候才能歇一歇？幸好，每過一段時期，我就會受點小傷，正好歇一歇，讓狂躁厭煩之心安頓一下，才不至發(fā)瘋。更妙的是，受傷靜養(yǎng)，正好有空閑盤點前階段辦案得失，思考對敵手段中的疏漏。”

由式(2)可以看出基于L1范數(shù)多次波自適應(yīng)減方法的目標(biāo)函數(shù)E(s)的一階導(dǎo)數(shù)不是處處連續(xù)，因此不能運用高效的Levinson遞推算法確定濾波器，故在求解的過程中會帶來一定的困難。公式(2)中的最優(yōu)濾波器s需要用非線性反演方法來確定[1]。

由于L1范數(shù)的目標(biāo)函數(shù)一階導(dǎo)數(shù)不連續(xù)，為此人們提出了求解L1范數(shù)最小值的線性方法，比如Huber范數(shù)方法[8-9]。本文采用基于迭代重加權(quán)最小二乘算法(IRLS算法)的混合L1/ L2范數(shù)的方法來近似得到L1范數(shù)解。IRLS算法采用迭代過程求解一系列加權(quán)最小二乘問題，并在每次迭代過程中按照一定規(guī)則調(diào)整權(quán)重，以逐步逼近最優(yōu)解，其具有收斂速度快的特點；同時，采用此方法可將非線性問題的求解，轉(zhuǎn)化為分段線性問題的求解[6,10]。此方法是將L1范數(shù)轉(zhuǎn)化為加權(quán)后的L2范數(shù)，其目標(biāo)函數(shù)為式(3)

(3)

式中W為加權(quán)矩陣，是對角矩陣，主對角線上的元素根據(jù)式(4)確定[17]：

(4)

其中ri為第i個采樣點去除多次波后的殘差，即

r=d-Ms

(5)

在式(4)中，ε是一個非常重要的先驗值參數(shù)，其大小對加權(quán)矩陣主對角線上元素數(shù)值有影響，進而影響求解濾波算子s時方程組的系數(shù)矩陣。由式(4)可知，加權(quán)矩陣W中元素的數(shù)值與殘差ri有關(guān)，由式(5)知殘差ri中又包含算子s，故式(3)是非線性的，不能用一般線性方法求解。

由式(4)和式(5)可知：給定ε，若殘差ri小于ε，二者相差越大，wii越大且wii越接近“1”，則矩陣W接近于單位對角矩陣；當(dāng)殘差ri趨于零時，即剩余能量最小，式(3)等價于式(1)，即此時的解相當(dāng)于L2范數(shù)解；若殘差ri大于ε，二者相差越大，式(4)近似等于(ri/ε)-1/2，即此時的解為L1范數(shù)解。通過合適的ε的可以實現(xiàn)L1范數(shù)與L2范數(shù)的平滑過渡[10-11]。ε通常取值如下[11-12]：

(6)

對式(3)中的s求導(dǎo)并令其等于零并化簡，根據(jù)基于IRLS算法的混合L1/L2范數(shù)準(zhǔn)則，求解使目標(biāo)函數(shù)E(s)最小可轉(zhuǎn)化為求解以下含濾波算子s的方程組：

MTWTWMs=MTWTWd

(7)

記矩陣A=MTWTWM，向量b=MTWTWd。為防止系數(shù)矩陣A成為奇異矩陣，在矩陣A的每個主對角線元素加上一個阻尼系數(shù)，則式(7)轉(zhuǎn)化為下式：

(A+λI)s=b

(8)

式中λ為阻尼系數(shù)，一般取0.001～0.002。

對于式(8)的求解，可采用牛頓迭代法，李學(xué)聰?shù)萚13]采用了直接迭代法，這里采用共軛梯度算法進行求解，因其收斂速度較快，且得到的解較為穩(wěn)定。基于IRLS算法的混合L1/L2范數(shù)單道多次波自適應(yīng)減方法的實現(xiàn)過程簡單描述如下：首先設(shè)定IRLS算法、共軛梯度法的求解精度、最大迭代次數(shù)(以防止解不收斂)以及濾波算子s的長度；根據(jù)SRME的迭代實現(xiàn)方法原理，得到預(yù)測的多次波模型；其次設(shè)定各初始值并代入式(8)由共軛梯度法求解得到濾波算子s；然后計算當(dāng)前剩余誤差r并判斷是否達(dá)到精度要求；最后如果滿足剩余誤差r滿足精度要求，則輸出該道處理結(jié)果做迭代處理，不滿足，則修正權(quán)矩陣的數(shù)值重復(fù)求解。整個實現(xiàn)過程如圖1所示。

圖1 混合L1/L2范數(shù)多次波自適應(yīng)減方法流程圖Fig．1 Flow chart of adaptive multiple subtraction using hybrid L1/L2 norm

圖2 模型一的介質(zhì)模型Fig．2 The medium model of model 1

2 實例分析

將基于L1范數(shù)的多次波自適應(yīng)減方法應(yīng)用于兩個模型數(shù)據(jù)和實際地震數(shù)據(jù)，分析多次波壓制效果，并與L2范數(shù)多次波自適應(yīng)減方法作比較。

2.1 模型一

模型一為一復(fù)雜海洋地質(zhì)模型，海底及其下伏地層界面橫向起伏較大(圖2)，采用聲波方程的有限差分法產(chǎn)生人工合成地震記錄，分別采用吸收和反射邊界條件處理模型邊界[14-16]。由于產(chǎn)生多次波的介質(zhì)系統(tǒng)較為復(fù)雜，多次波在橫向上的變化較大，故此模型可有效檢驗多次波壓制方法。本數(shù)據(jù)共141炮，每炮101道，每道501個采樣點，4 ms采樣。

圖3(a)為第1、第30、第60、第90、第125炮共五炮原始含多次波地震記錄，可以看出，多次波形式較復(fù)雜，且部分多次波與有效波重疊較大。為使L2范數(shù)和L1范數(shù)多次波自適應(yīng)減方法具有可比性，兩者算子長度均取15個采樣點，分相同的兩個時窗進行單道處理，在第一個時窗內(nèi)僅處理一階多次波；同時在L1范數(shù)多次波自適應(yīng)減方法中，共軛梯度法與IRLS算法的控制精度均設(shè)為10-3、最大迭代次數(shù)均為100。由圖3(b)、(c)、(d)可以看出，采用L2范數(shù)多次波自適應(yīng)減方法雖然壓制了部分多次波，但同時引入了一些干擾波；而采用L1范數(shù)多次波自適應(yīng)減方法在多炮記錄和零炮檢距剖面上多次波的大部分能量均已得到有效壓制，具有較強的適應(yīng)性，且有效波和多次波重疊在一起時，也能獲得較好地壓制效果。在零炮檢距剖面(圖4)可以看出，箭頭所指示的多次波L1范數(shù)自適應(yīng)減方法較L2范數(shù)自適應(yīng)減方法壓制效果好。

2.2 模型二

SMARRT Pluto1.5模型數(shù)據(jù)是檢驗多次波壓制的一套標(biāo)準(zhǔn)數(shù)據(jù)。作者選取部分Pluto1.5模型數(shù)據(jù)測試，共500炮，每炮含200道，每道601個采樣點。圖5(a)為Pluto 1.5模型的第1、第200、第400炮共三炮原始地震記錄。由于該模型較復(fù)雜，一階多次波能量較強且原始數(shù)據(jù)中存在有效波與多次波交叉的情況。對本數(shù)據(jù)處理時，L2范數(shù)與L1范數(shù)的多次波自適應(yīng)減方法的算子長度均取20個采樣點，分相同的兩個時窗進行單道處理，在第一個時窗內(nèi)只處理一階多次波，迭代3次；在L1范數(shù)的多次波自適應(yīng)減方法中，共軛梯度法和IRLS算法的控制精度均為10-3、最大迭代次數(shù)均設(shè)為100。圖5(b)為L2范數(shù)多次波自適應(yīng)減結(jié)果，可以看到，對一階多次波L2范數(shù)自適應(yīng)減方法壓制效果較好，但同時損害了有效波；而基于L1范數(shù)自適應(yīng)減方法減去后不僅多次波的壓制較為干凈，且有效波能量保護得較好。由圖6也可以看到，L2范數(shù)自適應(yīng)減方法與L1范數(shù)自適應(yīng)減方法相比有更大的多次波剩余。

2.3 實際地震數(shù)據(jù)

這里選取一個單邊放炮的二維實際海洋地震數(shù)據(jù)進行多次波的壓制處理，以進一步驗證此方法的正確性和適用性。本數(shù)據(jù)炮點距為50 m，檢波點距為25 m，最小炮檢距為 200 m，采用插值方法重構(gòu)近偏移距地震道，對所缺失的八個地震道進行外推重建。圖7(a)為預(yù)處理及近炮檢距數(shù)據(jù)重建后三炮原始記錄，由圖7(a)可見，此數(shù)據(jù)多次波能量較強，且其表現(xiàn)形式復(fù)雜，部分多次波和有效波疊加在一起，既有一階多次波，還有高階多次波；另外，還有層間多次波。在本次處理中，L2范數(shù)和L1范數(shù)多次波自適應(yīng)減方法的算子長度均取25個采樣點，分相同時窗進行單道處理，迭代3次；在L1范數(shù)多次波自適應(yīng)減方法中，共軛梯度法和IRLS算法的控制精度均設(shè)為10-3、最大迭代次數(shù)均設(shè)為100。圖7為壓制多次波前后炮集記錄對比圖，重建后的數(shù)據(jù)含有最小偏移距。圖8為壓制多次波前后最小炮檢距剖面對比圖，可以看到，箭頭所指示處的多次波由L2范數(shù)多次波自適應(yīng)減方法處理之后仍匹配不足，尤其是一階多次波；而基于L1范數(shù)自適應(yīng)減方法處理后大部分多次波被壓制，且有效波的大部分能量保護得較好。

圖3 模型一壓制多次波前后炮集數(shù)據(jù)對比圖Fig．3 Results of shot gather data by various multiple elimination method of model 1(a)原始五炮記錄；(b)L2范數(shù)自適應(yīng)減結(jié)果；(c)L1范數(shù)自適應(yīng)減結(jié)果；(d) L1范數(shù)自適應(yīng)減方法減掉的多次波

圖4 模型一壓制多次波前后零炮檢距 剖面對比圖Fig．4 Results of zero offset section of model 1 by various multiple elimination method(a) 原始含多次波零炮檢距記錄；(b) L2范數(shù)自適應(yīng)減結(jié)果；(c) L2范數(shù)自適應(yīng)減方法減掉的多次波； (d)L1范數(shù)自適應(yīng)減結(jié)果；(e) L1范數(shù)自適應(yīng)減方法減掉的多次波

圖5 模型二壓制多次波前后炮集數(shù)據(jù)對比圖Fig．5 Results of shot gather data by various multiple elimination method of model 2(a)原始三炮記錄；(b)L2范數(shù)自適應(yīng)減結(jié)果；(c)L1范數(shù)自適應(yīng)減結(jié)果；(d)L1范數(shù)自適應(yīng)減方法減掉的多次波

圖6 模型二壓制多次波前后零炮檢 距剖面對比圖Fig．6 Results of zero offset section of model 2 by various multiple elimination method(a)原始含多次波的零炮檢距記錄；(b)L2范數(shù)自適應(yīng)減結(jié)果；(c)L2范數(shù)自適應(yīng)減方法減掉的多次波； (d)L1范數(shù)自適應(yīng)減結(jié)果；(e) L1范數(shù)自適應(yīng)減方法減掉的多次波

3 結(jié)論

應(yīng)用基于IRLS算法的混合L1/L2范數(shù)方法，將L1范數(shù)多次波自適應(yīng)減方法具體實現(xiàn)，并進行了模型數(shù)據(jù)與實際地震數(shù)據(jù)的處理，實現(xiàn)了自由表面相關(guān)多次波的迭代自適應(yīng)減。結(jié)果表明，文中方法對于復(fù)雜情況下的多次波壓制也具有較好的適應(yīng)性。通過對理論模型數(shù)據(jù)和實際數(shù)據(jù)多次波壓制效果對比可知，基于L1范數(shù)自適應(yīng)減方法能較好地匹配預(yù)測多次波模型和地震記錄中的實際多次波，從而能更好地壓制多次波，并保護有效波的能量。但同時，該方法的計算量比L2范數(shù)自適應(yīng)減方法相比要大，計算過程也較L2范數(shù)自適應(yīng)減方法復(fù)雜，但該方法把L1范數(shù)與L2范數(shù)自適應(yīng)減方法很好地結(jié)合起來，適應(yīng)性較強，且共軛梯度法具有收斂速度快的特點，綜合考慮之下認(rèn)為該方法可以被接受。在實際應(yīng)用中為獲得更好的壓制效果，應(yīng)注意以下方面：

圖7 實際數(shù)據(jù)壓制多次波前后炮集數(shù)據(jù)對比圖Fig．7 Results of shot gather data by various multiple elimination method of real data(a)原始三炮含多次波炮集；(b)L2范數(shù)自適應(yīng)減結(jié)果；(c)L1范數(shù)自適應(yīng)減結(jié)果；(d)L1范數(shù)自適應(yīng)減方法減掉的多次波

圖8 實際數(shù)據(jù)壓制多次波前后零炮檢 距剖面圖Fig．8 Results of zero offset section of real data by various multiple elimination method(a)原始含多次波的零炮檢距記錄；(b)L2范數(shù)自適應(yīng)減結(jié)果；(c)L1范數(shù)自適應(yīng)減結(jié)果

(1)預(yù)測多次波模型精度。基于SRME方法原理，由反饋迭代法得到預(yù)測的多次波模型。如果預(yù)測的多次波模型誤差較大，將會影響多次波的壓制效果或者導(dǎo)致有效波受損。通常在測線兩端，由于預(yù)測多次波的數(shù)據(jù)信息有限，容易導(dǎo)致預(yù)測多次波不準(zhǔn)確，從而影響多次波的壓制效果。

(2)濾波算子長度。濾波算子s的長度不宜過大(一般不超過35個采樣點)，否則可能會導(dǎo)致有效波受損；但也不能過小(一般不小于8個采樣點)，否則可能會導(dǎo)致匹配不足而有較大的多次波剩余，其一般在10到30個采樣點之間取值。

(3)時窗長度及其控制。對各道地震記錄分時窗處理，考慮了多次波關(guān)于時間上的變化，可使多次波壓制效果更好。關(guān)于時窗長度的選擇不是越短越好(例如一個階次的多次波分多個時窗處理)，否則會導(dǎo)致有效波受損嚴(yán)重；時窗長度也不宜取得過長(例如一個時窗內(nèi)處理多個階次的多次波)，否則會導(dǎo)致多次波去除不干凈，一般是在一個時窗內(nèi)處理一個階次的多次波。

參考文獻(xiàn)：

[1] VERSCHUUR D J. 地震多次波去除技術(shù)的過去現(xiàn)在和未來[M]. 陳浩林，張保慶，劉軍,等,譯,北京：石油工業(yè)出版社，2010.

[2] VERSCHUUR D J. Surface-related multiple elimination: an inversion approach: Ph.D.thesis[M]. Delft: University of Technology, 1991.

[3] VERSCHUUR D J, BERKHOUT A J, WAPENAAR C P A. Adaptive surface-related multiple elimination [J]. Geophysics, 1992, 57(9):1166-1177.

[4] VERSCHUUR D J, BERKHOUT A J. Estimation of multiple scattering by iterative inversion ,partⅠ:theoretical considerations[J].Geophysics,1997,62(5):1586-1595.

[5] CLAERBOUT J F, MUIT F. Robust modeling with erratic data[J]. Geophysics,1973,38: 820-844.

[6] BUDE K P, LANGAN R T. Hybrid l1/l2 minimization with applications to tomography[J].Geophysics,1997,62:1183-1195.

[7] GUITTON A , VERSCHUUR D. J. Adaptive subtraction of multiples using the L1-norm[J]. Geophysical Prospecting, 2004, 52: 27-38.

[8] GUITTOH A .Implementation of a nonlinear solver for minimizing the Huber norm[J].SEP,2000,103:281-289.

[9] Antoine Guitton. Huber solver versus IRLS algorithm for quasi L1inversion[J]. Stanford Exploration Project, Report 103, April 27, 2000, 103:205-266.

[10] ADAM GERSZTENKOR, BEE BEDNAR J, LARRY LINES R. Robust iterative inversion for the one-dimensional acoustic wave equation[J]. Geophysics, 1986,51(2):357-368.

[11] DARCHE G. Iterative L1 deconvolution[J]. SEP Annual Report , 1989,61: 281-301.

[12] WEGLEIN A B. Multiple attenuation: an overview of recent advances and the road ahead[J].The Leading Edge, 1999, 18(1): 40- 44.

[13] 李學(xué)聰，劉伊克，常旭，等. 均衡多道1范數(shù)匹配多次波衰減的方法與應(yīng)用研究[J]. 地球物理學(xué)報，2010，53(4)：963-973.

[14] 劉光. 壓制自由表面相關(guān)多次波的自適應(yīng)減方法研究及應(yīng)用[D]. 北京: 中國地質(zhì)大學(xué)，2009.

[15] 黃新武，牛濱華，劉光,等. 基于數(shù)據(jù)一致性原理預(yù)測與壓制自由表面多次波的效果分析[J]. 石油地球物理勘探，2009， 44 (4)：409- 416.

[16] 牛濱華，沈操，黃新武. 波動方程多次波壓制技術(shù)的進展[J]. 地球物理學(xué)進展，2002，17 (3) ：480- 485.