利用稀疏約束非平穩多項式回歸去除地震噪聲及拾取初至

劉國昌 蔡加銘 閆海洋 李潔麗 陳小宏

(①中國石油大學(北京)地球物理學院，北京 102249； ②東方地球物理公司研究院，河北涿州 072751； ③東方地球物理公司海洋物探處，天津 300450)

0 引言

多項式擬合回歸被廣泛用于地震數據隨機噪聲衰減、多次波壓制、品質因子估計及初至拾取靜校正等方面[1-7]。早在1988年俞壽朋等[8]利用信號橫向相干性，提出了多項式擬合加強有效信號的方法，并在疊后地震資料處理中獲得較好應用效果。田小平等[9]提出利用多項式擬合模板法消除地震信號的低頻干擾，分析了模板性質并給出了模板與低頻隨機干擾最高頻率之間的關系，實驗證明了方法的有效性。Johansen等[10]利用正交多項式提取AVO特征，很好地保護了剖面特征。李鯤鵬等[11]提出了小波變換過零點匹配和基于奇異值分解總體最小二乘(SVD-TLS)算法的多項式擬合去噪方法，提高了地震資料分辨率和信噪比。夏洪瑞等[12]提出了二次多項式擬合與中值約束的矢量分解方法，解決了常規正交多項式擬合過程中出現的斷點模糊問題。Lu等[1]利用滑動窗多項式擬合壓制局部線性噪聲，并用于衰減隨機噪聲以保護邊緣信息[13]。多項式擬合不僅用于疊后去噪，也可應用于疊前CMP去噪。薛亞茹等[14]利用多項式變換壓制動校正后CMP道集的隨機噪聲，既提高了信號和噪聲的分離效果，又有效地保護了地震信號的AVO信息。結合時頻變換的多項式擬合也被用于地震數據處理。李向云等[15]提出了改進的正交多項式變換，利用奇異值分解確定有效信號正交多項式系數譜的階數，結合小波變換減弱了有效信號和噪聲在低階的混疊，改善了隨機噪聲壓制效果。陸文凱[16]沿多次波旅行時軌跡，利用L1范數多項式擬合技術估計多次波，有效地消除了一次波和隨機噪聲對預測多次波的影響。在初至拾取靜校正方面，王輝等[17]針對常規靜校正方法的不足，提出利用多項式擬合單炮記錄初至時間進行地表靜校正，取得了較好效果。

Fomel[18-19]利用整形正則化技術研究非平穩回歸，針對非平穩信號給出了時變非平穩回歸系數估計方法。非平穩回歸已被用于多次波自適應相減[19]、隨機噪聲壓制[20]和時頻分析[21-22]等領域，應用效果較好。在非平穩回歸基礎上，Liu等[23]提出了基于非平穩多項式擬合的地震隨機噪聲壓制方法，通過自適應估計非平穩信號的相干分量，在壓制噪聲的同時更好地保護了有效信號。

非平穩多項式擬合是L2范數下的優化問題，盡管考慮了信號的時變特征，但是仍然假設殘差呈隨機分布，當地震數據中存在較強非隨機噪聲時，常規的基于L2范數的非平穩多項式擬合不再適用。因此，本文針對復雜稀疏分布殘差問題，在反問題正則化理論框架下，結合非平穩多項式回歸和L1范數約束，研究了稀疏約束非平穩多項式回歸理論與方法。首先，回顧了非平穩多項式回歸的基本原理；然后提出了L1范數約束的估計時變光滑多項式回歸系數和稀疏分布回歸殘差的方法，給出了算法基本流程和參數分析；最后，應用稀疏約束非平穩多項式回歸衰減噪聲及拾取初至，獲得了較好效果。

1 方法原理

1.1 非平穩多項式回歸

假設信號d(p)可以通過M階多項式近似，有

(1)

式中：pi為多項式自變量；bi為多項式系數；n為擬合誤差。如果n為隨機噪聲，通過擬合誤差能量最小的方法得到bi，即其為L2范數約束的最小二乘優化問題

(2)

式中‖·‖2表示L2范數。在L2范數約束下可以將多項式擬合問題改進為非平穩形式，非平穩多項式回歸假設多項式系數是時變的，則式(2)變為

(3)

可以看出，由于待估計的bi變為bi(p)后未知數個數變多，因此式(3)為欠定反問題，需要對其強加約束限制才能求解。在反問題研究領域，存在多種正則化形式實現反問題約束求解。若考慮經典的Tikhonov正則化，則式(3)在正則化約束條件下的最小平方優化問題為

(4)

式中：σ為正則化參數；Γ為Tikhonov正則化算子。為了方便起見，采用矩陣或算子的形式描述式(4)，即式(4)變為聯立最小平方問題

(5)

式中d、P和b分別為式(4)中的d、p和b對應的矩陣形式。式(5)為線性優化問題，通過共軛梯度法求解。本文選取整形正則化作為正則化方式，其類似于反問題預條件方法，在計算效率和參數選取方面具有一定優勢[19]。

1.2 稀疏約束非平穩多項式回歸

在常規非平穩多項式擬合中，假設噪聲呈隨機分布，采用L2范數約束能夠得到較好的結果。在噪聲復雜情況下，如存在強野值噪聲、涌浪噪聲等非隨機分布的噪聲，L2范數約束不再適合。在存在非隨機噪聲情況下，式(5)修正為

min∶‖d-Pb‖1

(6a)

min∶‖Γb‖2

(6b)

式中‖·‖1表示L1范數，表明擬合噪聲誤差是稀疏的。為更好地求解式(6)，在式(6a)中將殘差稀疏噪聲作為變量，即

min∶‖d-Pb-n‖2

(7a)

min∶‖Γb‖2

(7b)

min∶‖n‖1

(7c)

式中n為稀疏分布非隨機噪聲。式(7a)表示經多項式擬合后稀疏分布噪聲的誤差在最小平方意義下最小； 式(7b)表示對多項式系數的約束項，是對多項式系數的假設； 式(7c)表示噪聲的L1范數最小，表明噪聲滿足稀疏約束。在不同的情況下，式(7)可以具有不同的含義：

(1)只考慮式(7a)和式(7b)，可以估計L2范數約束的非平穩多項式系數，此時假設噪聲n滿足L2范數下最小。則式(7a)和式(7b)退化為式(5)，在理論上具有最小平方解，即

(8)

=[I+S(PTP-I)]-1SPTd

(9)

式中：S為整形算子，一般取三角平滑或者高斯平滑算子；I為單位矩陣。需要注意的是，式(9)是式(7a)和式(7b)的正則化約束理論解，采用共軛梯度類算法求解。為了簡便，將式(9)記為

(10)

式中G代表綜合線性算子， 包含整形正則化的形式。

(2)綜合考慮式(7a)～式(7c)，可以估計L1范數約束的非平穩多項式擬合。此時式(7)變為

min∶‖d-Gb-n‖2

(11a)

min∶‖n‖1

(11b)

算子G包含在L2范數約束下對b的平滑約束。式中未知量為b和n， 為了求解非線性優化問題，將式(11a)寫為

(12)

此時未知量變為(nb)T。假設噪聲n服從稀疏分布、b沿著時間方向是平滑的，則式(12)簡化為

d-Fx≈0

(13)

式中：F=(IG)；x=(nb)T。采用Daubechies迭代[24]求解式(13)，即

xk+1=Tλk[FTd+(I-FTF)xk]

(14)

式中Tλk為噪聲n的閾值算子，k為迭代次數。

將式(14)展開，得

(15)

由于Tλk僅僅對噪聲n進行閾值約束，而不對多項式系數b約束，則式(15)變為

(16)

式(16)為最終的L1范數約束的非平穩多項式回歸迭代公式。由于采用共軛梯度類算法求解L2范數整形正則化約束的非平穩多項式回歸公式(式(9))，因此在式(16)中可以直接用線性共軛梯度迭代幾次后的結果代替GT，即

bk+1=GTd-GTnk=GT(d-nk)≈G-1(d-nk)

(17)

1.3 稀疏約束非平穩多項式回歸算法描述

通過上述分析，可以得到稀疏約束非平穩多項式回歸的具體實現算法。

(1)設置迭代初始值。設初始噪聲n0=0，初始系數b0=0，P為多項式，非線性最大迭代次數為K，共軛梯度線性迭代次數為L。

(2)設非線性迭代次數k為1～K，計算bk+1=G-1(d-nk)。值得注意的是，具體實現是利用整形正則化約束共軛梯度算法實現的[19]，其求解為線性算法，最大線性共軛梯度迭代次數為L。

(3)計算nk+1=Tλi(d-Gbk)。即將殘差做閾值處理，以滿足殘差噪聲服從稀疏分布，解決擬合殘差中含有強稀疏分布噪聲的問題。

(4)重復步驟(2)、步驟(3)。迭代終止條件可以通過最大迭代次數K控制，也可以通過約束殘差的能量控制，在實際應用中前者更方便。

1.4 稀疏約束非平穩多項式回歸關鍵參數分析

在稀疏約束非平穩多項式回歸算法流程中，整形正則化平滑參數和閾值參數是較重要的2個參數。正則化平滑參數控制多項式系數的平滑性，它是非平穩系數估計的重要度量。平滑參數設置越大，求取的多項式系數越平滑，則多項式回歸系數隨自變量變化不大、越平穩；平滑參數設置越小，說明多項式系數時變性強、非平穩性越強。

閾值參數是稀疏約束反演的重要參數之一，閾值一般分為軟閾值和硬閾值。軟閾值的形式為

(18)

硬閾值的形式為[22]

(19)

閾值參數λk有不同的選取方式，如常數閾值、線性下降閾值、指數下降閾值等[25-26]，本文采用容易實現且收斂較快的百分比閾值[27]。

2 理論模型算例

圖1為含非隨機噪聲數據的平穩、非平穩多項式回歸結果。由圖可見：①含有非高斯分布噪聲的散點(圖1a)是在直線y=b0+b1x上加隨機噪聲和稀疏分布的野值得到的，基本符合直線y=b0+b1x分布，但是在x∈[25,40]時存在強野值噪聲，這些噪聲偏離了直線，且不符合隨機分布。②由于L2范數約束的平穩多項式回歸結果(圖1b紅線)存在野值，導致回歸直線向上方偏離，這是因為在回歸過程中引入了遠離真解的采樣點所致；采用L1范數稀疏約束的平穩多項式回歸結果較理想(圖1b黑線)，這是由于在L1范數約束下假設回歸殘差符合稀疏分布。③對比L2范數和L1范數回歸結果表明，L1范數更好地處理了含有稀疏分布噪聲的情況(圖1c黑線)。

為了測試非平穩多項式回歸效果，設計了一條曲線，并在該曲線上加隨機噪聲和稀疏分布的噪聲(圖2a)。圖2為含非隨機噪聲數據的非平穩多項式回歸結果。由圖可見：①平穩回歸(即直線)考慮了所有點的影響，并且呈直線分布(圖2b紅線)；L2范數約束非平穩回歸(圖2b粉線)考慮了曲線的時變特征，在x∈[0,20]時由于不存在稀疏分布的噪聲，回歸效果較好，但是在x∈[20,45]時存在稀疏分布的噪聲，回歸效果變差，表明異常值干擾了回歸結果。②對比不同回歸方法結果的絕對誤差發現，對于含有稀疏分布噪聲的情況，L1范數約束非平穩回歸擬合的曲線絕對誤差最小(圖2c黑線)。

圖1 含非隨機噪聲數據的平穩多項式(y=b0+b1x)回歸結果(a)含有非高斯分布噪聲的散點； (b)回歸結果； (c)回歸曲線與真實曲線(不含噪聲)的絕對誤差紅線代表L2范數約束，黑線代表L1范數稀疏約束

圖2 含非隨機噪聲數據的非平穩多項式(y=b0(x)+b1(x)x)回歸結果(a)含有非高斯分布噪聲的散點； (b)回歸結果； (c)回歸曲線與真實曲線(不含噪聲)的絕對誤差紅線代表L2范數約束平穩回歸，粉線代表L2范數約束非平穩回歸，黑線代表L1范數稀疏約束非平穩回歸

3 地震數據處理算例

3.1 地震CMP道集數據去噪

動校正后CMP道集上某一反射點的反射系數曲線(AVO曲線)可以通過多項式近似表示[28]。薛亞茹等[14]探討了多項式變換壓制CMP道集隨機噪聲的方法； Liu等[23]將其擴展到非平穩多項式的情形，研究了非平穩多項式回歸CMP道集去噪方法。上述方法在噪聲呈隨機分布時去噪效果較好，但如果存在稀疏分布的噪聲，如陸地地震數據的強野值噪聲、海洋地震數據的涌浪噪聲等，L2范數約束的非平穩多項式回歸不再適用。圖3為含強稀疏分布噪聲的CMP道集非平穩多項式回歸去噪結果。由圖可見： ①動校正后的CMP道集中存在強稀疏分布噪聲及隨機噪聲(圖3a)； ②L2范數約束去噪效果欠佳(圖3b)，L1范數約束去噪效果很好(圖3c)； ③L2范數約束去除的噪聲呈隨機分布(圖3d)，強噪聲“平均”影響了周圍信號，L1范數約束去除的噪聲中強噪聲基本呈稀疏分布(圖3e)。

圖4為不同方法去噪結果在旅行時為4.5s時刻的振幅。由圖可見：在強稀疏噪聲處(箭頭處)，L2范數約束方法去噪結果的振幅與原始數據不一致，L1范數約束方法去噪結果幾乎摒棄了強振幅的影響；在沒有強噪聲處(如50～80道處)，L1范數約束和L2范數約束的去噪結果幾乎一致，這是由于在回歸過程中充分考慮了非平穩特征所致。

圖3 含強稀疏分布噪聲的CMP道集非平穩多項式回歸去噪效果(a)含噪CMP道集； (b)L2范數約束去噪結果； (c)L1范數稀疏約束去噪結果； (d)L2范數約束去除的噪聲； (e)L1范數約束去除的噪聲

圖4 不同方法去噪結果在旅行時為4.5s時刻的振幅黑線為原始數據，紅線為L2范數約束去噪結果，藍線為L1范數約束去噪結果

3.2 地震初至自動拾取

地震初至拾取對建立近地表地震速度模型和靜校正非常重要，當前有很多初至拾取方法，其中能量比法是較常用的方法之一[29-32]，該方法計算滑動窗口內能量與累計能量的比值

(20)

式中：Ai為地震振幅；l為滑動窗口長度。可以看出，能量比值Rj的第1個峰值在前幾個采樣點，因為在前幾個采樣點累計能量與滑動窗口內能量差別不大。能量比值的第2個峰值一般對應初至[29]，但是當存在噪聲時，較難拾取第2個峰值，尤其存在強野值噪聲時拾取更不準確。因此，考慮采用稀疏約束多項式回歸方法處理強野值噪聲。

圖5為模擬的炮記錄及能量比。由圖可見：假設炮點和檢波點不在一條直線上，因此模擬的炮記錄的初至是一條雙曲線(圖5a)；為了測試方法的有效性，在模擬炮記錄上加隨機噪聲和強稀疏分布噪聲(箭頭所示)，強稀疏分布噪聲導致能量比第2個峰值與初至時刻差距較大(圖5b)，因此難以自動拾取初至。圖6為不同方法自動拾取的初至對比。由圖可見，噪聲的存在使自動拾取的第2個峰值(粉色線)與真實初至時刻有一定誤差，尤其對于稀疏噪聲存在的位置，獲得的初至結果不精確。分別利用自動拾取的初至位置進行L2范數和L1范數約束的非平穩多項式回歸發現：前者的初至更平滑，但是異常抖動也嚴重影響了L2范數回歸結果；后者通過稀疏約束減小了稀疏分布野值噪聲的影響，提高了初至拾取精度。

圖5 模擬的炮記錄(a)及能量比(b)

圖6 不同方法自動拾取的初至對比粉線為自動拾取能量比的第2個峰值得到的初至，綠線為L2范數非平穩回歸得到的初至，藍線為L1范數非平穩回歸得到的初至

在VSP數據處理中，初至拾取是非常重要的一步，利用初至信息可以直接求取速度、衰減參數及進行靜校正等處理。圖7為VSP數據及能量比。由圖可見： 在井源距大于0.5km處存在一些噪聲，導致能量比法初至拾取效果不好(圖7b)； 通過非平穩回歸可以改善初至拾取精度，尤其采用L1范數約束的非平穩回歸(圖7a的紅線)基本可以消除噪聲影響，可以獲得較可靠的初至。

圖8、圖9分別為陸地地震數據及其滑動時窗能量比剖面。由圖可見，由于實際數據的復雜性，自動拾取的初至效果不理想(圖8箭頭處)，經L1范數約束的非平穩回歸改善了初至拾取效果(圖8紅線)。

圖7 VSP數據(a)及能量比(b)黃線為能量比直接拾取的初至，綠線為L2范數拾取的初至，紅線為L1范數拾取的初至

圖8 陸地地震數據最大炮檢距為3km，共100炮。黃線為能量比直接拾取的初至，綠線為L2范數拾取的初至，紅線為L1范數拾取的初至。上為時間切片，左下為共炮點道集，右下為共炮檢距剖面

圖9 對應圖8的滑動時窗能量比剖面

圖10為常規初至拾取結果、稀疏約束非平穩多項式回歸初至拾取結果。由圖可見： ①遠炮檢距處由于信號較弱，信噪比低，兩種方法拾取的初至均存在異常點，但稀疏約束非平穩多項式回歸拾取的初至更平滑(圖10b)。②近炮檢距處由于存在野值等噪聲導致常規初至拾取結果存在異常值，效果較差(圖10a)；稀疏約束非平穩多項式回歸摒棄了稀疏分布噪聲的影響，初至拾取結果穩定、可靠(圖10b)。

需要說明的是，稀疏約束非平穩多項式回歸的初至拾取誤差由不規則稀疏分布噪聲引起，對地表高程變化劇烈或采樣不規則引起的初至劇烈抖動，初至拾取效果會受到影響。

圖10 常規初至拾取結果(a)、稀疏約束非平穩多項式回歸初至拾取結果(b)

4 討論與結論

(1)針對地震數據非平穩特征及噪聲稀疏分布問題，在稀疏約束反問題正則化理論框架下，提出了稀疏約束非平穩多項式回歸方法，并用于地震噪聲壓制和初至拾取。該方法采用整形正則化和L1范數聯合約束策略，利用共軛梯度和投影算法求解聯合約束反問題，同時估計具有時變光滑特征的多項式回歸系數和具有稀疏分布特征的回歸殘差，克服了稀疏分布強噪聲對反演的影響。

(2)整形正則化平滑參數和閾值參數是稀疏約束非平穩多項式回歸方法的兩個重要參數。正則化平滑參數控制多項式系數的平滑性，參數越小表明參數非平穩特征越強；閾值參數選取與噪聲的稀疏特征有關，本文采用容易實現且收斂較快的百分比閾值參數。

(3)與常規的多項式回歸相比，稀疏約束非平穩多項式回歸方法計算量略大，但由于數據處理過程中一般采用一維多項式回歸，因此計算量在可接受范圍內。