2013, Vol. 56
2. 福建省地震局, 福州 350003
2. Earthquake Administration of Fujian Province, Fuzhou 350003, China
地震预警技术是近20年来发展起来的一种新的减轻地震损失、降低地震次生灾害、减少人员伤亡的有效手段[1-2]，其基本原理是：携带地震信息的P波波速（6~7 km/s)大于带来较大破坏的S波波速(约为3.5 km/s)和面波波速，且地震波速远远小于电磁波速度，地震预警系统可在破坏性地震发生后，基于近震源地震监测仪器得到的记录数据，快速侦测地震，以空间传播换取时间，在破坏性地震波到达特定目标区前，发布地震警报.
近年来，随着地震观测仪器、通信技术、计算机技术和自动处理能力的迅速发展，国际上已经建设了多个针对特定设施、单个城市甚至全国范围的地震预警系统[1-3]，很多已经取得了减灾实效，如日本铁路UrEDAS系统[4]、日本紧急地震速报系统[3, 5-6]、墨西哥SAS和SASO系统[7-8]、我国台湾地区试验性预警系统[9-10]、土耳其伊斯坦布尔市预警系统[11]及罗马尼亚布加勒斯特地震预警系统[12]等.此外美国加利福尼亚州及意大利等一些国家也在建设类似系统[1].我国也在研发基于实时传输测震台网和强震动台网的地震预警系统并拟先期在福建地区、首都圈地区和兰州地区进行布设.
地震预警是在真实的地震发生后，基于仪器记录，尽可能地利用先期获得的信息，较精确地估计地震发生的地点和大小，在破坏性地震波到达前预测特定目标区地震动场、破坏性地震波到达时间，并预测目标区可能的地震破坏情况，并决定是否发布预警信息.从技术上来说，地震事件的判别、震相的自动拾取是实现基于地震台网地震预警的首要条件，在传统地震学中，震相一般在地震图上人工拾取，震相的自动识别开始于对海量波形数据自动处理，随着实时地震学的发展，对震相的自动精确拾取越来越受到地震学家的重视.对地震预警重点考虑的地方震来说，目前主要考虑的震相为容易自动拾取的P波震相及S波震相，在本文中我们重点研究了可用于地震预警的P波震相到时自动拾取技术，关于S波的自动拾取问题将另文讨论.
目前绝大多数震相自动拾取方法是提取信号和噪声的不同特征来作为震相到来的判据，比如用幅值变化、频率组成变化、波形的相似性和动力学特征变化等来判断震相.归纳起来目前常用的方法有：能量变化分析、偏振分析和自回归方法（AR)等.能量方法中最常用的为长短时间平均（STA/LTA)方法[13-16]，长短时平均反映了幅值的瞬时变化.偏振分析方法是利用震相到来时质点运动的偏振方向会发生改变的特征来拾取和判别震相[17-19].自回归方法是基于可以把地震波划分为局部平稳段的假定，应用AIC或者AR-AIC方法来进行不同震相的特征拾取[20-24].
对于地震预警，其P波震相拾取的特殊性主要表现在：只能用到地震台网中有限台站的有限时间段信息，其方法能进行实时自动处理，且能剔除非地震事件干扰.本文依据地震预警技术的特点，发展了一套基于长短时平均（STA/LTA)、池赤准则（AIC)和Delaunay三角剖分的P波震相多步骤自动拾取技术，其可靠性和时效性能满足地震预警需求，方法可对地震台网多种实时传输记录（如测震速度记录和强震动加速度记录)进行实时处理.
2 基于STA/LTA和AIC算法的P波震相自动拾取技术 2.1 宽频带记录的实时仿真处理对于宽频带速度型地震仪器，其频带宽，记录到的地脉动噪声水平亦高，对小震及距离较远震来说，由于信噪比较低，实际地震信号往往湮没在噪声中而造成震相自动识别的困难，并影响拾取精度，为了更好地突出近震及地方震事件信号，并尽量压制长周期噪声，我们把宽频带仪器速度记录统一仿真为自振周期1 s，阻尼比0.707的短周期速度记录[25-26].图 1为2006年9月17日发生在福建省华安的ML2. 5级地震的尤溪台原始宽频带记录波形及仿真成短周期后时程，从图中可以看出，仿真后记录突出了近震及地方震的高频信息，压制了低频脉动噪声.
|
图 1 典型宽频带垂直分向速度记录(a)及仿真为周期为1s的速度时程(b) Fig. 1 Typical broadband velocity record of vertical component (a) and the simulation record with natural period 1 second (b) |
对P波震相，一般在垂直向幅值较大，对S波来说，一般水平向幅值较大, 对于要拾取的P波，我们希望其特征尽量放大.我们应用Allen[13-14]提出的特征函数来放大垂直向的记录特征.
P波拾取特征函数：
|
(1) |
式中，xud(k）为k时刻垂直向记录值.在以上特征函数中，考虑到了实际速度记录及向前差分(与加速度项有关）的影响，突出了近震及地方震记录中的高频信息.图 2表示了应用式(1)后的P波特征函数.
|
图 2 垂直向记录(a)，P波特征函数(b)及STA/LTA值(c) Fig. 2 Broadband velocity record (a), characteristic function of P-wave (b) and STA/LTA time history (c) |
STA/LTA方法是一种能量方法，广泛应用于信号检测中[13-16]，特别是对于地震弱信号.目前已经发展了很多应用不同时间窗的拾取方法，本文应用公式如下：
|
(2) |
式中CF为特征函数，i为当前时刻点，k1和k2为当前时刻i前某时刻点且k2 < k1 < i.
在P波到来后，记录的幅值会有较大的变化，在短窗内，其平均值STA变化快，刻画的主要是信号幅值的瞬时变化，其窗长一般取长于待测特征信号的几个周期左右，太短则对短周期的干扰更敏感，容易产生误触发，如果太长则显示不出待测信号的瞬时特征，容易产生漏触发；在长窗内其平均值LTA变化稍缓，刻画的是相对于待检信号的背景噪声平均大小，其窗长取值应该能反映背景噪声水平. STA/LTA进一步刻画了记录幅值的瞬时变化，在求STA/LTA后，其特征变得更为明显，由于记录的平均，其时程曲线变的更为平滑.设定一定的触发阈值，如果STA/LTA值超过设定阈值，则认为P波触发，对于阈值的选取，如果阈值过大，则有可能会出现漏触发，如果太小，则对很多微弱干扰会产生误触发.经过对大量记录的试验，对P波拾取中所用短时间及长时间窗长及阈值参数如表 1所示.图 2为对一典型记录(2006年9月17日在福建省华安ML2.5级地震龙岩台记录）求取特征函数并计算STA/LTA的结果，由图中可以看出，在求取特征函数后，其P波特征被放大，在求取STA/LTA后其P波段上升明显.STA/LTA方法的突出优点为适应性强、拾取效率高、稳定可靠，其缺点一是其触发点一般滞后于实际P波到时点，二是如果记录中存在干扰信号，对于设定阈值同样会产生地震事件触发.
|
表 1 STA/LTA窗长及阈值参数选取 Table 1 Time windows and threshold-trigger value |
对于STA/LTA阈值P波触发方法，只能粗略识别P震相的到时，且识别出的到时往往滞后于实际的到时，也就是说在达到阈值时，得到的只是P波到达的大致位置.20世纪70年代日本学者赤池弘次（Akaike)提出一个基本信息量的定阶准则，AIC准则[27]，对地震记录如假设震相到时前后的地震记录是两个不同的稳态过程，可应用自回归的AR-AIC方法来进行地震震相的判别.区别于AR-AIC方法，Maeda (1985)建议可由地震波形数据直接计算AIC函数，而不需要求出AR系数[28-30]，对地震记录x(i)(i=1，2, …，L)，AIC检测器定义为：
|
(3) |
其中，k的范围为地震图某窗长内所有的采样点，var表示方差.震相到时对应于窗长内AIC函数的最小值.
不同于在整条地震图上直接应用AIC准则，我们在用STA/LTA粗略拾取到P波后，在固定窗内用记录的垂直向特征函数来进行第二步P波震相精拾取.具体方法为：在STA/LTA触发点前推和后推一定时间窗，在窗内应用AIC准则，窗内AIC最小值认为是到时点.考虑到地震预警的实时性和时效性要求，在本文中对P波拾取的前推点数为100个采样点，后推点数为10个采样点，即对采样频率为50 Hz的地震记录，前推时间为2 s，后推时间为0.2 s.图 3为应用固定窗AIC准则对P波进行进一步拾取的情况（所用台站记录数据同图 2)，由图中可以看出，在STA/LTA粗略得到到时点后，应用AIC最小值能很好地判断出P波震相实际到时.
|
图 3 P波AIC捡拾结果与STA/LTA及人工捡拾结果对比 Fig. 3 P-arrival picks of AIC, STA/LTA and human |
地震记录上会记到非天然震动信号，对地震事件和P波来说相当于干扰，干扰的存在会降低地震信号的自动拾取的精度，还会造成天然地震事件的误拾取.对于非震动干扰，其各台站到时在顺序上是不相关的，即不符合地震震相的传播规律，此规律可作为是否为震动信号的逻辑判据.
Delaunay (以下简称D)三角剖分是计算几何学的重要基础，应用于许多学科领域中，其性质和定义可参考相关文献[31]，图 4为福建测震台网的D三角剖分，可以看出D三角剖分在平面上反映了各台站的相邻关系，可作为台站到时顺序的逻辑判别依据.假定地震P波传播速度恒定且地球介质均匀，台网中对同一次地震第i个接收到P波的台站，称为第i台，记为S(i).由D三角剖分性质和台站的几何分布，可以对误触发干扰信号进行剔除.具体操作如下，如果对于台网数据我们得到某台站S(i)的P波到时，则下一个P到时台站S(i+ 1)-定位于以S(i)台站为顶点的D三角的另外顶点.应用以上结论我们对判断出的P波到时错误信号进行排除，以前三台为例，方法如下：第1台触发后，如果触发的第2台不在第1台周围D三角顶点上，则等待第3台触发; 如果第3台在第1台周围D三角顶点上则认为第2台的第1个触发错误，进而判断第2个触发到时; 如果第3台在第2台周围D三角顶点上则认为第1台第1个到时错误，进而判断第2个触发；如果三台D三角顶点互不相关，则把第3台当成第1台并循环向下判断.如果某台发生多个不相关触发，则把此台直接剔除.图 5表示了三角剖分与干扰信号的判别流程图.应用以上D三角剖分逻辑判断虽然可能出现误相关情况，但我们对大量台网记录验证后发现，应用此方法后能剔除大部分不相关干扰触发情况，大大提高了拾取的可靠性.
|
图 4 福建数字地震台网台站Delaunay三角剖分图 Fig. 4 Delaunay triangulation of Fujian digital seismic net work |
|
图 5 Delaunay三角剖分与干扰信号的判别流程图(以前三台为例) Fig. 5 Discriminant flow chart of interference signal using Delaunay triangulation (Take the first three stations as an example) |
应用以上介绍的P波拾取方法，在地震预警中我们采取以下步骤对实时传输地震台网中每个台站的地震记录进行P波自动拾取，流程如图 6所示.
|
图 6 P波拾取流程图 Fig. 6 Flow chart of automatic P-arrival picking |
(1) 对实时传输的测震速度记录或加速度记录进行实时仿真，统一仿真为自振周期为1 s，阻尼比为0.707的短周期速度；
(2) 应用公式(1)实时计算仿真后记录的特征函数；
(3) 利用公式(2)等步长移动长窗和短窗（根据计算速度需要可取1个或者若干个时间间隔）并实时计算STA/LTA数值；
(4) 如超过触发阈值后，则在固定窗长内计算AIC值；
(5) 等待下一个台站触发并用Delaunay三角剖分判别是否是干扰信号.
5 计算实例 5.1 所用数据我们应用了福建数字地震台网自2000年2月至2007年3月，记录到的所有ML≥2. 5网内地震82个，最大震级ML=4.9, 2005年7月至2007年5月所有网外100 km范围内（网缘至网缘以外100 km) ML≥3.0地震67个，最大震级ML=4.8.地震台站及震中分布如图 7和图 8所示.
|
图 7 福建省地震观测台站分布图(截至2008年) Fig. 7 Distribution of Fujian digital seismic network (2008) |
|
图 8 地震台站及震中分布图 Fig. 8 Distribution of the stations and epicenters |
对网内82个地震164套三分向记录，由于选取的震级下限较低（P波信噪比太低或者波形不清），有3个台站人工亦不能分辨出P波到时，对其余161套记录，自动拾取有1套结果误差1. 5 s以上(近台小震干扰），其余160套记录均得出较好结果，P波拾取正确率达到99. 4%，自动拾取平均偏差0.00 s，标准差0.05 s.
对网外100 km范围内67个地震，其中3个地震前10 s左右叠加小震，无法自动处理出结果，对于震前叠加小震干扰这种情况不在本文讨论的范围内.对其余64个地震128套三分向记录，对P波到时自动拾取结果如下：拾取正确率达到100%，自动拾取平均偏差-0.01s，标准差0.13 s.
相对于0.02 s的数据采样间隔，其结果是可靠的.拾取结果如图 9, 图 10和表 2所示.
|
图 9 网内地震P波自动拾取偏差 Fig. 9 P-arrival picking results of the earthquakes inside the network |
|
图 10 网外100km地震P波自动拾取偏差 Fig. 10 P-arrival picking results of the earthquakes 100 km outside the network |
|
|
表 2 网内地震及网外100 km地震拾取结果 Table 2 P-arrival picking results of inside and 100 km outside seismic network |
针对地震预警技术中对p波震相拾取的要求，综合应用了实时地震记录仿真，P波特征函数计算、STA/LTA方法和AIC方法，发展了一种多步骤的拾取方法，提出了一种基于有限先到台的Delaunay三角剖分的非震动干扰信号剔除方法，并应用福建省地震数字台网记录对方法进行了验证，目前方法已经用到了福建省地震预警试验系统中，通过对结果的分析可得到如下结论：
(1) 所用方法均为递归算法，可对地震台网实时数据流进行实时计算，适合地震预警对数据的实时自动处理，对单台其理论拾取时间延迟小于应用STA/LTA的短窗时间窗长，试验表明计算需时可忽略不计；
(2) SLA/LTA方法具有适应性强、拾取效率高、稳定性好的特点，适合于地震弱信号判断，但其也会带误拾取和难以精确捡拾到时点的缺点，AIC方法具有在粗略到时点附近进行到时精确捡拾的特点.本文方法综合了SLA/LTA方法和AIC方法的特点，根据地震预警要求对方法进行了细节处理，采取多步骤对实时传输台网记录进行处理，提高了拾取的稳定性和精确性；
(3) 针对地震预警只能用到有限台站记录的特点，提出了一种基于有限先到台的Delaunay三角剖分干扰信号剔除方法，大大降低了误触发情况，提高了拾取的可靠性；
(4) 应用福建省地震台网近震源台站分析结果表明，其自动拾取率在99%以上，拾取标准差小于0.2 s，方法可用于地震预警：
(5) 本文选取的震例震级下限较低（ML2. 5~4.9)，其记录信噪比也相对较低，地震预警主要针对的是有破坏性大震，可以预见，对于大震其P波拾取效果会更好.
需要指出的是，地震预警具有自动处理和实时处理的特点，P波震相拾取的可靠性基于对地震事件的判别，实际地震记录上经常记到非天然地震干扰信号，具体可以分为环境干扰及仪器干扰.环境干扰包括气旋和海浪、周边人为活动、矿井塌陷及爆破等，仪器自身干扰包括仪器噪声、尖峰信号、直流分量及仪器脉检等.干扰的存在会降低初动地震信号的自动拾取的精度，还会造成天然地震事件的误拾取.对于干扰信号的其它判别方法，我们将另文讨论.
致谢感谢福建省地震局监测中心为本文的撰写提供了测震记录数据及编目数据.
| [1] | Satriano C, Wu Y M, Zollo A, et al. Earthquake early warning: Concepts, methods and physical grounds. Soil Dynamics and Earthquake Engineering , 2011, 31(2): 106-118. DOI:10.1016/j.soildyn.2010.07.007 |
| [2] | 马强.地震预警技术研究及应用[博士论文].哈尔滨:中国地震局工程力学研究所, 2008. Ma Q. Study and application on earthquake early warning[Ph. D. thesis](in Chinese). Harbin: Institute of Engineering Mechanics, China Earthquake Administration, 2008. http://cdmd.cnki.com.cn/Article/CDMD-85406-2009057325.htm |
| [3] | Allen R M, Gasparini P, Kamigaichi O, et al. The status of earthquake early warning around the world: an introductory overview. Seism. Res. Lett. , 2009, 80(5): 682-693. DOI:10.1785/gssrl.80.5.682 |
| [4] | Nakamura Y. UrEDAS, urgent earthquake detection and alarm system, now and future. 13th World Conference on Earthquake Engineering, Vancouver, 2004. |
| [5] | Kamigaichi O. JMA earthquake early warning. J. Japan. Assoc. Earthq. Eng. , 2004, 4(3): 134-137. |
| [6] | Wu Y M, Kanamori H. Development of an earthquake early warning system using real-time strong motion signals. Sensor , 2008, 8(1): 1-9. DOI:10.3390/s8010001 |
| [7] | Espinosa-Aranda J M, Jiménez A, Ibarrola G, et al. Mexico City seismic alert system. Seism. Res. Lett. , 1995, 66(6): 42-53. DOI:10.1785/gssrl.66.6.42 |
| [8] | Suarez G, Novelo D, Mansilla E. Performance evaluation of the seismic alert system (SAS) in Mexico City: a seismological and a social perspective. Seism. Res. Lett. , 2009, 80(5): 707-716. DOI:10.1785/gssrl.80.5.707 |
| [9] | Wu Y M, Teng T L. A virtual subnetwork approach to earthquake early warning. Bull. Seism. Soc. Am. , 2002, 92(5): 2008-2018. DOI:10.1785/0120010217 |
| [10] | Hsiao N C, Wu Y M, Shin T C, et al. Development of earthquake early warning system in Taiwan. Geophys. Res. Lett. , 2009, 36(5): L00B02. DOI:10.1029/2008GL036596 |
| [11] | Bose M, Wenzel F, Erdik M. PreSEIS: a neural network-based approach to earthquake early warning for finite faults. Bull. Seism. Soc. Am. , 2008, 98(1): 366-382. DOI:10.1785/0120070002 |
| [12] | B?se M, Ionescu C, Wenzel F. Earthquake early warning for Bucharest, Romania: novel and revised scaling relations. Geophys. Res. Lett. , 2007, 34(7): L07302. DOI:10.1029/2007GL029396 |
| [13] | Allen R E. Automatic earthquake recognition and timing from single traces. Bull. Seism. Soc. Am. , 1978, 68(5): 1521-1532. |
| [14] | Allen R E. Automatic phase pickers: Their present use and future prospects. Bull. Seism. Soc. Am. , 1982, 72(6B): 225-242. |
| [15] | Baer M, Kradolfer U. An automatic phase picker for local and teleseismic events. Bull. Seism. Soc. Am. , 1987, 77(4): 1437-1445. |
| [16] | Earle P S, Shearer P M. Characterization of global seismograms using an automatic-picking algorithm. Bull. Seism. Soc. Am. , 1994, 84(2): 366-376. |
| [17] | Cichowicz A R. An automatic S-phase picker. Bull. Seism. Soc. Am. , 1993, 83(1): 180-189. |
| [18] | Earle P S. Polarization of the Earth's teleseismic wavefield. Geophys. J. Int. , 1999, 139(1): 1-8. DOI:10.1046/j.1365-246X.1999.00908.x |
| [19] | Vidale J E. Complex polarization analysis of particle motion. Bull. Seism. Soc. Am. , 1986, 76(5): 1393-1405. |
| [20] | Sleeman R, Van Eck T. Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings. Phys. Earth Planet. Interiors , 1999, 113(1-4): 265-275. DOI:10.1016/S0031-9201(99)00007-2 |
| [21] | Leonard M, Kennett M B L N. Multi-component autoregressive techniques for the analysis of seismograms. Phys. Earth Planet. Interiors , 1999, 113(1-4): 247-264. DOI:10.1016/S0031-9201(99)00054-0 |
| [22] | Leonard M. Comparison of manual and automatic onset time picking. Bull. Seism. Soc. Am. , 1999, 90(6): 1384-1390. |
| [23] | Takanami T, Kitagawa G. A new efficient procedure for the estimation of onset times of seismic waves. J. Phys. Earth. , 1988, 36(6): 267-290. DOI:10.4294/jpe1952.36.267 |
| [24] | Takanami T, Kitagawa G. Multivariate time-series model to estimate the arrival times of S-waves. Computers and Geosciences , 1993, 19(2): 295-301. DOI:10.1016/0098-3004(93)90127-Q |
| [25] | 金星, 马强, 李山有, 等. 宽频带强震仪与地震仪同一台基上记录仿真对比研究. 地震工程与工程振动 , 2004, 24(5): 7–12. Jin X, Ma Q, Li S Y, et al. Comparison research on the records of wide band strong motion seismograph and seismometer at same station. Earthquake Engineering and Engineering Vibration (in Chinese) , 2004, 24(5): 7-12. |
| [26] | Jin X, Ma Q, Li S Y. Real-time simulation of ground displacement by digital accelerograph record. Acta Seismologica Sinica , 2005, 18(1): 82-88. DOI:10.1007/s11589-005-0009-9 |
| [27] | Akaike H. Information theory and an extension of the maximum likelihood principle.//Petrov B, Csaki F eds. 2nd International Symposium on Information Theory. Budapest: Akademiai Kiado, 1973: 267-281. |
| [28] | Maeda N. A method for reading and checking phase times in auto-processing system of seismic wave data. Zisin (J. Seismol. Soc. Jpn.) , 1985, 38(3): 365-380. |
| [29] | Zhang H J, Thurber C, Rowe C. Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single-component recordings. Bull. Seism. Soc. Am. , 2003, 93(5): 1904-1912. DOI:10.1785/0120020241 |
| [30] | 王继, 陈九辉, 刘启元, 等. 流动地震台阵观测初至震相的自动检测. 地震学报 , 2006, 28(1): 42–51. Wang J, Chen J H, Liu Q Y, et al. Automatic onset phase picking for portable seismic array observation. Acta Seismologica Sinica (in Chinese) , 2006, 28(1): 42-51. |
| [31] | O'Rourke J. Computational Geometry in C (2nd ed). London: University Press, 1998 : 181 -226. |