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

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

0 引言

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

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

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

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

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

(b) 工作示意圖

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

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

1 初至拾取方法

1.1 最大振幅法

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

1.2 能量比法

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

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

(1)

(2)

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

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

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

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即峰度赤池信息量準(zhǔn)則，由日本統(tǒng)計(jì)學(xué)家赤池弘次創(chuàng)立和發(fā)展； var表示序列的方差函數(shù)。

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

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

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

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ù)的長短時(shí)均比值等特征拾取初至[18-19]。當(dāng)信號到達(dá)時(shí)，STA要比LTA變化更快，相應(yīng)的STA/LTA值會出現(xiàn)明顯的增加，當(dāng)比值大于設(shè)定的閾值時(shí)，就可判定初至波到達(dá)[20-21]。該方法容易受到時(shí)窗長度和滑動步長的影響[22-23]。

(4)

(5)

(6)

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

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

圖4 STA/LTA算法流程圖

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

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

2 對比分析

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

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

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

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

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

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

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種不同的地震波初至自動拾取算法對同一組隧道超前探測地震記錄事件進(jìn)行初至拾取，可以得到19道地震數(shù)據(jù)的初至到時(shí)，見表2。

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

表2 地震波初至到時(shí)拾取表

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

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

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

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

選取不同隧道超前探測項(xiàng)目中記錄到的數(shù)據(jù)，以1次隧道超前探測產(chǎn)生的所有地震數(shù)據(jù)合稱為1次地震事件。事件1至事件4表示4次超前探測產(chǎn)生的地震數(shù)據(jù)，其原始數(shù)據(jù)如圖12所示。利用幾種算法繼續(xù)進(jìn)行初至拾取，將拾取的各個(gè)道到時(shí)與手動拾取的結(jié)果進(jìn)行差值，并計(jì)算所有地震事件拾取結(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ì)量，突出初至波到來時(shí)的各種震相變化幅度。

2)初至自動拾取方法減小了人工拾取P波到時(shí)的不穩(wěn)定性及個(gè)體因素的差異性，自動拾取初至算法拾取的初至到時(shí)穩(wěn)定、耗時(shí)小，效率高。通過對地震事件中19道地震數(shù)據(jù)測試，人工拾取法平均耗時(shí)為12 min，自動拾取法平均耗時(shí)為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ù)，引入抗干擾能力更強(qiáng)的算法獲得更加高效、準(zhǔn)確的拾取結(jié)果。