Relocation and focal mechanism inversion for the Changning，Sichuan，earthquake on 7 June 2019
-
摘要: 采用双差重定位和W震相波形反演方法分析 “地震编目系统” 和中国地震台网中心提供的地震观测报告及区域地震波形数据，对2019年四川长宁地震序列进行了重定位，反演获取了M＞4.5地震的震源机制解。地震序列重定位结果显示，长宁地震序列沿NW优势方向呈条带状分布，集中分布于5—10 km深度范围，且发震断层面呈高倾角。震源机制反演结果表明，2019年6月17日四川长宁MS6.0主震的两个可能发震断层面参数分别为：节面Ⅰ走向12°，倾角50°，滑动角139°；节面Ⅱ走向131°，倾角59°，滑动角48°，最优矩心深度为7.5 km，矩震级MW5.74。此外几个M＞4.5余震的震源机制也基本与主震类似，均为以逆断为主外加少量走滑的地震破裂事件。综合分析长宁地震序列的重定位、震源机制反演结果以及震中和附近区域的地质构造背景信息推断，本次长宁主震的发震破裂面呈NW−SE走向，发震断层为长宁—双河背斜东北翼发育的逆冲断层。Abstract: In this study, we apply the double-difference relocation technique and W-phase method to relocate the 2019 Changning, Sichuan, earthquake sequence, and to determine the focal mechanisms of the M＞4.5 aftershocks using the seismic data provided by the Earthquake Catalog System and China Earthquake Networks Center. The relocated epicentral locations of Changning earthquake sequence show a NW-striking dominant narrow belt in the predominant epicentral depth range of 5−10 km with a near vertical seismogenic fault plane. The focal parameters for the main shock of 17 June, 2019 Sichuan Changning earthquake are: strike 12°, dip 50°, strike 139° for the nodal plane Ⅰ , and strike 131°, dip 59°, strike 48° for the nodal plane Ⅱ , focal depth is 7.5 km and moment magnitude is MW5.74. Moreover, the M＞4.5 aftershocks show a similar focal mechanism to the main shock. Combined analyses of the relocated epicentral location, focal mechanisms and tectonic settings suggest that the seismogenic fault of 17 June 2019 Changning main shock could be related to the thrust fault developed in north-eastern wing of NW-SE-striking Changning-Shuanhe anticline.
-
引言
地震震相信息是地震学研究的重要基础数据，广泛应用于地震定位、震源机制分析、走时层析成像等研究，震相自动识别也是地震预警技术的重要基础. 随着地震观测技术的不断发展，利用密集地震台网和大规模流动地震台阵开展地震监测研究受到广泛重视，与此同时观测产生的地震波形数据急剧增加. 为了获取观测区内的微小地震，需要对大量的连续地震观测数据进行仔细分析，而传统方法是依靠人工一帧一帧地浏览地震波形识别地震，并逐一分析读取震相，工作量巨大，整个过程单调枯燥且效率极低，因此地震震相自动拾取和地震自动检测技术的发展，一直备受关注.
震相拾取方法众多，每种方法均有其优势和局限性，而有效性和计算效率是衡量震相拾取方法优劣的两个重要指标，而且这两个指标会受到诸如数据质量、算法简洁性等因素的影响（Chen，Holland，2016）. 应用于实时监测的震相自动识别方法还要求在保证一定效果的情况下，注重算法的计算效率. 与之区别的常规测震分析则需注重有效性和精度，常用于微震识别和地震定位等分析.
震相到时拾取算法可分为在时域、频域中进行两类，包括能量瞬时方法（Allen，1978，1982；Withers et al，1998 ；高淑芳等，2008）、偏振分析（Flinn，1965；Vidale，1986；Bai，2000；Jurkevics，2002）、自回归赤池信息准则（Akaike information criterion，简写为AIC）方法（Akaike，1974；Kitagawa，Akaike，1978；Maeda，1985； Leonard，Kennett，1999；Leonarda，2000；Zhang et al，2003 ；赵大鹏等，2012）、人工神经元（Wang，Teng，1997；Gentili，Michelini，2006；Dai，MacBeth，2007）、高阶统计（Saragiotis et al，2002 ；Panagiotakis et al，2008 ； Baillard et al，2014 ；Hibert et al，2014 ）、小波变换（Anant，Dowla，1997；刘希强等，1998；刘代志等，2005；Akansu et al，2010 ；Vassallo et al，2012 ；Bogiatzis，Ishii，2015）、奇异值分解（Rosenberger，2010；Kurzon et al，2014 ）等方法. 刘希强等（2009）的研究表明，根据信号与噪声的单一差异特征而提出的方法均具有其应用的局限性，特别是低噪声水平条件下震相自动识别的准确率仍不理想，而综合多个信号与噪声的差异特征有利于提高震相自动识别的准确性.
在震相拾取时，利用多台信息能够降低噪声的干扰，提高震相识别的可靠性. 马强等（2013）将德洛内（Delaunay）三角用于地震预警的P波震相拾取，其基本原理是利用短长时窗均值之比（short term average to long term average，简写为STA/LTA）方法得到初步到时，再用AIC方法得到P波的精确到时，通过德洛内三角判别该到时是否为干扰信号. 该方法现已应用于福建省地震预警试验系统中，在剔除干扰信号方面具有良好的效果，但并未将各台站触发到时归类到相同的地震事件中.
针对以上问题，本文拟对拾取的到时信息，利用台间距与台站相邻关系作为约束来减少噪声干扰，并在快速地震定位的基础上，根据台站理论到时选取合理的时间窗，采用STA/LTA方法、偏振分析与基于方差的AIC （var-AIC）的组合方法进行更精细的分析，以期通过该方法获取更可靠的P波震相到时信息.
1. 方法
目前震相拾取主要分为3步：① 数据预处理. 主要对波形进行滤波或时频转换，包括小波变换、奇异值分解、傅里叶变换等；② 地震信号检测. 以STA/LTA方法为主，但由于STA/LTA方法具有一定的延迟性，因此只能得到地震信号所在的大致到时；③ 拾取震相到时. 主要通过AIC方法、偏振分析方法以及人工神经元等方法拾取P波或S波到时.
以上步骤首先基于STA/LTA的检测到时结果，但是地震观测数据经常会受到台站周边人文环境、气候等因素的干扰，例如车辆、工厂、动物活动以及大风、雷电等产生的干扰，由此检测到虚假震相；但是由于这些干扰的传播距离较小，往往仅能被邻近台站观测到，利用相同地震事件中震相到时具有的同源性，可滤除大部分这类噪声干扰. 利用震源位置信息获得预测的P波到时选取窗口，再使用偏振分析和AIC函数等方法进行更精细的震相到时拾取. 同时，对于一些STA/LTA方法获得的信噪比较低的台站，在P波理论到时窗口内进一步分析，从而判定是否能获得有效的到时信息.
1.1 地震事件检测
首先，本文使用STA/LTA方法（Allen，1978，1982）快速获得低精度的地震波初至到时. 计算短时窗N2与长时窗N1内的特征函数CF的平均值比值STA（t）/LTA（t），其中：
${{CF}}(t) = {Z^2}(t),$
(1) ${{LTA}}(i) = {{LTA}}(i - 1)+ \frac{{{{CF}}(i) - {{CF}}(i - {N_1})}}{{{N_1}}},$
(2) ${{STA}}(i) = {{STA}}(i - 1) + \frac{{{{CF}}(i) - {{CF}}(i - {N_2})}}{{{N_2}}},$
(3) 式中，Z表示地震记录的垂直分量. 当STA（t）/LTA（t）超过阈值C1时，该时刻tP记为地震波初至到时；当STA（t）/LTA（tP）低于阈值C2时，该时刻记为信号恢复时刻te，此时获得的tP到时相较P波实际到时有一定滞后. Maeda（1985）提出基于方差的AIC方法，利用该方法来计算AIC函数的最小值，能得到更精确的P波到时，其中AIC函数为
${ AIC}(t) = t\lg [{\mathop{ var}} (D[1,t])]+ (N - t - 1)\lg [{\mathop{ var}} (D[t + 1,N])],$
(4) 式中，N为信号长度，D为截取信号. var-AIC方法需要预先设置好窗长才能计算，在STA/LTA初步拾取后，本文采用一个长度为N3的小窗对低精度的初至到时进行适当校正.
将各台站记录到的触发到时按时间先后顺序排列，并将震相到时和台站信息按照地震事件进行分组归类. 与台站附近的局部噪声干扰不同，一个地震事件若能被多个地震台站观测到，则相邻台站记录的初至震相到时应该具有较好的相关性. 根据地震波传播的特点可知，相邻台站间震相的最大到时差应小于台间距除以波速并考虑一定的误差量，而台站附近的局部噪声干扰具有随机性，相关性较弱. 根据这一特征，通过引入台网德洛内三角划分，并考虑台间距的影响，即可将震相到时差大于约束条件的台间距排除在外. 对于相继触发的两个台站1和2，根据以下3个条件判断二者触发的到时信号是否可能属于同一事件：① 台站2与台站1记录的到时差是否小于台间距除以波速vP；② 台站2与台站1的台间距是否小于台间距阈值C3；③ 台站2与台站1是否在同一德洛内三角内. 满足条件①且满足条件②或③的两台站所触发的到时认为在同一事件内，其它情况则排除在同一事件外.
通过以上步骤，将时空上具有相关性的相邻到时记录暂时归类为同一事件. 但考虑到后续触发的台站，其震中距不断增大，同一事件触发的台站虽然在时间上相邻，但是其既可能在空间上不相邻，也可能不在台间距阈值内，此时可能会造成同一事件记录的重复拾取. 为避免这种情况，将各个相邻并且互相有所关联的事件合并：两个相邻的事件A和B，取事件B触发的前5个台站与事件A触发的最后5个台站，分别计算前一组中台站与后一组每个台站的距离和到时差（事件A和B不足5台时均计算全部），若有两个台站满足上述3个条件，则将事件A与事件B合并为同一个事件. 最后提取其中触发台站数大于一定数值的事件记为地震事件，其它事件则作为干扰噪声予以剔除.
滤除无相关性的噪声并获得地震事件对应的各台站到时后，利用偏振特性计算各台站的线性偏振度P，从而获取S波震相到时. 首先选取一个长度为N的滑动视窗，计算三分向协方差矩阵的3个特征值λ1，λ2和λ3 （λ1>λ2>λ3），则可得到线性偏振度P为
$P(t) = 1 - \frac{{{\lambda _2} + {\lambda _3}}}{{{\lambda _1}}},$
(5) This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
$D{_{{ b}t}} = \frac{{{\lambda _{t + }}}}{{{\lambda _{t - }}}};$
(6) 同样取一短窗长Ls，计算t短窗的极化度Dst. 拾取震相的特征函数Pol为
${ Pol}(t) = D{_{{ b}t}}D{_{{ s}t}}.$
(7) 当Pol取最大值时，t点即为震相初至. 由于P波和S波均具有偏振特性，实际波形中P波的尾波会影响到S波的偏振性，因此若Pol函数的值低于阈值C4，则认为S波到时的计算结果不可靠，不记录该S波的到时.
1.2 地震定位判断
首先计算某事件所有触发台站的平均坐标，取平均坐标为中心的正方形区域，按照经纬度和深度进行三维网格化，并确定网格步长.
This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
${t_0} = \sum\limits_{i = 1}^n {\frac{{t_2^i - t_1^i}}{n}} ,$
(8) ${t_3} = \sqrt {\sum\limits_{i = 1}^n {\frac{{{{(t_2^i - t_1^i - {t_0})}^2}}}{n}} } ,$
(9) 式中n为事件所包含的震相到时数. 取所有网格的标准差t3的最小值，其所处位置为震源位置的可能性最大，此时各台站的残差为
$t_4^i = \left| {t_2^i - t_1^i - {t_0}} \right|.$
(10) This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
1.3 震相到时拾取的组合函数方法
在时窗内拾取震相到时之前，须判断该时窗内是否存在比较清晰的震相. 为此，本文设计了一个简单的震相存在判别方法：首先查看该台站是否为地震检测过程中触发的台站，若是，则认为该台站记录到震相并进行下一步计算，否则计算时窗内的最大偏振度以及窗内长短时触发位置处的信噪比，若二者之一大于各自的阈值C6或C7，则判定该时窗内有震相.
考虑到AIC函数和STA/LTA拾取到时的结果随信噪比的变化较为强烈，当信噪比较低时，AIC函数或STA/LTA可能表现不稳定，从而使得二者拾取到时的差值较大；而偏振分析方法虽然受到信噪比的影响程度相对较小，但其在高信噪比时拾取精度不如前两者高，因此我们设计了一种新的组合方法，提高震相拾取的可靠性，即利用AIC函数与STA/LTA两方法拾取到时的差值作为信噪比高低的判别依据，当差值较小即信噪比较高时，选择较准确的AIC函数的极小值作为可靠结果，信噪比较低时设计新的拾取函数
${ Pick}\left( {{i}} \right) = \left[ {1 - {{AIC}}\left( {{i}} \right)} \right] { Pol}\left( {{i}} \right),$
(11) 该函数的最大值对应P波到时. Pick函数利用归一化后的AIC函数压制了P波到达后其它信号对Pol函数造成的影响. 通过上述方法在时窗内进一步分析P波、S波到时，主要步骤为：① 先在预测P波、S波到时前后取窗长N4，根据震相存在判别方法判断该时窗内是否具有震相，若有则在该时窗内同时使用STA/LTA方法和var-AIC方法检测到时；② 若AIC函数与STA/LTA检测到的到时差小于C8，则将该到时记为P波到时；若大于C8或STA/LTA未检测到任何到时信号，则对时窗内的波形数据使用偏振分析方法并对AIC函数归一化，取拾取函数Pick的最大值为P波到时；③ 以P波到时与AK135计算得到的S波到时的中间值为时窗的起始点，S波到时向后N4范围时窗结束，在该时窗内计算Pol函数的最大值，若最大值大于C4，则记为S波到时.
2. 计算实例
2013年1月，房立华等（2014）在安宁河断裂带附近布设了由32个流动地震台站组成的密集流动地震台阵，即西昌地震观测台阵，开展针对断层带的地震观测研究. 台阵的台间距介于17.7—60.5 km，平均台间距为27.8 km. 地震观测仪的频带为50 Hz—1 200 s，采用GPS授时，观测使用的采样率为100 Hz.
本文以西昌地震台阵记录到的2013年1月28日21时4分ML1.2地震为例，此时共有28个台站正常运行. 首先使用德洛内三角剖分算法将台站划分为若干三角区（图1），接着使用本文提出的方法对该地震事件进行震相到时自动拾取试验，采用的相关参数列于表1. 首先对波形数据进行1—30 Hz的带通滤波处理，然后计算STA/LTA，并拾取震相到时. 在STA/LTA的计算过程中，长时窗N1设为30 s，短时窗N2设为0.05 s，触发阈值C1设为10，考虑到地震波到来后短时间内无法恢复至触发前的噪声水平，因此STA/LTA的恢复阈值C2设为1.5. 在利用STA/LTA方法拾取震相到时的基础上，进一步通过var-AIC和偏振分析方法拾取P波和S波到时. 图2给出了针对XC03台站的地震波形记录获得的分析结果. 由STA/LTA，var-AIC，偏振分析方法以及组合函数方法计算得到的P波初至结果分别为7.52，7.49，7.48和7.48 s，人机交互拾取结果为7.44 s. 与人工方法拾取的到时相比，var-AIC，偏振分析方法和组合函数方法拾取的到时分别相差0.05 s和0.04 s，STA/LTA方法的拾取到时则相差0.08 s.
![]()
图 1 西昌地震台阵的德洛内三角剖分图正方形为快速定位时网格搜索区域，星形为地震实例发生位置，三角形为台站位置Figure 1. Delaunay triangulation of Xi’chang seismic arrayThe square represents the search area，the star stands for the hypocentral location，and triangles for stations表 1 本文方法采用的相关参数Table 1. Related parameters used in the method proposed in this paper参数名 数值 单位 意义 C1 10 长短时方法获得到时的触发阈值 C2 1.5 长短时方法信号恢复时的触发阈值 C3 30 km 多台检测事件台间距的上限 C4 20 偏振分析Pol函数的触发阈值 C5 5 s 定位时的残差阈值 C6 10 震相存在判别的信噪比阈值 C7 1 000 震相存在判别的Pol函数阈值 C8 0.5 s AIC与STA/LTA结果接近的到时差范围 N1 30 s 长时平均值计算的滑动窗长 N2 0.5 s 短时平均值计算的滑动窗长 N3 1 s 为纠正STA/LTA拾取到时的AIC函数计算窗长 N4 4 s 代入AK135后拾取P波和S波到时的窗长 Lb 0.5 s 通过Pol函数计算S波到时设置的长窗长 Ls 0.2 s 通过Pol函数计算S波到时设置的短窗长 vP 5.8 km/s 检测地震事件时所取的P波速度 ![]()
图 2 西昌地震台阵记录的2013年1月28日21时4分ML1.2地震（a） 垂直分量波形；（b） 短长时窗均值之比；（c） 归一化的var-AIC函数；（d） Pol函数拾取P波；（e） Pol函数拾取S波；（f） 改良后的组合函数Figure 2. Event of ML1.2 recorded by Xi’chang seimic array at 21:04 on January 28，2013（a） Vertical component of waveform；（b） STA/LTA；（c） Normalized var-AIC；（d） P-wave picking by Pol-function；（e） S-wave picking by Pol-function；（f） Improved compound function图3给出了2013年1月28日21时4分ML1.2地震在15个观测台站的垂直向地震记录. 由于噪声干扰等因素的影响，使用STA/LTA方法与var-AIC方法初步拾取的震相到时出现了许多与地震信号无关的虚假震相到时（图3a）. 利用地震信号在多台记录中满足同源性的特点，采用德洛内三角和台间距与传播时间关系的约束，在很大程度上消除了台站附近局部噪声干扰的影响，很好地保留了基本满足同源传播特征的地震震相到时（图3b），成功地检测出地震事件. 值得注意的是，本文设置判断到时相关性的台间距阈值C3略大于平均台间距，主要是考虑到台间距阈值C3的主要作用是约束远处台站干扰信号的影响，而平均台间距能够作为台站远近的一项评价指标.
在初步完成地震事件检测后，以同一地震事件所有触发台站的平均坐标为中心，南北向、东西向的长度均取为1°，深度取0—60 km，以0.05°×0.05°×5 km的三维网格进行快速地震定位. 图1中正方形区域为该地震事件的搜索区域，图4给出了该事件的定位结果. 定位结果显示，与人工分析的正式定位结果相比，震中仅相差5.39 km. 在地震定位过程中，对于一些走时残差过大的震相数据进行适当剔除，可以进一步消除各种干扰的影响.
![]()
图 3 地震观测记录及拾取的震相到时图（a） STA/LTA方法拾取的各种可能的P波震相到时；（b） 利用本文地震检测方法分析处理后获得的P波震相到时；（c） 用AK135模型计算的理论到时，其中紫线为P波预测到时，黄线为S波预测到时，蓝线为地震检测方法得到的P波震相到时，红线和绿线分别为取窗拾取的P波到时和S波到时，其中剔除了Pol函数小于C4的S波到时Figure 3. Waveform records and arrival picking in different steps（a） Probable P-waves picking by STA/LTA；（b） P-waves processing by the earthquake detection method in this paper；（c） Arrival time calculated by AK135 model，where purple lines denote arrival time of P-wave，yellow lines denote arrival time of S-wave，blue lines denote arrival time of P-wave processing by earthquake detection method，red and green lines denote P-wave and S-wave picking by setting windows in waveforms，S-wave of which Pol-function below C4 is deleted根据地震定位结果，本文采用AK135模型计算各个台站的P波和S波预测到时（如图3c中的紫色和黄色竖线所示），并在预测到时附近设定适当宽度的时间窗，进一步拾取更可靠的震相到时. 由于采用的速度模型与真实模型存在差异、震相到时拾取存在误差以及定位结果本身也存在误差等多种因素的影响，理论到时与实际震相的到时之间存在一定的差异. 在实际资料的处理中，以P波预测到时为中心前后取4 s时窗，拾取P波到时. 取P波到时与S波预测到时的中间值与S波预测到时后4 s作为时窗计算Pol函数的最大值，即为S波到时. 信噪比较高时（图2），AIC函数与STA/LTA之差很小，此时较准确的AIC函数的极小值可作为可靠结果；若信噪较低或受到干扰信号影响，AIC函数或短长时之比可能表现得并不稳定，二者差值也会较大，此时可通过本文拾取函数拾取P波到时，如图5所示. 图3c为AK135模型计算的预测到时与取窗分析到时对比图，可见相较于预测到时，取窗分析的结果已有很大提高.
![]()
图 4 实际震源位置与算法计算结果的对比图中数字表示事件触发顺序Figure 4. Comparison of hypocentral location with calculation resultThe numbers represent the triggered order of events![]()
图 5 各震相拾取方法精度对比（a） 波形计算窗口；（b） 垂直分量波形；（c） 短长时比值；（d） var-AIC函数；（e） Pol函数；（f） 改良后的组合函数Figure 5. Accuracy comparison of different phase picking methods（a） Calculation window；（b） Vertical component of waveform；（c） STA/LTA；（d） var-AIC；（e） Pol-function；（f） Improved compound function3. 实际资料处理试验与讨论
为验证本文方法的有效性，本文选用2013年西昌地震观测台阵1—8月记录到的304个地震事件，震级分布于ML−0.5—2.2，震源深度最大为37 km. 这些地震是通过人工一帧一帧地浏览观测波形仔细分析后得到的，对许多微弱的地震信号均进行了分析. 根据人工分析获得的地震目录，将每一个地震截取为180 s的单事件波形数据，利用本文方法进行震相自动拾取，并将结果与人工拾取结果进行对比，以确认其有效性.
首先采用台间距阈值约束方法、德洛内三角判别的方法和本文方法对304个地震事件进行检测，并以人工拾取结果作为参照. 表2列出了3种检测方法的结果对比，可以看出，使用本文提出的改进方法检测，遗漏的地震事件明显少于其它两种方法，误拾取的地震事件数与德洛内三角判别检测方法相近.
由于待检测的地震大多震级较小，对于平均台间距为28 km的流动地震台阵而言，一些微小地震事件的地震波形往往只能被少量台站清晰地记录到. 图6给出了检测到的地震数与地震震级及设定的震相触发最少台站数的关系，其中触发震相要满足改进方法的相关约束条件. 可以看出：设定的最少触发台站数越小，能检测到的地震数越大；当ML>0.9时，采用最少3台触发能检测到所有的地震；当ML<0.6时，设定的最小触发台站数对遗漏的地震事件数具有明显的影响；对于ML>0.6地震，设定的最少触发台站数为3，4和5时，地震的正确识别率分别为95%，88%和81%；对于ML>0.3地震，相应的正确识别率则分别为90%，77%和69%. 为了检验震源深度可能对地震检测效果的影响，本文将不同深度的地震检测遗漏情况进行统计，结果如图7所示，可以看出不同深度的遗漏事件占比并无显著的相关性.
检测完地震事件后，为尽可能地保证定位效果，选取触发4台或4台以上的地震事件用于定位，共包含186个事件. 其中，有2个事件在定位过程中因剔除较大残差的到时而导致触发台站少于4个而放弃，最终共有184个事件用于计算. 地震定位采用AK135模型作为理论走时计算的速度模型. 由于本文在通过定位结果获取的预测到时前后4 s的时窗内分析实际到时，因此定位误差小于20 km足以满足进一步分析的要求，将定位结果与人工定位结果作比较，共有80%的地震事件定位误差小于20 km，89%的ML>1.0地震事件的定位误差小于20 km.
根据地震定位结果，采用AK135速度模型计算各台站的预测到时，并以预测到时为中心设定时窗，然后根据震相存在判别方法判断各台站是否存在震相，若存在则拾取该台站时窗内的震相. 为验证震相存在判别方法的准确性，用该方法判别184个地震事件中存在震相的台站数并与地震检测过程以及人工拾取结果进行对比分析，结果如图8所示，可以看出：该方法判别有震相记录的台站与人工拾取到震相记录的台站高度重合，人工拾取震相的85%能够被检测到，而之前不考虑理论走时的震相检测方法仅有77%的震相能被检测出来；对于ML>0.6地震，两种方法的检测能力分别是人工拾取震相的87%和79%；对于ML>1.5地震，分别能检测到89%和73%的震相. 基于理论到时的震相识别，更容易检测到一些人工识别遗漏的震相.图9给出了新方法震相检测情况与人工震相识别结果对比的两个实例. 可以看出，由于受多种因素的影响，人工识别既可能遗漏掉一些微弱的震相，也可能遗漏掉个别清晰的震相，而这些震相有可能被新的自动识别方法很好地识别出来（图9a）；自动识别方法也可能受短时强噪声干扰的影响错误地将噪声干扰识别为震相（图9b），需要通过其它方法尽可能排除这些干扰.
表 2 3种地震检测方法对304个地震事件的检测效果对比Table 2. Comparison of detection results for 304 events by three event detection methods台站数 拾取正确事件数 拾取错误事件数 拾取重复事件数 遗漏事件数 拾取事件数 本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别3 239 186 214 19 3 21 32 19 29 65 118 90 290 208 264 4 186 135 167 8 0 7 12 5 14 118 169 137 206 140 188 5 149 107 135 1 0 2 10 3 10 155 197 169 160 110 147 ![]()
图 6 触发台站数不同时所拾取事件数与震级的关系Figure 6. Relationship between the amount of events with magnitude with different amount of triggered station![]()
图 7 遗漏事件占比与深度关系图Figure 7. Relationship between the proportion of omitted events with focal depth![]()
图 9 本文方法与人工拾取结果对比图横坐标为相对预测到时的时间，纵坐标已对齐预测到时。红线为本文方法拾取到的P波震相，蓝线为人工拾取P波震相Figure 9. Comparison of picking results by manual and the method proposed in this paperHorizontal ordinate is relative time. Vertical ordinate has been lined up with forecasted arrival time. Red lines represent P-waves picking by the method in this paper，and blue lines represent P-waves picking by manual基于地震定位后的震相预测到时设定的时窗，通过震相检测判别后，分别采用var-AIC方法、偏振分析方法和本文组合函数方法进行震相自动拾取，并将结果与人工的拾取结果进行对比分析，拾取效果如图10所示，可见：与人工拾取的震相到时相比，误差在0.2 s内的震相数量，本文方法要优于var-AIC方法和偏振分析方法；在误差大于0.2 s且信噪比低于4 dB时，本文方法优于var-AIC方法但稍弱于偏振分析方法.
![]()
图 10 不同方法震相拾取率与信噪比SNR的关系图（a） 拾取误差小于0.1 s；（b） 拾取误差小于0.2 s；（c） 拾取误差小于0.3 s；（d） 拾取误差小于0.4 sFigure 10. Relationship between phase picking rate with SNR by different methods（a） Deviation below 0.1 s；（b） Deviation below 0.2 s；（c） Deviation below 0.3 s；（d） Deviation below 0.4 s4. 讨论与结论
本文基于var-AIC和偏振分析方法的特点，提出了一种改进的震相拾取方法. 实际资料处理试验揭示，该方法总体上明显优于原有的两种方法，只是在信噪比低于4 dB左右时稍弱于偏振分析方法.
采用STA/LTA方法、var-AIC以及偏振分析方法拾取可能的震相到时，利用地震具有同源性的特征，通过台站的德洛内三角形划分、与台间距有关的最大走时约束以及地震定位约束等，可以较好地排除噪声干扰的影响，提高地震事件检测的可靠性. 西昌台阵2013年1—8月的实际资料检验表明，当作为地震事件判别的最小触发台站数选为4时，本文提出的方法对于ML>0.5地震的识别率为82%，ML>0.8地震的识别率为89%，ML>1.0地震的识别率为96%，比单独使用台间距阈值约束方法或使用德洛内三角判别方法得到的地震识别率分别提高21%和3%—4%.
以地震定位结果为基础，在计算的震相理论到时前后设定合理的震相搜索时间窗，可以从观测资料中拾取更多的震相，提高震相到时拾取的可靠性. 本研究仅采用简单的AK135模型进行地震定位，对于观测震相数量较少、台站分布较差的地震，定位精度可能不高，会在一定程度上影响最终震相拾取的可靠性. 在今后的研究中，可以结合观测区域的实际情况，建立合适的区域速度模型或三维速度模型开展地震定位和地震震相到时计算，有可能获得更加可靠的地震检测和震相拾取结果.
-
![]()
图 1 2019年6月17日四川长宁地震序列震中及附近区域地质构造（邓起东等，2002；唐永等，2018）
Figure 1. Regional tectonic settings and location of epicenters of 17 June，2019 Changning， Sichuan，earthquake sequence （Deng et al，2002；Tang et al，2018）
![]()
图 2 （a） 地震序列重定位中所使用的地震台站及长宁MS6.0主震震中（坐标来自于下文的矩张量反演），图中活动断层引自邓起东等 （2002）；（b） P波和S波观测走时曲线
Figure 2. （a） Distribution of seismic stations used to relocation and location of the MS6.0 Changning main shock （coordinate from following moment tensor inversion），where dash red lines indicate active faults （Deng et al，2002）； （b） Observed P and S wave travel time curves
![]()
图 3 经重新定位后四川长宁地震序列的平面（a）及深度剖面图（b，c）
主震坐标来自于下文的矩张量反演
Figure 3. Map view （a） and vertical profiles （b，c） of the relocated Changning，Sichuan，earthquake sequence
The star indicates the main shock，its coordinate comes from following moment tensor inversion
![]()
图 4 （a） 2019年6月17日四川长宁MS6.0地震的矩张量反演结果；（b） 部分参加反演台站的波形拟合及台站分布图，台站名称标示于波形图上方
Figure 4. （a） Centroid moment tensor solution for the MS6.0 Changning，Sichuan，earthquake on 17 June 2019；（b） Comparison of the observed waveforms （black lines） and synthetic ones （red lines） for moment tensor solution，where the station name is marked above each waveform
![]()
图 5 2019年6月17日四川长宁MS6.0地震震源深度及震中位置空间网格搜索结果
（a） 深度搜索与波形拟合均方差变化图；（b） 经纬度搜索与波形拟合均方差分布图，图中黑色十字标示最佳震中位置
Figure 5. Focal depth and grid search results of optimal centroid location for 17 June 2019 Changning earthquake
（a） Variations of waveform fitting root mean square （RMS） with depth；（b） Spatial variation of waveform fitting RMS，where best-fitting location is marked with a black cross
表 1 2019年四川长宁地震序列重定位所采用的速度模型
Table 1 Crustal velocity model used to relocate the 2019 Changning，Sichuan，earthquake sequence
顶层深度/km vP/（km·s−1） vP/vS 0.0 4.58 1.72 4.0 5.65 13.0 6.12 26.0 6.55 32.0 7.32 46.0 8.10 
下载: 导出CSV
表 2 国内外相关机构发布的2019年6月17日四川长宁MS6.0地震的震源机制参数
Table 2 The focal mechanism parameters for the MS6.0 Changning，Sichuan，earthquake on 17 June 2019 released by related institutions
机构 节面Ⅰ 节面Ⅱ 震源
深度
/kmMW 走向/° 倾角/° 滑动角/° 走向/° 倾角/° 滑动角/° 美国地质调查局（USGS） 308 45 40 187 63 127 11.5 5.79 哈佛大学全球矩张量（GCMT） 323 57 65 184 40 123 12.0 5.7 中国地震台网中心（CENC） 350 67 93 162 23 82 16.0 6.0 中国科学院地质与地球物理研究所（IGGCAS） 119 58 9 24 82 147 5.0 5.8 中国地震局地球物理研究所（IGP-CEA） 307 58 40 191 57 141 16.0 6.0 中国地震局地震预测研究所（IEF-CEA） 296 81 −13 28 77 −170 19.0 5.8
下载: 导出CSV
表 3 四川长宁地震序列六次M＞4.5余震的震源参数
Table 3 Focal parameters for six M＞4.5 aftershocks of the Changning earthquake sequence
发震时刻 节面Ⅰ 节面Ⅱ 震源
深度
/kmMW 年−月−日 时:分:秒 走向/° 倾角/° 滑动角/° 走向/° 倾角/° 滑动角/° 2019−06−17 23:36:01 175 51 117 316 46 60 9.5 5.08 2019−06−18 07:34:33 156 42 118 300 53 66 7.5 4.72 2019−06−22 22:29:56 162 44 75 3 47 104 9.5 5.19 2019−06−23 08:28:17 146 59 20 45 72 147 7.5 4.44 2019−07−03 12:26:53 135 51 93 309 38 85 7.5 4.61 2019−07−04 10:17:58 192 43 82 22 46 96 9.5 5.16
下载: 导出CSV
-
国家测震台网数据备份中心. 2017. 国家测震台网地震波形数据[DB/OL]. [2019-07-28]. 中国地震局地球物理研究所. doi: 10.11998/SeisDmc/SN. Data Management Centre of China National Seismic Network. 2017. Waveform data of China National Seismic Network[DB/OL]. [2019-07-28]. Institute of Geophysics, China Earthquake Administration. doi: 10.11998/SeisDmc/SN （in Chinese）.
邓起东,张培震,冉勇康,杨晓平,闵伟,楚全芝. 2002. 中国活动构造基本特征[J]. 中国科学:D辑,<bold>32</bold>(12):1020–1030. Deng Q D,Zhang P Z,Ran Y K,Yang X P,Min W,Chu Q Z. 2003. Basic characteristics of active tectonics of China[J]. <italic>Science in China</italic>:<italic>Series D</italic>,<bold>46</bold>(4):356–372.
郭志,陈立春,李通,高星. 2018. 2017年8月8日四川九寨沟<italic>M</italic>7.0和9日新疆精河<italic>M</italic>6.6地震震源机制解[J]. 地震地质,<bold>40</bold>(6):1294–1304. Guo Z,Chen L C,Li T,Gao X. 2018. Focal mechanism of the August 8,<italic>M</italic>7.0 Sichuan Jiuzhaigou and August 9,<italic>M</italic>6.6 Xinjiang Jinghe earthquakes of 2017[J]. <italic>Seismology and Geology</italic>,<bold>40</bold>(6):1294–1304 (in Chinese).
郝金来,姚振兴. 2012. 小波域反演确定区域地震震源机制[J]. 中国科学:地球科学,<bold>42</bold>(2):191–201. Hao J L,Yao Z Y. 2013. Determination of regional earthquake source parameters in wavelet domain[J]. <italic>Science China Earth Sciences</italic>,<bold>55</bold>:296–305 (in Chinese).
何登发,鲁人齐,黄涵宇,王晓山,姜华,张伟康. 2019. 长宁页岩气开发区地震的构造地质背景[J]. 石油勘探与开发,<bold>46</bold>(5):993–1006. He D F, Lu R Q, Huang H Y, Wang X S, Jiang H, Zhang W K. 2019. Tectonic and geological background of the earthquake hazards in Changning shale gas development zone[J]. Petroleum Exploration and Development, 46（5）: 993−1006 （in Chinese）.
黄浩,付虹,沙成宁,王培玲. 2017. 2016年青海门源<italic>M</italic><sub>S</sub>6.4地震重定位[J]. 地震学报,<bold>39</bold>(2):176–187. doi: 10.11939/jass.2017.02.002 Huang H,Fu H,Sha C N,Wang P L. 2017. Relocation of the 2016 <italic>M</italic><sub>S</sub>6.4 Menyuan,Qinghai earthquake[J]. <italic>Acta Seismologica Sinica</italic>,<bold>39</bold>(2):176–187 (in Chinese). doi: 10.11939/jass.2017.02.002（inChinese）
唐永,周立夫,陈孔全,董晓霞,唐文军. 2018. 川东南构造应力场地质分析及构造变形成因机制讨论[J]. 地质论评,<bold>64</bold>(1):15–28. Tang Y,Zhou L F,Chen K Q,Dong X X,Tang W J. 2019. Analysis of tectonic stress field of southeastern Sichuan and formation mechansim of tectonic deformation[J]. <italic>Geological Review</italic>,<bold>64</bold>(1):15–28 (in Chinese).
韦伟,谢超,周本刚,郭志,尹昕忠,李保华,王萍,董绍鹏,魏费翔,王禹钦. 2018. 西藏米林<italic>M</italic>6.9级地震及其余震序列地震定位[J]. 科学通报,<bold>63</bold>(15):1493–1501. Wei W,Xie C,Zou B G,Guo Z,Yin X Z,Li B H,Wang P,Dong S P,Wei F X,Wang Y Q. 2018. Location of the mainshock and aftershock sequences of the <italic>M</italic>6.9 Mainling earthquake,Tibet[J]. <italic>Chinese Science Bulletin</italic>,<bold>63</bold>(15):1493–1501 (in Chinese). doi: 10.1360/N972017-01286
吴建平,明跃红,王椿镛. 2006. 川滇地区速度结构的区域地震波形反演研究[J]. 地球物理学报,<bold>49</bold>(5):1369–1376. doi: 10.3321/j.issn:0001-5733.2006.05.016 Wu J P,Ming Y H,Wang C Y. 2006. Regional waveform inversion for crustal and upper mantle velocity structure below ChuanDian region[J]. <italic>Chinese Journal of Geophysics</italic>,<bold>49</bold>(5):1369–1376 (in Chinese).
郑秀芬,欧阳飚,张东宁,姚志祥,梁建宏,郑洁. 2009. “国家数字测震台网数据备份中心” 技术系统建设及其对汶川大地震研究的数据支撑[J]. 地球物理学报,<bold>52</bold>(5):1412–1417. doi: 10.3969/j.issn.0001-5733.2009.05.031 Zheng X F,Ouyang B,Zhang D N,Yao Z X,Liang J H,Zheng J. 2009. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake[J]. <italic>Chinese Journal of Geophysics</italic>,<bold>52</bold>(5):1412–1417 (in Chinese).
翟亮,张晓东,王伟君. 2019. 2017年8月9日精河<italic>M</italic><sub>S</sub>6.6地震余震序列精定位及发震构造分析[J]. 地震学报,<bold>41</bold>(3):314–328. doi: 10.11939/jass.20180124 Zhai L,Zhang X D,Wang W J. 2019. Precise location and seismogenic structure analysis of aftershock sequence of Jinghe <italic>M</italic><sub>S</sub>6.6 earthquake on August 9,2017[J]. <italic>Acta Seismologica Sinica</italic>,<bold>41</bold>(3):314–328 (in Chinese).
中国地震局. 2019. 2019年 6月17日22时55分四川宜宾市长宁县发生6.0级地震[EB/OL]. [2019-07-28]. https://www.cea.gov.cn/cea/dzpd/dzzt/5437173/5437174/5477399/index.html. China Earthquake Administration. 2019. An earthquake with M6.0 occurred in Changning county, Yibin City, Sichuan Province at 22:55 on 17 June 2019[EB/OL]. [2019−07−28]. https://www.cea.gov.cn/cea/dzpd/dzzt/5437173/5437174/5477399/index.html （in Chinese）.
Dreger D S,Helmberger D V. 1993. Determination of source parameters at regional distances with three‐component sparse network data[J]. <italic>J Geophys Res</italic>:<italic>Solid Earth</italic>,<bold>98</bold>(B5):8107–8125. doi: 10.1029/93JB00023
Duputel Z,Rivera L,Kanamori H,Hayes G. 2012. W phase source inversion for moderate to large earthquakes (1990−2010)[J]. <italic>Geophys J Int</italic>,<bold>189</bold>(2):1125–1147. doi: 10.1111/j.1365-246X.2012.05419.x
Ekström G,Nettles M,Dziewoński A. 2012. The global CMT project 2004−2010:Centroid-moment tensors for 13 017 earthquakes[J]. <italic>Phys Earth Planet Inter</italic>,<bold>200/201</bold>:1–9.
Guo Z,Gao X,Li T,Wang W. 2018. Crustal and uppermost mantle structures of the South China from joint analysis of receiver functions and Rayleigh wave dispersions[J]. <italic>Phys Earth Planet Inter</italic>,<bold>278</bold>:16–25. doi: 10.1016/j.pepi.2018.03.001
Hunter J D. 2007. Matplotlib:A 2D graphics environment[J]. <italic>Comput Sci </italic>&<italic>Eng</italic>,<bold>9</bold>(3):90–95.
Jost M U,Herrmann R. 1989. A student’s guide to and review of moment tensors[J]. <italic>Seism Res Lett</italic>,<bold>60</bold>(2):37–57. doi: 10.1785/gssrl.60.2.37
Kanamori H. 1993. W phase[J]. <italic>Geophys Res Lett</italic>,<bold>20</bold>(16):1691–1694. doi: 10.1029/93GL01883
Kanamori H,Rivera L. 2008. Source inversion of W phase:Speeding up seismic tsunami warning[J]. <italic>Geophys J Int</italic>,<bold>175</bold>(1):222–238. doi: 10.1111/j.1365-246X.2008.03887.x
Masters G, Woodhouse J H, Freeman G. 2011. Mineos V1.0.2 [CP/OL]. [2019−07−28]. https://geodynamics.org/cig/software/mineos/.
Michelini A,Lomax A. 2004. The effect of velocity structure errors on double-difference earthquake location[J]. <italic>Geophys Res Lett</italic>,<bold>31</bold>(9):L09602.
Waldhauser F,Ellsworth W L. 2000. A double-difference earthquake location algorithm:Method and application to the northern Hayward fault,California[J]. <italic>Bull Seismol Soc Am</italic>,<bold>90</bold>(6):1353–1368. doi: 10.1785/0120000006
Wessel P,Smith W H F. 1991. Free software helps map and display data[J]. <italic>Eos</italic>,<italic>Trans Am Geophys Union</italic>,<bold>72</bold>(41):441–446.
Zhu L P,Helmberger D V. 1996. Advancement in source estimation techniques using broadband regional seismograms[J]. <italic>Bull Seismol Soc Am</italic>,<bold>86</bold>(5):1634–1641.
Zhu L P,Rivera L A. 2002. A note on the dynamic and static displacements from a point source in multilayered media[J]. <italic>Geophys J Inter</italic>,<bold>148</bold>(3):619–627. doi: 10.1046/j.1365-246X.2002.01610.x
Zhu L P. 2003. Recovering permanent displacements from seismic records of the June 9,1994 Bolivia deep earthquake[J]. <italic>Geophys Res Lett</italic>,<bold>30</bold>(14):1740.
-
期刊类型引用(12)
1. 杨敬松，廖府，李智涛，卢泽葳. 一种基于Delaunay三角剖分的天然地震波初至实时拾取算法. 地球物理学进展. 2024(05): 1788-1799 . 
百度学术
2. 邱磊，李彩华. STA/LTA方法拾取天然地震波初至及其改进分析. 地球物理学进展. 2023(04): 1497-1506 . 百度学术
3. 裴玮来，周仕勇. 基于最大顺序统计量的P波初动极性自动读取方法的小震震源机制反演及在云南小江地区的应用研究. 地球物理学报. 2022(03): 992-1005 . 百度学术
4. 颜利君，庞瑶，唐淋. 四川地震编目自动处理系统产出结果分析. 四川地震. 2022(01): 20-23 . 百度学术
5. 张珂，张帆，张晖，王鑫，郝美仙，翟浩. 内蒙古敖汉旗震群的微震匹配定位及其发震断层研究. 地震学报. 2022(04): 594-607 . 本站查看
6. 刘磊，王怀富，王志峰. 基于长短时窗比的地震检测方法应用. 电声技术. 2022(09): 93-96 . 百度学术
7. 金星，林彬华，蔡辉腾，韦永祥，李军，郭阳，于伟恒. 基于气枪信号偏振特征的震相自动判别研究. 地球物理学进展. 2021(01): 8-18 . 百度学术
8. 王旭林，李玉海，王鹏，王晓阳，吕明慧. 基于改进Delaunay三角剖分算法的微地震数据SRV估算与应用. 地球物理学进展. 2021(02): 689-695 . 百度学术
9. 钟媛，胡星星，滕云田，陈波，何朝博，沈晓宇. STA/LTA算法在MEMS烈度计中的适用性分析. 地震研究. 2021(02): 208-215 . 百度学术
10. 周本伟，范莉苹，张龙，李珀任，房立华. 利用卷积神经网络检测地震的方法与优化. 地震学报. 2020(06): 669-683 . 本站查看
11. Haichao Ma，Risheng Chu，Minhan Sheng，Ziye Yu. Template matching for simple waveforms with low signalto-noise ratio and its application to icequake detection. Earthquake Science. 2020(Z1): 256-263 . 必应学术
12. XU Zhen，WANG Tao，XU Shanhui，WANG Baoshan，FENG Xuping，SHI Jing，YANG Minghan. Active Source Seismic Identification and Automatic Picking of the P-wave First Arrival Using a Convolutional Neural Network. Earthquake Research in China. 2019(02): 288-304 . 必应学术
其他类型引用(10)
计量
- 文章访问数: 1137
- HTML全文浏览量: 502
- PDF下载量: 142
- 被引次数: 22
