2019, Vol. 62
2. 中国地震局地震观测与地球物理成像重点实验室, 北京 100081
2. Key Laboratory of Seismic Observation and Geophysical Imaging, China Earthquake Administration, Beijing 100081, China
随着地震监测台网的全面覆盖和日益密集，产生了海量的地震波形数据，极大增加了人工识别震相的工作量.准确的拾取震相到时，是诸如震源机制分析、层析成像、地震定位等研究分析工作的基础(Bai and Kennett, 2000；Bogiatzis and Ishii, 2015；刘劲松等，2013；王彩霞等，2013；何先龙等，2016).
地震震相自动拾取算法应运而生，依据各算法的特点可将其归纳为两类：一是依据震相的特征变换识别初至.通过波形振幅、能量、偏振、分形维数等特征变换判断震相到时，比较经典的算法有：长短时窗平均比(STA/LTA)算法(Stevenson，1976；Allen, 1978, 1982；Baer and Kradolfer, 1987；刘晗和张建中，2014)、AIC算法(Maeda，1985；St-Onge，2011；张唤兰等，2013)、高阶统计方法(Saragiotis et al., 2002；Küperkoch et al., 2010)等.二是综合判别分析.该类方法先通过某种算法粗略拾取震相到时，再结合其他方法得到精细初至到时.如Ross和Ben-Zion(2014)联合偏振分析、STA/LTA、峰度等3种方法进行多震相综合检测.
利用前述算法拾取震相到时过程中，自动拾取参数的选择尤为关键，当参数选取不恰当时往往会导致震相误拾、漏拾.对拾取参数进行优化选择可减少人为选择参数导致的震相不当拾取，为此Nippress等(2010)提出了基于参数优化的自动拾取算法，优化了三种传统算法(STA/LTA方法、卓越周期Tpd方法、峰度方法)的拾取参数.
另一方面，对拾取结果的可靠性进行评估是震相自动拾取过程中十分重要的环节.拾取结果的准确性是自动拾取算法优劣程度的重要评判依据之一.如Gentili和Bragato(2006)、Gentili和Michelini(2006)以信噪比为判据，评估了拾取结果的质量.在构建评估方案时，Nippress等(2010)主要以拾取点前后原始波形、特征函数的最大值之比为评判依据，但是当初动十分明显时会因幅值过大产生极值导致计算结果不稳定.
本文改进了Nippress等(2010)提出的参数优化算法和质量评估方案，引入了拾取到时优化方案，优化了频带-带宽拾取法、AICD拾取算法、峰度拾取算法的拾取参数，并将该方法应用到云南腾冲地区地震P、S波到时自动识别工作中.进一步以人工识别的初至为参照，检验了本文方法的可行性和可靠性.
1 自动拾取方法简介本文使用了Python软件包PhasePApy(Chen and Holland, 2016)中提供的3种自动拾取算法(频带-带宽拾取法、AICD拾取算法、峰度拾取算法)，由滑动窗口内特征函数(CF)的均方根值乘以阈值系数得到动态阈值，当CF超过阈值则判定为地震波初至信号.
1.1 频带-带宽拾取算法频带-带宽拾取法(FBpicker)是一种基于STA/LTA的改进算法，用最小二乘法去除非零均值和线性趋势，并采用滤波器组法对信号进行倍频程分析，在滤波后的序列上使用Cosine Tapered窗以减少Gibbs现象的影响.数据的瞬时能量为(Lomax et al., 2012；Chen and Holland, 2016)
|
(1) |
其中，BFn[i]是数据点i经过第n个滤波器后的幅值.
每一组滤波数据的CF为
|
(2) |
式中l是滑动窗口的长度, rms(En[i-1-l:i-1])是数据点i之前的滑动长窗口中能量的均方根值.取各组滤波数据中CFnrms[i]的最大值为该点的特征函数.CF体现了短时窗相对于长时窗的能量变化，值越高则表明短时窗中能量变化越大.
1.2 AICD拾取算法AR-AIC依据自回归过程将地震记录划分成背景噪声和地震信号两个局部模型，两个模型的分界点即为震相到时，也是AIC函数的极小值点(Leonard and Kennett, 1999).若P为分界点，地震记录x(i)(i = 1, 2, 3, …, N)的方差AIC表达式为(Maeda，1985)：
|
(3) |
其中，N为信号长度，var(x[1, P])和var(x[P+1, N])分别为分界点前后两个阶段的方差.因为AIC函数一阶后向差分(AICD)对原函数的变化十分敏感，因此将其作为AICD拾取算法(AICD picker)的CF.
1.3 峰度拾取算法四阶累积量的中心值称为“峰度”(田优平和赵爱华，2016).将数据的峰度作为CF，可得四阶统计量峰度的表达式(Baillard et al., 2014)：
|
(4) |
式中x是样本值，E是期望，μ是样本的平均值，μ4是四阶中心矩，σ是标准差.
波形中的背景噪声一般服从高斯分布(K=3)，发生地震后信号由高斯分布逐渐向非高斯分布过渡，地震记录的峰度值变高.峰度拾取算法(KTpicker)基于样本的分布形态进行判别，同样适用于信噪比较低的数据.
2 基于参数优化的自动拾取及质量评估 2.1 拾取参数优化算法为了避免人为参数设置不当而引起的震相误拾，Nippress等(2010)提出了拾取参数优化算法，优化了3种传统算法的相关窗长、静态阈值，判定优化函数(TC)为0的参数为最佳参数，但产生震相漏拾的参数因避开所有判别过程维持初值0而同最佳参数相混淆.本文首先在Nippress等(2010)基础上增加了震相漏拾参数判别环节，然后优化了基于动态阈值触发机制拾取算法(FBpicker、KTpicker、AICDpicker)的参数(表 1).
|
表 1 待优化参数含义及取值范围 Table 1 Definitions and values of picking parameters for three automatic pickers |
改进后的TC由残差精度常数(Ca)、误拾常数1(Cm)、误拾常数2(Cop)、漏拾常数(Ce)相加组成：
|
(5) |
(1) 以accuracy表示自动拾取和人工拾取的时间差值，用于划分P、S波残差精度常数值(CaP、CaS)，accuracy越大则CaP、CaS越大.
|
(6) |
|
(7) |
(2)Cm随误拾个数的增加而增加：
|
(8) |
|
(9) |
(3) 如果在垂直分量上仅拾取P波初至则Cop为0；如果同时拾取P、S波初至则值为1；如果只拾取S波初至则值为2.在水平分量上的定义类似.
|
(10) |
|
(11) |
(4) 如果产生漏拾则Ce为10：
|
(12) |
|
(13) |
改进后的方案给予震相漏拾参数高TC值，从而将其与最佳拾取参数明显区分.由上可知，TC值越小说明自动拾取结果越准确.
2.2 P、S波拾取结果优化在实际波形记录中，不可避免的存在干扰因素.针对该问题，采用时窗限定法去除干扰，只保留时窗内的拾取结果(图 1).如果出现漏拾则降低nsigma(见表 1)或自动替换相应台站TC值小的新参数.对于单个台站，新参数的针对性更强.
|
图 1 地震记录中存在多个地震(a)及干扰(b)实例 绿色虚线是筛选窗口位置，落在窗口内的拾取结果作为该条记录的P波到时.窗口外的红色虚线视为拾取到的干扰. Fig. 1 Examples of earthquakes (a) and interference (b) in seismic records Green dash lines represent the lower and upper ends of the search window. The P-phase arrival time is mostly located in the interval between the two green dash lines. The red dash lines outside the two green dash lines represent the questionable picks. |
P波筛选窗口的上限B1、下限B2为
|
(14) |
|
(15) |
S波筛选窗口的上限B3、下限B4为
|
(16) |
|
(17) |
式中，T是发震时刻，S是S波理论到时.C1、C2源于数据的实测结果，与速度模型的偏差、震中距呈正相关.P是经加权方案评估后质量较高的P波到时，D是理论P、S波到时差.
2.3 质量评估方案根据Diehl等(2009b)的方案，可由accuracy评判自动拾取结果的等级(表 2).
|
|
表 2 自动拾取结果的权重和准确度关系表 Table 2 Relationship between weight and accuracy of automatic picking results |
Di等(2006)制定了拾取到时加权方案，旨在以初至信号的综合变化特征为判据实现拾取结果等级的自动分类.Nippress等(2010)在其基础上，根据拾取点前后1 s内原始波形、特征函数的最大值之比计算了加权因子值(PW).该方法中两个最大值之比很容易产生极值，极大制约了算法的稳健性.本文仅保留PW2(式(18)—(19))，引入了三个新PW.
另一方面，Nippress等(2010)对P、S波加权因子采用相同的计算窗长(1 s)和权重方案.实际上对于水平分量，1 s的窗口太长会因包含了过多P波信息而压制S波初至信息，并且S波PW值通常要低于P波.所以本研究中将S波PW窗长减至0.4 s并依据实际测试结果对P、S波分别建立权重方案.改进后的4个PW分别为：
PW1：拾取点前、后窗长t_ma(s)内特征函数的信噪比.
PW2：
|
(18) |
|
(19) |
式(18)是P波PW2(N取100)，式(19)是S波PW2(N取40).CFM、TM分别是拾取点前滑动窗口N中特征函数、动态阈值的平均值.PW2反应了拾取点前特征函数的扰动与偏差.
PW3：拾取点前、后窗口(P波1 s、S波0.4 s)内原始波形的信噪比.
PW4：拾取点后窗口(P波1 s、S波0.4 s)内特征函数平均值除以拾取点前相同窗长内特征函数平均值.
经反复测试得到P、S波总权重Weight的分类区间：
|
(20) |
|
(21) |
本文收集了腾冲火山区附近2个中国国家数字地震台网固定台站(郑秀芬等，2009)和中国地震局地球物理研究所布设的5个流动台站2015年22325条波形记录用于震相拾取研究.到时拾取过程中，先通过对部分数据的测试得到各台站最佳拾取参数、筛选窗口常数值、权重分隔区间，再运用于全部波形的到时拾取.
3.1 数据集划分依据中国地震台网中心编目系统提供的地震目录，对23°N—26°N，96°E—100°E范围内，ML≥1.5的1533个地震事件的三分量波形记录进行截取，波形采样率为100 Hz.将所有记录划分成2个数据子集，数据集1(2015年5—8月)用于确定相关参数，数据集2(2015年1—4月、9—12月)为参数优化后的应用对象.
3.2 优化结果及分析通过对数据集1的反复测试，得到各台站3种算法P、S波最佳拾取参数和C1、C2常数值(图 3—图 4).计算中使用的P波速度模型如表 3所示(王椿镛等，2000)，采用VP/VS=1.72(王椿镛等，2002；李永华等，2009).由于7个台站波形背景噪声水平相近，所以其P、S波最佳拾取参数比较接近，主要是在阈值控制系数n sigma上有所调整.
|
|
表 3 P波速度模型 Table 3 P-wave velocity model |
|
图 2 研究区台站分布及震中位置分布图 黑色三角表示地震台站；红色圆点表示地震事件. Fig. 2 Map showing the distributions of seismic stations and epicenters of earthquakes in study region The black triangles denote the stations. The red dots are earthquakes. |
|
图 3 P波到时最优拾取参数统计图 图中1~3行分别是FBpicker, KTpicker和AICDpicker方法确定的地震台站(三角形)最优拾取参数.第1~3列为不同的拾取参数, 具体拾取参数见表 1和公式(14). Fig. 3 Statistics of P- wave optimized picking parameters for each picker at each station From upper to lower: FBpicker, KTpicker, and AICDpicker method. All parameters in each column are shown in Table 1 and Formula (14). |
|
图 4 S波到时最优拾取参数统计图 图中1~3行分别是FBpicker, KTpicker和AICDpicker方法确定的地震台站(三角形)最优拾取参数.第1~3列为不同的拾取参数, 具体拾取参数见表 1和公式(15). Fig. 4 Same as Fig. 3 but for S waves |
优化参数后的拾取结果如表 4所示，结果表明经优化参数、拾取结果优化后准确率得到明显提升.P波到时与人工识别结果绝对差值在0.2 s以内的占比(>94.74%)与田优平和赵爱华(2016)方法的效果相近，S波拾取的准确率(>90.59%)高于Castellazzi等(2015)得到的73%的准确率.
|
|
表 4 初至拾取情况统计 Table 4 Picking results of three methods |
加权方案中的4个Weight应当和以人工识别为基准的4个Weight(表 2)一一对应，即图 5和图 6中对角线的网格表示正确的分类，白色背景的网格表示降级，深灰色背景的网格表示提升1~2级，浅灰色背景的网格表示至少提升2级.
|
图 5 P波拾取到时加权方案统计结果 对角线的网格表示正确的分类, 白色背景的网格表示降级, 深灰色背景的网格表示提升1~2级, 浅灰色背景的网格表示至少提升2级.N代表每一列(行)的拾取个数总和. μt和σt分别为对应Weight到时残差的平均值、标准偏差. Fig. 5 Statistics of weighting schemes for automatic P-wave pickers Diagonal elements (bold) represent correctly weighted automatic picks. Off-diagonal elements represent picks that were either downgraded (white cells) or upgraded (grey) by the automatic picker. N represents the number of picks in each column or row. μt denotes the mean difference between the automatic pick and manual pick for all automatic picks classified as each weight t. σt represents the standard deviation of the difference between the automatic pick and manual pick for all automatic picks classified as each weight t. |
图 5和图 6结果表明，3种方法P波、S波拾取到时权重Weight0-1的准确率分别高于96.4%、95.9%，拾取精度高于同类研究方法(Di Stefano et al., 2006；Diehl et al., 2009a)，也体现了本文加权方案的合理性和可靠性.实际数据测试中，计算PW时并没有出现极值从而稳定的反应了初至波拾取点前、后质点振动性质的变化.虽然质量评估方案相对保守造成了很多优质拾取的降级，但也保证了两种判别方法的分类结果相互对应.改进后的加权方案中P、S波拾取结果连续提升2级以上的个数均低于0.62%，表明该加权方案的分类结果是可靠的，实现了拾取到时精度的自动判别.
|
图 6 S波拾取到时加权方案统计结果(具体说明同图 5) Fig. 6 Performance of the automatic S-wave weighting scheme (Notes as in Fig. 5) |
本文改进了Nippress等(2010)提出的拾取参数优化算法和自动拾取到时质量评估方案，引入了拾取结果的优化方案.改进后的优化算法将震相漏拾参数与最优拾取参数加以区分，拾取结果优化方案剔除了各种干扰，改进后的加权方案消除了极值的影响、分别体现了P、S波初至的变化特征.使用改进后的算法得到了云南腾冲地区7个地震台站FBpicker、KTpicker、AICDpicker最佳拾取参数，并对各台站2015年记录的地震事件进行了震相自动拾取、拾取结果优化、拾取结果质量评估.
分析表明，每个台站不同自动拾取算法的最佳参数需分别调整，阈值控制系数的高低取决于台站背景噪声的强弱.结果还表明，经参数优化、拾取优化后，3种方法自动识别的P、S波到时与人工拾取到时的绝对差值在0.1 s以内的分别高于74.66%、70.98%，准确度较高.
本文中的TC函数基于4个常函数的线性叠加，原理简单容易实现；加权方案仅限于时域内的计算较之于频域省时便捷，因而对于不同自动拾取算法的参数确定和拾取到时等级判定，该方法具有一定的普遍性.
致谢 感谢中国地震局地球物理研究所“国家数字测震台网数据备份中心”(doi:10.11998/SeisDmc/SN)提供的固定台地震波形数据.感谢审稿专家提出的宝贵修改意见.
Allen R V. 1978. Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America, 68(5): 1521-1532. |
Allen R V. 1982. Automatic phase pickers:their present use and future prospects. Bulletin of the Seismological Society of America, 72(6B): S225-S242. |
Baer M, Kradolfer U. 1987. An automatic phase picker for local and teleseismic events. Bulletin of the Seismological Society of America, 77(4): 1437-1445. |
Bai C Y, Kennett B L N. 2000. Automatic phase-detection and identification by full use of a single three-component broadband seismogram. Bulletin of the Seismological Society of America, 90(1): 187-198. DOI:10.1785/0119990070 |
Baillard C, Crawford W C, Ballu V, et al. 2014. An automatic kurtosis-based P- and S-phase picker designed for local seismic networks. Bulletin of the Seismological Society of America, 104(1): 394-409. DOI:10.1785/0120120347 |
Bogiatzis P, Ishii M. 2015. Continuous wavelet decomposition algorithm for automatic detection of compressional- and shear-wave arrival times. Bulletin of the Seismological Society of America, 105(3): 1628-1641. DOI:10.1785/0120140267 |
Castellazzi C, Savage M K, Walsh E, et al. 2015. Shear wave automatic picking and splitting measurements at Ruapehu volcano, New Zealand. Journal of Geophysical Research:Solid Earth, 120(5): 3363-3384. DOI:10.1002/2014JB011585 |
Chen C, Holland A A. 2016. PhasePApy:A robust pure Python package for automatic identification of seismic phases. Seismological Research Letters, 87(6): 1384-1396. DOI:10.1785/0220160019 |
Di Stefano R, Aldersons F, Kissling E, et al. 2006. Automatic seismic phase picking and consistent observation error assessment:application to the Italian seismicity. Geophysical Journal International, 165(1): 121-134. DOI:10.1111/j.1365-246X.2005.02799.x |
Diehl T, Deichmann N, Kissling E, et al. 2009a. Automatic S-wave picker for local earthquake tomography. Bulletin of the Seismological Society of America, 99(3): 1906-1920. DOI:10.1785/0120080019 |
Diehl T, Kissling E, Husen S, et al. 2009b. Consistent phase picking for regional tomography models:application to the greater Alpine region. Geophysical Journal International, 176(2): 542-554. DOI:10.1111/j.1365-246X.2008.03985.x |
Gentili S, Bragato P. 2006. A neural-tree-based system for automatic location of earthquakes in Northeastern Italy. Journal of Seismology, 10(1): 73-89. DOI:10.1007/s10950-005-9001-z |
Gentili S, Michelini A. 2006. Automatic picking of P and S phases using a neural tree. Journal of Seismology, 10(1): 39-63. DOI:10.1007/s10950-006-2296-6 |
He X L, She T L, Gao F. 2016. A new method for picking up arrival times of seismic P and S waves automatically. Chinese Journal of Geophysics (in Chinese), 59(7): 2519-2527. DOI:10.6038/cjg20160717 |
Küperkoch L, Meier T, Lee J, et al. 2010. Automated determination of P-phase arrival times at regional and local distances using higher order statistics. Geophysical Journal International, 181(2): 1159-1170. |
Leonard M, Kennett B L N. 1999. Multi-component autoregressive techniques for the analysis of seismograms. Physics of the Earth and Planetary Interiors, 113(1-4): 247-263. DOI:10.1016/S0031-9201(99)00054-0 |
Li Y H, Wu Q J, Tian X B, et al. 2009. Crustal structure in the Yunnan region determined by modeling receiver functions. Chinese Journal of Geophysics (in Chinese), 52(1): 67-80. |
Liu H, Zhang J Z. 2014. STA/LTA algorithm analysis and improvement of microseismic signal automatic detection. Progress in Geophysics (in Chinese), 29(4): 1708-1714. DOI:10.6038/pg20140429 |
Liu J S, Wang Y, Yao Z X. 2013. On micro-seismic first arrival identification:A case study. Chinese Journal of Geophysics (in Chinese), 56(5): 1660-1666. DOI:10.6038/cjg20130523 |
Lomax A, Satriano C, Vassallo M. 2012. Automatic picker developments and optimization:Filter Picker-A robust, broadband picker for real-time seismic monitoring and earthquake early warning. Seismological Research Letters, 83(3): 531-540. DOI:10.1785/gssrl.83.3.531 |
Maeda N. 1985. A method for reading and checking phase time in auto-processing system of seismic wave data. Zisin (Journal of the Seismological Society of Japan. 2nd ser.), 38(3): 365-379. DOI:10.4294/zisin1948.38.3_365 |
Nippress S E J, Rietbrock A, Heath A E. 2010. Optimized automatic pickers:application to the ANCORP data set. Geophysical Journal International, 181(2): 911-925. |
Ross Z E, Ben-Zion Y. 2014. Automatic picking of direct P, S seismic phases and fault zone head waves. Geophysical Journal International, 199(1): 368-381. DOI:10.1093/gji/ggu267 |
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. DOI:10.1109/TGRS.2002.800438 |
Stevenson P R. 1976. Microearthquakes at Flathead Lake, Montana:A study using automatic earthquake processing. Bulletin of the Seismological Society of America, 66(1): 61-80. |
St-Onge A. 2011. Akaike information criterion applied to detecting first arrival times on microseismic data. //81st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts: 1658-1662. |
Tian Y P, Zhao A H. 2016. Automatic identification of P-phase based on wavelet packet and Kurtosis-AIC method. Acta Seismologica Sinica (in Chinese), 38(1): 71-85. DOI:10.11939/jass.2016.01.007 |
Wang C X, Bai C Y, Wang X. 2013. Review of automatic onset time picking for seismic arrivals. Progress in Geophysics (in Chinese), 28(5): 2363-2375. DOI:10.6038/pg20130518 |
Wang C Y, Huangfu G, Wan D B, et al. 2000. Deep seismic sounding of crustal structure in Tengchong volcanic area. Journal of Seismological Research (in Chinese), 23(2): 148-156. |
Wang C Y, Mooney W D, Wang X L, et al. 2002. Study on 3-D velocity structure of crust and upper mantle in Sichuan-Yunnan region, China. Acta Seismologica Sinica (in Chinese), 24(1): 1-16. |
Zhang H L, Zhu G M, Wang Y H. 2013. Automatic microseismic event detection and picking method. Geophysical and Geochemical Exploration (in Chinese), 37(2): 269-273. DOI:10.11720/j.issn.1000-8918.2013.2.17 |
Zheng X F, Ouyang B, Zhang D N, et al. 2009. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake. Chinese Journal of Geophysics (in Chinese), 52(5): 1412-1417. DOI:10.3969/j.issn.0001-5733.2009.05.031 |
何先龙, 佘天莉, 高峰. 2016. 一种地震P波和S波初至时间自动拾取的新方法. 地球物理学报, 59(7): 2519-2527. DOI:10.6038/cjg20160717 |
李永华, 吴庆举, 田小波, 等. 2009. 用接收函数方法研究云南及其邻区地壳上地幔结构. 地球物理学报, 52(1): 67-80. |
刘劲松, 王赟, 姚振兴. 2013. 微地震信号到时自动拾取方法. 地球物理学报, 56(5): 1660-1666. DOI:10.6038/cjg20130523 |
刘晗, 张建中. 2014. 微震信号自动检测的STA/LTA算法及其改进分析. 地球物理学进展, 29(4): 1708-1714. DOI:10.6038/pg20140429 |
田优平, 赵爱华. 2016. 基于小波包和峰度赤池信息量准则的P波震相自动识别方法. 地震学报, 38(1): 71-85. DOI:10.11939/jass.2016.01.007 |
王彩霞, 白超英, 王馨. 2013. 地震震相初至自动检测技术综述. 地球物理学进展, 28(5): 2363-2375. DOI:10.6038/pg20130518 |
王椿镛, 皇甫岗, 万登堡, 等. 2000. 腾冲火山区地壳结构的人工地震探测. 地震研究, 23(2): 148-156. DOI:10.3969/j.issn.1000-0666.2000.02.009 |
王椿镛, Mooney W D, 王溪莉, 等. 2002. 川滇地区地壳上地幔三维速度结构研究. 地震学报, 24(1): 1-16. DOI:10.3321/j.issn:0253-3782.2002.01.001 |
张唤兰, 朱光明, 王云宏. 2013. 基于时窗能量比和AIC的两步法微震初至自动拾取. 物探与化探, 37(2): 269-273. DOI:10.11720/j.issn.1000-8918.2013.2.17 |
郑秀芬, 欧阳飚, 张东宁, 等. 2009. "国家数字测震台网数据备份中心"技术系统建设及其对汶川大地震研究的数据支撑. 地球物理学报, 52(5): 1412-1417. DOI:10.3969/j.issn.0001-5733.2009.05.031 |