引言

随着现代通信手段和数字地震观测技术的发展，地震数据自动处理技术已经广泛应用于地震日常监测分析中，其中地震信号到时的自动拾取是自动处理过程中的重要一环。到时估算的准确性越高，事件关联和定位的可靠性也越高，所以研究地震信号到时的自动拾取技术具有十分重要的意义。

震相到时拾取方法主要依据地震记录中信号到达前后振幅、频率、偏振特性等特征的不同将信号从噪声中识别出来，最基础的震相到时拾取方法为长短时平均能量比（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）识别震相的基础是假设地震波可以分成若干段平稳过程，每个过程均满足自回归模型。假定一个包含信号的时间序列{x_i} （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 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

下载: 全尺寸图片 幻灯片

This page contains the following errors:

Below is a rendering of the page up to the first error.

本文采用改进后的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 method

下载: 全尺寸图片 幻灯片

1.2   HOS方法

HOS方法（Küperkoch et al，2010 ）是通过计算地震波形的偏斜度和峰度来判断信号初动，由于噪声为平稳过程而信号为非平稳过程，因此当有信号到达时，统计量会发生明显的改变。偏斜度S和峰度K分别定义为

This page contains the following errors:

Below is a rendering of the page up to the first error.

递推计算为

根据式（5）和式（6）得到偏斜度和峰度，令偏斜度或峰度超过某一阈值的位置为信号到时；但是为了更准确地拾取到时，Küperkoch等（2010）以偏斜度或峰度为特征函数，再加入AIC方法的思想，有

式中F为特征函数，L为特征函数的长度，0＜k＜L，P波到时对应AIC函数最大值前的局部最小值。本文选择峰度作为特征函数，计算kur-AIC曲线。但是在实际处理中发现，kur-AIC曲线最大值之前可能不存在局部最小值，因此本文还加入另一种拾取思路，即令峰度超过某一阈值的点对应震相初至，例如阈值为αmax（K′），K′为峰度的一阶导数，α为比例系数。HOS方法拾取到时的结果如图2d所示。

1.3   CUSUM方法

累积和（CUSUM）方法（Inclán，Tiao，1994）是一种序贯分析方法，可以检测出开始发生异常的数据点位置，如整个数列的平均值或者均方差开始发生改变的位置。假定{x_i}为一段包含信号的记录，长度为N，则归一化的CUSUM曲线可定义为

式中

另外，假定噪声和信号的幅值均服从均值为0的正态分布，则可以求得k点处有信号和无任何信号两种假设之间的似然函数比为

此时，震相到时对应于归一化曲线和似然比曲线的极小值点，该方法自动拾取到时的结果如图2e所示。

2.   到时拾取方法的参数优化

2.1   参数变化的影响

AR-AIC，HOS和CUSUM这3种到时拾取方法所需要确定的参数不尽相同（表1），但均包含拾取时窗长度参数项，这是因为：进行到时拾取时，需要首先在地震记录上用STA/LTA方法检测到信号到达的大致位置，再在以此为中心的时窗内用上述方法进一步估算震相到时。参数不同会影响到时估计的结果，尤其是对于低信噪比P波信号而言。为了说明参数的影响，选择新疆地震台网富蕴台2015年10月1日一段包含远震P波的记录，以叠加地震噪声的方式改变该段数据的信噪比，然后使用3种自动拾取方法确定震相初至时间，参数取值如表2所示。信噪比的计算方法是取人工到时前10 s、到时后5 s地震波形的噪声峰峰值n_pp与信号峰峰值s_pp之比，即

表  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 

| 显示表格

表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

下载: 导出CSV 

| 显示表格

另外，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 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 

| 显示表格

图  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

下载: 全尺寸图片 幻灯片

2.2   到时拾取方法的参数确定

本文选择新疆地区台站2012年期间的波形记录，共450个远震P波震相（震中距大于30°），m_b介于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范围内变化，T_n和T_s的取值均小于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，T_n=2，T_s=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

下载: 导出CSV 

| 显示表格

表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波震相，并且事件数有限，所确定的最优参数存在局限性。此外，方法评估时，样本集的震相有所筛选，未能包含实际处理中出现的各种情况，不足之处，有待于进一步改进。

感谢新疆地震局为本文提供台站波形数据。