地震記錄的 P波自動撿拾＊

毛 燕,崔建文,鄭定昌,李正光,盧吉高

(1.昆明理工大學(xué)建筑工程學(xué)院,昆明 650093;2.云南省地震局,昆明 650224)

地震記錄的 P波自動撿拾＊

毛 燕1,2,崔建文2,鄭定昌2,李正光2,盧吉高2

(1.昆明理工大學(xué)建筑工程學(xué)院,昆明 650093;2.云南省地震局,昆明 650224)

震相到時的精確撿拾是地震定位的關(guān)鍵所在,是進(jìn)行地震預(yù)警的前提。對云南測震臺網(wǎng)的觀測數(shù)據(jù)進(jìn)行 P波自動撿拾試驗。用基于幅值和頻率的 P波識別方法和 STA/LTA方法撿拾到的 P波到時,與人工撿拾的結(jié)果比較接近,取得較好的結(jié)果;用該方法對云南強(qiáng)震臺網(wǎng)的部分強(qiáng)震記錄的豎向資料進(jìn)行 P波到時自動識別,也獲得了較好的結(jié)果。

自動撿拾;P波;地震記錄

0 前言

自然災(zāi)害是不可避免的,但是災(zāi)害造成的損失是可以減輕的。在實時地震監(jiān)測臺網(wǎng)的基礎(chǔ)上建立地震預(yù)警系統(tǒng)和地震應(yīng)急控制系統(tǒng)是近年來國際上的潮流,此類系統(tǒng)可以大大減輕地震災(zāi)害。然而,地震預(yù)警中的地震定位依賴于對震相到時的精確撿拾。震相識別是現(xiàn)代地震學(xué)研究中的重要課題,它是地球內(nèi)部構(gòu)造、地震定位、震源機(jī)制等一系列研究的基礎(chǔ)。當(dāng)人工撿拾震相滿足不了人們的需求時,震相的自動識別就愈顯重要。此工作開始于對大量地震記錄的自動處理,隨著實時地震學(xué)的發(fā)展,震相自動識別的精確性越來越引起地震學(xué)家的重視(Mao,Gubbins,1995;李山有等,2004)。

目前,大多數(shù)的震相識別方法是將信號和噪聲的不同特征作為震相到來的判據(jù)。常用的方法有能量變化分析、偏振分析和自回歸方法(AR)等。其中,能量分析方法中最常用的是長短時間平均方法(STA/LT A方法)(Baer,Kradolfer,1987;Earle,Shearer,1994),該方法反應(yīng)了幅值的瞬時變化,具有算法簡單、速度快、便于實時處理等特點(diǎn),被廣泛應(yīng)用于地震波的初動識別。此外,還有瞬時頻率方法(武東坡,2004)。地震波的偏振分析可以提取地震波的偏振特征,因而常應(yīng)用于地震波撿拾和震相判別。自回歸方法則是假定可以把地震波劃分為局部平穩(wěn)段,并且在觸發(fā)點(diǎn)前后是不相同的平穩(wěn)過程(Sleeman,Eck,1999)。

近 30年來我國數(shù)字地震觀測臺網(wǎng)迅速發(fā)展,特別是隨著“十五”“中國數(shù)字地震觀測網(wǎng)絡(luò)項目”的完成以及在云南地區(qū)建成的兩百多個數(shù)字強(qiáng)震動臺的正式運(yùn)行(崔建文等,2006),越來越多的地震觀測數(shù)據(jù)可以用于預(yù)警相關(guān)內(nèi)容的研究。云南地區(qū)的多震優(yōu)勢以及新建立的兩百多個強(qiáng)震臺站使得云南的地震預(yù)警系統(tǒng)也進(jìn)入了試驗研究階段,毛燕等(2009)進(jìn)行了強(qiáng)震數(shù)據(jù)的仿真處理。

1 理論

1.1 基于幅值和頻率的 P波識別

P波到達(dá)時在記錄波形中會產(chǎn)生幅值和頻率的變化,所以用于識別 P波震相的特征參量必須能夠反映這些變化趨勢,才能較為準(zhǔn)確地識別 P波震相。

對于數(shù)字信號 x(t)=a(t)cosφ(t),在微小的時間間隔內(nèi),a(t)和φ(t)可以看作常量,即可以寫為 x(n)=Acos(ωnΔt+φ)的形式。式中,n=1,2,…,n,Δt為信號的采樣間隔,顯然有下列式子成立:

由三角變換得

則令

由于上式是在 a(t)和φ(t)為常量的假設(shè)下得出的,所以在實際應(yīng)用中作如下改進(jìn),令

式中 0≤α＜1。

當(dāng) P波到達(dá)時,在 Pk曲線上會產(chǎn)生一個跳躍點(diǎn),第一個跳躍點(diǎn)位即為 P波到時點(diǎn)。

圖1 2007年6月 3日 10時 48分寧洱 MS5.1地震滄源臺記錄的 P波撿拾(a)仿真波形;(b)特征函數(shù)Fig.1 Pick-up of the P-wave of the recordings of the Ning'erMS5.1 earthquake at 10:48,June 3,2007 recorded by Cangyuan station(a)emulated waveform;(b)eigenvalue function

圖2 2007年6月 3日 10時 48分寧洱 MS5.1地震滄源臺記錄的 P波撿拾(用未仿真波形)(a)仿真波形;(b)特征函數(shù)Fig.2 Pick-up of the P-wave of the recordings of the Ning'erMS5.1 earthquake at 10:48,June 3,2007 recorded by Cangyuan station(no s imulation)(a)emulated wavefo rm;(b)eigenvalue function

圖3 2007年6月 3日 05時 34分寧洱 MS6.4地震祿勸臺記錄的 P波撿拾(a)仿真波形;(b)特征函數(shù)Fig.3 Pick-up of the P-wave of the recordings of the Ning'erMS6.4 earthquake at 05:34,June 3,2007 recorded byLuquan station(a)originalwaveform;(b)eigenvalue function

滄源臺記錄儀器的采樣率是 50 Hz,由圖 1、圖 2可以看出,P波初動點(diǎn)出現(xiàn)在 Pk曲線的跳躍點(diǎn)處,即 86.9 s處。祿勸臺記錄儀器的采樣率是50 Hz,由圖 3可以看出,P波初動點(diǎn)出現(xiàn)在 Pk曲線的跳躍點(diǎn)處,即 35 s處。

1.2 長短時間平均方法(STA/LTA方法)

1.2.1 原理

STA/LTA算法原理是用 STA(信號短時平均值)和LTA(信號長時平均值)的比值來反映信號水平或能量的變化。當(dāng)信號到達(dá)時,STA要比LTA變化的快,相應(yīng)的 STA/LTA值會有一個明顯的增加,當(dāng)比值大于設(shè)定的閾值時,就可判定有一震相到達(dá)。

在計算中,選用的特征函數(shù)如下式所示:

STA/LTA采用遞歸計算公式:

其中:STAi和 LTAi分別為信號在 i時刻的短時平均值和長時平均值,CF(i)為信號在 i時刻的特征函數(shù)值,Nsta和 Nlta分別為短時平均值和長時平均值時間窗包含的記錄點(diǎn)數(shù)。

特征函數(shù)的觸發(fā)閾值設(shè)定為 10,若大于此閾值,即可判定地震事件的發(fā)生。短窗取 0.2 s,長窗取 10 s。

利用該方法撿拾到的 P波也在 86.90 s處觸發(fā)(圖4)。

圖4 2007年6月 3日 10時 48分寧洱 MS5.1地震滄源臺 STA/LTA方法撿拾到的 P波(a)原始記錄;(b)特征函數(shù);(c)特征波形Fig.4 Pick-up through STA/LTA method of the P-wave of the recordings of the Ning'erMS5.1 earthquake at 10:48,June 3,2007 recorded by Cangyauan station(a)originalwaveform;(b)eigenvalue function;(c)eigenvalue wavefo rm

1.2.2 精確撿拾

(1)線性最小二乘法精確估計 P波到時

將 P波觸發(fā)點(diǎn)回推一個時間窗,在此取長為1 s的時間窗,設(shè)定閾值標(biāo)準(zhǔn) THR(計算時設(shè)為10)。確定大于 2×ASN(ASN即該窗內(nèi)的平均STA/LTA值)小于 THR的擬合目標(biāo)點(diǎn);

該事件目標(biāo)點(diǎn)間存在極大值點(diǎn),進(jìn)行目標(biāo)點(diǎn)處理,保證擬合目標(biāo)點(diǎn)在波峰一側(cè)。利用最小二乘法對這些目標(biāo)點(diǎn)進(jìn)行一元線性擬合,得到的擬合直線和噪聲均值線的交點(diǎn)就是近似的 P波到時(圖5)。

圖5 精確確定 P波觸發(fā)點(diǎn)Fig.5 Accurate determination of the P-wave triggerpoint

由圖 5可知,通過最小二乘法精確確定的 P波觸發(fā)點(diǎn)在 86.94 s處,人工確定的觸發(fā)點(diǎn)在 86.98 s處,誤差為 0.04 s。

(2)A IC準(zhǔn)則精確確定 P波到時

日本學(xué)者 Akaike(1973)提出一個基本信息量的定階準(zhǔn)則——A IC準(zhǔn)則。王繼等(2006)提出可由地震波形數(shù)據(jù)直接計算 A IC函數(shù),即對地震記錄 x(i)(i=1,2,….L)來說,A IC檢測器定義為:

其中,k的范圍是窗口內(nèi)所有的采樣點(diǎn)。震相到時對應(yīng)于 AIC的最小值。先對 P波進(jìn)行 STA/LTA粗略撿拾,然后在找到的觸發(fā)點(diǎn)向前推 100個點(diǎn),后推 10個點(diǎn),即對采樣率為 50 Hz的該記錄前推2 s,后推 0.2 s。用 A I C準(zhǔn)則計算的觸發(fā)點(diǎn)為第86.96 s處,人工確定的觸發(fā)點(diǎn)在 86.98 s處,誤差為 0.02 s。

2 實際應(yīng)用

筆者將以上識別方法應(yīng)用到在云南地區(qū)獲取的部分地震的強(qiáng)震記錄中,使用豎向記錄進(jìn)行地震 P波初動的自動識別,取得了較好的結(jié)果。表 1給出了地震 P波到時的自動識別結(jié)果和人工識別結(jié)果。從表中可以看出,人工識別和自動識別差值的平均值為 0.030 s,誤差很小,識別結(jié)果較好。圖 6、圖 7給出了 2007年6月 3日寧洱MS6.4地震正興臺和 2008年8月 20日盈江MS5.0地震梁河臺記錄的 STA/LTA曲線。

表1 P波到時的人工識別和自動識別結(jié)果比較Tab.1 Automatic and manual identified results of arrival time of P wave

圖6 2007年6月 3日寧洱 MS6.4地震正興臺記錄波形自動識別結(jié)果(a)原始記錄;(b)特征函數(shù)Fig.6 Automatic identification of the recordings of the Ning'erMS6.4 earthquake at June 3,2007 recorded by Zhengxing Station(a)originalwaveform;(b)eigenvalue function

圖7 2008年8月 20日盈江 MS5.0地震梁河臺記錄波形自動識別結(jié)果(a)原始記錄;(b)特征函數(shù)Fig.7 Automatic identification of the Yingjiang MS5.0 earthquake atAugust 20,2008recorded byLianghe Station(a)originalwaveform;(b)eigenvalue function

3 結(jié)論

通過對滄源臺的地震波記錄進(jìn)行基于幅值和頻率的 P波識別和 STA/LTA算法的 P波到時初步確定,再用最小二乘法和A IC準(zhǔn)則進(jìn)行精確確定,可以看出這兩種方法均具有較準(zhǔn)確的判定結(jié)果。

為了驗證這兩種方法的實用性,筆者運(yùn)用這兩種方法對云南強(qiáng)震臺站記錄到的豎向觀測資料進(jìn)行 P波初動的自動識別,從表 1可以看出,這兩種方法對強(qiáng)震記錄的識別比較準(zhǔn)確,人工識別和自動識別的差別較小,平均差值為 0.030 s。隨著云南省越來越多觀測臺站的投入運(yùn)行,在以后的地震預(yù)警中,可以用該方法來進(jìn)行 P波撿拾,推進(jìn)地震預(yù)警系統(tǒng)的發(fā)展。

崔建文,高東,李世成,等 .2006.新的云南數(shù)字強(qiáng)震動觀測臺網(wǎng)[J].地震研究,29(增刊):453-458.

李山有,金星,馬強(qiáng),等 .2004.地震預(yù)警系統(tǒng)與智能應(yīng)急控制系統(tǒng)研究[J].世界地震工程,20(4):21-26.

毛燕,崔建文,虎雄林,等 .2009.地震數(shù)據(jù)的仿真處理 [J].地震研究,32(增刊):461-463.

王繼,陳九輝,劉啟元,等 .2006.流動地震臺陣觀測初至震相的自動檢測[J].地震學(xué)報,28(1):42-51.

武東坡 .2004.震相識別的實時方法研究[D].哈爾濱:中國地震局工程力學(xué)研究所 .

Akaike H.1973.Information theory and an extension of the maximum likelihood principle[C]//the 2nd International Symposium on Information Theory.Akademia Kiado,Budapest,Hungary:267-281.

Baer M,Kradolfer U.1987.An automatic phase picker for local and teleseismic events[J].BSSA.77(4):1437-1445.

Earle P S,Shearer P M.1994.Characterization of global seismograms using an automatic-picking algorithm[J].BSSA.84(2):366-376.

MaoW,GubbinsD.1995.Simultaneous deter mination of time delays and stackingweights in seis mic array beamfor ming[J].Geophysics,60(2),491-502.

Sleeman R,Van E T.1999.Robust automatic P-phase picking:an on-line implementation in the analysis of broadband seis mogram recordings[J].Phys.Earth Planet.Interiors,113,265-275.

Automatic P-wave Detection of the Earthquake Recordings

MAO Yan1,2,CU I Jian-w en2,ZHENG D ing-chang2,L I Zheng-guang2,LU Ji-gao2
(1.The Architectural engineering Institute of Kunm ing University of Science and Technology,Kunm ing650093,Yunnan,China)(2.Earthquake Adm inistration of Yunnan Province,Kunm ing650224,Yunnan,China)

The autom atic,accurate detection of the arrival tim e of the seism ic phase is crucial for the earthquake location,and is a prerequisite for earthquake early w arning.W e use the data acquired by Yunnan Seism ic N etw orks to detect P-w ave autom atically by the tw o m ethods described in the article and get better results w hich are close to the artificial results.B etter results are also gotten w hen w e detect P-w ave autom atically for som e strong,seism ic,vertical component recordings by the sam e m ethod.

autom atic detection;P-w ave;earthquake records

P315.3+1

A

1000-0666(2011)01-0047-05

2010-03-19.

地震行業(yè)科研專項(20708031-5)、云南省科技廳運(yùn)用基金面上項目(2008ZC160M)共同資助 .