-
中国位于欧亚、太平洋及印度洋三大板块交汇地带，新构造运动活跃，地震活动十分频繁。20世纪以来，全球7级以上强震之中，中国约占35%。地震灾害危害性大、涉及面广，对社会经济发展和人们的正常生产生活都有巨大的影响，仅2008年的汶川地震就造成了69 227人遇难，直接经济损失8 452 亿元。为有效减少人员伤亡和经济损失，2010年1月“国家地震烈度速报与预警工程”项目启动，投入20亿元建设5 000个预警台站，组成覆盖全国的国家地震烈度速报与预警系统，以期最大限度减小地震灾害。
地震是地层受到应力作用形成破裂，以地层为介质向四周传播地震波，引起的地震动。根据地震波在介质内部振动方向和传播方向的关系不同，可分为纵波和横波，纵波（P波）的波速约为7.0 km/s，横波（S波）约为3.5 km/s，因此纵波（P波）率先到达台站，也称其为初至波[1]。通过捡拾P波震相计算得到震中位置的信息，推导出更大破坏波横波的到达时间，及时通过电视、广播、短信等方式向社会发布，提高预警意识，及时做好应急措施，紧急制动高速行驶的列车、快速关闭天然气管道、自动关闭家庭用电用水会等。因此，如何快速、准确地捡拾P波震相就显得尤为重要。
-
当震源在A点产生破裂时，向各个方向传播地震波，设传播介质是均匀的，纵波的波速为VP，横波的波速为VS，纵波率先到达台站B点，并被强震仪记录（图1）。
图 1 地震预警示意图
台站接收到初至波的时间可表达为：
$$ {t}_{AB}=\frac{\sqrt{{s}_{2}^{2}+{h}^{2}}}{{V}_{\mathrm{p}}} $$ (1) 同理，可得到震源的S波到达预警目标区的时间，其计算式为：
$$ {t}_{AC}=\frac{\sqrt{{s}_{1}^{2}+{h}^{2}}}{{V}_{\mathrm{s}}} $$ (2) 而经过计算机对初至波的自动捡拾，计算出S波到预警目标区的时刻，最后告知民众。假设这个过程所需要时间为tC，得到预警目标区的避险时间的表达式为：
$$ {t}_{AC}-{t}_{AB}-{t}_{C}=\frac{\sqrt{{s}_{1}^{2}+{h}^{2}}}{{V}_{\mathrm{s}}}-\frac{\sqrt{{s}_{2}^{2}+{h}^{2}}}{{V}_{\mathrm{p}}}-{t}_{C} $$ (3) 由式(3)可知，避险时间不仅与台站B的距离S2和预警目标区C的震中距S1有关，还与预警系统中计算时间有关，而初至波自动捡拾的时间是唯一人为可控的因素。
地震预警又分为单台预警和多台预警，单台预警是指只用一个台站的纵波震相信息来计算地震的主要参数达到预警的作用，而多台预警是综合多台对初至波的捡拾来计算，多台定位的结果更为精细，但所需要的的时间要更长[2]。因此对于台站建设不够密级的人烟稀少的地区一般选择使用单台定位，争分夺秒地提供预警信息。如何有效地实现这一过程，最关键的还是初至波P初动震相的自动捡拾[3]。
-
震相是反映地震波到时、波形、振幅和质点运动特征的总和，地震的震中、震级和发震时刻可通过震相捡拾和计算得到。距离震中最近的强震仪最先记录到的是P波初动，并对其做震相分析，可快速确定震级、震中方位角和震中距等信息。最常用的P波震相捡拾方法主要有4种，分别是利用能量分析的长短时平均（STA/LTA）方法、通过自回归分析的AIC方法、长短时平均（STA/LTA）和AIC综合择优法、以及基于小波系数的W-AIC方法[4-6]。
-
长短时平均比值法是通过短时窗内振幅特征函数的平均值（STA）与长时窗内特征函数的平均值（LTA）之比（STA/LTA）来反映初至波幅值的瞬时变化[7]。为使计算简便，短时窗一般为地震波的一个周期，也可选择半个周期，而长时窗的长度为周期的整数倍[8]。受地脉动和背景噪声的干扰，需要设定一个阈值，当平均值的比值会大于阈值时，才认为捡拾的初动点是准确的。
又因长时窗内特征函数的平均值（LTA）反映的是地脉动和背景噪声，因此受地脉动和背景噪声干扰较大。在地脉动和背景噪声一定的条件下，短时窗和长时窗的时间越长，P波初动点越容易识别，但也意味着需要花费更多的时间，因此需要根据所研究地区背景噪声特点，选取合适的长短时窗长度[9]。
-
AIC即最小信息量准则（an information criterion），反映了P波到时前后2个不同的稳态过程，一个是地脉动与背景噪声的稳态过程，另一个是地震波到达后幅值变化的稳态过程[10]。2个稳态的过程中就有一个分界点，该分界点就是初至波到时点，在分界点前后创建最大似然函数，AIC的核心表达式为
$${K_{{\rm{AIC}}}} = - 2\log \left( {\text{最大似然函数}} \right) + C\left( {\text{参数}} \right)。$$ (4) 由地震记录的数据计算AIC函数，通过多次拟合对比，得到最大似然函数的最大值，该最大值对应的点即AIC函数最小值的点，AIC最小值点处捡拾P波到时点。但时间序列越长，相关信息就会越分散，使拟合精度比较高就需要多个自变量进行复杂拟合，因此需要选择合适的时间序列[11]。
-
长短时平均（STA/LTA）方法捡拾到的P波初动点，往往会因为初动幅值较小导致捡拾点而滞后。而AIC方法，计算精度高，但计算过程费的时间更长。为提高自动捡拾的精度和效率，综合采用了长短时平均（STA/LTA）方法捡拾快和AIC方法捡拾准确两种方法的优点[12]。首先，通过长短时平均（STA/LTA）方法捡拾P波初动点，然后在初步捡拾点附近选择一个合适的窗长，在此窗内应用AIC方法精确捡拾，最后得到更为精细的捡拾点。
-
初至波到达强震仪后，信号的幅值发生在地脉动和背景噪声下，突变引起信号的不连续，即初动点就是一个间断点，使信号在初至点处存在奇异性。基于小波系数的W-AIC法就是利用了间断点的奇异性，通过对地震信号进行一维多尺度的离散小波分解得到小波系数[13-14]，不同尺度下的小波系数反映出初动的主要特征，然后应用固定时窗AIC方法捡拾P波到时，最后分析不同尺度下捡拾点差值大小来择优捡拾P波到时点[15]。
-
丹江地震台地处丹江口市东郊，距丹江水库大坝3.6 km。在地质构造上位于南阳-襄樊盆地西缘与南秦岭褶皱带接壤部位，NWW向两郧断裂东段帚状构造内,南秦岭NW向构造带是较大中强震活动带，主要监测丹江口水库及周边地区的地震活动。本文选取湖北丹江地震台强震仪实际记录到的6条地震为例（表1），分别对4种P波震相捡拾方法进行分析。
表 1 运用到的地震记录统计表
地震记录 震中距/km 震级M 发震时刻 震中位置 1 259.6 4.9 2019-12-26T18:36:34 湖北应城 2 216.7 4.5 2018-10-11T15:06:30 湖北秭归 3 95.6 4.0 2014-05-27T21:57:09 湖北房县 4 135.9 4.2 2008-03-24T23:24:37 湖北竹山 5 30.5 3.6 2019-11-30T04:28:48 河南淅川 6 32.9 4.5 2018-02-09T19:00:37 河南淅川 -
根据研究区的特点，选择捡拾短时窗长为0.1 s，长时窗长为2 s，阈值取6，得到6条地震动记录的长短时平均值的比值（STA/LTA）曲线图，并自动捡拾到初动点（图2）。
图 2 地震记录垂直向波形与长短时平均比值原图及局部放大图
表2是6条地震动记录通过STA/LTA方法得到的捡拾点，与人机交互捡拾的初动点的误差统计表。由表2可知，6条地震记录的初动点与人机交互捡拾的初动点相比分别滞后0.13 s、0.25 s、0.21 s、0.18 s、0.22 s、0.27 s，捡拾点均滞后于人机交互捡拾的初动点，而地震记录2和6因为初动幅值小，捡拾误差较大。
表 2 应用STA/LTA震相自动捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.47 0.13 2 36.18 36.43 0.25 3 15.95 16.16 0.21 4 22.71 22.89 0.18 5 5.41 5.63 0.22 6 5.76 6.03 0.27 -
根据地震记录和研究区的特点，选取AIC法固定窗长5 s，移动步长1 s，对6条地震记录进行捡拾，得到AIC值的曲线图，并捡拾到初动点（图3）。
图 3 地震记录垂直向地震波形图以及应用固定时窗AIC震相自动捡拾曲线图
由AIC法捡拾结果的误差分析（表3）可知，地震记录1、3、4、5、6的P波到时捡拾结果分别滞后于人机交互捡拾的初动点0.11 s、0.16 s、0.11 s、0.31 s、0.13 s，地震记录2则提前了0.06 s，并没有明显规律，但比STA/LTA方法捡拾精度高。
表 3 应用固定时窗AIC震相自动捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.45 0.11 2 36.18 36.12 −0.06 3 15.95 16.11 0.16 4 22.71 22.82 0.11 5 5.41 5.74 0.31 6 5.76 5.89 0.13 -
首先，长短时平均（STA/LTA）选择捡拾短时窗长为0.1 s，长时窗长为2 s，阈值取6，找到大于阈值6的初步捡拾点；然后，以该初步捡拾点为原点，往前步长取0.4 s，往后步长取0.2 s，在该时窗范围内应用AIC方法，自动捡拾P波到时。综合长短时平均（STA/LTA）法和固定时窗AIC法的捡拾结果见图4。
图 4 地震记录垂直向应用STA/LTA与AIC震相捡拾图
由STA/LTA法结合 AIC法的捡拾结果（表4）可知，地震记录1、3、4、5、6的P波到时捡拾结果分别滞后于人机交互捡拾到时0.09 s、0.12 s、0.10 s、0.20 s、0.12 s，地震记录2则提前了0.06 s，捡拾误差都在0.20 s以下，进一步提高了捡拾精度。
表 4 结合应用STA/LTA与AIC震相捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.43 0.09 2 36.18 36.12 −0.06 3 15.95 16.07 0.12 4 22.71 22.81 0.10 5 5.41 5.61 0.20 6 5.76 5.88 0.12 -
首先，通过Dbn作为小波滤波器的函数，而这里的n越大，表示消失矩数字越大，小波就越平滑，根据地震记录的特点，本文选使n取2，也就是说结果长度为4[16]；然后用wbmpen作为一维小波降噪的函数，得到降噪的全局正阈值，由小波系数阈值得到降噪后的信号和小波分解结构；再由小波分解结构提取小波系数，并重构原始地震波信号；最后在3个尺度下应用AIC法，固定窗长5 s，移动步长1 s，对重构信号自动捡拾，分析对比不同尺度下捡拾点的差值大小，对捡拾点限制得到最优捡拾点[17]。
限制条件如下：首先，相邻尺度的捡拾点到时相差在0.2 s之内，即尺度1与尺度2、尺度2与尺度3捡拾到时差值要在0.2 s以内；其次，3个尺度下的捡拾点必须距离捡拾窗两端边界的0.4 s之外，若在0.4 s以内，则向前移动1 s重新捡拾，直至满足要求；最后，在尺度2下做最初的捡拾，往前步长取0.4 s，往后步长取0.2 s，在该时窗范围内应用AIC方法自动捡拾P波到时，并以此作为最终的捡拾点[18-19]。
由W-AIC方法对6条地震记录的捡拾结果（表5）可知，地震记录1、3、4、5、6的P波到时捡拾结果分别滞后于人机交互捡拾到时0.11 s、0.11 s、0.09 s、0.15 s、0.05 s，地震记录2则提前了0.09 s。该方法对于6条地震记录精度都有提高，特别是对于地震记录6效果显著。
表 5 应用基于小波变换的AIC值震相捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.45 0.11 2 36.18 36.09 −0.09 3 15.95 16.06 0.11 4 22.71 22.80 0.09 5 5.41 5.56 0.15 6 5.76 5.81 0.05 分析以上4种方法的捡拾结果（图5）可以看出，长短时平均方法捡拾的初至波到时往往滞后于实际初动点。对于地震动初动幅值较大的地震，长短时平均方法捡拾误差较小；对于地震动初动幅值较小的地震，长短时平均方法捡拾结果误差稍大。其他3种方法 P波到时捡拾效果较好，但都需要根据所研究震相的特点，合理设置模型参数，并且基于小波函数的AIC值，还要通过对比不同尺度的误差来选择最优解。
图 5 6次地震记录垂直向4种方法捡拾误差对比图
-
介绍了长短时平均（STA/LTA）方法、固定时窗AIC法、长短时平均（STA/LTA）与固定时窗AIC结合法以及基于小波函数的W-AIC法等4种P波自动捡拾方法，结合所研究地震记录的特点，设置相应方法的合适参数，以6次地震记录为例，运用上述4种方法捡拾初至波到时并分析对比捡拾误差，得到以下结论：
1）相较于AIC准则法，长短时平均（STA/LTA）与固定时窗AIC结合法、基于小波系数的W-AIC法、长短时平均值之比（STA/LTA）法3种方法算法简单、运算速度较快，但对背景噪声大或初动幅值较小的地震捡拾误差大；
2）固定时窗AIC法，需设置固定窗长和每次移动步长，计算多自变量的复杂拟合，对地震动记录的捡拾精度较高，但对于初动幅值较小的P波捡拾效果不佳；
3）基于小波函数的W-AIC法，通过不同尺度下的小波系数反映地震信号在初至点的奇异性，但噪声在大的尺度下会逐渐消失，误差的大小随尺度选取而变化，因此需要设置捡拾的限制条件，对比不同尺度下的误差大小择优捡拾，捡拾精度高，但运算过程需要更多时间；
4）长短时平均值之比（STA/LTA）和固定时窗AIC结合法，根据不同地区的地震波特点，合理设置计算参数，一方面可以避免长短时平均（STA/LTA）法滞后的误捡拾，另外一方面解决了固定时窗AIC法对初动幅值较小时捡拾效果不佳的问题，是最为理想的捡拾方法。
Automatic Identification of Direct P-wave First Motion in Earthquake Early Warning
-
摘要: 通过长短时平均（STA/LTA）方法、固定时窗AIC法、长短时平均（STA/LTA）与固定时窗AIC结合法、基于小波系数的W-AIC法等4种方法研究地震初至波P波的自动捡拾，结合所研究的地震记录的特点选择合适的计算参数，通过对比误差发现：长短时平均（STA/LTA）和固定时窗AIC结合法，既避免了长短时平均（STA/LTA）法的将背景噪声误捡拾为初至波的误捡拾，又解决了固定时窗AIC法对初动幅值较小捡拾效果不佳的问题，是较为理想的初至波捡拾方法。Abstract: Through the combination of the STA/LTA method, the fixed time window AIC method, the STA/LTA method and the AIC criterion, and the W-AIC method, based on the wavelet coefficient, the automatically identification of direct P-wave first motion was studied. Based on the characteristics of the seismic records studied, the appropriate calculation parameters were selected and the error was compared. The results show that: the combination of STA/LTA and AIC criterion, on the one hand, can avoid the error identification of STA/LTA method, and on the other hand, it can solve the problem that the combination of AIC has poor effect on the identification of small initial amplitude, so as to achieve more accurate seismic phase identification.
-
表 1 运用到的地震记录统计表
地震记录 震中距/km 震级M 发震时刻 震中位置 1 259.6 4.9 2019-12-26T18:36:34 湖北应城 2 216.7 4.5 2018-10-11T15:06:30 湖北秭归 3 95.6 4.0 2014-05-27T21:57:09 湖北房县 4 135.9 4.2 2008-03-24T23:24:37 湖北竹山 5 30.5 3.6 2019-11-30T04:28:48 河南淅川 6 32.9 4.5 2018-02-09T19:00:37 河南淅川 
下载: 导出CSV
表 2 应用STA/LTA震相自动捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.47 0.13 2 36.18 36.43 0.25 3 15.95 16.16 0.21 4 22.71 22.89 0.18 5 5.41 5.63 0.22 6 5.76 6.03 0.27
下载: 导出CSV
表 3 应用固定时窗AIC震相自动捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.45 0.11 2 36.18 36.12 −0.06 3 15.95 16.11 0.16 4 22.71 22.82 0.11 5 5.41 5.74 0.31 6 5.76 5.89 0.13
下载: 导出CSV
表 4 结合应用STA/LTA与AIC震相捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.43 0.09 2 36.18 36.12 −0.06 3 15.95 16.07 0.12 4 22.71 22.81 0.10 5 5.41 5.61 0.20 6 5.76 5.88 0.12
下载: 导出CSV
表 5 应用基于小波变换的AIC值震相捡拾误差的统计表
地震记录 人机交互到时/s 自动捡拾到时/s 捡拾误差/s 1 40.34 40.45 0.11 2 36.18 36.09 −0.09 3 15.95 16.06 0.11 4 22.71 22.80 0.09 5 5.41 5.56 0.15 6 5.76 5.81 0.05
下载: 导出CSV
-
[1] 马亮. 用于地震预警的单台定位技术研究[D]. 哈尔滨: 中国地震局工程力学研究所, 2013. [2] 江汶乡, 于海英, 赵晓芬, 等. 用于地震预警系统的单台站防误触发算法研究[J]. 自然灾害学报, 2015, 24(2): 23-31. [3] 王庆民, 刘希强, 沈得秀, 等. 地震预警中的动态地震定位方法研究及应用[J]. 地震研究, 2013, 36(4): 485-489. doi: 10.3969/j.issn.1000-0666.2013.04.012 [4] 刘希强, 孙亚强, 赵冰, 等. 网内地震震中近实时动态定位方法研究[C]//中国地球物理学会第二十七届年会论文集. 长沙: 中国地球物理学会, 2011: 483. [5] Nakayachi K, Becker J S, Potter S H, et al. Residents' reactions to earthquake early warnings in Japan[J]. Risk Analysis, 2019, 39(8): 1723-1740. doi: 10.1111/risa.13306 [6] Sheen D H. A robust maximum-likelihood earthquake location method for early warning[J]. Bulletin of the Seismological Society of America, 2015, 105(3): 1301-1313. doi: 10.1785/0120140188 [7] 周彦文. 基于单台P波记录的早期地震预警方法研究[D]. 兰州: 中国地震局兰州地震研究所, 2008. [8] Minson S E, Baltay A S, Cochran E S, et al. The limits of earthquake early warning accuracy and best alerting strategy[J]. Scientific Reports, 2019, 9(1): 2478. doi: 10.1038/s41598-019-39384-y [9] 尹得余. P波震相自动捡拾与地震预警震级实时测定技术研究[D]. 哈尔滨: 中国地震局工程力学研究所, 2012. [10] 毛燕, 崔建文, 郑定昌, 等. 地震记录的P波自动捡拾[J]. 地震研究, 2011, 34(1): 47-51. doi: 10.3969/j.issn.1000-0666.2011.01.008 [11] 姚开一, 李英玉. 基于神经网络的地震震相自动拾取方法[J]. 电子设计工程, 2018, 26(22): 1-5. doi: 10.3969/j.issn.1674-6236.2018.22.001 [12] 郭铁龙, 张雪梅, 邹立晔. STA/LTA-AIC算法对地震P波震相拾取稳定性影响[J]. 地震地磁观测与研究, 2017, 38(3): 13-17. doi: 10.3969/j.issn.1003-3246.2017.03.003 [13] 王喜珍. 小波变换在地震数据压缩和震相到时拾取中的应用研究[D]. 北京: 中国地震局地球物理研究所, 2004. [14] 刘希强, 周彦文, 曲均浩, 等. 应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别[J]. 地震学报, 2009, 31(3): 260-271. doi: 10.3321/j.issn:0253-3782.2009.03.003 [15] Kodera Y. Real‐time detection of rupture development: earthquake early warning using P waves from growing ruptures[J]. Geophysical Research Letters, 2018, 45(1): 156-165. doi: 10.1002/2017GL076118 [16] 赵大鹏, 刘希强, 刘尧兴, 等. 高阶统计量及AIC方法在区域地震事件和直达P波初动识别中的应用[J]. 地震地磁观测与研究, 2013, 34(5): 61-69. doi: 10.3969/j.issn.1003-3246.2013.05/06.010 [17] Allen R M, Melgar D. Earthquake early warning: advances, scientific challenges, and societal needs[J]. Annual Review of Earth and Planetary Sciences, 2019, 47: 361-388. doi: 10.1146/annurev-earth-053018-060457 [18] 时伟, 林春华, 王维红, 等. 双约束变换时窗统计能量比地震波初至拾取方法[J]. 物探与化探, 2019, 43(5): 1064-1073. doi: 10.11720/wtyht.2019.0065 [19] Mărmureanu A, Ionescu C, Cioflan C O. Advanced real-time acquisition of the Vrancea earthquake early warning system[J]. Soil Dynamics and Earthquake Engineering, 2011, 31(2): 163-169. doi: 10.1016/j.soildyn.2010.10.002 -
点击查看大图
图(5) / 表 (5)
计量
- 文章访问数: 2492
- HTML全文浏览量: 1735
- PDF下载量: 30
