低信噪比微震P波震相初至自動(dòng)拾取方法

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

(1.山東科技大學(xué)礦山災(zāi)害預(yù)防控制省部共建國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地,山東青島 266590;2.山東科技大學(xué)信息科學(xué)與工程學(xué)院,山東青島266590;3.山東科技大學(xué)礦業(yè)與安全工程學(xué)院,山東青島 266590)

?

低信噪比微震P波震相初至自動(dòng)拾取方法

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

(1.山東科技大學(xué)礦山災(zāi)害預(yù)防控制省部共建國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地,山東青島 266590;2.山東科技大學(xué)信息科學(xué)與工程學(xué)院,山東青島
266590;3.山東科技大學(xué)礦業(yè)與安全工程學(xué)院,山東青島 266590)

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

關(guān)鍵詞:微震P波;震相初至拾取;經(jīng)驗(yàn)?zāi)B(tài)分解;希爾伯特-黃變換;AIC準(zhǔn)則

責(zé)任編輯:常 琛

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

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

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

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

1 微震P波震相自動(dòng)拾取算法分析

1.1 長(zhǎng)短時(shí)均值比(STA/ LTA)法

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

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

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

其中,t0為當(dāng)前時(shí)刻;t1,t2為t0之前某時(shí)刻,且滿足t2

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

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

1.2 AIC方法

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

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

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

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

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

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

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

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

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

2 對(duì)震相自動(dòng)拾取算法的改進(jìn)

2.1 EMD變換及信號(hào)降噪處理

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

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

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

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

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

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

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

(5)剔除高頻IMF分量,對(duì)剩余的IMF分量進(jìn)行重構(gòu),得到降噪后的微震信號(hào):

圖5　一個(gè)含噪微震信號(hào)的EMD分解Fig.5 Empirical mode decomposition of a micro seismic signal with noise

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

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

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

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

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

式中,a(t)為實(shí)信號(hào)x(t)的包絡(luò)。

將降噪后的微震信號(hào)進(jìn)行Hilbert變換并進(jìn)行包絡(luò)分析,即按時(shí)序?qū)⑸舷掳j(luò)的絕對(duì)值相加并將其歸一化處理,則圖6(c)的Hilbert包絡(luò)信號(hào)波形如圖7所示,對(duì)同一類微震信號(hào)使用同一閾值(本文設(shè)為0.3)沿時(shí)間軸順序搜索,當(dāng)包絡(luò)值大于設(shè)定閾值時(shí),即可得到P波震相初至的大致位置t0,本例中t0= 1 243 ms。在經(jīng)過HHT降噪后的微震信號(hào)波形上以t0為基準(zhǔn)向前及向后分別取500個(gè)采樣點(diǎn)作為計(jì)算時(shí)窗,在該時(shí)窗內(nèi)應(yīng)用AIC函數(shù)求解P波震相初至,求解結(jié)果如圖8所示。

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

圖7　Hilbert包絡(luò)信號(hào)波形Fig.7 Waveform of Hilbert envelope signal

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

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

3 算法檢驗(yàn)及討論

3.1 實(shí)驗(yàn)數(shù)據(jù)來源

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

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

分別應(yīng)用AIC方法及本文方法對(duì)上述2組不同類型的微震信號(hào)進(jìn)行處理,并將其與人工拾取的震相初至結(jié)果對(duì)比,從拾取結(jié)果的正確率、算法耗時(shí)等方面進(jìn)行評(píng)價(jià)。

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

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

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

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

從表1可以看出,對(duì)于H組中的高信噪比微震信號(hào),3種方法均獲得了準(zhǔn)確率很高的震相初至拾取結(jié)果,對(duì)比它們的算法平均耗時(shí),可以看出HT-AIC算法耗時(shí)最少。

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

表1　高信噪比微震信號(hào)震相初至拾取結(jié)果Table 1 Result of detection on microseismic P-arrival time

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

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

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

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

表2　低信噪比微震信號(hào)震相初至拾取結(jié)果Table 2 Result of detection on microseismic P-arrival time

3.3 關(guān)于震相自動(dòng)拾取方法的討論

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

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

4 結(jié) 論

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

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

參考文獻(xiàn):

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

[2]鞏思園,竇林名,馬小平,等.煤礦礦震定位中異向波速模型的構(gòu)建與求解[J].地球物理學(xué)報(bào),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]朱權(quán)潔,姜福興,王存文,等.微震波自動(dòng)拾取與多通道聯(lián)合定位優(yōu)化[J].煤炭學(xué)報(bào),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]常 旭,劉伊克.地震記錄的廣義分維及其應(yīng)用[J].地球物理學(xué)報(bào),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]馬 強(qiáng),金 星,李山有,等.用于地震預(yù)警的P波震相到時(shí)自動(dòng)拾取[J].地球物理學(xué)報(bào),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]王 繼,陳九輝,劉啟元,等.流動(dòng)地震臺(tái)陣觀測(cè)初至震相的自動(dòng)檢測(cè)[J].地震學(xué)報(bào),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]劉勁松,王 赟,姚振興.微地震信號(hào)到時(shí)自動(dòng)拾取方法[J].地球物理學(xué)報(bào),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

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

基金項(xiàng)目:山東省自然科學(xué)基金資助項(xiàng)目(ZR2013EEM019);國(guó)家“十二五”科技支撐計(jì)劃資助項(xiàng)目(2012BAK04B06)

收稿日期:2014-08-30

中圖分類號(hào):TD32

文獻(xiàn)標(biāo)志碼:A

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