胡瑞卿 王彥春 尹志恒 王 偉

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

0 引言

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

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

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

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

1 方法原理

1.1 CEEMDAN算法原理

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

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

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

(1)生成含噪信號集

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

(1)

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

(2)

(3)計算一階殘差量

(3)

(4)計算二階IMF

(4)

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

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

(5)

(6)計算k+1階分量

(6)

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

最終殘差滿足

(7)

信號可表示為

(8)

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

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

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

(a)及其EEMD

(b)和CEEMDAN

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

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

圖3 實際微地震信號(上)及其頻譜(下)

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

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

1.2 主成分分析與特征提取

通常采用三分量檢波器采集微地震數(shù)據(jù)，有效信號在各分量上能量不同，但初至點一致、有效成分波形高度相似，故可通過對三分量信號先進行PCA，再進行特征提取。微地震信號中常出現(xiàn)的噪聲為近似滿足獨立高斯分布的自然噪聲與低頻震蕩。通過CEEMDAN分解，不同階次的IMF承載原始信號中的不同信息。由于原信號中噪聲能量強于有效成分，故在各階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進行PCA，統(tǒng)計各主成分所占的能量百分比(圖6)，可見4、5、6階IMF的第一主成分具有明顯的能量優(yōu)勢。其原因在于有效成分能量主要落于這三階IMF的第一主成分中，而其他各階IMF中，前三項主成分能量差異不大，主要由高能量的噪聲構(gòu)成。通過將高能量噪聲等特征較弱的成分分配到各階IMF中，而特征明顯的有效成分能量則集中分配于前幾項主成分中，剔除后幾項主成分，可對有效成分進行提取。

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

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

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

圖7 實測含噪信號及其各階IMF

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

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

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

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

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

1.4 處理流程

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

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

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

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

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

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

2.1 模型Ⅰ

設(shè)計模型Ⅰ合成信號為

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

(16)

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

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

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

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

2.2 模型Ⅱ

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

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

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

圖11 模型Ⅱ不同噪聲能量下的合成信號頻譜

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

通過對2000組不同信噪比信號的測試，本文方法與常規(guī)STA/LTA-AIC法的拾取精度誤差統(tǒng)計結(jié)果如圖12所示?？梢姳疚姆椒ǖ氖叭【仍赟NR為-2～-20dB時，遠(yuǎn)優(yōu)于常規(guī)STA/LTA-AIC方法。

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

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

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

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

3 結(jié)論

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