基于地震预警参数阈值评估内蒙古地区地震潜在破坏区

郝美仙 王鑫 张珂 刘颖 尹战军 张建中 郑钰

1)内蒙古自治区地震局，呼和浩特 010010 2)中国地震局地球物理研究所，北京 100081

0 引言

地震是地球上主要的自然灾害之一，其突发性极强，破坏性巨大，一旦在城市或人口稠密地区附近发生破坏性地震，将会对重大基础设施、生命线工程以及人民的生命财产造成严重危害。因此，科学有效地做好地震监测和灾害防御工作，可以最大限度地减轻地震灾害所造成的损失(蒋长胜等，2016)。

地震预警是有效减轻地震灾害损失的有效手段之一，其原理是在地震发生后，抢在破坏性地震波到达之前，向可能遭受地震影响的目标区发出警报(Zhang et al，2016；马强，2008)，为相关行业紧急处置和人们逃生避险预留时间。地震预警信息发布的核心之一就是在破坏性地震发生过程中，能够快速准确地测定出潜在破坏区的范围(张红才等，2013)。

关于破坏性地震发生后如何快速估算预警目标区的潜在破坏区，国内外多名学者已基于地震动参数展开了相关研究(Teng et al，1997；马强，2008；马东，2013；宋晋东等，2017；王玉石等，2010；陈鲲等，2012)。地震动参数主要包括峰值加速度PGA、峰值位移PGD和峰值速度PGV等，由这些参数产出的潜在破坏区空间分布图依赖于地震基本参数和地震事件持续时间，过程耗时几分钟到十几分钟，时效性不强，因此有必要采用一种时效性强的方法来快速获取破坏性地震的潜在破坏区范围，以便更好地服务于社会公众及城市重要设施。

目前我国多地已建立了地震预警示范区，这些示范区依托于现有的地震观测台网，采用将实时数据流传输到台网中心进行统一处理的方式，对于地震预警算法的研究主要集中在震级估算与目标区的破坏情况研究2个方面。对于震级估算，目前主要使用2种预警参数(位移幅值Pd和特征周期τc)，国内外众多学者对此展开了研究，并获取到了适应于特定地区的经验关系式。对于目标区的破坏情况研究，Zollo等(2010)曾提出一种综合的地震预警方法，该方法实时联合测量P波到达后3s窗内的Pd与τc值，通过与预设的阈值相比较，达到快速区分是否为破坏性地震的目的；彭朝勇等(2019)利用该方法对四川庐山和云南鲁甸2次地震的潜在破坏区范围进行计算，计算结果与实际烈度调查结果基本吻合。地震预警参数的经验关系式具有区域相依性，与内蒙古地区地壳速度结构、断裂带分布及发震地区附近的台站分布情况等因素有关。至今适合内蒙古地区的地震潜在破坏区预警参数的经验关系式的相关研究工作还未开展。

基于此，本文采用内蒙古地区测震台站记录到的地震事件数据，拟合出适应于内蒙古地区的预警参数关系式，包括特征周期τc与震级ML、位移幅值Pd与τc和震源距R的统计关系；根据评估内蒙古地区地震潜在破坏区范围时定义的相应预警参数阈值，利用2次已有地震开展线下模拟，对拟合关系式进行验证。

1 数据选取与预处理

本研究选取2016—2019年内蒙古及周边地区ML≥2.0地震，共计251次，所选地震事件参数信息直接采用全国统一地震编目(1)中国地震台网中心地震编目系统. 内部资料.结果，波形数据来源于内蒙古测震台网波形数据库，文中所用震级标度均为ML，所选事件和记录台站如图1 所示。

图1 本文选取的地震事件与台站分布

由于所选该时间段内发生的中强震事件较少，仅有1次ML5.4地震和1次ML6.1地震，因此在对获取到的内蒙古地区的测震波形数据进行筛选处理时，采用了与彭朝勇等(2019)不同的规则，限制了震级下限、记录台站个数和震源距范围。具体筛选条件为：①所选地震事件震级ML≥2.0；②有3个以上台站可以清晰记录到该事件；③震源距在200km范围内。通过以上条件共筛选出251个地震事件和3075条台站记录，本文所用震例分布情况见图2。

图2 本文所用震例分布情况统计(a)研究所用事件的震级-频度分布；(b)研究所用记录的震级-震源距分布

对于筛选后的台站记录数据，首先使用Matlab逐一打开各条记录并手动选取P波和S波到时，保存至文本文件中，然后对各台站垂直向速度记录进行积分后得到位移值，再利用Butterworth高通滤波器对位移值进行连续滤波处理，最后分别根据P波触发后3s窗段和全P波段垂直向位移时程和速度时程计算得出τc和Pd值，以进行对比分析。另外，在处理过程中分别使用了2阶和4阶Butterworth高通滤波器进行滤波处理，以便于考量不同阶数滤波器对结果的影响。

2 相关参数的拟合关系式

2.1 τc与ML的统计关系式拟合

利用选取的所有地震事件记录，采用下式进行线性拟合

lg(τc)=A×ML+B

(1)

式中，A和B为待拟合系数，τc的单位为s。

在τc与ML的关系式拟合过程中，根据震级范围和选取的P波不同窗长分别使用2阶和4阶Butterworth高通滤波器进行数据拟合计算。对比结果显示，ML2.0～6.1、ML2.5～6.1和ML3.0～6.1的记录数据拟合获得的结果，无论使用的滤波器阶数和P波窗长为多少，标准偏差相对较大，相关性微弱，甚至出现了负相关的情况。通过对ML<3.5的地震事件记录进行分析发现，部分记录在P波到达前存在较大的波形扰动情况，从而对计算获得的τc值造成了极大的影响，由个别记录算得的τc值超过了5s。这说明在利用τc参数开展震级估算时，记录的数据需要具有较高的质量。

震级选取在ML3.5～6.1和ML4.0～6.1的记录数据，无论使用的滤波器阶数和P波窗长为多少，均具备较好的相关性，尤其是对于震级范围在ML4.0～6.1的记录，其相关系数均在0.9以上。此外，使用4阶Butterworth高通滤波器处理的结果要优于使用2阶滤波器，这主要归因于高通滤波器的过滤带更窄，可以更好地滤除滤波频带范围外的能量。图3 为震级范围为ML4.0～6.1时，不同窗长的P波段经4阶滤波器处理后所对应的拟合曲线。

图3 特征周期τc与震级ML的拟合曲线实心圆为每次事件所对应各个台站τc的平均值，空心圆为每次事件中各记录台站计算得出的τc值；黑色实线为τc与ML拟合曲线；黑色虚线为1倍标准方差；SDV为标准差；R为相关系数；Pole为所用的滤波器阶数

2.2 Pd与τc和震源距R的统计关系

由于较小的地震不会引起破坏，一般不需要对其进行预警处理，因此本文只分析ML4.0～6.1地震事件波形。对于所选用的震级范围ML4.0～6.1的记录数据，采用不同的垂直向P波窗长和4阶Butterworth高通滤波器进行处理后，经过多变量线性拟合，所获得的Pd与τc和震源距R的相关性分别为

lg(Pd)=1.5122+2.4665×lg(τc)-2.0336×lg(R)

(2)

lg(Pd)=0.3338+2.2531×lg(τc)-1.4228×lg(R)

(3)

式(2)、(3)分别为采用垂直向P波3s窗和全P波窗得出的拟合关系式，式中Pd的单位为cm；τc的单位为s；R为震源距。关系式对应的拟合曲线如图4，从图中可以看出，Pd值明显呈现出随τc值增加而稳步增加的趋势，且利用垂直向全P波段数据得到的Pd值变化范围更为稳定。

图4 Pd与τc和震源距R的相关性拟合曲线(a)P波3s窗；(b)全P波窗

3 台站地震预警阈值设定

在开展潜在破坏区范围估算时，需要定义相应的预警参数阈值，以便实时跟踪每个触发台站的变化情况。按照Wu等(2008)、Zollo等(2010)、彭朝勇等(2019)的阈值定义模式，需要同时获取Pd和τc的阈值。对于Pd阈值，一般根据Pd与PGV的相关性以及PGV与仪器地震烈度的对应关系来获得。采用震级范围ML4.0～6.1的记录数据，同时采用4阶Butterworth滤波器，针对不同的垂直向P波窗长拟合得到的相关性分别为：

垂直向P波3s窗

lg(PGV)=0.5977×lg(Pd)+0.6000±0.3717

(4)

垂直向全P波窗

lg(PGV)=0.7491×lg(Pd)+0.9736±0.3137

(5)

其中，PGV表示台站三分向合成的地震动速度峰值，其单位为cm/s。

图 5为Pd与PGV的相关性拟合曲线，由图可以看出，无论是在标准偏差还是相关系数方面，利用垂直向全P波窗获得的Pd与PGV的相关性要明显优于P波3s窗。因此，在后续的数据处理与分析中，所有的拟合参数关系式均选取震级范围ML4.0～6.1，且采用4阶Butterworth高通滤波器对垂直向全P波窗进行分析。

图5 Pd与PGV的相关性曲线(a)P波3s窗；(b)全P波窗

根据《仪器地震烈度计算暂行规程》(2)《仪器地震烈度计算暂行规程》(中震测发[2015]18号).内部资料.，同时考虑到现有地震事件对应的最大震级较低(ML6.1)，因此设定Ⅳ度对应的三分向合成PGV下限值(1.77cm/s)作为对应的阈值。根据式(5)可得到相对应的Pd阈值为0.1075cm。当设定ML=5时，得出相应的τc值为0.686s。依据Pd和τc的阈值，可定义4类台站地震预警级别：①3级：Pd≥0.1075cm且τc≥0.686s，说明正在发生的地震震级可能很大，并且发生在记录台站附近，预计会对记录台站附近和较远地区造成较大破坏，此为最高警报级别；②2级：Pd≥0.1075cm且τc<0.686s，说明地震发生在记录台站附近且震级较小，只可能对台站周边地区造成影响；③1级：Pd<0.1075cm 且τc≥0.686s，说明地震震级大但震中远离记录台站，仅在震中附近区域会引起破坏；④0级：Pd<0.1075cm且τc<0.686s，地震不会造成破坏。

在地震发生过程中，获取到多个触发台站所对应的Pd和τc实际测量值后，就可以采用与彭朝勇等(2013、2019)一致的方法，通过对各个台站的τc值取平均值，同时将Pd值设定在阈值0.1075cm，根据式(3)可以计算出潜在破坏区的半径。

在实际应用时，首先对台站覆盖区进行固定大小的单元网格划分，单元格内如果有触发台站，则直接利用记录台站Pd和τc的实测值；如果没有触发台站，则使用式(3)和各台站τc测量平均值来推算Pd值，对震中区域进行内插，以填补没有记录台站的区域。生成的等值线Pd=0.1075cm即为评估破坏程度大小的临界线。在中强地震发生时，利用该方法可以根据一定的时间间隔来实时估算破坏区范围，并可随着记录台站的增多逐步更新。

在使用式(3)对Pd值进行计算时，需要先根据P波信息确定地震震中位置与震源深度，以此来计算震源距R值。目前国际和国内用于地震预警的实时地震定位技术已有多种，如“着未着”定位方法(Horiuchi et al，2005)、RTLoc定位方法(Satriano et al，2008)、前三/四台P波到时信息实时定位方法(张红才等，2011)及四台连续定位方法(李军等，2014)等。由于地震定位方法不是本文主要的研究内容，因此在进行线下模拟时，将直接采用地震目录中的震中信息和震源深度，以3台触发的1s时间点作为起始时间。

4 线下模拟与结果分析

利用内蒙古测震台网记录到的阿拉善左旗2015年4月15日ML6.1和2017年6月3日ML5.4地震开展线下模拟(表1)。对记录台站覆盖区域按照0.25°大小划分为均匀的单元网格，每条事件记录数据从P波触发后开始，每隔1s计算一次Pd和τc值，直到S波到达后停止计算，线下模拟的2次地震事件相应的潜在破坏区估算结果如图6 所示，每个事件对应的4张图分别为发震后不同时刻计算得出的Pd分布图。

表1 线下模拟的阿拉善左旗2次中强地震事件参数

图6 2次地震在P波触发后不同时间段的潜在破坏区实验结果五角星为震中位置；三角形为触发台站

从模拟结果可以看出，由于该地区测震台站分布非常稀疏，无法获取到地震发生后完整的潜在破坏区范围。因此，在地震发生过程中采用该方法实时估算潜在破坏区时，需要台网具有一定密度，如能均匀分布，效果更好，否则就会出现结果仅由经验关系式决定、而实际观测值没有发挥作用的问题。

5 讨论与结论

本文利用内蒙古测震台网记录到的2016—2019年地震事件，按照筛选条件选取波形记录，人工拾取P波和S波到时后，基于P波参数阈值开展评估内蒙古地区地震潜在破坏区的研究。该方法所拾取的预警参数包括特征周期τc和位移幅值Pd，首先利用P波3s段和全P波段垂直向位移时程和速度时程计算τc和Pd值，研究中选用不同震级范围、不同垂直向P波窗长和不同阶数Butterworth高通滤波器分别得到内蒙古地区τc与震级ML、Pd与τc和震源距R的统计关系式。

结果表明，在特征周期τc与震级ML的相应关系拟合过程中，震级选取范围在ML4.0～6.1 的事件记录，对其P波3s窗和全P波窗波形数据进行拟合，τc与ML均体现出了较好的相关性，且采用4阶Butterworth高通滤波器处理的结果要优于2阶滤波器；而采用小于ML3.5的地震事件记录拟合结果标准偏差均较大，相关性微弱，说明在采用τc参数开展震级估算时，需要波形记录具有较高的数据质量，且利用垂直向全P波窗数据获取到的Pd值变化范围更为稳定。

震源距的计算受到震源深度不确定的影响，进而会影响经验关系的可靠性。利用已拟合得到的τc与震级ML、Pd与τc和震源距R的统计关系式，同时考虑到内蒙古地区现有地震事件最大震级为ML6.1，按照震级ML=5.0，烈度为Ⅳ度定义了4类台站地震预警级别，利用式(3)获得了计算内蒙古地区潜在破坏区的半径范围的关系式。需要注意的是，这里仅以Ⅳ度对应的单水平方向PGV下限值计算Pd阈值，主要还是由内蒙古地区现有地震记录对应的震级较小导致。今后在实际使用时，如果发生了震级较高的地震，则应将PGV下限值提高，这样会与实际计算获得的潜在破坏区范围更加吻合。

在地震发生过程中，依据地震预警的实时定位算法获得震中位置、震级和震源深度，并通过获取到多个触发台站所对应的Pd和τc平均值后，计算潜在破坏区的半径范围。在计算时如只选取少量震中附近触发台站信息，则误差较大，但时效性强；而选取更远、更多的触发台站信息，估算结果会更为准确，但预警时间也会增加。

本文利用P波全波段的计算结果和经验关系仅限于理论研究及线下模拟，效果较好，但在预警阶段的实际作用有限，可主要用于在烈度速报阶段配合实际仪器估算烈度的应用，便于提高预警参数估算结果的准确性。

本文研究获得的内蒙古地区预警参数经验关系、阈值，由于区域台站分布稀疏和选取地震事件的震级分布范围不均匀且缺少大地震事件记录数据，计算结果与其他地区存在差异。

随着“国家地震烈度速报与预警工程”内蒙古子项目的建设，在未来2～3年内，台网密度将会得到较大程度的提高，内蒙古地震预警项目新建或改造共计61个基准站、133个基本站和265个一般站。项目建成后将在呼和浩特、包头、鄂尔多斯和乌兰察布地区形成台间距为14km的预警监测网，在其他一般预警区形成台站间距为60～70km的监测网(图7)。

图7 “国家地震烈度速报与预警工程”内蒙古子项目建成后台网(a)和现有台站分布(b)对比

将该方法应用到内蒙古近年来发生的阿拉善左旗2015年4月15日ML6.1和2017年6月3日ML5.4地震，结果显示，由于该地区台站分布稀疏，未得到潜在破坏区估算范围，结果仅由经验关系式决定，而实际观测值没有发挥作用。随着“国家地震烈度速报与预警工程”内蒙古子项目的建设，利用该方法实时开展潜在破坏区范围估算所能取得的效果将会大大提高，精度将会得到明显的改善。

致谢：感谢中国地震局地球物理研究所彭朝勇研究员为本文提供计算程序。