基于CEEMDAN的時(shí)頻峰值濾波微地震隨機(jī)噪聲壓制方法

陳毅軍 程 浩 鞏恩普 薛 林

(東北大學(xué)深部金屬礦山安全開(kāi)采教育部重點(diǎn)實(shí)驗(yàn)室，遼寧沈陽(yáng) 110819)

0 引言

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

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

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

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

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

根據(jù)微地震數(shù)據(jù)具有隨機(jī)性、非平穩(wěn)性和時(shí)頻耦合等特性，兼顧EMD呈現(xiàn)的模態(tài)混疊及低信噪比去噪能力欠缺問(wèn)題，本文提出基于CEEMDAN的時(shí)頻峰值濾波微地震數(shù)據(jù)去噪方法。模型測(cè)試和實(shí)際數(shù)據(jù)處理結(jié)果驗(yàn)證了方法的有效性與實(shí)用性。

1 基本原理

1.1 自適應(yīng)白噪聲完整集合經(jīng)驗(yàn)?zāi)B(tài)分解

CEEMDAN通過(guò)在各階段添加有限次的自適應(yīng)白噪聲，實(shí)現(xiàn)在較小平均次數(shù)下重構(gòu)誤差趨近于零。該算法在程序上實(shí)現(xiàn)了添加自適應(yīng)白噪聲，可很好地抑制模態(tài)混淆且避免原始信號(hào)失真。具體步驟如下。

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

(1)

(2)構(gòu)造新的待分解信號(hào)S(t)，S(t)=R1(t)+αqBq(t)，進(jìn)行到第50次分解，得到CEEMDAN的二階IMF分量

(2)

余項(xiàng)R2(t)=S(t)-IMF2。

(3)重復(fù)步驟(1)和步驟(2)，到程序終止，共產(chǎn)生w個(gè)IMF，最后的余項(xiàng)為

(3)

1.2 樣本熵分析

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

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

(1)按序號(hào)組成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個(gè)點(diǎn)開(kāi)始的m個(gè)連續(xù)的值；

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

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)設(shè)相似容限(Tolerance for accepting matches)為r，通常取值為0.10～0.25個(gè)時(shí)間序列標(biāo)準(zhǔn)差。向量Xm(i)與Xm(j)的距離小于或等于r的數(shù)目為Bi，i≠j。其均值定義為

(5)

(6)

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

(6)當(dāng)N為有限值時(shí)，定義樣本熵

(7)

1.3 時(shí)頻峰值濾波(TFPF)

1.3.1 基本原理

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

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

(8)

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

(1)對(duì)含噪數(shù)據(jù)進(jìn)行編碼，將其變?yōu)榻馕鲂盘?hào)的形式； 對(duì)含噪數(shù)據(jù)s(t)進(jìn)行頻率調(diào)制，得到單位幅度的解析信號(hào)

(9)

式中μ是頻率調(diào)制指數(shù)。

(2)取解析信號(hào)z(t)的偽Wigner-Ville分布(PWVD)的峰值，對(duì)解析信號(hào)進(jìn)行瞬時(shí)頻率估計(jì)，作為有效信號(hào)x(t)的估計(jì)值

(10)

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

1.3.2 窗長(zhǎng)選擇對(duì)TFPF濾波的影響

為了說(shuō)明窗長(zhǎng)選擇對(duì)TFPF算法在信號(hào)保幅和噪聲壓制上存在的問(wèn)題，通過(guò)合成理論波形，用不同的窗長(zhǎng)濾波，對(duì)比分析濾波結(jié)果。圖1a為有效信號(hào)，對(duì)其加入-1dB的噪聲，可見(jiàn)在噪聲背景下，該含噪數(shù)據(jù)(圖1b)中原始有效信號(hào)已極難辨別。

選用不同窗長(zhǎng)的TFPF算法對(duì)s(t)進(jìn)行濾波，從處理結(jié)果(圖2)可看出，各窗長(zhǎng)的該算法均削減了噪聲，驗(yàn)證了TFPF算法對(duì)噪聲壓制的有效性。進(jìn)一步對(duì)比分析發(fā)現(xiàn)，長(zhǎng)窗長(zhǎng)具有更強(qiáng)的噪聲壓制能力，但信號(hào)幅度損失極大。

圖1 有效信號(hào)(a)及其含噪數(shù)據(jù)(b)

圖2 窗長(zhǎng)分別為9(a)和13(b)的濾波結(jié)果

2 基于CEEMDAN的TFPF

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

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

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

(3)在選定需進(jìn)行TFPF的IMFs后，根據(jù)地震數(shù)據(jù)的經(jīng)驗(yàn)公式對(duì)各個(gè)IMF選擇合適窗長(zhǎng)

(11)

式中：fs為地震波的采樣頻率；fd為主頻。從該式可見(jiàn)，若fs增大或fd減小，則窗長(zhǎng)相應(yīng)變大。對(duì)于高頻IMFs，采用長(zhǎng)窗長(zhǎng)； 對(duì)于低頻IMFs，適用短窗長(zhǎng)。這樣，就能達(dá)到保持有效信號(hào)幅度與壓制隨機(jī)噪聲間的平衡。

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

圖3 本文方法處理流程圖

3 數(shù)值實(shí)驗(yàn)結(jié)果與分析

3.1 理論模型試算

選用兩個(gè)主頻分別為30Hz和25Hz的Ricker子波組成的波形進(jìn)行實(shí)驗(yàn)。采樣率為1ms，共計(jì)1000個(gè)采樣點(diǎn)。向該信號(hào)中加入高斯白噪聲，使信噪比(SNR)達(dá)到約5dB，有效信號(hào)基本被噪聲淹沒(méi)。

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

圖4 理論數(shù)據(jù)模型

傳統(tǒng)EMD去噪方法是舍棄前三個(gè)含大量噪聲的IMFs，以其余的IMFs重構(gòu)信號(hào)。而固定窗長(zhǎng)TFPF去噪方法，根據(jù)式(11)選定窗口長(zhǎng)度為13。

圖5 理論模型數(shù)據(jù)經(jīng)CEEMDAN分解后的IMFs

圖6 圖5各個(gè)IMFs對(duì)應(yīng)的樣本熵

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

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

圖7 理論模型數(shù)據(jù)不同方法去噪結(jié)果

3.2 實(shí)際微地震數(shù)據(jù)應(yīng)用

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

圖8 不同方法去噪結(jié)果波形對(duì)比

圖9 對(duì)應(yīng)圖8波峰(a)、波谷(b)波形的放大顯示

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

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

圖10 針對(duì)實(shí)際微地震數(shù)據(jù)應(yīng)用不同方法去噪所得波形、頻譜及時(shí)頻對(duì)比

圖11 實(shí)際微地震數(shù)據(jù)經(jīng)CEEMDAN分解后的IMFs

圖12 實(shí)際微地震數(shù)據(jù)各IMF對(duì)應(yīng)的樣本熵

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

4 結(jié)論

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

(1)CEEMDAN-TFPF去噪方法是自適應(yīng)的，對(duì)參數(shù)的選擇無(wú)嚴(yán)格限制，能提高數(shù)據(jù)分析的效率。

(2)該方法能有效克服EMD的模態(tài)混疊問(wèn)題，且可更精準(zhǔn)地判別噪聲分量，對(duì)需濾波的IMF靈活選擇窗長(zhǎng)，達(dá)到信號(hào)保幅與噪聲壓制間的平衡，提高微地震數(shù)據(jù)信噪比。

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