結(jié)合CEEMDAN和主成分分析的低信噪比微地震初至信號(hào)檢測

胡瑞卿 王彥春 尹志恒 王 偉

(①中國地質(zhì)大學(xué)(北京)地下信息探測技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室，北京 100089；②中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院，北京 100089；③中石化石油工程地球物理有限公司，北京 100020；④東方地球物理公司培訓(xùn)中心，河北涿州 072750)

0 引言

從微地震數(shù)據(jù)中提取微地震事件，可實(shí)時(shí)監(jiān)測裂縫的空間位置與尺度，為后續(xù)壓裂參數(shù)的調(diào)整與控制提供數(shù)據(jù)支撐，微地震監(jiān)測對油氣開發(fā)具有重要的指導(dǎo)意義[1]。

當(dāng)前微地震監(jiān)測面臨的主要問題是資料信噪比較低[2-4]，甚至有效信號(hào)被噪聲完全掩蓋。這對以初至拾取為基礎(chǔ)的震源定位、壓裂預(yù)測、機(jī)制分析等后續(xù)工作造成了嚴(yán)重影響。因此，針對低信噪比或超低信噪比的微地震資料，選擇合理有效的初至拾取方法是微地震資料處理的關(guān)鍵。

通常，初至拾取以信號(hào)與噪聲在振幅、頻率、偏振等方面的差異為基礎(chǔ)。常采用長短時(shí)窗平均能量比(STA/LTA)法[5-7]、基于自回歸模型的Aksike信息準(zhǔn)則(AIC)法[8]、相關(guān)法[9-12]、偏振分析法[13]、二分法[14]、分維相關(guān)法、時(shí)頻分析法、神經(jīng)網(wǎng)絡(luò)法[15]、獨(dú)立分量法[16]、基于高階統(tǒng)計(jì)量的PAI-S/K法[17]以及基于小波分解的高階統(tǒng)計(jì)初至拾取方法[18]。Coppens[19]對不同時(shí)窗內(nèi)不同炮檢距道集的能量進(jìn)行統(tǒng)計(jì)相關(guān)分析，提出一種自動(dòng)定位初至點(diǎn)的方法，有效提高了拾取精度；鄒紅星[20]提出一種基于自回歸模型的判讀法，該方法具有一定自適應(yīng)能力，可對不同噪聲影響程度的接受信息進(jìn)行初至拾取。賈瑞生等[21]、蔡劍華等[22]基于經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和本征模態(tài)函數(shù)(IMF)重構(gòu)對低信噪比信號(hào)去噪，然后應(yīng)用正交性指標(biāo)與能量保存度指標(biāo)對傳統(tǒng)EMD中的端點(diǎn)效應(yīng)進(jìn)行壓制，從而實(shí)現(xiàn)震相初至的拾取，其應(yīng)用效果優(yōu)于常規(guī)基于AIC準(zhǔn)則的拾取方法。譚玉陽等[23]提出了基于信號(hào)信噪比構(gòu)造不同檢測函數(shù)的初至拾取算法(SLPEA)，通過綜合考慮有效信號(hào)與環(huán)境噪聲間的多種特征差異，提高拾取過程的抗噪能力，但時(shí)窗的選擇依據(jù)并不明確，需依賴經(jīng)驗(yàn)或反復(fù)調(diào)試。

在前人的基礎(chǔ)上，本文基于三分量間有效成分的特征相似性，采用CEEMDAN算法對各分量信號(hào)進(jìn)行分解，然后用各階IMF構(gòu)建三分量IMF矩陣，再進(jìn)行主成分分析(PCA)，提取主成分特征，達(dá)到對低信噪比微地震資料有效波初至拾取的目的。通過對不同信噪比合成信號(hào)與實(shí)測信號(hào)的處理，證實(shí)該方法具有很強(qiáng)的抗噪能力，適用于有效信號(hào)被噪聲掩蓋的低信噪比微地震資料。

1 方法原理

1.1 CEEMDAN算法原理

EMD可將信號(hào)x(t)分解為一系列IMF，但信號(hào)需要滿足兩個(gè)條件：①極值點(diǎn)數(shù)與過零點(diǎn)數(shù)相等或相差不超過一個(gè)；②上、下包絡(luò)線的均值在各處均為0。

對于實(shí)際非平穩(wěn)時(shí)變信號(hào)，常規(guī)EMD會(huì)受到模態(tài)混疊的影響。通過引入正態(tài)分布的白噪聲，形成噪聲輔助信號(hào)處理方法(NADA)，即為常用的EEMD[24]。而本文采用基于改進(jìn)EEMD方法——CEEMDAN算法，通過在原始信號(hào)上加入有限方差約束的多組獨(dú)立同分布自適應(yīng)白噪聲，得到含噪信號(hào)的集合。該改進(jìn)方法可進(jìn)一步降低迭代次數(shù)，壓縮頻率混疊區(qū)域，提高收斂性能，對于非平穩(wěn)信號(hào)不同頻率成分具有更高的分辨能力。

對信號(hào)x(n)的CEEMDAN的實(shí)現(xiàn)步驟如下。

(1)生成含噪信號(hào)集

xi(n)=x(n)+wi(n)

(1)

式中wi(n)(i=1，2，…，I)為滿足高斯分布的噪聲，I為集合樣本數(shù)；

(2)

(3)計(jì)算一階殘差量

(3)

(4)計(jì)算二階IMF

(4)

式中：Ej(·)表示信號(hào)的j階IMF；εj為控制白噪能量的參數(shù)，Wu等[24]認(rèn)為對于主頻較高的信號(hào)，取值應(yīng)較小，反之較大，本文采用較小常量進(jìn)行試算。

(5)對于k(k=2，3，…，K，K為設(shè)置的最高IMF階次)階分量，計(jì)算k階殘差

(5)

(6)計(jì)算k+1階分量

(6)

(7)重復(fù)步驟(5)、步驟(6)直到殘差不可再分解或達(dá)到最高IMF階次。

最終殘差滿足

(7)

信號(hào)可表示為

(8)

與傳統(tǒng)EEMD方法相比，CEEMDAN算法最終得到的各階IMF分量對原始信號(hào)頻率成分的表征不同。如圖1所示，對3～6階IMF分量進(jìn)行頻譜對比，可見傳統(tǒng)EEMD各分量間能量差異較大，頻率混置區(qū)較大。而CEEMDAN方法得到的各IMF分量間能量均衡，頻率混疊區(qū)較窄，對于非平穩(wěn)信號(hào)，不同頻率成分具有更高的分辨能力。通過對心電圖信號(hào)(圖1a)的測試考察其計(jì)算性能，該信號(hào)具有非平衡性、非線性與隨機(jī)性，具有豐富的頻率成分，適合用于時(shí)頻分析方法測試。該測試信號(hào)的各階分量的迭代次數(shù)箱形圖如圖2所示，可見CEEMDAN各階分量的迭代次數(shù)均顯著降低。對于主要承載有效信息的3階至6階IMF分量，CEEMDAN方法的計(jì)算量約為EEMD方法的1/3。

圖3為實(shí)測微地震信號(hào)及其頻譜，信號(hào)采樣間隔為0.25ms。圖4為實(shí)際微地震數(shù)據(jù)傳統(tǒng)EMD方法分解的6～9階IMF及其頻譜。各階IMF在靠近端點(diǎn)處產(chǎn)生較大幅度的振蕩，并且在6階和7階IMF對應(yīng)的頻譜上均顯示了兩個(gè)明顯峰值(圖4a、圖4b中的藍(lán)色區(qū)域)，明顯包含兩組頻率成分，表明存在明顯的模態(tài)混疊現(xiàn)象[25-26]。

圖5為實(shí)際微地震數(shù)據(jù)CEEMDAN方法分解的6～9階IMF分量及其頻譜，端點(diǎn)附近的異常振蕩得到有效壓制，且各階IMF對應(yīng)的頻譜均無明顯的多峰現(xiàn)象，模態(tài)混疊現(xiàn)象得到有效消除。這是因?yàn)镃EEMDAN方法利用高斯白噪聲頻率均勻分布的統(tǒng)計(jì)特征，向原始信號(hào)中加入不同的白噪聲，使信號(hào)在不同尺度上具有連續(xù)性，并通過在分解時(shí)的每一階段添加特定白噪聲，并計(jì)算一個(gè)唯一殘差以得到每個(gè)IMF分量，使原始信號(hào)的模態(tài)頻譜分離效果更好，端點(diǎn)效果的影響顯著降低。

(a)及其EEMD

(b)和CEEMDAN

(c)方法分解的3～7階IMF頻譜對比

圖2 不同分解方法各階IMF迭代次數(shù)統(tǒng)計(jì)結(jié)果

圖3 實(shí)際微地震信號(hào)(上)及其頻譜(下)

圖4 實(shí)際微地震信號(hào)傳統(tǒng)EMD方法6～9階IMF分解結(jié)果(上)及其頻譜(下)

圖5 實(shí)際微地震信號(hào)CEEMDAN方法6～9階IMF分解結(jié)果(上)及其頻譜(下)

1.2 主成分分析與特征提取

通常采用三分量檢波器采集微地震數(shù)據(jù)，有效信號(hào)在各分量上能量不同，但初至點(diǎn)一致、有效成分波形高度相似，故可通過對三分量信號(hào)先進(jìn)行PCA，再進(jìn)行特征提取。微地震信號(hào)中常出現(xiàn)的噪聲為近似滿足獨(dú)立高斯分布的自然噪聲與低頻震蕩。通過CEEMDAN分解，不同階次的IMF承載原始信號(hào)中的不同信息。由于原信號(hào)中噪聲能量強(qiáng)于有效成分，故在各階IMF中，仍存在噪聲能量，可采用PCA對三分量IMF矩陣提取主成分。三分量數(shù)據(jù)的k階IMF矩陣可表示為

(9)

Ek=M(k)·W

(10)

其中第一主成分對應(yīng)的特征向量，應(yīng)該滿足

(11)

寫成矩陣形式

W1=arg max‖W‖=1[(WTW)-1WTMT(k)M(k)W]

(12)

W1使第一主成分具有最高能量，第一主成分能量可表示為

(13)

(14)

(15)

通過對相鄰三道三分量微地震記錄的各階IMF進(jìn)行PCA，統(tǒng)計(jì)各主成分所占的能量百分比(圖6)，可見4、5、6階IMF的第一主成分具有明顯的能量優(yōu)勢。其原因在于有效成分能量主要落于這三階IMF的第一主成分中，而其他各階IMF中，前三項(xiàng)主成分能量差異不大，主要由高能量的噪聲構(gòu)成。通過將高能量噪聲等特征較弱的成分分配到各階IMF中，而特征明顯的有效成分能量則集中分配于前幾項(xiàng)主成分中，剔除后幾項(xiàng)主成分，可對有效成分進(jìn)行提取。

1.3 IMF主成分加權(quán)重構(gòu)

對三分量數(shù)據(jù)進(jìn)行CEEMD后，對各階IMF進(jìn)行PCA ，選取目標(biāo)主成分，最后加權(quán)重構(gòu)。通過CEEMD分解得到的各階IMF，反映出信號(hào)的不同成分。以圖7含噪信號(hào)為例，可以看出，1、2階IMF主要反映信號(hào)的大體趨勢，7～9階IMF則主要承載高頻信息。

圖6 各階IMF主成分能量統(tǒng)計(jì)

圖7 實(shí)測含噪信號(hào)及其各階IMF

在低信噪比情況下，高頻信息主要為噪聲，那么有效信號(hào)的時(shí)域特征則主要由中間的3～6階IMF反映。因此，重構(gòu)時(shí)，可將1、2階和7～9階IMF 的權(quán)重減少，將3～6階IMF的權(quán)重加大，以便更加突出有效信號(hào)的特征。

對于重構(gòu)時(shí)各級(jí)主成分權(quán)重系數(shù)的選取，本文的思路如下。

(1)歸一化處理：將各階IMF主成分單獨(dú)重構(gòu)回時(shí)域波形，對時(shí)域波形進(jìn)行振幅歸一化處理；

(2)計(jì)算各階時(shí)域波形與原始信號(hào)所提取的子波間的相關(guān)系數(shù)；

(3)以相關(guān)系數(shù)絕對值作為各階IMF主成分的重構(gòu)權(quán)值。

1.4 處理流程

針對低信噪比條件下的微地震信號(hào)，首先對原始信號(hào)進(jìn)行CEEMDAN；然后基于同一檢波器中三個(gè)分量間的有效成分相關(guān)性，在IMF域內(nèi)從高能量噪聲中提取出有效成分特征，并剔除其余成分；最后將剩余的IMF主成分加權(quán)重構(gòu)，得到體現(xiàn)有效信號(hào)特征的時(shí)域波形。具體處理流程如下：

(1)對三分量微地震資料進(jìn)行CEEMDAN；

(2)對三分量中各階IMF分量進(jìn)行PCA，以75%為門檻值，剔除次要成分；

(3)按不同權(quán)值進(jìn)行IMF主成分重構(gòu)，得到主成分在時(shí)域內(nèi)波形。

2 模型測試與實(shí)際應(yīng)用

采用兩組模型對本文方法進(jìn)行可行性測試與適應(yīng)性測試。模型Ⅰ采用合成信號(hào)與白噪聲，測試該方法對噪聲掩蓋下的具有內(nèi)在固有特征成分檢測能力。模型Ⅱ采用實(shí)際有效信號(hào)，疊加白噪聲，通過調(diào)整噪聲能量，測試該方法適用的信噪比下限。通過對2000組信號(hào)測試，對比基于本文方法與常規(guī)基于AIC信息準(zhǔn)則的長短視時(shí)窗(STA/LTA-AIC)方法的拾取誤差。最后對具有不同品質(zhì)的低信噪比實(shí)測資料進(jìn)行了處理。

2.1 模型Ⅰ

設(shè)計(jì)模型Ⅰ合成信號(hào)為

x(t)=λy(t)?[exp(-4.5t)sin(100πt)]+w(t)

(16)

式中：y(t)為標(biāo)準(zhǔn)Ricker子波；“?”為褶積運(yùn)算符；w(t)為白噪聲；λ為振幅因子，控制合成信號(hào)的信噪比。初至?xí)r間設(shè)置為500ms，采樣間隔為0.25ms，用式(16)生成能量隨時(shí)間呈指數(shù)衰減的正弦信號(hào)與標(biāo)準(zhǔn)Ricker子波的褶積作為有效信號(hào)(圖8中紅線)，調(diào)整λ得到一系列不同信噪比(SNR)的測試信號(hào)(圖8中黑線)。

應(yīng)用本文方法對測試信號(hào)進(jìn)行處理，結(jié)果如圖9所示。圖9a中由于信噪比過低，在500～650ms可以看到持續(xù)的低能量偽響應(yīng)，隨著SNR的提高，該偽響應(yīng)的相對能量降低。經(jīng)過本文方法處理后，初至事件的響應(yīng)清晰度與分辨率都有顯著提升，應(yīng)用常規(guī)方法在該結(jié)果上均能得到高精度的初至拾取結(jié)果。需要指明，有效事件的檢測響應(yīng)，需要一致地出現(xiàn)在三分量信號(hào)中的至少兩個(gè)分量上，否則不應(yīng)被視為微地震事件的有效響應(yīng)。如圖9e中所示，在約210ms處出現(xiàn)單一響應(yīng)，而在其他分量上則并沒有同時(shí)出現(xiàn)同步響應(yīng)，因此將其作為偽響應(yīng)剔除。該模型測試結(jié)果表明，本文方法最低可適用于SNR=-20dB的信號(hào)。當(dāng)SNR低于-20dB，有效信號(hào)因能量過低而無法識(shí)別。

圖8 模型Ⅰ不同SNR的合成信號(hào)(紅線代表有效信號(hào)、黑線代表加噪信號(hào))

圖9 模型Ⅰ不同SNR數(shù)據(jù)本文方法處理結(jié)果

2.2 模型Ⅱ

江漢油田水力壓裂實(shí)驗(yàn)中采集到高信噪比微地震信號(hào)(圖10a)，采樣頻率為4000Hz，有效頻帶為50～220Hz。利用常規(guī)方法(如STA/LTA-AIC等)拾取的P波初至為287.50ms，S波初至為527.25ms。加入不同能量的白噪聲，獲得模型Ⅱ不同SNR的合成數(shù)據(jù)(圖10b～圖10f)，當(dāng)SNR降至-6.9787dB時(shí)(圖10c)，STA/LTA方法無法有效拾取初至。不同SNR的測試信號(hào)頻譜如圖11所示，隨著噪聲能量的加強(qiáng)，有效頻段的能量優(yōu)勢逐漸消失。本文方法對不同SNR的合成數(shù)據(jù)(圖10b～圖10f)的處理結(jié)果如圖10g～圖10k所示，可以看出，即使SNR降到-19.0199dB，處理后仍可對微地震事件進(jìn)行有效檢測。

圖10 模型Ⅱ合成信號(hào)及本文方法處理結(jié)果

(a)原始微地震信號(hào)；(b)SNR為-0.9581dB的合成信號(hào)；(c)SNR為-6.9787dB的合成信號(hào)；(d)SNR為-12.9993dB的合成信號(hào)；(e)SNR為-16.5211dB的合成信號(hào)；(f)SNR為-19.0199dB的合成信號(hào)；(g)圖10b合成信號(hào)的檢測響應(yīng)；(h)圖10c合成信號(hào)的檢測響應(yīng)；(i)圖10d合成信號(hào)的檢測響應(yīng)；(j)圖10e合成信號(hào)的檢測響應(yīng)；(k)圖10f合成信號(hào)的檢測響應(yīng)

圖11 模型Ⅱ不同噪聲能量下的合成信號(hào)頻譜

圖12 模型Ⅱ不同方法初至拾取誤差統(tǒng)計(jì)

通過對2000組不同信噪比信號(hào)的測試，本文方法與常規(guī)STA/LTA-AIC法的拾取精度誤差統(tǒng)計(jì)結(jié)果如圖12所示。可見本文方法的拾取精度在SNR為-2～-20dB時(shí)，遠(yuǎn)優(yōu)于常規(guī)STA/LTA-AIC方法。

2.3 實(shí)際資料應(yīng)用

實(shí)測井下壓裂微地震數(shù)據(jù)(圖13a)采樣頻率為4000Hz，部分記錄數(shù)據(jù)品質(zhì)較差，常規(guī)方法拾取效果不理想。圖13a中藍(lán)色顯示的記錄上可肉眼觀測到疑似P波與S波初至波峰，約在1080ms與1130ms處(激發(fā)時(shí)刻不為0)，在紅色顯示的紀(jì)錄上僅在1230ms附近可見單一疑似初至波峰，而在黑色顯示的紀(jì)錄上，無任何突變波形體現(xiàn)。本文方法處理結(jié)果如圖13b所示，藍(lán)、紅色記錄均可見明確的檢測響應(yīng)。

圖13 實(shí)際井下微地震數(shù)據(jù)(a)及本文方法處理結(jié)果(b)

而對于黑色記錄，本文方法沒有得到有效的處理結(jié)果，原因在于有效成分能量過低，處理后仍被噪聲殘留成分所淹沒。在藍(lán)、紅色記錄的檢測結(jié)果基礎(chǔ)上，采用常規(guī)STA/LTA-AIC方法拾取初至，以圖13b中不同底色(藍(lán)—黃、黃—灰)的分界線標(biāo)識(shí)P波、S波初始時(shí)刻。

3 結(jié)論

本文針對噪聲能量較大、有效信號(hào)被淹沒的初至拾取難題，提出在CEEMDAN的基礎(chǔ)上，利用微地震記錄中三分量信號(hào)有效成分具有高度相似性，對三分量信號(hào)的各階IMF進(jìn)行主成分分析，保留有效成分的特征，剔除干擾成分，有效實(shí)現(xiàn)了低信噪比資料中的初至特征檢測。實(shí)際資料的處理結(jié)果表明該方法對于在有效信號(hào)被噪聲掩蓋的情況下仍能有效、精確地拾取初至。