基于L1范數的多次波自適應減方法及應用分析

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

(中國地質大學 工程技術學院，北京 100083)

0 引言

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

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

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

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

1 方法原理

1.1 基于L2范數的多次波自適應減方法

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

(1)

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

1.2 基于L1范數的多次波自適應減方法

對于L1范數解，則不需要滿足有效波與多次波正交?；贚1范數的多次波自適應減方法的目標函數為式(2)[7]：

E(s)=|d-Ms|1

(2)

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

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

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

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

(3)

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

(4)

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

r=d-Ms

(5)

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

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

(6)

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

MTWTWMs=MTWTWd

(7)

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

(A+λI)s=b

(8)

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

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

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

圖2 模型一的介質模型Fig．2 The medium model of model 1

2 實例分析

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

2.1 模型一

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

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

2.2 模型二

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

2.3 實際地震數據

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

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

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

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

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

3 結論

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

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

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

(1)預測多次波模型精度?；赟RME方法原理，由反饋迭代法得到預測的多次波模型。如果預測的多次波模型誤差較大，將會影響多次波的壓制效果或者導致有效波受損。通常在測線兩端，由于預測多次波的數據信息有限，容易導致預測多次波不準確，從而影響多次波的壓制效果。

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

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

參考文獻：

[1] VERSCHUUR D J. 地震多次波去除技術的過去現在和未來[M]. 陳浩林，張保慶，劉軍,等,譯,北京：石油工業出版社，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] 李學聰，劉伊克，常旭，等. 均衡多道1范數匹配多次波衰減的方法與應用研究[J]. 地球物理學報，2010，53(4)：963-973.

[14] 劉光. 壓制自由表面相關多次波的自適應減方法研究及應用[D]. 北京: 中國地質大學，2009.

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

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