微震信號初至拾取的AIC算法及其分析

趙震華，張杏莉，盧新明

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

在沖擊地壓微震監測技術中，微震信號P波初至到時對于震源準確定位、沖擊地壓預測預警具有重要的意義[1]。微震信號初至拾取的高準確性是實現震源定位可靠性的前提。

目前，國內外針對微震信號初至拾取的方法主要有赤地信息量準則(Akaike information criterion,AIC)算法[2-3]、長短時窗平均能量比法(short-term average/long-term average,STA/LTA)[4-5]、高階統計量法[6]、人工神經網絡法[7-8]等。STA/LTA方法是目前常用的微震事件識別方法，利用時域信號短時窗能量平均值和長時窗能量平均值之比動態反映信號振幅變化，通過識別微震事件的初至時刻達到自動檢測微震事件的目的。該方法實現簡單，運算速度快，但誤判率較高。高階統計量法通過從非高斯信號的高階統計量中獲取有效信息，根據高階統計量中的峰度或偏斜度曲線的變化自動拾取震相到時，是研究非高斯、非線性和非平穩信號的有力工具，但在信噪比低時識別精度較差[9]。人工神經網絡法在拾取震相初至時，一般將表征地震波的特征量如偏振特征、信噪比、頻譜特征、振幅等作為輸入參數，將震相類別作為輸出，并使用樣本集訓練該網絡使之獲得模式分類識別功能，最后將地震波某一特征量輸入訓練好的網絡，若輸出值大于設定好的閾值，就認為獲得了震相的類型。人工神經網絡法具有良好的非線性映射功能，可并行處理大量信息，但存在工作量大、難以選擇隱含層節點數和傳遞函數的缺點。

AIC算法是目前應用最廣泛、拾取精度較高的震相初至拾取算法之一，具有運算速度快、拾取準確率高、噪聲魯棒性強等特點。王李管等[10]利用長短時窗法初步定位P波S波的到時區間，在到時區間內應用自回歸赤池信息準則(auto regression-Akaike information criterion,AR-AIC)算法，計算出對應P波S波的初至到時，并結合時差閾值判別和時頻圖譜對拾取結果進行檢測；賈瑞生等[11]利用Hilbert變換計算出包絡信號，設置包絡閾值，搜索出震相初至的大致位置，并設置合適的時窗，利用VAR-AIC算法計算P波震相初至；尚雪義等[12]利用經驗模態分解(empirical mode decomposition,EMD)自適應地將地震數據和噪聲分離為不同的本征模式函數(intrisic modcel function,IMF)，可靠地保留了P波震相到達幅值和相位譜，利用VAR-AIC算法識別P波初至。本研究對AR-AIC、VAR-AIC算法及其影響因素進行全面分析對比，并結合實測煤礦井下微震信號的處理，分析AR模型階數、特征函數、時窗長度對算法效果的影響。

1 AIC算法簡介

AIC[13]算法是由日本統計學家赤池弘次創立和發展的。AR-AIC[2]算法是一種基于AIC信息準則的到時拾取方法，該方法基于微震信號的非平穩特征，將信號劃分成若干個固定長度的波形段，然后進行自回歸(auto regression,AR)處理，求解AR模型階數和系數，AR-AIC算法表示為[14]

(1)

Maeda[3]對AR-AIC算法進行改進，提出VAR-AIC算法，該算法直接由微震波形數據計算AIC函數，求解AIC函數的最小值，該最小值對應的位置即為震相初至。VAR-AIC算法表示為

VAR-AIC(k)=kln{VAR(CF(Xt[1,k]))}+(N-k-1)ln{VAR(CF(Xt[k+1,N]))}。

(2)

式中：Xt[1,k]為微震波形數據Xt中采樣點1～k的振幅值，k的范圍為時間窗口內所有的微震信號采樣點，N為微震信號長度，VAR()為方差函數，CF()為特征函數。VAR-AIC算法計算簡單，拾取精度高，因此得到廣泛應用。

2 AIC算法的影響因素分析及參數確定

由式(1)和式(2)可知，影響兩種AIC算法拾取精度的因素有AR模型階數、特征函數和時窗長度。本部分對兩種AIC算法的影響因素進行分析討論并確定各參數的取值。

2.1 階數M的選取

AR模型階數M最常用的確定方法是AIC準則[15]，AIC準則的表達式為

AIC(M)=Nlnσ2+2M。

(3)

式中：M為AR模型階數；N為微震信號長度；σ2為AR模型系數。取式(3)最小值時的M，即為AR模型階數。利用式(3)得到某微震信號的AIC準則曲線如圖1所示。

圖1 AIC準則曲線圖Fig. 1 AIC criterion curve

由圖1可見，AIC(M)值在M=100時已經基本趨于穩定，在M=265時到達最小值。M取25、50、100、200、265和300時，AR-AIC算法對某微震信號的拾取結果如圖2所示。

(a)為手工拾取結果；(b)～(g)分別為AR-AIC算法在M為25、50、100、200、265和300時的拾取結果。圖2 不同AR模型階數M影響計算結果Fig. 2 Calculation results of difference AR model order M

由圖1和圖2可知，當M=50和100時，M值偏小，AIC(M)值還未穩定，AR-AIC算法拾取結果與手工拾取結果差距較大；當M=200和265時，AIC(M)值已基本穩定，拾取結果與手工拾取結果差異較小；當M=300時，AIC(M)值已經穩定，拾取結果與手工拾取結果差異相對較小。結果表明AIC準則確定的AR模型階數能夠使AR-AIC算法取得相對良好的效果。

AR-AIC算法的計算時間會隨著AR模型階數M的增大而增加。以長度N=4 000的微震信號為例，使用AR-AIC算法在不同階數M下的計算時間見表1。由表1可知，當M值增大時，算法耗時會顯著增加。因此，在選擇AR模型階數時，建議選擇AIC(M)最小值時的階數M或者選擇AIC(M)最小值稍前的階數M，但還需不斷計算，以保證該M值時拾取的初至到時與AIC(M)最小值時的M值拾取的初至到時相差不大，既保證拾取結果，又能適當降低算法運行時間。

表1 不同AR模型階數M的計算時間Tab. 1 Calculation time of different AR model order M

2.2 特征函數的選取

特征函數對微震信號的初至到時拾取至關重要，選取對微震信號振幅和頻率變化反應靈敏的特征函數，會提高算法對微震信號初至到時的拾取精度。

如式(2)所示，VAR-AIC算法的主要影響因素是能夠影響微震波形的特征函數(characteristic function，CF)，特征函數的選取會直接影響初至到時拾取的結果[16]。目前常用的特征函數有：

CF1(i)=|Xt(i)|，

(4)

(5)

CF3(i)=Xt(i)-Xt(i-1)，

(6)

(7)

CF5(i)=SUM(|Xt(i)|)，

(8)

(9)

(10)

利用7種特征函數分別對余弦函數進行分析，對振幅變化的響應特征如圖3所示。

(a)和(i)為振幅變化的余弦信號；(b)和(j)、(c)和(k)、(d)和(l)、(e)和(m)、(f)和(n)、(g)和(o)、(h)和(p)分別為特征函數CF1～CF7的拾取結果。圖3 7種特征函數對振幅變化的響應Fig. 3 Response of 7 characteristic functions to amplitude variation

由圖3可見，特征函數CF1～CF7對圖3(a)的余弦信號振幅變化皆有響應，其中CF4在分段點處的VAR-AIC值變化最大，響應最靈敏；CF7相比其他特征函數的響應程度稍弱一點。特征函數CF1～CF6對圖3(i)的余弦信號振幅變化皆有響應，其中CF4的響應最靈敏；CF7的響應最弱，并且拾取到錯誤的分段點。

對頻率變化的響應特征如圖4所示。特征函數CF1、CF2、CF5和CF6對圖4(a)和圖4(i)余弦信號頻率變化沒有響應，CF3和CF4對頻率變化有響應。相比CF3，CF4在分段點處的VAR-AIC值變化更大，響應更為靈敏，CF7雖對頻率變化有響應，但響應程度較弱。

(a)和(i)為頻率變化的余弦信號；(b)和(j)、(c)和(k)、(d)和(l)、(e)和(m)、(f)和(n)、(g)和(o)、(h)和(p)分別為特征函數CF1～CF7的拾取結果。圖4 7種特征函數對頻率的響應Fig. 4 Response of seven characteristic functions to frequency

由此得出的7種特征函數對振幅和頻率變化的響應特性見表2。

表2 7種特征函數對振幅和頻率的變化響應特性Tab. 2 The response characteristics of seven characteristic functions to the variations of amplitude and frequency

圖5為VAR-AIC算法在7種特征函數作用下的某實測微震信號初至拾取結果。可見，VAR-AIC算法對 7種特征函數作用下的微震信號都能進行初至拾取，但拾取結果受到一定影響。圖5(d)和圖5(e)中，CF3和CF4作用下的微震信號拾取結果比CF1、CF2、CF5、CF6和CF7更接近手工拾取結果，拾取結果更優，故建議選用特征函數CF3或CF4來進行微震信號初至拾取。

(a)為手工拾取結果；(b)～(h)分別為特征函數CF1～CF7作用下VAR-AIC算法拾取結果(虛線)圖5 不同特征函數VAR-AIC算法拾取結果Fig. 5 VAR-AIC algorithm picks the first arrival results using different characteristic functions

2.3 時窗長度

在VAR-AIC算法實際計算過程中，會出現AIC函數最小值不在信號初至位置的情況，如圖6所示。

(a)為某微震信號；(b)為VAR-AIC算法拾取錯誤初至圖6 VAR-AIC算法拾取錯誤初至Fig. 6 VAR-AIC algorithm picks up error first arrival

為避免此類情況的發生，將微震信號的振幅最大值時刻p作為時窗的結束時刻(如圖7)，從p向前取長度為l個采樣點為初至拾取時窗，并在時窗內計算AIC函數最小值。為了分析時窗長度對拾取結果的影響，分別取l為200、400、600、1 000和1 200個采樣點，且特征函數選擇CF4，此時式(2)變為

圖7 振幅最大值時刻pFig. 7 Maximum amplitude moment p

VAR-AIC(k)=kln{VAR(CF4(Xt[p-l,k]))}+
(N-k-1)ln{VAR(CF4(Xt[k+1,p]))}。

(11)

式中，p為振幅最大值時刻，l為時窗長度。

VAR-AIC算法對圖7所示的微震信號拾取結果如圖8所示。圖8(a)～8(f)分別為時窗長度為200、400、600、800、1 000和1 200 個采樣點的拾取結果，矩形框表示時窗長度。

圖8 不同時窗長度的VAR-AIC算法拾取結果Fig. 8 VAR-AIC algorithm picks up signals with different time window lengths

由圖8可見，時窗長度對VAR-AIC算法的拾取結果有較大影響。圖8(a)和圖8(b)的拾取結果與手工拾取結果相差較大，而圖8(c)～8(f)的拾取結果與手工拾取結果誤差僅為8個采樣點。可見，當時窗長度過小時，時窗內手工拾取點后的微震數據占總時窗長度的25%以上，計算初至點前的方差偏大，從而導致拾取的初至點稍靠后。為保證拾取結果的準確性，建議在選擇時窗長度時，開始時刻為信號最大振幅時刻向前移動600～1 000個采樣點，結束時刻選擇信號最大振幅時刻，且時窗內需要包含信號初至時刻。

3 影響因素實例驗證

3.1 AR模型階數M的選取

選取10條實測微震信號記錄編號1～10，如圖9所示。

圖9 10條實測微震信號記錄Fig. 9 10 micro-seismic signal records

利用AIC準則得到10條微震信號記錄的AIC(M)的最小值，并利用AR-AIC算法得到各階數M下的初至時刻結果如表3所示，AR-AIC算法的計算用時如表4所示。

表3 各階數M下的初至時刻Tab. 3 The first arrival time of each order ms

由表3可知，在AIC(M)最小值時的階數M的拾取結果相比于其他階數的結果更接近手工拾取，效果更好。由表4可知，運行時間會隨著階數的增大而增加，因此階數M應選擇AIC(M)最小值時的階數M或者選擇AIC(M)最小值稍前的階數M，既保證拾取效果，也盡可能減少算法的計算時間。

表4 AR-AIC算法在各階數M下的計算用時Fig. 4 Running time of AR-AIC algorithm in different orders M s

3.2 特征函數的選取

選取圖9中的10條實測微震信號，利用VAR-AIC算法分別計算10條微震信號記錄在不同特征函數下的拾取結果，如表5所示。

表5 不同特征函數下的初至時刻Tab. 5 First arrival time under different characteristic functionsms ms

由表5可知，與手工拾取相比，特征函數CF4作用下的微震信號更接近于手工拾取結果，魯棒性更強，因此特征函數選擇CF4。此外，其他特征函數作用的微震信號在拾取初至時會出現拾取錯誤的情況，與圖6所示結果一樣，原因為初至時刻的AIC值不是AIC函數的全局最小值，導致拾取錯誤。

3.3 時窗長度的選取

同樣選取圖9所示的10條微震信號，利用式(11)計算10條微震信號記錄在不同時窗長度下的拾取結果，特征函數選擇CF4。時窗的開始位置為結束位置向前移動的時窗大小，結束位置為最大振幅值時刻，拾取結果如表6所示。

由表6可知，與特征函數CF4的拾取結果相比，時窗長度在200～400 個采樣點之間的拾取效果與CF4的拾取效果相差較大，600～1 200 個采樣點之間的拾取效果更接近于特征函數CF4的拾取結果，魯棒性更強，故時窗長度應選擇在最大振幅值之前的600～1 200個采樣點。

表6 不同時窗長度下的初至時刻Tab. 6 First arrival time under different window lengthms

通過對以上10條實測微震信號記錄的研究及分析，得出的結論與第2部分中給出的建議基本一致，驗證了所得結論的正確性。

4 AIC算法應用分析

利用兩種AIC算法及其改進算法拾取微震信號初至應用較多。由于AR-AIC算法本身存在計算量巨大、拾取精度不高、計算時間長等缺點，決定了該算法在微震信號初至拾取中應用相對較少。Leonard[17]和王海軍等[18]分別對AR-AIC算法進行改進，提出混合初至估計的方法，改善了低信噪比震相的拾取效果。而在改進VAR-AIC算法時，一種思路是對待拾取的信號波形進行改變，例如本研究利用特征函數作用于信號波形，使其對信號初至拾取產生影響。劉希強等[19]提出利用三階累積量進行精細震相拾取的TOC-AIC(third-order cumulant AIC)算法，即用兩個時間段數據的三階累積量代替VAR-AIC算法中的方差。TOC-AIC算法無論是在信噪比較低還是較高的情況下，得到的震相初至結果與實際結果都非常接近。

由于VAR-AIC算法存在拾取錯誤的問題，衍生出另外一種思路，即對微震信號長度進行限制，避免拾取錯誤。崔云潔等[20]應用小波閾值降噪方法對微震信號進行降噪，利用STA/LTA方法進行預拾取，獲取初至大致時刻，并以該時刻為基準，向前向后各截取500個采樣點構成時間窗口，在時窗內利用VAR-AIC算法拾取初至；朱權潔等[16]對原有的小波包去噪法模型進行改造，建立基于多去噪規則的多閾值去噪法，在VAR-AIC算法的基礎上，提出了具有時窗限制的I-AIC(improved AIC)算法，精確拾取水力壓裂微震信號P波初至，與其他4種方法相比，I-AIC算法對去噪結果更為敏感，拾取精度、計算速度更優。

通過對上述幾種AIC改進算法的應用分析，說明對兩種AIC算法的影響因素進行研究很有必要，可以為AIC算法在特征函數或時窗長度上進一步改進，以獲得精度更高的初至拾取結果提供參考。

5 結論

通過對AIC算法的實驗和應用分析，得出以下結論：

1) AR-AIC算法中，AR模型階數M對拾取初至結果有較大影響。M值越大，算法復雜度越高，建議選擇AIC(M)最小值時的階數M或者選擇AIC(M)最小值稍前的階數M，該M值下算法拾取的初至點與AIC(M)最小值時的階數M值拾取的初至點相差較小。這樣既能保證拾取效果，又能適當降低算法計算用時。

3) 為避免拾取錯誤，應確定合適的時窗長度，建議時窗長度為600～1 000個采樣點，時窗內需包含信號初至，時窗結束時刻通常選擇信號最大振幅值時刻，時窗開始時刻為信號最大振幅時刻向前移動600～1 000個采樣點。這樣既能保證拾取效果、減少拾取誤差，又縮短了計算時長。