基于改進(jìn)長短時(shí)窗比值及優(yōu)化變分模態(tài)分解的微震初至拾取算法＊

孟 娟 吳燕雄 李亞南

（中國河北三河 065201 防災(zāi)科技學(xué)院）

引言

微震監(jiān)測是通過對(duì)微小地震事件的實(shí)時(shí)監(jiān)測來實(shí)現(xiàn)對(duì)地壓穩(wěn)定性或地下施工過程的安全監(jiān)控與管理，在開采、礦藏勘探、大壩監(jiān)測、礦山安全等工程領(lǐng)域，尤其是井下水力壓裂、救援等方面應(yīng)用廣泛。準(zhǔn)確可靠的P 波初至到時(shí)拾取是微震監(jiān)測系統(tǒng)的關(guān)鍵技術(shù)之一，但井下工作環(huán)境復(fù)雜，實(shí)時(shí)監(jiān)測的微震數(shù)據(jù)中經(jīng)常混入較多的噪聲干擾，加之微震信號(hào)本身能量微弱，頻率范圍廣，且持續(xù)時(shí)間短，噪聲干擾的增加無疑給P 波到時(shí)拾取帶來了巨大挑戰(zhàn)，因此，快速、準(zhǔn)確、可靠地拾取微震P 波到時(shí)成為業(yè)界研究的重點(diǎn)。

微震P 波到時(shí)拾取算法，目前應(yīng)用較廣泛的有長短時(shí)窗比值（short term averaging/long term averaging，縮寫為STA/LTA）（劉晗，張建中，2014；劉曉明等，2017）、赤池信息準(zhǔn)則（Akaike Information Criterion，縮寫為AIC）（趙大鵬等，2013；朱權(quán)潔等，2018）。STA/LTA 快速簡單，但在信噪比（signal-to-noise ratio，縮寫為SNR）較低或微震不明顯時(shí)準(zhǔn)確度低；AIC 方法較穩(wěn)定，但在低SNR 時(shí)準(zhǔn)確度有所下降。針對(duì)STA/LTA 算法的不足，劉晗和張建中（2014）引入長短時(shí)窗內(nèi)信號(hào)幅度的標(biāo)準(zhǔn)差之比的加權(quán)系數(shù)來放大信號(hào)，增加檢測靈敏度，但該系數(shù)對(duì)噪聲亦有放大作用，使檢測準(zhǔn)確度降低；劉曉明等（2017）基于信號(hào)一階導(dǎo)數(shù)引入權(quán)重因子K構(gòu)建新特征函數(shù)，從而快速反映信號(hào)振幅及頻率的變化，提高檢測靈敏度，但該因子為全局性值，未根據(jù)信號(hào)幅度變化作靈活調(diào)整，存在一定的缺陷。

單一方法進(jìn)行初至拾取各有適用條件和局限性，近年，利用各種算法優(yōu)勢的綜合分析方法備受關(guān)注，尤其是將地震信號(hào)先進(jìn)行變換或分解，再用AIC 準(zhǔn)則或高階統(tǒng)計(jì)量提取P 波到時(shí)，能大大提高識(shí)別準(zhǔn)確度。例如：希爾伯特-黃變換與AIC 結(jié)合（賈瑞生等，2015）；小波分解與AIC 結(jié)合（Gaci，2014）；小波包分解AIC 的P 波初至拾?。ㄌ飪?yōu)平，趙愛華，2016）；基于離散小波變換（discrete wavelet transform，縮寫為DWT）與峰度進(jìn)行P 波到時(shí)精確估計(jì)（Liet al，2016）；利用多尺度離散平穩(wěn)小波變換（discrete stationary wavelet transform，縮寫為DSWT）和峰度/峭度實(shí)現(xiàn)P 波自動(dòng)拾取（Charles，Ge，2018）等。

經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition，縮寫為EMD）是近年發(fā)展起來的自適應(yīng)信號(hào)分解方法，非常適合處理非線性非平穩(wěn)的地震信號(hào)，其被運(yùn)用于地震波初至拾取也處于不斷探索和研究中。例如：Liu 等（2017）利用EMD 方法分解信號(hào)，找到本征模函數(shù)（intrinsic mode functions，縮寫為IMF）分量1（IMF1）的瞬時(shí)頻率突增位置，再對(duì)該段記錄利用AIC 準(zhǔn)則來確認(rèn)信號(hào)初至?xí)r間；Li 等（2017）對(duì)地震信號(hào)EMD 后，以IMF5 的P 波到時(shí)為原信號(hào)的P 波到時(shí)；Shang 等（2018）對(duì)地震信號(hào)進(jìn)行EMD，對(duì)IMF1 和IMF2 去噪后重構(gòu)原信號(hào)，再利用余弦函數(shù)去除工頻噪聲，最后基于AIC 計(jì)算P 波到時(shí)；Kirbas 和Peker （2017）對(duì)地震信號(hào)進(jìn)行EMD，選擇IMF3 計(jì)算滑動(dòng)時(shí)間窗口長度的Teager 能量算子（Teager-Kalser energy operator，縮寫為TKEO），以TKEO 大于設(shè)定閾值且為最大值的位置為P 波到時(shí)；李偉等（2018）利用集合經(jīng)驗(yàn)?zāi)B(tài)分解（ensemble empirical mode decomposition，縮寫為EEMD）進(jìn)行信號(hào)分解后再用Hankel 矩陣奇異值分解（singular value decomposition with Hankle matrix，縮寫為Hankle SVD）法降噪，最后利用AIC 進(jìn)行微震初至拾取。這些算法通過EMD 分解原記錄，利用某些IMF 的細(xì)節(jié)特征進(jìn)行到時(shí)拾取，取得了一定的效果，但EMD 算法本身存在模態(tài)混疊及偽分量的問題，這無疑會(huì)影響拾取的準(zhǔn)確度。

基于以上背景，為解決低信噪比下初至拾取準(zhǔn)確率低的問題，本文擬利用STA/LTA 快速和AIC 識(shí)別穩(wěn)定的特點(diǎn)，同時(shí)結(jié)合變分模態(tài)分解（variational mode decomposition，縮寫為VMD）可有效減少偽分量及抑制模態(tài)混疊的優(yōu)勢，提出一種基于改進(jìn)STA/LTA 和自適應(yīng)VMD-峰度AIC 的兩步識(shí)別P 波初至的算法。首先，改進(jìn)傳統(tǒng)STA/LTA 算法，引入反映信號(hào)幅度實(shí)時(shí)變化的權(quán)重因子，構(gòu)建能反映信號(hào)能量和頻率變化的特征函數(shù)，提高P 波初次識(shí)別的準(zhǔn)確度；其次，對(duì)初次判別的P 波到時(shí)前后2—3 s 信號(hào)進(jìn)行優(yōu)化VMD，使IMF 分量既能準(zhǔn)確地描述原信號(hào)的細(xì)節(jié)特征，又可防止出現(xiàn)過分解和模態(tài)混疊；最后，對(duì)分解后的IMF 進(jìn)行峰度AIC 拾取，并根據(jù)各IMF 的能量比進(jìn)行綜合加權(quán)，實(shí)現(xiàn)P 波到時(shí)的精確拾取。

1 基本算法原理

1.1 STA/LTA 算法原理

STA/LTA 算法根據(jù)噪聲與有效地震信號(hào)的振幅、能量或頻率等特征的變化趨勢來拾取P 波到時(shí)，具體方法為分別計(jì)算長、短時(shí)窗內(nèi)定義的特征函數(shù)值的STA 和LTA，其中STA 刻畫的是地震信號(hào)的變化趨勢，而LTA 刻畫的是背景噪聲的變化趨勢。當(dāng)?shù)卣鹗录絹頃r(shí)，引起STA 的變化快于LTA，即STA/LTA 會(huì)出現(xiàn)突增。當(dāng)STA/LTA 大于預(yù)先設(shè)定的閾值時(shí)，則判定為有效地震P 波到達(dá)。特征函數(shù)是基于原始地震記錄構(gòu)建的能靈敏反映信號(hào)幅度或頻率變換的新序列，目前廣泛使用的有

STA/LTA 的計(jì)算方法為

式中，CF（i）為i時(shí)刻地震信號(hào)的特征函數(shù)值，LSTA為短時(shí)窗長度，LLTA為長時(shí)窗長度。STA/LTA 拾取結(jié)果的準(zhǔn)確性受特征函數(shù)、時(shí)窗長度及閾值的影響，需要構(gòu)建靈敏的特征函數(shù)，選取合理的閾值和時(shí)窗長度。

1.2 AIC 算法基本原理

AIC 算法假定地震信號(hào)可分為若干局部平穩(wěn)部分，每部分均可模擬為自回歸過程。將地震到來前后的時(shí)間序列假設(shè)為兩個(gè)不同的平穩(wěn)序列，則其AIC 值的最小值點(diǎn)為兩個(gè)序列的分界點(diǎn)，以此確定P 波的初達(dá)時(shí)間。對(duì)于長度為N的序列，其AIC 值可表示為可表示為

式中，k為數(shù)據(jù)窗口內(nèi)的采樣點(diǎn)數(shù)，var 為序列方差，則AIC 最小值對(duì)應(yīng)P 波初至到時(shí)。根據(jù)趙大鵬等（2013）的方法，用峰度（Kurtosis）作為特征函數(shù)替代式（3）中的方差能更好地識(shí)別微震時(shí)間，即Kur-AIC 為

峰度AIC 方法不依賴于時(shí)窗長度，計(jì)算穩(wěn)定，對(duì)微弱P 波亦有較好的分辨能力。AIC 算法雖穩(wěn)定，但當(dāng)信號(hào)信噪比較低時(shí)其準(zhǔn)確度會(huì)大大降低，因?yàn)榧词故窃肼曈涗浀腁IC 值也會(huì)出現(xiàn)一個(gè)局部最小值，故AIC 不能準(zhǔn)確判斷一段記錄中是否存在有效微震信號(hào)，只能用于拾取微震P 波的初至?xí)r間。因此，一般先利用STA/LTA 初步確定微震P 波到時(shí)區(qū)間，再利用AIC進(jìn)行拾取，降低誤差。

1.3 變分模態(tài)分解的原理

變分模態(tài)分解（VMD）是一種自適應(yīng)的信號(hào)分解方法（Dragomiretskiy，Zosso，2014），能將信號(hào)分解為若干從高頻到低頻，帶限、準(zhǔn)正交的IMF 分量之和，且每個(gè)IMF 都是一個(gè)窄帶的調(diào)幅調(diào)頻信號(hào)。VMD 算法包括構(gòu)造約束變分模型和變分模型的求解，其實(shí)質(zhì)是在“各IMF 之和為原信號(hào)，且IMF 的估計(jì)帶寬之和最小”的約束條件下搜尋L個(gè)IMF 分量。不同于EMD 分解分量個(gè)數(shù)不可人為設(shè)定，VMD 的IMF 分量個(gè)數(shù)可設(shè)定。通過迭代搜尋變分模型的最優(yōu)解，VMD 確定L個(gè)IMF 的中心頻率和帶寬，從而實(shí)現(xiàn)信號(hào)從低頻到高頻的分解。相比EMD，VMD 模態(tài)穩(wěn)定性好，收斂速度快，具有更強(qiáng)的抗噪聲能力和局部分解能力，且能有效克服EMD 模態(tài)混疊及存在偽分量的不足，比較適合進(jìn)行地震資料處理（Xueet al，2016；Liet al，2018）。

VMD 能準(zhǔn)確地分解出每一個(gè)IMF 分量，分解后的IMF 瞬時(shí)頻譜和幅值具有較高的分辨率，不僅能刻畫原信號(hào)的時(shí)頻特征，也能更精細(xì)地分析信號(hào)，更好地描述信號(hào)的局部突變或漸變特征（鄭小霞等，2017）。因此可先基于VMD 對(duì)地震信號(hào)進(jìn)行分解，再利用AIC 分析各IMF 的初至?xí)r間。

VMD 分解有兩個(gè)重要參數(shù)需要預(yù)先設(shè)定，一個(gè)是分解層數(shù)L，一個(gè)是用于確保準(zhǔn)確重構(gòu)信號(hào)的懲罰因子α。隨著分解層數(shù)L的不同，VMD 會(huì)出現(xiàn)不同分解結(jié)果，L過大會(huì)使信號(hào)斷斷續(xù)續(xù)無規(guī)律，過小則易導(dǎo)致模態(tài)混疊，直接影響分解的準(zhǔn)確性，因此需要選擇合理的分解層數(shù)。本文擬對(duì)VMD 進(jìn)行優(yōu)化，根據(jù)分解后IMF 的特點(diǎn)以及IMF 與原信號(hào)的相關(guān)性自適應(yīng)確定分解層數(shù)，使分解后的IMF 不混疊、較純凈，能精細(xì)刻畫原信號(hào)的時(shí)頻特性，從而進(jìn)行初至識(shí)別。

2 改進(jìn)STA/LTA 及優(yōu)化VMD 的峰度AIC 微震初至拾取

2.1 算法基本思想

針對(duì)信噪比低導(dǎo)致P 波拾取精確度低的問題，本文對(duì)傳統(tǒng)STA/LTA 算法進(jìn)行改進(jìn)，基于地震信號(hào)幅度的實(shí)時(shí)變化引入權(quán)重因子，可快速、靈敏地從微震監(jiān)測數(shù)據(jù)中提取微地震事件，實(shí)現(xiàn)P 波的粗略提取；然后選取初次拾取時(shí)間點(diǎn)前后2—3 s 的記錄，根據(jù)互相關(guān)系數(shù)和排列熵準(zhǔn)則確定VMD 分解層數(shù)，使分解后的IMF 無偽分量和模態(tài)混疊，能較好地表征微震信號(hào)特征；最后計(jì)算各IMF 分量的能量比，并進(jìn)行峰度AIC 識(shí)別，以各IMF 分量的拾取結(jié)果及能量比進(jìn)行綜合加權(quán)，得到精確的P 波到時(shí)，算法流程圖如圖1 所示。

圖1 基于改進(jìn)STA/LTA 及優(yōu)化VMD 的微震初至拾取算法流程圖Fig. 1 The flow chart of the first break picking of micro-seismic based on improved STA/LTA and adaptive VMD

2.2 基于改進(jìn)STA/LTA 的初步拾取

特征函數(shù)的選取決定了初至拾取的精度，因此要構(gòu)建盡可能靈敏反映信號(hào)頻率或振幅變化的特征函數(shù)，使微震到來時(shí)，根據(jù)特征函數(shù)計(jì)算所得的STA 發(fā)生明顯變化，STA/LTA 出現(xiàn)較大幅度的增大。本文根據(jù)信號(hào)幅度的實(shí)時(shí)變化引入加權(quán)系數(shù)K，其計(jì)算方法為

式中y為地震信號(hào)序列。當(dāng)微震事件未到來時(shí)，信號(hào)幅度變化較小，K較小（接近于0），而當(dāng)微震到來，信號(hào)幅度發(fā)生較大變化時(shí)，K將增大，可見K能實(shí)時(shí)反映信號(hào)的幅度變化。利用K構(gòu)造新的特征函數(shù)，即

y（i）2反映信號(hào)幅度特性，［y（i）－y（i－1）］2反映信號(hào)頻率特性。K能根據(jù)信號(hào)采樣頻率對(duì)信號(hào)幅度和頻率進(jìn)行加權(quán)分配，K的引入將增加特征函數(shù)的敏感度。當(dāng)微震事件到來時(shí)，K較大，對(duì)特征函數(shù)有放大作用，使特征函數(shù)隨之發(fā)生較大變化，STA/LTA 可迅速反映出信號(hào)幅度或頻率的變化，從而更易識(shí)別微震事件。相比劉曉明等（2017）的研究中權(quán)重因子不能根據(jù)信號(hào)幅度變化靈活調(diào)整的不足，本文特征函數(shù)能更好地反映并放大地震信號(hào)的實(shí)時(shí)變化，從而提高微震信號(hào)檢測的靈敏度。

改進(jìn)的STA/LTA 可快速可靠地初步識(shí)別微震，為了提高識(shí)別精度，取初次拾取時(shí)間點(diǎn)前后2—3 s 的記錄進(jìn)行二次拾取。由于VMD 算法能較好地描述信號(hào)的局部突變或漸變特征，本文對(duì)VMD 算法進(jìn)行優(yōu)化，對(duì)分解后的IMF 進(jìn)行基于峰度-AIC 的二次拾取。

2.3 基于互相關(guān)系數(shù)和排列熵優(yōu)化的VMD

經(jīng)VMD 后，地震信號(hào)s（t）可表示為

式中，sIMF（t）為從高頻到低頻的IMF 分量，rn（t）為殘差分量。為確保分解后IMF 既不包含偽分量，又能較好地表征原信號(hào)的特征，本文根據(jù)互相關(guān)系數(shù)和排列熵準(zhǔn)則來確定VMD 分解層數(shù)。互相關(guān)系數(shù)的計(jì)算方法為

互相關(guān)系數(shù)衡量的是兩個(gè)信號(hào)的相關(guān)程度，根據(jù)互相關(guān)系數(shù)理論，當(dāng)兩個(gè)信號(hào)間的互相關(guān)系數(shù)大于0.3 時(shí)，表明相關(guān)性較大（鄭近德等，2013），計(jì)算分解后IMF 與原信號(hào)的互相關(guān)系數(shù)，若互相關(guān)系數(shù)不小于0.3，則表明不含偽分量，反之則表明出現(xiàn)了過分解。

為了進(jìn)一步確定分解層數(shù)，本文利用排列熵（Bandt，Pompe，2002）進(jìn)行分解后IMF 的隨機(jī)性和異常檢測。排列熵能較好地反映時(shí)間序列細(xì)節(jié)的微小變化，有效進(jìn)行異?；蛲蛔儥z測。其計(jì)算方法為：

① 設(shè)原時(shí)間序列為x（i）（i＝1，2，···，N），利用相空間重構(gòu)得到原序列的重構(gòu)矩陣，即

排列熵通過時(shí)間序列的相鄰數(shù)據(jù)進(jìn)行對(duì)比而不是考慮數(shù)據(jù)的具體值來度量序列隨機(jī)性，其變化反映并放大了原序列x（i）的微弱變化，從而易于檢測信號(hào)異常。對(duì)歸一化后的排列熵，其值越趨近于1，表明信號(hào)越隨機(jī)，據(jù)此可判斷分解后IMF 是否為異常信號(hào)。根據(jù)王志堅(jiān)等（2018）的研究，分解到某一層后的IMF 分量排列熵不大于0.6，則表明分解后的IMF 無異常分量，能較好地刻畫原信號(hào)的特征，否則表明出現(xiàn)過分解。基于互相關(guān)系數(shù)和排列熵準(zhǔn)則優(yōu)化VMD，確定分解層數(shù)L的算法步驟為：① 設(shè)定初始分解層數(shù)i＝2；② 基于VMD 對(duì)原信號(hào)進(jìn)行分解，得到一系列IMFj（j＝1，2，···，i），計(jì)算IMFj與原信號(hào)的互相關(guān)系數(shù)corj，以及IMFj的排列熵Hp；③ 若corj≤0.3，且Hp≥0.6，則表明出現(xiàn)過分解，存在異常分量，則停止分解，分解層數(shù)L＝i－1，否則令i＝i＋1 進(jìn)入第2 步繼續(xù)VMD 分解。經(jīng)優(yōu)化VMD 分解后的IMF 既不包含異常分量，又能刻畫原信號(hào)的起伏特征和細(xì)節(jié)信息。

2.4 基于優(yōu)化VMD 的峰度AIC 精確識(shí)別P 波初至

為了提高拾取精度，對(duì)改進(jìn)STA/LTA 拾取的P 波到時(shí)tP前后2—3 s 記錄進(jìn)行二次拾取，以獲得精確的P 波到時(shí)。根據(jù)互相關(guān)系數(shù)和排列熵確定分解層數(shù)，對(duì)這段記錄進(jìn)行優(yōu)化VMD 分解，實(shí)現(xiàn)信號(hào)分解。如前所示，VMD 中參數(shù)懲罰因子α需提前設(shè)定，一般為采樣頻率的0.25—2 倍。α越小，分解后IMF 分量帶寬越大，反之分解后IMF 分量帶寬越小 （唐貴基，王曉龍，2015）。鑒于峰度AIC 的穩(wěn)定性和精度高，對(duì)分解出的L個(gè)IMF 根據(jù)式（4）計(jì)算峰度AIC 值，以最小值為IMFi（i≤L）的拾取時(shí)間tAi。VMD 分解后IMF 分量的中心頻率由低至高分布，且具有不同的能量，其中低頻IMF 是包含有效微震信號(hào)自身特征信息的主模態(tài)分量，占原信號(hào)總能量中的主要部分。因此，對(duì)L個(gè)IMF 計(jì)算能量比，其定義為

根據(jù)各IMF 的拾取時(shí)間及計(jì)算所得能量比，按式（11）進(jìn)行加權(quán)，得到最終的P 波拾取到時(shí)tAi，即

3 實(shí)驗(yàn)結(jié)果及分析

3.1 改進(jìn)STA/LTA 的模擬微震拾取

設(shè)雷克（Ricker）子波主頻為20 Hz，采樣頻率為1 000 Hz，基于卷積模型得到的理論地震記錄如圖2 所示，其中加入不同SNR 的隨機(jī)噪聲后可得模擬地震記錄，當(dāng)SNR 為-5 dB 時(shí)，由于噪聲嚴(yán)重，起震不明顯，且有效微震信號(hào)被噪聲嚴(yán)重干擾，增加拾取難度。

圖2 基于雷克子波和卷積模型的理論和模擬地震記錄Fig. 2 Theoretical and simulated seismogram based on Ricker wavelet and convolution model

設(shè)STA＝0.1 s，LTA＝0.5 s，微震閾值為1.5，分別采用本文改進(jìn)STA/LTA，傳統(tǒng)STA/LTA，劉曉明等（2017）的改進(jìn)STA/LTA 進(jìn)行初次識(shí)別，結(jié)果如圖3 所示，人工拾取時(shí)間為3.065 s，本文改進(jìn)STA/LTA 識(shí)別時(shí)間為3.125 s，傳統(tǒng)STA/LTA 識(shí)別時(shí)間為3.218 s，劉曉明等的改進(jìn)S T A/L T A 拾取時(shí)間為3.203 s?？梢?，STA/LTA 算法拾取時(shí)間有延遲，相比傳統(tǒng)STA/LTA，劉曉明等的改進(jìn)STA/LTA 和本文的改進(jìn)STA/LTA 算法能更快地反映信號(hào)變化，從而快速檢測微震事件。相比劉曉明等的改進(jìn)SAT/LTA，本文改進(jìn)STA/LTA 權(quán)重因子K能實(shí)時(shí)反映信號(hào)幅度變化，而非固定值，特征函數(shù)能更快速、真實(shí)地反映信號(hào)變化，可快速增加超過閾值，從而能快速識(shí)別微震。以人工拾取結(jié)果為參考，本文改進(jìn)STA/LTA 拾取誤差為0.06 s，劉曉明等的改進(jìn)STA/LTA 和傳統(tǒng)STA/LTA 的誤差分別為0.138 s 和0.153 s，可見本文改進(jìn)STA/LTA 算法能減小初次拾取誤差。

圖3 三種STA/LTA 算法拾取結(jié)果Fig. 3 Picking up results of three kinds of STA/LTA algorithms

為了驗(yàn)證改進(jìn)STA/LTA 的性能及抗噪性，將此理論地震記錄隨機(jī)添加不同SNR 噪聲進(jìn)行拾取（SNR 為-5—20 dB），并隨機(jī)仿真100 次，計(jì)算與人工拾取誤差的均值。圖4 給出了3 種STA/LTA 算法在不同SNR 下的拾取誤差曲線，可見，SNR 越低，拾取誤差越大。但在較低SNR 時(shí)，本文改進(jìn)STA/LTA 拾取誤差更小，在5—-5 dB 時(shí)拾取結(jié)果與人工拾取的誤差為0.05—0.28 s，比傳統(tǒng)STA/LTA 和劉曉明等的改進(jìn)STA/LTA 算法的拾取誤小差約0.02—0.1 s；當(dāng)SNR 在5—-20 dB 時(shí)，本文與人工拾取的誤差為0.02—0.05 s，比傳統(tǒng)STA/LTA 和劉曉明等的改進(jìn)STA/LTA 算法的拾取誤差小約0.01—0.02 s，以上說明本文的改進(jìn)STA/LTA 有效，比傳統(tǒng)STA/LTA 方法降低誤差0.01 s 以上，適用于強(qiáng)噪聲微震初至拾取，且性能較好，尤其是在低信噪比時(shí)，能有效降低初次拾取誤差。

圖4 不同SNR 下三種STA/LTA 算法拾取誤差曲線Fig. 4 Picking up error curves of three STA/LTA algorithms under different SNRs

3.2 優(yōu)化VMD 峰度AIC 的模擬微震拾取

為了提高拾取精度，本文利用優(yōu)化VMD對(duì)信號(hào)進(jìn)行分解，并利用峰度AIC 再次識(shí)別。根據(jù)優(yōu)化VMD 確定分解層數(shù)，并將初次拾取時(shí)間向前后各推2—3 s，精確確定該記錄P 波到時(shí)。圖2 中的模擬地震記錄，經(jīng)優(yōu)化VMD 確定的分解層數(shù)為2，得到IMF1 和IMF2。由圖5 可見，分解后的IMF 不包含噪聲異常分量，一定程度上消除噪聲影響，且不同頻率段的各IMF 能較準(zhǔn)確地反映原信號(hào)中不同成分的振蕩信息，保持信號(hào)細(xì)節(jié)特征。根據(jù)式（4）計(jì)算對(duì)應(yīng)的峰度AIC 曲線，以曲線最小值點(diǎn)為初至?xí)r間，圖5 給出了峰度AIC 曲線拾取時(shí)間分別為3.06 s 和3.087 s 的IMF1 和IMF2，按式（10）計(jì)算出的IMF1 和IMF2 的能量比分別為0.55，0.45，最終加權(quán)拾取時(shí)間為3.07 s，與人工拾取時(shí)間3.065 s 的誤差為0.005 s，可見優(yōu)化VMD 峰度AIC 能有效提高拾取精度。

圖5 優(yōu)化VMD 峰度AIC 為3.07 s 時(shí)的初至拾取結(jié)果Fig. 5 First arrival picking based on improved VMD when Kurtosis-AIC is 3.07 s

為了驗(yàn)證算法性能，將小波包峰度AIC （田優(yōu)平，趙愛華，2016）、EMD-AIC 算法（Liuet al，2017）與本文優(yōu)化VMD 峰度AIC 算法進(jìn)行對(duì)比。如圖6，7 所示，對(duì)于圖5 的原信號(hào)地震記錄，基于EMD-AIC 的初至拾取時(shí)間為3.077 s （3 個(gè)IMF 的拾取時(shí)間分別為3.098 s，3.052 s，3.081 s）；小波包峰度AIC 的初至拾取時(shí)間為3.075 s （三個(gè)分量的拾取時(shí)間分別為3.073 s，3.078 s，3.075 s），與參考結(jié)果的誤差分別為0.012 s 和0.01 s?？梢?，EMD-AIC 算法、小波包峰度AIC 和優(yōu)化VMD 峰度AIC 算法均能較準(zhǔn)確地拾取到時(shí)，且本文算法拾取誤差更小，精度較高。

圖6 基于EMD-AIC 算法的初至拾取Fig. 6 First arrival picking based on EMD-AIC

圖7 基于小波包峰度AIC 的初至拾取Fig. 7 First arrival picking based on wavelet packet Kurtosis-AIC

為了分析本文優(yōu)化VMD 峰度AIC 算法在降低二次拾取誤差方面的性能，取不同SNR（-5—20 dB）的模擬地震信號(hào)，同一SNR 添加隨機(jī)噪聲仿真100 次，并對(duì)比初次與二次拾取時(shí)間與人工結(jié)果的誤差均值。如圖8 所示，當(dāng)SNR 較低（-5—3 dB）時(shí)，初次拾取誤差偏大，約為0.05—0.28 s；經(jīng)二次拾取后，誤差明顯降低，最終結(jié)果與人工拾取結(jié)果平均相差在0.023 s以內(nèi)；隨著SNR 的提高，初次拾取誤差逐漸減小且再次拾取誤差大大降低，與人工拾取結(jié)果非常接近，誤差在0.01 s 以內(nèi)。

圖8 初次與二次P 波拾取的誤差對(duì)比Fig. 8 Comparison of initial and second P-wave picking errors

為了驗(yàn)證VMD 峰度AIC 算法在降低二次拾取誤差方面的性能，添加不同SNR 的噪聲并隨機(jī)仿真100 次取均值，對(duì)比小波包峰度AIC 和EMD-AIC 算法與人工拾取結(jié)果的誤差。如圖9 所示，三種算法都能較好地拾取到初至P 波，本文算法在低SNR 下表現(xiàn)穩(wěn)定，平均誤差在0.023 s 以內(nèi)。相比小波包分解和EMD，經(jīng)優(yōu)化VMD 后IMF 分量能更好地保留原信號(hào)的細(xì)節(jié)和突變特征，且有效消除噪聲影響，經(jīng)峰度AIC 拾取后能有效降低拾取誤差。

圖9 3 種不同算法的拾取誤差對(duì)比Fig. 9 Picking errors of three different algorithms

3.3 實(shí)際微震記錄拾取

為測試本文算法拾取真實(shí)微震記錄的準(zhǔn)確性，對(duì)2013—2019 年間廣西平果、武鳴、貴港、樂業(yè)等地的若干次ML1.3—2.0 的礦藏勘探、開采等人工爆破微震記錄共計(jì)1 137 條進(jìn)行拾取，圖10 給出了平果地區(qū)2013 年7 月4 日15 時(shí)3 分在（23.343°N，107.583°E）發(fā)生的ML＝1.3 爆破的某條監(jiān)測記錄，其采樣率為100 Hz，爆破持續(xù)時(shí)間較短，且有較多外部環(huán)境噪聲，SNR 較低，該記錄人工拾取的參考時(shí)間為14.76 s。本文改進(jìn)STA/LTA 初次拾取時(shí)間為14.78 s，傳統(tǒng)STA/LTA 拾取時(shí)間為14.81 s，劉曉明等的改進(jìn)STA/LTA 拾取時(shí)間為14.79 s（圖11）?？梢姡瑯拥拈撝?，本文改進(jìn)STA/LTA 更敏感快速反映信號(hào)變化，減少初次拾取誤差。

圖10 廣西平果2013 年7 月4 日ML1.3 爆破的一條實(shí)際微震記錄Fig. 10 Real micro-seismic record of some coal mine with ML1.3 at Pingguo，Guangxi on July 4th 2013

圖11 本文改進(jìn)STA/LTA 初次拾取結(jié)果Fig. 11 First picking results of improved STA/LTA in this paper

在二次拾取中，本文優(yōu)化VMD 自適應(yīng)確定的分解層數(shù)為3 （圖12），各分量能較好地刻畫原地震信號(hào)的細(xì)節(jié)特征，各分量拾取時(shí)間分別為14.8 s，14.73 s 和14.77 s，能量比為0.31，059，0.1，則最終拾取時(shí)間為14.76 s，與人工拾取結(jié)果一致。

圖12 本文改進(jìn)STA/LTA 及優(yōu)化VMD峰度AIC 拾取結(jié)果Fig. 12 First arrival picking by improved STA/LTA and adaptive VMD Kurtosis-AIC

EMD-AIC 3 個(gè)IMF 分量的拾取時(shí)間分別如圖13 所示，為14.76 s，14.76 s，14.83 s，最終拾取時(shí)間為14.78 s，誤差為0.02 s。小波包峰度AIC 算法分解的3 個(gè)IMF 拾取時(shí)間分別為圖14 所示的14.78 s，14.76 s，14.76 s，最終拾取時(shí)間為14.77 s，誤差為0.01 s。3 種算法都能較準(zhǔn)確拾取初至波，誤差在0.02 s內(nèi)，但本算法誤差更小。

圖13 EMD-AIC 拾取結(jié)果Fig. 13 Picking results of EMD-AIC

圖14 小波包峰度AIC 的拾取結(jié)果Fig. 14 Picking results of wavelet packet Kurtosis-AIC

對(duì)其余記錄分別用本算法、EMD-AIC，小波包峰度AIC 算法進(jìn)行拾取，以人工拾取結(jié)果為參考，3 種算法的拾取結(jié)果列于表1。

表1 實(shí)際微震記錄初至拾取結(jié)果Table 1 First arrival picking of real micro-seismics

本文改進(jìn)STA/LTA 及優(yōu)化VMD 峰度AIC 算法拾取誤差在30，20，10 ms 的比例分別為98.2%，95.9%，90.7%，將誤差在0.02 s （20 ms）以內(nèi)的視為準(zhǔn)確拾取，結(jié)果顯示，本文算法有效，拾取準(zhǔn)確率為95.9%，分別高于EMD-AIC 和小波包峰度AIC 約3.4%和4.7%，且本文拾取誤差在10 ms 的比例比其它兩種算法的要高4.5%以上，表明本文算法拾取誤差更小，具有更高的精度。

4 討論與結(jié)論

本文對(duì)傳統(tǒng)STA/LTA 的特征函數(shù)進(jìn)行了改進(jìn)，提高P 波的初次拾取精度；再對(duì)VMD 進(jìn)行了優(yōu)化，能自適應(yīng)確定分解層數(shù)，對(duì)分解后IMF 基于峰度AIC 二次拾取后，基于各IMF 的能量比綜合加權(quán)得到最終二次拾取到時(shí)。實(shí)驗(yàn)表明：

1） 改進(jìn)STA/LTA 利用信號(hào)幅度設(shè)計(jì)權(quán)重因子構(gòu)建特征函數(shù)，能快速反映信號(hào)變化，提高初至P 波的初次檢測靈敏度；

2） 優(yōu)化多尺度VMD 分解，能根據(jù)分解IMF 分量的互相關(guān)系數(shù)和排列熵確定VMD 分解層數(shù)，既能解決模態(tài)混疊，又能防止出現(xiàn)過分解，且充分保留微震信號(hào)的細(xì)節(jié)特征；

3） 仿真實(shí)驗(yàn)表明，在較低信噪比時(shí),算法仍能較好識(shí)別初至P 波，比傳統(tǒng)STA/LTA 和劉曉明等（2017）改進(jìn)STA/LTA 減小初次識(shí)別誤差約0.01 s 以上；基于優(yōu)化VMD 分解后，能再次降低拾取誤差，相比小波包峰度AIC 和EMD-AIC 拾取算法，本文最終誤差基本在0.023 s以內(nèi)，表明本文算法的有效性、準(zhǔn)確性和抗噪性；

4） 實(shí)際微震拾取表明，本文拾取結(jié)果與人工拾取誤差在20 ms 內(nèi)的記錄占比95.9%，準(zhǔn)確率高于EMD-AIC 和小波包峰度AIC，表現(xiàn)出較強(qiáng)性能。