迭代累积平方和算法与自回归赤池准则算法在地震信号到时估算中的比较研究

何 燕 靳 平 肖卫国 王红春

（中国西安710024西北核技术研究所）

引言

可靠、高效的地震信号自动处理技术对于日常地震监测以及灾后地震速报具有非常重要的意义（张诚鎏，2010），它不仅可以提高庞大地震数据的处理效率，减轻分析人员工作量，同时可以提高检测定位地震事件的时效性和完整性，从而提高整个台网系统的监测能力以及快速响应能力，而其中的震相自动提取技术是自动处理过程中不可或缺的关键技术环节之一.

在地震信号自动处理技术中，震相的自动提取主要包含两方面的内容：一是信号检测（初步的震相初至估计）；二是到时估算（精确的震相初至估计）.针对这两个方面，国内外诸多学者进行了不同程度的研究，其中信号检测方面主要以长短时平均能量比检测法STA／LTA（Stevenson，1976；Allen，1978，1982）为主，其它的检测方法还包括F检测法（Melton，Bailey，1957；Blandford，1974；唐明帅，王海涛，2011）、互相关检测法（Gibbons，Frode，2006）、偏振检测法（Kedrov，Ovtchinnikon，1990）、高阶统计量检测法（Longbottomet al，1988；Gentili，Bragato，2006）、人工神经网络检测法（Wang，Teng，1997；Zhao，Kiyoshi，1999）、卓越周期Tpd检测法（Hildyard et al，2008）等.到时估算方面则主要以赤池信息准则（Akaike’s information criterion，简写为AIC）算法（Akaike，1973）为基础.例如，自回归赤池准则算法（AR-AIC）（Leonard，Kennett，1999；王海军等，2003；杨配新等，2004；王继等，2006），无需计算AR系数的VRA-AIC方法（Maeda，1985），基于三阶统计量的TOC-AIC方法（刘希强等，2009）；另外还有基于神经网络及基于小波变换的AIC方法（Sleeman，van Eck，1999；Zhang，Clifford，2003）.其中，ARAIC方法由于具有较高精确度的特性（Küperkoch et al，2012），是目前国内外普遍采用的到时估算方法，但这种算法存在着运算速度较慢以及对于常见的S，Lg震相的估算误差相对较大的问题.本文将介绍不同于AR-AIC方法的另外一种到时估算方法——累积和法.虽然目前基于这种方法所开展的信号到时研究不是到时估算技术的主要研究方向，但在Inclán和Tiao（1994）、Der等（1998）以及Der和Shumway（1999）等人的研究基础上，我们发现在震相初至的精确估算方面，累积和算法在一定程度上可与AR-AIC算法相媲美.

本文重点介绍迭代累积平方和算法（iterative cumulative sums of squares algorithm，简写为ICSS）的基本原理（Inclán，Tiao，1994）及应用实例，并将其与AR-AIC算法的到时估算结果进行分析比较.

1 ICSS算法简述

1.1 算法起源

累积和（cumulative sum，简写为CUSUM）是一种序贯分析方法.该算法主要用于某个相对稳定的数据序列中，以检测开始发生异常的数据点.所谓异常的数据点是指从某点开始，整个数列的平均值或者均方差开始发生改变，进而影响到整组数据稳定的点.该算法假设时间数列起初呈一定稳定状态直至数列中的数据点突然发生改变为止，其基本原理如下.

令时间数列xt服从一个均值为0、方差为σ2的正态分布，其长度为T，为了估计这一长度区间内异常数据点的位置，须计算时间数列xt的累积平方和，即

并定义统计量Dk

当k为异常数据点所在位置时，Dk具有极小值.如果在所关注的区间内只有单个异常数据点时，上述算法能够得到较为准确的结果；但当区间内同时存在多个异常数据点时，潜在的掩蔽效应使Dk不具有充分性，结果往往存在较大误差.为解决这一问题，Inclán和Tiao（1994）提出利用改进后的累积和算法（即ICSS）来系统寻找该数列在不同区间内异常数据点的位置.

1.2 ICSS算法基本原理

ICSS算法具体来讲就是将时间长度T划分为多个不同的样本区间，即异常数据点的区间集合可表示为1＜t1＜t2…＜tNT＜T，NT表示从T个观察值中所检测到的数据点发生异常改变的总数.假设在所估计的样本区间内数据点没有发生异常改变，那么统计量Dk将在水平横轴为0的轴线附近随机波动；而如果该数列在所观测的区间内有一个或两个数据点发生异常改变时，Dk的值则会由0开始增加或是减少.在均值为0的假设条件下，根据Dk的分布可以导出临界值来检测数据点是否存在显著的改变.而当Dk绝对值的最大值大于临界值时，0假设条件也随之不成立，也即当｜Dk｜取得最大值时，k所在的位置即为发生异常改变的时间点.根据Inclán和Tiao（1994）、Der等（1998）以及Der和Shumway（1999）的描述并结合具体的到时估算窗口，ICSS算法大致可分为以下几个步骤：

1）令到时估算窗口取为［t1，t2］，由公式（1）可知这一时间段内的累积平方和函数为

再由公式（2）得出经规范化后的迭代累积和统计量为

理论上，当t为实际的信号初动时，D（t1，t）取得最大值.

2）截取一段事件记录的原始波形并作滤波处理.

3）将事件持续时间范围内的起止时间［t1，tn］作为到时估算的时间窗口.

4）计算这段时间长度内的迭代累积和统计量.

5）取起止时间［t1，t2］，计算D（t1，t2）并找出最大值所对应的时间点.

6）同样找出位于区间［t2，tn］范围内的D（t2，tn）最大值所对应的时间点.

7）重复上述过程直至没有新的时间点被找出.

由上述过程可以看出，运用ICSS算法估算信号到时需经过多次迭代计算，这是由于为了精确寻找在事件持续时间范围内的一个或多个震相到时，需要在不同区间内反复查找.下面将结合一些具体的例子来说明如何运用ICSS算法估算信号到时.

2 运用ICSS算法估算信号到时

选取发生在新疆地区的一次区域地震作为分析实例.首先将台站所记录到的地震信号进行滤波处理，然后截取一段波形数据进行到时估算.经过多次对比计算发现，当时间窗的起止时间以及滤波器的频带满足一定条件时，在事件持续时间范围内的所有震相到时（如P波、S波）均能被正确检出，且只需较少的迭代计算次数.而若改变估算的起止时间或者滤波器的频带范围，则有可能需要更多的迭代计算才能最终得到准确的结果.限于篇幅，本文仅以其中一个台站上的到时估算结果为例进行分析.

图1a，b分别给出了台站记录的原始地震波形及经滤波器相关频带滤波后的地震波形，图1c－－f则分别对应用ICSS算法经4次迭代后的累积和统计量曲线.可以看出图1a中短红线所在的位置为经人工分析后所确定的正确到时位置（包括Pn，Sn与Lg），虚线所在的位置为统计量最大值所对应的时间点（即实际估算时的信号到时位置）.从图1c所展现的一次迭代结果不难看出，虚线所标出的位置即是Sn波的到时所在.以此类推，图1d与图1f的虚线位置则分别对应Lg波和Pn波的到时.从图中可以看出，ICSS算法所估算出的3个到时位置都是比较准确的，与人工分析所确定的到时基本一致（虚线与短红线基本重合）.图2所展示的是改变了到时估算窗口的起止时间及滤波器频带后的到时估算结果.可明显看出经1次迭代后Sn波的到时仍较为准确，而经过5次迭代计算后，Lg波和Pn波的到时位置与正确位置（用短红线标出）相比却显得误差较大.实验证明，若在5次迭代基础上继续进行重复计算，最终会得到准确结果（限于篇幅，5次迭代后的迭代曲线不在图上画出），但这无疑将使计算过程变得十分繁琐且耗时.

图1 某区域震记录及ICSS算法到时估算结果（a）原始地震信号波形；（b）滤波后的地震信号波形；（c）经1次迭代后的ICSS曲线；（d）经2次迭代后的ICSS曲线；（e）经3次迭代后的ICSS曲线；（f）经4次迭代后的ICSS曲线Fig.1 A regional earthquake record and onset time estimation results by ICSS（a）Original seismic waveform；（b）Filtered seismic waveform；（c）ICSS curve after the first iteration；（d）ICSS curve after the second iteration；（e）ICSS curve after the third iteration；（f）ICSS curve after the fourth iteration

图2 改变估算条件后的某区域震记录及ICSS算法到时估算结果（a）滤波后的地震信号波形；（b）经1次迭代后的ICSS曲线；（c）经2次迭代后的ICSS曲线；（d）经3次迭代后的ICSS曲线；（e）经4次迭代后的ICSS曲线；（f）经5次迭代后的ICSS曲线Fig.2 A regional earthquake record and onset time estimation results by ICSS after the estimation conditions are modified（a）Filtered seismic waveform；（b）ICSS curve after the first iteration；（c）ICSS curve after the second iteration；（d）ICSS curve after the third iteration；（e）ICSS curve after the fourth iteration；（f）ICSS curve after the fifth iteration

由以上实例可以看出，运用ICSS算法估算信号到时的传统做法是在事件的持续时间范围内对这一时间段内的所有震相进行到时估算，这样做虽可一次性估算出所有信号到时，但其结果受滤波器频带及所选时间窗起止时间影响较大，有时为了得到更为准确的到时往往需要重复多次迭代计算，使程序执行效率变得低下.而在目前常用的地震信号自动处理技术中，其通常做法是在进行到时估算前通过STA／LTA信号检测算法得出初步的信号触发位置，并在此基础上进一步精确估算信号到时.由此我们认为，若能结合地震信号自动处理技术在进行到时估算前用到这一处理方法，那么不仅可以弥补ICSS算法在频带和时间窗选取方面的不足，同时也大大缩短了估算时间，提高了结果的准确性.有了初步的信号触发位置作为基础，相应地我们不再需要以一段事件波形的起始和结束作为到时估算的时间窗口，而只需以触发检测的时间点为中心来进行时间窗口的选取，这样做的好处可使迭代计算次数大为减少.图3给出了这两种时间窗口的对比示意图，具体的时间窗口长度可根据需要进行调整.

由图3可见，在结合地震信号自动处理方法后，传统意义上的ICSS算法所选取的到时估算窗口被分割为3个独立的部分，它们分别用来估算P波、S波以及Lg波的到时.图4给出了以上提到的区域震相（包括Pn，Sn和Lg震相）经独立的时间窗口估算后的到时结果，可以清晰地看出原本需经多次迭代计算才能得到的到时结果，现在只需1次迭代便可做到，大大提高了程序的执行效率.

图3 到时估算时间窗选取对比Fig.3 Comparison of time windows selection for onset time estimation

图4 结合地震信号自动处理方法后的ICSS算法到时估算结果（a）原始地震信号波形；（b）P波估算的ICSS曲线；（c）S波估算的ICSS曲线；（d）Lg波估算的ICSS曲线Fig.4 The onset time estimation results by ICSS combined with automatic seismic signal processing method（a）Original seismic waveform；（b）ICSS curve of P wave onset time estimation；（c）ICSS curve of S wave onset time estimation；（d）ICSS curve of Lg wave onset time estimation

3 ICSS算法与AR-AIC算法的估算结果对比及分析

AR-AIC算法经过长期的实践证明是比较准确和有效的，而为了验证ICSS算法估算信号到时的准确率到底如何，有必要将其与AR-AIC算法作一对比分析.首先对AR-AIC算法作一简要介绍.

3.1 AR-AIC算法简介

AR-AIC算法是在信号触发检测的基础上，以检测的触发时间为中心截取一段数据作为到时估算窗口，假定窗口的前若干秒为噪声，最后若干秒为信号，分别计算噪声和信号J阶的AR模型.假定信号初动到时对应第k个点，用噪声和信号的AR模型分别拟合前k－1个点和后M－k＋1个点，计算拟合误差，并计算相应的AIC值，即

式中，¯en，¯es分别为拟合误差序列en（i）和es（i）的均方根值.理论上，当k为实际的信号初动时AIC（k）取最小值.

当信噪比较大、信号初动比较明显时，该算法可以较好地确定信号初动；但当信噪比很低、信号初动也不明显时，则可能造成较大的误差.为此本文实际采用的是经过改进后的AR-AIC算法（王海军等，2003）.

3.2 ICSS算法改进

在应用ICSS算法估算信号到时时发现，当采用绝对幅值取代平方幅值进行累积和的计算时，有时可以得到更为准确的结果且计算时间也更短（Der，Shumway，1999）.即原有的累积平方和函数累积和统计量的公式则保持不变.

需要说明的是，由于通过STA／LTA检测得到的结果与实际到时相比通常偏晚，因此实际应用时最多进行两次迭代计算即可得到准确的结果.

3.3 估算结果对比与分析

选取新疆数字台网记录到的330个震相记录作为样本数据（包括P，S和Lg震相），在同等估算条件下（相同的触发检测时间及同样的滤波通道和频带）分别采用两种估算方法对各震相进行到时估算.两种方法所得到的到时误差统计结果见表1，到时误差分布情况见图5.这些记录涉及18个三分向地震台站以及20次近震和区域震事件，具体的台站分布及震中位置如图6所示.所选震相记录的震中距范围在0°—12°之间，其中位于2°—10°的区域震震相占绝大多数.此外所选地震事件中震级ML＜3的有8次，ML3—4之间的有9次，ML＞4的有3次，最高为ML4.7.

表1 ICSS和AR-AIC两种算法所得到的震相到时估算误差结果统计Table 1 Statistic results of the errors in onset time estimation using ICSS and AR-AIC methods

图5 ICSS（a）和AR-AIC（b）算法所得到的震相到时估算误差分布以及各震相（P，S，Lg）的误差分布对比情况（c）Fig.5 Distribution of the errors in onset time estimation using ICSS（a）and AR-AIC（b）methods as well as comparison of error distributions for onset time estimation of all phases（including P，S，Lg）（c）

图6 台站及震中位置分布Fig.6 Location of stations（triangles）and epicenters（circles）

由图5可见，ICSS与AR-AIC两种算法所给出的到时估算误差相对都不大，基本都在5s以内，且两种算法均呈现P震相误差较小，而S和Lg震相误差较大的情况.进一步结合表1中的统计结果进行分析发现，对于P震相，虽然两种算法估算误差小于2s的震相数非常接近，但在误差小于0.8s的震相数方面，可明显看出AR-AIC算法所占比例要大于ICSS算法，由此可见AR-AIC算法估算P震相到时的精度要优于ICSS算法，而在中位数方面，前者也明显优于后者；对于S和Lg震相，从表1及图5的统计情况可明显看出结果正好相反，ICSS算法的各项指标均优于AR-AIC算法.图7为到时误差与信噪比关系图.由图7可见，P震相的信噪比明显高于S及Lg震相的信噪比，且P震相的信噪比主要分布于1—20之间，而S及Lg震相的信噪比则主要分布于1—5之间.通过分析误差大小与信噪比高低之间的变化关系可见，两种算法均呈现到时误差随信噪比增大而减小的情况.对于P震相，两种算法无明显差别；对于S及Lg震相，可看出当信噪比大于5dB时，两者之间差别不大；而当信噪比小于5dB时，AR-AIC算法的到时估算误差比ICSS算法明显偏大，可见ICSS算法在估算低信噪比信号方面具有一定的优势.

图7 ICSS（a）和AR-AIC（b）算法所得到的震相到时估算误差与信噪比关系Fig.7 Relationship between errors in onset time estimation and SNRs by using ICSS（a）and AR-AIC（b）methods

此外，在对比过程中我们发现，ICSS算法的运算速度要优于AR-AIC算法.以处理100个震相为例，ICSS算法所需的运算时间比AR-AIC算法要快大约36s.由此可见当需要处理大量的震相数据时，ICSS算法的速度优势将更为明显.

4 讨论与结论

将迭代累积平方和算法应用于地震信号到时估算的传统做法是在事件的持续时间范围内对这一时间段内的所有震相进行到时估算，这样做虽可一次性估算出所有信号到时，但其结果受滤波器频带及所选时间窗的起止时间影响较大，为了得到更为准确的到时结果往往需要重复多次迭代计算，使程序执行效率变得低下.为了克服这一缺点，本文提出将ICSS算法与地震信号自动处理方法相结合，并将改进后的算法与目前所普遍采用的到时估算方法——AR-AIC算法进行了估算结果对比，初步结果表明：

1）就总的误差分布范围而言，两种算法相差不大，基本都在5s以内.对于P震相，ICSS算法在估算误差小于2s的震相数方面与AR-AIC算法非常接近，但在到时精度等方面则不如AR-AIC算法；而对于S和Lg震相，ICSS算法的各项指标均优于AR-AIC算法.

2）对于信噪比较高的信号，两种算法无太大差别；但对于信噪比较低的信号，ICSS算法所得到的S和Lg震相的到时误差明显小于AR-AIC算法.

3）ICSS算法的运算速度优于AR-AIC算法.

综上所述，我们认为与AR-AIC算法相比，ICSS算法在对区域震S和Lg震相的到时估算方面具有估算误差小、计算时间短的优势，具有一定的应用价值.

尽管运用ICSS算法估算信号到时具备一定的实用价值，但这只是在有限的数据测试得出的结果，并未接受大量数据的长期检验；同时本文所选震相多为区域震震相，并未涵盖所有震相，且ICSS算法与其它的AIC算法相比有可能会有不同的结果，本文只是选用了目前所普遍采用的AR-AIC到时估算方法来验证ICSS算法的有效性，其不足之处，有待我们进一步研究和改进.

新疆地震局为本文提供了台站波形数据；西北核技术研究所刘文学副研究员对本文修改给予了很大帮助.作者在此一并表示感谢.

刘希强，周彦文，曲均浩，石玉燕，李铂.2009.应用单台垂向记录进行区域地震事件实时检测和直达P波初动自动识别［J］.地震学报，31（3）：260－－271.

Liu X Q，Zhou Y W，Qu J H，Shi Y Y，Li B.2009.Real-time detection of regional events and automatic P-phase identification from the vertical component of a single station record［J］.Acta Seismologica Sinica，31（3）：260－－271（in Chinese）.

唐明帅，王海涛.2011.F检测算法及其在识别地震地表反射波震相中的初步应用［J］.地震学报，33（6）：776－－787.

Tang M S，Wang H T.2011.F-detection algorithm and its preliminary application to seismic surface reflected wave phase identification［J］.Acta Seismologica Sinica，33（6）：776－－787（in Chinese）.

王海军，靳平，刘贵忠.2003.低信噪比地震记录中信号初至的估计［J］.西安交通大学学报，37（6）：659－－660.

Wang H J，Jin P，Liu G Z.2003.Onset estimation for signals with low signal to noise ratio［J］.Journal of Xi’an Jiao tong University，37（6）：659－－660（in Chinese）.

王继，陈九辉，刘启元，李顺成，郭飚.2006.流动地震台阵观测初至震相的自动检测［J］.地震学报，28（1）：42－－51.

Wang J，Chen J H，Liu Q Y，Li S C，Guo B.2006.Automatic onset phase picking for portable seismic array observation［J］.Acta Seismologica Sinica，28（1）：42－－51（in Chinese）.

杨配新，邓存华，刘希强，任勇，颜其中.2004.数字化地震记录震相自动识别的方法研究［J］.地震研究，27（4）：308－－313.

Yang P X，Deng C H，Liu X Q，Ren Y，Yan Q Z.2004.Studies on auto-distinguishing phase of digital seismic recordings［J］.Journal of Seismological Research，27（4）：308－－313（in Chinese）.

张诚鎏.2010.地震信号识别和关联技术研究［D］.西安：西北核技术研究所：10－－15.

Zhang C L.2010.Study on seismic phases identification and association［D］.Xi’an：Northwest Institute of Nuclear Technology：10－－15（in Chinese）.

Akaike H.1973.Information theory and an extension of the maximum likelihood principle［C］∥Petrov B N，Csaki F eds.International Symposium on Information Theory.Budapest：Akademiai Kiado：267－－281.

Allen R V.1978.Automatic earthquake recognition and timing from single trace［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）：225－－242.

Blandford R R.1974.An automatic event detector at the Tonto Forest seismic observatory［J］.Geophysics，39（5）：633－－643.Der Z A，Shumway R H，McGarvey M W.1998.Automatic interpretation of regional seismic signals using the CUSUMSA algorithms［G］∥21st Seismic Research Symposium.Davis：University of California：393－－403.

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／2／3／4）：227－－246.

Gentili S，Bragato P.2006.A neural-tree-based system for automatic location of earthquakes in Northeastern Italy［J］.Seismology，10（1）：73－－89.

Gibbons S J，Frode R.2006.The detection of low magnitude seismic events using array-based waveform correlation［J］.Geophys J Int，165（1）：149－－166.

Hildyard M W，Nippress S E J，Rietbrok A.2008.Event detection and phase picking using a time-domain estimate of predominate period Tpd［J］.Bull Seismol Soc Am，98（6）：3025－－3032.

Inclán C，Tiao G C.1994.Use of cumulative sums of squares for retrospective detection of changes of variance［J］.Journal of the American Statistical Association，89（427）：913－－923.

Kedrov O K，Ovtchinnikov V M.1990.An online analysis system for three-component seismic data：Method and preliminary results［J］.Bull Seismol Soc Am，80（6）：2053－－2071.

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：23－－30.

Leonard M，Kennett B L N.1999.Multi-component autoregressive techniques for the analysis of seismograms［J］.Phys Earth Planet Inter，113（1／2／3／4）：247－－263.

Longbottom J，Walden A T，White R E.1988.Principles and application of maximum kurtosis phase estimation［J］.Geophysical Prospecting，36（2）：115－－138.

Maeda N.1985.A method for reading and checking phase times in autoprocessing system of seismic data［J］.Zisin，38（3）：365－－379.

Melton R S，Bailey L F.1957.Multi signal correlators［J］.Gepphysics，22（3）：565－－588.

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／2／3／4）：265－－275.

Stevenson R.1976.Microearthquakes at Flathead Lake，Montana：A study using automatic earthquake processing［J］.Bull Seismol Soc Am，66（1）：61－－79.

Wang J，Teng T.1997.Identification and picking of S phase using an artificial neural network［J］.Bull Seismol Soc Am，87（5）：1140－－1149.

Zhang H，Clifford T.2003.Automatic P-wave arrival detection and picking with multiscale wavelet analysis for single component recordings［J］.Bull Seismol Soc Am，93（5）：1904－－1912.

Zhao Y，Kiyoshi T.1999.An artificial neural network approach for broadband seismic phase picking［J］.Bull Seismol Soc Am，89（3）：670－－680.