趙震華，張杏莉，盧新明

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

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

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

AIC算法是目前應(yīng)用最廣泛、拾取精度較高的震相初至拾取算法之一，具有運(yùn)算速度快、拾取準(zhǔn)確率高、噪聲魯棒性強(qiáng)等特點(diǎn)。王李管等[10]利用長(zhǎng)短時(shí)窗法初步定位P波S波的到時(shí)區(qū)間，在到時(shí)區(qū)間內(nèi)應(yīng)用自回歸赤池信息準(zhǔn)則(auto regression-Akaike information criterion,AR-AIC)算法，計(jì)算出對(duì)應(yīng)P波S波的初至到時(shí)，并結(jié)合時(shí)差閾值判別和時(shí)頻圖譜對(duì)拾取結(jié)果進(jìn)行檢測(cè)；賈瑞生等[11]利用Hilbert變換計(jì)算出包絡(luò)信號(hào)，設(shè)置包絡(luò)閾值，搜索出震相初至的大致位置，并設(shè)置合適的時(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算法識(shí)別P波初至。本研究對(duì)AR-AIC、VAR-AIC算法及其影響因素進(jìn)行全面分析對(duì)比，并結(jié)合實(shí)測(cè)煤礦井下微震信號(hào)的處理，分析AR模型階數(shù)、特征函數(shù)、時(shí)窗長(zhǎng)度對(duì)算法效果的影響。

1 AIC算法簡(jiǎn)介

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

(1)

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

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

由式(1)和式(2)可知，影響兩種AIC算法拾取精度的因素有AR模型階數(shù)、特征函數(shù)和時(shí)窗長(zhǎng)度。本部分對(duì)兩種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為微震信號(hào)長(zhǎng)度；σ2為AR模型系數(shù)。取式(3)最小值時(shí)的M，即為AR模型階數(shù)。利用式(3)得到某微震信號(hào)的AIC準(zhǔn)則曲線(xiàn)如圖1所示。

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

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

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

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

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

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

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

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

如式(2)所示，VAR-AIC算法的主要影響因素是能夠影響微震波形的特征函數(shù)(characteristic function，CF)，特征函數(shù)的選取會(huì)直接影響初至到時(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ù)分別對(duì)余弦函數(shù)進(jìn)行分析，對(duì)振幅變化的響應(yīng)特征如圖3所示。

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

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

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

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

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

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

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

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

2.3 時(shí)窗長(zhǎng)度

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

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

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

圖7 振幅最大值時(shí)刻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為振幅最大值時(shí)刻，l為時(shí)窗長(zhǎng)度。

VAR-AIC算法對(duì)圖7所示的微震信號(hào)拾取結(jié)果如圖8所示。圖8(a)～8(f)分別為時(shí)窗長(zhǎng)度為200、400、600、800、1 000和1 200 個(gè)采樣點(diǎn)的拾取結(jié)果，矩形框表示時(shí)窗長(zhǎng)度。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 時(shí)窗長(zhǎng)度的選取

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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