峰度和AIC方法在区域地震事件和直达P波初动自动识别方面的应用*

赵大鹏，刘希强，李 红，周彦文

(1.中国地震局兰州地震研究所，甘肃兰州730000;2.山东地震局，山东济南250014)

峰度和AIC方法在区域地震事件和直达P波初动自动识别方面的应用*

赵大鹏1，刘希强2，李 红2，周彦文2

(1.中国地震局兰州地震研究所，甘肃兰州730000;2.山东地震局，山东济南250014)

提出了一种基于直达P波信号的峰度和Kurtosis-AIC方法进行区域地震事件实时检测和直达P波初动精细识别的新方法，并应用于山东地震台网记录的地震波资料处理.结果表明:(1)应用峰度方法能够有效识别出地震事件，可有效减低地震事件的错误报警率和漏报率;(2)与人工识别震相到时的结果相比，根据Kurtosis-AIC震相自动识别方法得到的震相到时的平均绝对值误差为 (0.09±0.08)s。

峰度方法;Kurtosis-AIC方法;地震事件识别;震相识别

0 引言

日本是第一个建立地震预警系统的国家，日本预警系统在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波震相初至识别的技术和方法。

1 地震波形资料的椭圆型带通滤波方法

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

图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 bandpass 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。

2 区域地震事件触发判断的峰度方法

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

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

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

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

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

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

3 初至波震相自动精确识别的Kurtosis-AIC方法

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

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

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

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

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

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

4 资料和结果分析

本文对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

采用逐点滑动的方式对数据的峰度进行实时计算和分析，计算时间窗长为10 s，步长为1个采样点。为了总结地震事件触发阈值的分布规律，选择大于噪声段最大峰度值作为每次地震事件触发的最低阈值，峰度超过阈值所对应的时间初定为P波震相初至。然后应用时间段的峰度数据进行Kurtosis-AIC方法分析，其最小值对应时刻即为精细识别的P波震相。图3给出了梁山台记录的2009年6月15日22时27分范县2.4级地震垂直向波形、事件触发时刻及应用峰度和Kurtosis-AIC方法对初至震相到时的拾取结果。当选取峰度阈值等于3时，由峰度方法得到的地震事件触发时刻为49.69 s，应用Kurtosis-AIC方法得到的P波震相初至结果为49.68 s，与人机交互拾取结果49.64 s相比，震相识别的误差为0.04 s，识别结果非常理想。

图4给出了数据滤波前后信噪比的比较及信噪比随震中距的变化。从图4可以看出，由于对信号进行了椭圆型带通滤波，有效提高了信噪比。图4a中数据滤波前信噪比的变化范围为－3.8～47.6 dB，均值为8.7 dB，标准差为9.4 dB，滤波后信噪比的变化范围为5.1～49.8 dB，均值为22.0 dB，标准差为9.0 dB，可见对观测数据进行滤波后信噪比得到了明显提高。图4b表明对发生在山东地区2级左右小地震进行滤波处理后的信噪比随震中距的变化并未出现明显衰减的特点。

图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

图4 数据滤波前后信噪比的比较 (a)及信噪比随震中距的变化 (b)Fig.4 Signal-to-noise ratio comparison of data before and after the filtering(a)and signal-tonoise 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)

图5给出了地震事件触发判断的峰度最低阈值和最大峰度值分布图。从图5a中可以看出，在对数据样本处理时，保证地震事件峰度触发阈值变化范围为2.5～8，均值为3.35，标准差为0.70。与之相对应，图5b中的峰度值的变化范围为4.2～691.9，其中位于4.2～8.0区间有3次，均值为137.0，标准差为118.5。为了保证地震事件既不漏报又不误报，需要在积累的大量样本中选择1个经验阈值。如果本文选定峰度阈值为8，研究样本中有3个单台记录出现漏报，这3次漏报地震的震级都为1.3级，漏报率为1.5%，误报率为0。上述出现漏报资料的信噪比都不高，分别是3.9 dB、1.08 dB和－1.7 dB。如果设置峰度阈值为8，并且地方震的单台垂向记录数据满足信噪比大于4 dB，则较大地震不会出现漏报和误报。

图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)

5 结论

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

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

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

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

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

常旭，刘伊克.2002.地震记录的广义分维及其应用［J］.地球物理学报，45(6):839－846.

陈会忠，侯燕燕，何加勇，等.2011.日本地震预警系统日趋完善［J］.国际地震动态，(4):10－15.

李慧婷，黄文辉.2000.应用人工神经网络方法识别近震与远震［J］.华南地震，20(4):71－75.

刘希强，周蕙兰，曹文海，等.2002.高斯线调频小波变换极其在地震震相识别中的应用［J］.地震学报，24(6):607－616.

刘希强，周蕙兰，沈萍，等.2000.用于三分向记录震相识别的小波变换方法［J］.地震学报，22(2):125－131.

刘希强，周彦文，曲均浩，等.2009.应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别［J］.地震学报，31(3):260－271.

毛燕，崔建文，郑定昌，等.2011.地震记录的P波自动捡拾［J］.地震研究，34(1):47－51.

王海军，靳平，刘贵忠，等.2003.区域震相初至估计［J］.西北地震学报，25(4):298－303.

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

王娟，刘俊民，范万春.2004.神经网络在震相识别中的应用［J］.现代电子技术，27(8):35－37.

杨配新，邓存华，刘希强，等.2004.数字化地震记录震相自动识别的方法研究［J］.地震研究，27(4):308－313.

张范民，李清河.1998.利用人工神经网络理论对地震信号及地震震相进行识别［J］.西北地震学报，20(4):43－49.

周彦文，刘希强，李铂，等.2010.基于单台P波记录的快速自动地震定位方法研究［J］.地震研究，33(2):183－188.

朱元清，佟玉霞，于海英，等.2002.数字化台网的近震震相自动识别［J］.西北地震学报，24(1):5 －12.

Allen R V.1978.Automatic earthquake recognition and timing from single traces［J］.BSSA，68(5):1521 － 1532.

Allen R V.1982.Automatic phase pickers:their present use and future prospects［J］.BSSA，72(68):225 －242.

Bai C Y，Kennett B L N.2000.Automatic phase-detection and identification by full use of single three-component broadband seismogram［J］.BSSA，90(1):187 －198.

Bear M，Kradolfer U.1987.An automatic phase picker for local and teleseismic events［J］.BSSA，77(4):1437 －1445.

Diehl T，Kissling E，Husen S，et al.2009.Consistent phase picking for regional tomography models:application to the Greater Alpine region［J］.Geophys J Int，176(2):542 － 554.

Haykin S.1996.Adaptive Filter Theory［M］.New Jersey:Upper Saddle River:1－989.

Kuperkoch L，Meier T，Lee J，et al.2010.Automated determination of P-phase arrival times at regional and local distances using higher order statistics［J］.Geophys J Int，181(2):1159 － 1170.

Maeda N.1985.A method for reading and checking phase times in autoprocessing system of seismic wave data［J］.Zisin，38(3):385－379.

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 Interi，113(1 － 4):265 － 275.

Wang J，Teng T L.1995.Artificial neural network-based seismic detector［J］.BSSA，85(1):308 －319.

Zhang H J，Thurber C，Rowe C.2003.Automatic P-wave arrival detection and picking with muhiscale wavelet analysis for single-component recordings［J］.BSSA，93(5):1904 －1912.

Detection of Regional Seismic Events by Kurtosis Method and Automatic Identification of Direct P-wave First Motion by Kurtosis-AIC Method

ZHAO Da-peng1，LIU Xi-qiang2，LI Hong2，ZHOU Yan-wen2
(1.Lanzhou Institute of Seismology，CEA，Lanzhou 73000，Gansu，China)
(2.Earthquake Administration of Shandong Province，Jinan 250014，Shandong，China)

Basing on Kurtosis-AIC and kurtosis methods of P-wave signal，we put forward a new method for real-time detection of regional earthquake event and automatic identification of direct P-wave first motion and apply it to process seismic data recorded by Shandong Seismic Network.The results show as follows:(1)The kurtosis method effectively detect earthquake events，and may effectively reduce false alarm and missing report rates;(2)Compared with phase arrival time results in manual identification，average absolute error of phase arrival time in automatic identification based on Kurtosis-AIC method is(0.09±0.08)s.

kurtosis method;Kurtosis-AIC method;earthquake identification;phase identification

P315.6

A

1000－0666(2012)02－0220－06

2011－02－14.

国家科技支撑计划课题 (2012BAK19B04)和中国地震局地震科技星火计划攻关项目 (XH12029)联合资助.