賈瑞生,譚云亮,孫紅梅,洪永發(fā)

(1.山東科技大學礦山災害預防控制省部共建國家重點實驗室培育基地,山東青島 266590;2.山東科技大學信息科學與工程學院,山東青島266590;3.山東科技大學礦業(yè)與安全工程學院,山東青島 266590)

?

低信噪比微震P波震相初至自動拾取方法

賈瑞生1,2,譚云亮1,3,孫紅梅2,洪永發(fā)2

(1.山東科技大學礦山災害預防控制省部共建國家重點實驗室培育基地,山東青島 266590;2.山東科技大學信息科學與工程學院,山東青島
266590;3.山東科技大學礦業(yè)與安全工程學院,山東青島 266590)

摘 要:針對微震信號固有的低信噪比、非平穩(wěn)性、隨機性等特征,發(fā)展了一種基于希爾伯特-黃變換(HHT)和AIC準則相結合的震相初至自動拾取方法,該方法首先對含噪信號進行經驗模態(tài)分解(EMD)及內蘊模態(tài)函數(shù)(IMF)重構,在保留微震信號固有特征基礎上實現(xiàn)降噪;其次,基于Hilbert變換計算出歸一化包絡信號,通過設置包絡閾值搜索震相初至的大致位置,并以該位置為基礎為AIC函數(shù)選擇合適的計算時窗以降低計算復雜度;最后,在選擇的時窗內應用AIC準則計算出P波震相初至。應用該方法對隨機選取的高、低信噪比兩類共80個微震信號進行震相初至拾取,若以人工拾取結果為基準、時差在10 ms以內視為準確拾取,結果顯示,高信噪比信號拾取準確率為100%,低信噪比信號拾取準確率為92%。

關鍵詞:微震P波;震相初至拾取;經驗模態(tài)分解;希爾伯特-黃變換;AIC準則

責任編輯:常 琛

賈瑞生,譚云亮,孫紅梅,等.低信噪比微震P波震相初至自動拾取方法[J].煤炭學報,2015,40(8):1845-1852.doi:10.13225/ j.cnki.jccs.2014.1122

微震監(jiān)測是預警沖擊地壓、煤與瓦斯突出、礦井突水等煤礦災害事故的主要手段之一,目前已在微震信號頻譜特征提取、災害危險預警及防治等方面取得許多研究成果[1],這些成果的取得均建立在震源準確定位的基礎之上,而震源準確定位技術的關鍵環(huán)節(jié)之一則是微震P波初至到時的準確拾取[2-3],因此微震P波震相初至自動拾取是微震監(jiān)測系統(tǒng)的關鍵技術之一,是震源定位及震源機制解釋的前提與基礎。

微震監(jiān)測系統(tǒng)拾取的震動信號具有瞬態(tài)性、多樣性及不確定性等特點,信號中夾雜有機械震動、環(huán)境噪聲、電磁噪聲等多種干擾信號,因此,如何在低信噪比信號中準確辨識微震P波到時十分困難。對于地震事件,目前已發(fā)展多種震相自動識別及拾取方法, 如Allen等基于時間域信號能量變化構建特征函數(shù),提出了長短時均值比方法(STA/ LTA)[4-6];Akaike等基于地震波可劃分為局部平穩(wěn)段的假設,依據(jù)到達前后地震波形數(shù)據(jù)統(tǒng)計的差別,提出了判別震相初至的AIC準則[7-11];Saragiotis等基于地震波形的偏斜度和峰度提出了PAI-S/ K方法[12];常旭等應用分形理論對地震記錄進行了有效的初至自動拾取[13];馬強等綜合應用STA/ LTA方法及AIC準則,基于Delaunay三角剖分剔除干擾信號,提出了多步驟的P波自動拾取方法[14];王繼等應用單臺Akaike信息準則和多臺AIC最小二乘互相關方法,發(fā)展了震相自動精確檢測技術,實現(xiàn)了流動地震臺陣觀測震相初至的自動拾取[15];劉勁松等通過分析STA/ LTA,AIC,PAI-S/ K等幾種方法的原理及特點,提出了移動時窗峰度的快速算法和改進的峰度拾取初至算法[16]。以上方法均已在不同場合得到實際應用,能夠實現(xiàn)對較高信噪比地震波的到時拾取,但當信號信噪比較低時可能會出現(xiàn)拾取結果精度不夠甚至錯誤。對于微地震事件,其P波震相到時自動拾取的特殊性在于:相較于天然地震而言,微震信號頻率較高、信噪比低,自動拾取震相初至易受外部噪聲干擾。

為此,本文從信號分析理論入手,引入希爾伯特-黃變換對含噪微震信號進行降噪處理及包絡分析,通過設置包絡閾值為AIC算法選擇合適的計算時窗,再在時窗內應用AIC方法,實現(xiàn)了低信噪比微震P波震相初至的自動拾取。

1 微震P波震相自動拾取算法分析

1.1 長短時均值比(STA/ LTA)法

設x(t)為t時刻P波垂向幅值,f(t)為高通濾波后的垂向幅值數(shù)據(jù),則有

其中,α為高通濾波因子,由于P波震相在垂直方向幅值較大,可以通過構造信號特征函數(shù)來放大垂直方向幅值,常用的P波拾取特征函數(shù)為

STA/ LTA方法的基本原理是根據(jù)微震P波特征函數(shù)值的長短時均值比判斷其初至[8-10],應用STA/ LTA的基本形式如下:

其中,t0為當前時刻;t1,t2為t0之前某時刻,且滿足t2

圖1　STA/ LTA法拾取P波到時Fig.1 Pick up P-arrival time based on STA/ LTA

STA/ LTA方法的突出優(yōu)點是算法穩(wěn)定可靠、拾取效率高,但其觸發(fā)點一般滯后于實際初至點(圖1),當微震信號中存在干擾信號時,通常會產生震相初至點的錯誤判斷[14],因而適用于高信噪比微震信號的震相初至拾取。

1.2 AIC方法

AIC方法的基本原理是求取地震信號AIC函數(shù)的局部最小值,Sleema提出了AR-AIC準則[10],它根據(jù)自回歸過程將地震波形數(shù)據(jù)分成2個局部統(tǒng)計時段(圖2),AR-AIC函數(shù)表示為

其中,k為兩個局部統(tǒng)計時段分界點;p為AR過程階數(shù);l為地震波形數(shù)據(jù)長度;分別為2個局部統(tǒng)計時段的擬合誤差;C為一個常數(shù)。為了求出震相初至,必須求出該函數(shù)中AR模型的階數(shù)和系數(shù),該方法計算復雜度較高,不利于震相初至的實時拾取。

圖2　AR-AIC法拾取P波初至Fig.2 Pick up P-arrival time based on AR-AIC

不同于AR-AIC模型,Maeda提出直接由地震波形數(shù)據(jù)計算AIC函數(shù)[11],求取AIC函數(shù)的局部最小值(圖3),該值對應的位置即為震相初至,AIC函數(shù)表示為

其中,x(i)(i=1,2,…,l)為地震波形離散數(shù)據(jù);k的取值范圍是數(shù)據(jù)窗口內所有采樣點。

圖3　微震波形和時窗[500,2 500]內的AIC函數(shù)曲線Fig.3 A microseismic signal and AIC function curve whentime window between 500 ms to 2 500 ms

對比以上2種AIC方法,后者不需要計算AR模型的階數(shù)即可直接求取AIC值,在震相初至拾取實時性要求很高的情況下,是一種更加高效實用的算法,但是該算法需要在震相初至的附近尋找一個合適的時窗來計算AIC值,這是因為不同的時窗可能使AIC函數(shù)局部最小值出現(xiàn)的位置不同。圖3(a),4(a)為同一個微震波形數(shù)據(jù)在不同時窗內的波形,圖3(b),4(b)為對應的AIC函數(shù)曲線,其中圖4由于時窗設置不合理導致震相初至拾取錯誤,如何合理選擇時窗是AIC法準確拾取震相初至的關鍵問題之一。

圖4　微震波形和時窗[1 000,3 000]內的AIC函數(shù)曲線Fig.4 A micro-seismic signal and AIC function figure whentime window between 1 000 ms to 3 000 ms

此外,微震監(jiān)測過程中拾取的信號易受外界噪聲干擾,使得AIC法震相初至拾取的精度大打折扣,有時還會造成微震事件的誤拾取[10-11],因此,如何降低環(huán)境噪聲對震相初至拾取的影響成為目前該領域研究的關鍵問題。由于微震信號及外部噪聲具有隨機性、非平穩(wěn)性等特征,常規(guī)降噪方法(如Fourier方法)無法有效消除微震信號中夾雜的外部噪聲,而HHT則能對非平穩(wěn)信號進行EMD分解并重構,在保留信號固有非平穩(wěn)特征基礎上實現(xiàn)降噪。

2 對震相自動拾取算法的改進

2.1 EMD變換及信號降噪處理

Huang等提出了分析處理非平穩(wěn)信號的EMD方法[17],即Huang變換,其核心是:任何非平穩(wěn)信號都是由一些互不相同的、非正弦的、簡單的IMF函數(shù)組成的,因此可將信號分解成從高頻到低頻若干個基本時間序列的組合,每個基本時間序列即是一個IMF分量。基于HHT的核心思想,含噪微震信號的EMD分解及消噪算法如下:

(1)確定微震信號x(t)時間曲線上所有局部極大、極小值點,分別將其用3次樣條函數(shù)擬合為原數(shù)據(jù)序列的上、下包絡線,上、下包絡線對應的坐標均值連接成均值線m1(t),由式(8)計算出h1(t),即

(2)判斷h1(t)是否滿足IMF函數(shù)所需條件,令

當SD值介于[0.2,0.3]之間時,篩選過程終止,得到信號x(t)的第1個IMF分量;若SD值不在[0.2,0.3]之間,此時把h1(t)看作新的時間序列曲線,重復以上步驟繼續(xù)分解,則有

經過k次迭代后,當SD值介于[0.2,0.3]之間時,就得到了第1階IMF分量h1k(t),將其記為c1(t),它表示微震信號的最高頻部分。

(3)從原信號x(t)中減去該信號的高頻部分c1(t),得到頻率較低的信號殘差r1(t):

(4)將r1(t)看作一個新的信號序列,重復以上步驟,得到一系列的ci(t)和最后一個不可分解的信號殘差r(t),則原信號可表達如下:

(5)剔除高頻IMF分量,對剩余的IMF分量進行重構,得到降噪后的微震信號:

圖5　一個含噪微震信號的EMD分解Fig.5 Empirical mode decomposition of a micro seismic signal with noise

圖5顯示,含噪信號經EMD分解后得到7個IMF分量,這些IMF分量從高頻到低頻依次排列,由于微震頻率范圍一般介于0~200 Hz之間,觀察IMF分量的頻譜分布,可以判斷出imf1為高頻噪聲成分,故可將其剔除,重構剩余的IMF即得到降噪后的微震有效信號。圖6為EMD降噪前后微震信號的波形及頻譜,對比圖6(a),(c)可知,經過EMD降噪后較好地保留了信號的尖峰、突變及波形特征,說明EMD法能充分保留信號的瞬態(tài)非平穩(wěn)特征;對比圖6(b),(d)可知,圖6(d)可清晰辨別出微震信號的時頻特征,無用高頻干擾噪聲被有效濾除,為進一步識別震相初至奠定基礎。

2.2 Hilbert變換及震相初至拾取算法

Hilbert變換是非平穩(wěn)信號分析與處理的有效方法,實信號x(t)的Hilbert變換可定義[18]如下

Hilbert變換即是信號和時間倒數(shù)的卷積,實信號x(t)和它的Hilbert變換可組成解析信號:

其中, a(t)為解析信號的幅值;φ(t)為解析信號的相位,分別表示為

式中,a(t)為實信號x(t)的包絡。

將降噪后的微震信號進行Hilbert變換并進行包絡分析,即按時序將上下包絡的絕對值相加并將其歸一化處理,則圖6(c)的Hilbert包絡信號波形如圖7所示,對同一類微震信號使用同一閾值(本文設為0.3)沿時間軸順序搜索,當包絡值大于設定閾值時,即可得到P波震相初至的大致位置t0,本例中t0= 1 243 ms。在經過HHT降噪后的微震信號波形上以t0為基準向前及向后分別取500個采樣點作為計算時窗,在該時窗內應用AIC函數(shù)求解P波震相初至,求解結果如圖8所示。

圖6　去噪前后的微震信號波形及其頻譜Fig.6 Microseismic signal waveform and its spectrum before and after de-noising

圖7　Hilbert包絡信號波形Fig.7 Waveform of Hilbert envelope signal

圖8　時窗[744,1 743]內的微震波形和對應AIC函數(shù)曲線Fig.8 Microseismic signal when time window between 744 msand 1 743 ms and the corresponding AIC function curve

綜上,基于HHT及AIC求解震相初至的改進算法(簡稱HHT-AIC法)如下:①對含噪微震信號進行EMD分解,按照頻率從高到低獲得一系列IMF分量,根據(jù)外部噪聲特征對IMF分量進行剔除,并把剩余的IMF分量進行重構,實現(xiàn)微震信號的降噪;②對降噪后的微震信號進行Hilbert變換,求出其包絡信號,并對其進行歸一化處理;③設置包絡閾值,沿時間軸查找第1個大于包絡閾值的時刻,該時刻即為震相初至的粗略值,以該時刻為基準分別向前、向后各取500個采樣點作為計算時窗;④在選取的時窗內應用式(7)計算AIC函數(shù)值,AIC函數(shù)局部最小值對應的時刻即為震相初至。

3 算法檢驗及討論

3.1 實驗數(shù)據(jù)來源

為檢驗算法的有效性,我們分別選取一組高信噪比(標記為H組)及一組低信噪比(標記為L組)微震信號進行對比試驗,其中H組信號來源于我國西部某水電站微震監(jiān)測數(shù)據(jù),該水電站所在地主要為硬巖結構,巖石破裂產生的微震信號起跳明顯(圖9(a)),監(jiān)測環(huán)境較為理想,信號中噪聲干擾較少,從中隨機抽取30組微震信號用于實驗;L組微震信號來源于我國西北某煤礦井下采場,屬軟巖地質構造,且機械振動、運輸、放炮等產生的外部噪聲較多,所采集的微震信號信噪比低,波形起跳不明顯(圖9(b)),從中隨機抽取50組微震信號用于實驗。

圖9　高信噪比和低信噪比微震波形Fig.9 Micro seismic waveform with high signal to noise ratio and Micro seismic waveform with low signal to noise ratio

分別應用AIC方法及本文方法對上述2組不同類型的微震信號進行處理,并將其與人工拾取的震相初至結果對比,從拾取結果的正確率、算法耗時等方面進行評價。

3.2 實驗結果及分析

應用Matlab編程實現(xiàn)AIC算法及HHT-AIC算法,對圖9(a)所示的高信噪比微震波形進行震相初至拾取,與人工拾取的結果對比如圖10所示。

圖10所示的AIC及HHT-AIC兩種拾取方法都獲得了同樣的震相初至時刻739 ms,在相同軟硬件配置的PC機上,AIC法運算耗時0.641 s,HHT-AIC法耗時1.892 s,AIC法耗時最少,但AIC法在選擇時窗大小及位置不合理時可能會出現(xiàn)圖4(b)出現(xiàn)的錯誤拾取結果,故AIC法是不穩(wěn)定的;HHT-AIC法包括EMD降噪(耗時1.689 s)、Hilbert包絡分析及AIC函數(shù)計算(耗時0.203 s)等過程,計算復雜度較高,但當處理高信噪比信號時可以省去EMD降噪過程,僅保留Hilbert包絡分析及AIC函數(shù)計算過程(以下稱HT-AIC法),此時算法僅耗時0.203 s,大大低于AIC法的計算耗時,且能保持震相初至拾取精度不變,這是因為HHT-AIC法利用Hilbert變換進行包絡分析時為AIC函數(shù)選擇了合適的計算時窗(圖10(c)),時窗內僅有1 000個采樣點數(shù)據(jù),大大降低了AIC函數(shù)的計算強度。

對H組中30個高信噪比微震信號進行處理,以人工拾取的震相初至作為參照,震相初至拾取結果見表1。其中AIC法的計算時窗為P波初至前750個至P波初至后2 250個采樣點數(shù)據(jù),時窗內共3 000個采樣點數(shù)據(jù);HHT-AIC時窗為P波初至前后各500個采樣點數(shù)據(jù),時窗內共1 000個采樣點數(shù)據(jù); HT-AIC法所采用的時窗位置及大小與HHT-AIC方法所采用的時窗相同。

從表1可以看出,對于H組中的高信噪比微震信號,3種方法均獲得了準確率很高的震相初至拾取結果,對比它們的算法平均耗時,可以看出HT-AIC算法耗時最少。

圖10　使用AIC法及HHT-AIC法拾取高信噪比微震P波初至結果對比Fig.10 Detection on micro seismic P-arrival time by using AIC and HHT-AIC method under high SNR

表1　高信噪比微震信號震相初至拾取結果Table 1 Result of detection on microseismic P-arrival time

分別應用AIC,HHT-AIC方法對圖9(b)所示的低信噪比微震信號進行震相初至拾取,與人工拾取的結果對比如圖11所示。圖11(a)為原始微震信號,起跳位置較為模糊,即使人工辨識信號起跳點也存在一定困難,通過信號局部放大,人工拾取的震相初至t=631 ms;圖11(b)為未經消噪處理直接應用AIC法求取的震相初至t=2 326 ms,拾取結果出現(xiàn)錯誤; 圖11(c)為經過EMD降噪后的微震信號,相較于圖11(a)而言信號波形較為清晰,且充分保留了微震信號的隨機非平穩(wěn)特征;圖11(d)為在消噪基礎上應用HHT-AIC法計算得到的震相初至t=617 ms,計算結果與人工拾取結果相比誤差在20 ms之內,說明HHT-AIC法能有效處理低信噪比信號震相初至的自動拾取,算法具有較強的抗噪性能。

從算法耗時來看,AIC方法耗時0.635 s,但自動拾取的震相初至出現(xiàn)錯誤; HHT - AIC法耗時1.912 s,雖然算法耗時較長,但獲得了正確的拾取結果。若在圖11(a)所示的微震波形上直接應用HTAIC法,算法耗時0.204 s,但拾取的震相初至為927 ms,與人工拾取結果631 ms相比產生了較大誤差,這是因為該含噪信號波形起跳不明顯,HT-AIC算法通過包絡閾值搜尋震相初至的大致位置出現(xiàn)了偏差,導致為AIC函數(shù)的選取計算時窗不合理,從而使AIC函數(shù)計算得到的震相初至出現(xiàn)較大誤差。因此,HT-AIC算法在處理低信噪比微震信號震相初至拾取時也是不穩(wěn)定。

對L組中50個低信噪比微震信號進行處理,以人工拾取的震相初至作為參照,震相初至拾取結果見表2。表2中HHT-AIC法震相初至的拾取準確率最高,而HT-AIC法耗時最少。綜合表1和2的拾取結果可以看出:對于高信噪比微震信號,HT-AIC法利用Hilbert變換對信號震相初至位置進行了粗略判定,為AIC函數(shù)選擇了位置及大小合適的計算時窗,有效降低了AIC函數(shù)的計算強度,因而能在保持較高識別準確率的前提下使算法耗時最少;但HT-AIC方法在處理低信噪比微震信號時,如果不對含噪信號進行降噪處理,仍會導致計算時窗位置選擇不合理,最終出現(xiàn)震相拾取結果錯誤;HHT-AIC法則首先應用EMD法對低信噪比微震信號進行降噪,提高了信號的信噪比,然后再應用HT-AIC算法實現(xiàn)震相初至的自動拾取,識別準確率比AIC法提高約10%,比HT-AIC法提高約6%。

圖11　使用AIC法及HHT-AIC法拾取低信噪比微震P波初至結果對比Fig.11 Detection on micro seismic P-arrival time by using AIC and HHT-AIC method under low SNR

表2　低信噪比微震信號震相初至拾取結果Table 2 Result of detection on microseismic P-arrival time

3.3 關于震相自動拾取方法的討論

AIC算法在時窗恰好包含有震相初至前后各一段波形的前提下,能夠得到較好的拾取結果,故AIC算法適用于已知震相初至大致位置的情況,工程上通常與STA/ LTA法組合使用,即首先使用STA/ LTA法獲取震相初至的大致位置,再以該位置為基礎向前及向后各取一段數(shù)據(jù)作為AIC函數(shù)的時窗計算出震相初至位置,但在處理低信噪比微震信號時,應用AIC 及STA/ LTA方法時震相初至辨識的準確率均顯著下降,算法抗噪性能較差。

HHT-AIC算法通過信號的EMD分解與IMF分量重構實現(xiàn)降噪,充分保留了微震信號固有的隨機非平穩(wěn)特征,并在消噪濾波的上應用Hilbert變換求得歸一化包絡信號波形,通過設定包絡閾值搜索震相初至的大致位置后,再應用AIC法求解震相初至,因而HHT-AIC算法具有較強的抗噪性能,在相同信噪比情況下,震相初至辨識的準確率顯著提高。

4 結 論

(1)對于高信噪比微震信號,HT-AIC方法能依據(jù)Hilbert變換及包絡閾值為AIC函數(shù)確定計算時窗位置及大小,包絡閾值參數(shù)選取簡單,自動拾取的震相初至與人工拾取結果一致,誤差在10 ms內拾取占比為100%,且算法實時性強,克服了傳統(tǒng)AIC方法由于時窗選擇不合理造成初至錯誤拾取的不足。

(2)對于低信噪比微震信號,HHT-AIC方法基于EMD分解及IMF重構兩個過程對含噪信號進行降噪,充分保留了微震信號固有的瞬態(tài)非平穩(wěn)特征,在此基礎上再應用HT-AIC方法實現(xiàn)震相初至的自動拾取,通過實驗與人工拾取結果相比,時差在10 ms內的微震信號占比為92%,自動拾取準確率高于AIC方法及HT-AIC方法,表現(xiàn)出較強的抗噪性能;不足之處是算法復雜度較高,但以犧牲時間為代價換取精度也是慣常的選擇。

參考文獻:

[1]齊慶新,竇林名.沖擊地壓理論與技術[M].徐州:中國礦業(yè)大學出版社,2008:32-38.

[2]鞏思園,竇林名,馬小平,等.煤礦礦震定位中異向波速模型的構建與求解[J].地球物理學報,2012,55(5):1757-1763.

Gong Siyuan,Dou Linming,Ma Xiaoping,et al.Study on the construction and solution technique of anisotropic velocity model in the location of coal mine tremor[J].Chinese Journal of Geophysics, 2012,55(5):1757-1763.

[3]朱權潔,姜福興,王存文,等.微震波自動拾取與多通道聯(lián)合定位優(yōu)化[J].煤炭學報,2013,38(8):397-403.

Zhu Quanjie, Jiang Fuxing, Wang Cunwen, et al.Automated micro-seismic event arrival picking and multi-channel recognition and location[J].Journal of China Coal Society,2013,38(8):397-403.[4]Allen R V.Automatic earthquake recognition and timing from single traces[J].Bull.Seismol.Soc.Amer.,1978,68 (5):1521 -1532.

[5]Allen R V.Automatic phase pickers:Their present use and future prospects[J].Bull.Seismol.Soc.Amer.,1982,72(6):225-242.[6]Baer M,Kardolfer U.An automatic phase picker for local and tele-seismic events[J].Bull.Seismol.Soc.Amer.,1987,77 (4): 1437-1445.

[7]Akaike H.Information theory and an extension of the maximum likelihood principle[A].2ndInternational Symposium on Information Theory (Tsahkadsor)[C].1971:267-281.

[8]Takanami T,Kitagawa G.Multivariate time-series model to estimate the arrival times of S-waves[J].Computers and Geosciences, 1993,19(2):295-301.

[9]Leonard M,Kennett M B L N.Multi-component autoregressive techniques for the analysis of seismograms[J].Phys.Earth Planet.Interiors,1999,113(1-4):247-264.

[10]Sleeman R,Orild V E.Robust automatic P-phase picking:an online implementation in the analysis of broadband seismogram recordings [J].Phys.Earth Planet.Interiors,1999,113(1-4):265-275.

[11]Maeda N.A method for reading and checking phase times in autoprocessing system of seismic wave data[J].J.Seismol.Soc.Jpn., 1985,38(3):365-379.

[12]Saragiotis, Christos D.PAI-S/ K: A robust automatic seismic P phase arrival identification scheme[J].IEEE Transactions on Geosciences and Remote Sensing,2002,40(6):1395.

[13]常 旭,劉伊克.地震記錄的廣義分維及其應用[J].地球物理學報,2002,11(6):839-846.

Chang Xu,Liu Yike.The generalized fractal dimension of seismic records and its application[J].Chinese Journal of Geophysics, 2002,11(6):839-846.

[14]馬 強,金 星,李山有,等.用于地震預警的P波震相到時自動拾取[J].地球物理學報,2013,56(7):2313-2321.

Ma Qiang,Jin Xing,Li Shanyou,et al.Automatic P-arrival detection for earthquake early warning[J].Chinese Journal of Geophysics,2013,56(7):2313-2321.

[15]王 繼,陳九輝,劉啟元,等.流動地震臺陣觀測初至震相的自動檢測[J].地震學報,2006,28(1):42-51.

Wang Ji,Chen Jiuhui,Liu Qiyuan,et al.Automatic onset phase picking for portable seismic array observation[J].Acta Seismologica Sinica,2006,28(1):42-51.

[16]劉勁松,王 赟,姚振興.微地震信號到時自動拾取方法[J].地球物理學報,2013,56(5):1660-1666.

Liu Jinsong,Wang Yun,Yao Zhenxing.On micro-seismic first arrival identification:A case study[J].Chinese Journal of Geophysics, 2013,56(5):1660-1666.

[17]Huang N E,Zheng S,Steven R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear non-stationary time series analysis[A].Proceedings:Mathematical,Physical and Engineering Sciences [C].London: The Royal Society Press,1998: 903-995.

[18]Rashid A.IIR discrete-time Hilbert transformers[J].IEEE Transactions on Acoustics,Speech,and Signal Processing,1987,35(8): 1116-1119.

Jia Ruisheng,Tan Yunliang,Sun Hongmei,et al.Method of automatic detection on micro-seismic P-arrival time under low signal to noise ratio[J].Journal of China Coal Society,2015,40(8):1845-1852.doi:10.13225/ j.cnki.jccs.2014.1122

Method of automatic detection on micro-seismic P-arrival time under low signal-to-noise ratio

JIA Rui-sheng1,2,TAN Yun-liang1,3,SUN Hong-mei2,HONG Yong-fa2

(1.State Key Laboratory of Mining Disaster Prevention and Control Co-founded by Shandong Province and Ministry of Science and Technology,Shandong University of Science and Technology,Qingdao 266590,China;2.College of Information Science and Engineering,Shandong University of Science and Technology,Qingdao 266590,China;3.College of Mining and Safety Engineering,Shandong University of Science and Technology,Qingdao 266590,China)

Abstract:According to the low signal-to-noise ratio,non-stationary,and inherent randomness of micro-seismic signals, a new method based on Hilbert-Huang Transform (HHT) and AIC criterion was proposed.Firstly,to implement denoising and filtering,the signals with noise were decomposed with empirical mode decomposition (EMD) and then reconstructed with intrinsic mode function (IMF);Secondly,based on the Hilbert transform,the normalized envelope signal was calculated,the coarse position of P-arrival was searched by setting the envelope threshold,and on the basis of the position,the time window of AIC function was determined;Finally,the P-arrival time was calculated by using AIC criterion in the selected time window.This method was applied to 80 random micro-seismic signals with both high SNR and low SNR.Considering the manual picking results as a benchmark and the time difference within 10 ms as an accurate result,the result shows that for the signals with high SNR,the picking accuracy is 100%.For the signals with low SNR,the accuracy is 92%.

Key words:micro-seismic P-arrival;seismic phase first arrival identification;empirical mode decomposition;Hilbert-Huang transform;Akaike information criterion

作者簡介:賈瑞生(1972—),男,安徽碭山人,教授,博士。E-mail:jrs716@163.com

基金項目:山東省自然科學基金資助項目(ZR2013EEM019);國家“十二五”科技支撐計劃資助項目(2012BAK04B06)

收稿日期:2014-08-30

中圖分類號:TD32

文獻標志碼:A

文章編號:0253-9993(2015)08-1845-08