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

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

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

0 引言

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

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

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

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

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ù)測(cè)多次波模型的矩陣形式；s為自適應(yīng)濾波算子。由于實(shí)際的地震數(shù)據(jù)并不是都能滿足輸出能量E(s)最小，即有效波與多次波正交，因此應(yīng)用此方法有可能會(huì)導(dǎo)致大的多次波剩余，或者衰減有效波的能量。

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

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

E(s)=|d-Ms|1

(2)

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

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

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

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

(3)

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

(4)

其中ri為第i個(gè)采樣點(diǎn)去除多次波后的殘差，即

r=d-Ms

(5)

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

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

(6)

對(duì)式(3)中的s求導(dǎo)并令其等于零并化簡(jiǎn)，根據(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的每個(gè)主對(duì)角線元素加上一個(gè)阻尼系數(shù)，則式(7)轉(zhuǎn)化為下式：

(A+λI)s=b

(8)

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

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

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

2.1 模型一

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

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

2.2 模型二

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

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

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

圖3 模型一壓制多次波前后炮集數(shù)據(jù)對(duì)比圖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 模型一壓制多次波前后零炮檢距 剖面對(duì)比圖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ù)對(duì)比圖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 模型二壓制多次波前后零炮檢 距剖面對(duì)比圖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)減方法具體實(shí)現(xiàn)，并進(jìn)行了模型數(shù)據(jù)與實(shí)際地震數(shù)據(jù)的處理，實(shí)現(xiàn)了自由表面相關(guān)多次波的迭代自適應(yīng)減。結(jié)果表明，文中方法對(duì)于復(fù)雜情況下的多次波壓制也具有較好的適應(yīng)性。通過(guò)對(duì)理論模型數(shù)據(jù)和實(shí)際數(shù)據(jù)多次波壓制效果對(duì)比可知，基于L1范數(shù)自適應(yīng)減方法能較好地匹配預(yù)測(cè)多次波模型和地震記錄中的實(shí)際多次波，從而能更好地壓制多次波，并保護(hù)有效波的能量。但同時(shí)，該方法的計(jì)算量比L2范數(shù)自適應(yīng)減方法相比要大，計(jì)算過(guò)程也較L2范數(shù)自適應(yīng)減方法復(fù)雜，但該方法把L1范數(shù)與L2范數(shù)自適應(yīng)減方法很好地結(jié)合起來(lái)，適應(yīng)性較強(qiáng)，且共軛梯度法具有收斂速度快的特點(diǎn)，綜合考慮之下認(rèn)為該方法可以被接受。在實(shí)際應(yīng)用中為獲得更好的壓制效果，應(yīng)注意以下方面：

圖7 實(shí)際數(shù)據(jù)壓制多次波前后炮集數(shù)據(jù)對(duì)比圖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í)際數(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ù)測(cè)多次波模型精度?；赟RME方法原理，由反饋迭代法得到預(yù)測(cè)的多次波模型。如果預(yù)測(cè)的多次波模型誤差較大，將會(huì)影響多次波的壓制效果或者導(dǎo)致有效波受損。通常在測(cè)線兩端，由于預(yù)測(cè)多次波的數(shù)據(jù)信息有限，容易導(dǎo)致預(yù)測(cè)多次波不準(zhǔn)確，從而影響多次波的壓制效果。

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

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

參考文獻(xiàn)：

[1] VERSCHUUR D J. 地震多次波去除技術(shù)的過(guò)去現(xiàn)在和未來(lái)[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é)報(bào)，2010，53(4)：963-973.

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

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

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