地震学报  2012, Vol. 34 Issue (5): 593-610
引用本文
金星, 张红才, 李军, 韦永祥, 马强. 2012. 地震预警震级确定方法研究. 地震学报, 34 (5): 593-610.
Jin Xing, Zhang Hongcai, Li Jun, Wei Yongxiang, Ma Qiang. 2012. Research on earthquake early warning magnitude estimate. Acta Seismologica Sinica , 34 (5): 593-610.
地震预警震级确定方法研究
金星1,2, 张红才1,2 , 李军1, 韦永祥1, 马强2    
1. 中国福州350003福建省地震局;
2. 中国哈尔滨150008中国地震局工程力学研究所
收稿日期 2011-08-03收到初稿， 2011-12-19决定采用修改稿.
基金项目： 国家科技支撑计划(2009BAK55B02， 2009BAK55B03)和福建省地震局青年基金课题(Y201002) 共同资助.
通讯联系人,E-mail:zhanghongcai521@163.com
摘要：地震预警技术是减轻地震灾害损失的有效手段之一. 地震预警系统中， 地震震级计算是最重要也是最困难的部分之一. 利用日本KiK-net台网和四川汶川余震共142次地震事件的记录， 分别采用τc和Pd方法统计得到了地震预警震级的计算公式， 震级计算的方差分别为0.62和0.56个震级单位. 为消除震级计算过程中出现的震级饱和现象， 作者拓展了Pd方法， 提出了一套对位移幅值连续追踪测定的算法. 当时间窗长度为10 s时， 采用该方法的震级计算方差仅为0.37个震级单位， 充分满足地震预警系统的需求. 同时， 该方法也实现了信息的连续过渡， 提高了对现有信息的利用率. 最后， 还对位移幅值Pd用于地震动峰值PGV的估计以及不同特征参数间的相容性等内容进行了讨论.
关键词： 地震预警     地震震级     τc方法     Pd方法    
Research on earthquake early warning magnitude estimate
Jin Xing1,2, Zhang Hongcai1,2 , Li Jun1, Wei Yongxiang1, Ma Qiang2    
1. Earthquake Administration of Fujian Pronvince, Fuzhou 350003, China;
2. Institute of Engineering Mechanics, China Earthquake Administration, Harbin 150008, China
Abstract: Earthquake early warning (EEW) is one of the effective ways to reduce earthquake damage. Earthquake magnitude estimation is the most important and also most difficult part in the whole EEW system. Based on 142 earthquake records from KiK-net in Japan and aftershock records of the great Wenchuan earthquake in Sichuan, China, we worked out the magnitude estimation relationships by using τc and Pd methods, with corresponding variances of 0.62 and 0.56, respectively. In order to remove the effect of magnitude saturation, we expanded the Pd method and proposed a new method of continuously calculating the peak displacement. When the time window is 10 s, the variance of magnitude estimation is reduced to 0.37, fully meeting the need of EEW systems. Meanwhile, our method realized the continuous transition of different information, and improved the utilization of seismic records. We also discussed the relation between Pd and PGV, and the compatibility between different characteristic parameters.
Key words: earthquake early warning     earthquake magnitude     τc method     Pd method    

引言

我国作为世界上地震活动最强烈的国家之一遭受着严重的地震灾害． 20世纪内1/3的大陆地震都发生在我国， 这些地震造成的死亡人数占全球因地震死亡人数的50%以上． 与此同时， 我国有超过一半以上国土都位于高烈度区域， 其中就包括全国23个省会城市及2/3的百万人口以上的大城市． 一旦破坏性地震在这些城市周围发生， 必将造成难以估量的财产损失和人员伤亡． 我国所面临的防灾减灾形势依然十分严峻． 1976年唐山7.8级地震造成的死亡人数超过24万， 全城倾毁； 2008年5·12汶川8.0级特大地震造成近7万人死亡2万人失踪； 2010年4月14日青海玉树7.1级地震也造成了近3 000人的伤亡． 因此， 我们必须高度重视地震灾害的防御工作， 同时采取有效措施切实减轻地震所造成的灾害损失．

自然灾害本身不可避免， 但如果能预先采取一些合理的防御措施就可以有效地减少这些灾害造成的损失． 为了达到减轻地震灾害的目的， 除了加强城市工程结构抗震设计外， 人们最先想到的就是地震预报． 虽然我国曾有地震短临预报的成功实例， 但现有科技水平还无法彻底攻克这一难题， 地震的预测预报必将长期处于探索和研究阶段． 地震预警是目前世界上公认的能够有效减轻地震灾害的新手段之一．

地震预警系统利用布设在潜在震源区周围的实时传输地震观测台站， 在破坏性地震发生后很短时间内， 根据距离震中较近的若干个触发台站的信息， 迅速测定地震发生的地点和规模， 并利用地震动P波传播速度大于更具破坏性的S波和面波传播速度的特点及地震波传播速度远远小于电磁波的原理， 在破坏性地震动到达之前向周围可能遭受地震影响的地区发出警报． 按照工作方式的不同， 地震预警系统又可分为现地预警、 异地预警及混合预警等3种模式． 由于预警目标区内各地点距离震中的远近不同， 可以利用的反应时间从几秒至几十秒不等， 从而人们可以及时地采取一些必要的紧急避险逃生措施. 例如， 及时告知人们撤离危险场所， 以减少不必要的伤亡； 使正在行驶中的高速列车减速或采取制动措施， 从而防止列车在地震中发生出轨事故； 自动切断燃气供应， 以防止火灾发生； 通知正在从事危险生产的工人停止作业， 以消除潜在威胁， 从而达到减轻地震灾害， 防止次生灾害的目的． 目前， 世界上多个国家或地区已经建设了多个针对特定设施、 单个城市甚至大区域的地震预警系统， 很多已经取得了一定的减灾实效(Kamigaichi et al, 2009; Allen et al, 2007; Iglesias et al, 2007; Erdik et al, 2003; Iannaccone et al, 2010; Hsiao et al, 2009)．

1 预警震级的相关研究

一个完整的地震预警系统至少应该包括实时地震定位， 实时震级计算， 预警目标区烈度估计， 以及预警信息发布等几个重要功能模块． 其中实时震级计算是最重要、 也是最困难的一部分． 由于地震预警应用中对于信息的高度时效性要求， 计算预警震级时往往只有震中附近的少数几个触发台站、 有限时间长度内的信息可以利用， 因此不可能采用常规震级算法进行计算． 同时， 地震预警系统还要求所发布的信息具有足够可靠度， 因此必须发展一些非常规的、 稳定可靠的实时震级计算方法(张红才等， 2012)． 目前， 国际上也已发展了一些实用的实时震级计算方法， 这些方法大致可以归纳为三大类: 与周期/频率相关算法(如τmaxp， τc方法等)、 与幅值相关算法(如Pd方法等)及与强度相关算法(如MI方法等)． 本文分别就主要采用的τc和Pd两种方法介绍如下.

1.1 τc方法

Nakamura (1988)提出了利用实时速度记录计算地震动卓越周期的算法， 即τmaxp方法． Allen和Kanamori (2003), Kanamori等(1997), Kanamori(2005), Olson和Allen (2005), Wu和 Kanamori (2005), Wu和Li(2006)及Wu和Kanamori(2008)采用该方法进行了一系列相关研究． 卓越周期计算公式如下：

式中， Xi=αXi-1+x2i， Di=αDi-1+(dx/dt)2i， τpi为i秒时测定的卓越周期， xi是记录的地面运动速度值， Xi是平滑后地面运动速度导数的平方值， α为平滑参数， 它决定了平滑的速度， 一般取为0.999． Wolfe (2006)研究认为， τmaxp参数是幅值和频率的非线性函数， 对不同的幅值和频率分量有着不同的权重． 此外， 由式(1)也可以看出， 由于采用了单边加权， 因而算法的稳定性也会受到较大影响． Wu和Li (2006)的研究结果表明， τmaxp方法的准确性和稳定性不仅受采样率影响, 还与对记录的预处理过程密切相关． 选用不同的低通滤波器或是不同长度的时间窗, 都会对特征周期的计算结果产生较大影响． 当选用较大的高频截止频率对记录低通滤波处理后， 特征周期计算结果会明显偏小．

针对τmaxp方法存在的局限性， Kanamori (2005)提出了一种改进后的特征周期计算方法， 即τc方法． τc的计算公式如下：

式中， 积分区间[0, τ0]即从台站触发后开始计． 时间窗长度τ0如果选择太长， 那么可以用于预警信息发布的时间就会变少； 而如果τ0选择过短， 又不能比较准确估计地震的规模. Kanamori等(1997)综合考虑各种因素， 建议一般τ0取为3 s． Wu和Kanamori (2005)分析认为， 在该时间段内， 对于6.5级以下的地震， 断层破裂已经基本上完成了， 因此就可以根据τc参数对地震的规模做出比较准确的估计． 而对于更大规模的地震， 由τc参数也能够对地震是否具有破坏性快速做出判断．

运用巴什瓦定律， 对式(2)进一步分析可以得到

其中， 〈f2〉为位移谱关于f2的平均频率． 因此， 计算得到的τc值本质上对应了位移谱重心位置处的周期． 如果假定前3 s内的波形全部为P波信息， 那么采用中小地震的Brune震源模型可以从理论上证明， 此时计算得到的周期参数τc本质上也就是P波的拐角周期； 对于破裂尺度较大的破坏性大震， τc参数也是衡量地震规模的有效参数． 这也就能够比较合理地解释τc方法的物理意义． τc方法是本文主要采用的方法之一．

1.2 Pd方法

由地震动幅值计算地震预警震级的方法似乎与传统的震级计算方法更接近， 也更容易理解． 由于预警系统应用时效性的要求， 可以利用的只是少数几秒钟的信息， 相对于整个波列的信息而言这是相当有限的， 因此采用幅值参数的预警震级计算方法并不简单等同于传统的震级计算方法．

地震动在地球介质中传播时由于几何扩散、 介质吸收、 反射折射等因素的影响， 幅值会随着震中距的增大而衰减， 地震震级的测定过程中就是利用了这种衰减关系． 类似地， 地震预警中若已知P波段(或是P波段前若干秒)的幅值， 那么也可以根据衰减关系推算得到地震的震级． 但正如前文提到的一样， 由于此时可用信息的有限性， 由幅值参数计算预警震级与传统的地震震级计算之间还存在着本质的区别． 相关研究结果表明(马强， 2008)， 相同情况下位移幅值参数的预警震级估计结果优于速度、 加速度参数， 国际上也大都采用位移参数进行预警震级的计算． 同样由于垂直向分量中的P波最为发育， 因此一般也都只从垂直向记录提取位移幅值参数．

Wu等(2007) 将Pd定义为采用2阶高通巴特沃斯滤波器(低频截止频率为0.075 Hz)滤波后的位移记录中P波触发3 s时间内的位移幅值． 利用发生在美国南加州地区地震的强震记录， 他们统计得到了Pd与地震震级间存在的关系并用于地震预警系统中． 他们的研究还表明， Pd与地震动峰值速度PGV、 峰值加速度PGA之间都存在着良好的相关性， 因此Pd参数也可以作为预警目标区地震烈度估计的有利判据之一． Lancieri和Zollo (2008)利用震中距50 km内的强震记录对Wu等(2007)的方法进行了进一步研究， 认为除了采用3 s时间窗内的Pd参数外， 也可以用捡拾到P波后2 s， 4 s时间窗内的P波初始位移的最大值Pd参数， 或是捡拾到S波后1 s， 2 s时间窗内的S波初始位移的最大值Sd参数来预测地震震级． 其他一些学者也都对该方法进行了研究， 结果都表明Pd方法也是稳定性和可靠性都非常好的算法之一， 甚至比τc方法的稳定性还要高一些. 同时Pd的计算也比较简便并容易实现， 因此本研究中也将该方法选为主要采用的方法之一． 但采用Pd方法计算预警震级时， 必须注意到由于对断层的欠采样而引起的震级饱和现象. 下文中作者还将针对这一问题展开探讨．

预警震级计算中实际上包含了一个十分重要的基本物理问题， 即能否由初始破裂的、 极其有限的信息估计整个地震的规模． 目前， 国际上对于这个问题还没有比较统一的认识． Olson和Allen (2005)在《自然》杂志上发表了他们的相关研究成果， 认为采用τc方法在有限的时间内(即使断层破裂过程还未结束)就可以得到整个地震的最终震级， 即使是对于大震这个结论也是成立的． 而Rydelek和Horiuchi (2006)随后对其研究结果表示质疑， 因为他们所进行的研究显示τc与震级之间并不存在明显的相关性． Kanamori (2005)对相关的研究进行了归纳总结， 认为至少有两种理论可以支持由少量初始信息计算地震震级的可行性： 其一是成核震相的观点； 其二是P波、 S波分别携带不同信息的观点． Iio(1992, 1995), Umeda (1990, 1992), Ellsworth和Heaton(1994), Ellsworth和Beroza (1995), Nakatani等(2000)研究认为， 地震在一定长度和时间域内的成核过程决定了初始破裂的形态以及地震的最终规模， 至少在统计意义上这一观点是成立的． 而Kanamori(2005)则认为， 地震产生的P波携带着地震本身的信息(P波段的波形能反映断层是怎样滑动的)， 而S波则携带着地震能量的信息， 因而可以分别加以利用． 总而言之， 利用初始破裂的有限信息估计整个地震的规模是有一定理论依据的． 而其他一些研究者， 如Mori和Kanamori (1996), Kilb和Gomberg (1999)则持相反的观点． 他们认为不同规模的地震， 初始破裂过程与最终震级之间并不存在关联性． 关于这一问题的争论还在继续， 但近几年来发表的一系列相关论文却从一个侧面证明， 由初始破裂信息估计整个地震的规模是可能的， 采用一些特殊、 稳定的算法是能够获得足够精度的震级估计结果的．

初始P， S波段的位移幅值与地震震级间的相关性， 可以用基本的震源理论进行解释． 假定位移幅值Pd仅取决于记录信号中相对高频的成分， 地震计与破裂断层之间具有一定距离， 断层破裂的方向性和辐射方式等因素的影响被处于不同方位的台站平均化， 此时， 地震辐射的形式即可一阶近似为点源模型远场效应． 因此， 震中距为R处的P， S波的位移场u(t)与地震矩速率 的关系式可写成

式中， c为波速， Δ 为断层平均滑动速率， Σ为断层中初始阶段时的滑动断层面积， L为断层线性尺度， C为一阶几何参数， const为常数． 根据断层动力学理论模型， 滑动速率幅值Δ 与动应力降Δσ线性相关．

另一方面， 地震破裂的发育、 传播由弹性能能流密度控制, 即

式中， f是由破裂速度νr和荷载情况决定的无量纲函数， μ为刚度系数．

根据以上分析， 远场位移u(t)和能流密度G都由应力降Δσ和断裂的几何尺寸(参数L)决定． 因而大部分的观点认为， 具有较大初始能量的破裂能够传播到更远的距离． 因此， 初始P， S波位移幅值与震级间的相关性可以认为是地震震级与初始阶段应力降水平和(或)断层初始破裂面积间的关系． 当然， 断层破裂的传播还受到断层区域内相对强度等其它因素的影响， 因此， 上述结论只是建立于概率意义上(并非确定性的)．

如果假设断层破裂面积Σ随震级的增大而单调增加， 那么也就容易得到滑动时间(或上升时间)也与震级的大小密切相关． 为了解释观测得到的特征周期τmaxp与震级间存在的关系性， Olson和Allen(2005)进一步深化了该假设， 并认为特征周期与初始破裂阶段的滑动持时相关． τmaxp参数与地震震级间的统计关系也能够进一步作为滑动面积与震级之间存在关联性的证据， 即使是在破裂的初始阶段这种关联性也是存在的．

2 本文统计所用震例及记录挑选

为了能够准确地由初始P波段的特征信息估计整个地震的规模， 需要有一定数量的地震记录进行统计. 统计样本中的震级区间应尽量完备， 每个震级段内的事件数目也应尽量 均匀． 本研究中， 作者一共收集了55个日本KiK-net台网记录到的事件(Mjma4.0—Mjma7.3) 和87个汶川余震事件(MS3.5—ML6.3). 本研究忽略不同震级标度间的区别， 统称为M． 所用地震事件的空间分布及发震时刻、 震级大小、 震中位置等其它详细资料分别如图1—3所示． 对强震记录进行基线校正等基本处理后， 本文采用金星等(2003)和马强等(2003)提出的实时仿真算法， 分别将加速度型记录仿真为速度时程和位移时程等， 以满足计算需要． 随后， 又采用人工识别方法对各台的记录进行P波到时捡拾．

图1 本文选用的KiK-net台网记录的地震事件 Fig.1 KiK-net event records used in this study

图2 本文选用的汶川余震地震事件 Fig.2 Wenchuan aftershocks used in this study

图3 本文所用事件和记录情况统计 (a) 所用记录的震中距情况； (b) 所用事件的震级情况； (c) 所用数据的震中距-震级分布情况 Fig.3 Events and records used in this study (a) Number of records versus epicenter distance; (b) Number of earthquakes versus magnitude; (c) Earthquake magnitudes versus epicenter distances

众所周知， 垂直向地震动分量的P波最发育、 震相最清楚， 因此, 通常从垂直向地震动记录中读取P波到时及计算相关参数(如方位角、 初动方向等)． 现有的地震预警中所采用的一些震级算法也大都是从垂直向地震动记录中提取能反映地震规模的特征参数， 并据此计算地震预警震级． 因此， 本文也主要利用垂直向地震动记录进行相关研究． 相关理论分析也表明， 当震中距超过100 km后, 由于反射、 折射等原因地震记录中震相会变得异常复杂， Pg波也不再是首先到达的震相， 因此就不能比较可靠地提取出相应的特征参数， 震级的估计结果也会因此产生较大误差． 鉴于上述原因, 并充分考虑地震预警对于信息的高度时效性要求， 本文只选用震中距30 km内的台站记录参与统计． 在该震中距范围内， P波段的波形成分较为单纯, 同时又能充分满足时效性要求， 因此也就能够保证地震震级估计结果的可靠性和准确性． 由于地震成因的复杂性， 即使对于同一次地震， 由不同台站的观测记录中提取出的特征参数之间也会存在较大的差异． 但如果对同一次地震中多个观测台站特征值结果取平均， 再用该平均值进行相关统计分析， 震级估计结果的离散性将得到大大改善， 这样也就能够有效地排除单个台站特征值计算结果跳动引起的震级估计结果的不确定性． 实际工作中作者还注意到， 所使用的记录还需满足信噪比的要求， 即所选记录必须具备足够的可用性． 实际上， 尤其对于大震而言， 近场台站的信噪比是肯定能够满足要求的， 信噪比的挑选主要还是为了保证中小地震记录的可靠性． 因此， 本文参考Wu和Li (2006)的记录挑选标准对记录再次进行筛选， 只选用P波触发后3 s时间内加速度记录的幅值超过5×10-2 m/s2， 同时该时间窗长度内的平均幅值与触发前平均幅值之比超过10的强震记录参与相关统计． 根据上述记录挑选原则， 本文共挑选出253条垂直向强震记录参与统计(图3)． 3 预警震级统计结果

本研究基于这些地震事件记录， 采用τc和Pd两种方法提取出地震动P波段的周期和幅值特征参数， 然后再分别采用M=a+blg(τc)和M=a+blg(D)+clg(Pd)两种形式拟合， 得到了特征参数与地震震级间的关系． 统计结果分别如下.

3.1 τc方法

从图4和式(6)的结果中可以看出， 周期参数τc能够较好地反映地震规模， 由τc参数计算地震震级的方差为0.62个震级单位． 总体而言， 由汶川余震和KiK-net台网记录中提取得到的τc结果有着明显的规律性， 但个别记录的计算结果却明显地偏离了这一趋势． 正如前面提到的一样， 即使是同一次地震， 由不同台站计算得到的τc值之间也会存在着一些差异性． 因此， 本文对每个震级段中各个τc计算结果求取平均， 并再次统计得到了τc平均值与地震震级间存在的关系， 结果分别如图5和式(7)所示．

图4 τc参数与地震震级间的统计关系

Fig.4 Statistical relationship between τc and magnitude

图5 τc平均值与震级间的统计关系 Fig.5 Statistical relation between mean τc value and magnitude

从图5统计结果中可以看出， 当采用τc平均值计算地震震级时， 统计关系的离散性明显得到了改善， 此时震级计算的方差减小为0.50个震级单位． 由于震级为6.2时只有两个台站的记录， 所以该平均值偏离得较多． 图4， 5中作者还分别添加了四川汶川主震中首先触发的PXZ台(震中距36 km)和DXY台(震中距47 km)的τc计算结果和台湾集集主震中震中距30 km内的13个触发台站的τc计算结果， 以及由这些计算结果得到这两次地震事件的τc平均值(分别为2.04 s和2.51 s， 均没有参与相关公式的统计)． 从这些对比分析中可以发现， 当震级超过7.0时(包括KiK-net台网记录的一次震级为7.3的事件， τc平 均值为1.87 s)， 会存在一个比较明显的震级 “饱和”现象． 此时， 如果再用周期参数τc计算震级， 就会明显地低估震级． Kanamori (2005)采用Sato圆盘破裂模型对τc方法进行数值模拟后认为， τc方法只对震级6.5以下的地震有效， 而一旦超过这一界限， 仅靠这3 s时间窗内的信息将无法反映整个破裂断层的规模， 因此， 不可能利用τc参数对震级做出准确的估计， 即对断层的欠采样是造成震级“饱和”现象的主要原因．

从本文相关研究中也可以看出， 周期参数τc在不同区域间的差异性也较小． Wu和Li (2006)及Wu和Kanamori (2008)即利用发生在台湾地区、 南加州以及日本等3个地区的54次地震事件(震级4.1—8.3)， 同样统计得到每次地震中首先触发的前6个台站τc计算结果的平均值与矩震级MW间的关系

图6 不同统计关系式的对比Fig.6 Comparison between different statistical relationships

将本文统计得到的式(6)和(7)与式(8)进行对比， 如图6所示.

从图6的对比中可以看出， 由于研究过程中所使用的样本不同， 3个关系间还是存在着少许的区别. 相同情况下， 由本文关系计算地震震级时结果将稍小于Wu (2008)推荐关系式的结果． 同时从该图中也可以明显观察到， 上述3条线的斜率基本是一致的， 这也说明τc参数受到地区性的影响较小， 因此τc参数是一个十分稳定可靠的参数． 综上所述， 由周期参数τc能够对地震震级做出较好的估计．

3.2 Pd方法

本文借鉴Wu 和Kanamori(2005)的相关研究， 选用与之相同的定义计算Pd参数． 但与τc方法不同的是， 采用Pd方法计算预警震级时还需用到地震预警定位的结果， 因此， 为保证本文结果的可靠性， 作者首先仔细核对了每次地震事件的定位结果， 随后再进行相关统计分析. 统计结果如图7和式(9)所示．

图7中， 纵轴P10 kmd表示将不同震中距处的位移幅值Pd统一校正到震中距为10 km的“参考”场地时的位移幅值， 这样处理主要是为了画图方便． 采用Pd方法的震级估计方差为0.61个震级单位， 与τc方法的计算精度相当． 同时， 从图7中也可以明显观察到， 当震级小于6.5时， 绝大多数参数都处于1倍方差之内， 而当震级超过6.5后也会出现显著的“饱和”现象， 此时采用Pd算法的预警震级估计结果也会明显“低估”实际震级． Zollo等(2006)研究结果表明， 不仅采用Pd参数时会产生饱和现象， 采用S波前2 s的幅值Sd计算预警震级时同样也会产生类似的现象， 只不过产生饱和现象时的震级会更大些． 利用P, S波段的位移幅值计算预警震级时所出现的不同饱和现象, 可以用等时线理论进行解释． 通常而言， 前若干秒内辐射出P, S波的断层面积要大于等于相同时间内破裂断层的面积(这取决于地震仪与断层间的相对位置)． 而P, S波产生后的很短的时间内， 等时线即覆盖了辐射出高频地震动的断层区域． 但由于S波速度小于P波， 相同持续时间内由S波等时线包围的断层面积要大于P波等时线包围的断层面积， Sd的饱和震级也因此较Pd更大些．

图7 Pd与震级间的统计关系 Fig.7 Statistical relation between Pd and magnitude

Murphy和Nielsen (2009)研究认为， 对于饱和震级(小于6.5级)以下的地震， 由P波触发3 s内的等时线包围的面积要大于或等同于最终的断层尺度. 这也就能够解释Pd与震级间存在的统计关系, 而且也从另一方面说明, 当采用窗长为3 s的P波段信息估计震级较大的地震时产生的饱和效应主要是由于对断层的欠采样造成的． 当时间窗长度扩展为4 s或者更长或是选用更长时间窗内的S波幅值信息时， 将能对更大的断层进行采样， 此时饱和震级下限甚至能扩展到7.5级． 然而， 对于震级更大的地震， 由于饱和效应更为显著， 因此仅由P, S波初始阶段的信息还是会低估实际震级．

因此， 为消除震级饱和现象的影响， 本研究从台站触发开始即对Pd进行连续追踪测定， 从而实现了震级连续测定． 图8为捡拾到P波后时间窗长分别为4—10 s时的Pd与震级间的统计关系. 同样采用M=ai+bilg(D)+cilg(Pd, i)的形式进行回归拟合. 其中Pd, i为各时间窗长内的位移幅值Pd值; ai, bi, ci为系数(图中集集主震和汶川主震的记录并未参与相关统计)． 各时间窗长内统计关系式的系数如表1所示．

表1 不同时间窗长内回归公式的系数 Table 1 Coefficients of statistical relations for each time window

图8 采用Pd方法对震级进行连续测定. (a)—(g)分别为时间窗长4—10 s时的统计关系 Fig.8 Successive magnitude estimation using Pd method. (a)—(g) show the relations with time windows being 4—10 s, respectively

由以上统计结果可以看出， 随着时间窗长度的增加， 位移幅值Pd与地震震级之间的相关性明显增强， 震级饱和现象越来越不明显， 震级估计的方差也随之逐渐减小． 当时间窗的长度达到10 s时， 采用Pd方法的震级估计方差已经减小为0.37个震级单位， 震级饱和现象也基本上观察不到． 采用本文方法即可对震级进行连续追踪计算， 并在10 s左右时获得稳定可靠的震级估计结果． 由于本研究中选用记录的震中距都小于30 km， 因此可忽略震中距的影响， 对各个时段内的Pd取平均， 然后再采用M=ai×lg(Pd, m)+bi±std的形式回归统计该平均值与震级间的关系． 结果分别如图9、 表2所示．

表2 Pd平均值与震级间统计关系式的系数 Table 2 Coefficients of different statistical relations between mean Pd value and magnitude

图9 不同时间窗长内Pd平均值与地震震级间的统计关系. (a)—(h)分别为窗长3—10 s时的统计关系 Fig.9 Statistical relations between mean Pd value for different time windows and earthquake magnitudes. (a)—(h) show the relations with time windows being 3—10 s, respectively

分析以上统计结果可见， 由各震级段中Pd平均值估计地震震级时， 数据的离散性明显减小， Pd平均值与震级间的相关性也随着时间窗长度的增加而增强． 当时间窗长度为10 s时， 震级估计的方差仅为0.28个震级单位．

前文中也已经提到， Pd参数还与地震动峰值速度PGV、 峰值加速度PGA都存在着良好的相关性， 其中与PGV之间的相关程度最高． 通过Pd与PGV之间关系即可以对记录台站处地震动强度迅速做出估计， 据此即可研发“现地预警”模式的相关算法． 本文采取与震级类似处理思路， 分别统计得到了时间窗长度分别为3—10 s时位移幅值Pd与峰值速度PGV之间的统计关系， 结果分别如图10和表3所示．

表3 不同窗长的Pd与PGV统计关系式的系数 Table 3 Coefficients of different statistical relations between Pd and PGV for different time windows

图10 不同时间窗长内Pd与PGV间的统计关系. (a)—(h)分别为窗长3—10 s时的统计关系 Fig.10 Statistical relations between Pd value for different time windows and PGV. (a)—(h) show relations with time windows being 3—10 s, respectively

由以上统计结果可以看出， 无论是日本KiK-net台网记录还是汶川余震记录， Pd与PGV之间都存在着良好的相关性， 而且这种相关性也不存在区域性的特点． 即使时间窗的长度仅为3 s， 由位移幅值Pd也能够对峰值速度PGV做出比较准确的估计（统计关系式方差为0.37）． 同时， 随着时间窗长度的增加， Pd与PGV的相关性越来越高， 由Pd估计PGV的方差也越来越小， 当时间窗长度为10 s时统计关系式的方差仅为0.26． 因此， Pd参数也可用于对地震动峰值PGV的持续跟踪预测． 但从上图中也能明显观察到， 对于集集主震和汶川主震等大地震而言， 位移幅值Pd与PGV之间的关系与其它中小地震的记录之间也存在着显著的不同. 当时间窗长度较短时， 整个断层破裂过程还远没有完成， 因此也不可能由Pd参数对地震动的强度做出准确估计， 且采用本文上述关系对大震级时的PGV进行估计时会造成明显的低估， 因此， 本文统计关系的适用性还需进一步探讨， 大震级时Pd与PGV之间的关系还需要收集更多资料进行统计分析． 综上所述， Pd方法有着多重的优越性， 不仅可以用于震级的测定， 而且也能够对地震动强度进行连续预测． 4 特征参数间相容性检测

由于公众对于地震信息普遍存在恐慌心理， 一旦发布地震预警警报， 势必会对社会秩序和人们的生产生活产生重要影响， 因此, 必须充分确保预警信息的可靠性． 这就需要在技术上采取一些特别考虑， 保证每次发布地震预警信息的准确性， 尽量避免因为错报而使公众产生的恐慌．

Wu 和Kanamori (2008)研究认为， τc和Pd参数都可以作为地震预警信息发布的粗略判据． 他们发现， 如果一次地震中计算得到的τc超过1 s， 同时Pd又超过0.5 cm， 那么就可以大致断定这次地震很可能是一次破坏性地震， 据此就可以发布预警信息． 同时τc和Pd也可以结合在一起对地震事件进行判定， 如果它们之间满足相关的统计关系， 那么就能够比较可靠地鉴别出需要发布地震预警信息的破坏性地震事件．

从本文相关分析中可以看出， τc参数实际上反映了地震动的周期（频率）特征， 而Pd参数则反映了地震动的幅值（强度）特征． 本文结合两个参数的不同属性特征， 认为对于一次真正需要发布警报信息的地震事件， 它的周期（频率）特征和幅值（强度）特征应该是充分相容的， 即该事件的频谱应足够宽， 同时强度也应足够大． 经过τcPd相容性检测也就能够排除相当大一部分Pd很大而τc很小（爆破或干扰信息）、 或τc很大而Pd很小（远震）的事件的影响． 本文统计得到了τc与Pd之间的关系， 如图11和式（10）所示．

图11 Pd与τc间的统计关系

Fig.11 Statistical relation between Pd and τc

从图11中可以看出， τc参数和Pd参数具有良好的协调性， 绝大部分的记录都处于2倍方差之内． 根据方差分析理论， 本文认为对于处于1倍方差之内的散点即可以判定为一次“确定”的地震事件； 对于处于1倍方差与2倍方差之间的， 则认为是一次“可能”的地震事件； 而一旦超出2倍方差则认为是“不可能”的地震事件． 通过τc与Pd相容性检测即可以剔除一些异常信息的干扰， 保证信息的可靠性． 同时注意到， 对于绝大部分震级超过6.0的地震， τc值都超过1 s， 但Pd值只超过0.1 cm， 这与Wu和Kanamori（2005）的研究结论不太一致. 或许这些差异是由本文统计中选用的大震级事件偏少而造成的， 作者还将继续针对该问题进行研究．

由Pd参数的定义中可知， Pd参数仅能反映P波前3 s内地震动的最大强度而不是平均强度． 为保证发布预警信息的可靠性， 作者提出如下地震动平均强度参数速度均方根vrms． 速度均方根vrms采用下式计算

其中vEW， vNS， vUD分别为东西、 南北、 垂直三向速度时程， 即vrms反映了P波触发3 s内的地震动强度的平均水平． 随后， 本文又统计得到了Pd与速度均方根vrms之间的关系， 如图12和式(12)所示．

图12 Pd与vrms间统计关系

Fig.12 Statistical relation between Pd and vrms

从图12中可以看出， 由于vrms与Pd都是强度参数， 它们之间的相关性明显好于τc与Pd间的相关性． 同理， 本文认为处于1倍方差之内的是一次“确定”的地震事件； 处于1倍方差与2倍方差之间的则认为是一次“可能”的地震事件； 而一旦超出2倍方差则认为是“不可能”的地震事件． 经过vrms与Pd的相容性检测， 也能有效地排除一大批异常信息的干扰． 能够通过上述两个相容性检测后发布的预警信息才具有足够的可信度和准确度．

5 讨论与结论

地震预警作为一种新颖的、 能够减轻地震灾害的有效手段， 过去20几年间已经受到越来越多国家的重视， 很多国家和地区也都已经建立了自己的地震预警系统， 或是正在对这些系统进行测试研究． 本文即简要介绍了地震预警技术在防震减灾中的重要作用以及地震预警系统在世界各地的建设情况， 并主要针对地震预警系统中预警震级的计算方法进行了研究． 借鉴国际上成熟的τc和Pd两种算法， 作者利用日本KiK-net台网和四川汶川余震 的强震观测记录， 分别统计得到了上述两个特征参数与地震震级之间的关系， 并对两种方法进行详细对比分析． 为保证预警震级计算结果的稳定性和可靠性， 作者还根据两个参数的不同属性， 提出了τc与Pd参数相容性检测以及Pd与速度均方根vrms参数的相容性检测等． 从本文相关分析结果中可以看出， 利用以上两个参数都可以对地震震级做出比较准确的估计． 本文还拓展了Pd方法的应用范围， 并从3 s开始连续对其跟踪计算， 这样也就能够实现地震预警震级的连续测定， 当时间窗长度达到10 s时即可获得稳定可信的震级估计结果． 同时， 采用Pd参数也能够对速度峰值PGV进行连续估算， 并可以根据PGV估计结果及时发布地震预警信息．

必须特别强调的是， 地震震级确定中还涉及到震源过程、 传播介质、 场地条件、 仪器响应等多个方面， 而地震预警中可用的信息十分有限， 仅仅依靠幅值相对较小的P波段单一的周期参数或幅值参数很难得到可靠的、 稳定的地震震级估计结果． 实际的地震记录是非常复杂的， 单一的周期参数或幅值参数在一定程度上能够反映地震的规模， 但不可能反映地震的全部特征， 这也决定了本文方法用于预警震级估计时其结果的离散性． 同时， 采用两种不同参数估计得到的预警震级间可能还会存在一些差异性， 因此， 如何更合理地综合运用两种参数， 或者发展新的特征参数， 也是非常值得研究的问题．

参考文献
[1] 金星, 马强, 李山有. 2003. 四种计算地震反应数值方法的比较研究[J]. 地震工程与工程振动, 23(1): 18--30.(1)
[2] 马强, 金星, 李山有. 2003. 单自由度系统地震动力反应的实时计算方法[J]. 地震工程与工程振动, 23(5): 61--68.(0)
[3] 马强. 2008. 地震预警技术研究与应用[D]. 哈尔滨: 中国地震局工程力学研究所: 124--134.(1)
[4] 张红才, 金星, 李军, 韦永祥, 马强. 2012. 地震预警震级计算方法研究综述[J]. 地球物理学进展, 27(2): 464--474.(1)
[5] Allen R M, Kanamori H. 2003. The potential for earthquake early warning in southern California[J]. Science, 300: 786--789.(1)
[6] Allen R M, Brown H, Hellweg M, Kireev A, Neuhauser D. 2007. Real-time test of the ElarmS earthquake early warning methodology[J]. J Geophys Res, 112(B08311), doi:10.1029/2008GL036366.(1)
[7] Ellsworth W L, Beroza G C. 1995. Seismic evidence for an earthquake nucleation phase[J]. Science, 268: 851--855.(1)
[8] Ellsworth W L, Heaton T. 1994. Real-time analysis of earthquakes: Early-warning systems and rapid damage assessment[J]. Sensors, (April): 27--33.(1)
[9] Erdik M, Fahjan Y, Ozel O, Alcik H, Mert A, Gul M. 2003. Istanbul earthquake rapid response and the early warning system[J]. Bull Earthq Eng, 1(1): 157--163.(1)
[10] Hsiao N C, Wu Y M, Shin T C, Li Z, Teng T L. 2009. Development of an earthquake early warning system in Taiwan[J]. Geophys Res Lett, 36: L00B02, doi:10.1029/2008 GL036596.(1)
[11] Iannaccone G, Zollo A, Elia L. 2010. A prototype system for earthquake early-warning and alert management in southern Italy[J]. Bull Earthq Eng, 8(5): 1105--1129.(1)
[12] Iglesias A, Singh S K, Santoyo M A, Pacheco J, Ordaz M. 2007. The seismic alert system for Mexico City: An evaluation of its performance and a strategy for its improvement[J]. Bull Seism Soc Amer, 97(5): 1718--1729. (1)
[13] Iio Y. 1992. Slow initial phase of the P-wave velocity pulse generated by micro-earthquakes[J]. Geophys Res Lett, 19(5): 477--80.(1)
[14] Iio Y. 1995. Observations of the slow initial phase generated by micro-earthquakes: Implications for earthquake nucleation and propagation[J]. J Geophys Res, 100(B8): 15333--15349.(1)
[15] Kamigaichi O, Saito M, Doi K, Matsumori T, Tsukada S, Takeda K, Shimoyama T, Nakamura K, Kiyomoto M, Watanabe Y. 2009. Earthquake early warning in Japan: Warning the general public and future prospects[J]. Seism Res Lett, 80(5): 717--726.(1)
[16] Kanamori H, Hauksson E, Heaton T. 1997. Real-time seismology and earthquake hazard mitigation[J]. Nature, 390: 461--464.(2)
[17] Kanamori H. 2005. Real-time seismology and earthquake damage mitigation[J]. Annu Rev Earth Pl Sci, 33: 195--214.(5)
[18] Kilb D, Gomberg J. 1999. The initial subevent of the 1994 Northridge, California, earthquake: Is earthquake size predictable?[J]. J Seismol, 3(4): 409--420.(1)
[19] Lancieri M, Zollo A. 2008. A Bayesian approach to the real-time estimation of magnitude from the early P and S wave displacement peaks[J]. J Geophys Res, 113: B12302, doi:10.1029/2007 JB005386.(1)
[20] Mori J, Kanamori H. 1996. Rupture initiations of microearthquakes in the 1995 Ridgecrest, California, sequence[J]. Geophys Res Lett, 23: 2437--2440.(1)
[21] Murphy S, Nielsen S. 2009. Estimating earthquake magnitude with early arrivals: A test using dynamic and kinematic models[J]. Bull Seism Soc Amer, 99(1): 1--23.(1)
[22] Nakamura Y. 1988. On the urgent earthquake detection and alarm system (UrEDAS)[C]//Proceedings of Ninth World Conference on Earthquake Engineering, 7: 673--678.(1)
[23] Nakatani M, Kaneshima S, Fukao Y. 2000. Size dependent microearthquake initiation inferred from high-gain and low-noise observations at Nikko district, Japan[J]. J Geophys Res, 105(B12): 28095--28109.(1)
[24] Olson E L, Allen R M. 2005. The deterministic nature of earthquake rupture[J]. Nature, 438: 212--215.(3)
[25] Rydelek P, Horiuchi S. 2006. Is earthquake rupture deterministic?[J]. Nature, 442: E5--E6, doi:10.1038/nature04963.(1)
[26] Umeda Y. 1990. High-amplitude seismic waves radiated from the bright spot of an earthquake[J]. Tectonophysics, 175(1-3): 81--92.(1)
[27] Umeda Y. 1992. The bright spot of an earthquake[J]. Tectonophysics, 211(1-4): 13--22.(1)
[28] Wolfe C J. 2006. On the properties of predominant-period estimators for earthquake early warning[J]. Bull Seism Soc Amer, 96(5): 1961--1965. (1)
[29] Wu Y M, Kanamori H. 2005. Experiment on an onsite early warning method for the Taiwan early warning system[J]. Bull Seism Soc Amer, 95(1): 347--353.(3)
[30] Wu Y M, Li Z. 2006. Magnitude estimation using the first three seconds P-wave amplitude in earthquake early warning[J]. Geophys Res Lett, 33: L16312, doi:10.1029/2006GL026871.(4)
[31] Wu Y M, Kanamori H, Richard M A, Hauksson E. 2007. Determination of earthquake early warning parameters, τc and Pd, for southern California[J]. Geophys J Int, 170: 711--717.(2)
[32] Wu Y M, Kanamori H. 2008. Development of an earthquake early warning system using real-time strong motion signals[J]. Sensors, 8(1): 1--9.(5)
[33] Zollo A, Lancieri M, Nielsen S. 2006. Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records[J]. Geophys Res Lett, 33: L23312, dio:10.1029/2006GL027795.(1)