趙震華，張杏莉，盧新明

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

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

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

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

1 AIC算法簡介

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

(1)

Maeda[3]對AR-AIC算法進(jìn)行改進(jìn)，提出VAR-AIC算法，該算法直接由微震波形數(shù)據(jù)計(jì)算AIC函數(shù)，求解AIC函數(shù)的最小值，該最小值對應(yīng)的位置即為震相初至。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]為微震波形數(shù)據(jù)Xt中采樣點(diǎn)1～k的振幅值，k的范圍為時間窗口內(nèi)所有的微震信號采樣點(diǎn)，N為微震信號長度，VAR()為方差函數(shù)，CF()為特征函數(shù)。VAR-AIC算法計(jì)算簡單，拾取精度高，因此得到廣泛應(yīng)用。

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

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

2.1 階數(shù)M的選取

AR模型階數(shù)M最常用的確定方法是AIC準(zhǔn)則[15]，AIC準(zhǔn)則的表達(dá)式為

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

(3)

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

圖1 AIC準(zhǔn)則曲線圖Fig. 1 AIC criterion curve

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

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

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

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

表1 不同AR模型階數(shù)M的計(jì)算時間Tab. 1 Calculation time of different AR model order M

2.2 特征函數(shù)的選取

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

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

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

(4)

(5)

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

(6)

(7)

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

(8)

(9)

(10)

利用7種特征函數(shù)分別對余弦函數(shù)進(jìn)行分析，對振幅變化的響應(yīng)特征如圖3所示。

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

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

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

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

由此得出的7種特征函數(shù)對振幅和頻率變化的響應(yīng)特性見表2。

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

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

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

2.3 時窗長度

在VAR-AIC算法實(shí)際計(jì)算過程中，會出現(xiàn)AIC函數(shù)最小值不在信號初至位置的情況，如圖6所示。

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

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

圖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所示的微震信號拾取結(jié)果如圖8所示。圖8(a)～8(f)分別為時窗長度為200、400、600、800、1 000和1 200 個采樣點(diǎn)的拾取結(jié)果，矩形框表示時窗長度。

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

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

3 影響因素實(shí)例驗(yàn)證

3.1 AR模型階數(shù)M的選取

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

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

利用AIC準(zhǔn)則得到10條微震信號記錄的AIC(M)的最小值，并利用AR-AIC算法得到各階數(shù)M下的初至?xí)r刻結(jié)果如表3所示，AR-AIC算法的計(jì)算用時如表4所示。

表3 各階數(shù)M下的初至?xí)r刻Tab. 3 The first arrival time of each order ms

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

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

3.2 特征函數(shù)的選取

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

表5 不同特征函數(shù)下的初至?xí)r刻Tab. 5 First arrival time under different characteristic functionsms ms

由表5可知，與手工拾取相比，特征函數(shù)CF4作用下的微震信號更接近于手工拾取結(jié)果，魯棒性更強(qiáng)，因此特征函數(shù)選擇CF4。此外，其他特征函數(shù)作用的微震信號在拾取初至?xí)r會出現(xiàn)拾取錯誤的情況，與圖6所示結(jié)果一樣，原因?yàn)槌踔習(xí)r刻的AIC值不是AIC函數(shù)的全局最小值，導(dǎo)致拾取錯誤。

3.3 時窗長度的選取

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

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

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

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

4 AIC算法應(yīng)用分析

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

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

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

5 結(jié)論

通過對AIC算法的實(shí)驗(yàn)和應(yīng)用分析，得出以下結(jié)論：

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

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