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

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

0 引言

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

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

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

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

1 基本理論

1.1 壓縮感知基本理論回顧

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

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

y=Φx，

(1)

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

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

y=ΦΨc=Θc，

(2)

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

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

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

(3)

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

1.2 壓縮感知框架下的地震數(shù)據(jù)隨機(jī)噪聲壓制問題

信號中含有隨機(jī)噪聲干擾時，公式(1)變?yōu)?/p>

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

(4)

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

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

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

(5)

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

(6)

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

(7)

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

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

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

(8)

其中：

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

(9)

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

(10)

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

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

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

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

(11)

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

1.4 壓縮感知迭代噪聲壓制方法的實現(xiàn)流程

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

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

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

圖1 實現(xiàn)流程框架Fig.1 Processing framework

2 模型測試

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

(12)

圖2a是不含噪聲的模型數(shù)據(jù)，圖2b是加入較強(qiáng)振幅的非平穩(wěn)隨機(jī)噪聲數(shù)據(jù)，含噪聲數(shù)據(jù)信噪比為4.4 dB，與噪聲方差為7.5×10-7的平穩(wěn)噪聲信噪比相當(dāng).模型數(shù)據(jù)中的非平穩(wěn)噪聲被定義為時間和空間方向上性質(zhì)不穩(wěn)定的高斯白噪.將時間方向上平穩(wěn)的高斯隨機(jī)噪聲序列記為n1(t)(圖3a)，對平穩(wěn)隨機(jī)噪聲序列中的部分取值進(jìn)行放大或改變?nèi)≈档姆植继卣骺梢缘玫椒瞧椒€(wěn)高斯隨機(jī)噪聲(鐘鐵，2016)，則通過平穩(wěn)的噪聲序列n1(t)可以構(gòu)造并表示出非平穩(wěn)的噪聲序列(圖3b)，非平穩(wěn)噪聲添加在時間方向上，記為n2(t)：

(13)

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

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

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

圖4 不同去噪方法結(jié)果對比(a) f-x預(yù)測濾波去噪結(jié)果； (b) f-x預(yù)測濾波壓制的噪聲； (c) 本文迭代方法小波變換下去噪結(jié)果(未除偏)； (d) 本文迭代方法小波變換下壓制的噪聲(未除偏)； (e) 本文迭代方法小波變換下去噪結(jié)果(除偏)； (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).

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

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

圖6 模型數(shù)據(jù)f-k譜(a) 原始信號f-k譜； (b) 含噪聲數(shù)據(jù)f-k譜； (c) f-x預(yù)測濾波結(jié)果f-k譜； (d) 本文迭代方法小波變換下未除偏結(jié)果f-k譜； (e) 本文迭代方法小波變換下除偏結(jié)果f-k譜； (f) 本文迭代方法seislet變換下去噪結(jié)果f-k譜； (g) f-x域流式預(yù)測濾波方法去噪結(jié)果f-k譜； (h) f-x域RNA方法去噪結(jié)果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 實際數(shù)據(jù)處理

3.1 實際疊后數(shù)據(jù)處理

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

圖7 實際數(shù)據(jù)處理結(jié)果(a) 實際數(shù)據(jù)； (b) f-x預(yù)測濾波方法的去噪結(jié)果； (c) f-x預(yù)測濾波方法壓制的噪聲； (d) 本文方法的去噪結(jié)果； (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 處理結(jié)果f-k譜對比(a) 疊后實際數(shù)據(jù)f-k譜； (b) f-x預(yù)測濾波方法的去噪結(jié)果f-k譜； (c) 本文方法的去噪結(jié)果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 實際疊前數(shù)據(jù)處理

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

圖9 疊前實際數(shù)據(jù)處理結(jié)果(a) 疊前實際數(shù)據(jù)； (b) f-x預(yù)測濾波方法的去噪結(jié)果； (c) 本文方法的去噪結(jié)果.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 不同方法處理結(jié)果疊加剖面(a) f-x預(yù)測濾波方法去噪結(jié)果疊加剖面； (b) 本文方法的去噪結(jié)果疊加剖面.Fig.10 Stacked profiles using different methods(a) f-x prediction filtering； (b) The proposed method with seislet transform.

4 結(jié)論與展望

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

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