基金项目:地震行业科研专项(20708031-5)、云南省科技厅运用基金面上项目(2008ZC160M)共同资助.
(1.昆明理工大学 建筑工程学院,昆明 650093; 2.云南省地震局,昆明 650224)
(1.The Architectural engineering Institute of Kunming University of Science and Technology,Kunming 650093,Yunnan,China)(2.Earthquake Administration of Yunnan Province,Kunming 650224,Yunnan,China)
备注
基金项目:地震行业科研专项(20708031-5)、云南省科技厅运用基金面上项目(2008ZC160M)共同资助.
震相到时的精确捡拾是地震定位的关键所在,是进行地震预警的前提。对云南测震台网的观测数据进行P波自动捡拾试验。用基于幅值和频率的P波识别方法和STA/LTA方法捡拾到的P波到时,与人工捡拾的结果比较接近,取得较好的结果; 用该方法对云南强震台网的部分强震记录的竖向资料进行P波到时自动识别,也获得了较好的结果。
The automatic,accurate detection of the arrival time of the seismic phase is crucial for the earthquake location,and is a prerequisite for earthquake early warning.We use the data acquired by Yunnan Seismic Networks to detect P-wave automatically by the two methods described in the article and get better results which are close to the artificial results.Better results are also gotten when we detect P-wave automatically for some strong,seismic,vertical component recordings by the same method.
引言
自然灾害是不可避免的,但是灾害造成的损失是可以减轻的。在实时地震监测台网的基础上建立地震预警系统和地震应急控制系统是近年来国际上的潮流,此类系统可以大大减轻地震灾害。然而,地震预警中的地震定位依赖于对震相到时的精确捡拾。震相识别是现代地震学研究中的重要课题,它是地球内部构造、地震定位、震源机制等一系列研究的基础。当人工捡拾震相满足不了人们的需求时,震相的自动识别就愈显重要。此工作开始于对大量地震记录的自动处理,随着实时地震学的发展,震相自动识别的精确性越来越引起地震学家的重视(Mao,Gubbins,1995; 李山有等,2004)。
目前,大多数的震相识别方法是将信号和噪声的不同特征作为震相到来的判据。常用的方法有能量变化分析、偏振分析和自回归方法(AR)等。其中,能量分析方法中最常用的是长短时间平均方法(STA/LTA方法)(Baer,Kradolfer,1987; Earle,Shearer,1994),该方法反应了幅值的瞬时变化,具有算法简单、速度快、便于实时处理等特点,被广泛应用于地震波的初动识别。此外,还有瞬时频率方法(武东坡,2004)。地震波的偏振分析可以提取地震波的偏振特征,因而常应用于地震波捡拾和震相判别。自回归方法则是假定可以把地震波划分为局部平稳段,并且在触发点前后是不相同的平稳过程(Sleeman,Eck,1999)。
近30年来我国数字地震观测台网迅速发展,特别是随着“十五”“中国数字地震观测网络项目”的完成以及在云南地区建成的两百多个数字强震动台的正式运行(崔建文等,2006),越来越多的地震观测数据可以用于预警相关内容的研究。云南地区的多震优势以及新建立的两百多个强震台站使得云南的地震预警系统也进入了试验研究阶段,毛燕等(2009)进行了强震数据的仿真处理。
1 理论
2 实际应用
笔者将以上识别方法应用到在云南地区获取的部分地震的强震记录中,使用竖向记录进行地震P波初动的自动识别,取得了较好的结果。表1给出了地震P波到时的自动识别结果和人工识别结果。从表中可以看出,人工识别和自动识别差值的平均值为0.030 s,误差很小,识别结果较好。图6、图7给出了2007年6月3日宁洱MS6.4地震正兴台和2008年8月20日盈江MS5.0地震梁河台记录的STA/LTA曲线。
图6 2007年6月3日宁洱MS6.4地震正兴台记录波形自动识别结果(a)原始记录;(b)特征函数
Fig.6 Automatic identification of the recordings of the Ning'er MS6.4 earthquake at June 3,2007 recorded by Zhengxing Station(a)original waveform;(b)eigenvalue function图7 2008年8月20日盈江MS5.0地震梁河台记录波形自动识别结果(a)原始记录;(b)特征函数
Fig.7 Automatic identification of the Yingjiang MS5.0 earthquake at August 20,2008 recorded by Lianghe Station(a)original waveform;(b)eigenvalue function3 结论
通过对沧源台的地震波记录进行基于幅值和频率的P波识别和STA/LTA算法的P波到时初步确定,再用最小二乘法和AIC准则进行精确确定,可以看出这两种方法均具有较准确的判定结果。
为了验证这两种方法的实用性,笔者运用这两种方法对云南强震台站记录到的竖向观测资料进行P波初动的自动识别,从表1可以看出,这两种方法对强震记录的识别比较准确,人工识别和自动识别的差别较小,平均差值为0.030 s。随着云南省越来越多观测台站的投入运行,在以后的地震预警中,可以用该方法来进行P波捡拾,推进地震预警系统的发展。
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为信号的采样间隔,显然有下列式子成立:
x(n-1)=Acos[ω(n-1)Δt+φ],(1)
x(n+1)=Acos[ω(n+1)Δt+φ].(2)
由三角变换得
x(n-1)+x(n+1)=2x(n)cosωΔt,(3)
则令
ω(n)=1/(Δt)arccos(x(n-1)+x(n+1))/(2x(n)).(4)
由于上式是在a(t)和φ(t)为常量的假设下得出的,所以在实际应用中作如下改进,令
Xg(n)=αXg(n-1)+[x(n-1)+x(n+1)]2,(5)
X(n)=αX(n-1)+[2x(n)]2.(6)
式中0≤α<1。
瞬时频率f(n)=1/(Δt)arccos((Xg(n))/(X(n)))1/2.(7)
再计算给定时间窗内(这里取为0.5 s)的峰值平均值X^-pk和瞬时频率平均值fk,最后利用下式计算特征函数Pk:
Pk={X^-pk·(f^-k)/(f^-k-1)f^-k/f^-k-1>1,
X^-pk·(f^-k-1)/(f^-k)f^-k/f^-k-1<1.(8)
当P波到达时,在Pk曲线上会产生一个跳跃点,第一个跳跃点位即为P波到时点。
图1 2007年6月3日10时48分宁洱MS 5.1地震沧源台记录的P波捡拾(a)仿真波形;(b)特征函数
Fig.1 Pick-up of the P-wave of the recordings of the Ning'er MS 5.1 earthquake at 10:48,June 3,2007 recorded by Cangyuan station(a)emulated waveform;(b)eigenvalue function图2 2007年6月3日10时48分宁洱MS 5.1地震沧源台记录的P波捡拾(用未仿真波形)(a)仿真波形;(b)特征函数
Fig.2 Pick-up of the P-wave of the recordings of the Ning'er MS 5.1 earthquake at 10:48,June 3,2007 recorded by Cangyuan station(no simulation)(a)emulated waveform;(b)eigenvalue function图3 2007年6月3日05时34分宁洱MS 6.4地震禄劝台记录的P波捡拾(a)仿真波形;(b)特征函数
Fig.3 Pick-up of the P-wave of the recordings of the Ning'er MS 6.4 earthquake at 05:34,June 3,2007 recorded by Luquan station(a)original waveform;(b)eigenvalue function1.2 长短时间平均方法(STA/LTA方法)1.2.1 原理STA/LTA算法原理是用STA(信号短时平均值)和LTA(信号长时平均值)的比值来反映信号水平或能量的变化。当信号到达时,STA要比LTA变化的快,相应的STA/LTA值会有一个明显的增加,当比值大于设定的阈值时,就可判定有一震相到达。
在计算中,选用的特征函数如下式所示:
CF(i)=X(i)2+[X(i)-X(i-1)]2.(9)
STA/LTA采用递归计算公式:
STAi=STAi-1+(CF(i)-STAi-1)/(Nsta),(10)
LTAi=LTAi-1+(CF(i-Nsta-1)-LTAi-1)/(Nlta).(11)
其中: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分宁洱MS 5.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'er MS5.1 earthquake at 10:48,June 3,2007 recorded by Cangyauan station(a)original waveform;(b)eigenvalue function; (c)eigenvalue waveform1.2.2 精确捡拾(1)线性最小二乘法精确估计P波到时
将P波触发点回推一个时间窗,在此取长为1 s的时间窗,设定阈值标准THR(计算时设为10)。确定大于2×ASN(ASN即该窗内的平均STA/LTA值)小于THR的拟合目标点;
该事件目标点间存在极大值点,进行目标点处理,保证拟合目标点在波峰一侧。利用最小二乘法对这些目标点进行一元线性拟合,得到的拟合直线和噪声均值线的交点就是近似的P波到时(图5)。
(2)AIC准则精确确定P波到时
日本学者Akaike(1973)提出一个基本信息量的定阶准则——AIC准则。王继等(2006)提出可由地震波形数据直接计算AIC函数,即对地震记录x(i)(i=1,2,….L)来说,AIC检测器定义为:
AIC(k)=k·log{var(x[1,k])}+(L-k-1)·
log{var(x[k+1,L])}.(12)
其中,k的范围是窗口内所有的采样点。震相到时对应于AIC的最小值。先对P波进行STA/LTA粗略捡拾,然后在找到的触发点向前推100个点,后推10个点,即对采样率为50 Hz的该记录前推2 s,后推0.2 s。用AIC准则计算的触发点为第86.96 s处,人工确定的触发点在86.98 s处,误差为0.02 s。
- 崔建文,高东,李世成,等.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.
- Mao W,Gubbins D.1995.Simultaneous determination of time delays and stacking weights in seismic array beamforming[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 seismogram recordings[J].Phys.Earth Planet.Interiors,113,265-275.