强震预警中P波到时STA/LTA和贝叶斯BIC双步骤捡拾研究
杨黎薇+邱志刚+林国良+王玉石
摘要：以云南地区的地震数据为基础，借鉴国内外P波震相自动识别相关研究，提出一套可实时处理P波震相的方法，即STA/LTA和贝叶斯BIC双步骤捡拾法。应用此方法对所选取的云南强震动台网观测记录进行P波自动精确识别，并与人工捡拾方法结果进行对比，确定STA/LTA和贝叶斯BIC双步骤捡拾法的识别精度能满足地震预警快速准确的要求。
关键词：地震预警；P波；强震记录；STA/LTA；贝叶斯准则；双步骤捡拾
中图分类号：P315.91文献标识码：A文章编号：1000-0666（2017）04-0619-10
0引言
地震要素有3个：发震时间、发震地点和震级大小，地震预警的核心内容就是对地震三要素进行快速估计。作为一种有效的地震防灾减灾手段，地震预警必须满足2个特性：可靠性和时效性。首先，预警系统所提供的信息必须具有相对准确性和稳定性；其次，预警系统应能够提供尽可能长的应急反应时间。换言之，地震预警系统就是要满足地震警报发布速度与地震警报内容准确度之间的最大平衡关系。可见，在最短时间内确定震中位置和震级大小是预警的重点环节，而地震定位和震级估算又直接依赖于精确的震相捡拾结果，因此，深入研究地震动震相自动捡拾问题很有必要。
震相自动捡拾是基于识别地震波到达后能量、频率成分、极化等特征上的变化来实现的，多数震相自动捡拾方法是通过提取信号和噪声中的不同特征来区分真实地震信号与噪声（尹得余，2012）。自动捡拾作为地震预警系统高效准确进行震相判别的唯一选择，其准确度对地震定位和震级确定结果有非常重要的影响。近几年，国内外很多学者对此进行了多方位研究，也取得了一些进展（陈少波，2015）。
常用P波震相捡拾方法主要有能量分析、地震波形偏振分析、自回归分析、神经网络、小波变换、高阶统计、复合方法等（史洪山，2012）。到目前为止，还没有任何一种单一方法可以识别所有震相，每一种方法都有一定局限性和应用范围。能量分析是通过选取合适特征函数去突出地震震相主要特征，再使用长短时平均方法捡拾震相。长短时平均（STA/LTA）算法是能量分析中最常用的捡拾P波到时方法，由Stevenson（1976）最早提出并应用于判别地震波初至。随后，Allen（1978，1982）、Baer和Kradolfer（1987）、Earle 和Shearer（1994）对长短时间方法进行不同程度的改进。Ambuter和Solomon（1974）、Anderson（1981）、Mcevilly和Majer（1982）曾利用地震记录幅值的绝对值作为特征函数，Swindell和Snell（1977）利用地震记录幅值的平方作为特征函数，Earle和Shearer（1994）用Hilbert变换求取特征函数，在确定相应特征函数之后，通过长短时平均方法去判定P波震相初动（马强，2008）。这类方法适用于震相清晰的地震事件，当信噪比低或初动不明显时，捡拾效果较差（Stevenson，1976；Allen，1978，1982；Baer，Kradolfer，1987；高淑芳等，2008）。地震波形偏振分析主要是综合三分向地震资料后开展奇异值解析，提出能反映偏振特性的量去进行震相识别（Vernon，1987；Jukevics，1988；刘建华等，2006；马强，2008）。Flinn（1965）、Vidale（1986）、Cichowicz（1993）、Mao和Gubbins（1995）、Earle（1999）发现震相到来时质点运动的偏振方向会发生改变，于是他们提取出使用地震波的偏振特征去捡拾和判别震相到时。这类方法主要应用于检测原地应力和地震预报等方面，其对资料窗口长度和信噪比变化特别敏感，会产生同向轴重叠问题。自回归分析先假设震相到时前后地震状态为2个不同的稳态过程，再找到震相到时点，即2个稳态过程划分点（王继等，2006；周彦文，2008；Sleeman，1999）。Maeda（1985）提出一种时域赤池信息量准则，即直接通过地震波形记录数据求取AIC函数，跳过对AR模型阶数的求取过程。Sleeman（1999）提出了自回归方法，在此基础上，Leonard and Kennett（1999）和Leonard（2000）将自回归的假定AR-AIC用作P波震相的自动捡拾。这类方法属于目前应用最为广泛的P波震相识别方法，其精确度和稳定性能满足地震预警的工作需要，算法高效但模型结构复杂，计算量也较大。神经网络把与地震资料相关的可用信息作为网络输入去识别震相（Wang，Teng，1995；Zhao，Takano 1999）。Zhao和Takano（1999）、Dai和MacBeth（1995，1997）、Murat和Rudman（1992）、Wang 和Teng（1995）都對神经网络方法进行震相识别做过研究。Tong和Kemett（1996）利用该方法进行震相识别，建立有关地震波特性的专家知识库，再根据检测到的震相特性去推理判别震相的类型。这种理论化的人脑神经网络数学模型，是基于模仿大脑神经网络结构和功能而建立的一种信息处理系统，具有高度非线性，能够进行复杂的逻辑操作和非线性实现系统，捡拾精确度高。小波分析把地震记录P波到时的特征表现在不同小波变换系数上，Chakraborty和Okaya（1995）、Kanawaldip和Farid（1997）提出小波变换的主要成分分析，刘希强等（1998，2000a，b）对该算法做出进一步探索；Zhang等（2003）将AIC准则与小波变换结合起来进行震相捡拾，小波变换虽然在震相识别方面取得了一定成效，但其本质依然依托于Fourier变换，摆脱不了Fourier分析的局限性。这类方法会限制信号的精细分辨，随着尺度增大，小波变换相应正交基函数的频谱局部性就越差。高阶统计是引入偏度和峰度2个量去自动识别地震资料从高斯段变为非高斯段的点，该点即被认为对应地震信号的初至（Saragiotis et al，2002）。Saragiotis等（2002）为了消除人为因素和噪声对震相捡拾的影响，提出运用高阶统计学的方法来捡拾震相。这类方法对波形清晰程度的要求高，捡拾结果稳定性不算理想。复合方法则是把2种或2种以上的方法联合使用，进而获取更好的结果（Bai，Kennett，2000）。Bai和Kennett（2000）提出将能量分析方法、瞬时频率分析方法、自回归分析方法等结合起来捡拾震相。这类方法在一定程度上弥补了单一算法的缺陷，权衡考虑各方面因素，得到较完善精确的结果。综上所述，本文采用复合方法自动识别P波到时。
本文借鉴国内外P波震相自动识别的相关研究，以云南地区的地震数据为基础，提出一套可实时处理P波震相的复合方法，即STA/LTA和贝叶斯BIC双步骤捡拾法。应用此方法对所选取的云南强震动台网观测记录进行P波自动精确识别，并通过比较此方法与人工捡拾方法的实测结果，以确定STA/LTA和贝叶斯BIC双步骤捡拾法的识别精度是否满足地震预警对时效性与准确度的要求。
1基于STA/LTA和BIC算法的P波震相自动捡拾研究
1.1P波震相特征函数震相在地震图上显示为性质不同或传播路径不同的地震波组。各种震相在到时、波形、振幅、周期和质点运动方式等方面都各有特征。在P波震相中，质点沿着波的传播方向运动，通常在震中距105°的范围内，P波震相是地震图上的初至震相。
对地震记录来说，P波在垂直向的特征明显，一般在地震波P波初动点的幅值和频率会明显发生变化，在大多数情况下，可以观察到一个明显的向上或向下的脉动，由此可判定P波初动点。所以用于识别P波震相的特征参量必须能够反映这些变化趋势，才能准确地识别P波震相（毛燕等，2011）。选取特征函数是为了使用STA/LTA震相识别方法时，引入一个新的时间序列CF（i）来反映原始记录信号的变化特征，以此为基础再计算STA/LTA值。特征函数的选取对识别结果精度至关重要，它必须能灵敏地反映信号到达时的频率或幅值特征变化，甚至能增强这些变化。为了将P波的特征尽量放大，本文使用Allen（1978）所提出的特征函数进行波相特征放大处理：CFP=xud（k）2+[xud（k）-xud（k-1）]2（1）式中：CFP表示P波捡拾的特征函数；xud（k）为第k时刻竖向的速度记录。使用（1）式主要是考虑到P波作为纵波，在传递的过程中应该保持其传播方向与振动方向一致，故而采用垂直记录上的变化去体现P波的特征。
1.2STA/LTA+BIC双步骤捡拾法
STA/LTA+BIC双步骤捡拾法主要是将常用的长短时平均算法（STA/LTA）与目前国际上较受关注的贝叶斯分割算法（BIC）结合起来对地震P波震相进行精确捡拾，确定出一种性能稳定、准确度高、可应用于地震预警系统中P波震相自动判别的方法。
1.2.1STA/LTA捡拾
长短时平均（STA/LTA）算法是地震预警中最常用的捡拾P波到时方法，这种方法反应了幅值的瞬时变化。STA表征信号短时平均值，主要反映信号幅值瞬时的变化；LTA表征信号长时平均值，主要反映相对于待检信号的背景噪声平均水平。STA比LTA变化快很多，该算法通过STA与LTA之比反映信号水平或者能量的变化。当地震信号到达时，STA/LTA值会出现明显增大，当该比值大于预先所设定的阈值，其对应的时刻点则被认为是P波的初动点。
STA/LTA的计算值体现了长短时窗内的能量比，本文应用时间窗的捡拾公式为（Allen，1982）：STA（i）/LTA（i）=∑ik1CFP（i）/（i-k1+1）∑ik2CFP（i）/（i-k2+1）（2）式中：STA（i）和LTA（i）分别表示P波信号在i时刻的短时和长时平均值；CFP（i）为P波信号在i时刻的特征函数值；i为变化时刻点；k1和k2均为变化时刻前某一时刻点（k2 <k1<i）。
P波到达前后，记录的幅值变化较大，STA平均值的窗长取值与待测信号的周期相关，决定了记录事件是否能够正确触发，取值太短会误触发，取值太长会漏触发；而LTA平均值的窗长取值仅是与信号的背景噪声水平相关，取值范围相对固定。以2014年10月7日发生的云南景谷M6.6主震中景谷永平台的UD向记录为例，STA和LTA时间窗长的选取对计算结果的精度有较大的影响，其值的选取需结合特征函数的特点与STA/LTA计算方法来确定，同时为了准确判定地震事件初动点，特征函数的阈值通常都会定为10。在本文中，设定特征函数阈值参数为10，长窗长取15 s，短窗长取0.5 s，即可达到较为理想的捡拾效果。当特征函数值大于10时，可判定有地震事件发生，如图1所示，该记录的P波初动点出现在曲线跳跃点处，即13.49 s处；用STA/LTA方法捡拾的到P波也在13.49 s处对应触发。这种以能量变化实现特征捡拾的方法对于地震弱信号适应性较强，捡拾效果高，适用于震相清晰、规则型的地震事件。
1.2.2貝叶斯（BIC）信息准则
当前有许多模型选择方法，如似然法、贝叶斯法和信息准则法等。1973年统计学家Hirotugu Akaike提出Akaike 信息准则（AIC）用来进行模型选择，对模型信息的发展产生了深远影响，Schwarz（1978）在此基础上提出了贝叶斯（BIC）信息准则。
最大似然准则是一种具有理论性的点估计法，此方法的基本思想是：从模型总体随机抽取n组样本观测值后，最合理的参数估计量以观测值出现的概率最大为最大估计准则。而模型选择问题即是从K个候选模型M={Mi；i=1，2，…，K}中选出一个最能表示给定N个数据X= {Xi；i=1，2，…，N}分布的模型选择问题。BIC准则是一种基于模型复杂度惩罚机制的渐进最大似然准则，用来解决模型选择问题（薛昊，2010）。
假设数据集X对于每个模型Mj的最大似然式为L（x，Mj），kj为模型Mj的参数个数，那么模型Mj的BIC计算式可用（3）式表示：BIC（Mj）=lgL（x，Mj）-12λkjlg（N）（3）式中：对数似然函数lgL（x，Mj）反映了模型精确度，它可以测量效率参数化模型的预测数据，也可用来选择集群数量内出现的一个特定数据集；λ为惩罚因子；λkjlg（N）反映了对模型复杂度的惩罚，它可提高模型的灵活性，满足规则内的惩罚需求（薛昊，2010）。
由于样本数据的似然值会随着模型参数的增加而增加，单纯采用极大似然估计会造成模型参数过多、维数过大，所以我们在BIC准则中引入了对模型复杂度的惩罚。由（3）式使BIC计算取得最大值得模型，即为基于贝叶斯信息准则的最优模型： </k1
1.2.3BIC算法流程
在波形图中，并不是只会出现单一的分割点，很多时候，会有多个分割点出现，为了得到最优分割点，本文采用贝叶斯准则中顺序检测算法逐个找出分割点，最后得到分割点的最大值，确定P波震相的精确到时。具体流程如图2所示。
图2中，首先是初始化检测窗[a，b]，且a=1；b=1。这一步骤表明视窗一开始仅包含2个样本点a和b。接着开始进行BIC检测[a，b]窗内是否存在分割点，不存在分割点时，对检测窗窗长增加1，即令b=b+1，重新开始BIC检测分割点；存在分割点时，则移动检测窗，即a=i+1；b=a+1，i为分割点，读取该分割点值并按顺序逐步增加检测窗长度，直至整个完整波频段检测结束。
从整个算法的流程可看出，算法的时间复杂度和分割点数有关，虽然BIC的检测精度很高，但其在检测过程中计算量非常大。
从长短时平均（STA/LTA）方法的捡拾效果来看，其算法简单、稳定快速，主要缺点是捡拾到的P波到时往往比实际到时滞后。而贝叶斯（BIC）方法捡拾精度高、稳定可靠，但其捡拾P波到时需要对整个完整波频按顺序逐步检测分割点，最后才能判定出所有分割点的最大值，确定出P波到时。因此，如果采用STA/LTA粗略捡拾到P波到时后，缩小检测范围，再用BIC准则进行P波震相的精确拾取，比在整条地震记录图上直接应用BIC准则要更快，精度与直接应用BIC准则捡拾相同。
ΔBIC（i）>0时，表明BIC（M2）>BIC（M1），即模型M2优于模型M1，当ΔBIC（i）的取值为最大值，则认为这是最优分割点i，最终使ΔBIC（i）取得最大值得帧i为P波震相到时点。i=argmax1<i0ΔBIC（i）（11）使用BIC准则方法捡拾2014年云南景谷M6.6主震中景谷永平台站UD向记录的P波到时，获得结果如图3所示。图3a描述了在该加速度时程曲线中，P波到时相对于记录起始时刻是13.37 s；
人工捡拾景谷永平台UD向记录P波到时为13.40 s，STA/LTA自动捡拾结果为13.49 s，比人工捡拾结果滞后0.09 s；BIC自动捡拾结果为13.37 s，比人工捡拾结果提前0.03 s。
1.3STA/LTA+BIC双步骤捡拾的工作流程
STA/LTA+BIC双步骤捡拾法的工作流程如图4所示，强震仪的阈值产生触发事件时，实时数据的特征函数可直接采用（1）式计算特征函数值；STA/LTA粗略捡拾P波到时采用（2）式计算；BIC精确捡拾P波到时采用（11）式计算。首先根据特征函数计算结果，对记录波形展开P波粗略捡拾，应用STA/LTA长短时平均法检测或报警触发事件，当STA/LTA刻画出记录幅值的瞬时变化时，立刻判别P波阈值是否触发，若是STA/LTA计算结果超过常数阈值THR，则认为是地震事件，进入下一步精确捡拾；反之则是干扰事件，捡拾工作及时结束。本方法对地震事件P波的精确捡拾采用BIC自动捡拾法，其捡拾结果作为最终的P波到时成为地震预警的重要参数。
归纳STA/LTA+BIC双步骤捡拾法的基本步骤主要有2点：（1）先用STA/LTA法粗略捡拾P波到时；（2）以粗略到时点为基准，前后各推0.5 s，在此范围内使用BIC准则进行精确捡拾。
2实际算例分析
根据本文对震相捡拾方法的原理介绍与工作流程，对云南省数字强震动台网记录的17次ML≥5.0地震事件进行STA/LTA+BIC双步骤自动捡拾，捡拾记录均来自震中距在100 km范围内的强震动台站记录。
2.1数据选取
本文所用数据来自于云南省数字强震动台网获取的强震记录，包括了2008—2014年17次M5.0～7.0地震事件，涵盖了云南盈江、腾冲、彝良、洱源、景谷等地震多发区域。对P波震相而言，其在竖直向的幅值变化较大，我们对地震动记录进行初次人工筛选，剔除丢头、叠加等波形异常的强震记录，选取震中距在100 km以内、P波在竖直向相对突出的竖向原始加速度时程记录173个，并对数据进行预处理。筛选后的各地区的加速度时程记录数量分布如表1所示。
2.2实际捡拾情况分析
如表1所示，台网内选取这17个地震事件的实际记录，有9个震中距在100 km内的记录出现漏检，其它记录捡拾工作正常。分析原因发现，其中有4个记录人工捡拾也无法分辨P波时，2个记录中的地震波形有其他小震的干扰叠加，其余3个原因不明。
表2中，STA/LTA捡拾与人工捡拾的P波到时平均差值为0.698 s，标准差值为0.531 s；STA/LTA+BIC捡拾与人工捡拾的P波到时平均差值为0.091 s，标准差值为0.156 s。从地震系统的工作要求而言，STA/LTA捡拾结果误差相对较大，其精确度可达到地震烈度速报的工作要求；STA/LTA+BIC捡拾结果误差相对很小，其精确度能够达到地震预警系统的工作要求。
图5显示了STA/LTA自动捡拾与人工捡拾P波到时结果的偏差分布，可以看出，STA/LTA自动捡拾P波到時的偏差基本以正数为主，最大偏差将近2.0，其捡拾结果通常大于人工捡拾结果，即其对P波到时的自动捡拾往往存在滞后的现象。图6表示了STA/LTA+BIC自动捡拾与人工捡拾P波到时结果的偏差分布，可以看出，STA/LTA+BIC自动捡拾P波到时的偏差均匀分布在-0.5～0.5之间，以0～0.2之间的分布最广。该方法对P波到时的自动捡拾结果与人工捡拾结果吻合。综上所述，采用STA/LTA+BIC自动捡拾P波到时的结果与人工捡拾P波到时结果偏差很小，误差不超过0.5。
2.3震中距对P波到时准确度及时效性的影响
为了研究震中距是否会影响地震P波震相到时的准确程度，我们统计了人工判别、STA/LTA和STA/LTA+BIC自动判别3种捡拾结果，将记录台站的震中距范围扩大至0～160 km，具体变化如图7所示。由图7可见，采用STA/LTA+BIC自动图5STA/LTA自动捡拾与人工捡拾判别方法与人工捡拾方法的捡拾结果几乎是完全吻合的，而STA/LTA自动判别方法捡拾的结果有延时滞后现象产生。同样，将STA/LTA+BIC自动捡拾P波到时的偏差与震中距的变化分布绘制成图8，其结果与图6相符，无论震中距如何变化，采用STA/LTA+BIC自动捡拾P波到时的偏差都很小。由此可见，震中距远近对P波到时的准确程度没有任何影响，采用合适的捡拾方法才能提高捡拾结果的精确度。为了得到不同震中距下采用STA/LTA+BIC自动捡拾P波到时的速度，本文引入速度节奏来显示捡拾速度，即规定一个节奏值，每个数据捡拾的用时以规定节奏值相应的倍数表示，这样就能通过不同的节奏数值看出每个数据的捡拾用时长短，进而推断出P波传播于不同介质时的速度变化。如图9所示，采用STA/LTA+BIC自动捡拾P波到时速度节奏随震中距的变化是逐渐增大的。究其原因，可能是随着震中距逐渐增大，不同场地介质会对P波传播速度及方向产生改变，影响后续对P波震相识别的速度。 </i
綜上所述，震中距的远近对P波震相识别的精度没有影响，P波到时捡拾结果的准确程度取决于采用方法是否适合；但震中距的远近对P波震相识别用时长短有一定影响，随着震中距逐渐增大，自动捡拾P波初至的耗时也会呈逐渐增长趋势。
3讨论与结论
本文借鉴了国内外P波震相自动识别的相关研究，以云南地区的地震数据为基础，提出一套可实时处理P波震相的方法，即STA/LTA和贝叶斯BIC双步骤捡拾法。应用此方法对所选取的云南强震动台网观测记录进行P波自动精确识别，并与人工捡拾方法的实测结果进行比较，得到以下结论：
（1）STA/LTA自动捡拾与人工捡拾P波到时的偏差基本以正数为主，最大偏差将近2.0 s，其捡拾结果通常大于人工捡拾结果，即其对P波到时的自动捡拾往往存在滞后的现象。STA/LTA+BIC自动捡拾P波到时的偏差在0.0值上下均匀分布，主要分布在-0.2～0.2 s之间，最大偏差也超不过0.5 s，该方法其对P波到时的自动捡拾结果与人工捡拾结果接近吻合。
（2）长短时平均（STA/LTA）方法捡拾速度快；贝叶斯（BIC）方法精确度高。因此，在用STA/LTA粗略捡拾到P波后，再用BIC准则进行P波震相的精确拾取，从地震预警的实时性与时效性要求来看，效果比在整条地震记录图上直接应用BIC准则要好。
（3）通过大量数据实际验证，P波特征函数触发阈值设为10较合适。阈值小于10误判率会增高，阈值大于10漏检率也会增高。进行P波到时捡拾时，信噪比高捡拾结果理想，若是其中有噪声等干扰信号影响，捡拾质量会受到很大影响。
（4）采用P波震相自动捡拾方法，结果虽然具有较高的精度，但对S波震相自动捡拾并未尝试，可在下一步工作中继续研究S波震相自动捡拾。
通过对比分析，STA/LTA与BIC相结合的P波到时捡拾方法不但利用了STA/LTA稳定可靠的特点，在初步捡拾点的附近选取一个恰当的时间窗包含P波震相初至过程，再以这个合适的时间窗进一步开展BIC准则精确捡拾，捡拾结果稳定可靠，准确性高。STA/LTA和贝叶斯BIC双步骤捡拾法的识别精度能够达到地震预警的要求，可作为云南地区预警系统建设的技术参考。
参考文献：
陈少波.2015.地震预警系统的P波震相自动识别方法研究[D].哈尔滨：中国地震局工程力学研究所.
高淑芳，李山有，武东坡，等.2008.一种改进的STA/LTA震相自动识别方法[J].世界地震工程，24（2）：37-41.
刘建华，刘福田，胥颐.2006.三分量地震资料的偏振分析[J].地球物理学进展，21（1）：6-10.
刘希强，周蕙兰，李红.2000a.基于小波包变换的地震数据时频分析方法[J].西北地震学报，22（2）：143-146.
刘希强，周蕙兰，沈萍，等.2000b.用于三分向记录震相识别的小波变换方法[J].地震学报，32（2）：125-131.
刘希强，周蕙兰，郑治真，等.1998.基于小波包变换的弱震相识别方法[J].地震学报，20（4）：373-380.
马强.2008.地震预警技术研究及应用[D].哈尔滨：中国地震局工程力学研究所.
毛燕，崔建文，郑定昌，等.2011.地震记录的P波自动捡拾[J].地震研究，34（1）：47-51.
史洪山.2012.地震预警震级测定技术研究[D].哈尔滨：中国地震局工程力学研究所.
王继，陈九辉，刘启元，等.2006.流动地震台站观测初至震相的自动检测[J].地震学报，28（1）：42-51.
薛昊.2010.基于BIC的通用音频分割方法研究[D].哈尔滨：哈尔滨工业大学.
尹得余.2012.P波震相自动捡拾与地震预警震级实时测定技术研究[D].哈尔滨：中国地震局工程力学研究所.
周彦文.2008.基于单台P波记录的早期地震预警方法研究[D].兰州：中国地震局兰州地震研究所.
ALLEN R E.1978.Automatic earthquake recognition and timing from single traces[J].Bulletin of the Seismological Society of America，68（5）：1521-1532.
ALLEN R E.1982.Automatic phase pickers：Their present use and future prospects[J].Bulletin of the Seismological Society of America，72（6B）：225-242.
AMBUTER B P，SOLOMON.1974.An event-recording system for monitoring small earthquakes[J].Bulletin of Seismological Society of America，64（4）：1181-1188.
ANDERSON K R.1981.Epicentral location using arrival time order[J].Bulletin of the Seismological Society of America，70（2）：541-545.
BAER M，KRADOLFER U.1987.An Automatic phase picker for local and teleseismic[J].Bulletin of the Seismological Society of America，72（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[J].Bulletin of Seismological Society of America，90（1）：187-198.
CHAKRABORT Y A，OKAYA D.1995.Frequency-time decomposition of seismic data using wavelet-based methods[J].Geophysics，60（6）：1906-1916.
TONG C，KENNETT B L N.1996.Auto seismic event recognition and later phase identification for broadband seismograms[J].Bulletin of the Seismological Society of America，86（6）：1896-1909.
CICHOWICZ A R.1993.An automatic S-phase picker[J].Bulletin of the Seismological Society of America，83（1）：180-189.
DAI H C，MACBETH C.1995.Automatic picking of seismic arrivals in local earthquake data using an artificial neural network[J].Geophysical Journal International，120（3）：758-774.
DAI H C，MACBETH C.1997.The application of back-propagation neural network to automatic picking seismic arrivals from single-component recordings[J].J Geophys Res，102（B7）：15105-15113.
EARLA P.1999.Polarization of the Earths teleseismic wave field，Geophys[J].J Int，139（1）：1-8.
EARLE P S，SHEARER P M.1994.Characterization of glob seismograms using an automatic-picking algorithm[J].Bulletin of the Seismological Society of America，84（2）：366-376.
FLINN E.1965.Signal analysis using recti linearity and direction of particle motion[J].Proc IEEE 53：1874-1876.
JUKEVICS A.1988.Polarization analysis of three-component array data[J].Bulletin of the Seismological Society of America，78（5）：1725-1743.
KANWALDIP A S，FARID U D.1997.Wavelet transform methods for phase identification in three-component seismograms[J].Bull Seism Soc Amer，87（6）：1598-1612.
LEONARD M，KENNETT B L N.1999.Multi-component autoregressive techniques for the analysis of seismograms[J].Phys Earth Planet Interiors，113（1-4）：247-264.
LEONARD M.2000.Comparison of manual and automatic onset time picking[J].Bull Seism Soc Am，90：1384-1390.
MAEDA N.1985.A method for reading and checking phase times in auto-processing system of seismic wave data[J].The Seismological Society of Japan，38：365-379.
MAO W，GUBBINS D.1995.Simultaneous determination of time delays and stacking weights in seismic array beam forming[J].Geophysics，60（2）：491-502.
MCEVILLY T V，MAJER E L.1982.An automated seismic processor for micro earthquake network[J].Bulletin of the Seismological Society of America，72（1）：303-325.
MURAT M E，RUDMAN A J.1992.Automated first arrival picking：a neural network approach[J].Geophysics Prospecting，40（6）：587-604.
SARAGIOTIS C H D，HADJILEONTIADIS L J，PANAS S M P.2002.AI-S/K：a robust automatic seismic P-phase arrival identification scheme[J].IEEE Sransactions on Geoscience and Remote Sensing，40（6）：1395-1404.
SCHWARZ E G.1978.Estimating the dimension of a model[J].Annals of Statistics，6（2）：461-464.
SLEEMAN R.1999.Robust automatic P-phase picking：an on-line implementation in the analysis of broadband Seismogram recordings[J].Phys Earth Planet Interiors，113：265-275.
STEVENSON P S.1976.Micro earthquakes at flathead lake，Montana：a study using automatic earthquake processing[J].Bulletin of the Seismological Society of America，66（1）：1-80.
SWINDELL H，SNELL N.1977.Station Processor Automatic Signal Detection System，Phase Ι：Final Report，Station Processor Software Development[J].Texas Instruments Report，ALEX（01）-FR-77-01.
VERNON F L.1987.Frequency dependent polarizatianalysis of high-frequency seismograms[J].Journal of geophysical research，92（B12）：664-674.
VIDALE T.1986.Complex polarization analysis of particle motion[J].Bull Seism Soc Am，76（5）：1393-1405.
WANG J，TENG T L.1995.Artificial neural network-based seismic detector[J].Bulletin of the Seismological Society of America，85（1）：308-319.
ZHAO Y，TAKANO K.1999.An artificial neural network approach for broadband seismic phase Picking[J].Bulletin of the Seismological Society of America，89（3）：70-680.
ZHANG H，THURBER C，ROWE C.2003.Automatic P-Wave Arriva Detection and Picking with Multiscale Wavelet Analysis for Single-Component Recordings[J].Bulletin of the Seismological Society of America，93（5）：1904-1912.
摘要：以云南地区的地震数据为基础，借鉴国内外P波震相自动识别相关研究，提出一套可实时处理P波震相的方法，即STA/LTA和贝叶斯BIC双步骤捡拾法。应用此方法对所选取的云南强震动台网观测记录进行P波自动精确识别，并与人工捡拾方法结果进行对比，确定STA/LTA和贝叶斯BIC双步骤捡拾法的识别精度能满足地震预警快速准确的要求。
关键词：地震预警；P波；强震记录；STA/LTA；贝叶斯准则；双步骤捡拾
中图分类号：P315.91文献标识码：A文章编号：1000-0666（2017）04-0619-10
0引言
地震要素有3个：发震时间、发震地点和震级大小，地震预警的核心内容就是对地震三要素进行快速估计。作为一种有效的地震防灾减灾手段，地震预警必须满足2个特性：可靠性和时效性。首先，预警系统所提供的信息必须具有相对准确性和稳定性；其次，预警系统应能够提供尽可能长的应急反应时间。换言之，地震预警系统就是要满足地震警报发布速度与地震警报内容准确度之间的最大平衡关系。可见，在最短时间内确定震中位置和震级大小是预警的重点环节，而地震定位和震级估算又直接依赖于精确的震相捡拾结果，因此，深入研究地震动震相自动捡拾问题很有必要。
震相自动捡拾是基于识别地震波到达后能量、频率成分、极化等特征上的变化来实现的，多数震相自动捡拾方法是通过提取信号和噪声中的不同特征来区分真实地震信号与噪声（尹得余，2012）。自动捡拾作为地震预警系统高效准确进行震相判别的唯一选择，其准确度对地震定位和震级确定结果有非常重要的影响。近几年，国内外很多学者对此进行了多方位研究，也取得了一些进展（陈少波，2015）。
常用P波震相捡拾方法主要有能量分析、地震波形偏振分析、自回归分析、神经网络、小波变换、高阶统计、复合方法等（史洪山，2012）。到目前为止，还没有任何一种单一方法可以识别所有震相，每一种方法都有一定局限性和应用范围。能量分析是通过选取合适特征函数去突出地震震相主要特征，再使用长短时平均方法捡拾震相。长短时平均（STA/LTA）算法是能量分析中最常用的捡拾P波到时方法，由Stevenson（1976）最早提出并应用于判别地震波初至。随后，Allen（1978，1982）、Baer和Kradolfer（1987）、Earle 和Shearer（1994）对长短时间方法进行不同程度的改进。Ambuter和Solomon（1974）、Anderson（1981）、Mcevilly和Majer（1982）曾利用地震记录幅值的绝对值作为特征函数，Swindell和Snell（1977）利用地震记录幅值的平方作为特征函数，Earle和Shearer（1994）用Hilbert变换求取特征函数，在确定相应特征函数之后，通过长短时平均方法去判定P波震相初动（马强，2008）。这类方法适用于震相清晰的地震事件，当信噪比低或初动不明显时，捡拾效果较差（Stevenson，1976；Allen，1978，1982；Baer，Kradolfer，1987；高淑芳等，2008）。地震波形偏振分析主要是综合三分向地震资料后开展奇异值解析，提出能反映偏振特性的量去进行震相识别（Vernon，1987；Jukevics，1988；刘建华等，2006；马强，2008）。Flinn（1965）、Vidale（1986）、Cichowicz（1993）、Mao和Gubbins（1995）、Earle（1999）发现震相到来时质点运动的偏振方向会发生改变，于是他们提取出使用地震波的偏振特征去捡拾和判别震相到时。这类方法主要应用于检测原地应力和地震预报等方面，其对资料窗口长度和信噪比变化特别敏感，会产生同向轴重叠问题。自回归分析先假设震相到时前后地震状态为2个不同的稳态过程，再找到震相到时点，即2个稳态过程划分点（王继等，2006；周彦文，2008；Sleeman，1999）。Maeda（1985）提出一种时域赤池信息量准则，即直接通过地震波形记录数据求取AIC函数，跳过对AR模型阶数的求取过程。Sleeman（1999）提出了自回归方法，在此基础上，Leonard and Kennett（1999）和Leonard（2000）将自回归的假定AR-AIC用作P波震相的自动捡拾。这类方法属于目前应用最为广泛的P波震相识别方法，其精确度和稳定性能满足地震预警的工作需要，算法高效但模型结构复杂，计算量也较大。神经网络把与地震资料相关的可用信息作为网络输入去识别震相（Wang，Teng，1995；Zhao，Takano 1999）。Zhao和Takano（1999）、Dai和MacBeth（1995，1997）、Murat和Rudman（1992）、Wang 和Teng（1995）都對神经网络方法进行震相识别做过研究。Tong和Kemett（1996）利用该方法进行震相识别，建立有关地震波特性的专家知识库，再根据检测到的震相特性去推理判别震相的类型。这种理论化的人脑神经网络数学模型，是基于模仿大脑神经网络结构和功能而建立的一种信息处理系统，具有高度非线性，能够进行复杂的逻辑操作和非线性实现系统，捡拾精确度高。小波分析把地震记录P波到时的特征表现在不同小波变换系数上，Chakraborty和Okaya（1995）、Kanawaldip和Farid（1997）提出小波变换的主要成分分析，刘希强等（1998，2000a，b）对该算法做出进一步探索；Zhang等（2003）将AIC准则与小波变换结合起来进行震相捡拾，小波变换虽然在震相识别方面取得了一定成效，但其本质依然依托于Fourier变换，摆脱不了Fourier分析的局限性。这类方法会限制信号的精细分辨，随着尺度增大，小波变换相应正交基函数的频谱局部性就越差。高阶统计是引入偏度和峰度2个量去自动识别地震资料从高斯段变为非高斯段的点，该点即被认为对应地震信号的初至（Saragiotis et al，2002）。Saragiotis等（2002）为了消除人为因素和噪声对震相捡拾的影响，提出运用高阶统计学的方法来捡拾震相。这类方法对波形清晰程度的要求高，捡拾结果稳定性不算理想。复合方法则是把2种或2种以上的方法联合使用，进而获取更好的结果（Bai，Kennett，2000）。Bai和Kennett（2000）提出将能量分析方法、瞬时频率分析方法、自回归分析方法等结合起来捡拾震相。这类方法在一定程度上弥补了单一算法的缺陷，权衡考虑各方面因素，得到较完善精确的结果。综上所述，本文采用复合方法自动识别P波到时。
本文借鉴国内外P波震相自动识别的相关研究，以云南地区的地震数据为基础，提出一套可实时处理P波震相的复合方法，即STA/LTA和贝叶斯BIC双步骤捡拾法。应用此方法对所选取的云南强震动台网观测记录进行P波自动精确识别，并通过比较此方法与人工捡拾方法的实测结果，以确定STA/LTA和贝叶斯BIC双步骤捡拾法的识别精度是否满足地震预警对时效性与准确度的要求。
1基于STA/LTA和BIC算法的P波震相自动捡拾研究
1.1P波震相特征函数震相在地震图上显示为性质不同或传播路径不同的地震波组。各种震相在到时、波形、振幅、周期和质点运动方式等方面都各有特征。在P波震相中，质点沿着波的传播方向运动，通常在震中距105°的范围内，P波震相是地震图上的初至震相。
对地震记录来说，P波在垂直向的特征明显，一般在地震波P波初动点的幅值和频率会明显发生变化，在大多数情况下，可以观察到一个明显的向上或向下的脉动，由此可判定P波初动点。所以用于识别P波震相的特征参量必须能够反映这些变化趋势，才能准确地识别P波震相（毛燕等，2011）。选取特征函数是为了使用STA/LTA震相识别方法时，引入一个新的时间序列CF（i）来反映原始记录信号的变化特征，以此为基础再计算STA/LTA值。特征函数的选取对识别结果精度至关重要，它必须能灵敏地反映信号到达时的频率或幅值特征变化，甚至能增强这些变化。为了将P波的特征尽量放大，本文使用Allen（1978）所提出的特征函数进行波相特征放大处理：CFP=xud（k）2+[xud（k）-xud（k-1）]2（1）式中：CFP表示P波捡拾的特征函数；xud（k）为第k时刻竖向的速度记录。使用（1）式主要是考虑到P波作为纵波，在传递的过程中应该保持其传播方向与振动方向一致，故而采用垂直记录上的变化去体现P波的特征。
1.2STA/LTA+BIC双步骤捡拾法
STA/LTA+BIC双步骤捡拾法主要是将常用的长短时平均算法（STA/LTA）与目前国际上较受关注的贝叶斯分割算法（BIC）结合起来对地震P波震相进行精确捡拾，确定出一种性能稳定、准确度高、可应用于地震预警系统中P波震相自动判别的方法。
1.2.1STA/LTA捡拾
长短时平均（STA/LTA）算法是地震预警中最常用的捡拾P波到时方法，这种方法反应了幅值的瞬时变化。STA表征信号短时平均值，主要反映信号幅值瞬时的变化；LTA表征信号长时平均值，主要反映相对于待检信号的背景噪声平均水平。STA比LTA变化快很多，该算法通过STA与LTA之比反映信号水平或者能量的变化。当地震信号到达时，STA/LTA值会出现明显增大，当该比值大于预先所设定的阈值，其对应的时刻点则被认为是P波的初动点。
STA/LTA的计算值体现了长短时窗内的能量比，本文应用时间窗的捡拾公式为（Allen，1982）：STA（i）/LTA（i）=∑ik1CFP（i）/（i-k1+1）∑ik2CFP（i）/（i-k2+1）（2）式中：STA（i）和LTA（i）分别表示P波信号在i时刻的短时和长时平均值；CFP（i）为P波信号在i时刻的特征函数值；i为变化时刻点；k1和k2均为变化时刻前某一时刻点（k2 <k1<i）。
P波到达前后，记录的幅值变化较大，STA平均值的窗长取值与待测信号的周期相关，决定了记录事件是否能够正确触发，取值太短会误触发，取值太长会漏触发；而LTA平均值的窗长取值仅是与信号的背景噪声水平相关，取值范围相对固定。以2014年10月7日发生的云南景谷M6.6主震中景谷永平台的UD向记录为例，STA和LTA时间窗长的选取对计算结果的精度有较大的影响，其值的选取需结合特征函数的特点与STA/LTA计算方法来确定，同时为了准确判定地震事件初动点，特征函数的阈值通常都会定为10。在本文中，设定特征函数阈值参数为10，长窗长取15 s，短窗长取0.5 s，即可达到较为理想的捡拾效果。当特征函数值大于10时，可判定有地震事件发生，如图1所示，该记录的P波初动点出现在曲线跳跃点处，即13.49 s处；用STA/LTA方法捡拾的到P波也在13.49 s处对应触发。这种以能量变化实现特征捡拾的方法对于地震弱信号适应性较强，捡拾效果高，适用于震相清晰、规则型的地震事件。
1.2.2貝叶斯（BIC）信息准则
当前有许多模型选择方法，如似然法、贝叶斯法和信息准则法等。1973年统计学家Hirotugu Akaike提出Akaike 信息准则（AIC）用来进行模型选择，对模型信息的发展产生了深远影响，Schwarz（1978）在此基础上提出了贝叶斯（BIC）信息准则。
最大似然准则是一种具有理论性的点估计法，此方法的基本思想是：从模型总体随机抽取n组样本观测值后，最合理的参数估计量以观测值出现的概率最大为最大估计准则。而模型选择问题即是从K个候选模型M={Mi；i=1，2，…，K}中选出一个最能表示给定N个数据X= {Xi；i=1，2，…，N}分布的模型选择问题。BIC准则是一种基于模型复杂度惩罚机制的渐进最大似然准则，用来解决模型选择问题（薛昊，2010）。
假设数据集X对于每个模型Mj的最大似然式为L（x，Mj），kj为模型Mj的参数个数，那么模型Mj的BIC计算式可用（3）式表示：BIC（Mj）=lgL（x，Mj）-12λkjlg（N）（3）式中：对数似然函数lgL（x，Mj）反映了模型精确度，它可以测量效率参数化模型的预测数据，也可用来选择集群数量内出现的一个特定数据集；λ为惩罚因子；λkjlg（N）反映了对模型复杂度的惩罚，它可提高模型的灵活性，满足规则内的惩罚需求（薛昊，2010）。
由于样本数据的似然值会随着模型参数的增加而增加，单纯采用极大似然估计会造成模型参数过多、维数过大，所以我们在BIC准则中引入了对模型复杂度的惩罚。由（3）式使BIC计算取得最大值得模型，即为基于贝叶斯信息准则的最优模型： </k1
1.2.3BIC算法流程
在波形图中，并不是只会出现单一的分割点，很多时候，会有多个分割点出现，为了得到最优分割点，本文采用贝叶斯准则中顺序检测算法逐个找出分割点，最后得到分割点的最大值，确定P波震相的精确到时。具体流程如图2所示。
图2中，首先是初始化检测窗[a，b]，且a=1；b=1。这一步骤表明视窗一开始仅包含2个样本点a和b。接着开始进行BIC检测[a，b]窗内是否存在分割点，不存在分割点时，对检测窗窗长增加1，即令b=b+1，重新开始BIC检测分割点；存在分割点时，则移动检测窗，即a=i+1；b=a+1，i为分割点，读取该分割点值并按顺序逐步增加检测窗长度，直至整个完整波频段检测结束。
从整个算法的流程可看出，算法的时间复杂度和分割点数有关，虽然BIC的检测精度很高，但其在检测过程中计算量非常大。
从长短时平均（STA/LTA）方法的捡拾效果来看，其算法简单、稳定快速，主要缺点是捡拾到的P波到时往往比实际到时滞后。而贝叶斯（BIC）方法捡拾精度高、稳定可靠，但其捡拾P波到时需要对整个完整波频按顺序逐步检测分割点，最后才能判定出所有分割点的最大值，确定出P波到时。因此，如果采用STA/LTA粗略捡拾到P波到时后，缩小检测范围，再用BIC准则进行P波震相的精确拾取，比在整条地震记录图上直接应用BIC准则要更快，精度与直接应用BIC准则捡拾相同。
ΔBIC（i）>0时，表明BIC（M2）>BIC（M1），即模型M2优于模型M1，当ΔBIC（i）的取值为最大值，则认为这是最优分割点i，最终使ΔBIC（i）取得最大值得帧i为P波震相到时点。i=argmax1<i0ΔBIC（i）（11）使用BIC准则方法捡拾2014年云南景谷M6.6主震中景谷永平台站UD向记录的P波到时，获得结果如图3所示。图3a描述了在该加速度时程曲线中，P波到时相对于记录起始时刻是13.37 s；
人工捡拾景谷永平台UD向记录P波到时为13.40 s，STA/LTA自动捡拾结果为13.49 s，比人工捡拾结果滞后0.09 s；BIC自动捡拾结果为13.37 s，比人工捡拾结果提前0.03 s。
1.3STA/LTA+BIC双步骤捡拾的工作流程
STA/LTA+BIC双步骤捡拾法的工作流程如图4所示，强震仪的阈值产生触发事件时，实时数据的特征函数可直接采用（1）式计算特征函数值；STA/LTA粗略捡拾P波到时采用（2）式计算；BIC精确捡拾P波到时采用（11）式计算。首先根据特征函数计算结果，对记录波形展开P波粗略捡拾，应用STA/LTA长短时平均法检测或报警触发事件，当STA/LTA刻画出记录幅值的瞬时变化时，立刻判别P波阈值是否触发，若是STA/LTA计算结果超过常数阈值THR，则认为是地震事件，进入下一步精确捡拾；反之则是干扰事件，捡拾工作及时结束。本方法对地震事件P波的精确捡拾采用BIC自动捡拾法，其捡拾结果作为最终的P波到时成为地震预警的重要参数。
归纳STA/LTA+BIC双步骤捡拾法的基本步骤主要有2点：（1）先用STA/LTA法粗略捡拾P波到时；（2）以粗略到时点为基准，前后各推0.5 s，在此范围内使用BIC准则进行精确捡拾。
2实际算例分析
根据本文对震相捡拾方法的原理介绍与工作流程，对云南省数字强震动台网记录的17次ML≥5.0地震事件进行STA/LTA+BIC双步骤自动捡拾，捡拾记录均来自震中距在100 km范围内的强震动台站记录。
2.1数据选取
本文所用数据来自于云南省数字强震动台网获取的强震记录，包括了2008—2014年17次M5.0～7.0地震事件，涵盖了云南盈江、腾冲、彝良、洱源、景谷等地震多发区域。对P波震相而言，其在竖直向的幅值变化较大，我们对地震动记录进行初次人工筛选，剔除丢头、叠加等波形异常的强震记录，选取震中距在100 km以内、P波在竖直向相对突出的竖向原始加速度时程记录173个，并对数据进行预处理。筛选后的各地区的加速度时程记录数量分布如表1所示。
2.2实际捡拾情况分析
如表1所示，台网内选取这17个地震事件的实际记录，有9个震中距在100 km内的记录出现漏检，其它记录捡拾工作正常。分析原因发现，其中有4个记录人工捡拾也无法分辨P波时，2个记录中的地震波形有其他小震的干扰叠加，其余3个原因不明。
表2中，STA/LTA捡拾与人工捡拾的P波到时平均差值为0.698 s，标准差值为0.531 s；STA/LTA+BIC捡拾与人工捡拾的P波到时平均差值为0.091 s，标准差值为0.156 s。从地震系统的工作要求而言，STA/LTA捡拾结果误差相对较大，其精确度可达到地震烈度速报的工作要求；STA/LTA+BIC捡拾结果误差相对很小，其精确度能够达到地震预警系统的工作要求。
图5显示了STA/LTA自动捡拾与人工捡拾P波到时结果的偏差分布，可以看出，STA/LTA自动捡拾P波到時的偏差基本以正数为主，最大偏差将近2.0，其捡拾结果通常大于人工捡拾结果，即其对P波到时的自动捡拾往往存在滞后的现象。图6表示了STA/LTA+BIC自动捡拾与人工捡拾P波到时结果的偏差分布，可以看出，STA/LTA+BIC自动捡拾P波到时的偏差均匀分布在-0.5～0.5之间，以0～0.2之间的分布最广。该方法对P波到时的自动捡拾结果与人工捡拾结果吻合。综上所述，采用STA/LTA+BIC自动捡拾P波到时的结果与人工捡拾P波到时结果偏差很小，误差不超过0.5。
2.3震中距对P波到时准确度及时效性的影响
为了研究震中距是否会影响地震P波震相到时的准确程度，我们统计了人工判别、STA/LTA和STA/LTA+BIC自动判别3种捡拾结果，将记录台站的震中距范围扩大至0～160 km，具体变化如图7所示。由图7可见，采用STA/LTA+BIC自动图5STA/LTA自动捡拾与人工捡拾判别方法与人工捡拾方法的捡拾结果几乎是完全吻合的，而STA/LTA自动判别方法捡拾的结果有延时滞后现象产生。同样，将STA/LTA+BIC自动捡拾P波到时的偏差与震中距的变化分布绘制成图8，其结果与图6相符，无论震中距如何变化，采用STA/LTA+BIC自动捡拾P波到时的偏差都很小。由此可见，震中距远近对P波到时的准确程度没有任何影响，采用合适的捡拾方法才能提高捡拾结果的精确度。为了得到不同震中距下采用STA/LTA+BIC自动捡拾P波到时的速度，本文引入速度节奏来显示捡拾速度，即规定一个节奏值，每个数据捡拾的用时以规定节奏值相应的倍数表示，这样就能通过不同的节奏数值看出每个数据的捡拾用时长短，进而推断出P波传播于不同介质时的速度变化。如图9所示，采用STA/LTA+BIC自动捡拾P波到时速度节奏随震中距的变化是逐渐增大的。究其原因，可能是随着震中距逐渐增大，不同场地介质会对P波传播速度及方向产生改变，影响后续对P波震相识别的速度。 </i
綜上所述，震中距的远近对P波震相识别的精度没有影响，P波到时捡拾结果的准确程度取决于采用方法是否适合；但震中距的远近对P波震相识别用时长短有一定影响，随着震中距逐渐增大，自动捡拾P波初至的耗时也会呈逐渐增长趋势。
3讨论与结论
本文借鉴了国内外P波震相自动识别的相关研究，以云南地区的地震数据为基础，提出一套可实时处理P波震相的方法，即STA/LTA和贝叶斯BIC双步骤捡拾法。应用此方法对所选取的云南强震动台网观测记录进行P波自动精确识别，并与人工捡拾方法的实测结果进行比较，得到以下结论：
（1）STA/LTA自动捡拾与人工捡拾P波到时的偏差基本以正数为主，最大偏差将近2.0 s，其捡拾结果通常大于人工捡拾结果，即其对P波到时的自动捡拾往往存在滞后的现象。STA/LTA+BIC自动捡拾P波到时的偏差在0.0值上下均匀分布，主要分布在-0.2～0.2 s之间，最大偏差也超不过0.5 s，该方法其对P波到时的自动捡拾结果与人工捡拾结果接近吻合。
（2）长短时平均（STA/LTA）方法捡拾速度快；贝叶斯（BIC）方法精确度高。因此，在用STA/LTA粗略捡拾到P波后，再用BIC准则进行P波震相的精确拾取，从地震预警的实时性与时效性要求来看，效果比在整条地震记录图上直接应用BIC准则要好。
（3）通过大量数据实际验证，P波特征函数触发阈值设为10较合适。阈值小于10误判率会增高，阈值大于10漏检率也会增高。进行P波到时捡拾时，信噪比高捡拾结果理想，若是其中有噪声等干扰信号影响，捡拾质量会受到很大影响。
（4）采用P波震相自动捡拾方法，结果虽然具有较高的精度，但对S波震相自动捡拾并未尝试，可在下一步工作中继续研究S波震相自动捡拾。
通过对比分析，STA/LTA与BIC相结合的P波到时捡拾方法不但利用了STA/LTA稳定可靠的特点，在初步捡拾点的附近选取一个恰当的时间窗包含P波震相初至过程，再以这个合适的时间窗进一步开展BIC准则精确捡拾，捡拾结果稳定可靠，准确性高。STA/LTA和贝叶斯BIC双步骤捡拾法的识别精度能够达到地震预警的要求，可作为云南地区预警系统建设的技术参考。
参考文献：
陈少波.2015.地震预警系统的P波震相自动识别方法研究[D].哈尔滨：中国地震局工程力学研究所.
高淑芳，李山有，武东坡，等.2008.一种改进的STA/LTA震相自动识别方法[J].世界地震工程，24（2）：37-41.
刘建华，刘福田，胥颐.2006.三分量地震资料的偏振分析[J].地球物理学进展，21（1）：6-10.
刘希强，周蕙兰，李红.2000a.基于小波包变换的地震数据时频分析方法[J].西北地震学报，22（2）：143-146.
刘希强，周蕙兰，沈萍，等.2000b.用于三分向记录震相识别的小波变换方法[J].地震学报，32（2）：125-131.
刘希强，周蕙兰，郑治真，等.1998.基于小波包变换的弱震相识别方法[J].地震学报，20（4）：373-380.
马强.2008.地震预警技术研究及应用[D].哈尔滨：中国地震局工程力学研究所.
毛燕，崔建文，郑定昌，等.2011.地震记录的P波自动捡拾[J].地震研究，34（1）：47-51.
史洪山.2012.地震预警震级测定技术研究[D].哈尔滨：中国地震局工程力学研究所.
王继，陈九辉，刘启元，等.2006.流动地震台站观测初至震相的自动检测[J].地震学报，28（1）：42-51.
薛昊.2010.基于BIC的通用音频分割方法研究[D].哈尔滨：哈尔滨工业大学.
尹得余.2012.P波震相自动捡拾与地震预警震级实时测定技术研究[D].哈尔滨：中国地震局工程力学研究所.
周彦文.2008.基于单台P波记录的早期地震预警方法研究[D].兰州：中国地震局兰州地震研究所.
ALLEN R E.1978.Automatic earthquake recognition and timing from single traces[J].Bulletin of the Seismological Society of America，68（5）：1521-1532.
ALLEN R E.1982.Automatic phase pickers：Their present use and future prospects[J].Bulletin of the Seismological Society of America，72（6B）：225-242.
AMBUTER B P，SOLOMON.1974.An event-recording system for monitoring small earthquakes[J].Bulletin of Seismological Society of America，64（4）：1181-1188.
ANDERSON K R.1981.Epicentral location using arrival time order[J].Bulletin of the Seismological Society of America，70（2）：541-545.
BAER M，KRADOLFER U.1987.An Automatic phase picker for local and teleseismic[J].Bulletin of the Seismological Society of America，72（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[J].Bulletin of Seismological Society of America，90（1）：187-198.
CHAKRABORT Y A，OKAYA D.1995.Frequency-time decomposition of seismic data using wavelet-based methods[J].Geophysics，60（6）：1906-1916.
TONG C，KENNETT B L N.1996.Auto seismic event recognition and later phase identification for broadband seismograms[J].Bulletin of the Seismological Society of America，86（6）：1896-1909.
CICHOWICZ A R.1993.An automatic S-phase picker[J].Bulletin of the Seismological Society of America，83（1）：180-189.
DAI H C，MACBETH C.1995.Automatic picking of seismic arrivals in local earthquake data using an artificial neural network[J].Geophysical Journal International，120（3）：758-774.
DAI H C，MACBETH C.1997.The application of back-propagation neural network to automatic picking seismic arrivals from single-component recordings[J].J Geophys Res，102（B7）：15105-15113.
EARLA P.1999.Polarization of the Earths teleseismic wave field，Geophys[J].J Int，139（1）：1-8.
EARLE P S，SHEARER P M.1994.Characterization of glob seismograms using an automatic-picking algorithm[J].Bulletin of the Seismological Society of America，84（2）：366-376.
FLINN E.1965.Signal analysis using recti linearity and direction of particle motion[J].Proc IEEE 53：1874-1876.
JUKEVICS A.1988.Polarization analysis of three-component array data[J].Bulletin of the Seismological Society of America，78（5）：1725-1743.
KANWALDIP A S，FARID U D.1997.Wavelet transform methods for phase identification in three-component seismograms[J].Bull Seism Soc Amer，87（6）：1598-1612.
LEONARD M，KENNETT B L N.1999.Multi-component autoregressive techniques for the analysis of seismograms[J].Phys Earth Planet Interiors，113（1-4）：247-264.
LEONARD M.2000.Comparison of manual and automatic onset time picking[J].Bull Seism Soc Am，90：1384-1390.
MAEDA N.1985.A method for reading and checking phase times in auto-processing system of seismic wave data[J].The Seismological Society of Japan，38：365-379.
MAO W，GUBBINS D.1995.Simultaneous determination of time delays and stacking weights in seismic array beam forming[J].Geophysics，60（2）：491-502.
MCEVILLY T V，MAJER E L.1982.An automated seismic processor for micro earthquake network[J].Bulletin of the Seismological Society of America，72（1）：303-325.
MURAT M E，RUDMAN A J.1992.Automated first arrival picking：a neural network approach[J].Geophysics Prospecting，40（6）：587-604.
SARAGIOTIS C H D，HADJILEONTIADIS L J，PANAS S M P.2002.AI-S/K：a robust automatic seismic P-phase arrival identification scheme[J].IEEE Sransactions on Geoscience and Remote Sensing，40（6）：1395-1404.
SCHWARZ E G.1978.Estimating the dimension of a model[J].Annals of Statistics，6（2）：461-464.
SLEEMAN R.1999.Robust automatic P-phase picking：an on-line implementation in the analysis of broadband Seismogram recordings[J].Phys Earth Planet Interiors，113：265-275.
STEVENSON P S.1976.Micro earthquakes at flathead lake，Montana：a study using automatic earthquake processing[J].Bulletin of the Seismological Society of America，66（1）：1-80.
SWINDELL H，SNELL N.1977.Station Processor Automatic Signal Detection System，Phase Ι：Final Report，Station Processor Software Development[J].Texas Instruments Report，ALEX（01）-FR-77-01.
VERNON F L.1987.Frequency dependent polarizatianalysis of high-frequency seismograms[J].Journal of geophysical research，92（B12）：664-674.
VIDALE T.1986.Complex polarization analysis of particle motion[J].Bull Seism Soc Am，76（5）：1393-1405.
WANG J，TENG T L.1995.Artificial neural network-based seismic detector[J].Bulletin of the Seismological Society of America，85（1）：308-319.
ZHAO Y，TAKANO K.1999.An artificial neural network approach for broadband seismic phase Picking[J].Bulletin of the Seismological Society of America，89（3）：70-680.
ZHANG H，THURBER C，ROWE C.2003.Automatic P-Wave Arriva Detection and Picking with Multiscale Wavelet Analysis for Single-Component Recordings[J].Bulletin of the Seismological Society of America，93（5）：1904-1912.
