何 燕 靳 平 肖衛(wèi)國 王紅春

（中國西安710024西北核技術(shù)研究所）

引言

可靠、高效的地震信號自動處理技術(shù)對于日常地震監(jiān)測以及災(zāi)后地震速報具有非常重要的意義（張誠鎏，2010），它不僅可以提高龐大地震數(shù)據(jù)的處理效率，減輕分析人員工作量，同時可以提高檢測定位地震事件的時效性和完整性，從而提高整個臺網(wǎng)系統(tǒng)的監(jiān)測能力以及快速響應(yīng)能力，而其中的震相自動提取技術(shù)是自動處理過程中不可或缺的關(guān)鍵技術(shù)環(huán)節(jié)之一.

在地震信號自動處理技術(shù)中，震相的自動提取主要包含兩方面的內(nèi)容：一是信號檢測（初步的震相初至估計）；二是到時估算（精確的震相初至估計）.針對這兩個方面，國內(nèi)外諸多學(xué)者進行了不同程度的研究，其中信號檢測方面主要以長短時平均能量比檢測法STA／LTA（Stevenson，1976；Allen，1978，1982）為主，其它的檢測方法還包括F檢測法（Melton，Bailey，1957；Blandford，1974；唐明帥，王海濤，2011）、互相關(guān)檢測法（Gibbons，F(xiàn)rode，2006）、偏振檢測法（Kedrov，Ovtchinnikon，1990）、高階統(tǒng)計量檢測法（Longbottomet al，1988；Gentili，Bragato，2006）、人工神經(jīng)網(wǎng)絡(luò)檢測法（Wang，Teng，1997；Zhao，Kiyoshi，1999）、卓越周期Tpd檢測法（Hildyard et al，2008）等.到時估算方面則主要以赤池信息準(zhǔn)則（Akaike’s information criterion，簡寫為AIC）算法（Akaike，1973）為基礎(chǔ).例如，自回歸赤池準(zhǔn)則算法（AR-AIC）（Leonard，Kennett，1999；王海軍等，2003；楊配新等，2004；王繼等，2006），無需計算AR系數(shù)的VRA-AIC方法（Maeda，1985），基于三階統(tǒng)計量的TOC-AIC方法（劉希強等，2009）；另外還有基于神經(jīng)網(wǎng)絡(luò)及基于小波變換的AIC方法（Sleeman，van Eck，1999；Zhang，Clifford，2003）.其中，ARAIC方法由于具有較高精確度的特性（Küperkoch et al，2012），是目前國內(nèi)外普遍采用的到時估算方法，但這種算法存在著運算速度較慢以及對于常見的S，Lg震相的估算誤差相對較大的問題.本文將介紹不同于AR-AIC方法的另外一種到時估算方法——累積和法.雖然目前基于這種方法所開展的信號到時研究不是到時估算技術(shù)的主要研究方向，但在Inclán和Tiao（1994）、Der等（1998）以及Der和Shumway（1999）等人的研究基礎(chǔ)上，我們發(fā)現(xiàn)在震相初至的精確估算方面，累積和算法在一定程度上可與AR-AIC算法相媲美.

本文重點介紹迭代累積平方和算法（iterative cumulative sums of squares algorithm，簡寫為ICSS）的基本原理（Inclán，Tiao，1994）及應(yīng)用實例，并將其與AR-AIC算法的到時估算結(jié)果進行分析比較.

1 ICSS算法簡述

1.1 算法起源

累積和（cumulative sum，簡寫為CUSUM）是一種序貫分析方法.該算法主要用于某個相對穩(wěn)定的數(shù)據(jù)序列中，以檢測開始發(fā)生異常的數(shù)據(jù)點.所謂異常的數(shù)據(jù)點是指從某點開始，整個數(shù)列的平均值或者均方差開始發(fā)生改變，進而影響到整組數(shù)據(jù)穩(wěn)定的點.該算法假設(shè)時間數(shù)列起初呈一定穩(wěn)定狀態(tài)直至數(shù)列中的數(shù)據(jù)點突然發(fā)生改變?yōu)橹?，其基本原理如?

令時間數(shù)列xt服從一個均值為0、方差為σ2的正態(tài)分布，其長度為T，為了估計這一長度區(qū)間內(nèi)異常數(shù)據(jù)點的位置，須計算時間數(shù)列xt的累積平方和，即

并定義統(tǒng)計量Dk

當(dāng)k為異常數(shù)據(jù)點所在位置時，Dk具有極小值.如果在所關(guān)注的區(qū)間內(nèi)只有單個異常數(shù)據(jù)點時，上述算法能夠得到較為準(zhǔn)確的結(jié)果；但當(dāng)區(qū)間內(nèi)同時存在多個異常數(shù)據(jù)點時，潛在的掩蔽效應(yīng)使Dk不具有充分性，結(jié)果往往存在較大誤差.為解決這一問題，Inclán和Tiao（1994）提出利用改進后的累積和算法（即ICSS）來系統(tǒng)尋找該數(shù)列在不同區(qū)間內(nèi)異常數(shù)據(jù)點的位置.

1.2 ICSS算法基本原理

ICSS算法具體來講就是將時間長度T劃分為多個不同的樣本區(qū)間，即異常數(shù)據(jù)點的區(qū)間集合可表示為1＜t1＜t2…＜tNT＜T，NT表示從T個觀察值中所檢測到的數(shù)據(jù)點發(fā)生異常改變的總數(shù).假設(shè)在所估計的樣本區(qū)間內(nèi)數(shù)據(jù)點沒有發(fā)生異常改變，那么統(tǒng)計量Dk將在水平橫軸為0的軸線附近隨機波動；而如果該數(shù)列在所觀測的區(qū)間內(nèi)有一個或兩個數(shù)據(jù)點發(fā)生異常改變時，Dk的值則會由0開始增加或是減少.在均值為0的假設(shè)條件下，根據(jù)Dk的分布可以導(dǎo)出臨界值來檢測數(shù)據(jù)點是否存在顯著的改變.而當(dāng)Dk絕對值的最大值大于臨界值時，0假設(shè)條件也隨之不成立，也即當(dāng)｜Dk｜取得最大值時，k所在的位置即為發(fā)生異常改變的時間點.根據(jù)Inclán和Tiao（1994）、Der等（1998）以及Der和Shumway（1999）的描述并結(jié)合具體的到時估算窗口，ICSS算法大致可分為以下幾個步驟：

1）令到時估算窗口取為［t1，t2］，由公式（1）可知這一時間段內(nèi)的累積平方和函數(shù)為

再由公式（2）得出經(jīng)規(guī)范化后的迭代累積和統(tǒng)計量為

理論上，當(dāng)t為實際的信號初動時，D（t1，t）取得最大值.

2）截取一段事件記錄的原始波形并作濾波處理.

3）將事件持續(xù)時間范圍內(nèi)的起止時間［t1，tn］作為到時估算的時間窗口.

4）計算這段時間長度內(nèi)的迭代累積和統(tǒng)計量.

5）取起止時間［t1，t2］，計算D（t1，t2）并找出最大值所對應(yīng)的時間點.

6）同樣找出位于區(qū)間［t2，tn］范圍內(nèi)的D（t2，tn）最大值所對應(yīng)的時間點.

7）重復(fù)上述過程直至沒有新的時間點被找出.

由上述過程可以看出，運用ICSS算法估算信號到時需經(jīng)過多次迭代計算，這是由于為了精確尋找在事件持續(xù)時間范圍內(nèi)的一個或多個震相到時，需要在不同區(qū)間內(nèi)反復(fù)查找.下面將結(jié)合一些具體的例子來說明如何運用ICSS算法估算信號到時.

2 運用ICSS算法估算信號到時

選取發(fā)生在新疆地區(qū)的一次區(qū)域地震作為分析實例.首先將臺站所記錄到的地震信號進行濾波處理，然后截取一段波形數(shù)據(jù)進行到時估算.經(jīng)過多次對比計算發(fā)現(xiàn)，當(dāng)時間窗的起止時間以及濾波器的頻帶滿足一定條件時，在事件持續(xù)時間范圍內(nèi)的所有震相到時（如P波、S波）均能被正確檢出，且只需較少的迭代計算次數(shù).而若改變估算的起止時間或者濾波器的頻帶范圍，則有可能需要更多的迭代計算才能最終得到準(zhǔn)確的結(jié)果.限于篇幅，本文僅以其中一個臺站上的到時估算結(jié)果為例進行分析.

圖1a，b分別給出了臺站記錄的原始地震波形及經(jīng)濾波器相關(guān)頻帶濾波后的地震波形，圖1c－－f則分別對應(yīng)用ICSS算法經(jīng)4次迭代后的累積和統(tǒng)計量曲線.可以看出圖1a中短紅線所在的位置為經(jīng)人工分析后所確定的正確到時位置（包括Pn，Sn與Lg），虛線所在的位置為統(tǒng)計量最大值所對應(yīng)的時間點（即實際估算時的信號到時位置）.從圖1c所展現(xiàn)的一次迭代結(jié)果不難看出，虛線所標(biāo)出的位置即是Sn波的到時所在.以此類推，圖1d與圖1f的虛線位置則分別對應(yīng)Lg波和Pn波的到時.從圖中可以看出，ICSS算法所估算出的3個到時位置都是比較準(zhǔn)確的，與人工分析所確定的到時基本一致（虛線與短紅線基本重合）.圖2所展示的是改變了到時估算窗口的起止時間及濾波器頻帶后的到時估算結(jié)果.可明顯看出經(jīng)1次迭代后Sn波的到時仍較為準(zhǔn)確，而經(jīng)過5次迭代計算后，Lg波和Pn波的到時位置與正確位置（用短紅線標(biāo)出）相比卻顯得誤差較大.實驗證明，若在5次迭代基礎(chǔ)上繼續(xù)進行重復(fù)計算，最終會得到準(zhǔn)確結(jié)果（限于篇幅，5次迭代后的迭代曲線不在圖上畫出），但這無疑將使計算過程變得十分繁瑣且耗時.

圖1 某區(qū)域震記錄及ICSS算法到時估算結(jié)果（a）原始地震信號波形；（b）濾波后的地震信號波形；（c）經(jīng)1次迭代后的ICSS曲線；（d）經(jīng)2次迭代后的ICSS曲線；（e）經(jīng)3次迭代后的ICSS曲線；（f）經(jīng)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 改變估算條件后的某區(qū)域震記錄及ICSS算法到時估算結(jié)果（a）濾波后的地震信號波形；（b）經(jīng)1次迭代后的ICSS曲線；（c）經(jīng)2次迭代后的ICSS曲線；（d）經(jīng)3次迭代后的ICSS曲線；（e）經(jīng)4次迭代后的ICSS曲線；（f）經(jīng)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算法估算信號到時的傳統(tǒng)做法是在事件的持續(xù)時間范圍內(nèi)對這一時間段內(nèi)的所有震相進行到時估算，這樣做雖可一次性估算出所有信號到時，但其結(jié)果受濾波器頻帶及所選時間窗起止時間影響較大，有時為了得到更為準(zhǔn)確的到時往往需要重復(fù)多次迭代計算，使程序執(zhí)行效率變得低下.而在目前常用的地震信號自動處理技術(shù)中，其通常做法是在進行到時估算前通過STA／LTA信號檢測算法得出初步的信號觸發(fā)位置，并在此基礎(chǔ)上進一步精確估算信號到時.由此我們認為，若能結(jié)合地震信號自動處理技術(shù)在進行到時估算前用到這一處理方法，那么不僅可以彌補ICSS算法在頻帶和時間窗選取方面的不足，同時也大大縮短了估算時間，提高了結(jié)果的準(zhǔn)確性.有了初步的信號觸發(fā)位置作為基礎(chǔ)，相應(yīng)地我們不再需要以一段事件波形的起始和結(jié)束作為到時估算的時間窗口，而只需以觸發(fā)檢測的時間點為中心來進行時間窗口的選取，這樣做的好處可使迭代計算次數(shù)大為減少.圖3給出了這兩種時間窗口的對比示意圖，具體的時間窗口長度可根據(jù)需要進行調(diào)整.

由圖3可見，在結(jié)合地震信號自動處理方法后，傳統(tǒng)意義上的ICSS算法所選取的到時估算窗口被分割為3個獨立的部分，它們分別用來估算P波、S波以及Lg波的到時.圖4給出了以上提到的區(qū)域震相（包括Pn，Sn和Lg震相）經(jīng)獨立的時間窗口估算后的到時結(jié)果，可以清晰地看出原本需經(jīng)多次迭代計算才能得到的到時結(jié)果，現(xiàn)在只需1次迭代便可做到，大大提高了程序的執(zhí)行效率.

圖3 到時估算時間窗選取對比Fig.3 Comparison of time windows selection for onset time estimation

圖4 結(jié)合地震信號自動處理方法后的ICSS算法到時估算結(jié)果（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算法的估算結(jié)果對比及分析

AR-AIC算法經(jīng)過長期的實踐證明是比較準(zhǔn)確和有效的，而為了驗證ICSS算法估算信號到時的準(zhǔn)確率到底如何，有必要將其與AR-AIC算法作一對比分析.首先對AR-AIC算法作一簡要介紹.

3.1 AR-AIC算法簡介

AR-AIC算法是在信號觸發(fā)檢測的基礎(chǔ)上，以檢測的觸發(fā)時間為中心截取一段數(shù)據(jù)作為到時估算窗口，假定窗口的前若干秒為噪聲，最后若干秒為信號，分別計算噪聲和信號J階的AR模型.假定信號初動到時對應(yīng)第k個點，用噪聲和信號的AR模型分別擬合前k－1個點和后M－k＋1個點，計算擬合誤差，并計算相應(yīng)的AIC值，即

式中，ˉen，ˉes分別為擬合誤差序列en（i）和es（i）的均方根值.理論上，當(dāng)k為實際的信號初動時AIC（k）取最小值.

當(dāng)信噪比較大、信號初動比較明顯時，該算法可以較好地確定信號初動；但當(dāng)信噪比很低、信號初動也不明顯時，則可能造成較大的誤差.為此本文實際采用的是經(jīng)過改進后的AR-AIC算法（王海軍等，2003）.

3.2 ICSS算法改進

在應(yīng)用ICSS算法估算信號到時時發(fā)現(xiàn)，當(dāng)采用絕對幅值取代平方幅值進行累積和的計算時，有時可以得到更為準(zhǔn)確的結(jié)果且計算時間也更短（Der，Shumway，1999）.即原有的累積平方和函數(shù)累積和統(tǒng)計量的公式則保持不變.

需要說明的是，由于通過STA／LTA檢測得到的結(jié)果與實際到時相比通常偏晚，因此實際應(yīng)用時最多進行兩次迭代計算即可得到準(zhǔn)確的結(jié)果.

3.3 估算結(jié)果對比與分析

選取新疆?dāng)?shù)字臺網(wǎng)記錄到的330個震相記錄作為樣本數(shù)據(jù)（包括P，S和Lg震相），在同等估算條件下（相同的觸發(fā)檢測時間及同樣的濾波通道和頻帶）分別采用兩種估算方法對各震相進行到時估算.兩種方法所得到的到時誤差統(tǒng)計結(jié)果見表1，到時誤差分布情況見圖5.這些記錄涉及18個三分向地震臺站以及20次近震和區(qū)域震事件，具體的臺站分布及震中位置如圖6所示.所選震相記錄的震中距范圍在0°—12°之間，其中位于2°—10°的區(qū)域震震相占絕大多數(shù).此外所選地震事件中震級ML＜3的有8次，ML3—4之間的有9次，ML＞4的有3次，最高為ML4.7.

表1 ICSS和AR-AIC兩種算法所得到的震相到時估算誤差結(jié)果統(tǒng)計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以內(nèi)，且兩種算法均呈現(xiàn)P震相誤差較小，而S和Lg震相誤差較大的情況.進一步結(jié)合表1中的統(tǒng)計結(jié)果進行分析發(fā)現(xiàn)，對于P震相，雖然兩種算法估算誤差小于2s的震相數(shù)非常接近，但在誤差小于0.8s的震相數(shù)方面，可明顯看出AR-AIC算法所占比例要大于ICSS算法，由此可見AR-AIC算法估算P震相到時的精度要優(yōu)于ICSS算法，而在中位數(shù)方面，前者也明顯優(yōu)于后者；對于S和Lg震相，從表1及圖5的統(tǒng)計情況可明顯看出結(jié)果正好相反，ICSS算法的各項指標(biāo)均優(yōu)于AR-AIC算法.圖7為到時誤差與信噪比關(guān)系圖.由圖7可見，P震相的信噪比明顯高于S及Lg震相的信噪比，且P震相的信噪比主要分布于1—20之間，而S及Lg震相的信噪比則主要分布于1—5之間.通過分析誤差大小與信噪比高低之間的變化關(guān)系可見，兩種算法均呈現(xiàn)到時誤差隨信噪比增大而減小的情況.對于P震相，兩種算法無明顯差別；對于S及Lg震相，可看出當(dāng)信噪比大于5dB時，兩者之間差別不大；而當(dāng)信噪比小于5dB時，AR-AIC算法的到時估算誤差比ICSS算法明顯偏大，可見ICSS算法在估算低信噪比信號方面具有一定的優(yōu)勢.

圖7 ICSS（a）和AR-AIC（b）算法所得到的震相到時估算誤差與信噪比關(guān)系Fig.7 Relationship between errors in onset time estimation and SNRs by using ICSS（a）and AR-AIC（b）methods

此外，在對比過程中我們發(fā)現(xiàn)，ICSS算法的運算速度要優(yōu)于AR-AIC算法.以處理100個震相為例，ICSS算法所需的運算時間比AR-AIC算法要快大約36s.由此可見當(dāng)需要處理大量的震相數(shù)據(jù)時，ICSS算法的速度優(yōu)勢將更為明顯.

4 討論與結(jié)論

將迭代累積平方和算法應(yīng)用于地震信號到時估算的傳統(tǒng)做法是在事件的持續(xù)時間范圍內(nèi)對這一時間段內(nèi)的所有震相進行到時估算，這樣做雖可一次性估算出所有信號到時，但其結(jié)果受濾波器頻帶及所選時間窗的起止時間影響較大，為了得到更為準(zhǔn)確的到時結(jié)果往往需要重復(fù)多次迭代計算，使程序執(zhí)行效率變得低下.為了克服這一缺點，本文提出將ICSS算法與地震信號自動處理方法相結(jié)合，并將改進后的算法與目前所普遍采用的到時估算方法——AR-AIC算法進行了估算結(jié)果對比，初步結(jié)果表明：

1）就總的誤差分布范圍而言，兩種算法相差不大，基本都在5s以內(nèi).對于P震相，ICSS算法在估算誤差小于2s的震相數(shù)方面與AR-AIC算法非常接近，但在到時精度等方面則不如AR-AIC算法；而對于S和Lg震相，ICSS算法的各項指標(biāo)均優(yōu)于AR-AIC算法.

2）對于信噪比較高的信號，兩種算法無太大差別；但對于信噪比較低的信號，ICSS算法所得到的S和Lg震相的到時誤差明顯小于AR-AIC算法.

3）ICSS算法的運算速度優(yōu)于AR-AIC算法.

綜上所述，我們認為與AR-AIC算法相比，ICSS算法在對區(qū)域震S和Lg震相的到時估算方面具有估算誤差小、計算時間短的優(yōu)勢，具有一定的應(yīng)用價值.

盡管運用ICSS算法估算信號到時具備一定的實用價值，但這只是在有限的數(shù)據(jù)測試得出的結(jié)果，并未接受大量數(shù)據(jù)的長期檢驗；同時本文所選震相多為區(qū)域震震相，并未涵蓋所有震相，且ICSS算法與其它的AIC算法相比有可能會有不同的結(jié)果，本文只是選用了目前所普遍采用的AR-AIC到時估算方法來驗證ICSS算法的有效性，其不足之處，有待我們進一步研究和改進.

新疆地震局為本文提供了臺站波形數(shù)據(jù)；西北核技術(shù)研究所劉文學(xué)副研究員對本文修改給予了很大幫助.作者在此一并表示感謝.

劉希強，周彥文，曲均浩，石玉燕，李鉑.2009.應(yīng)用單臺垂向記錄進行區(qū)域地震事件實時檢測和直達P波初動自動識別［J］.地震學(xué)報，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檢測算法及其在識別地震地表反射波震相中的初步應(yīng)用［J］.地震學(xué)報，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］.西安交通大學(xué)學(xué)報，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］.地震學(xué)報，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.數(shù)字化地震記錄震相自動識別的方法研究［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.地震信號識別和關(guān)聯(lián)技術(shù)研究［D］.西安：西北核技術(shù)研究所：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，F(xiàn)rode 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.