毛 燕,崔建文,郑定昌,李正光,卢吉高

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

地震记录的 P波自动捡拾＊

毛 燕1,2,崔建文2,郑定昌2,李正光2,卢吉高2

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

震相到时的精确捡拾是地震定位的关键所在,是进行地震预警的前提。对云南测震台网的观测数据进行 P波自动捡拾试验。用基于幅值和频率的 P波识别方法和 STA/LTA方法捡拾到的 P波到时,与人工捡拾的结果比较接近,取得较好的结果;用该方法对云南强震台网的部分强震记录的竖向资料进行 P波到时自动识别,也获得了较好的结果。

自动捡拾;P波;地震记录

0 前言

自然灾害是不可避免的,但是灾害造成的损失是可以减轻的。在实时地震监测台网的基础上建立地震预警系统和地震应急控制系统是近年来国际上的潮流,此类系统可以大大减轻地震灾害。然而,地震预警中的地震定位依赖于对震相到时的精确捡拾。震相识别是现代地震学研究中的重要课题,它是地球内部构造、地震定位、震源机制等一系列研究的基础。当人工捡拾震相满足不了人们的需求时,震相的自动识别就愈显重要。此工作开始于对大量地震记录的自动处理,随着实时地震学的发展,震相自动识别的精确性越来越引起地震学家的重视(Mao,Gubbins,1995;李山有等,2004)。

目前,大多数的震相识别方法是将信号和噪声的不同特征作为震相到来的判据。常用的方法有能量变化分析、偏振分析和自回归方法(AR)等。其中,能量分析方法中最常用的是长短时间平均方法(STA/LT A方法)(Baer,Kradolfer,1987;Earle,Shearer,1994),该方法反应了幅值的瞬时变化,具有算法简单、速度快、便于实时处理等特点,被广泛应用于地震波的初动识别。此外,还有瞬时频率方法(武东坡,2004)。地震波的偏振分析可以提取地震波的偏振特征,因而常应用于地震波捡拾和震相判别。自回归方法则是假定可以把地震波划分为局部平稳段,并且在触发点前后是不相同的平稳过程(Sleeman,Eck,1999)。

近 30年来我国数字地震观测台网迅速发展,特别是随着“十五”“中国数字地震观测网络项目”的完成以及在云南地区建成的两百多个数字强震动台的正式运行(崔建文等,2006),越来越多的地震观测数据可以用于预警相关内容的研究。云南地区的多震优势以及新建立的两百多个强震台站使得云南的地震预警系统也进入了试验研究阶段,毛燕等(2009)进行了强震数据的仿真处理。

1 理论

1.1 基于幅值和频率的 P波识别

P波到达时在记录波形中会产生幅值和频率的变化,所以用于识别 P波震相的特征参量必须能够反映这些变化趋势,才能较为准确地识别 P波震相。

对于数字信号 x(t)=a(t)cosφ(t),在微小的时间间隔内,a(t)和φ(t)可以看作常量,即可以写为 x(n)=Acos(ωnΔt+φ)的形式。式中,n=1,2,…,n,Δt为信号的采样间隔,显然有下列式子成立:

由三角变换得

则令

由于上式是在 a(t)和φ(t)为常量的假设下得出的,所以在实际应用中作如下改进,令

式中 0≤α＜1。

当 P波到达时,在 Pk曲线上会产生一个跳跃点,第一个跳跃点位即为 P波到时点。

图1 2007年6月 3日 10时 48分宁洱 MS5.1地震沧源台记录的 P波捡拾(a)仿真波形;(b)特征函数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)特征函数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)特征函数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波初动点出现在 Pk曲线的跳跃点处,即 86.9 s处。禄劝台记录仪器的采样率是50 Hz,由图 3可以看出,P波初动点出现在 Pk曲线的跳跃点处,即 35 s处。

1.2 长短时间平均方法(STA/LTA方法)

1.2.1 原理

STA/LTA算法原理是用 STA(信号短时平均值)和LTA(信号长时平均值)的比值来反映信号水平或能量的变化。当信号到达时,STA要比LTA变化的快,相应的 STA/LTA值会有一个明显的增加,当比值大于设定的阈值时,就可判定有一震相到达。

在计算中,选用的特征函数如下式所示:

STA/LTA采用递归计算公式:

其中:STAi和 LTAi分别为信号在 i时刻的短时平均值和长时平均值,CF(i)为信号在 i时刻的特征函数值,Nsta和 Nlta分别为短时平均值和长时平均值时间窗包含的记录点数。

特征函数的触发阈值设定为 10,若大于此阈值,即可判定地震事件的发生。短窗取 0.2 s,长窗取 10 s。

利用该方法捡拾到的 P波也在 86.90 s处触发(图4)。

图4 2007年6月 3日 10时 48分宁洱 MS5.1地震沧源台 STA/LTA方法捡拾到的 P波(a)原始记录;(b)特征函数;(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波触发点回推一个时间窗,在此取长为1 s的时间窗,设定阈值标准 THR(计算时设为10)。确定大于 2×ASN(ASN即该窗内的平均STA/LTA值)小于 THR的拟合目标点;

该事件目标点间存在极大值点,进行目标点处理,保证拟合目标点在波峰一侧。利用最小二乘法对这些目标点进行一元线性拟合,得到的拟合直线和噪声均值线的交点就是近似的 P波到时(图5)。

图5 精确确定 P波触发点Fig.5 Accurate determination of the P-wave triggerpoint

由图 5可知,通过最小二乘法精确确定的 P波触发点在 86.94 s处,人工确定的触发点在 86.98 s处,误差为 0.04 s。

(2)A IC准则精确确定 P波到时

日本学者 Akaike(1973)提出一个基本信息量的定阶准则——A IC准则。王继等(2006)提出可由地震波形数据直接计算 A IC函数,即对地震记录 x(i)(i=1,2,….L)来说,A IC检测器定义为:

其中,k的范围是窗口内所有的采样点。震相到时对应于 AIC的最小值。先对 P波进行 STA/LTA粗略捡拾,然后在找到的触发点向前推 100个点,后推 10个点,即对采样率为 50 Hz的该记录前推2 s,后推 0.2 s。用 A I C准则计算的触发点为第86.96 s处,人工确定的触发点在 86.98 s处,误差为 0.02 s。

2 实际应用

笔者将以上识别方法应用到在云南地区获取的部分地震的强震记录中,使用竖向记录进行地震 P波初动的自动识别,取得了较好的结果。表 1给出了地震 P波到时的自动识别结果和人工识别结果。从表中可以看出,人工识别和自动识别差值的平均值为 0.030 s,误差很小,识别结果较好。图 6、图 7给出了 2007年6月 3日宁洱MS6.4地震正兴台和 2008年8月 20日盈江MS5.0地震梁河台记录的 STA/LTA曲线。

表1 P波到时的人工识别和自动识别结果比较Tab.1 Automatic and manual identified results of arrival time of P wave

图6 2007年6月 3日宁洱 MS6.4地震正兴台记录波形自动识别结果(a)原始记录;(b)特征函数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地震梁河台记录波形自动识别结果(a)原始记录;(b)特征函数Fig.7 Automatic identification of the Yingjiang MS5.0 earthquake atAugust 20,2008recorded byLianghe Station(a)originalwaveform;(b)eigenvalue function

3 结论

通过对沧源台的地震波记录进行基于幅值和频率的 P波识别和 STA/LTA算法的 P波到时初步确定,再用最小二乘法和A IC准则进行精确确定,可以看出这两种方法均具有较准确的判定结果。

为了验证这两种方法的实用性,笔者运用这两种方法对云南强震台站记录到的竖向观测资料进行 P波初动的自动识别,从表 1可以看出,这两种方法对强震记录的识别比较准确,人工识别和自动识别的差别较小,平均差值为 0.030 s。随着云南省越来越多观测台站的投入运行,在以后的地震预警中,可以用该方法来进行 P波捡拾,推进地震预警系统的发展。

崔建文,高东,李世成,等 .2006.新的云南数字强震动观测台网[J].地震研究,29(增刊):453-458.

李山有,金星,马强,等 .2004.地震预警系统与智能应急控制系统研究[J].世界地震工程,20(4):21-26.

毛燕,崔建文,虎雄林,等 .2009.地震数据的仿真处理 [J].地震研究,32(增刊):461-463.

王继,陈九辉,刘启元,等 .2006.流动地震台阵观测初至震相的自动检测[J].地震学报,28(1):42-51.

武东坡 .2004.震相识别的实时方法研究[D].哈尔滨:中国地震局工程力学研究所 .

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.

地震行业科研专项(20708031-5)、云南省科技厅运用基金面上项目(2008ZC160M)共同资助 .