First arrival time picking algorithm of micro-seismic based on improved STA/LTA and adaptive VMD
-
摘要: 针对低信噪比条件下微震初至拾取准确度低的问题，基于信号幅度变化引入权重因子，对传统长短时窗比值（STA/LTA）算法进行改进，提高初次拾取精度。为了进一步降低拾取误差，对变分模态分解（VMD）算法进行优化，基于互相关系数和排列熵准则自适应确定VMD分解层数，对初次拾取结果前后2—3 s的记录进行优化VMD，并计算分解后各本征模函数（IMF）的峰度赤池信息准则值，得到各IMF的到时，以各IMF的拾取结果及能量比综合加权得到二次拾取到时。仿真实验表明：改进后的STA/LTA在较低信噪比下可降低初次拾取误差约0.01 s以上；相比经验模态分解（EMD)和小波包分解，自适应VMD分解后能再次降低误差，最终与人工拾取结果平均误差在0.023 s以内。实际微震信号初至拾取结果表明，本算法能快速有效地识别初至P波，与人工拾取结果相比误差小，准确率高。Abstract: Accurate and reliable picking of the first arrival time is one of the critical steps in micro-seismic monitoring. Aiming at the problem of low accuracy of first arrival picking for micro-seisms under low signal-to-noise ratio, the traditional short term averaging/long term averaging algorithm is improved by introducing weight factor according to the change of signal amplitude to improve the accuracy of initial pickup. In order to further reduce the pickup error, variational mode decomposition （VMD） is optimized based on cross-correlation coefficient and permutation entropy criterion, and decomposition layers are determined adaptively. Then, the signals of 2−3 s before and after the initial pickup are decomposed by VMD, and the Kurtosis-Akaike information criterion （AIC） values of the decomposed intrinsic mode functions （IMF） are calculated to get the arrival time of each IMF, and the secondary arrival time is obtained by comprehensively weighting the picking results and energy ratios of each IMF. Simulation results show that the improved STA/LTA can reduce the initial picking error by more than 0.01 s at low SNR; compared with empirical mode decomposition （EMD） and wavelet packet decomposition, the adaptive VMD decomposition can reduce the picking error again, and the finalaverage picking error is less than 0.023 s. The first arrival time picking results of real micro-seismic signals show that the proposed algorithm can identify the first break of P-wave quickly and effectively, and the error is smaller than that of manual picking, which shows that the algorithm is effective and the picking accuracy is high.
-
Keywords:
- STA/LTA /
- variational mode decomposition /
- first arrival picking /
- micro-seismic /
- AIC
-
引言
近年来，从强震之间的应力交互角度开展强震间相互作用的研究，已成为震源物理学的一个热点课题（Harris，1998；Perfettini et al，1999 ；Stein，1999；King，Cocco，2001；Steacy et al，2005 ）。通过计算强震序列中前一次事件对后续事件的库仑应力改变量，可以定量地判断前者对后者的触发/延迟作用，帮助研究人员理解强震序列的孕育过程（Wan，Shen，2010；Wang et al，2014 ；Jia et al，2018 ）。同时，大地震的发生将显著影响其周缘地区的应力场状态，并可能造成周围断层提前或推迟发生地震，影响区域地震危险性（Parsons et al，2008 ；Toda et al，2008 ；Hu et al，2017 ）。此外，以统计学的方法研究区域地震活动性的时空演化规律已成为刻画、了解地震事件之间的因果关系及区域构造活动的有力手段（Zhuang et al，2005 ；Peng et al，2012 ；Jia et al，2014 ，2018）。自1970年我国国家地震台网建立以来，地震目录越来越完备，这为统计分析地震活动特征提供了数据基础。在速率-状态本构关系框架下（Dieterich，1979；Ruina，1983），强震引起的库仑应力变化会直接导致地震活动性的时空变化，从而将库仑应力变化和地震活动性变化从理论上衔接起来，为强震之间的触发机制以及未来区域地震危险性的评估提供了有效方法（Aron，Hardebeck，2009；Cocco et al，2010 ；Toda，Enescu，2011；Ishibe et al，2015 ）。
青藏高原是我国最主要的地震活动区之一，尤其是青藏高原中部的巴颜喀拉块体周缘，近20年来发生了一系列强震事件，显示了其构造运动正处于活跃期（程佳等，2011；闻学泽等，2011；Jia et al，2012 ；徐锡伟等，2017）。自本世纪初开始，巴颜喀拉块体周缘相继发生了2001年可可西里MW7.8地震、2008年于田MW7.1地震、2008年汶川MW7.9地震、2010年玉树MW6.9地震、2013年芦山MW6.6地震、2014年于田MW6.9地震和2017年九寨沟MW6.5地震（图1）。这7次强震在时间和空间上的密集发生，造成了巨大的人员伤亡和财产损失，引起了学术界和公众的普遍关注。从动力学的角度看，这些强震的发震机制均与巴颜喀拉块体东南向的运动相吻合，可能是块体构造运动的产物（Jia et al，2012 ；徐锡伟等，2017）。另一方面，一些研究表明，强震间的触发/延迟作用也可能主导了某些强震的发生（Jia et al，2014 ，2018）。因此，深入地细致地研究该强震序列之间可能存在的触发/延迟机制，有助于理解该强震序列的孕震机理和巴颜喀拉块体及周缘构造区的强震危险性。
![]()
图 1 巴颜喀拉块体构造背景及7次强震事件（活动断层数据引自张培震等，2003）Figure 1. The tectonic background of the Bayan Har block and the seven major earthquakes（The active faults data is from Zhang et al，2003 ）多项研究表明，部分强震之间存在应力交互现象，并且强震的发生会对周围其它发震断层产生触发或抑制的作用，影响余震丛集、强震孕育等地震活动（Stein et al，1992 ；Harris，1998；Stein，1999；King，Cocco，2001；Pollitz et al，2003 ；Steacy et al，2005 ）。King等（1994）研究了强震引起的同震库仑应力改变，并运用于1992年MW7.3 Landers地震，成功地解释了Landers地震的余震空间分布，即大多数余震发生在库仑应力增加的区域，只有少量地震发生在库仑应力减小的区域。Freed和Lin （2001）首次使用黏弹性应力触发假说解释了1992年Landers MW7.3地震对1999年Hector Mine MW7.1地震的延迟触发作用，使得黏弹性应力触发机制在地震触发中的作用被广泛关注。1999年土耳其MW7.4 Izmit地震提高了其后续MW7.1 Duzce地震的震源破裂起始点的黏弹性应力，增加量达70% （Hearn et al，2002 ）。上述研究表明各个时空尺度的地震触发已成为一种普遍存在的现象，可为地震孕震机理提供更多的信息。
另一方面，使用统计学的方法分析地震活动的时空特征，进而判断地震之间的触发关系，也成为地震触发关系研究的热点内容之一。其中，结合了大森-宇津（Omori-Utsu）定律和古登堡-里查德（Gutenberg-Richter）定律的传染型余震序列（epidemic type aftershock sequence，简写为ETAS模型）（Ogata，1988，1998；Zhuang et al，2002 ；Ogata，Zhuang，2006），因其提供了一种有效区分背景地震和丛集（触发）地震的方法，而被广泛应用于不同构造背景的地震活动性研究中（Zhuang et al，2005 ；蒋长胜，庄建仓，2010；Peng et al，2012 ；Jia et al，2014 ，2018）。Dieterich （1994）从理论上将地震活动性与库仑应力改变直接联系起来，证实了两者之间存在定量的正相关关系，从而为结合地震活动性分析和库仑应力改变来探究地震触发关系提供了理论依据。因此，同时考虑地震活动性变化和库仑应力改变，能够互相印证地甄别地震之间的触发关系，有助于深入理解区域构造活动和地震危险性。
巴颜喀拉块体周缘强震序列之间是否存在触发关系是科学界普遍关注的问题。尤其是关于2008年汶川地震与2013年芦山地震及2017年九寨沟地震之间，以及2008年和2014年两次于田地震之间的关系，科学界存在不同的看法。刁桂苓等（2010）研究发现，于田地区的地震活动性在汶川地震发生之后显著增强，进而推断2008年于田地震和汶川地震之间可能存在触发关系。而关于芦山地震是否是汶川地震的余震这一科学问题，学术界仍存在争论：既有学者认为是余震的（陈运泰等，2013；王为民等，2013），也有研究推断是构造事件的（杜方等，2013；Xu et al，2013 ）；同时，一些研究人员认为芦山地震的发生既有汶川地震的触发作用，也有构造应力加载的贡献（Jia et al，2014 ；Wang et al，2014 ）；Jia等（2014）定量分析了汶川地震对芦山地区背景地震活动性的改变，认为芦山地震是汶川地震余震的概率是62%；Wang等（2014）以构造应力和库仑应力改变的比例为判断依据，推断芦山地震是汶川地震余震的概率是85%。针对2017年九寨沟地震与汶川地震的关系，也存在不同的声音：一些学者认为前者是后者触发引起的（单斌等，2017；Hu et al，2017 ），另一些研究则表明后者延迟了前者的发震时间（Wang et al，2014 ；Jia et al，2018 ）。关于两次于田地震之间的关系，多项研究均表明，2008年于田地震所引起的静态和黏弹性应力的增加提前了2014年于田地震的发生（Li et al，2015b ；Wang et al，2017 ）。程佳等（2011）考虑黏弹性应力作用，计算了玛尼地震、可可西里地震和汶川地震所形成的同震和震后形变场的变化过程和特征，认为各地震之间的影响很微弱。以上针对巴颜喀拉块体周缘强震之间触发关系的研究，对理解该强震序列的发震机理是构造运动还是触发作用，或者两者皆有，均具有重要的科学意义。
本文拟对自本世纪初起发生在巴颜喀拉块体边界的7次强震之间的触发关系进行系统研究（表1）。以库仑应力改变量和地震活动性的改变为依据，甄别地震之间的交互作用。根据发震时间和震中距离的不同，本文对汶川地震−芦山地震−九寨沟地震、两次于田地震进行单独研究，将其余强震作为一个序列进行分析讨论，再根据这些地震之间的触发关系，推断巴颜喀拉块体现今的构造运动和未来的地震危险性。
表 1 巴颜喀拉块体周缘7次强震的震源机制解Table 1. The focal mechanism solutions of seven major earthquakes occurred around the Bayan Har block序号 发震时间
年-月-日东经/° 北纬/° MW 走向/° 倾角/° 滑动角/° 深度/km 1 2001−11−14 92.91 35.80 7.8 94 51 −12 15.0 2 2008−03−21 81.37 35.43 7.1 203 52 −74 12.0 3 2008−05−12 104.10 31.44 7.9 231 35 138 12.8 4 2010−04−14 96.79 33.05 6.9 300 88 23 15.7 5 2013−04−20 103.12 30.22 6.6 212 42 100 21.9 6 2014−02−12 82.57 36.22 6.9 242 86 −5 18.3 7 2017−08−08 103.90 33.21 6.5 150 78 −13 14.9 注：数据参考Ekström等 （2012）。 1. 库仑应力改变
由于各强震间发震时刻的间隔（数年）远大于动态应力触发的特征时间（数小时至数周），因此本文不考虑该强震序列的动态应力改变，仅关注同震时的静态应力改变和震后的黏弹性应力改变。本研究根据巴颜喀拉块体深部结构的研究结果，将地球模型简化为包含弹性上地壳、黏弹性下地壳和黏弹性上地幔的层状模型，其中下地壳和上地幔的粘性系数分别为1×1018 Pa·s和1×1020 Pa·s，模型具体细节请参考Xiong等（2010）和Jia等（2012）。由于汶川地震与芦山地震、九寨沟地震以及两次于田地震之间的空间距离较近，且相关研究较多（Wang et al，2014 ；赵立波等，2016；单斌等，2017）。为了方便对比研究，分别将汶川地震对芦山地震和九寨沟地震的影响，芦山地震对九寨沟地震的影响，以及2008年于田地震对2014年于田地震的影响进行单独分析。同时定量地计算可可西里地震、2008年于田地震、汶川地震和玉树地震之间的库仑应力改变量。
1.1 库仑应力改变的计算方法
库仑应力改变的计算公式为（Robinson，2005；Jia et al，2014 ）
${\rm d}{\rm{CFS}} = {\rm d}{\tau _{\rm{s}}} + \mu \left( {{\rm d}{\tau _{\rm{n}}} - \left( {\frac{\beta }{3}} \right)\sum {{\tau _{ii}}} } \right),$
(1) 式中，dτs为剪切应力改变量，μ为摩擦系数，dτn为正应力改变量（拉张为正），β为Skempton系数。相比于部分研究使用视摩擦系数综合孔隙压的变化，我们更倾向于使用Skempton系数来更准确地刻画孔隙压力的变化，因为前者可能导致库仑应力计算结果不准确（Beeler et al，2000 ）。Skempton系数根据实验室和野外观测的结果取0.5，而摩擦系数因不同发震断层所处的环境不同而取不同数值，取值范围介于0—1.0之间 （Robinson，2005；Jia et al，2014 ）。使用PSGRN/PSCMP程序（Wang et al，2006 ），给定地震滑动分布模型，接受断层的发震机制，即可根据黏弹性位错模型计算出强震引起的静态及黏弹性应力改变量。
1.2 可可西里地震、2008年于田地震、汶川地震和玉树地震之间的库仑应力改变
可可西里地震、2008年于田地震和汶川地震的破裂模型分别参考Lasserre等（2005）、张国宏等（2011）以及Ji和Hayes （2008）的结果。摩擦系数和Skempton系数根据实验室和野外观察的结果，分别选取0.75和0.5 （Robinson，2005）。库仑应力的计算结果列于表2。其中由于2008年于田地震和汶川地震的震中相距非常远（2 500 km），而同震应力改变与两次地震的震中距离呈二次反比的关系，因此这两次地震之间的同震应力改变可以忽略；同时这两次地震的发震时间间隔较短（52天），远小于震后黏弹性应力释放发挥作用的特征时间，因此震后黏弹性应力改变也可以忽略，所以表2中并未列出这两次地震之间的库仑应力改变量。从表2可见，可可西里地震对后续2008年于田、汶川和玉树地震的库仑应力改变非常小，远小于可能引起地震活动性变化的阈值0.01 MPa （Reasenberg，Simpson，1992；Stein et al，1992 ）。类似的，2008年于田地震和汶川地震对玉树地震的库仑应力改变量也非常小。由上述结果可以推断，这4次地震事件之间并无触发作用。
表 2 巴颜喀拉强震序列在其后续地震发震位置所引起的库仑应力变化Table 2. The dCFS results induced by the strong earthquake sequence occurred around the Bayan Har block on following earthquakes2001年可可西里地震 2010年玉树地震 同震应力改变
/MPa震后应力改变
/MPa同震应力改变
/MPa震后应力改变
/MPa于田地震 7.00×10−5 1.30×10−4 7.28×10−7 8.00×10−7 汶川地震 1.14×10−5 1.20×10−5 −1.32×10−4 −2.50×10−4 玉树地震 −7.15×10−4 −1.35×10−3 注：震后应力改变指同震和震后黏弹性松弛共同导致的结果。 1.3 汶川地震、芦山地震和九寨沟地震之间的库仑应力改变
为了讨论汶川地震对芦山地震和九寨沟地震的影响，以及芦山地震对九寨沟地震的影响，本文使用不同的滑动分布模型和摩擦系数计算强震之间的库仑应力改变量，计算结果如图2所示。汶川地震的滑动模型使用Ji和Hayes （2008）以及Wang等（2011）的反演结果，芦山地震使用Jiang等（2014）和Zhang等（2014）的反演结果。由图2a可知，分别使用Wang等（2011）、Ji和Hayes （2008）的破裂模型，汶川地震对芦山地震的同震库仑应力改变分别达到了0.006—0.009 MPa和0.010—0.012 MPa，而总体的库仑应力改变（同震与震后库仑应力改变之和）分别为0.031—0.034 MPa和0.041—0.049 MPa。显然，这种程度的库仑应力改变量表明汶川地震对芦山地震具有一定的触发作用，也与其它研究结果相似（Parsons et al，2008 ；Toda et al，2008 ；单斌等，2009；Wan，Shen，2010；Jia et al，2014 ；Wang et al，2014 ）。由图2b可知，除了使用Wang等（2011）的模型在μ=0.8的条件下，库仑应力改变为正值外，在μ=0.4以及使用Ji和Hayes （2008）的模型时，库仑应力改变量均为负值。这是由库仑应力改变的计算结果对摩擦系数有较高的敏感度造成的，也能解释为何针对汶川地震对九寨沟地震的库仑应力改变，不同研究可能得到相反的结果（Parsons et al，2008 ；Toda et al，2008 ；Luo，Liu，2010；Wan，Shen，2010；Nalbant，McCloskey，2011；Wang et al，2014 ；Hu et al，2017 ）。通常，断层面上摩擦系数的精确值很难获得，因而仅假设其取值范围介于0—1.0之间。Li等（2015a）通过汶川地震断裂带科学钻探项目直接观测到断层面上的温度数据，推断龙门山断裂带处于低摩擦环境（μ<0.2）。因此，使用低摩擦系数计算得到的负库仑应力改变量更为可信，表明汶川地震可能对九寨沟地震具有一定的延迟作用。由图2c可知，芦山地震对九寨沟地震的库仑应力改变量非常微弱，约为10−6 MPa，可以推断芦山地震对九寨沟地震并无影响。
![]()
图 2 强震之间的库仑应力改变随时间的变化（a） 汶川地震对芦山地震的库仑应力改变；（b） 汶川地震对九寨沟地震的库仑应力改变；（c） 芦山地震对九寨沟地震的库仑应力改变Figure 2. Coulomb stress changes dCFS of strong earthquakes versus time（a） The dCFS induced by the Wenchuan earthquake on the Lushan earthquake；（b） The dCFS induced by the Wenchuan earthquake on the Jiuzhaigou earthquake；（c） The dCFS induced by the Lushan earthquake on the Jiuzhaigou earthquake1.4 两次于田地震之间的库仑应力改变
使用同样的方法计算两次于田地震之间的库仑应力改变量（图3），解析深度为10 km。由图3可知，2008年于田地震对2014年于田地震的库仑应力改变量从同震时的约−0.001 MPa至2014年于田地震前的约−0.003 MPa，虽然是负值，但其强度较小。该结果与前人的研究结果并不完全一致（Li et al，2015b ；王辉等，2016），但各项计算结果均显示库仑应力改变量非常小（赵立波等，2016；Wang et al，2017 ），表明这两次地震事件之间触发关系很微弱。
![]()
图 3 2008年于田地震对2014年于田地震的库仑应力改变量Figure 3. The dCFS induced by the 2008 Yutian earthquake on the 2014 Yutian earthquake2. 地震活动性改变
2.1 ETAS模型和随机除丛法
ETAS模型结合了古登堡-里查德定律、大森-宇津定律及地震发生率与震级的关系等统计学规律，将某一地区某一时段内大于某震级的地震活动性描述为点过程，并将整体地震活动性看作不随时间变化的背景活动性和余震活动性之和。最重要的是，ETAS模型中的每次地震（包括余震）均有能力触发其余震，这种触发能力与其震级大小有关。ETAS模型以概率的形式来定义每次地震，即每次地震均有一定概率是其之前地震的余震，余下的概率作为背景地震的概率。然而在现实中无法实际观测到每次地震与其余震之间的关系，而只能记录到地震的时间序列。但ETAS模型能够通过对某一地区的地震目录建模，拟合其模型参数，从而分辨出各地震的余震概率和背景概率。
ETAS模型的数学含义为：给定时刻 t之前的地震目录，在单位时间内，位置（x，y）处震级大于完备震级Mc的地震数目，其表达式为（Ogata，1998；Zhuang et al，2005 ）
$\lambda \left( {t,x,y} \right) = \mu \left( {x,y} \right) + \sum\limits_{\left( {i:{t_i} < t} \right)} {\kappa \left( {{M_i}} \right)g\left( {t - {t_i}} \right)f\left( {x - {x_i},y - {y_i};{M_i}} \right)} ,$
(2) $\kappa \left( M_i \right) = A{{\rm{e}}^{\alpha \left( {M_i - {M_{\rm{c}}}} \right)}},\quad M_i {\text{≥}} {M_{\rm{c}}},$
(3) $g\left( t \right) = \frac{{p - 1}}{c}{\left( {1 + \frac{t}{c}} \right)^{ - p}},$
(4) $f\left( {x,y;M_i} \right) = \frac{{q - 1}}{{{\rm{\pi }}D{{\rm{e}}^{\gamma \left( {M_i - {M_c}} \right)}}}}{\left[ {1 + \frac{{{x^2} + {y^2}}}{{D{{\rm{e}}^{\gamma \left( {M_i - {M_{\rm{c}}}} \right)}}}}} \right]^{ - q}},$
(5) 式中：μ（x，y）为背景地震活动率，假设其不随时间变化；κ（Mi）为地震发生率，即震级为Mi的主震能够触发的余震数目；g（t）为大森-宇津定律的概率密度函数；f（x，y；Mi）为空间概率密度函数；A，α，p，c，q，D，γ均为模型参数，表示给定区域地震活动性的统计特征。给定完备的地震目录和合适的模型参数初始值，通过最大似然法即可估计出各模型参数的值。
2.2 ETAS模型拟合
将整个巴颜喀拉块体看作一个整体进行ETAS模型拟合，空间范围为（29°N—37°N，79°E—105°E），时间范围为1980年1月1日至2017年8月27日。前人的研究表明，强震发生之后存在小地震无法被记录到的现象（Enescu et al，2007 ；Peng et al，2012 ），导致地震目录不完备，进而引起ETAS模型拟合参数的不稳定（Utsu et al，1995 ；Ogata，1998）。这一问题可以通过选用高于完备震级的地震目录来解决。由于观测到汶川地震之后大量小震（M≤4.0）缺失，因此本研究选取的完备震级为4.0 （Jia et al，2014 ，2018）。利用最大似然法估计ETAS模型参数，最终结果为μ=0.754，A=0.474，c=0.025，α=0.839，p=1.134，D=2.848×10−2，q=1.765，γ=0.46。
2.3 背景地震活动性变化
得到ETAS模型参数之后，可以计算地震事件j处背景地震活动率与整体地震活动率的比值，即地震事件j的背景概率，其表达式为
${\varphi _j} = \frac{{\mu \left( {{x_j},{y_j}} \right)}}{{\lambda \left( {{t_j},{x_j},{y_j}} \right)}}.$
(6) 为了描述背景地震活动率随时间的变化，Zhuang等（2005）定义了累积背景地震概率S（t）为
$S(t) = \sum\limits_{{t_i} < t} {{\varphi _i}} .$
(7) 对于构造运动稳定的地区，背景地震活动性不随时间变化，因此S（t）曲线应为随时间稳定上升的直线。当S（t）曲线发生突然变化时，预示着该区域背景地震活动性的突变。
根据活动构造与地震活动性空间分布的特征，在7次强震周围分别选取6个研究区（图4），展示其各自的地震活动性随时间的变化规律。由于2008年与2014年两次于田地震的震中位置较近，故将这两次地震的发震区域作为同一个研究区来分析地震活动性的变化。
![]()
图 4 青藏高原北部地震活动性空间分布及各研究区空间范围Figure 4. Seismicity of the northern Tibetan Plateau and the spatial region of the studied area图5给出了各个研究区的总体地震活动性和背景地震概率累积曲线，可以根据后者随时间的变化推断各地区构造活动的变化，进而推断该研究区被强震影响的可能性。从图5a可以看出，2001年可可西里地震之后，该地区背景地震活动性显著升高，表明该地区构造运动活跃；从图5b可以看出，于田地区背景地震活动性在每次强震发生时均无明显变化，表明其它地震对两次于田地震并无明显影响，两次于田地震之间也无显著触发/延迟关系。同时也可以看出，两次于田地震也并未提升该地区的背景地震活动性；类似的，汶川地区背景地震活动性在可可西里地震时并未出现显著变化，而2008年汶川地震显著提高了该地区的背景地震活动性（图5c）；玉树地区在可可西里地震、2008年于田地震和汶川地震前后背景地震活动性也均未发生显著变化（图5d），指示了玉树地区并未受这3次地震事件的影响；而芦山地区的背景地震活动性也仅仅在2008年汶川地震之后有明显升高，表明汶川地震对芦山地震有触发作用（图5e）；九寨沟地区在汶川地震后背景地震活动性出现明显降低，而在其它地震发生时则无显著变化（图5f），表明汶川地震对九寨沟地区地震活动性有压制作用。
综合库仑应力改变和地震活动性变化的结果可以看出，各强震的库仑应力变化结果与其造成的背景地震活动性变化相一致。在该强震序列中，汶川地震对芦山地震和九寨沟地震均有显著的库仑应力改变以及背景地震活动性变化，而其它地震之间均无明显的触发/延迟作用以及显著的背景地震活动性变化。
![]()
图 5 累积地震次数和背景地震概率（a） 可可西里；（b） 于田；（c） 汶川；（d） 玉树；（e） 芦山；（f） 九寨沟Figure 5. The cumulative total earthquakes and the cumulative background earthquake probabilities（a） Hoh Xil；（b） Yutian；（c） Wenchuan；（d） Yushu；（e） Lushan；（f） Jiuzhaigou3. 讨论与结论
强震之间的应力交互是一个复杂的物理过程。由于计算参数的不确定性，例如缺乏约束的摩擦系数或滑移分布模型的不一致性等，往往很难准确地计算出强震引起的库仑应力改变。此外，观测到的地震活动性不仅包含来源于构造运动的背景地震事件，还包含了强震触发引起的触发地震。利用ETAS模型能够有效地区分出背景地震和触发地震，从而为分析背景地震活动性的时空变化提供基础。库仑应力改变与地震活动性改变呈正相关，通过速率-状态摩擦率有机地联系在一起，为研究强震之间的触发机制提供有力的证据。
总体来说，本文针对巴颜喀拉块体周缘的强震序列，从库仑应力改变和地震活动性这两个角度进行探讨分析，试图辨别这一备受关注的强震序列是构造运动的产物还是触发作用的结果，并进一步讨论了巴颜喀拉块体未来的地震危险性。库仑应力改变和地震活动性分析的结果均表明，汶川地震对芦山地震的触发作用和对九寨沟地震的延迟作用较为明显，展现出局部地区也存在强震之间互相影响的特征。除此之外，其它地震之间均无明显的触发或延迟作用。这表明该强震序列的孕震机制主要来源于巴颜喀拉块体不断向东南挤压的构造运动。因此，可以推断巴颜喀拉块体目前的构造运动较为活跃，该区域未来的地震危险性不可忽视。同时，强震触发/延迟周围活动断层的破裂也是值得关注的现象。
中国地震局地球物理研究所蒋长胜研究员提供了整理和修正后的地震目录，德国波茨坦地学研究中心汪荣江研究员提供了PSGRN/PSCMP的程序，日本数理统计研究所庄建仓教授提供了ETAS模型和随机除丛法的程序，两位匿名评审专家为本文的撰写提出了宝贵的意见和建议，作者在此一并表示衷心的感谢。
-
![]()
图 1 基于改进STA/LTA及优化VMD的微震初至拾取算法流程图
Figure 1. The flow chart of the first break picking of micro-seismic based on imp-roved STA/LTA and adaptive VMD
![]()
图 2 基于雷克子波和卷积模型的理论和模拟地震记录
Figure 2. Theoretical and simulated seismogram based on Ricker wavelet and convolution model
![]()
图 4 不同SNR下三种STA/LTA算法拾取误差曲线
Figure 4. Picking up error curves of three STA/LTA algorithms under different SNRs
![]()
图 5 优化VMD峰度AIC为3.07 s时的初至拾取结果
Figure 5. First arrival picking based on improved VMD when Kurtosis-AIC is 3.07 s
![]()
图 7 基于小波包峰度AIC的初至拾取
Figure 7. First arrival picking based on wavelet packet Kurtosis-AIC
![]()
图 10 广西平果2013年7月4日ML1.3爆破的一条实际微震记录
Figure 10. Real micro-seismic record of some coal mine with ML1.3 at Pingguo，Guangxi on July 4th 2013
![]()
图 11 本文改进STA/LTA初次拾取结果
Figure 11. First picking results of improved STA/LTA in this paper
![]()
图 12 本文改进STA/LTA及优化VMD峰度AIC拾取结果
Figure 12. First arrival picking by improved STA/LTA and adaptive VMD Kurtosis-AIC
表 1 实际微震记录初至拾取结果
Table 1 First arrival picking of real micro-seismics
拾取方法 误差 30 ms 误差 20 ms 误差 10 ms 记录条数 占比 记录条数 占比 记录条数 占比 本文 1 117 98.2% 1 090 95.9% 1 031 90.7% EMD-AIC 1 090 95.9% 1 052 92.5% 980 86.2% 小波包峰度AIC 1 075 94.5% 1 037 91.2% 968 85.1% 
下载: 导出CSV
-
贾瑞生,谭云亮,孙红梅,洪永发. 2015. 低信噪比微震P波震相初至自动拾取方法[J]. 煤炭学报,40(8):1845–1852. Jia R S,Tan Y L,Sun H M,Hong Y F. 2015. Method of automatic detection on micro-seismic P-arrival time under low signal-to-noise ratio[J]. Journal of China Coal Society,40(8):1845–1852 (in Chinese).
李伟,江晓林,陈海波,金珠鹏,刘志军,李兴伟,林井祥. 2018. 基于EEMD_Hankel_SVD的矿山微震信号降噪方法[J]. 煤炭学报,43(7):1910–1917. Li W,Jiang X L,Chen H B,Jin Z P,Liu Z J,Li X W,Lin J X. 2018. Denosing method of mine microseismic signal based on EEMD_Hankel_SVD[J]. Journal of China Coal Society,43(7):1910–1917 (in Chinese).
刘晗,张建中. 2014. 微震信号自动检测的STA/LTA算法及其改进分析[J]. 地球物理学进展,29(4):1708–1714. doi: 10.6038/pg20140429 Liu H,Zhang J Z. 2014. STA/LTA algorithm analysis and improvement of microseismic signal automatic detection[J]. Progress in Geophysics,29(4):1708–1714 (in Chinese).
刘晓明,赵君杰,王运敏,彭平安. 2017. 基于改进的STA/LTA方法的微地震P波自动拾取技术[J]. 东北大学学报(自然科学版),38(5):740–745. Liu X M,Zhao J J,Wang Y M,Peng P A. 2017. Automatic picking of microseismic events P-wave arrivals based on improved method of STA/LTA[J]. Journal of Northeastern University (Natural Science)
,38(5):740–745 (in Chinese). 唐贵基,王晓龙. 2015. 参数优化变分模态分解方法在滚动轴承早期故障诊断中的应用[J]. 西安交通大学学报,49(5):73–81. doi: 10.7652/xjtuxb201505012 Tang G J,Wang X L. 2015. Parameter optimized variational mode decomposition method with application to incipient fault diagnosis of rolling bearing[J]. Journal of Xi’an Jiaotong University,49(5):73–81 (in Chinese).
田优平,赵爱华. 2016. 基于小波包和峰度赤池信息量准则的P波震相自动识别方法[J]. 地震学报,38(1):71–85. doi: 10.11939/jass.2016.01.007 Tian Y P,Zhao A H. 2016. Automatic identification of P-phase based on wavelet packet and kurtosis-AIC method[J]. Acta Seismologica Sinica,38(1):71–85 (in Chinese).
王志坚,常雪,王俊元,杜文华,段能全,党长营. 2018. 排列熵优化改进变模态分解算法诊断齿轮箱故障[J]. 农业工程学报,34(23):59–66. doi: 10.11975/j.issn.1002-6819.2018.23.007 Wang Z J,Chang X,Wang J Y,Du W H,Duan N Q,Dang C Y. 2018. Gearbox fault diagnosis based on permutation entropy optimized variational mode decomposition[J]. Transactions of the Chinese Society of Agricultural Engineering,34(23):59–66 (in Chinese).
赵大鹏,刘希强,刘尧兴,王志铄,赵晖,张亚琳. 2013. 高阶统计量及AIC方法在区域地震事件和直达P波初动识别中的应用[J]. 地震地磁观测与研究,34(5/6):61–69. Zhao D P,Liu X Q,Liu Y X,Wang Z S,Zhao H,Zhang Y L. 2013. Detection of regional seismic events by high order statistics method and automatic identification of direct P-wave first motion by AIC method[J]. Seismological and Geomagnetic Observation and Research,34(5/6):61–69 (in Chinese).
郑近德,程军圣,杨宇. 2013. 改进的EEMD算法及其应用研究[J]. 振动与冲击,32(21):21–26. doi: 10.3969/j.issn.1000-3835.2013.21.004 Zheng J D,Cheng J S,Yang Y. 2013. Modified EEMD algorithm and its applications[J]. Journal of Vibration and Shock,32(21):21–26 (in Chinese).
郑小霞,周国旺,任浩翰,符杨. 2017. 基于变分模态分解和排列熵的滚动轴承故障诊断[J]. 振动与冲击,36(22):22–28. Zheng X X,Zhou G W,Ren H H,Fu Y. 2017. A rolling bearing fault diagnosis method based on variational mode decomposition and permutation entropy[J]. Journal of Vibration and Shock,36(22):22–28 (in Chinese).
朱权洁,姜福兴,魏全德,王博,刘金海,刘晓辉. 2018. 煤层水力压裂微震信号P波初至的自动拾取方法[J]. 岩石力学与工程学报,37(10):2319–2333. Zhu Q J,Jiang F X,Wei Q D,Wang B,Liu J H,Liu X H. 2018. An automatic method determining arrival times of microseismic P-phase in Hydraulic fracturing of coal seam[J]. Chinese Journal of Rock Mechanics and Engineering,37(10):2319–2333 (in Chinese).
Bandt C,Pompe B. 2002. Permutation entropy:A natural complexity measure for time series[J]. Phys Rev Lett,88(17):174102. doi: 10.1103/PhysRevLett.88.174102
Charles M,Ge M C. 2018. Enhancing manual P-phase arrival detection and automatic onset time picking in a noisy microseismic data in underground mines[J]. Int J Min Sci Technol,28(4):691–699. doi: 10.1016/j.ijmst.2017.05.024
Dragomiretskiy K,Zosso D. 2014. Variational mode decomposition[J]. IEEE Trans Signal Proc,62(3):531–544. doi: 10.1109/TSP.2013.2288675
Gaci S. 2014. The use of wavelet-based denoising techniques to enhance the first-arrival picking on seismic traces[J]. IEEE Trans Geosci Remote Sens,52(8):4558–4563. doi: 10.1109/TGRS.2013.2282422
Kirbas I,Peker M. 2017. Signal detection based on empirical mode decomposition and Teager-Kaiser energy operator and its application to P and S wave arrival time detection in seismic signal analysis[J]. Neural Comput Appl,28(10):3035–3045. doi: 10.1007/s00521-016-2333-5
Li F Y,Zhang B,Verma S,Marfurt K J. 2018. Seismic signal denoising using thresholded variational mode decomposition[J]. Explora Geophys,49(4):450–461. doi: 10.1071/EG17004
Li X B,Shang X Y,Wang Z W,Dong L J,Weng L. 2016. Identifying P-phase arrivals with noise:An improved Kurtosis method based on DWT and STA/LTA[J]. J Appl Geophys,133:50–61. doi: 10.1016/j.jappgeo.2016.07.022
Li X B, Shang X Y, Morales-Esteban A, Wang Z W. 2017. Identifying P phase arrival of weak events: The Akaike information criterion picking application based on the empirical mode decomposition[J] Comput Geosci, 100: 57–66.
Liu M Z,Yang J X,Cao Y P,Fu W N,Cao Y L. 2017. A new method for arrival time determination of impact signal based on HHT and AIC[J]. Mech Syst Signal Proc,86:177–187. doi: 10.1016/j.ymssp.2016.10.003
Shang X Y,Li X B,Morales-Esteban A,Dong L J. 2018. Enhancing micro-seismic P-phase arrival picking:EMD-cosine function-based denoising with an application to the AIC picker[J]. J Appl Geophys,150:325–337. doi: 10.1016/j.jappgeo.2017.09.012
Xue Y J,Cao J X,Wang D X,Du H K,Yao Y. 2016. Application of the variational-mode decomposition for seismic time-frequency analysis[J]. IEEE J Select Top Appl Earth Observ Remote Sens,9(8):3821–3831. doi: 10.1109/JSTARS.2016.2529702
-
期刊类型引用(12)
1. 杨磊，陈双贵，闫万生，张瑜，雷光，杨龙. 南北地震带北段及周边地区岩石圈磁场时空变化特征对地震活动的指示意义. 大地测量与地球动力学. 2024(03): 304-309 . 
百度学术
2. 宋成科，张鹏涛，陈斌. 基于地磁场重复观测资料分析玛多M_W7.3地震前地磁场变化. 防灾科技学院学报. 2024(01): 38-46 . 百度学术
3. 陈政宇，倪喆，周思远，金云华，杨薪俊. 岩石圈磁场与地质构造和地震活动性之间的关系——以漾濞地震为例. 地震地质. 2024(02): 449-461 . 百度学术
4. 王朝景，李博，苏树朋. 基于多期累积岩石圈磁场变化分析唐山M_S5.1地震震磁异常. 地震研究. 2024(04): 517-527 . 百度学术
5. 李晨阳，池成全. 机器学习在地震观测异常数据提取中的应用. 海南师范大学学报(自然科学版). 2024(03): 348-356 . 百度学术
6. 张瑜，陈双贵，闫万生，雷光，杨磊，杨龙，马辉源，肖世堂，董兴洲，岳敏. 甘肃青海部分地区流动地磁场时空分布特征. 大地测量与地球动力学. 2023(01): 65-70 . 百度学术
7. 张瑜，陈双贵，马辉源，杨磊，雷光，肖世堂，闫万生，杨龙. 2019—2021年门源M_S6.9地震岩石圈磁场异常回溯分析. 华南地震. 2023(01): 46-54 . 百度学术
8. 张瑜，陈双贵，马辉源，闫万生，杨磊，雷光，杨龙，杜建清. 2021年5月22日玛多M_S 7.4地震周边地区岩石圈磁场变化及地震前后异常特征分析. 地震地磁观测与研究. 2023(03): 64-72 . 百度学术
9. 张海洋，苏树朋，赵慧琴. 2022年青海门源6.9级地震前岩石圈磁场异常变化分析. 地震工程学报. 2022(03): 735-743 . 百度学术
10. 蔡苏苏，陈斌. 中国大陆岩石圈震磁异常年变统计分析. 地震研究. 2022(04): 592-598 . 百度学术
11. 马永，张海江，高磊，陈志刚. 2021年玛多Ms7.4地震三维地壳速度结构与活动特征（英文）. Applied Geophysics. 2022(04): 590-602+605 . 百度学术
12. 文丽敏，康国发，白春华，高国明. 南北地震带南段地壳磁异常与强震活动关系研究（英文）. Applied Geophysics. 2021(03): 408-419+434 . 百度学术
其他类型引用(1)
计量
- 文章访问数: 730
- HTML全文浏览量: 326
- PDF下载量: 125
- 被引次数: 13
