A comparative study on three methods of onset-time determination for teleseismic P arrivals and parameters optimization
-
摘要: 分析了自回归赤池信息准则（AR-AIC）、高阶统计量（HOS）和累积和（CUSUM）等3种到时拾取方法中参数对远震P波到时估计的影响，以450个远震P波信号为样本集，参考人工拾取到时，以网格搜索方式确定了每种方法的最优参数。之后重新选取信噪比处于[2，20]区间的100个远震P波信号，用确定最优参数后的3种到时拾取方法估计其P波初至时间，并比较了3种方法对低信噪比远震P波的拾取准确度。结果表明，AR-AIC方法和CUSUM方法对低信噪比远震P波的拾取准确度要优于HOS方法，CUSUM方法的计算速度最快，HOS方法由于其原理的限制更适用于信噪比较大、初动较尖锐的信号。
-
关键词:
- 到时拾取 /
- 自回归赤池信息准则（AR-AIC）方法 /
- 高阶统计量（HOS）方法 /
- 累积和（CUSUM）方法
Abstract: This paper analyzed the effect of parameters on estimation of onset times for teleseismic P waves by AR-AIC （auto-regression Akaike information criterion）, HOS （higher order statistics） and CUSUM （cumulative sums） methods. The parameters of each method were optimized in the way of grid searching by taking 450 teleseismic P phases as sample set and the manual picked arrival times as reference. And then 100 teleseismic P phases were reselected with SNR between 2 and 20, and the onset times were estimated by the three methods with optimized parameters, so as to evaluate the accuracy of estimation for low-SNR teleseismic P phases. The results show that the accuracy of estimated phase arrival times by the AR-AIC and CUSUM methods is better than by the HOS method, the calculation speed by the CUSUM method is the fastest, and the HOS method is better applied to high-SNR and sharp onset signals.-
Keywords:
- automated picking /
- AR-AIC method /
- HOS method /
- CUSUM method
-
引言
随着现代通信手段和数字地震观测技术的发展，地震数据自动处理技术已经广泛应用于地震日常监测分析中，其中地震信号到时的自动拾取是自动处理过程中的重要一环。到时估算的准确性越高，事件关联和定位的可靠性也越高，所以研究地震信号到时的自动拾取技术具有十分重要的意义。
震相到时拾取方法主要依据地震记录中信号到达前后振幅、频率、偏振特性等特征的不同将信号从噪声中识别出来，最基础的震相到时拾取方法为长短时平均能量比（short term average to long term average，简写为STA/LTA）算法，该算法由Allen （1978, 1982）提出，利用长短时间窗的振幅能量比来识别信号；Baer和Kradolfer （1987）改进了该方法，并引入动态门限。然而这两种方法所确定的初至时间精度不高，多用于信号的检测。目前应用较多的是以赤池信息准则（Akaike information criterion，简写为AIC）（Akaike，1973）为基础而发展起来的基于自回归模型的AIC （auto-regression AIC，简写为AR-AIC）方法（Sleeman，van Eck，1999；王海军等，2003）、计算方差的AIC （variance-AIC，简写为VAR-AIC）方法（Maeda，1985）和计算三阶统计量的AIC （three-order cumulative AIC，简写为TOC-AIC）方法（周银兴，2009），此外还有偏振分析方法（Jurkevics，1988；Kedrov，Ovtchinnikov，1990）、分形分维方法（曾富英等，2002）、基于时频分析的方法（岳龙等，2015）、高阶统计量（higher order statistics，简写为HOS）方法（Saragiotis et al，2002 ）和累积和（cumulative sums，简写为CUSUM）方法（Iclán，Tiao，1994；Der，Shumway，1999，2003；何燕等，2014）等。
上述方法均采用单一特征估算震相到时，各有其局限性和适用条件，后来发展起来的综合分析方法，例如，基于峰度的AIC （kurtosis-AIC，简写为kur-AIC）方法（Küperkoch et al，2010 ）、STA/LTA和峰度或者卓越周期和峰度联合的方法（Nippress et al，2010 ）、基于小波分解与高阶统计量的方法（盛冠群等，2015）、基于小波包和峰度AIC的拾取方法（田优平，赵爱华，2016）、基于STA/LTA与多道互相关结合的拾取方法（林凡生，邹志辉，2016）和STA/LTA、偏振分析与AR-AIC相结合的微地震事件初至拾取方法（谭玉阳等，2016），则综合了多种单一方法的优点。这些综合分析方法虽然拾取震相到时的精度较高，但计算量也相对较大。
在多种震相到时自动拾取方法中，AR-AIC方法、HOS方法和CUSUM方法均以振幅变化为特征，且具有快速简单的特点。AR-AIC方法基于信息理论，精度较高（Küperkoch et al，2012 ），全面禁核试条约国际数据处理中心的自动处理系统（Greg et al，2001 ）便是使用该方法估算信号到时；HOS方法的灵敏性好，Küperkoch等（2010）将该方法应用于爱琴海南部的区域地震台网监测中，结果表明该方法适用于实时处理系统；Der和Shumway （1999，2003）将CUSUM方法用于区域震震相拾取中，其结果优于自回归模型的测试结果（Der，Shumway，2003；何燕等，2014）。
为提高包括低震级事件在内的所有地震事件自动检测定位结果的准确性，靳平等曾提出同时基于信号整体和局部特征的区域震相识别关联（靳平等，2014；Jin et al，2014 ）和远震事件渐进关联（Jin et al，2015 ）两种区域台网地震数据自动处理新技术，并在此基础上，自主开发了区域台网数据自动处理系统DRSN，其中在到时估算环节采用一种改进的AR-AIC方法（王海军等，2003）。DRSN系统在本单位使用的区域台网上已应用了十余年，取得了较好的实际应用效果，但相对于有经验的分析员的分析结果，尚存在诸如信号到时拾取不够精确、关联信号不够完整、部分事件定位误差较大等问题，有时甚至会出现虚假事件。类似问题在目前各种地震数据自动处理系统包括全面禁核试条约国际数据中心的地震数据自动处理系统中也都存在，甚至非常严重（靳平等，2014）。针对现有系统存在的这些问题，着眼于进一步提高地震数据自动处理结果的准确性和完整性，实现日常地震监测的全过程自动处理，现阶段我们正在进一步开展地震监测计算机分析员关键技术研究和原型系统研制（Jin，2015）。而作为支持该系统研制的一部分，本文拟对AR-AIC方法、HOS方法和CUSUM方法这3种常用的到时拾取方法进行优化，并比较其拾取结果，给出其优点和适用范围，以期进一步提高震相到时拾取的精度。
1. 3种到时自动拾取方法
1.1 AR-AIC方法
AR-AIC方法（Sleeman，van Eck，1999；王海军等，2003）识别震相的基础是假设地震波可以分成若干段平稳过程，每个过程均满足自回归模型。假定一个包含信号的时间序列{xi} （i＝1，2，…，N），已知大概的初至时间，该序列的前若干秒为噪声段，最后若干秒为信号段，分别用M阶自回归模型来拟合这两段，示意图如图1所示。自回归模型的系数可由Yule-Walker方法、Burgs算法或最小二乘法计算。假设自变量k对应于噪声与信号的分离点，利用噪声自回归模型拟合k点之前的序列，利用信号自回归模型拟合k点之后的序列，得到拟合误差，最终的AIC函数为
![]()
图 1 AR-AIC模型示意图N为序列的长度，M为拟合的模型阶数，k为噪声与信号序列的分离点Figure 1. The schematic diagram of AR-AIC modelN is the length of the sequence，M is the order of the AR-AIC model，and k is the divided point of noise and signal parts${\rm{AIC}}(k) \text{＝} (k \text{－} 1){\rm{lg}} \; \sigma _1^2 \text{＋} (N \text{－} k \text{＋} 1){\rm{lg}} \; \sigma _2^2\text{，}$
(1) 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.
$w(k) \text{＝} \frac{{(k \text{－} 1)\sum\limits_{i \text{＝} 1}^{k \text{－}1} {\left| {x(i)} \right|} }}{{(N - k \text{＋} 1)\sum\limits_{i \text{＝} k}^N {\left| {x(i)} \right|} }}$
(2) $f'(k) \text{＝} \frac{{f(k) \text{－} {f_{\min }}}}{{{f_{\max }} \text{－} {f_{\min }}}}$
(3) $f''(k) \text{＝} \text{±} \frac{{\left| {[f'(k) \text{－} f'(1)]\displaystyle\frac{{N \text{－} 1}}{N} \text{－} [f'(N) \text{－} f'(1)]\displaystyle\frac{k}{N}} \right|}}{{\sqrt {\displaystyle\frac{{{{(N \text{－} 1)}^2}}}{{{N^2}}} \text{＋} {{[f'(N) - f'(1)]}^2}} }}$
(4) 本文采用改进后的AR-AIC方法（王海军等，2003）拾取震相到时。图2a给出了新疆地震台网清河台2012年1月1日某地震事件的垂向记录，图2b为1—2 Hz滤波后的波形，图2c为使用改进AR-AIC方法拾取到时的结果。
![]()
图 2 实际地震的垂向记录及基于3种方法的到时拾取结果（a） 实际地震信号的垂向记录；（b） 1—2 Hz滤波后的地震波形；（c） AR-AIC方法自动拾取的到时；（d） HOS方法自动拾取的到时；（e） CUSUM方法自动拾取的到时Figure 2. The vertical component record of a real earthquake and the results of arrival time picked by the three estimation methods（a） The vertical component record of a real earthquake；（b） The seismic waveform after 1—2 Hz bandpass filtering；（c） The estimated result with AR-AIC method；（d） The estimated result with HOS method；（e） The estimated result with CUSUM method1.2 HOS方法
HOS方法（Küperkoch et al，2010 ）是通过计算地震波形的偏斜度和峰度来判断信号初动，由于噪声为平稳过程而信号为非平稳过程，因此当有信号到达时，统计量会发生明显的改变。偏斜度S和峰度K分别定义为
$S \text{＝} \frac{{E[{{(X \text{－} E[X])}^3}]}}{{E{{[{{(X \text{－} E[X])}^2}]}^{{3}/{2}}}}}{\rm{ \text{＝} }}\frac{{{m_3}}}{{{m_2^{3/2}}}}\text{，}$
(5) $K \text{＝} \frac{{E[{{(X \text{－} E[X])}^4}]}}{{E{{[{{(X \text{－} E[X])}^2}]}^{4/2}}}} \text{＝} \frac{{{m_4}}}{{m_2^2}}\text{，}$
(6) 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.
${\hat m_k}(j) \text{＝} \frac{1}{N}\sum\limits_{l \text{＝} 0}^{N \text{－} 1} {x_{j \text{－} l}^k},$
(7) 递推计算为
${\hat m_k}(j) \text{＝} {\hat m_k}(j \text{－} 1) \text{－} \frac{{{x^k}(j \text{－} N) \text{－} {x^k}(j)}}{N}.$
(8) 根据式（5）和式（6）得到偏斜度和峰度，令偏斜度或峰度超过某一阈值的位置为信号到时；但是为了更准确地拾取到时，Küperkoch等（2010）以偏斜度或峰度为特征函数，再加入AIC方法的思想，有
${\rm{AIC}}(k) \text{＝} (k \text{－} 1){\rm{lg}}\left({\frac{1}{{k \text{－} 1}}\sum\limits_{j \text{＝} 1}^{k \text{－} 1} {F_j^2} } \right) \text{＋} (L \text{－} k \text{＋} 1){\rm{lg}}\left({\frac{1}{{L \text{－} k \text{＋} 1}}\sum\limits_{j \text{＝} k}^L {F_j^2} } \right)\text{，}$
(9) 式中F为特征函数，L为特征函数的长度，0＜k＜L，P波到时对应AIC函数最大值前的局部最小值。本文选择峰度作为特征函数，计算kur-AIC曲线。但是在实际处理中发现，kur-AIC曲线最大值之前可能不存在局部最小值，因此本文还加入另一种拾取思路，即令峰度超过某一阈值的点对应震相初至，例如阈值为αmax（K′），K′为峰度的一阶导数，α为比例系数。HOS方法拾取到时的结果如图2d所示。
1.3 CUSUM方法
累积和（CUSUM）方法（Inclán，Tiao，1994）是一种序贯分析方法，可以检测出开始发生异常的数据点位置，如整个数列的平均值或者均方差开始发生改变的位置。假定{xi}为一段包含信号的记录，长度为N，则归一化的CUSUM曲线可定义为
${D_k} \text{＝} \frac{{{C_k}}}{{{C_N}}} \text{－} \frac{k}{N}, \quad \quad 0 < k < N\text{，}$
(10) 式中
${C_k} \text{＝} \sum\limits_{i \text{＝} 1}^k {x_i^2} .$
(11) 另外，假定噪声和信号的幅值均服从均值为0的正态分布，则可以求得k点处有信号和无任何信号两种假设之间的似然函数比为
$\lg L(k) \text{＝} \text{－} \frac{k}{2}\lg \left({1 \text{＋} \frac{N}{k}{D_k}} \right) \text{－} \frac{{N \text{－} k}}{k}\lg \left({1 \text{－} \frac{N}{{N \text{－} k}}{D_k}} \right)\text{，}$
(12) 此时，震相到时对应于归一化曲线和似然比曲线的极小值点，该方法自动拾取到时的结果如图2e所示。
2. 到时拾取方法的参数优化
2.1 参数变化的影响
AR-AIC，HOS和CUSUM这3种到时拾取方法所需要确定的参数不尽相同（表1），但均包含拾取时窗长度参数项，这是因为：进行到时拾取时，需要首先在地震记录上用STA/LTA方法检测到信号到达的大致位置，再在以此为中心的时窗内用上述方法进一步估算震相到时。参数不同会影响到时估计的结果，尤其是对于低信噪比P波信号而言。为了说明参数的影响，选择新疆地震台网富蕴台2015年10月1日一段包含远震P波的记录，以叠加地震噪声的方式改变该段数据的信噪比，然后使用3种自动拾取方法确定震相初至时间，参数取值如表2所示。信噪比的计算方法是取人工到时前10 s、到时后5 s地震波形的噪声峰峰值npp与信号峰峰值spp之比，即
表 1 3种到时拾取方法中待确定的参数Table 1. The parameters used in the three methods for determining arrival times方法 待确定参数 AR-AIC方法 拾取时窗长度T，噪声段长度Tn，信号段长度Ts，模型的阶数M HOS方法 拾取时窗长度T，滑移窗长度W，一阶导数门限系数α CUSUM方法 拾取时窗长度T ${\rm{SNR}} \text{＝} \frac{{{s_{{\rm{pp}}}}}}{{{n_{{\rm{pp}}}}}}.$
(13) 表3给出了当信噪比降低时，用表2所示参数取值组合得到的自动到时与人工到时的误差，可见：无论信噪比取高值还是低值，AR-AIC方法和CUSUM方法的误差均小于HOS方法的误差，说明AR-AIC方法和CUSUM的拾取精度高于HOS方法；当信噪比降低时，3种方法的到时误差均增大，说明对于低信噪比信号，到时拾取的精度较低；当信噪比大于19时，同一方法到时结果的最大、最小值之差基本处于0.03 s内，而当信噪比小于5时，同一方法到时结果的最大、最小值之差大于0.10 s，有些甚至大于0.50 s，说明对于低信噪比信号，参数变化对到时结果的影响较大。所以为了更准确地估计低信噪比P波信号的到时，需要确定每种到时拾取方法的最优参数。
表 2 用于合成地震波的3种到时拾取方法的参数取值组合Table 2. The combination of parameter selection of the three methods for arrival times picking in synthesizing seismic wave序号 AR-AIC方法 HOS方法 CUSUM方法 T Tn Ts M T W α T 组合① 20 2 3 17 30 9 0.36 22 组合② 20 3 4 17 26 7 0.36 12 组合③ 30 3 4 8 24 4 0.24 32 另外，HOS方法的到时均晚于AR-AIC方法和CUSUM方法的到时，这与HOS方法的原理有关，该方法是以峰度开始增加的位置为信号到时，而计算峰度的滑移窗有一定的平滑作用，使得峰度有明显增加的位置滞后于信号初动位置，尤其是对于初动较为平缓的信号。图3和图4分别给出了表3中SNR为31.1和4.5时的波形拾取到时过程，可以看出对于低信噪比信号，其峰度曲线中峰度增加的位置晚于信号初动位置。
表 3 不同参数取值时自动到时与人工到时的误差Table 3. The difference between auto-onset times and manual onset times with different parameter combination selectionSNR AR-AIC HOS CUSUM 组合① 组合② 组合③ 组合① 组合② 组合③ 组合① 组合② 组合③ 31.1 −0.03 −0.02 0 0.12 0.09 0.11 0.03 0.05 0.03 19.3 −0.02 −0.05 0.06 0.20 0.19 0.18 0.06 0.08 0.07 9.6 0.09 0.09 0.10 0.33 0.43 0.28 0.13 0.15 −0.34 4.5 0.16 0.22 0.31 0.46 0.57 0.56 −0.39 0.15 −0.39 2.4 0.19 0.18 −0.39 0.44 0.66 0.63 −0.39 0.18 −0.39 ![]()
图 3 信噪比为31.1时基于HOS方法的P波信号拾取到时过程（a） 0.75—1.5 Hz滤波后的地震波形；（b） 峰度曲线；（c） 峰度-AIC曲线Figure 3. The process of onset time estimation using HOS method for P signal with SNR 31.1（a) The seismic wave after 0.75-1.5 Hz bandpass filtering；（b） The kurtosis curve；（c） The kur-AIC curve![]()
图 4 信噪比为4.5时基于HOS方法的P波信号拾取到时过程（a） 0.75—1.5 Hz滤波后的地震波形；（b） 峰度曲线；（c） 峰度-AIC曲线Figure 4. The process of onset time estimation using HOS method for P signal with SNR 4.5（a) The seismic wave after 0.75-1.5 Hz bandpass filtering；（b） The kurtosis curve；（c） The kur-AIC curve2.2 到时拾取方法的参数确定
本文选择新疆地区台站2012年期间的波形记录，共450个远震P波震相（震中距大于30°），mb介于3.5—5.0，以人工拾取的到时作为参考样本。样本数据中，SNR＜5.0的P波震相210个，5≤SNR＜10的166个，10≤SNR＜15的46个，15≤SNR＜20的12个，20≤SNR＜30的11个，SNR≥30的5个。由于特别针对低信噪比信号，所以样本集中低信噪比信号所占的比例最高。
对样本数据，同样先用STA/LTA方法检测信号，再在拾取时窗内采用3种方法拾取震相到时，由于前两种方法中待确定的参数不只一种，因此以网格搜索方式逐步缩小参数取值范围。拾取时窗的长度是首要确定的因素，AR-AIC方法中T在12—40范围内变化，Tn和Ts的取值均小于T，M从1变化至20，逐步缩小参数取值范围，所用参数组合共280种；HOS方法中T从12变化至50，W小于T，α从0.2变化至0.5，所用参数组合共131种；CUSUM方法中T从6变化至40，参数组合共18种。
拟合每次搜索中自动拾取结果与人工到时的误差，目标函数为样本集误差标准差的最小值。经网格搜索后可知AR-AIC方法中T=20，Tn=2，Ts=3，M=17时，HOS方法中T=16，W=4，α=0.38时，CUSUM方法中T=22时，样本集的误差标准差最小，此即为样本集P波信号的最优参数组合。用此最优参数组合拾取样本集的信号到时，其误差分布直方图如图5所示，可以看到AR-AIC方法的到时误差均值和标准差最小，CUSUM方法其次，HOS方法最大，初步说明对于远震P波到时的拾取，AR-AIC方法的稳定性最优，CUSUM方法次之，HOS方法较差。
![]()
图 5 基于AR-AIC方法（a），HOS方法（b）和CUSUM方法（c）的自动拾取到时与人工到时误差分布及其比较（d）Figure 5. Error between auto-onset times and manual onset times based on the methods AR-AIC （a），HOS （b） and CUSUM （c） and their comparison （d）3. 3种到时拾取方法的比较
本文涉及的3种到时拾取方法均能比较有效地拾取P波初至时间，对于高信噪比信号其拾取结果均较准确，但对低信噪比信号拾取精度各有不同。为了进一步说明3种到时拾取方法对低信噪比P波信号的拾取精度，重新选择信噪比在[2，20]区间之内的100个远震P波记录，同样以人工拾取的到时为参考样本，用得到最优参数的3种方法估计P波初至时间。表4给出了不同误差区间的统计结果，图6则给出了自动拾取结果与人工拾取结果的误差和信噪比的关系图。
![]()
图 6 3种方法拾取到时的误差与信噪比的关系图Figure 6. The relationship between errors in onset time estimation and SNR by using the three methods表 4 3种方法所得到的到时拾取误差统计Table 4. Error statistics of phase onset time estimation by using the three methods到时拾取方法 自动拾取与人工拾取的误差处于不同区间所占比例 误差均值/s 误差标准差/s 误差<0.1 s 误差<0.3 s 误差<0.5 s 误差<1 s AR-AIC方法 36% 79% 93% 100% 0.19 0.15 HOS方法 13% 59% 87% 97% 0.31 0.28 CUSUM方法 28% 82% 95% 100% 0.20 0.17 表4说明，AR-AIC方法和CUSUM方法的到时误差均在1 s 以内，二者的误差均值与误差标准差相当，而HOS方法有3个P震相的到时误差大于1 s，并且其误差均值和标准差最大；在误差小于0.1 s的震相个数方面，AR-AIC方法所占比例最高，误差小于0.3 s和0.5 s方面，AR-AIC方法所占比例略低于CUSUM方法。从图6则可看到：信噪比小于6时，3种方法的到时误差较大，且HOS方法中有3个信号的到时误差大于1 s；当信噪比大于10时，3种方法的到时误差均小于0.5 s。
上述结果说明，3种方法中：AR-AIC方法和CUSUM方法的拾取精度相当；HOS方法的拾取能力最差，特别是对于初动不够尖锐的信号，其到时误差可能更大；CUSUM方法的计算过程最为简单，所以其运行速度最快。
4. 讨论与结论
本文简要介绍了AR-AIC，HOS，CUSUM这3种到时拾取方法的基本原理，分析了变化的参数对到时结果的影响，然后利用新疆地区台站的地震数据确定每种方法的最优参数，并且针对低信噪比P波信号对3种方法进行评估，得到结论如下：
1） 3种方法均能有效地估计P波初至时间，信号的信噪比较高时，拾取结果比较准确，而信噪比较低时，到时拾取可能有较大误差，且3种方法中参数值的改变会影响到时估计的结果，尤其是对于微弱信号。
2） AR-AIC方法拾取到时的能力比较稳定，对不同信噪比P波信号的到时，其估计精度均较高；HOS方法适用于处理信噪比较高且初动尖锐的信号，对于初动平缓的信号，其估计的到时系统滞后；CUSUM方法计算简单，运行速度最快，拾取到时的精度与AR-AIC方法相当，处理大批量数据时可以考虑该方法。
本文对于参数优化仅考虑了远震P波震相，并且事件数有限，所确定的最优参数存在局限性。此外，方法评估时，样本集的震相有所筛选，未能包含实际处理中出现的各种情况，不足之处，有待于进一步改进。
感谢新疆地震局为本文提供台站波形数据。
-
![]()
图 1 AR-AIC模型示意图
N为序列的长度，M为拟合的模型阶数，k为噪声与信号序列的分离点
Figure 1. The schematic diagram of AR-AIC model
N is the length of the sequence，M is the order of the AR-AIC model，and k is the divided point of noise and signal parts
![]()
图 2 实际地震的垂向记录及基于3种方法的到时拾取结果
（a） 实际地震信号的垂向记录；（b） 1—2 Hz滤波后的地震波形；（c） AR-AIC方法自动拾取的到时；（d） HOS方法自动拾取的到时；（e） CUSUM方法自动拾取的到时
Figure 2. The vertical component record of a real earthquake and the results of arrival time picked by the three estimation methods
（a） The vertical component record of a real earthquake；（b） The seismic waveform after 1—2 Hz bandpass filtering；（c） The estimated result with AR-AIC method；（d） The estimated result with HOS method；（e） The estimated result with CUSUM method
![]()
图 3 信噪比为31.1时基于HOS方法的P波信号拾取到时过程
（a） 0.75—1.5 Hz滤波后的地震波形；（b） 峰度曲线；（c） 峰度-AIC曲线
Figure 3. The process of onset time estimation using HOS method for P signal with SNR 31.1
（a) The seismic wave after 0.75-1.5 Hz bandpass filtering；（b） The kurtosis curve；（c） The kur-AIC curve
![]()
图 4 信噪比为4.5时基于HOS方法的P波信号拾取到时过程
（a） 0.75—1.5 Hz滤波后的地震波形；（b） 峰度曲线；（c） 峰度-AIC曲线
Figure 4. The process of onset time estimation using HOS method for P signal with SNR 4.5
（a) The seismic wave after 0.75-1.5 Hz bandpass filtering；（b） The kurtosis curve；（c） The kur-AIC curve
![]()
图 5 基于AR-AIC方法（a），HOS方法（b）和CUSUM方法（c）的自动拾取到时与人工到时误差分布及其比较（d）
Figure 5. Error between auto-onset times and manual onset times based on the methods AR-AIC （a），HOS （b） and CUSUM （c） and their comparison （d）
![]()
图 6 3种方法拾取到时的误差与信噪比的关系图
Figure 6. The relationship between errors in onset time estimation and SNR by using the three methods
表 1 3种到时拾取方法中待确定的参数
Table 1 The parameters used in the three methods for determining arrival times
方法 待确定参数 AR-AIC方法 拾取时窗长度T，噪声段长度Tn，信号段长度Ts，模型的阶数M HOS方法 拾取时窗长度T，滑移窗长度W，一阶导数门限系数α CUSUM方法 拾取时窗长度T 
下载: 导出CSV
表 2 用于合成地震波的3种到时拾取方法的参数取值组合
Table 2 The combination of parameter selection of the three methods for arrival times picking in synthesizing seismic wave
序号 AR-AIC方法 HOS方法 CUSUM方法 T Tn Ts M T W α T 组合① 20 2 3 17 30 9 0.36 22 组合② 20 3 4 17 26 7 0.36 12 组合③ 30 3 4 8 24 4 0.24 32
下载: 导出CSV
表 3 不同参数取值时自动到时与人工到时的误差
Table 3 The difference between auto-onset times and manual onset times with different parameter combination selection
SNR AR-AIC HOS CUSUM 组合① 组合② 组合③ 组合① 组合② 组合③ 组合① 组合② 组合③ 31.1 −0.03 −0.02 0 0.12 0.09 0.11 0.03 0.05 0.03 19.3 −0.02 −0.05 0.06 0.20 0.19 0.18 0.06 0.08 0.07 9.6 0.09 0.09 0.10 0.33 0.43 0.28 0.13 0.15 −0.34 4.5 0.16 0.22 0.31 0.46 0.57 0.56 −0.39 0.15 −0.39 2.4 0.19 0.18 −0.39 0.44 0.66 0.63 −0.39 0.18 −0.39
下载: 导出CSV
表 4 3种方法所得到的到时拾取误差统计
Table 4 Error statistics of phase onset time estimation by using the three methods
到时拾取方法 自动拾取与人工拾取的误差处于不同区间所占比例 误差均值/s 误差标准差/s 误差<0.1 s 误差<0.3 s 误差<0.5 s 误差<1 s AR-AIC方法 36% 79% 93% 100% 0.19 0.15 HOS方法 13% 59% 87% 97% 0.31 0.28 CUSUM方法 28% 82% 95% 100% 0.20 0.17
下载: 导出CSV
-
何燕,靳平,肖卫国,王红春. 2014. 迭代累积平方和算法与自回归赤池准则算法在地震信号到时估算中的比较研究[J]. 地震学报,36(2):209–219. He Y,Jin P,Xiao W G,Wang H C. 2014. A comparative study between ICSS and AR-AIC algorithms on onset time estimation for seismic signals[J]. Acta Seismologica Sinica,36(2):209–219 (in Chinese)
靳平,张诚鎏,沈旭峰,王红春,潘常周,严峰,王电源. 2014. 基于信号整体与局部特征的地震数据自动处理方法研究[J]. 地震学报,36(3):464–479 Jin P,Zhang C L,Shen X F,Wang H C,Pan C Z,Yan F,Wang D Y. 2014. A novel technique for automatic seismic data processing using both integral and local features of seismograms[J]. Acta Seismologica Sinica,36(3):464–479 (in Chinese)
林凡生, 邹志辉. 2016. 基于STA/LTA与多道互相关结合的地震波初至自动拾取方法[G]//国家安全地球物理丛书（十二）: 地球物理与信息感知. 西安: 西安地图出版社: 74–80. Lin F, Zou Z H. 2016. Automatic seismic first-arrival picking method based on the combination of STA/LTA and multi-channel cross correlation[G]//National Security Geophysics Series （12）: Geophysics and Information Perception. Xi’an: Xi’an Cartogra-phic Publishing House: 74–80 （in Chinese）.
盛冠群,李振春,王维波,高澜. 2015. 基于小波分解与高阶统计量的微地震初至拾取方法研究[J]. 石油物探,54(4):388–395 Sheng G Q,Li Z C,Wang W B,Gao L. 2015. A new automatic detection method of microseismic events based on wavelet decomposition and high-order statistics[J]. Geophysical Prospecting for Petroleum,54(4):388–395 (in Chinese)
谭玉阳,于静,冯刚,何川. 2016. 微地震事件初至拾取SLPEA算法[J]. 地球物理学报,59(1):185–196 Tan Y Y,Yu J,Feng G,He C. 2016. Arrival picking of microseismic events using the SLPEA algorithm[J]. Chinese Journal of Geophysics,59(1):185–196 (in Chinese)
田优平,赵爱华. 2016. 基于小波包和峰度赤池信息量准则的P波震相自动识别方法[J]. 地震学报,38(1):71–85 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)
王海军,靳平,刘贵忠,王晓明. 2003. 区域震相初至估计[J]. 西北地震学报,25(4):298–303 Wang H J,Jin P,Liu G Z,Wang X M. 2003. Accurate estimation for arrival time of seismic wave[J]. Northwestern Seismologi-cal Journal,25(4):298–303 (in Chinese)
岳龙,刘怀山,刘凯,徐秀刚,邢磊. 2015. 基于时频分析的初至拾取方法研究[J]. 石油物探,54(5):508–520 Yue L,Liu H S,Liu K,Xu X G,Xing L. 2015. First-break picking based on time-frequency analysis[J]. Geophysical Prospecting for Petroleum,54(5):508–520 (in Chinese)
曾富英,李敏锋,申维. 2002. 地震波初至拾取的分形研究[J]. 现代地质,16(2):209–213 Zeng F Y,Li M F,Shen W. 2002. The fractal study on detecting arrival time of the first break[J]. Geoscience,16(2):209–213 (in Chinese)
周银兴. 2009. 微震事件检测及震相自动识别研究[D]. 北京: 中国地震局地震预测研究所: 49–51. Zhou Y X. 2009. Research on the Micro-Earthquake Detection and Seismic Phase Automatic Identification[D]. Beijing: Institute of Earthquake Prediction, China Earthquake Administration: 49–51 （in Chinese）.
Akaike H. 1973. Information theory and an extension of the maximum likelihood principle[C]//Proceedings of the 2nd International Symposium on Information Theory. Budapest: Akadémiai Kiadó: 267–281.
Allen R V. 1978. Automatic earthquake recognition and timing from single traces[J]. Bull Seismol Soc Am,68(5):1521–1532
Allen R V. 1982. Automatic phase pickers:Their present use and future prospects[J]. Bull Seismol Soc Am,72(6B):S225–S242
Baer M,Kradolfer U. 1987. An automatic phase picker for local and teleseismic events[J]. Bull Seismol Soc Am,77(4):1437–1445
Der Z A,Shumway R H. 1999. Phase onset time estimation at regional distances using the CUSUM algorithm[J]. Phys Earth Planet Inter,113(1/4):227–246
Der Z A,Shumway R H. 2003. Automatic interpretation of regional short period seismic signals using CUSUM-SA algorithms[J]. Lecture Notes in Earth Sciences,98:41–59
Greg W B, David J B, Jerry A C. 2001. IDC Processing of Seismic, Hydroacoustic, and Infrasonic Data[R]. Reston, Virginia: Science Applications International Corporation Research, SAIC-99/3023: 22–26.
Inclán C,Tiao G C. 1994. Use of cumulative sums of squares for retrospective detection of changes of variance[J]. J Am Stat Associat,89(427):913–923
Jin P,Zhang C L,Shen X F,Wang H C,Pan C Z,Lu N,Xu X. 2014. A novel technique for automatic seismic data processing using both integral and local feature of seismograms[J]. Earthquake Science,27(3):337–349
Jin P. 2015. Toward high-confidence and full automation of seismic data processing for CTBT monitoring[C/OL]//Science and Technology 2015 Conference. [2015-12-18]. http://www.ctbto.org/fileadmin/user_upload/SnT2015/SnT2015_Orals/T3.3-O8.pdf.
Jin P,Pan C Z,Zhang C L,Shen X F,Wang H C,Lu N. 2015. A novel progressive signal association algorithm for detecting teleseismic/network-outside events using regional seismic networks[J]. Geophys J Int,201(3):1950–1960
Jurkevics A. 1988. Polarization analysis of three-component array data[J]. Bull Seismol Soc Am,78(5):1725–1743
Kedrov O K,Ovtchinnikov V M. 1990. An on-line analysis system for three-component seismic data:Method and preliminary results[J]. Bull Seismol Soc Am,80(6B):2053–2071
Küperkoch L,Meier T,Lee J,Friederich W,EGELADOS Working Group. 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics[J]. Geophys J Int,181(2):1159–1170
Küperkoch L, Meier T, Diehl T. 2012. Chapter 16: Automated event and phase identification [G]//Bormann P ed. New Manual of Seismological Observatory Practice 2 （NMSOP-2）. Potsdam: IASPEI, GFZ German Research Centre for Geosciences: 17–23.
Maeda N. 1985. A method for reading and checking phase time in auto-processing system of seismic wave data[J]. Zisin,38(3):365–380
Nippress S E J,Rietbrock A,Heath A E. 2010. Optimized automatic pickers:Application to the ANCORP data set[J]. Geophys J Int,181(2):911–925
Saragiotis C D,Hadjileontiadis L J,Panas S M. 2002. PAI-S/K:A robust automatic seismic P phase arrival identification scheme[J]. IEEE Trans Geosci Remote Sens,40(6):1395–1404
Sleeman R,van Eck T. 1999. Robust automatic P-phase picking:An on-line implementation in the analysis of broadband seismogram recordings[J]. Phys Earth Planet Inter,113(1/4):265–275
-
期刊类型引用(3)
1. 王树旺，王文才，杨晓鹏，安昭，张卫东，石文兵，李亮，张蓉. 2023年甘肃白银平川区4.9级地震强地面运动特征分析. 地震工程学报. 2025(01): 207-216 . 
百度学术
2. 周守军，张国正，刘臣，林旭东，曹长红. 基于CUSUM方法的集中供热管网泄漏故障诊断. 山东建筑大学学报. 2023(05): 48-57 . 百度学术
3. 赵震华，张杏莉，卢新明. 微震信号初至拾取的AIC算法及其分析. 山东科技大学学报(自然科学版). 2022(01): 44-53 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 1181
- HTML全文浏览量: 520
- PDF下载量: 62
- 被引次数: 6
