幾種初至自動拾取算法在地震波法隧道超前探測中的對比分析

王 偉， 高 星， 劉孝衛(wèi)

(中國科學(xué)院地理科學(xué)與資源研究所 資源與環(huán)境信息系統(tǒng)國家重點實驗室， 北京 100101)

0 引言

隧道超前地質(zhì)預(yù)報進行地震數(shù)據(jù)處理時，地震波初至拾取精度會直接影響隧道掘進前方異常地質(zhì)體預(yù)測結(jié)果的準確度[1]，這就需要選擇初至拾取精度高的算法來進行拾取，而目前的地震初至波自動拾取的方法有很多。

在20世紀70年代以前，初至的拾取主要是靠人工根據(jù)波形特征及經(jīng)驗來完成的，過程枯燥、耗時且拾取結(jié)果容易因為拾取者的主觀因素出現(xiàn)偏差。在20世紀80年代至90年代末，隨著計算機的普及和相關(guān)技術(shù)的發(fā)展，大量的處理軟件實現(xiàn)了在窗口顯示地震波形，通過人機交互的方式來進行初至的拾取，拾取效率得到了很大的提高，但這并沒有解決耗時和主觀因素干擾的問題。進入21世紀后，大量的初至自動拾取算法被提出并實現(xiàn)應(yīng)用，包括Chen 和 Stewart于2005年提出的長短時窗算法、Sabbiobe 和 Velis于2010年改進的能量比值計算公式以及神經(jīng)網(wǎng)絡(luò)、小波變換等算法，極大地促進了初至自動拾取算法的發(fā)展。

事實上，地震波初至拾取是一項繁重復(fù)雜但非常重要的工作，現(xiàn)階段還沒有一種算法可以完美地解決如噪聲干擾大、信噪比低等環(huán)境下地震數(shù)據(jù)的初至拾取精度低的問題，至今研究人員仍在尋找精細準確的自動拾取方法[2]。

本文所使用的數(shù)據(jù)來源于TSP-SK(tunnel seismic prediction)型隧道地質(zhì)超前預(yù)報儀，如圖1所示。

(a) 隧道地質(zhì)超前預(yù)報儀

(b) 工作示意圖

Fig. 1 TSP-SK tunnel advance forecaster and its working sketch

該儀器運行穩(wěn)定，數(shù)據(jù)記錄清晰，能夠滿足本次初至拾取對比分析對數(shù)據(jù)質(zhì)量的要求。目前，隧道超前探測常用的初至自動拾取的方法主要有4種，分別為最大振幅法(TST，tunnel scatter technology使用的方法為振幅比法)、能量比法(TSP203、TSP303使用該方法)、STA/LTA(short time average /long time average，短時間的均值/長時間的均值)法、AIC(akaike information criterino)法。初至拾取的方法雖多，但其檢測震相初至的理論依據(jù)卻相似，即根據(jù)信號與噪音以及各震相在某領(lǐng)域(時域、頻域或時頻域)或某些屬性(運動學(xué)或動力學(xué))方面的差異來區(qū)分信號和噪聲，進而得到震相初至[3-4]。由于各個算法依據(jù)的地震屬性有很大的差異，因此對各初至拾取算法的拾取結(jié)果進行對比是非常必要的。本文通過對常用的4種初至拾取的算法進行對比，分析各種算法的拾取精度，力求在現(xiàn)階段找到一種簡單、高效的初至拾取算法。

1 初至拾取方法

1.1 最大振幅法

最大振幅法(最大能量法)，即選取地震道上能量最大的點作為地震波初至。優(yōu)點是算法操作簡單，計算速度快；缺點是受噪聲影響大，且容易產(chǎn)生錯誤[5]。

1.2 能量比法

該方法首先選取某一特定的地震屬性道，如振幅絕對值、波形長度、振幅包絡(luò)等，然后根據(jù)相鄰前、后窗口的屬性值比值特征確定地震波初至的位置。選取合適的地震屬性道是該方法成功的關(guān)鍵[6-7]。

本文利用能量比法拾取地震波初至?xí)r，采用變換時窗統(tǒng)計法記錄地震初至前后時窗內(nèi)的能量差異特征，其中時窗為固定大小的時間量[8-9]。能量比法的基本公式為：

(1)

(2)

式(1)—(2)中：A′為時窗內(nèi)后、前能量比；x(t)為地震記錄振幅值；T1為時窗起點；T0為時窗中點；T2為時窗終點；A為1個地震道的相對能量；α為穩(wěn)定系數(shù)，本文取值為3[10-11]。

利用該原理，提取能量比值函數(shù)A的最大值出現(xiàn)的位置作為初至到時拾取點[12]，如圖2所示。

圖2 第2道地震記錄能量比法初至到時拾取圖

Fig. 2 No.2 seismic record′s first arrival time pick-up by energy ratio method

1.3 AIC法

該方法的基本原理是求取地震信號AIC函數(shù)局部最小值的位置[13-14]。對于長度為N的信號x，其AIC函數(shù)可表示為：

AIC(k)=klog{var(x[1,k])}+(N-k-1)log{var(x[k+1,N])}。

(3)

式中： AIC即峰度赤池信息量準則，由日本統(tǒng)計學(xué)家赤池弘次創(chuàng)立和發(fā)展； var表示序列的方差函數(shù)。

AIC函數(shù)在初至點峰值尖銳，準確率較高；但AIC函數(shù)受時窗的影響較大，在不同的時窗下，其局部最小值出現(xiàn)的位置往往不同[7,15]。

利用該原理，在初步選定的初至到時點附近建立1個新的時窗，計算出這個時窗內(nèi)的AIC函數(shù)值，并提取能量函數(shù)AIC的最小值出現(xiàn)的位置作為初至到時拾取點[3,16]，如圖3所示。

圖3 第2道地震記錄AIC法初至到時拾取圖

Fig. 3 No.2 seismic record′s first arrival time pick-up by AIC method

1.4 STA/LTA法

Allen于1978年提出地震波STA/LTA 法[17]，其主要原理是根據(jù)地震波形特征函數(shù)的長短時均比值等特征拾取初至[18-19]。當信號到達時，STA要比LTA變化更快，相應(yīng)的STA/LTA值會出現(xiàn)明顯的增加，當比值大于設(shè)定的閾值時，就可判定初至波到達[20-21]。該方法容易受到時窗長度和滑動步長的影響[22-23]。

(4)

(5)

(6)

式(4)—(6)中： STA、LTA為短、長時窗的振幅均值；B為信號的振幅;X為閾值，經(jīng)過不斷地對比拾取分析，本文取值為11[24-25]。

利用該原理，設(shè)計滑動時窗并依次計算STA/LTA的值，短時窗的值為固定值，而長時窗的值不斷依次累加，這樣可以減少突變點對初至提取的干擾。計算的同時不斷判斷STA/LTA的值是否超過閾值，判斷流程如圖4所示。若該值大于閾值X(X=1，取該閾值是通過對數(shù)據(jù)多次運算分析得到)，則將該點作為初至到時拾取點，如圖5所示[1,26]。

圖4 STA/LTA算法流程圖

圖5 第2道地震記錄STA/LTA法初至到時拾取圖

Fig. 5 No.2 seismic record′s first arrival time pick-up by STA/LTA method

2 對比分析

2.1 數(shù)據(jù)采集

數(shù)據(jù)來源于云南某一隧道超前探測實例，采用TSP-SK三分量加速度檢波器單點接收、多點激發(fā)的方式，檢波點遠離隧道掌子面，采用1—19號有效激發(fā)數(shù)據(jù)，單炮藥量70 g，炮間距1.5 m，最小炮檢距20 m，0.062 5 ms采樣，采樣點數(shù)為7 218[27]。現(xiàn)場布置如圖1(b)所示，有效激發(fā)19炮記錄的x方向分量地震數(shù)據(jù)如圖6所示。由圖6可以看出： 各個道地震記錄88 ms之前的記錄占主要成分的為地震波能量，88 ms后的記錄參雜了較多成分的高頻聲波，這部分為干擾波，必須在后續(xù)的處理中予以消除；第11道的地震記錄受到的干擾比較大，因此后續(xù)分析應(yīng)適當減少各個初至拾取方法對第11道數(shù)據(jù)初至拾取結(jié)果的參考權(quán)重；初至到時從第1道到第19道是逐漸減小的，大體上呈線性減小，因此在確定第1道的初至后，可以確定在首道初至到時的時窗長度內(nèi)基本可以拾取到所有地震記錄道的初至到時，這樣可以大大減少算法的數(shù)據(jù)計算量，提高拾取效率。

圖6 TSP-SKA地震記錄圖

2.2 數(shù)據(jù)預(yù)處理

2.2.1 地震記錄去噪處理

由圖6可以看出，由于爆炸聲波或者環(huán)境噪音等干擾因素使得波形記錄中參雜了許多噪聲。本文利用巴特沃茲帶通濾波器，濾去主頻帶范圍外的噪聲，即將80～1 000 Hz內(nèi)的數(shù)據(jù)保留，有效消除由于噪聲引起的高低頻干擾。濾波后的部分道數(shù)據(jù)記錄如圖7和圖8所示，可以看出，經(jīng)過濾波去除了高頻干擾信號。

(a) 地震記錄原始圖

(b) 地震記錄濾波圖

2.2.2 提取特征函數(shù)

提取特征函數(shù)是描述地震波某類特征差異的函數(shù)，可以用CFi表示。不同學(xué)者提出了不同形式的特征函數(shù)，目的是為了反映地震波到達時信號的變化規(guī)律[28-29]。本文使用的特征函數(shù)為：

(7)

(a) 地震記錄原始圖

(b) 地震記錄濾波圖

選用此特征函數(shù)的原因是，其既可以反映波形信息的振幅變化，又可以反映頻率變化。推導(dǎo)過程如下。

設(shè)

x(n)=Bcos[ω(n)Δt+Φ]，n=1,2,3,…。

(8)

式中：ω為圓頻率；Φ為相位角。

則

x(n)=Bcos[ω(n)Δt+Φ]；

(9)

x(n-1)=Bcos[ω(n-1)Δt+Φ]；

(10)

x(n+1)=Bcos[ω(n+1)Δt+Φ]。

(11)

根據(jù)三角恒等式：

cos(2α)=2cos2α-1=1-2sin2α。

(12)

將式(10)和式(11)相乘得：

x(n-1)·x(n+1)=B2cos2[ω(n)Δt+Φ]-B2sin2ω。

(13)

將式(9)帶入式(13)可得：

B2sin2ω=x(n)2-x(n-1)·x(n+1)。

(14)

部分道數(shù)據(jù)的特征函數(shù)如圖9和圖10所示，通過與圖7和圖8對比可知，特征函數(shù)可以更加突出、準確地反映出地震記錄波形的各個屬性的變化，對提高初至拾取精度效果非常明顯，論證如下。

圖9 第8道特征函數(shù)圖

圖10 第16道特征函數(shù)圖

為了方便，本文選STA/LTA法對特征函數(shù)的作用予以對比說明，采用10道地震記錄數(shù)據(jù)。首先使用濾波后的波形數(shù)據(jù)作為能量比法的基礎(chǔ)分析數(shù)據(jù)進行初至拾取，再利用特征函數(shù)作為基礎(chǔ)分析數(shù)據(jù)進行初至拾取，可以得到表1中的結(jié)果。經(jīng)計算可以得到，以原始數(shù)據(jù)為基礎(chǔ)數(shù)據(jù)的初至拾取平均差值為0.547 ms，而以特征函數(shù)為基礎(chǔ)的平均差值只有0.292 ms，因此使用特征函數(shù)對于提高初至拾取算法的精度有非常明顯的作用。

表1以原始記錄為基礎(chǔ)和以特征函數(shù)為基礎(chǔ)的初至拾取到時對比表

Table 1 First arrival time pick-up comparison based on original record and characteristic function ms

2.3 4種初至波拾取方法對比

利用MATLAB數(shù)據(jù)處理平臺，結(jié)合4種算法的原理編寫其自動拾取函數(shù)，并利用4種不同的地震波初至自動拾取算法對同一組隧道超前探測地震記錄事件進行初至拾取，可以得到19道地震數(shù)據(jù)的初至到時，見表2。

由表2可以得知，最大振幅法拾取對信噪比的要求很高，容易受到環(huán)境雜波的干擾，在下文的對比中將不再對最大振幅法進行過多的分析。為了保障手動拾取初至到時結(jié)果的科學(xué)性，通過不同人進行10多次的拾取，然后取其平均值，這樣可以較大程度地減少人為因素的干擾。

表2 地震波初至到時拾取表

不同算法拾取結(jié)果的折線統(tǒng)計圖見圖11。

圖11 不同算法拾取結(jié)果的折線統(tǒng)計圖

Fig. 11 Polyline statistical graph of different pick-up consequences

由圖11可以看出，AIC法拾取的初至到時大部分較手動拾取的到時提前，能量比法相對于手動拾取結(jié)果也大部分提前，而STA/LTA法的拾取結(jié)果與手動拾取結(jié)果的波形較為一致。 2—15道STA/LTA法與手動拾取的結(jié)果極為接近；第1道AIC法和能量比法比較接近手動拾取結(jié)果；16—19道幾種方法的初至拾取點均比手動拾取的點要延后，其主要原因之一是炮檢距縮小使得算法空間較小，從而影響拾取精度。從折線圖可以初步判斷，STA/LTA法對于超前探測地震波初至自動拾取的精度較高。下文通過增加地震事件數(shù)，繼續(xù)利用這幾種方法進行拾取精度對比分析，驗證初步結(jié)論的準確性。

選取不同隧道超前探測項目中記錄到的數(shù)據(jù)，以1次隧道超前探測產(chǎn)生的所有地震數(shù)據(jù)合稱為1次地震事件。事件1至事件4表示4次超前探測產(chǎn)生的地震數(shù)據(jù)，其原始數(shù)據(jù)如圖12所示。利用幾種算法繼續(xù)進行初至拾取，將拾取的各個道到時與手動拾取的結(jié)果進行差值，并計算所有地震事件拾取結(jié)果的均方差值，結(jié)果見表3。

(a) 事件1

(b) 事件2

(c) 事件3

(d) 事件4

Table 3 Deviation of four first arrival time pick-up and manual pick-up ms

由表3可知，4種初至拾取精度最高的是STA/LTA法，最低的是最大振幅法，精度由高到低依次為STA/LTA法、能量比法、AIC法、最大振幅法。

3 結(jié)論與討論

在地震數(shù)據(jù)質(zhì)量及信噪比一定的條件下，地震波初至?xí)r間拾取算法影響初至?xí)r間拾取精度，通過對比分析人工拾取法與4種自動拾取算法提取的初至?xí)r間，得出如下結(jié)論：

1)特征函數(shù)可以很大程度上提高初至拾取所依賴的數(shù)據(jù)的質(zhì)量，突出初至波到來時的各種震相變化幅度。

2)初至自動拾取方法減小了人工拾取P波到時的不穩(wěn)定性及個體因素的差異性，自動拾取初至算法拾取的初至到時穩(wěn)定、耗時小，效率高。通過對地震事件中19道地震數(shù)據(jù)測試，人工拾取法平均耗時為12 min，自動拾取法平均耗時為5.68 s。

3)在自動拾取算法中，STA/LTA算法拾取的初至?xí)r間精度最高。通過對STA/LTA法、AIC法、能量比法、最大振幅法4種地震波初至?xí)r間對比分析，STA/LTA法與人工拾取法的結(jié)果平均差值為0.46 ms，能量比法為0.89 ms，AIC法為1.03 ms，最大振幅法為4.89 ms。

幾種初至自動拾取方法對于干擾大、信噪比低的數(shù)據(jù)，其精度均會降低。有待針對低信噪比數(shù)據(jù)，引入抗干擾能力更強的算法獲得更加高效、準確的拾取結(jié)果。