基于CEEMDAN的時頻峰值濾波微地震隨機噪聲壓制方法

陳毅軍 程 浩 鞏恩普 薛 林

(東北大學深部金屬礦山安全開采教育部重點實驗室，遼寧沈陽 110819)

0 引言

世界油氣需求量的持續增長與常規油氣產量的不斷下降使得非常規油氣成為近年全球油氣資源勘探開發新亮點[1]。在油氣開采過程中，采出、注水、注氣、水力壓裂等作業都會誘發地震。目前很多油田采用壓裂方法以提高油氣采收率，該方法同樣適用于低滲透性的非常規油氣田[2]。為了促進油氣藏開采的高產、高效，可利用現今廣泛應用的“微地震監測”技術研究壓裂誘發的微地震效應，繪制壓裂裂縫空間圖像獲取裂縫產狀等相關信息[3-4]。

通常由壓裂導致的微地震事件產生的地震波較弱，且壓裂施工現場噪聲較強，因此微地震數據一般呈現“弱信號、強干擾”特征，影響初至拾取、發震時刻計算、震源定位等分析。微地震數據中有效信號能量極其微弱、噪聲類型復雜，其信噪比遠低于常規地震記錄，如何壓制其中噪聲是微地震數據處理的關鍵[5]。經過去噪處理后，微地震數據的信噪比、分辨率及保真度等都能得到不同程度的提高，有利于后續的地震資料綜合解釋，對油氣及其他礦藏的開采具有重要的指導意義。

微地震數據去噪方法很多，根據微地震監測方式可分為井中和地面兩種微地震數據去噪。前者的信噪比相對較高，常用去噪方法為偏振濾波、F-K變換[6]等； 后者的信噪比較低，現行去噪方法主要包括基于小波變換、曲波變換、剪切波變換等變換域去噪方法、單道奇異值分解(SVD)及改進的K-L變換法。其中： 傳統傅里葉變換對非平穩信號去噪效果不理想[7]； 變換域方法能很好地區分非平穩信號的突變部分[8]，但去噪結果受基函數影響，且會出現吉布斯效應； 而SVD及K-L變換限制條件多，僅適用于水平及傾斜同相軸。由于地面微地震數據具有噪聲能量強、有效信號能量弱、非平穩等特點，很多常規去噪方法難以獲得較好去噪效果，因此亟待研發地面微地震數據去噪方法。

Huang等[9-11]的經驗模態分解(EMD)方法是一種信號自適應分解方法，無需先驗信息就能突顯系統的物理特性。CEEMDAN(Complete ensemble EMD with adaptive noise)[12]方法是對EMD方法的改進，利用白噪聲均勻分布的統計特性，通過對目標信號多次添加一定幅值的白噪聲來克服模態混疊的影響。但基于EMD的閾值去噪方法對于低信噪比地震記錄濾波效果不佳[13]。即使在信噪比為-9dB情況下，時頻峰值濾波(TFPF)算法也可有效壓制噪聲[14]。這兩種方法現廣泛應用于機器故障監測、信號處理、地震探測等領域。

武安緒等[15]通過EMD對實際地震波形信號進行時頻分析。吳琛等[16]利用EMD方法提取地震信號動力特征。方江雄等[17]提出基于VMD的地震隨機噪聲壓制方法，提高了計算效率和保幅性能。胡瑞卿等[18]結合CEEMDAN和主成分分析法檢測低信噪比微地震初至信號。金雷等[19-20]利用TFPF消除合成地震數據隨機噪聲。林紅波等[21-22]指出TFPF時窗長度隨時間變化能在壓制噪聲的同時更好地保護有效信號。李月等[23-24]系統地分析了TFPF時窗長度和窗型對壓制隨機噪聲的影響，并給出TFPF時窗長度的選取規則。秦桓等[25]對微地震數據進行EMD，利用同步壓縮提取高頻段有效信號，并與低頻部分重構，成功消除了微地震中的混疊噪聲。瞿明岳等[26]通過EMD改進TFPF方法，有效地消除電力線通信噪聲，降低了誤碼率。

根據微地震數據具有隨機性、非平穩性和時頻耦合等特性，兼顧EMD呈現的模態混疊及低信噪比去噪能力欠缺問題，本文提出基于CEEMDAN的時頻峰值濾波微地震數據去噪方法。模型測試和實際數據處理結果驗證了方法的有效性與實用性。

1 基本原理

1.1 自適應白噪聲完整集合經驗模態分解

CEEMDAN通過在各階段添加有限次的自適應白噪聲，實現在較小平均次數下重構誤差趨近于零。該算法在程序上實現了添加自適應白噪聲，可很好地抑制模態混淆且避免原始信號失真。具體步驟如下。

(1)在原始數據S(t)中分別多次添加自適應白噪聲Bq(t)，其中q表示添加噪聲次數，一般取10～50，本文取q=50，則第q次信號可表示為Sq(t)=S(t)+αqBq(t)(q=1，2，…，50)，其中αq為第q次加入白噪聲的標準差。CEEMDAN的一階IMF分量為

(1)

(2)構造新的待分解信號S(t)，S(t)=R1(t)+αqBq(t)，進行到第50次分解，得到CEEMDAN的二階IMF分量

(2)

余項R2(t)=S(t)-IMF2。

(3)重復步驟(1)和步驟(2)，到程序終止，共產生w個IMF，最后的余項為

(3)

1.2 樣本熵分析

樣本熵是一種與近似熵類似但精度更高的數據復雜度衡量指標，可對數據的復雜度進行量化分析。樣本熵的物理意義是表征信號中產生新模式的概率大小，以衡量時間序列復雜性。即新模式產生概率越小——樣本熵越小，序列自我相似度就越高； 樣本熵越大，樣本序列就越復雜。本文中樣本熵越大，IMFs中各個頻率越多，具體體現為噪聲含量高； 樣本熵越小，IMFs越趨向于有效信號。

通常，N個數據點組成的時間序列{x(n)=x(1)，x(2)，…，x(n)}的樣本熵的計算包括如下步驟：

(1)按序號組成1組m維的向量序列：Xm(1)，Xm(2)，…，Xm(i)，…，Xm(N-m+1)，其中Xm(i)={x(i)，x(i+1)，…，x(i+m-1)}，1≤i≤N-m+1，這些向量代表從第i個點開始的m個連續的值；

(2)將向量Xm(i)與Xm(j)之間的距離定義為對應元素中最大差值的絕對值

d[Xm(i)，Xm(j)]=max[|x(i+k)-x(j+k)|]

(4)

其中：1≤j≤N-m；j≠i；k=0，1，…，m-1。

(3)設相似容限(Tolerance for accepting matches)為r，通常取值為0.10～0.25個時間序列標準差。向量Xm(i)與Xm(j)的距離小于或等于r的數目為Bi，i≠j。其均值定義為

(5)

(6)

(5)將維數m增至m+1，重復以上步驟，可得Bm+1(r)；

(6)當N為有限值時，定義樣本熵

(7)

1.3 時頻峰值濾波(TFPF)

1.3.1 基本原理

TFPF的本質是基于Wigner-Ville分布(WVD)的瞬時頻率估計。待處理的含噪數據可表示為

s(t)=x(t)+n(t)

(8)

式中：x(t)是有效信號；n(t)是加性隨機噪聲。濾波的目的就是從含噪數據s(t)中恢復有效信號x(t)。TFPF可在不需假設條件的情況下很好地恢復有效信號，其具體步驟為：

(1)對含噪數據進行編碼，將其變為解析信號的形式； 對含噪數據s(t)進行頻率調制，得到單位幅度的解析信號

(9)

式中μ是頻率調制指數。

(2)取解析信號z(t)的偽Wigner-Ville分布(PWVD)的峰值，對解析信號進行瞬時頻率估計，作為有效信號x(t)的估計值

(10)

式中Wz(t，f)是解析信號z(t)的PWVD。

1.3.2 窗長選擇對TFPF濾波的影響

為了說明窗長選擇對TFPF算法在信號保幅和噪聲壓制上存在的問題，通過合成理論波形，用不同的窗長濾波，對比分析濾波結果。圖1a為有效信號，對其加入-1dB的噪聲，可見在噪聲背景下，該含噪數據(圖1b)中原始有效信號已極難辨別。

選用不同窗長的TFPF算法對s(t)進行濾波，從處理結果(圖2)可看出，各窗長的該算法均削減了噪聲，驗證了TFPF算法對噪聲壓制的有效性。進一步對比分析發現，長窗長具有更強的噪聲壓制能力，但信號幅度損失極大。

圖1 有效信號(a)及其含噪數據(b)

圖2 窗長分別為9(a)和13(b)的濾波結果

2 基于CEEMDAN的TFPF

本文提出基于CEEMDAN的TFPF去噪方法，可靈活選擇窗口長度，精準識別噪聲位置。詳細流程(圖3)如下。

(1)利用CEEMDAN分解信號，能使信號按從高到低頻率分解，得到若干個IMF分量。

(2)通過計算IMFs的樣本熵判斷含噪臨界點，高于此臨界值的IMFs需進行濾波。

(3)在選定需進行TFPF的IMFs后，根據地震數據的經驗公式對各個IMF選擇合適窗長

(11)

式中：fs為地震波的采樣頻率；fd為主頻。從該式可見，若fs增大或fd減小，則窗長相應變大。對于高頻IMFs，采用長窗長； 對于低頻IMFs，適用短窗長。這樣，就能達到保持有效信號幅度與壓制隨機噪聲間的平衡。

(4)最后將濾波后的IMFs與保留的IMFs重構，即可得到有效信號。

圖3 本文方法處理流程圖

3 數值實驗結果與分析

3.1 理論模型試算

選用兩個主頻分別為30Hz和25Hz的Ricker子波組成的波形進行實驗。采樣率為1ms，共計1000個采樣點。向該信號中加入高斯白噪聲，使信噪比(SNR)達到約5dB，有效信號基本被噪聲淹沒。

分別對原始的和加噪的理論模型數據(圖4)做CEEMDAN處理，得到分解后的多個IMFs(圖5)及其對應的樣本熵(圖6)。通過此樣本熵，可選擇需進行TFPF的IMFs； 再根據各個IMFs的頻率特性靈活地選擇窗長。選取該信號的前五個IMFs進行濾波。根據式(11)，對前五個IMFs分量分別選擇窗長為13、11、9、9、7，再將經TFPF濾波的信號與殘留IMFs構成完整濾波信號。

圖4 理論數據模型

傳統EMD去噪方法是舍棄前三個含大量噪聲的IMFs，以其余的IMFs重構信號。而固定窗長TFPF去噪方法，根據式(11)選定窗口長度為13。

圖5 理論模型數據經CEEMDAN分解后的IMFs

圖6 圖5各個IMFs對應的樣本熵

分別應用傳統EMD、固定窗長TFPF和本文CEEMDAN-TFPF三種去噪方法針對上述理論模型進行隨機噪聲壓制處理。可見傳統EMD去噪方法所得去噪結果(圖7a)中隨機噪聲得到一定程度壓制(相對于圖4b)，但仍有較多殘留噪聲； 固定窗長TFPF方法的去噪結果(圖7b)中噪聲得到較好壓制，但振幅與圖4a相差較大； 而CEEMDAN-TFPF方法的去噪結果(圖7c)不僅更徹底地壓制了隨機噪聲，且很好地保持了有效信號的振幅，兼顧了“保幅”與“去噪”。

為更充分、直觀地展示實驗結果，將原始數據、加噪信號和上述三種方法去噪結果置于同一坐標系下做波形對比(圖8)，并將其波峰、波谷局部放大(圖9)顯示，可見本文方法去噪結果更貼近原始數據。

圖7 理論模型數據不同方法去噪結果

3.2 實際微地震數據應用

選自M地區實際礦山微地震數據(圖10a)中包含大量隨機噪聲，嚴重影響后續監測分析與微地震源定位。對該數據做傅里葉變換，分析頻譜特征得知其主頻集中于10Hz附近(圖10a中)，信號在高頻部分分布較寬。因礦山微地震數據具低頻特性，故高頻部分的貢獻基本為隨機噪聲，應予壓制及去除。

圖8 不同方法去噪結果波形對比

圖9 對應圖8波峰(a)、波谷(b)波形的放大顯示

對上述數據做CEEMDAN處理，得到分解后的各個IMFs(圖11)及其對應的樣本熵(圖12)，據此樣本熵可選擇需進行TFPF的IMFs，且能由各個IMFs的頻率特性靈活地選擇窗長。選取該信號前六個IMFs進行濾波，根據式(11)對前六個IMFs分量分別選定窗長為9、9、7、7、5、5，最后將經TFPF濾波后的信號與殘留IMFs重構為完整濾波信號。對于傳統EMD去噪方法，也是舍棄前三個含大量噪聲的IMFs，以剩余IMFs重構信號； 固定窗長的TFPF去噪方法，據式(11)選定窗口長度為9。

分別應用傳統EMD、固定窗長TFPF和CEEMDAN-TFPF三種去噪方法對實際數據進行隨機噪聲壓制處理。可見傳統EMD方法所得去噪結果(圖10b上)中隨機噪聲得到一定程度的壓制，但仍殘留較多噪聲； 固定窗長TFPF方法所得去噪結果(圖10c上)中可看出噪聲得到了較好的壓制，但振幅與圖10a相差較大； 而本文方法去噪結果(圖10d上)中，不僅更徹底地(對比前兩種方法)壓制了隨機噪聲，且充分地保持了有效信號振幅，能達到保持有效信號幅度和壓制隨機噪聲間的平衡。

圖10 針對實際微地震數據應用不同方法去噪所得波形、頻譜及時頻對比

圖11 實際微地震數據經CEEMDAN分解后的IMFs

圖12 實際微地震數據各IMF對應的樣本熵

對比分析三種方法去噪后頻譜，可見傳統EMD去噪方法(圖10b中)能很好壓制高頻噪聲，但幾乎未壓制低頻噪聲； 固定窗長TFPF去噪方法(圖10c中)的中低頻噪聲未能很好地去除； 而本文方法對于全頻段噪聲都有很強壓制能力。從其時頻分布圖(圖10下)可見，去噪后數據中干擾信息減弱。

4 結論

傳統TFPF方法是選擇一個固定時窗進行濾波，不能同時兼顧信號去噪和幅值保持，過度壓制噪聲會造成有效信號幅值的嚴重衰減。而EMD存在模態混疊問題，且傳統EMD方法通過舍棄噪聲分量的處理方式不能有效壓制噪聲。本文提出一種基于CEEMDAN-TFPF的壓制隨機噪聲方法，利用EMD的尺度分解特性將原數據按頻率高低分解成一系列IMFs，再計算這些模態分量的樣本熵，確定濾波模態分量，然后對這些分量做自適應時窗長度的TFPF處理，最后將濾波后的分量與剩余分量相加得到最終濾波信號。通過對模擬數據和實際微地震數據的處理，得到以下認識和結論：

(1)CEEMDAN-TFPF去噪方法是自適應的，對參數的選擇無嚴格限制，能提高數據分析的效率。

(2)該方法能有效克服EMD的模態混疊問題，且可更精準地判別噪聲分量，對需濾波的IMF靈活選擇窗長，達到信號保幅與噪聲壓制間的平衡，提高微地震數據信噪比。

(3)與傳統EMD和定窗長TFPF兩種去噪方法相比，本文方法能在有效壓制隨機噪聲的同時，也很好地保護有效信號幅值，準確識別有效微地震信號，提高定位精度，保留原始微地震事件的數據特征，給后續微地震數據研究提供可靠基礎數據。