一種微震震相到時自動拾取方法

，，，

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

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

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

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

1 震相拾取原理

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

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

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

STA/LTA是通過STA(short-term average)和LTA(long-term average)的比值，反映信號的幅度、頻率等特征的變化。其工作原理圖如圖1所示，當地震信號到達時，由于短時窗STA內幅值平均值變化比長時窗LTA內變化劇烈，此時STA與LTA的比值會有一個突變，當其比值大于某一設定閾值λ時，則判定為有效信號[18]，且此劇變點被判定為震相初至點。

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

(1)

(2)

(3)

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

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

1.2 池赤準則(AIC)法

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

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

(4)

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

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

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

。

(5)

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

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

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

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

1.3 小波閾值降噪法

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

小波閾值降噪[21]是在小波變換的基礎上發展起來的一種降噪方法。小波變換對于信號具有良好的去相關性，在實際環境下，有意義的信號一般是比較平穩的或者低頻信號，但是含噪信號則為高頻。基于含有噪聲的信號在小波域的分布特性，可以通過將時域信號轉換到小波域內，通過設定適當的閾值濾除噪聲，從而達到對含有噪聲的信號進行降噪處理的目的。

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

(6)

(7)

(8)

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

(9)

2 微震震相初至拾取過程

2.1 算法流程

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

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

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

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

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

算法具體步驟如下：

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

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

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

(10)

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

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

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

(11)

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

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

(12)

3 實驗

實驗數據來源于微震監測系統，從監測數據中隨機選取一組含噪微震信號進行到時拾取實驗。

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

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

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

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

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

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

用STA/LTA法對降噪后的信號進行到時拾取，短時窗選取為150 ms，長時窗為900 ms，閾值設置為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方法的窗長進行到時拾取，得到最低點即到時點為1 231 ms。

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

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

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

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

4 結論與未來工作

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

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

[1]張文東,馬天輝,唐春安,等.錦屏二級水電站引水隧洞巖爆特征及微震監測規律研究[J].巖石力學與工程學報,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]田優平.近震P波震相自動識別方法研究[D].北京：中國地震局地球物理研究所,2015.

[3]ANDREW H RANKIN.作為賦存于花崗巖中鎢錫礦化勘探指南的流體包裹體異常-前景展望[J].巖石學報,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].地震學報,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].煤炭學報,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]常旭,劉伊克.地震記錄的廣義分維及其應用[J].地球物理學報,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]馬強,金星,李山有,等.用于地震預警的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.

[14]劉勁松,王赟,姚振興.微地震信號到時自動拾取方法[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.

[15]劉晗,張建中.微震信號自動檢測的STA/LTA算法及其改進分析[J].地球物理學進展,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].地球物理學進展,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].巖石力學與工程學報,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]姜福興,尹永明,朱權潔,等.單事件多通道微震波形的特征提取與聯合識別研究[J].煤炭學報,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]王宏強,尚春陽,高瑞鵬,等.基于小波系數變換的小波閾值去噪算法改進[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]王維,張英堂,任國全.小波閾值降噪算法中最優分解層數的自適應確定及仿真[J].儀器儀表學報,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]臧玉萍,張德江,王維正.小波分層閾值降噪法及其在發動機振動信號分析中的應用[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]鐘建軍,宋健,由長喜,等.基于信噪比評價的閾值優選小波去噪法[J].清華大學學報(自然科學版),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.