基于SSEC-EWT的地震資料噪聲壓制算法

孟 娟 高 琴 李亞南②

(①防災科技學院電子科學與控制工程學院，河北三河 065201；②河北省地震災害監(jiān)測儀器與監(jiān)測技術(shù)重點實驗室，河北三河 065201)

0 引言

地震勘探是油氣勘探開發(fā)的重要手段，且依賴于高質(zhì)量的處理成果數(shù)據(jù)。受人類活動、環(huán)境、天氣、儀器等多方面因素影響，實際采集的地震數(shù)據(jù)中往往存在各種噪聲干擾，這無疑會影響地震資料處理和解釋的準確性。為了提高地震勘探精度，使地震資料能更真實地反映地下地質(zhì)情況，壓制噪聲干擾在地震數(shù)據(jù)處理中就顯得極為重要了。

地震資料噪聲包括規(guī)則噪聲和隨機噪聲兩大類。面波是其中一種能量較強且廣泛存在的規(guī)則干擾，其頻率低、衰減慢，大大降低了地震資料信噪比；而隨機噪聲則顯得雜亂無規(guī)律，頻率分布寬泛，空間上無確定的視速度。因此，可根據(jù)面波和隨機噪聲各自的特點探尋相應的適用去噪方法。

在面波壓制領(lǐng)域，S變換[1]、廣義S變換[2-3]、Curvelet變換[4-5]、小波變換[6-7]和f-x域濾波[8]等方法廣泛使用，利用面波與有效波在變換域上的特征差異實施分離。這些方法有一定效果但也有局限性，且會不同程度地損害有效信號。

在隨機噪聲去除方面，經(jīng)驗模態(tài)分解(Empirical mode decomposition，EMD)[9]、小波閾值去噪[10]、深度學習[11]等方法被引入，但這些方法同樣有各自優(yōu)缺點。如EMD存在模態(tài)混疊、偽模態(tài)等問題，難以兼顧增強去噪效果與降低有效信號損傷；小波閾值去噪效果依賴于小波基的選?。簧疃葘W習需大量訓練數(shù)據(jù)，計算量大。這些算法一般只針對單一類型噪聲，難以同時處理面波和隨機噪聲這兩種最普遍存在的干擾。

作為一種自適應的信號處理方法，EWT(經(jīng)驗小波變換)[12]結(jié)合小波變換和EMD技術(shù)優(yōu)勢，通過分割信號的傅里葉譜，并在每個分割區(qū)間構(gòu)造正交小波基，將信號分解為若干分量，從而準確提取信號中固有模態(tài)分量，適用于非線性非平穩(wěn)信號。

EWT推出后被廣泛應用于噪聲壓制及故障診斷等方面[13-15]，在地震數(shù)據(jù)處理與分析領(lǐng)域，也倍受青睞。Liu等[16]基于EWT對多通道地震資料進行時頻分析，分解出不同頻率分量，并根據(jù)計算邊界提取地震信號分量，分辨率明顯高于傳統(tǒng)連續(xù)小波變換?；谟蜌獠貙е碌卣鸩ㄋp增大的原理，Hui等[17]利用EWT分析地震波衰減時頻特性，圏定衰減異常位置，以厘清含氣層厚度變化，實現(xiàn)油氣檢測。在地震數(shù)據(jù)去噪方面，覃發(fā)兵等[18]利用EWT將地震信號分解為若干分量，并計算各分量主頻，剔除閾值外分量后重構(gòu)實現(xiàn)隨機噪聲壓制。Chen等[19]利用EWT將地震信號分解為若干IMF，再利用聚類算法對含噪分量做閾值化處理。

以上應用表明利用EWT可有效分析地震信號，但傳統(tǒng)EWT根據(jù)信號頻譜幅值進行頻帶劃分和信號分解，易受噪聲影響，因此有必要根據(jù)地震信號特點改進EWT頻譜分割方法，使之更適用于地震信號。Liu等[20]基于尺度空間表示提取信號傅里葉譜的慢變分量，再根據(jù)尺度參數(shù)獲取地震信號中包含的頻率分量和邊界，自適應分割頻譜，但尺度參數(shù)通常依據(jù)經(jīng)驗獲取，使頻譜分割的準確性受限。孟娟等[21]對地震信號連續(xù)小波變換后計算各頻率點能量，得到小波譜能量曲線(WSEC)，以該曲線局部極大值點的頻率及ε鄰域法確定頻譜分割邊界進行改進EWT，實現(xiàn)依地震信號頻率和能量的自適應分解，取得良好效果，但連續(xù)小波變換中基小波函數(shù)的選取影響變換結(jié)果精度，且通過尺度參數(shù)調(diào)整小波函數(shù)形狀，獲得的并不是真實的時頻譜。

為同時分離地震資料中的面波和隨機噪聲，針對傳統(tǒng)EWT根據(jù)頻譜極值分割頻帶抗干擾性差的不足，本文提出基于S譜能量曲線的改進EWT地震資料噪聲壓制算法，實現(xiàn)面波和隨機噪聲同步壓制。通過對地震信號S變換后計算各頻點能量得到S譜能量曲線(SSEC)，根據(jù)曲線極大值點分布進行頻帶分割，將原信號自適應分解為若干按頻率和能量分布的IMF，并根據(jù)SSEC確定面波IMF；為保護面波IMF中的低頻有效信號，對其進行帶通濾波，以實現(xiàn)精準的面波剔除；最后根據(jù)有效波頻率閾值和各IMF主頻去除隨機噪聲，完成面波和隨機噪聲的壓制，有效提高地震資料信噪比。

1 EWT基本原理

EWT通過計算信號的傅里葉頻譜，并基于頻譜幅值的極值點進行頻譜分割，在每個分割區(qū)間內(nèi)構(gòu)建正交的小波濾波器組，提取不同頻帶內(nèi)的調(diào)幅、調(diào)頻分量，實現(xiàn)信號分解。其原理詳述如下：

(1)獲取原信號f(t)的傅里葉頻譜F(ω)，ω∈[0，π]。

(2)獲取F(ω)的M個按頻率升序排列的局部極大值點ω1，…，ωM，將[0，π]頻帶劃分成N(N≤M)個子頻帶，設以ωn為中心的分量子頻帶為[Ωn-1,Ωn]。獲取邊界Ω的常用方法有兩種：一是根據(jù)頻譜的M個極大值點，將相鄰兩個極大值中點定義為邊界，Ωn-1=(ωn-1+ωn)/2，n=2，…,N，Ω0=0,ΩN=π；或先找到信號頻譜的極大值點，再尋找兩相鄰極大值間的極小值點，以該極小值點為邊界。這樣將原信號頻譜劃分為N個連續(xù)、不重疊的子頻帶。

(1)

(2)

(3)

(4)

則得到的低頻分量和高頻分量分別為

(5)

(6)

須注意的是，傳統(tǒng)EWT分解信號需預先設定經(jīng)驗模態(tài)分解階數(shù)N，該參數(shù)選取會影響信號分解精度，過大或過小都不能較好表征原信號。

2 基于S譜能量曲線的改進EWT

2.1 S譜能量曲線SSEC

作為短時傅里葉變換和以Morlet小波為小波基的連續(xù)小波變換的發(fā)展，S變換采用帶有頻率變量的高斯窗函數(shù)截取信號，從而實現(xiàn)信號的局部分析。S變換免去了窗函數(shù)選擇，且其窗函數(shù)可隨頻率需求自適應地變化，具有多分辨率的特點。與連續(xù)小波變換相比，S變換分辨率依賴于頻率而不是尺度，因此更直觀；且每一頻率的絕對相位不變，與傅里葉譜直接關(guān)聯(lián)，能精確描述信號各時刻的頻譜，較好適應非平穩(wěn)信號頻率不規(guī)律變化的特點，因此在地震數(shù)據(jù)處理領(lǐng)域應用較廣。其原理如下：

根據(jù)短時傅里葉變換，定義高斯窗函數(shù)為

(7)

式中σ為尺度因子，可控制窗函數(shù)寬度。它與頻率直接聯(lián)系，從而改變時頻分辨率

(8)

代入短時傅里葉變換公式，即有

(9)

可得到其S變換

(10)

對信號進行S變換得到信號時頻分布的S譜后，根據(jù)式(11)對頻點i定義S譜能量，實現(xiàn)信號時頻域到頻率—能量域的轉(zhuǎn)換。表示為

(11)

式中：fi為S變換后頻點i的頻率；S(fi,τ)為該頻點對應的S變換值，即S譜。

對某信號按式(11)獲取所有頻點能量值，從而得到依頻率分布的SSEC，該曲線能清晰描繪不同頻率處信號的能量分布及不同頻率成分的能量強弱。

2.2 基于SSEC的改進EWT

對含面波的地震信號進行S變換并獲取SSEC，SSEC能直接反映不同頻率成分的能量大小。由于面波的頻率低且能量強，因此通過SSEC能區(qū)分面波與有效地震波。不同于傳統(tǒng)EWT根據(jù)傅里葉頻譜的局部極大值分割頻譜，本文根據(jù)SSEC的極大值點自適應地分割頻譜，即以SSEC極大值點對應的頻率為初始邊界，通過ε鄰域法搜索重新確定邊界，以進行改進EWT，根據(jù)信號不同頻率分量的能量大小自適應分解。算法如下：

(1)對地震信號s(t)進行傅里葉變換和S變換，得到傅里葉頻譜F(ω),ω∈[0,π]和S譜S(f,τ)。

(2)根據(jù)式(11)獲取各頻率點能量，繪制SSEC。

改進EWT可將原地震信號分解為一系列具有緊支撐頻譜的IMF分量，對各IMF解調(diào)可得到其瞬時頻率、幅度、相位等屬性，從而獲取原信號特征。

3 基于SSEC-EWT的噪聲壓制方法

與地震信號相比，面波頻率范圍低，頻帶一般為4～18Hz[2]，且能量較強。本文改進EWT是根據(jù)S變換后各頻點能量大小進行頻譜分割，各IMF是根據(jù)原信號不同頻率成分的能量大小分解的，從而可將低頻、高能量面波成分從原信號中分離出來。

而隨機噪聲在時域內(nèi)廣泛存在，且一般為高頻率，可根據(jù)有效地震信號頻率范圍設定閾值，剔除閾值范圍外的噪聲IMF，實現(xiàn)隨機噪聲壓制?；诖怂悸?，本文提出基于SSEC-EWT的面波和隨機噪聲壓制方法，其流程如圖1所示。

圖1 基于SSEC-EWT的面波和隨機噪聲壓制方法流程

3.1 面波壓制

根據(jù)有效波與面波在頻率和能量上的不同，本文利用改進EWT進行地震信號分解，提取面波IMF后進行面波壓制，流程如下：

(1)計算地震信號SSEC，以能量曲線極值點及ε鄰域法分割頻譜，進行改進的EWT，自適應得到依頻率和能量分解的各模態(tài)分量IMF1，…，IMFN。

(2)對各IMF進行希爾伯特變換，得到瞬時頻率、瞬時振幅，并據(jù)下式計算各IMF主頻

(12)

式中：fm為IMF主頻；f為瞬時頻率；s(f)為瞬時振幅譜。

(3)結(jié)合SSEC，找到頻率在4～18Hz范圍內(nèi)，且能量強于其他鄰近子頻帶的子頻段位置，確定面波所在的IMF。

(4)根據(jù)SSEC中面波頻率峰值，對該IMF進行帶通濾波。濾波器截止頻率為面波IMF頻譜峰值的6dB帶寬頻率，即以頻譜功率下降為功率峰值25%、幅度下降為峰值50%時的頻率為截止頻率，對面波IMF進行帶通濾波，以最大限度地保護面波IMF中的低頻有效波不被破壞。

(5)從原記錄剔除濾波后面波IMF，實現(xiàn)面波壓制。

3.2 隨機噪聲壓制

隨機噪聲一般以高頻率廣泛存在于地震資料中，常規(guī)地震勘探中檢波器接收的有效地震波頻率范圍一般在150Hz以內(nèi)[22]。因此，根據(jù)地震波與隨機噪聲的頻率范圍差異，可實現(xiàn)隨機噪聲衰減。本文基于改進EWT的隨機噪聲壓制具體算法如下：

(1)設定有效地震波最高頻率閾值為150Hz；

(2)計算完成面波壓制后的其余IMF主頻，剔除閾值范圍外IMF，進行隨機噪聲壓制；

(3)對剔除面波后IMF和隨機噪聲IMF 信號進行重構(gòu)，得到噪聲壓制后的地震信號。

至此，完成面波和隨機噪聲的同時壓制。對二維資料中所有道進行以上相同操作，即可得到最終去噪記錄。

4 實驗與分析

4.1 SSEC-EWT地震信號分解

先對模擬地震記錄進行SSEC-EWT分解，并與傳統(tǒng)的EWT及文獻[21]中的小波譜能量曲線EWT(WSEC-EWT)進行地震信號分解對比。利用Ricker子波模擬單道地震信號，設0.5s處的面波是主頻為10Hz的Ricker子波，幅值為5，能量較強，0.7、0.8、0.9s處為有效波，頻率分別為20、35、50Hz，幅值分別為2、2、3，其合成波及SNR=3dB的含噪合成記錄如圖2所示。

對該含噪記錄(圖2b)進行S變換，根據(jù)S譜對每個頻點按式(11)計算能量得到SSEC(圖3a)。可見SSEC和WSEC(圖3b)都能較好體現(xiàn)原信號不同頻率處的能量大小。但因CWT是通過尺度參數(shù)改變小波函數(shù)形狀，變換結(jié)果并不是真正的時頻譜，從而導致信號特征精度受影響，而S變換反映的是信號最真實的頻率構(gòu)成，得到的是真實時頻譜。因此，相比于WSEC，本文SSEC能更真實地體現(xiàn)原信號低頻和高頻部分各頻率成分的能量分布。

圖2 合成單道地震波(a)及含噪地震波(b)

圖3 圖2b數(shù)據(jù)的SSEC(a)和WSEC(b)

基于SSEC極大值點及ε鄰域法進行子頻帶劃分結(jié)果如圖4a所示，可見較好地將面波頻譜與鄰近頻率有效波頻譜分割開；WSEC-EWT分割結(jié)果(圖4b)中面波與有效波頻譜有少許混疊；而傳統(tǒng)EWT劃分的子帶(圖4c)更不能很好地反映面波與有效波的差異，故不能直接應用于面波壓制。相比WSEC，本文SSEC極大值點能更直接、準確地表示有效信號和面波的頻率和能量大小，因此本文頻譜分割方法能更顯著體現(xiàn)地震信號中面波與有效波的能量、頻率的不同，具有更好的頻譜分割性能。

圖4 SSEC-EWT(a)、WSEC-EWT(b)、傳統(tǒng)EWT(c)三種方法的頻譜分割結(jié)果

本文方法得到的IMF如圖5a所示，可見SSEC-EWT和WSEC-EWT(圖5b)都能較好地根據(jù)信號的頻率和能量進行模態(tài)分解。圖5a中IMF1主要為0.5s處的面波分量(含有微量有效波)，IMF2主要為0.7、0.8、0.9s處的有效波分量，而IMF4主要為高頻隨機噪聲，有利于實現(xiàn)面波與隨機噪聲的有效分離；圖5b中IMF1主要為面波，但其中包含部分有效波；而傳統(tǒng)EWT分解結(jié)果(圖5c)中，有效波與面波混疊不清，難以實現(xiàn)面波提取。因此，本文SSEC能完整描述原信號各頻率點的能量特征，其曲線極大值點能更直接、更真實地表示有效信號和面波的頻率和能量值，從而更好地進行頻譜分割和信號分解，實現(xiàn)面波與有效波的精確分離。

圖5 SSEC-EWT(a)、WSEC-EWT(b)、傳統(tǒng)EWT(c)三種分解方法所得的IMF

為防止面波IMF中有效波被同時壓制，根據(jù)SSEC面波頻率峰值，對頻段做帶通濾波，保留面波頻率附近的有效波。圖6中濾波前的面波IMF中含有少量有效波，濾波后大部分有效波被保留，二者與理想純面波的標準差分別為0.0225、0.0209，方差分別為0.1310、0.1154，可見能最大限度減少對有效波的損傷。從經(jīng)面波和隨機噪聲提取后重構(gòu)所得的最終去噪結(jié)果(圖7)可見，面波壓制較徹底，隨機噪聲也得到有效抑制，去噪后的SNR=8.5dB。

圖6 濾波前、后的面波IMF

圖7 理想合成波(a)及本文方法去噪結(jié)果(b)

4.2 SSEC-EWT抗噪性測試

為測試SSEC-EWT的抗噪性，定義

(13)

首先測試在SNR較低(圖8)情況下SSEC-EWT分解的準確性。選取有效波頻率分別為16、35、48Hz，面波頻率為6Hz，添加隨機噪聲后SNR=-8.8dB，有效波淹沒于噪聲中難以分辨。

圖8 模擬合成波(a)及其加噪信號(SNR=-8.8dB)(b)

從本文SSEC-EWT分解結(jié)果(圖9a)可見，IMF1是面波IMF，其中含有很少量有效波，IMF2、IMF3主要為有效波，IMF4、IMF5主要為高頻噪聲。顯然SSEC-EWT根據(jù)原信號的頻率成分和能量大小進行信號分解，在SNR較低時仍具有較好的信號分解能力，具有較好的抗噪性。面波和隨機噪聲壓制后所得結(jié)果(圖9b)中，SNR=9.3dB，面波壓制較徹底，隨機噪聲也得到有效壓制，且有效波得到較好保護。可見在較低SNR下，本文算法仍具有較好性能。

圖9 SNR較低時SSEC-EWT分解所得IMF(a)及去噪結(jié)果(b)

再針對不同SNR(-15～15dB)模擬信號，用本文算法進行處理(圖10)并做對比，可見經(jīng)本文算法處理后，SNR明顯提高(幅度高達17.9dB)，尤其是在低SNR時，本文算法也能較好壓制面波和隨機噪聲，具有較好抗噪性能。

圖10 相同SNR條件下本文去噪效果

進一步測試有效波與面波的頻率有重疊時，改進EWT的分解效果。圖11a中面波主頻設為6Hz(圖11a上)，最大峰值為5，有效波主頻為7、15、20Hz，最大峰值均為3，并添加隨機噪聲，面波與有效波的頻率幾近重疊(圖11a下)。按式(11)得到的SSEC如圖11b所示，該曲線較清晰地體現(xiàn)了不同頻率地震波成分的能量大小，有助于區(qū)分有效波與面波。從基于該SSEC曲線進行頻譜分割后所得IMF(圖11c)可見，IMF1主要為面波，IMF4為隨機噪聲，由于面波與有效波頻率非常接近，故面波IMF中尚含有少量的有效波。

圖11 面波和有效波幾近重疊的模擬地震波(a)及其SSEC(b)與分解所得IMF(c)

對IMF1進行帶通濾波(圖12a)后，除了極少部分有效波被損傷，絕大部分有效波得到有效保護。從最終去噪結(jié)果(圖12b下)也可見，面波得到充分壓制，有效波得到盡可能保護，顯著提高了資料的SNR。

圖12 面波IMF濾波前、后對比(a)及最終去噪結(jié)果(b)

4.3 合成二維地震道集

對同時含有面波和隨機噪聲的二維地震道集進行噪聲壓制測試。根據(jù)層狀速度模型和波場模擬得到圖13a所示地震道集。面波速度為200m/s，直達波速度為2000m/s，合成地震記錄采樣間隔為1ms，共計10道、1000個采樣點，面波主頻為8Hz，反射波頻率為40Hz。添加信噪比為15dB的隨機噪聲(圖13b)，面波與有效反射波疊加，導致同相軸模糊。

從應用本文算法壓制面波和隨機噪聲后結(jié)果道集(圖13c)可見，面波和隨機噪聲同時得到有效壓制，有效地震信號更光滑、更平穩(wěn)，且有效反射波與面波混疊的同相軸變得清晰。觀察最終被去除的噪聲殘差(圖13d)可知，其上基本為面波和隨機噪聲，僅損失極少量的有效波，且面波壓制徹底。

圖13 含面波(a)、且含隨機噪聲(b)地震道集、本文方法去噪結(jié)果(c)及去噪殘差(d)

對經(jīng)SSEC-EWT分解后所得面波IMF記錄(圖14a)進行濾波處理，可見濾波后記錄(圖14b)上大部分有效波被保留，僅極少部分有用信息被濾除。

圖14 本文方法濾波前(a)、后(b)的面波IMF記錄

將WSEC-EWT與本文SSEC-EWT提取面波的結(jié)果(圖15)進行對比，可見二者都能較好地提取面波，但相比WSEC-EWT中面波IMF會包含“少量”有效波，而本文方法提取面波則只含有“微量”有效波，面波提取更精準。這是因為相比連續(xù)小波變換，S變換的時頻分析精確度更高，而基于頻率和能量進行頻譜分割后能更好分解信號。

圖15 WSEC-EWT(a)與本文算法(b)提取面波對比

為對比隨機噪聲壓制效果，對圖13b地震道集應用文獻[18]算法去噪(圖16a)，可見雖然隨機噪聲被壓制，但面波未被壓制；而本文算法去噪后(圖13c)不僅能去除面波，還能有效壓制隨機噪聲，而且能最大限度保護有效波。從文獻[18]算法去噪殘差(圖16a)可知，部分有效波被損害；而本文算法去噪殘差(圖13d)中基本不含有效波。

圖16 文獻[18]去噪結(jié)果(a)及隨機噪聲殘差(b)

不同于傳統(tǒng)EWT按傅里葉頻譜極值點分割頻譜，本文算法是根據(jù)頻點能量進行頻譜分割，故可使面波、有效波、隨機噪聲根據(jù)各自能量大小自適應分割，進而能同時分離面波和隨機噪聲。

4.4 實際地震資料處理

圖17所示實際地震數(shù)據(jù)共有95道，每道采樣點數(shù)為500，采樣間隔為2ms，圖上可見典型的掃帚狀面波，降低了地震資料的SNR。

圖17 實際地震資料

觀察從該記錄抽取第17道信號進行S變換后根據(jù)SSEC所得子頻帶分割(圖18a)情況，可見低頻有效波和面波的頻譜被分割為不同子頻帶。在根據(jù)頻譜分割得到的IMF(圖18b)中，IMF2為面波，IMF1為低頻有效波，IMF3～IMF4為有效波，IMF5為隨機噪聲。對該記錄中每道都進行改進EWT，并根據(jù)各IMF的頻率和能量確定面波IMF，再做帶通濾波，分離出記錄中的面波。

圖18 第17道地震信號的子頻帶分割結(jié)果(a)與所得IMF(b)

對比本文方法與WSEC-EWT算法分離出的面波記錄，可見二者都能夠較好地分離出面波，但與WSEC-EWT(圖19b)相比，本文算法提取的面波記錄(圖19a)中僅含更少量的有效波成分。這是因為本文方法能更好地根據(jù)地震信號中各頻率成分的能量大小做信號分解，從而能更精準地提取面波記錄。

圖19 本文(a)與WSEC-EWT(b)算法分離出的面波記錄

本文方法提取的隨機噪聲記錄(圖20a)基本為雜亂無章的隨機噪聲，未見有效波殘留。從最終去除面波和隨機噪聲后的結(jié)果(圖20b)可見，有效波充分保留，面波和隨機噪聲被分離，提高了SNR。

圖20 本文分離出的隨機噪聲(a)與最終去噪結(jié)果(b)

實際地震資料處理結(jié)果表明，SSEC-EWT能根據(jù)地震信號的頻率、能量準確地分解信號，從而實現(xiàn)有效波、面波、隨機噪聲的分離，并最大限度地保護有效波，提高資料的SNR。

5 結(jié)論

(1)基于S變換的S譜能量曲線，能根據(jù)面波、有效波和隨機噪聲的頻率和能量進行頻譜分割和自適應分解，從而實現(xiàn)面波和隨機噪聲的同步壓制。

(2)為減少面波壓制過程中對低頻有效波的損傷，根據(jù)SSEC中的面波頻譜峰值頻率進行帶通濾波，從而達到有效提取面波的同時，最大限度地保留有效反射波。

(3)仿真實驗表明，SSEC能真實、完整地反映信號的能量、頻率特征，完成信號分解，相比WSEC-EWT，本文SSEC-EWT能更精準提取面波。

(4)對模擬地震道集的處理表明，本文方法對面波的壓制效果優(yōu)于WSEC-EWT面波壓制算法，能更充分保護有效波；對隨機噪聲壓制效果也優(yōu)于文獻[18]，表明所提算法的良好性能。

(5)實際地震資料處理結(jié)果表明，該算法不僅能有效抑制高頻隨機噪聲，而且能準確有效地分離面波，提高地震資料的信噪比，是一種能同時去除面波和隨機噪聲的可靠、實用方法。