地震隨機噪聲壓縮感知迭代壓制方法

劉璐， 劉洋， 劉財， 鄭植升

吉林大學地球探測科學與技術學院， 長春 130026

0 引言

隨機噪聲是地震勘探無法避免的噪聲之一，其頻率、視速度、方向、振幅均不固定，往往廣泛分布于有效信號頻帶內，嚴重影響著地震記錄的信噪比以及后續的處理和解釋工作.隨著地震勘探的目標逐漸轉向“雙復雜”(復雜地表條件和復雜介質)以及深層條件下的資源探查，隨機噪聲對于復雜信號和弱信號有效識別的影響更加突出，因此，如何有效壓制隨機噪聲，提升地震記錄信噪比，恢復更高質量的地震有效信號振幅信息尤為重要.隨機噪聲主要來源于環境噪聲、風動、記錄儀器以及檢波器與地面耦合差等情況，通常被認為是平穩、高斯隨機過程(李慶忠，1993)，而近年來的研究發現地震勘探中的隨機噪聲并不是傳統意義上的平穩隨機過程，其平穩性隨記錄時長增加而降低，且與采集環境也有關系，尤其在山區等地表條件復雜的地區常采集到非平穩隨機噪聲記錄(鐘鐵等，2017).因此，開發能夠保護復雜地震有效信號的非平穩隨機噪聲壓制方法具有重要的現實意義.

目前壓制地震隨機噪聲的方法有多種類型，包括：(1)數據域濾波類的噪聲壓制方法；(2)f-x域及t-x域預測濾波類方法；(3)基于深度學習的地震數據去噪方法；(4)變換類噪聲壓制方法；在數據域濾波的地震數據去噪方法中，中值濾波與均值濾波是最常見的方法，Liu等(2006)在傳統中值濾波方法的基礎上發展了二維多級中值濾波方法，相對于一維中值濾波，該方法能夠有效保護信號細節；Bonar 和Sacchi(2012)提出非局部均值濾波方法，利用數據點間的結構相似性以及隨機噪聲的非相關性進行去噪，具有良好的結構和細節保持能力.預測濾波方法中，Amba和Claerbout(1997)介紹并比較了f-x域及t-x域兩種同相軸預測方法，Liu和Li(2018)及國朧予等(2020)在傳統f-x域及t-x域預測濾波中引入流式處理框架來解決地震數據隨機噪聲的壓制問題，實現了計算效率的提升.基于深度學習的去噪方法在近年來發展迅速，并在應用中形成了多種網絡結構，地震數據去噪中常見的有U-net、Res-net等基于神經網絡的去噪方法(Liu et al.，2018；李海山等，2020).基于稀疏變換的地震數據去噪方法主要利用地震數據的稀疏特性進行信噪分離，自20世紀80年代，地震數據的稀疏特性就被應用于解決地震隨機噪聲的壓制問題，小波變換(Mallat and Zhong，1992)克服了傅里葉變換(Canales，1984；Vassiliou and Garossino，1998)在非平穩信號應用中的缺陷，實現了多頻率的分解，小波變換至今在地震數據隨機噪聲的壓制中仍是一種應用廣泛的變換方法(高靜懷等，2006；吳招財和劉天佑，2008)，近年來，針對傳統小波變換難以表征高維數據方向性等問題，發展出了多尺度幾何分析方法，如ridgelet(Donoho，2000)變換、curvelet變換(Herrmann et al.，2008)、seislet變換、shearlet變換等都在實際應用中取得了良好的效果，其中，seislet變換及shearlet變換以其良好的稀疏性在地震數據的去噪(劉洋等，2009；鄧盾等，2019；童思友等，2019)中得到了較多應用，以地震數據稀疏特性為基礎的去噪方法是近年來的研究方向之一，其能夠提供比傳統方法更加有效的弱信號恢復能力.隨著數學方法以及信號處理等相關領域的不斷發展，地震隨機噪聲壓制方法也在不斷地更新，目標是在時空變地震信號保真的前提下提高信噪比，為解決實際應用中的復雜問題提供更多的方法選擇.

壓縮感知理論由Candès(2006)提出，并在圖像處理領域得到了廣泛的應用.壓縮感知理論能夠突破奈奎斯特(Nyquist)采樣定理的限制，即使采樣頻率沒有達到信號最高頻率的二倍，也可以將原始信號精確地恢復出來，因此，基于壓縮感知理論的方法在地震數據的插值、去噪等方面能夠取得良好的應用效果.楊冠雨等(2020)在壓縮感知中結合shearlet變換，并使用雙正則化項進行約束，實現了地震數據的插值重建，南方舟等(2018)將壓縮感知與curvelet變換相結合，實現了海底地震儀采集數據的去噪.壓縮感知理論下隨機噪聲的壓制方法主要包括兩個方面：一是稀疏變換基的選擇，二是重構方法的選擇.由于地震數據具有特殊性，常規圖像處理領域的變換方法往往無法為地震數據提供理想的稀疏表征，因此Fomel和Liu(2010)提出了有效壓縮地震數據同相軸的seislet變換，seislet變換屬于類小波變換，是在小波提升算法的基礎上結合地震數據局部傾角等屬性構建的新方法，能夠針對地震數據特性實現更優的稀疏表征.劉洋等(2009)將低階seislet變換擴展到高階，Liu和Fomel(2010)將seislet變換與波動方程炮檢距連續算子結合構造了OC-seislet變換，Liu等(2015,2017)利用動校正方程構造了依賴速度的VD-seislet變換以及廣義VD-seislet變換，試圖降低復雜地質環境及隨機噪聲對地震數據屬性表征的影響，增強seislet變換對地震數據的稀疏表征能力.有效的稀疏表征是壓縮感知去噪效果的前提和保證，由于有效信號具有可壓縮性，而隨機噪聲不具有可壓縮性，經稀疏表征后的有效信號和隨機噪聲可以在重構過程中分離.常見的重構算法有貪婪類方法和凸優化方法.貪婪類方法的基本思想是將全局優化問題轉化為分階段優化，在每一次迭代中只考慮當前階段的最優解，算法時空復雜度低，但不能保證得到全局最優解，貪婪方法以匹配追蹤類(Mallat and Zhang，1993；陳興飛和孫紅梅，2019)為代表.凸優化方法主要包括內點法(Kim et al.，2007)、梯度投影法(Figueiredo et al.，2007)、全變差方法(Chen et al.，2001)、Bregman迭代(Yin et al.，2008)、迭代硬閾值和迭代軟閾值(Pope，2008；Daubechies et al.，2004)等.在凸優化方法中，迭代收縮閾值方法由于其計算復雜度低，適用于高維大規模數據的性質受到了廣泛的應用.盡管壓縮感知理論在地震數據隨機噪聲壓制方面有一些應用，但是其去噪數學反問題和匹配的實現方法并不明確.

本文針對地震數據隨機噪聲壓制問題，建立完善的壓縮感知框架下稀疏信號重建數學模型，通過求解基追蹤降噪模型下的去噪反演目標函數，提出基于迭代軟閾值算法的“采集-重建-修復”方案，對含隨機噪聲數據中的有效信號進行迭代分離.通過選取和對比不同的稀疏變換基，論證有效信號的稀疏表征能力對于提高地震隨機噪聲壓縮感知迭代壓制方法有效性的重要意義，同時，針對稀疏優化問題中正則化項使有效信號衰減的現象，利用除偏算法對損失的振幅進行恢復，進一步提高本文方法的弱信號保護能力.將本文方法與工業標準隨機噪聲壓制FXDECON方法進行比較，模型及實際數據測試結果表明，壓縮感知迭代去噪方法可以有效地壓制地震勘探記錄中的較強振幅隨機噪聲，有效保護復雜構造的同相軸信息.

1 基本理論

1.1 壓縮感知基本理論回顧

假設一組有限長的離散信號x∈RN,在某個變換域內可以線性表示為x=Ψc，式中c是x在Ψ域的變換域系數，Ψ是N×N的稀疏變換基，也用于表示稀疏反變換.當信號在某個變換基Ψ下的系數矩陣中絕大部分元素為零或約等于零時，就稱該信號是變換域內的稀疏信號.信號x本身或者變換域系數c具有稀疏性是壓縮感知理論應用的基礎.

當信號本身是稀疏的，壓縮感知的一般采樣問題可以描述為(Candès and Wakin，2008)：

y=Φx，

(1)

式中，x∈RN為原始信號，觀測矩陣Φ數據維度為M×N，當M?N時，可以獲得關于x的欠采樣信號y∈RM，壓縮感知理論主要用于實現利用y重建x的過程.

在地震數據中，有效信號本身一般不具有稀疏性，但是可以經過稀疏變換進行稀疏表征，即x=Ψc，此時壓縮感知問題為

y=ΦΨc=Θc，

(2)

式中，Θ=ΦΨ稱為傳感矩陣.

此時，由觀測信號y恢復原始信號x的最優化問題可以描述為

min{‖c‖0:y=Φx}，

(3)

信號恢復的過程利用壓縮感知中的重構算法.通過壓縮感知技術來觀測數據，可以節省大量采樣空間，且能夠準確重建不滿足采樣定理條件的數據.壓縮感知理論不僅可以有效解決稀疏采樣的恢復問題，還能夠用于壓制數據中的隨機噪聲干擾.

1.2 壓縮感知框架下的地震數據隨機噪聲壓制問題

信號中含有隨機噪聲干擾時，公式(1)變為

y=Φd=Φ(x+z)=Φx+z′，

(4)

式中，d為含隨機噪聲地震數據，x為不含噪聲的純信號，z為隨機噪聲，z′為誤差項，即壓縮感知采樣后的噪聲，y為壓縮感知的觀測信號.研究如何從觀測信號y中恢復有效信號x成為壓縮感知框架下的地震數據隨機噪聲壓制問題(Candès and Wakin，2008).

壓縮感知最優化問題可以描述為公式(3)所示的L0范數最小化的最優化問題，但該問題是一種典型的非凸NP-hard問題，在數值計算上難以有效求解，一種可行的方法是用L1范數替代問題中的L0范數.當信號滿足在變換域內稀疏，稀疏變換基Ψ和觀測矩陣Φ滿足非相干性或等距約束性時，最優化問題可以由L0范數最小化轉化為L1范數最小化，在基追蹤(Basis Pursuit，BP)問題(Chen et al.，2001)下進行求解：

min{‖c‖1∶y=Φx}.

(5)

基追蹤準則下的最優解意味著所得解的L1范數最小.在公式(4)數據中含有噪聲的條件下，求解壓縮感知框架下的隨機噪聲壓制問題需要對基追蹤問題(5)的約束條件進行松弛近似，得到基追蹤降噪(Basis Pursuit De-Noising，BPDN)問題(Chen et al.，2001)，該問題可以描述為

(6)

以及相應的非約束問題(Figueiredo et al.，2007)：

(7)

迭代軟閾值方法由非迭代的軟閾值方法發展而來，Donoho等(1995)提出軟閾值去噪方法，該方法通過一個簡單的閾值步驟將含噪聲數據中的有效信號與噪聲分離，且估計的去噪結果具有較好的平滑性，但是軟閾值去噪的目標函數與基追蹤降噪的更一般化目標函數并不相同，因此軟閾值方法不能解決基追蹤降噪問題.為了適應多種稀疏變換下從含噪聲數據中估計有效信號的線性反問題，Daubechies等(2004)提出1≤p≤2時解決Lp問題的迭代軟閾值方法，并且給出任意初值時的收斂性證明，為迭代軟閾值方法的應用提供了理論基礎.

迭代軟閾值將軟閾值函數與Landweber迭代(Landweber，1951)相結合，在迭代的每一步中應用軟閾值函數.迭代軟閾值方法可以處理多種Lp范數(1≤p≤2)問題，當取L1范數時，迭代軟閾值方法的目標函數與基追蹤降噪問題(7)一致，因此，迭代軟閾值方法可以用于解決壓縮感知框架下的隨機噪聲壓制問題.利用迭代軟閾值方法求解壓縮感知問題的迭代解如公式(8)所示，迭代初始值任意的情況下公式(8)最終收斂于公式(7)的解：

xn+1=Sλ[xn+ΦT(y-Φxn)]，

(8)

其中：

Sλ[·]=ΨsλΨ-1[·]，

(9)

[·]T表示為矩陣的轉置，Ψ-1為稀疏正變換，sλ是對信號中的每一個元素執行軟閾值，軟閾值函數定義為

(10)

壓縮感知的信號重構效果不僅與重構算法的選擇有關，還取決于稀疏變換及觀測矩陣的選擇.選擇的稀疏變換應該使有效信號經變換后的變換域系數足夠稀疏，實現有效信號和隨機噪聲的變換域系數振幅差異最大化，保證重構信號的精度.離散小波變換在科學和工程上有著廣泛的應用，其主要基于分片平滑假設，在地震數據隨機噪聲壓制方面有著廣泛的應用.但傳統小波變換在方向性上缺乏靈活性，同時地震數據往往不滿足其應用的分片平滑假設，因此，利用圖像中方向特性的類小波變換得到了持續的發展，seislet變換(Fomel and Liu，2010)可以根據地震數據的不同特征對地震數據進行更加有效地壓縮，因此本文選取離散小波變換和seislet變換作為壓縮感知迭代噪聲壓制方法中的稀疏變換基，并對比處理效果.

觀測矩陣的設計是另一個影響信號重建精度的重要環節.壓縮感知理論本質上是一個欠定問題，為使觀測信號y中包含有足夠多原始信號的信息以便進行信號的重構，觀測矩陣與稀疏變換基需要滿足非相干性，相干性越小，欠采樣的觀測信號y中包含的有效信息越多，才能夠保證信號重建的精確性.在觀測矩陣元素的選取方面，Candès和Wakin(2008)給出了隨機生成的方案.在地震數據的插值問題中，觀測矩陣與觀測系統的設計有關，單位矩陣是在全采樣的情況下對觀測系統的數學表示，欠采樣時，插值問題的觀測矩陣為對角線隨機缺失的單位陣.在去噪問題中，觀測矩陣與觀測系統無關，主要考慮計算上的優勢.高斯隨機矩陣作為觀測矩陣與正交稀疏變換基滿足非相干性，提供了信號恢復的前提，結構簡單易于構造且具有普適性，是壓縮感知方法中常用的觀測矩陣，因此本文選取高斯隨機矩陣作為壓縮感知迭代噪聲壓制方法中的觀測矩陣，且有Φ～N(0,1/M)，滿足均值為0，方差為1/M的高斯分布，其中M為Φ的行維度.

1.3 L1正則化振幅保真的除偏方法

最優化問題(6)和(7)中由于L1正則化項的存在，重構產生的結果往往在幅度上相比于原信號通常會有一定的損失，可以選擇除偏(debiasing)算法(Wright et al.，2008)來恢復L1正則項引起的信號幅度損失.

(11)

本文研究表明，當選擇的稀疏變換能夠提供高效壓縮結果時，L1正則化項的影響將變得較小，有效信號的損失也較小，此時除偏步驟可以被省略，因為除偏方法在恢復有效信號的同時會恢復出一定的噪聲，降低處理結果的信噪比.除偏步驟主要應用于對信號的保真性要求更高的處理任務中.

1.4 壓縮感知迭代噪聲壓制方法的實現流程

在壓縮感知框架中使用高斯隨機矩陣與稀疏變換的組合作為壓縮感知的傳感矩陣，對含噪聲信號進行“采集-重建-修復”的壓縮感知迭代噪聲壓制過程可以歸納如下3個步驟：

(1)帶有隨機噪聲的實際數據為原始信號d，根據公式(4)，利用高斯隨機矩陣Φ對信號d進行采樣得到觀測信號y=Φd，實現“采集”數據環節.

具體流程框架如圖1所示.本文方法的計算效率不僅取決于Landweber迭代次數，還取決于稀疏變換的計算速度，為了保證噪聲壓制效果，當選用seislet變換作為壓縮感知迭代噪聲壓制方法中的稀疏變換基時，計算效率比傳統的迭代去噪方法更低一些，但是可以通過并行計算改善處理效率.

圖1 實現流程框架Fig.1 Processing framework

2 模型測試

首先選取復雜的Sigmoid疊后地震數據模型(Claerbout，2008)測試本文方法，該模型中包括了傾斜地層、斷層和不整合面等構造特征，時間方向的采樣間隔為4 ms，空間方向的采樣間隔為8 m.為了定量分析噪聲壓制的結果，本文選取全局信噪比作為去噪效果的衡量(劉洋等，2017)：

(12)

圖2a是不含噪聲的模型數據，圖2b是加入較強振幅的非平穩隨機噪聲數據，含噪聲數據信噪比為4.4 dB，與噪聲方差為7.5×10-7的平穩噪聲信噪比相當.模型數據中的非平穩噪聲被定義為時間和空間方向上性質不穩定的高斯白噪.將時間方向上平穩的高斯隨機噪聲序列記為n1(t)(圖3a)，對平穩隨機噪聲序列中的部分取值進行放大或改變取值的分布特征可以得到非平穩高斯隨機噪聲(鐘鐵，2016)，則通過平穩的噪聲序列n1(t)可以構造并表示出非平穩的噪聲序列(圖3b)，非平穩噪聲添加在時間方向上，記為n2(t)：

(13)

其中k為噪聲序列中取值放大的開始位置，取值放大的長度為50，放大倍數取為2倍，不同地震道上噪聲放大的起始時間k隨機.將全頻帶的非平穩噪聲進行帶通濾波，僅保留合理頻率范圍內的噪聲，如圖3b所示，構造出的隨機噪聲具有非平穩的時變特征，更接近于真實的隨機噪聲.

圖2 模型數據(a) Sigmoid模型； (b) 含噪聲數據.Fig.2 Synthetic model data(a) Sigmoid model； (b) Noisy data.

圖3 不同平穩性的隨機噪聲(a) 平穩隨機噪聲； (b) 非平穩隨機噪聲.Fig.3 Random noise with different stationarities(a) Stationary random noise； (b) Non-stationary noise.

圖4 不同去噪方法結果對比(a) f-x預測濾波去噪結果； (b) f-x預測濾波壓制的噪聲； (c) 本文迭代方法小波變換下去噪結果(未除偏)； (d) 本文迭代方法小波變換下壓制的噪聲(未除偏)； (e) 本文迭代方法小波變換下去噪結果(除偏)； (f) 本文迭代方法小波變換下壓制的噪聲(除偏).Fig.4 Comparison of denoised results by different methods(a) Result using f-x prediction filtering； (b) The noise removed by f-x prediction filtering； (c) Result using the proposed method with wavelet transform(before debiasing)； (d) The noise removed by the proposed method with wavelet transform(before debiasing)； (e) Result using the proposed method with wavelet transform(after debiasing)； (f) The noise removed by the proposed method with wavelet transform (after debiasing).

為了進一步提升本文方法的有效性，使用seislet變換作為稀疏變換基，圖5a為seislet稀疏變換在本文迭代方法框架下的去噪結果，為了和小波變換進行對比，選取百分比閾值為16%，迭代次數為10，從去噪結果中可以看出本文方法的能夠有效壓制大部分隨機噪聲，尤其有效弱信號損失的恢復能力增強，處理之后的同相軸連續性提高，信噪比提升至11.7 dB(圖5a).在差剖面(圖5b)中可以看出，本文方法壓制的主要為隨機噪聲和少部分的斷層信息.由于seislet變換對有效信號的壓縮能力更強，因此L1正則化項的影響很小，故“修復”步驟被省略，此時除偏修復的幅度主要來自噪聲，會降低處理效果，因此“修復”步驟主要適用于壓縮能力一般的稀疏變換基.與f-x域流式預測濾波方法(圖5c)相比，本文方法計算結果信噪比更高，f-x域流式預測濾波方法處理結果信噪比為8.5 dB，差剖面(圖5d)中可以看出本文方法去掉的噪聲更多且損失有效信號更少，但流式方法的計算量要低于本文方法；f-x域RNA方法(圖5e)處理結果信噪比為10.7 dB，對比可見本文方法的信噪比更高，本文方法差剖面(圖5f)中去掉的有效信息更少，而且沒有頻率域預測濾波方法中常見的虛假同相軸問題.

圖5 本文迭代方法處理結果(a) seislet變換下去噪結果； (b) seislet變換下壓制的噪聲; (c) f-x域流式預測濾波方法去噪結果； (d) f-x域流式預測濾波方法壓制的噪聲； (e) f-x域RNA方法去噪結果； (f) f-x域RNA方法壓制的噪聲.Fig.5 Denoising results using the proposed method(a) The denoised result using the proposed method with seislet transform ； (b) The noise removed by the proposed method with seislet transform； (c) The denoised result using f-x streaming prediction filter； (d) The noise removed by f-x streaming prediction filter； (e) The denoised result using f-x RNA； (f) The noise removed by f-x RNA.

如圖6所示為Sigmoid模型不同方法處理結果的f-k譜.圖6a為不含噪聲模型的f-k譜，圖6b為含噪聲模型的f-k譜，對比各方法f-k譜可見，基于seislet變換的本文方法(圖6f)去噪效果最佳，基于wavelet變換的本文方法(圖6d)去除隨機噪聲不夠徹底，除偏(圖6e)還會引入少量噪聲，而f-x預測濾波方法(圖6c)處理后仍然存在部分高波數的隨機噪聲干擾.f-x域流式預測濾波方法(圖6g)處理后在有效信號的頻帶范圍內殘留噪聲較多，f-x域RNA方法(圖6h)去除噪聲較為徹底，僅在高波數部分殘留少量噪聲.各方法間去噪性能對比如表1所示.基于seislet變換的本文方法去噪結果信噪比(SNR)最高，均方誤差(MSE)最小，雖然迭代次數稍多于基于wavelet變換的本文方法，但其去噪效果好于后者，f-x預測濾波方法信噪比最低，均方誤差最大.一些其他的稀疏變換基，如shearlet(童思友等，2019)、contourlet(Li and Gao，2013)也可以用于本文去噪方法框架中，理論上也會實現較好的效果.

圖6 模型數據f-k譜(a) 原始信號f-k譜； (b) 含噪聲數據f-k譜； (c) f-x預測濾波結果f-k譜； (d) 本文迭代方法小波變換下未除偏結果f-k譜； (e) 本文迭代方法小波變換下除偏結果f-k譜； (f) 本文迭代方法seislet變換下去噪結果f-k譜； (g) f-x域流式預測濾波方法去噪結果f-k譜； (h) f-x域RNA方法去噪結果f-k譜.Fig.6 Comparison of different f-k spectras(a) Original signal； (b) Noisy data； (c) Denoised result using f-x prediction filtering； (d) Denoised result using the proposed method with wavelet before debiasing； (e) Denoised result using the proposed method with wavelet after debiasing； (f) Denoised result using the proposed method with seislet； (g) Denoised result using f-x streaming prediction filter； (h) Denoised result using f-x RNA method.

表1 各方法去噪性能對比Table 1 Comparison of denoising performance in different methods

3 實際數據處理

3.1 實際疊后數據處理

接下來，利用某地區實際疊后時間偏移數據對本文方法進行測試(圖7a)，該數據在時間方向上為700個采樣點，采樣間隔為1 ms，空間方向310個地震道.該數據淺層存在近平行的同相軸，在深層存在比較復雜的具有不同傾角的傾斜同相軸，數據中的較強隨機噪聲嚴重影響著復雜同相軸的連續性，進而影響高精度的地層解釋.圖7b為f-x預測濾波方法的去噪結果，噪聲水平得到一定程度的衰減，同時淺層同相軸連續性增強.將濾波后的結果與原始數據相減得到濾除的噪聲剖面，如圖7c所示.在噪聲剖面中可以看到標準f-x預測濾波方法對有效波振幅的保護效果差，即使淺層比較平穩的近水平同相軸也被部分衰減.圖7d為使用seislet變換作為稀疏變換的壓縮感知迭代噪聲壓制方法的去噪結果，迭代次數為10，百分比閾值選取為18%.從圖7b和圖7d對比可以看到，本文方法的處理結果優于f-x預測濾波方法，處理后數據同相軸清晰，深層復雜構造得到更好的保護，而f-x預測濾波的處理結果中同相軸被一定程度地模糊化，有效信號損失較大.對比兩種方法的差剖面(圖7c和圖7e)，本文方法去除的主要是隨機噪聲，很難直觀看到有效波的損失，而f-x預測濾波損失了大量有效波同相軸振幅.最后，通過處理結果的f-k譜來進一步對比兩種方法的處理效果，如圖8所示.圖8c表明本文能夠更好地壓制原始數據低頻和高波數范圍內(圖8a)的隨機噪聲干擾，而f-x預測濾波方法處理后仍然存在部分高波數的隨機噪聲干擾(圖8b).

圖7 實際數據處理結果(a) 實際數據； (b) f-x預測濾波方法的去噪結果； (c) f-x預測濾波方法壓制的噪聲； (d) 本文方法的去噪結果； (e) 本文方法壓制的噪聲.Fig.7 Processing result of real data(a) Real data; (b) Denoised result using f-x prediction filtering; (c) Noise removed by f-x prediction filtering； (d) Denoised result using the proposed method with seislet transform； (e) Noise removed by the proposed method with seislet transform.

圖8 處理結果f-k譜對比(a) 疊后實際數據f-k譜； (b) f-x預測濾波方法的去噪結果f-k譜； (c) 本文方法的去噪結果f-k譜.Fig.8 Comparison of f-k spectra for processing results(a) Poststack real data； (b) Denoised result using f-x prediction filtering； (c) Denoised result using the proposed method.

3.2 實際疊前數據處理

利用某地區陸上實際疊前數據對本文方法進行進一步的測試，圖9a為去面波后的CMP道集，可以看到該數據信噪比較低，強隨機噪聲的存在嚴重影響著該數據中的反射同相軸的連續性，有效信號的能量較弱.本文方法可以直接擴展為三維方式，三維處理方法能夠進一步提高效果，但是將會增加計算負擔，所以將三維數據按照二維切片進行獨立處理.圖9b為f-x預測濾波方法的去噪結果，可見由于強噪聲影響，濾波后的結果中也幾乎不可見同相軸的存在，數據信噪比仍然較低.圖9c為使用VD-seislet變換作為稀疏變換的壓縮感知迭代噪聲壓制方法的去噪結果，其中seislet變換的傾角模式由均方根速度計算(Liu et al.，2015)，迭代次數為11，百分比閾值選取為10%.從圖9c中可以看到，經本文方法處理后，結果中同相軸連續性明顯增強，同相軸形態恢復良好，數據信噪比得到提升.對比兩種方法處理后的疊加剖面(圖10a和圖10b)，本文方法的疊加剖面中同相軸連續性明顯增強，尤其中淺層結構恢復較為明顯，f-x預測濾波方法的疊加剖面中同相軸能量較弱，連續性較差，隨機噪聲的存在仍嚴重影響著該剖面的信噪比.

圖9 疊前實際數據處理結果(a) 疊前實際數據； (b) f-x預測濾波方法的去噪結果； (c) 本文方法的去噪結果.Fig.9 Processing results of pre-stack real data(a) Pre-stack real data; (b) Denoised result using f-x prediction filtering; (c) Denoised result using the proposed method with seislet transform.

圖10 不同方法處理結果疊加剖面(a) f-x預測濾波方法去噪結果疊加剖面； (b) 本文方法的去噪結果疊加剖面.Fig.10 Stacked profiles using different methods(a) f-x prediction filtering； (b) The proposed method with seislet transform.

4 結論與展望

本文系統地論述了壓縮感知框架下的隨機噪聲壓制問題，建立了壓縮感知、軟閾值去噪和基追蹤降噪的理論關聯，提出了基于迭代軟閾值算法的“采集-重建-修復”技術方案，為提升本文方法的時空變信號保護及非平穩隨機噪聲壓制能力，選用對地震信號具有更高壓縮能力的seislet變換作為稀疏變換，同時利用匹配的除偏方法，對最優化問題中由于較差壓縮能力的稀疏變換在L1正則化條件中所引起的振幅衰減實現有效的恢復.通過與工業標準的f-x預測濾波方法進行對比，本文方法能夠更好地處理同相軸傾角的時空變化，避免頻率域預測濾波類方法的虛假同相軸問題，能夠壓制強隨機噪聲且更好地平衡時空變有效信號保護和隨機噪聲壓制的關系，同時能夠適應隨機噪聲的非平穩分布特征，提供一種有效的隨機噪聲壓制方法.理論模型和實際數據測試結果表明本文提出方法能夠在保護時空變有效信息的前提下壓制較強幅度的非平穩隨機噪聲，提高同相軸連續性，為后續高精度的解釋提供數據基礎.為提升本文方法的性能，仍需要進一步研究，如尋求對數據更優的稀疏表示方法，例如Shearlet、Contourlet及超完備冗余字典可能會提供更加優秀的數據表征.

致謝感謝同濟大學王本鋒研究員的建議和討論，感謝兩位匿名審稿專家提出的寶貴意見.