2016, Vol. 59
地震P波、S波初至时间是精确计算震源位置的必要参数，因此地震P波、S波初至时间的拾取是地震数据处理的一项必要工作(周宝峰，2012).在地震发生后，大量地震记录需要快速处理来实现P波、S波初至时间的精确和快速拾取，特别是在地震预警方面，地震P波、S波初至时间的拾取速度尤其重要.因此，基于传统的人工拾取P波、S波初至时间的方法和手段已经难以满足社会发展需求，促使 相关科研人员研究地震P波、S波初至时间的自动 拾取技术(刘劲松等，2013；叶根喜等，2008；马强，2008).
地震P波、S波初至时间的自动拾取算法一直在不断改进，迄今为止已有多种较成熟的方法，如相关法(Bai and Kennett，2000；Allen，1978)、能量比 法(Allen，1982；Maeda，1985；王继等，2006；Saragiotis et al.，2002)、 最大振幅法、分形维法及神经网络法(刘伊克等，2001；Paige and Saunders，1982；Boschetti et al.，1996；Chang et al.，1999；Cao and Greenhalgh，1993)等.目前最常用的方法主要有以下几种(马强，2008)：
(1) 长短时间平均(STA/LTA)方法：Allen(1978，1982)提出的STA/LTA方法，是目前广泛使用的一种地震P波、S波初至时间的自动拾取方法，其主要原理是根据地震波形特征函数的长短时均比值等特征拾取初至.
(2) AIC方法：该方法的基本原理是求取地震信号AIC函数局部最小值的位置(Maeda，1985；王继等，2006).
(3) 基于高阶统计量偏斜度和峰度的PAI-S/K方法：Saragiotis等(2002)提出了基于地震波形峰度和偏斜度的拾取初至方法，称为PAI-S/K方法.
但是由于不同地震记录的特点存在较大差异，目前存在的地震P波、S波初至时间的自动拾取方法没有哪一种可以拾取出所有地震记录的P波、S波初至时间，即单种方法不能保证所有地震记录的P波、S波初至时间的100%的准确拾取(叶根喜等，2008；马强，2008).因此，为了提高地震P波、S波初至时间的自动拾取结果的可靠性，马强(2008)综合多种自动拾取方法的特点，提出了地震P波、S波初至时间的多步自动拾取法.
相比于地震P波初至时间的拾取，地震S波初至时间的拾取更难，也是地震初至时间自动拾取技术的主要研究内容和研究难点(刘劲松等，2013；叶根喜等，2008).P波初至波先于S波初至波到达，但是由于S波初至波常在P波未完时就到达，使S波初至波叠加在P波中而难以被区分.这也是基于现有的地震P波、S波初至时间的自动拾取方法难以 精确拾取S波初至时间的主要原因(Bai and Kennett，2000；Allen，1978). 因此，如何从P波、S波混叠部分识别S波初至波的到达时刻，是地震波初至时间的自动拾取技术的研究热点(刘劲松等，2013；叶根喜等，2008；马强，2008).
本文提出一种新的自动拾取P波和S波初至时间的方法，首先把地震波的三分量时程曲线变换为一组空间向的能量变化率时程曲线，然后对能量变化率时程曲线进行STA/LTA处理，自动拾取地震波的P波和S波大致初至时间，最后建立了一种二次方自回归模型并基于此模型对地震P波、S波初至波附近的能量变化率曲线进行二次方自回归处理，获得P波、S波初至波的二次方自回归时程曲线，最后自动拾取P波和S波的初至波的二次方自回归时程曲线的最大值所对应的时间，所拾取的两个时间值即为P波、S波的精确初至时间.
2 基于信号能量变化率的二次方自回归模型拾取初至原理介绍对某个地震信号x(t)进行等间隔Δt离散化后得序列{xi}，假设序列长度为N，可定义此信号的震动瞬态幅值变化差为
|
(1) |
则{pi}为一个等间隔时间为Δt的离散序列，它在物理意义上表示物体震动幅值的变化量值.
同时定义
|
(2) |
则序列 Ei 在物理意义上表示物体瞬态震动能量变化率.
假设由三分量传感器所记录到的某个地震信号的两个水平向和竖向震动可分别表示为{xi}，{yi}，{zi}，0≤i<N.则按照公式(2)可推得
|
(3) |
则序列{Qi}描述了地震信号在三维空间内震动能量变化率的时程特征.
当地震到达时，地震记录波形在三维空间内的振幅会突然增大；同样在地震到达时刻，地震记录波形在三维空间内的震动能量变化率也会突然增大.因此可根据序列{Qi}的幅值时程曲线来分析出地震到达时刻.为了方便描述，本文定义序列{Qi}为地震记录波形的震动能量变化率谱.
为了能准确地自动找到地震信号到达时刻，本文提出一种多步拾取的新方法.首先基于STA/LTA算法对序列{Qi}进行STA/LTA变换，且基 于STA/LTA的峰值特性确定信号到达时刻的狭窄时间段 2sΔt(s为截取序列长度，通常取0.2≤sΔt≤0.6).序列{Qi}通过STA/LTA变换后得到序列{ei}，其变换公式如下：
|
(4) |
N4为STA/LTA的短尺度，N3为STA/LTA的长尺度.由公式(4)可见序列{ei}相对序列{Qi}的起始 时间延迟了(N3-N4+1)Δt.序列{ei}的长度为N2=N-2-N3.
令
|
(5) |
en所在位置对应地震记录波上的时间为tn=(N3-N4+1+n)Δt.tn为地震波的大致初至时间，可以此值为中心确定一个包含了真实初至时间的狭窄时间段为[tn-sΔt，tn+sΔt]，则此时间段对应序列{Qi}的子序列为[QN3-N4+1+n-s，QN3-N4+1+n+s].如果定义此子序列为{qi}，则可得
|
(6) |
假设地震信号的准确到达时刻为T1=T0+T2，T0为开始记录波形的绝对时间，T2为地震信号准确到达时间与开始记录波形时刻的相对时间，则可得T2∈[tn-sΔt，tn+sΔt].
为了进一步精确自动拾取到达时刻，本文提出第二步拾取法：对序列{qi}基于二次方自回归模型法，精确地自动拾取初至.
此模型的理想二次方程假设为：
|
(7) |
为了求得此二次方的参数，本文提出基于公式(8)得到其逼近序列{c′i}.
|
(8) |
然后基于最小二乘法求解系数w，g，f.令
|
(9) |
|
(10) |
由式(10)可得：
|
(11) |
由式(11)可解得系数w，g，f.然后基于式(7)求取序列{ci}最大值对应的i，即为地震信号精确初至时间，得到地震信号到达的绝对精确时刻为：
|
(12) |
首先基于青川某斜坡地震监测站的一组实测场地脉动和某个水平向余震信号进行仿真实验来验证本文提出的地震初至自动拾取新方法的可靠性.余震信号为地震P波已经到达后截取的某段地震波，采样频率为256 Hz.为了方便分析，本文人为地在场地脉动信号的3.4 s时刻叠加余震信号，得到信号S3.因此，由信号S3可以看出地震到达时刻为3.4 s，即P波到达时刻为3.4 s，对S2进行人工拾取S波到达时刻为2.6006 s，则S3里的S波到达时刻实际值为6.006 s.
由图 2可得：对S3的绝对幅值按照长的平均取1 s、短的平均取0.05 s的参数进行STA/LTA分析可得到P波的到达时间大概在3.4123 s附近，而S波的到达时间比较难确定，可初步确定在4.9436、5.3195、6.0456、6.2547这四个时刻的某个时刻的附近.对S3的震动能量变化率按照相同参数进行STA/LTA分析，可得到P波的到达时间大概在3.4032 s 附近，而S波的到达时间可确定在6.0456 s附近. 因此，基于能量变化率的STA/LTA处理和基于绝对幅值的STA/LTA处理都能比较准确地分析出P波的到达时间，但是基于能量变化率的STA/LTA处理更能准确分析出S波的到达时刻.
|
图 1 模拟地震信号(a表示振动加速度) Fig. 1 Simulated seismic signals |
|
图 2 S3的能量变化率和长短时间平均处理 Fig. 2 Transforming S3 into its energy gradients and STA/LTA curves |
为了更准确地确定P波和S波的到达时间，对S3的能量变化率的时程波形分别以3.4032、6.0456 s 所对应的数据为中心，截取长度为200个的两个序列，采用本文建立的二次方自回归模型法来进一步确定P波和S波的初至，即初步确定P波初至在(3.0116，3.7948)s内，S波初至在(5.6450，6.4372)s内. 两个序列的二次方自回归模型曲线如图 3所示.读取两条回归曲线最大幅值对应的时间分别为 0.3892 s、0.3720 s.因此，P波到达时刻为3.4008 s，S波到达时刻为6.0170 s.
|
图 3 S3的能量变化率的二次方的自回归曲线 Fig. 3 Quadratic auto-regressive curves of energy gradients of S3 |
基于以上实验仿真分析可初步得出以下结论：
(1) 绝对幅值的STA/LTA处理结果的最大峰值对应P波的初至，第2个最大峰值不一定对应S波的初至.因此，可拾取绝对幅值的STA/LTA处理结果的最大峰值来确定P波的初至，但是难以依据此方法来拾取到S波的初至时间.
(2) 能量变化率的STA/LTA处理结果的最大峰值对应P波的初至，第2个最大峰值一定对应S波的初至.因此，可拾取能量变化率的STA/LTA处理结果的最大峰值来确定P波的初至，可拾取能量变化率的STA/LTA处理结果的第2个最大峰值来确定S波的初至时间.
(3) 基于能量变化率的STA/LTA处理结果和拾取的大致初至时间，分别对P波和S波初至附近的小段能量变化率曲线进行二次方自回归处理，分析出P波和S波的自回归曲线，然后分别拾取自回归曲线的最大幅值来确定P波和S波的精确初至时间.
此仿真实验初步证明了本文提出的基于地震信号的能量变化率分别进行STA/LTA处理和二次方自回归模型分析，可准确拾取出P波和S波初至时间.其拾取流程如图 4所示.
|
图 4 基于能量变化率的初至自动拾取流程 Fig. 4 Flow of picking up the arrival times of P and S waves automatically |
为了进一步论证噪声强度对此方法的拾取结果精度的影响，本文对前面仿真分析的噪声S1分别按照不同比例进行放大，直到其最大值达到地震波S2最大值的50%为止，然后把地震波S2在3.4 s时刻加入到这些经过放大后的噪声信号S1里，最后编写相关程序，采用本文提出的方法和传统的AIC方法及STA/LTA方法，共计三种方法对叠加后的信号进行地震波初至时间自动拾取.三种方法的分析结果如表 1所示.
|
表 1 三种方法的自动拾取误差统计 Table 1 Errors of picking-up the arrival times with three different methods |
由表 1可见：三种方法对P波和S波初至时间的拾取误差随着噪声放大倍数的增加而逐步增大，如果以0.1 s为判断拾取失效的阈值，则STA/LTA 方法在放大倍数达到8倍时的拾取结果失效，AIC方法在放大倍数达到12倍时的拾取结果失效，新方法在放大倍数达到12倍时的P波初至时间拾取结果失效，但是S波初至时间拾取结果还未失效.在失效前，新方法和AIC方法的拾取误差都明显小于STA/LTA方法的拾取误差；新方法对P波的拾取误差与AIC方法的拾取误差很接近，但是对S波初至时间的拾取误差要小于AIC方法的拾取误差.
通过以上仿真实验表明：本文提出的新方法和传统的AIC方法，二者相比于STA/LTA方法，具有较强的抗噪声的性能，自动拾取P波和S波的初至时间的误差也更小些，但是新方法对S波初至时间的拾取精度要更优于AIC方法.
如果以三种方法的计算时间长度来衡量计算效率，则STA/LTA方法的运行效率要明显优于新方法和AIC方法，所需时间不到其他两种方法的1/4；而新方法与AIC方法相比虽然稍微慢些，但是相差不到5%.
4 实测地震信号的初至分析为了进一步检验基于能量变化率和二次方自回归模型法自动拾取P波和S波初至的可靠性，本文采用此新方法、STA/LTA方法和AIC方法，分别 对2013年4月20日发生的芦山地震的10组余震 记录数据进行自动拾取P波和S波的相对初至，同时与手工拾取结果进行对比.图 5和图 6为第1组余震记录数据的详细分析结果.由图 5可得：手工拾取的P波相对初至为1.1835s、S波相对初至为2.5076 s；自动拾取能量变化率的STA/LTA分析结果的前两个最大峰值得到P波相对初至为1.0663 s、S波相对初至为2.5702 s，且根据此拾取结果，截取(0.6751，1.4575)s和(2.1790，2.9614)s 时间段内的能量变化率曲线进行二次方自回归分析.
|
图 5 芦山地震的余震记录数据的初至分析 Fig. 5 Picking-up approximate arrival times of Lushan earthquake aftershocks |
由图 6可得：P波的二次方自回归曲线的最大值对应时刻为0.4531 s、S波的二次方自回归曲线的最大值对应时刻为0.3571 s.因此，综合两步自动拾取结果可得P波的精确相对初至为：
|
图 6 能量变化率的二次方自回归曲线 Fig. 6 Quadratic auto-regressive curves of energy gradient |
|
S波的精确相对初至为：
|
基于本文提出的方法对这10组余震记录数据的相对初至的自动拾取结果如表 2所示.
|
|
表 2 对10组余震记录数据初至的自动拾取结果的统计 Table 2 Comparison of picking up the arrival times by different methods on 10 groups seismic data |
由表 2可得：
(1) 如以手工拾取的结果为基准，则对10组地震记录单独采用STA/LTA法拾取的10个P波的初至时间都接近手工拾取结果；对S波的拾取结果里有9组的拾取结果接近手工拾取结果，但对第8组地震记录S波的拾取结果与手工拾取结果相差了 0.5424 s，存在较大偏差.
(2) 如以手工拾取的结果为基准，分别采用AIC方法和新方法对这10组地震记录自动拾取的初至时间拾取结果，无论是P波还是S波的初至时间的拾取结果，都接近手工拾取的结果，因此，这两种方法的拾取结果精度都优于STA/LTA方法.
(3) 如以手工拾取的结果为基准，则新方法和AIC方法对P波的初至时间的拾取结果误差基本一致，但是新方法对S波的初至时间的拾取结果误差要小于AIC方法.
因此，通过以上比较表明了本文提出的新的地震波初至时间多步拾取方法具有较高精度和可靠性，特别是在S波的初至时间拾取上，其精度要优于AIC方法.
为了进一步论证本文提出的地震P波、S波自动多步拾取的新方法的可靠性，本文选择了150条汶川地震记录并编写了相关程序分别采用新方法、STA/LTA方法和AIC方法，进行了P波、S波初至时间的自动拾取，同时进行了手工拾取.如以手工 拾取结果为基准，则三种方法的自动拾取结果情况如表 3所示.
|
|
表 3 三种拾取方法对150组地震记录的拾取结果误差统计 Table 3 Errors of picking up the arrival times with different methods |
由表 3可见：STA/LTA方法拾取精度比较低，新方法和AIC方法的拾取精度都比较高.如果以0.1 s的误差作为拾取结果失效的阈值，则AIC和新方法对P波的初至时间拾取准确率都为99.333%，对S波的初至时间的拾取准确率分别为95.333%和98%.STA/LTA方法的计算效率明显优于AIC 方法和新方法，而新方法在此方面稍逊于AIC方法.
因此，通过汶川地震的150条记录数据的P波、S波初至时间的拾取结果误差的统计，进一步检验了本文提出的多步自动拾取地震P波S波初至时间的新方法，虽然在计算效率方面稍逊于AIC方法，但是其对S波的初至时间的拾取精度和准确率要稍微高于AIC方法，是一种较可靠的地震波初至时间自动拾取方法.
5 结论本文首先总结了常见的几种地震P波和S波初至时间的自动拾取方法，指出了S波的初至时间自动拾取的难点，然后提出了一种新的地震P波和S波初至时间的自动拾取方法.此方法主要分三步来实现.
第1步：把三个方向的地震记录时程曲线合成 为一组反映地震能量突变情况的能量变化率时程曲线.
第2步：地震P波和S波初至时间的初步拾取.对能量变化率时程曲线进行STA/LTA处理，自动拾取能量变化率的STA/LTA处理结果的前两个最大峰值，然后把这两个最大峰值按照对应时间的先后进行排序.首先出现的峰值对应的时间为P波初至，后面出现的为S波初至.依据拾取结果确定第3步分析数据范围.
第3步：地震P波和S波初至时间的精确拾取.对由第2步拾取结果确定的两小段能量变化率时程曲线进行二次方自回归分析，分别得到地震P波和S波到达时的二次方自回归时程曲线，然后自动拾取两条二次自回归曲线最大值对应的时间，分别作为地震P波和S波的精确初至时间.
基于实测地震数据的分析表明，能量变化率时程曲线的STA/LTA处理结果具有典型的特征：前两个最大峰值分别对应地震P波和S波的初至时间.地震P波和S波到达时的能量变化率时程曲线的二次方自回归曲线也具有典型特征：曲线的最大值对应P波或S波的精确初至时间.
基于此新方法、AIC方法和STA/LTA方法，分别对芦山地震的10组地震记录和汶川地震的150组地震记录数据进行了初至时间的自动拾取，且同时进行了人工拾取获得准确的参考初至时间值.分析结果和误差统计表明，新方法和AIC方法的拾取精度和可靠性明显高于STA/LTA方法.新方法与AIC方法相比，虽然计算效率稍微差些，但是对S波的拾取精度和可靠性都稍微高些.
由于本文提出的方法的实现步骤是首先对地震波能量变化率曲线采用STA/LTA方法初步拾取出P波、S波的初至时间，然后再采用二次自回归法进一步精确拾取，因此理论上可以理解为是对现有STA/LTA法的补充和提升.
致谢感谢中国地震局工程力学研究所“国家强震台网中心”对本文提供的数据支持.
| Allen R V. 1978. Automatic earthquake recognition and timing from single traces. Bull. Seismol. Soc. Amer. , 68(5): 1521–1532. | |
| Allen R V. 1982. Automatic phase pickers: their present use and future prospects. Bull. Seismol. Soc. Amer. , 72(6B): S225–S242. | |
| Bai C, Kennett B L N. 2000. Automatic phase-detection and identification by full use of a single three-component broadband seismogram. Bull. Seismol. Soc. Amer. , 90(1): 187–198. | |
| Boschetti F, Dentith M D, List R D. 1996. A fractal-based algorithm for detecting first arrivals on seismic traces. Geophysics , 61(4): 1095–1102. | |
| Cao S H, Greenhalgh S. 1993. Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm. Geophysical Journal International , 114(3): 593–600. | |
| Chang X, Liu Y K, Ashida Y. 1999. Hausdorff fractal algorithm for picking first break in seismic traces. Butsuri Tansa , 52(4): 316–322. | |
| Liu J S, Wang Y, Yao Z X. 2013. On micro-seismic first arrival identification: A case study. Chinese J. Geophys. (in Chinese) , 56(5): 1660–1666. doi: 10.6038/cjg20130523. | |
| Liu Y K, Chang X, Wang H, et al. 2001. Estimation of near-surface velocity and seismic tomographic static corrections. Chinese J. Geophys. (in Chinese) , 44(2): 272–278. | |
| Ma Q. 2008. Study and application on earthquake early warning [Ph. D. thesis] (in Chinese). Harbin: Institute of Engineering Mechanics, China Earthquake Administration. | |
| Maeda N. 1985. A method for reading and checking phase time in auto-processing system of seismic wave data. Zisin-JIsin , 38(3): 365–379. | |
| Paige C C, Saunders M A. 1982. LSQR: Sparse linear equations and least squares problems. ACM Transactions on Mathematical Software , 8(2): 195–209. | |
| Saragiotis C D, Hadjileontiadis L J, Panas S M. 2002. PAI-S/K: A robust automatic seismic P phase arrival identification scheme. IEEE Transactions on Geoscience and Remote Sensing , 40(6): 1395–1404. | |
| Wang J, Chen J H, Liu Q Y, et al. 2006. Automatic onset phase picking for portable seismic array observation. Acta Seismologica Sinica (in Chinese) (in Chinese) , 28(1): 42–51. | |
| Ye G X, Jiang F X, Yang S H. 2008. Possibility of automatically picking first arrival of microseismic wave by energy eigenvalue method. Chinese J. Geophys. (in Chinese) , 51(5): 1574–1581. | |
| Zhou B F. 2012. Some key issues on the strong motion observation [Ph. D. thesis] (in Chinese). Harbin: Institute of Engineering Mechanics, China Earthquake Administration. | |
| 刘劲松, 王赟, 姚振兴. 2013. 微地震信号到时自动拾取方法. 地球物理学报 , 56(5): 1660–1666. | |
| 刘伊克, 常旭, 王辉, 等. 2001. 三维复杂地形近地表速度估算及地 震层析静校正. 地球物理学报 , 44(2): 272–278. | |
| 马强. 2008. 地震预警技术研究及应用[博士论文]. 哈尔滨: 中国地震局工程力学研究所. | |
| 王继, 陈九辉, 刘启元, 等. 2006. 流动地震台阵观测初至震相的自动检测. 地震学报 , 28(1): 42–51. | |
| 叶根喜, 姜福兴, 杨淑华. 2008. 时窗能量特征法拾取微地震波初始到时的可行性研究. 地球物理学报 , 51(5): 1574–1581. | |
| 周宝峰. 2012. 强震观测中的关键技术研究[博士论文]. 哈尔滨: 中国地震局工程力学研究所. | |