日本是第一个建立地震预警系统的国家,日本预警系统在2011年3月11日日本发生9.0级特大地震时发挥了作用,东京地区民众在地震S波最大震动到达之前收到地震报警信息,得以避险(陈会忠等,2011)。随着经济的快速发展和社会的不断进步,社会公众对地震速报提出了越来越高的要求。地震预警具有极其重要的社会效益,它不仅可以减轻人员伤亡和降低次生地震灾害的发生,还可为震后紧急救援和抢修提供第一手的信息(刘希强等,2009)。

快速确定地震参数是地震预警系统的关键技术环节之一,关系到预警时间的长短。现场地震预警系统利用P波进行预警可以获得更长的预警时间,即在预警的目标区安装地震计,利用P波比S波传播速度快的原理,用P波到达后数秒的数据来确定地震参数,对S波还没到达的地区通过电磁波发出预警,从而达到预警的目的(周彦文等,2010)。目前国外已经涌现出大量关于地震事件自动实时检测、P波初动拾取、震相类型识别方面的研究方法和成果(Allen,1978,1982; Bear,Kradolfer,1987; Sleeman,Van,1999; Bai,Kennett,2000; Diehl et al,2009; Kuperkoch et al,2010)。国内学者提出了基于分形分维(常旭,刘伊克,2002)、自回归模型(Zhang et al.,2003; 王海军等,2003; 杨配新等,2004; 王继等,2006)、人工神经网络法(Wang,Teng,1995; 张范民,李清河,1998; 李慧婷,黄文辉,2000; 王娟等,2004)、波形变化增长法(朱元清等,2002)、小波变换的主成分分析法(刘希强等,2000)、定尺度比方法(刘希强等,2002)、地震事件自动检测的EFGLP方法(刘希强等,2009)、初至波震相自动精确识别的TOC-AIC方法(刘希强等,2009)、基于幅值和频率的P波识别(毛燕等,2011)等来识别震相,取得了较好的效果。

Kuperkoch等(2010)根据高阶统计量(HOS)中的峰度(kurtosis)和偏斜度(skewness)与Akaike信息准则提出了一种地震事件和震相的自动检测和识别算法,并将该算法与Allen(1978)、Bear和Kradolfer(1987)的检测算法对自动获得P波震相读取值进行了比较研究。研究表明,高阶统计量方法显示出了更佳的识别效果,尤其峰度识别效果最佳,且分别使用偏斜度和峰度对3 000多个人工测量的P波震相的到时识别其平均偏差分别为(0.38±0.75)s和(0.26±0.64)s。笔者使用峰度方法对山东地震台网记录的单台垂向记录进行了区域地震事件自动检测和震相识别探索研究,同时应用AIC方法对峰度方法所确定的包含地震P波在内的时间窗内记录进行自动分析,探索进一步提高P波震相初至识别的技术和方法。

椭圆型带通滤波器属于无限脉冲响应滤波器,它在有限频率范围内存在传输零点和极点; 其通带和阻带都具有等波纹特性,通带、阻带逼近特性良好; 对于同样的性能要求,该滤波器比巴特沃斯、切比雪夫这两种滤波器所需用的阶数都低,而且它的过渡带比较窄。

图1a给出了阶数为7的椭圆形带通滤波器,其基本参数如下:通带波动系数为0.1,阻带衰减系数为30、数据采样率为100 p/s、通带低端截止频率5.5 Hz、通带高端截止频率7.5 Hz、带通阻带低端截止频率为5 Hz、带通阻带高端截止频率为8 Hz。从幅频特性可以看出,椭圆型带通滤波器具有通带、阻带逼近特性良好的特点。图1b给出的是由1.5、6.5和10.5 Hz正弦波组成的合成信号,其中6.5 Hz的信号起始点在2 s。图1c给出了根据上述滤波器进行数据滤波的结果,可

图1 椭圆形滤波器对合成模拟仿真信号的带通滤波结果
(a)椭圆形带通滤波器的幅频特征;(b)由3种频率
合成的仿真信号;(c)带通滤波信号与指定频段仿真信号的对比
Fig.1 Synthesis analog simulation signal filtered by elliptical band-pass filter (a)Amplitude-frequency characteristics of the elliptical band-pass filter;(b)Analog simulation signal synthesized by sine wave of three different frequency;(c)Comparison between band-pass filter signal and simulation signal with designated frequency

在处理区域地震事件(震中距小于120 km)的数据时,考虑到其主频分布范围,对椭圆型滤波器的设计参数如下:通带低端截止频率2 Hz、通带高端截止频率15 Hz、带通阻带低端截止频率为1.5 Hz、带通阻带高端截止频率为16 Hz。

随机变量X的数学期望表示为

E[X]=∫^+∞∫_-∞xp(x)dx.(1)

式中,p(x)是随机变量X的概率密度。

利用期望函数,随机变量X的k阶统计量m_k可表示为

m_k=E[(X-E[X])^k], k>1.(2)

式中,当k=2时表示X的方差,m₂表示X的取值与数学期望的偏离程度。若X比较集中,则m₂较小; 反之,若取值比较分散,则m₂较大。因此,m₂是刻画X取值分散程度的一个量,它是衡量X取值分散程度的一个尺度。

峰度是指频数分布的集中程度,也就是分布曲线的尖峭程度。随机变量的峰度的表达式为

K=(m₄)/((m₂)²).(3)

峰度K是四阶中心动差m₄除以标准差四次方。将四阶动差m₄除以标准差四次方,同样也是为了得到相对数,以便进行不同频数分配的峰度比较。当峰度K>3时,表示频数分布比正态分布更集中,分布呈尖峰状态,平均数的代表性更大; K<3时表示频数分布比正态分布更分散,分布呈平坦峰,平均数的代表性较小。

常用的震相自动识别技术是AR-AIC方法,其基本原理是假定地震波记录可分为若干个局部稳态的过程,每个过程满足自回归模型,并认为新的地震波震相到达前后的信号是两个不同的稳态过程(Sleeman,Van,1999)。设长度为N的地震波记录包含有背景信号和地震波信号,通过滑动变量k的不断改变,得到长度分别为[1,k-M]和[k-M+1,N-M-k]数据段的最大似然函数的最大值,进而得到AIC函数:

AIC(k)=(k-M)lg(σ²_1,max)+(N-M-k)lgσ²_2,max)+C₂.(4)

式中,k表示滑动变量(下同),M表示AR拟合过程的阶数,σ²_1,max和σ²_2,max表示最大似然函数的最大值,C₂是常数。AR过程的阶数需要通过实际记录并经多次试验计算得到,AR系数根据Yule-Walker方程确定(Haykin,1996)。当所选数据中包含有背景噪声和地震波信号时,AIC函数的最小值对应的时间就是地震震相初至。

Maeda(1985)提出了无需计算AR系数,直接根据地震图记录而得到AIC函数的方法,称为VAR-AIC方法。AIC函数的最小值对应的时间就是地震震相初至。AIC函数表示为

AIC(k)=klg(var(x[1,k]))+(N-k-1)lg(varx[k+1,N])).(5)

式中,var(x[1,k])和var(x[k+1,N])分别表示两个时间段内数据的方差。

VAR-AIC方法属于二阶统计量方法。二阶统计量(如二阶矩、相关函数、功率谱密度函数等)方法需要假设系统是最小相位系统,噪声为高斯白噪声,然而实际信号常常不满足这个假设。高阶统计量与二阶统计量相比具有适于检测和描述系统的非线性,在低信噪比条件下,可有效抑制加性高斯有色噪声的影响,检测和识别较弱的地震信号。用峰度进行精细震相识别的具体思路是用两个时间段数据的峰度作为特征函数代替公式(5)中的方差,本文称为Kurtosis-AIC方法。Kurtosis-AIC函数表示为

Kurtosis-AIC(k)=klg(1/k∑^k_j=1CF²_j)+(L-k+1)lg(1/(L-k+1)∑^l_j=kCF²_j).(6)

式中,L是特征函数CF的长度,k取值为从1到L。

识别震相初至的具体思路是,以峰度方法自动确定的地震事件触发判断结果为参考,选择一个包含有噪声和地震信号的时间窗,Kurtosis-AIC函数的最小值对应时刻即为震相初至。

本文对50次地震的198个单台垂直向记录进行了分析与研究,涉及的台站数量有53个(图2)。所选取单台记录的最大震中距为110.0 km、最小为6.0 km,震中距小于100 km的单台记录有190个,占总数的96%。50次地震中的最小震级为0.6级,最大震级为4.2级。

图2 研究区地震和台站分布图
Fig.2 Distribution of earthquakes and stations in study area

图4给出了数据滤波前后信噪比的比较及信噪比随震中距的变化。从图4可以看出,由于对信号进行了椭圆型带通滤波,有效提高了信噪比。图4a中数据滤波前信噪比的变化范围为-3.8～47.6 dB,均值为8.7 dB,标准差为9.4 dB,滤波后信噪比的变化范围为5.1～49.8 dB,均值为22.0 dB,

图3 梁山台记录的2009年6月15日22时27分
范县2.4级地震垂直向波形及事件触发参数变化
(a)垂直向记录;(b)(a)中虚框内记录;(c)对(b)
中波形的滤波结果;(d)数据的峰度曲线;
(e)数据的Kurtosis-AIC变化曲线
Fig.3 Vertical component waveform of Fanxian ML2.4 earthquake at 22:27 on Jun.15,2009 recorded by Liangshan Station and variation of event-trigger parameters (a)Vertical component waveform;(b)Waveform within the dotted frame in graph(a);(c)Filter result of the waveform in graph(b);(d)Kurtosis variation curve of data; (e)Kurtosis-AIC variation curve of data

图5给出了地震事件触发判断的峰度最低阈值和最大峰度值分布图。从图5a中可以看出,在对数据样本处理时,保证地震事件峰度触发阈值变化范围为2.5～8,均值为3.35,标准差为0.70。与之相对应,图5b中的峰度值的变化范围为4.2～691.9,其中位于4.2～8.0区间有3次,均值为137.0,标准差为118.5。为了保证地震事件既不漏报又不误报,需要在积累的大量样本中选择1个

图4 数据滤波前后信噪比的比较(a)及
信噪比随震中距的变化(b)
Fig.4 Signal-to-noise ratio comparison of data before and after the filtering(a)and signal-to-noise ratio change with epicenter distance(b)

图5 用于地震时间触发判断的峰度最低阈值(a)
和峰度最大值分布图(b)
Fig.5 Minimum threshold of Kurtosis used to judge the seismic time trigger(a)and distribution of the maximum Kurtosis value(b)

图6给出了Kurtosis-AIC震相自动识别结果与人机交互识别结果对比的误差分布及误差分布与观测数据信噪比的关系。由图6a可以看出,震相识别误差范围为-0.19～0.5 s,平均绝对误差为0.09 s,标准差为0.08 s。198个单台垂直向记录的处理结果中有127个误差在0.1 s内,占总数的64.14%; 有174个误差在0.17 s内,占88%; 有184个误差在0.2 s内,占93%。从图6a和图6b对比可以看出,震相识别误差较大的样本相对应的垂向信噪比多数比较低,反应了在人工分析中同样面临的当P波信号较弱时震相初至较难识别这一难题。但总体上Kurtosis-AIC方法对较弱的P波地震信号,也具有较好的分辨能力,达到了比较高的精度。

图6 Kurtosis-AIC震相自动识别结果与人机交互
识别结果对比的误差分布(a)及误差分布与
观测数据信噪比的关系(b)
Fig.6 Error distribution of comparision between the result of Kurtosis-AIC automatic recognition of seismic phase and that of man-machine interaction recognition(a)and the relationship between error distribution and signal-to-noise ratio of data(b)

本文利用山东地震台网记录的单台垂向记录和峰度方法进行了区域地震事件自动检测和震相识别方法的探索研究,为了减小噪声干扰又不失地方震记录信号,本文设计了保持信号的振幅和相位畸变最小的地震波形资料的椭圆型带通滤波方法。通过研究取得了以下结果:

(1)椭圆型滤波器提高了含噪声P波信号的信噪比。

(2)提出了一种基于直达P波信号进行区域地震事件实时检测的峰度方法,并用该方法研究得到了山东地震台网记录判断地震事件触发的经验阈值。表明选择合适的峰度经验阈值可同时有效降低地震事件的漏报率和误报率。

(3)提出了峰度与AIC方法相结合的精细震相识别新方法对198个单台垂向记录的分析研究表明,与人工识别震相到时的结果相比,根据Kurtosis-AIC震相自动识别方法得到的震相到时的平均绝对值误差为0.09±0.08 s,精度非常高。

(4)本文研究样本仍偏少,经验阈值和识别精度有待进一步检验。