王子珺，赵伯明

(北京交通大学 土木建筑工程学院，北京　100044)

地震预警系指在破坏性强的地震抵达目标区域前的有限时间(几秒至数十秒)内，利用专用监测预警系统发出预警信息，进而启动防御灾害的紧急处置措施。作为近年来得到快速发展的抗震减灾的重要手段之一，日、美、意等地震多发国家已经在地震预警方面开展了大量的理论与应用研究[1-5]。尤其是对于可以控制其运行的生命线工程和高危行业，例如高铁、核电站、燃气等领域，地震预警系统(EEWS)已在多次地震中被证明能够有效地减轻地震灾害。

高速铁路地震预警主要包括2个重要过程：预警报警和紧急控制。预警报警是在破坏性地震动尚未到达沿线之前，基于预警判断过程发出报警信息；紧急控制是指在接收到报警信息后紧急启动控车措施，使高速运行中的列车实现减速、停车等目的。这种预警—控车机制，可以防止由于地震造成的列车脱轨或倾覆等灾难性事故的发生，减少次生灾害和人员伤亡，保障列车在运营过程中的安全。

根据《高速铁路自然灾害及异物侵限监测系统总体技术方案(暂行)》(铁科技[2013] 35号)，高铁采用基于沿线路布置的原位线型地震监测自备网络，地震监测点设置于牵引变电所、AT所、分区所内，布设间距平均为20 km左右，今后新建高铁也将同期建设此类沿线地震监测网络。线型网络依存于地震波的入射方位，其触发监测点的初始波到时具有较大的分散性，同时原位预警的方式也要求将预警信息进行快速实时处理，因此利用单台站的P波与S波初始段数据进行地震预警是高速铁路所必须采取的基本方式。事实上，日本是最早开展高速铁路地震预警研究的国家，在世界上率先提出了利用单台站记录快速进行地震预警的模式[2]。高铁预警系统是高铁行车安全的相关系统，一旦发生强震漏报误报将会导致高铁重大灾害事件。中国高铁运营速度高，采用以桥代路的里程所占比例较大，所处地震环境复杂，地震时的危险性更高，急需开发具有实时自动、精准快速的适合中国高铁特点的地震预警理论与方法。综上所述，基于单台站的初始先达记录，对地震波进行快速自动高精度地识别，是实现高铁单台站预警目标的首要环节和必要条件。目前，关于P波与S波的识别方法因其应用目的不同而具有不同的特点，但是均不能很好地同时达到自动、快速和准确这3项要求，难以满足中国高铁所需要的单台站预警模式。因此，针对我国高铁地震预警的实际需求，研究提出1套适合于单台站快速准确地震预警的三步骤P波与S波复合自动快速识别理论，并选取记录的2008年汶川Ms8.0级地震的主震及余震的地面运动加速度(以下简称为加速度)，对提出的理论方法进行综合验证。

1　地震数据与数据预处理

中国数字强震动台网在2008年汶川Ms8.0级地震中获得了大量的强震观测记录，观测台网所使用的仪器均为安装于自由地表的数字强震仪，采样频率为200 Hz，测量范围为±2g。除去420个台站获得的主震记录外，截至2008年9月30日12时，台网共记录余震383次；其中Ms4.0级及以上余震611次，Ms5.0级及以上余震56次，Ms6.0级及以上余震8次[6]。由于发震地区位于青藏高原向平原过渡的区域，具有复杂丰富的地质地貌特征，因此广泛分布于该地区的地震记录为开展预警方法的研究以及开发更具有适用性及稳健性的P波和S波识别方法提供了重要的数据保障。

近断层的浅源大震级地震往往造成严重的地震灾害[7]，为此本文选取了汶川地震的主震以及50个余震记录作为研究样本，地震记录的筛选条件为，震级大于Ms5.0且震源距小于150 km。除去其中没有完整得到初始先达P波的记录，最终使用并解析了1 548条满足要求的加速度，相关地震事件的震源位置以及观测台站的空间分布如图1所示。

图1　地震事件的震源位置以及观测台站的空间分布

数据预处理工作如下：首先将原始记录减去记录全长的平均值，完成基线调零；然后对调零后的加速度进行0.5～30 Hz的4阶Butterworth数字带通滤波。该频带适用于本研究区域内的地震记录，对于其他地区，其带通的取值范围则需要根据数据的特性进行调整。

2　P波与S波复合自动快速识别理论

提出的P波与S波复合自动快速识别理论由三步骤构成。首先通过极化分析方法将P波与S波进行分离；然后研究利用改进后的短时窗信号平均值(STA)与长时窗信号平均值(LTA)之比的方法初步确定P波与S波的到时范围；最后根据初步到时范围，通过改良的基于高阶统计量的AIC算法推导P波与S波的精确到时。

2.1　P波与S波的分离

P波与S波的分离主要基于Jurkevics提出的极化分析法[8]，首先通过1个移动的固定时间窗，建立预处理后的三分量加速度的协方差矩阵

(1)

式中：var(*)和cov(*, *)分别表示三分量加速度的自协方差和协方差；z表示竖直分量，n表示南北分量，e表示东西分量。

协方差矩阵C是半正定矩阵，其本质为椭球体方程的二次式系数矩阵。假设矩阵C的3个特征值分别是λ1，λ2和λ3(λ1≥λ2≥λ3)，其对应的特征向量为u=(u1，u2，u3)，表示椭球体的3个轴方向。则根据每一时间窗内的协方差矩阵可以计算得到相应的2组偏振特性参数，基于该偏振特性参数建立本文所使用的震相滤波器。

第1组偏振特性参数为偏振度，定义为

(2)

偏振度Rec的值介于0～1之间，由于纯体波具有完全线性偏振的特点，即偏振度等于1；而背景噪音则不具有相关特点，即偏振度为0；则该组参数可以用于地震波与噪音的分离。

第2组偏振特性参数是表面垂直入射角，可表示为

φinc=arccos|u11|

(3)

式中：u11为特征向量u1的垂直方向余弦。

表面垂直入射角φinc的值介于0°～90°之间，对于区域性地震，由于P波沿着近乎垂直的入射线发生偏振，因此，其入射角的余弦值接近于1；而S波的偏振方向趋于水平，则其入射角的正弦值接近于1。

基于上述这2组偏振特性参数，可以建立2个极化滤波器来分离P波和S波。采用Rosenberger[9]提出的滤波方法，将每一时间步的三分量加速度分别乘以相应的滤波器方程，从而获得极化后的地震波。

选取2008年5月12日14时43分汶川地震余震的理县木卡台站的原始加速度三分量记录，如图2所示，该余震震级为Ms6.3，震中距约为56 km；采用0.2，2.0，5.0 s不同长度时间窗进行P波和S波的滤波，如图3所示。由图3可知：采用较小的时间窗长度(0.2 s)不能很好地将地震信号与噪音信号分离；而随着时间窗长度的不断增大，滤波方程趋于平滑且对信号的变化更加灵敏，从而易于对不同类型的地震波进行判断识别；这是因为移动时间窗的长度直接关系到偏振度和入射角的计算，进而会对滤波器的滤波效果产生影响。根据Baillard等[10]的研究成果，在进行偏振分析时，时间窗长度应大于地震波最短周期的3倍，同时要小于P波与S波之间的时间差；通过分析发现，采用1个长度为2.0 s的移动时间窗即可达到理想的效果。通过极化滤波后的南北、东西、竖直三分量加速度如图4所示。

图2　原始加速度三分量

图3　不同长度移动时窗下P波与S波的滤波结果

图4　三分量加速度的P波与S波分离结果

由图4可知：利用建立的极化滤波器能够很好地完成P波与S波的分离，从而有效地保证震相的初至识别。

2.2　识别算法

首先对传统的STA/LTA方法进行优化改进研究，根据得到的极化后的P波与S波，利用该改进方法初步确定P波与S波的到时范围，在该范围内，利用基于高阶统计量的AIC算法分别推导P波与S波的精确到时，具体过程如下。

STA/LTA算法的基本原理是利用STA与LTA比值的变化来反应信号的变化。由于短时间窗反映的主要是信号幅值的瞬时变化，而长时间窗反映的是待检信号背景噪音的平均幅值；当震相到达时，STA会比LTA的变化频度高，相应的比值会有一个明显的增加，当比值超过某一触发阈值时，则判定有震相到达。为了更加明显地反映信号到达时频率及幅值等特征量的变化，通常采用Allen[11]提出的特征函数CF(characteristic function)代替原始信号进行STA/LTA的计算。该特征函数为

CF(i)=x(i)2+W(x(i)-x(i-1))2

(4)

式中：x(i)(i=1，2，…，L)为ti时刻的加速度，L为随机变量x的样本数；W为随着采样率和背景噪音性质而变化的加权常数，根据经验取值为3。

该特征函数CF考虑了实际地震的加速度以及向前一阶差分的影响，能够灵敏地反映信号到时的频率及幅值变化特征。虽然该特征函数能够有效降低震相识别的误判率，但是由于触发阈值是根据经验设定的固定值，在低信噪比的情况下可能会造成震相的漏判或误判。所以为了提高P波与S波的初始到时的识别率，提出导入1个调整系数δ，对传统的STA与LTA的比值进行优化改进，该调整系数的定义为

(5)

其中，

(6)

由于标准差可以衡量某一数据序列的数值离散程度，在地震事件发生前，系数δ趋近于1；当P波或S波的信号到达时，系数δ会有一个较大幅度的增长。因此这个改进后的STA/LTA方法可以明显地提高对信噪比低或初动信号不明显的记录的识别率，从而大幅度提高识别的精度和速度。

传统STA/LTA方法与改进后STA/LTA方法的识别效果如图5所示，可以看出，改进后的方法对震相初至识别有非常明显的提升效果。利用该方法可以确定震相的到时范围，其精确的到时判定将会在下一阶段完成。P波的拾取使用竖直向加速度分量，S波的拾取使用水平向加速度分量，并设初始到时的时间为T初。

图5传统STA/LTA方法与改进后STA/LTA方法的识别效果比较

移动时间窗的长度由采样频率和特征信号的周期所决定。对于STA，若时间窗的长度选择过大，则会导致特征函数过于平滑而无法检测到信号的瞬时变化，从而产生漏触发；若时间窗的长度选择过小，则会对短周期的干扰信号过于敏感从而造成误触发。对于LTA，其窗长取值应能反映背景噪音的平均水平。基于对大量记录的试算，本文取：短时窗tSTA=0.5 s，长时窗tLTA=3～5 s，用户定义参数Δ=35；P波与S波拾取的触发阈值分别设为5和10。

为了精确确定震相的到达时刻，在下一阶段中采用基于改良的高阶统计量的AIC函数对地震波进一步的识别。根据Maeda[12]，计算地震数据x(i)的AIC函数，其最小值对应的时间即为地震震相的精确到时，AIC函数为

AIC(k)=klog{var(x[1，k])}+(L-K-

1)log{var(x[k+1，L])}

(7)

式中：k为滑动变量，k∈[1，L-1]。

基于数学分析理论，高阶统计量与二阶统计量相比具有适于检测和描述系统的非线性的特点。在低信噪比的条件下，可有效抑制背景噪音的影响，易于对较弱的地震信号进行检测和识别，因此本文导入峰度函数代替二阶方差函数，并应用于AIC函数以确定震相的精确到时。峰度函数K为样本的四阶中心矩，即

(8)

在震相初始到时T初的基础上，通过确定1个时间窗[T初-0.5，T初+0.5]，应用峰度-AIC函数对P波与S波的到时进行精确识别。选择0.5 s既可以保证将触发时刻控制在该时间窗内，又能够满足地震预警所需的实时性和时效性要求。

图6表示应用峰度-AIC函数对震相进行二次拾取的结果。由图6可以看出，利用第1阶段改进后的STA/LTA方法初步确定P波与S波的到时后，再求取峰度-AIC函数的最小值，就能够准确地判断震相的实际到时。

图6　峰度-AIC函数的最小值对应震相的精确到时

由此可见，由于改进后的SLA/LTA方法具有适应性强、拾取率高的特点，在使用该方法进行震相初始识别的基础上，通过改良的基于高阶统计量的AIC函数进一步保证了识别的稳定性和精确性。因此，本文命名上述由3个步骤构成的识别方法为P波与S波复合自动快速识别方法。

3　汶川地震数据的应用与分析

将提出的P波与S波复合自动快速识别方法应用于记录的2008年5月12日14时43分汶川Ms8.0级地震的主震以及其后50个余震的1548条地震的加速度，识别震相的精确到时，并与人机交互处理结果和改进后的STA/LTA方法处理结果进行对比，验证提出的复合方法的准确性和有效性。在判别过程中，P波的拾取针对竖直方向的加速度分量，S波的拾取针对南北以及东西2个水平方向的加速度分量。以2.1节使用的余震为例，选取茂新南新、理县木卡和理县沙坝3个台站所获得的地震加速度。其P波的到时识别过程如图7所示，S波的到时识别过程如图8和图9所示。

图7不同台站的P波的到时识别过程(采用竖直方向加速度分量)

图8不同台站S波的到时识别过程(采用南北方向水平加速度分量)

为进一步研究P波与S波复合自动快速识别方法的准确性，以人机交互结果为标准，对使用的1 548个三分量加速度的自动识别结果进行误差计算，其偏差分布如图10所示。P波到时的自动识别的结果为：拾取正确率达到91%，其中89%的误差在0.1 s以内，平均偏差为-0.021 s，标准差为0.068 s。S波到时的自动识别结果为：按照识别误差小于0.5 s的精度要求，识别有效率高达85%，其中有92%的误差在0.2 s以内，平均偏差为-0.025 s，标准差为0.169 s。

图9不同台站S波的到时识别过程(采用东西方向水平加速度分量)

图10P波与S波的自动识别与人机交互识别结果的偏差分布

P波与S波自动识别的平均偏差基本一致，但P波识别的标准差以及识别率优于S波，因为P波属于纵向压缩波，其传播过程比横向S波相对单纯，而S波的产生及传播过程比较复杂，受传播路径的影响较大，对其识别相对有较大难度。另外，本文对P波的拾取只针对竖直方向的加速度，而S波的拾取则使用了水平2个方向的记录，S波的样本数量是P波的2倍，因此也会造成S波识别率要低于P波的结果。考虑到每个观测点都会受到复杂

的震源、传播路径以及场地效应的影响，根据既有研究无法实现百分百的识别率。

此外本文也对改进后的STA/LTA方法与人机交互的结果进行了比较：对于P波到时的识别，两者之间偏差的平均值为-0.217±0.326 s；对于S波到时的识别，两者偏差的平均值为-0.191±0.379 s。该结果表明，通过对传统的STA/LTA方法进行改进，大幅度提高了识别精度和速度，但由三步骤构成的P波与S波复合自动快速识别方法的效果更优。

因此，综合以上结果，本文提出的结合改进后的STA/LTA方法以及改良的基于高阶统计量的AIC函数得到的复合方法可以显著提高P波和S波的识别效果。

4　结　论

(1)高速铁路地震预警的研究就是保障预警过程的实时性、准确性和高效性。针对我国高铁的特点，以及高铁地震预警对时间和精度的要求，基于沿线路排列的原位地震监测网络所必需的单台站预警模式，本文提出了1套适合于高铁单台站快速准确地震预警的三步骤P波与S波复合自动快速识别理论。

(2)为了校验提案方法的泛用性，整理解析了具有代表性的地震——2008年汶川Ms8.0级地震的主震及余震记录1 548条，通过将自动识别方法与人机交互处理的结果进行对比，验证了该方法可以准确、有效地应用于具有复杂的震源、传播路径以及场地效应的地震进行快速预警。比之其他既有方法具有实时性、精度高、速度快的特点，能够满足高铁单台站P波或S波地震预警的严格要求。由于强震观测点条件的多样性，对P波与S波的识别率造成了一定影响，如果使用针对高铁地震预警建设的台站所获得的地震记录，那么成功率会进一步提高。

[1]NAKAMURA Y，SAITAJ，SATO T. On an Earthquake Early Warning System (EEW) and Its Applications [J]. Soil Dynamics and Earthquake Engineering，2011，31(2)：127-136.

[2]NAKAMURA Y. On the Urgent Earthquake Detection and Alarm System (UrEDAS)[C]//Proceedings of 9th World Conference on Earthquake Engineering. Tokyo-Kyoto，Japan：International Conference for Earthquake Engineering (IAEE)& the 9th WCEE Organizing Committee，1988，Ⅶ：673-678.

[3]BROWN HM，ALLEN RM，HELLWEG M，et al. Development of the ElarmS Methodology for Earthquake Early Warning: Realtime Application in California and Offline Testing in Japan[J]. Soil Dynamics and Earthquake Engineering，2011，31(2)：188-200.

[4]SATRIANO C，ELIA L，MARTINO C，et al. PRESTo, the Earthquake Early Warning System for Southern Italy: Concepts, Capabilities and Future Perspectives [J]. Soil Dynamics and Earthquake Engineering，2011，31(2)：137-153.

[5]IANNACCONE G，ZOLLO A，ELIA，et al. A Prototype System for Earthquake Early-Warning and Alert Management in Southern Italy [J]. Bulletin of Earthquake Engineering，2010，8(5)：1105-1129.

[6]国家强震动台网中心.2008年汶川地震数据集[DB/OL].[2015-11-01].ftp://www.csmnc.net.

[7]WU Y M，KANAMORI H. Experiment on an Onsite Early Warning Method for the Taiwan Early Warning System[J]. Bulletin of the Seismological Society of America，2005，95(1)：347-353.

[8]JURKEVICS A. Polarisation Analysis of Three-Component Array Data [J]. Bulletin of the Seismological Society of America，1988，78(5)：1725-1743.

[9]ROSENBERGER A. Real-Time Ground Motion Analysis: Distinguishing P and S Arrivals in a Noisy Environment[J]. Bulletin of the Seismological Society of America，2010，100(3)：1252-1262.

[10]BAILLARD C，CRAWFORD WC，BALLU V，et al. An Automatic Kurtosis-Based P- and S-Phase Picker Designed for Local Seismic Networks [J]. Bulletin of the Seismological Society of America，2014，104(1)：394-409.

[11]ALLEN R V. Automatic Phase Pickers: Their Present Use and Future Prospects [J]. Bulletin of the Seismological Society of America，1982，72(6)：S225-S242.

[12]MAEDA N. A Method for Reading and Checking Phase Times in Auto Processing System of Seismic Wave Data [J]. Zisin=Jishin，1985，38：365-379.