，，，

(1.山東科技大學(xué) 計算機(jī)科學(xué)與工程學(xué)院，山東 青島 266590；2.山東科技大學(xué) 山東省智慧礦山信息技術(shù)省級重點(diǎn)實驗室，山東 青島 266590)

微震的初至到時拾取對實現(xiàn)海量微地震數(shù)據(jù)的自動處理有重要意義[1]。由于微震監(jiān)測環(huán)境復(fù)雜，通常采集到的微震信號具有信噪比較低、非平穩(wěn)、非線性等特征。受煤礦井下機(jī)械振動的影響，微震信號常常夾雜大量的振動噪聲，這些噪聲對微震震相的初至到時拾取有很大干擾。因此，在機(jī)械振動噪聲干擾下研究微震震相初至的快速準(zhǔn)確拾取具有重要意義。

目前震相初至到時的自動拾取方法主要包括：長短時均值比方法(STA/LTA)、AIC法、PAI-S/K法、Hausdoff分形法等。其中，Stevenson等[2-4]提出了基于時間域的長短時均值比法(STA/LTA)； Akaike等[5-9]依據(jù)地震波到達(dá)前后數(shù)據(jù)統(tǒng)計的差別，提出了AIC準(zhǔn)則；Saragiotis等[10]基于地震波形的偏斜度和峰度提出了PAI-S/K方法；Chang等[11]提出了基于Hausdorff分形維識別震相的方法；Coppens[12]提出了在不同時窗內(nèi)進(jìn)行能量比較的方法；王繼等[4]應(yīng)用單臺Akaike信息準(zhǔn)則(AIC)和多臺最小二乘互相關(guān)方法，發(fā)展了震相自動精確檢測技術(shù)，實現(xiàn)了震相初至的自動拾?。?馬強(qiáng)等[13]綜合STA/LTA方法及AIC 準(zhǔn)則，提出了多步驟P波自動拾取方法；劉勁松等[14]通過分析STA/LTA，AIC，PAI-S/K等幾種方法的原理及特點(diǎn)，提出了移動時窗峰度的快速算法和改進(jìn)的峰度拾取初至算法。STA/LTA方法在拾取震相到時過程中受信噪比的影響較大，不易拾取低信噪比信號的準(zhǔn)確到時；而AIC算法拾取到時沒有選取時窗大小的規(guī)則，不易確定合適的窗長；分形維識別法、PAI-S/K方法、能量比法等都沒有對低信噪比信號進(jìn)行過多的研究。

由于煤礦井下環(huán)境復(fù)雜，采集到的信號還包含有大量機(jī)械振動噪聲，而上述拾取方法雖具有一定的抗噪能力，但對于低信噪比微震信號不能很好的拾取有效到時。目前工程上對于低信噪比微震信號的震相拾取主要是在專業(yè)軟件的支持下進(jìn)行人工拾取，拾取工作效率低，且震相初至的判識精度取決于操作人員的經(jīng)驗。為了解決低信噪比對到時拾取精度影響的問題，提出了一種基于低信噪比微震信號震相初至自動拾取方法。該方法首先應(yīng)用小波閾值降噪法對采集到的微震信號進(jìn)行降噪預(yù)處理，消除部分外部噪聲對到時拾取精度的干擾；然后用STA/LTA獲取信號的粗略到時，以該到時為基礎(chǔ)選取合適的時間窗口，在該事件窗口內(nèi)使用AIC方法計算震相的精確到時，從而實現(xiàn)低信噪比微震信號的震相初至的自動拾取。

1 震相拾取原理

本文的震相拾取方法是基于小波閾值降噪法的多方法聯(lián)合拾取，涉及到的震相識別方法有STA/LTA法及AIC方法。下面介紹兩種識別震相的基本理論基礎(chǔ)。

1.1 長短時平均(STA/LTA)方法

STA/LTA 方法是由 Stevenson 提出[2-4]，并應(yīng)用于地震初至到時的判別，后來Allen 等對該方法進(jìn)行了改進(jìn)并用于檢測地震事件和震相識別[15-17]。

STA/LTA是通過STA(short-term average)和LTA(long-term average)的比值，反映信號的幅度、頻率等特征的變化。其工作原理圖如圖1所示，當(dāng)?shù)卣鹦盘柕竭_(dá)時，由于短時窗STA內(nèi)幅值平均值變化比長時窗LTA內(nèi)變化劇烈，此時STA與LTA的比值會有一個突變，當(dāng)其比值大于某一設(shè)定閾值λ時，則判定為有效信號[18]，且此劇變點(diǎn)被判定為震相初至點(diǎn)。

STA/LTA拾取到時，是通過計算長短時窗特征函數(shù)的均值得到的STA、LTA的值，并通過STA和LTA的比值得到閾值。計算公式為:

(1)

(2)

(3)

其中，公式(1)～(3)中的n指微震信號的采樣時刻，sw為短時窗的長度，lw是長時窗的長度，α為設(shè)定的觸發(fā)閾值,CF(i)是在時刻i關(guān)于微震信號的特征函數(shù)值，表示微震數(shù)據(jù)的振幅、能量或其變化。

圖1 STA/LTA原理Fig.1 STA/LTA Principle

1.2 池赤準(zhǔn)則(AIC)法

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

圖2 AR-AIC模型Fig.2 AR-AIC Model

(4)

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

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

AIC(n)=n·log{var(x[1,n])}+(L-n-1)·log{var(x[n+1,L])}

。

(5)

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

圖3 使用AIC法拾取震相到時Fig.3 Onset time picking of seismic phase with AIC

通過以上2種AIC方法識別信號的對比，可見后者不需要計算AR模型的階數(shù)，即可直接求取AIC值。后者能滿足震相初至拾取的實時性較高的要求，但是該算法需要在震相初至的附近尋找一個合適的時窗來計算AIC值，這是因為不同的時窗可能使AIC函數(shù)局部最小值出現(xiàn)的位置不同。為同一個微震波形數(shù)據(jù)在不同時窗內(nèi)的波形，由于時窗設(shè)置不合理導(dǎo)致震相初至拾取錯誤。因此，如何選擇合理的時窗是AIC法拾取準(zhǔn)確震相初至的關(guān)鍵問題之一。

由于采集的微震信號包含有大量機(jī)械振動噪聲，而機(jī)械振動噪聲對微震震相到時的拾取是一項極大的挑戰(zhàn)，因此在進(jìn)行到時拾取之前要進(jìn)行機(jī)械振動噪聲的濾除。

1.3 小波閾值降噪法

為壓制機(jī)械噪聲對微震震相到時拾取的影響，使用小波閾值降噪法對低信噪比微震信號進(jìn)行降噪預(yù)處理。

小波閾值降噪[21]是在小波變換的基礎(chǔ)上發(fā)展起來的一種降噪方法。小波變換對于信號具有良好的去相關(guān)性，在實際環(huán)境下，有意義的信號一般是比較平穩(wěn)的或者低頻信號，但是含噪信號則為高頻?；诤性肼暤男盘栐谛〔ㄓ虻姆植继匦裕梢酝ㄟ^將時域信號轉(zhuǎn)換到小波域內(nèi)，通過設(shè)定適當(dāng)?shù)拈撝禐V除噪聲，從而達(dá)到對含有噪聲的信號進(jìn)行降噪處理的目的。

小波閾值降噪法可以采用硬閾值或軟閾值的處理方法對小波系數(shù)做閾值處理[22]。硬閾值是將指定閾值和信號的絕對值對比，把絕對值小于或等于閾值的置為零，對于絕對值大于閾值的信號則不變。軟閾值是將指定閾值和信號的絕對值對比，將小于或等于閾值的置為零，大于閾值的信號，設(shè)置為本身與閾值相減，這樣數(shù)據(jù)點(diǎn)就會向零靠攏。硬閾值、軟閾值的公式定義分別為公式(6)和公式(7)所示：

(6)

(7)

(8)

sgn(ωi,j)為符號函數(shù)，即：

(9)

2 微震震相初至拾取過程

2.1 算法流程

該微震震相初至拾取過程由三部分組成，分別為降噪、大致的到時拾取、精確拾取。

進(jìn)行第一部分降噪處理時，采用了小波閾值降噪法，對要進(jìn)行到時拾取的信號進(jìn)行預(yù)處理，消除大部分噪聲。第二部分進(jìn)行大致到時拾取，用的是STA/LTA方法。由于STA/LTA方法拾取到時的局限性，單純的用此方法得不到準(zhǔn)確的結(jié)果，將得到的到時作為下一步處理的一個中間值。第三部分為到時的精確拾取，在前兩部分的基礎(chǔ)上，用AIC方法進(jìn)行到時拾取，AIC的窗長根據(jù)第二部分得到的中間值確定。

微震信號到時拾取流程如圖4。

圖4 微震信號到時拾取流程Fig.4 Onset time picking process of microseismic signal

2.2 微震震相初至拾取算法描述

算法具體步驟如下：

Step 1: 獲取初始微震信號x(t)；

Step 2: 用小波閾值對信號x(t)進(jìn)行降噪預(yù)處理，得到降噪后的信號x′(t)；

Step 4: 利用設(shè)置好的STA/LTA算法對x′(t)進(jìn)行到時拾取，得到一個大致到時t1；

(10)

當(dāng)公式(10)的結(jié)果首次大于δ時，獲取t1的值；

Step 5： 在STA/LTA算法拾取到時的基礎(chǔ)上，確定AIC算法的時間窗長，以t1為對稱點(diǎn)，向前向后分別取500個采樣點(diǎn)；

L=[t1-500,t1+499]，

(11)

Step 6: 利用AIC算法對Step5中時窗內(nèi)的信號進(jìn)行微震信號的拾取，得到的到時為t

AIC(x′(t))=t1*log{var(x′[1,t1])}+(L0-t1-1)*log{var(x′[t1+1,Ln])}。

(12)

3 實驗

實驗數(shù)據(jù)來源于微震監(jiān)測系統(tǒng)，從監(jiān)測數(shù)據(jù)中隨機(jī)選取一組含噪微震信號進(jìn)行到時拾取實驗。

信號采集頻率設(shè)置為1 000 Hz，選取的采樣點(diǎn)數(shù)為4 000，該組信號的人工拾取到時點(diǎn)設(shè)置在1 229 ms處。原信號的波形圖與時頻譜如圖5。

圖5 原信號波形及時頻譜Fig.5 Waveform and frequency spectrum of the original signals

圖6為用STA/LTA以及AIC方法拾取的到時結(jié)果圖。用STA/LTA進(jìn)行到時拾取時，選取的短時窗為150 ms，長時窗為900 ms，設(shè)置的閾值為4，拾取的到時點(diǎn)為1 425 ms。如圖6所示，AIC方法理論拾取到時為1 231 ms，但是由于時窗選取太長，結(jié)果中的最低點(diǎn)可能不止一個。而且由于噪音的影響，STA/LTA、AIC方法拾取波形時，在拾取到時點(diǎn)附近分辨率不高，不能確定到時點(diǎn)。

圖6 不同方法對原信號進(jìn)行震相拾取Fig.6 Seismic phase picking of the original signals by different methods

針對上述出現(xiàn)的拾取到時不準(zhǔn)確以及分辨率不高的情況，提出了基于小波閾值降噪的低信噪比微震震相初至拾取方法。首先用小波閾值法對這組含噪的信號進(jìn)行降噪預(yù)處理，采用無偏風(fēng)險估計原則(rigrsure)[23]，設(shè)置小波分解層數(shù)為4層[24]，選取的小波基函數(shù)為Symlets小波系的sym8[25]。得到降噪后的信號波形圖(如圖7)。

圖7 降噪波形及頻譜圖Fig.7 De-noising waveform and frequency spectrum

用STA/LTA法對降噪后的信號進(jìn)行到時拾取，短時窗選取為150 ms，長時窗為900 ms，閾值設(shè)置為4，得到一個大致的到時，如圖8(b)所示。圖8(c)為用STA/LTA方法確定的時窗，該時窗是(b)中拾取到的到時(1 325±500) ms所得，即825～1 825 ms之間。

圖8 STA/LTA方法拾取降噪波形并確定時窗Fig.8 STA/LTA method of de-noising waveform picking and determination of the time window

圖9 本文方法拾取的到時Fig.9 Seismic phase picking of signals by the model

圖9為本文模型拾取的到時，即將圖8(c)得到的時窗作為AIC方法的窗長進(jìn)行到時拾取，得到最低點(diǎn)即到時點(diǎn)為1 231 ms。

表1為該實驗的性能分析，由本文方法拾取到時較傳統(tǒng)的STA/LTA方法運(yùn)行速度平均降低了34.28%；相對于AIC方法運(yùn)行平均速度提升了6.65%。

對50組微震信號進(jìn)行實驗分析，選取效果較好的實驗結(jié)果進(jìn)行不同方法拾取到時準(zhǔn)確率比對。以人工拾取到時為準(zhǔn)，誤差±10 ms作為本文所有拾取到時的準(zhǔn)確拾取，統(tǒng)計對比結(jié)果如表2所示。通過對比分析可見，本文方法拾取準(zhǔn)確率比傳統(tǒng)的STA/LTA方法平均提高9.78個百分點(diǎn)，比傳統(tǒng)的AIC方法平均提高10.08個百分點(diǎn)。

表1 一組信號不同方法拾取到時性能對比Tab.1 Performance comparison of different onset time picking methods for a set of signals

表2 不同方法拾取到時準(zhǔn)確率對比Tab. 2 Accuracy comparison of different onset time picking methods

4 結(jié)論與未來工作

通過對比STA/LTA、AIC以及本文方法在低信噪比環(huán)境下拾取微震震相拾取的準(zhǔn)確率，得出如下結(jié)論：STA/LTA算法簡單，拾取速度快，但在噪聲干擾下震相到時拾取誤差較大；AIC方法具有較強(qiáng)的抗噪能力，且能準(zhǔn)確拾取到時，但在時窗過長情況下算法運(yùn)行時間過長，不利于震相到時的實時拾取。本文方法使用STA/LTA法確定震相的大致到時，并以此為基礎(chǔ)為AIC算法確定一個盡可能小的時窗，大大縮短了AIC算法迭代的次數(shù)，不僅具有較強(qiáng)的抗噪能力，而且震相到時的拾取準(zhǔn)確率也較高。

使用小波閾值去噪能夠降低機(jī)械振動噪聲對微震信號的干擾，提高了震相到時拾取的準(zhǔn)確率，但由于微震信號采集環(huán)境復(fù)雜，尚未對其他類型的噪聲進(jìn)行深入研究，但初步實驗表明本文方法能夠?qū)﹄S機(jī)非平穩(wěn)噪聲進(jìn)行有效壓制，下一步我們將進(jìn)一步改進(jìn)小波閾值降噪算法，使其適應(yīng)更多類型的噪聲壓制，進(jìn)而提高震相初至拾取的精度及實時性。

[1]張文東,馬天輝,唐春安,等.錦屏二級水電站引水隧洞巖爆特征及微震監(jiān)測規(guī)律研究[J].巖石力學(xué)與工程學(xué)報,2014,33(2):339-348.

ZHANG Wendong,MA Tianhui,TANG Chun’an,et al.Research on characteristics of rockburst and rules of microseismic monitoring at diversion tunnels in Jinping II hydropower station[J].Chinese Journal of Rock Mechanics and Engineering,2014,33(2):339-348.

[2]田優(yōu)平.近震P波震相自動識別方法研究[D].北京：中國地震局地球物理研究所,2015.

[3]ANDREW H RANKIN.作為賦存于花崗巖中鎢錫礦化勘探指南的流體包裹體異常-前景展望[J].巖石學(xué)報,2007,23(1):3-14.

RANKIN H A.Fluid inclusion anomalies as exploration guides for granite hosted Sn-W mineralization:Prospects for the future[J].Acta Petrologica Sinica,2007，23(1):3-14.

[4]王繼,陳九輝,劉啟元,等.流動地震臺陣觀測初至震相的自動檢測[J].地震學(xué)報,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 (English Edition),2006，28(1):44-53,115.

[5]AKAIKE H.Information theory and an extension of the maximum likelihood principle[C]//2nd International Symposium on Information Theory.Tsahkadsor,1971:267-281.

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

JIA Ruisheng,TAN Yunliang,SUN Hongmei,et al.Method of automatic detection on micro-seismetic P-arrival time under low signal-to-noise ratio[J].Journal of China Coal Society,2015,40(8):1845-1852.

[7]LEONARD M,KENNETT B L N.Multi-component autoregressive techniques for analysis of seismograms[J].Physics of the Earth Planetary Interiors,1999,113(1-4):247-264.

[8]SLEEMAN R,ORILD V E.Robust automatic P-phase picking:An online implementation in the analysis of broadband seismogram recordings[J].Physics of the Earth planetary Interiors,1999,113(1/2/3/4):265-275.

[9]MAEDA N.A method for reading and checking phase times in autoprocessing system of seismic wave data[J].Journal of Seismological Society of Japan,1985,38(3):365-397.

[10]SARAGIOTIS,C.HADJILEONTIADIS L,PANAS S M.PAI-S/K:A robust automatic seismic P phase arrival identification scheme[J].IEEE Transactions on Geosciences and Remote Sensing,2002,40(6):1395.

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

CHANG Xu,LIU Yike.The generalized fractal dimension of seismic records and its applications[J].Chinese Journal of Geophysics,2002,45(6):839-846.

[12]COPPENS F.First arrival picking on common offset trace collections for automatic estimation of static correction[J].Geophysical Prospecting,1985,33(8):1212-1231.

[13]馬強(qiáng),金星,李山有,等.用于地震預(yù)警的P波震相到時自動拾取 [J].地球物理學(xué)報,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.

[14]劉勁松,王赟,姚振興.微地震信號到時自動拾取方法[J].地球物理學(xué)報,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.

[15]劉晗,張建中.微震信號自動檢測的STA/LTA算法及其改進(jìn)分析[J].地球物理學(xué)進(jìn)展,2014,29(4):1708-1714.

LIU Han,ZHANG Jianzhong.STA/LTA algorithm analysis and improvement of microseismic signal automatic detection[J].Progress in Geophysics,2014,29(4):1708-1714.

[16]吳治濤,李仕雄.STA/LTA算法拾取微地震事件P波到時對比研究[J].地球物理學(xué)進(jìn)展,2010,25(5):1577-1582.

WU Zhitao,LI Shixiong.Comparison of STA/LTA P-pickers for micro seismic monitoring [J].Progress in Geophysics,2010,25(5):1577-1582.

[17]白添羊,吳順川,王進(jìn)進(jìn),等.巖石破裂聲發(fā)射壓縮波到時拾取方法及其優(yōu)化改進(jìn)研究[J].巖石力學(xué)與工程學(xué)報,2016,35(9):1754-1766.

BAI Tianyang,WU Shunchuan,WANG Jinjin,et al.P-onset picking methods of acoustic emission compression waves and optimized improvement[J].Chinese Journal of Rock Mechanics and Engineering,2016,35(9):1754-1766.

[18]姜福興,尹永明,朱權(quán)潔,等.單事件多通道微震波形的特征提取與聯(lián)合識別研究[J].煤炭學(xué)報,2014,39(2):229-237.

JIANG Fuxing,YIN Yongming,ZHU Quanjie,et al.Feature extraction and classification of mining microseismic waveforms via multi-channels analysis[J].Journal of China Coal Society,2014,39(2):229-237.

[19]SLEEMAN R,ORILD V E.Robust automatic P-phase picking:An online implementation in the analysis of broadband seismogram recordings[J].Physics of the Earth Planetary Interiors,1999,113(1/2/3/4):265-275.

[20]MAEDA N.A method for reading and checking phase times in autoprocessing system of seismic wave data[J].Journal of Seismological Society of Japan,1985,38(3):365-379.

[21]DONOHO D L.Denoising by soft thresholding[J].IEEE Transactions on Information Theory,1995,41(3):613-627.

[22]王宏強(qiáng),尚春陽,高瑞鵬,等.基于小波系數(shù)變換的小波閾值去噪算法改進(jìn)[J].振動與沖擊,2011(10):165-168.

WANG Hongqiang,SHANG Chunyang,GAO Ruipeng,et al.An improvement of wavelet shrinkage denoising via wavelet coefficient transformation[J].Journal of Vibration and Shock,2011(10):165-168.

[23]王維,張英堂,任國全.小波閾值降噪算法中最優(yōu)分解層數(shù)的自適應(yīng)確定及仿真[J].儀器儀表學(xué)報,2009,30(3):526-530.

WANG Wei,ZHANG Yingtang,REN Guoquan.Adaptiveselection and simulation of optimal decomposition level in threshold de-noising algorithm based on wavelet transform[J].Chinese Journal of Scientific Instrument,2009,30(3):526-530.

[24]臧玉萍,張德江,王維正.小波分層閾值降噪法及其在發(fā)動機(jī)振動信號分析中的應(yīng)用[J].振動與沖擊,2009,28(8):57-60.

ZANG Yuping,ZHANG Dejiang,WANG Weizheng.Wavelet hierarchical threshold de-noising method and its application in vibration signal analysis of engine[J].Journal of Vibration and Shock,2009,28(8):57-60.

[25]鐘建軍,宋健,由長喜,等.基于信噪比評價的閾值優(yōu)選小波去噪法[J].清華大學(xué)學(xué)報(自然科學(xué)版),2014,54(2):259-263.

ZHONG Jianjun,SONG Jian,YOU Changxi,et al.Wavelet de-noising method with threshold selection rules based on SNR evaluations[J].Journal of Tsinghua University (Science and Technology),2014,54(2):259-263.