李力 張建國 李盛樂

摘要：基于余震識別自動化對地震監(jiān)測、災(zāi)后應(yīng)急、確定發(fā)震構(gòu)造、分析預(yù)報和科學(xué)研究的重要意義，總結(jié)和分析了余震自動識別技術(shù)的現(xiàn)狀和未來發(fā)展。首先，總結(jié)單臺震相到時自動拾取方法及其特征函數(shù)，包括單特征方法、互相關(guān)方法和機器學(xué)習(xí)方法。其次，分別總結(jié)基于震相走時的常規(guī)識別方法和基于特征函數(shù)的偏移疊加方法。最后，對余震識別技術(shù)的未來發(fā)展進行分析和展望，認(rèn)為偏移疊加方法的性能改進將成為未來余震識別的重點發(fā)展方向，深度學(xué)習(xí)方法作為一種新技術(shù)，將在包括震相拾取和余震識別在內(nèi)的地震數(shù)據(jù)處理領(lǐng)域發(fā)揮重要作用。

關(guān)鍵詞：余震；自動識別；震相自動拾?。蛔邥r定位；偏移疊加

中圖分類號：P315.61 文獻(xiàn)標(biāo)識碼：A 文章編號：1000-0666（2018）01-0001-13

0 引言

破壞性地震發(fā)生之后，隨之而來的是大量的余震，從連續(xù)波形中識別余震是開展余震震級測定、震源機制解反演等后續(xù)監(jiān)測工作的基礎(chǔ)。余震信息是評估災(zāi)情和指導(dǎo)救援力量部署（張?zhí)K平等，2013；鄭韻，2015）、確定發(fā)震構(gòu)造和揭示隱伏斷層（黃媛等，2008；徐錫偉等，2008，2013a，b）、判定震后趨勢（莊建倉，馬麗，2000；宋金，蔣海昆，2009；余娜等，2016）的重要依據(jù)，也是震后地下速度結(jié)構(gòu)變化、余震遷移等科研工作的基礎(chǔ)數(shù)據(jù)。人工識別余震是一項十分耗時且枯燥的工作，結(jié)果中不可避免地會引人主觀因素，隨著觀測臺網(wǎng)的建設(shè)，數(shù)據(jù)的產(chǎn)出量越來越大，完全由人工處理變得越來越不現(xiàn)實。大震之后，識別高頻度的余震需要持續(xù)投入大量的時間和人力，這讓有限的監(jiān)測資源更加捉襟見肘，數(shù)據(jù)產(chǎn)出時效也受到限制。

余震識別的自動化有重要的應(yīng)用和科研價值，但是現(xiàn)有技術(shù)仍然有很多問題亟待解決和完善，地震監(jiān)測領(lǐng)域也沒有專門的余震自動識別系統(tǒng)。為了提高對余震自動識別技術(shù)的關(guān)注和投入，2017年5月中國地震局地球物理研究所和阿里云聯(lián)合舉辦余震捕捉AI大賽，吸引來自國內(nèi)外的近千支隊伍報名參賽，一時間余震自動識別成為地震領(lǐng)域的熱門話題和技術(shù)熱點。

余震發(fā)生后，其能量以地震波的形式向外傳播，當(dāng)震相到達(dá)各個觀測臺站時，會引起波形記錄的變化。自動識別余震時，首先需要從連續(xù)觀測波形中提取特征，并自動拾取震相到時，然后聯(lián)合多個臺站的數(shù)據(jù)確定余震信息。對于震相清晰的獨立地震，可以使用單臺波形進行識別，并初步確定震源信息。但是余震在時間上高度密集，并且震級普遍較小，記錄到的信號受主震尾波和其他余震的強烈干擾，多臺聯(lián)合是提高余震識別可靠性的必要技術(shù)。同時，余震在構(gòu)造上與主震有密切聯(lián)系，在空間上沿發(fā)震構(gòu)造密集展布，發(fā)震機理有較高的相似性，這對于震相識別和震相聯(lián)合是可以利用的有利條件。傳統(tǒng)震相拾取方法所使用的特征函數(shù)在偏移疊加和人工智能等新技術(shù)中也有非常廣泛的應(yīng)用，因此本文從單特征方法出發(fā)，總結(jié)幾種重要的震相自動拾取方法及其特征函數(shù)。然后介紹基于震相走時定位的常規(guī)方法的發(fā)展，并分析其在余震自動識別中的適用性，最后對基于特征函數(shù)的偏移疊加方法進行綜述，并在此基礎(chǔ)上展望余震自動識別技術(shù)的未來發(fā)展方向，以期對開展相關(guān)研究工作有所幫助。震相特征的提取是走時定位和偏移疊加方法的基礎(chǔ)，為便于理解，對STA/LTA、VAR-AIC、峰度、極性分析等單臺特征提取算法進行了簡單實現(xiàn)，相關(guān)示例代碼發(fā)布于github版本庫https：//github.com/dualli/trigger。

1 震相自動拾取方法及其特征函數(shù)

當(dāng)?shù)卣鹫鹣嗟竭_(dá)觀測臺站時，會改變記錄信號的能量、頻譜、極性和整體形態(tài)等特征，這些變化可以通過不同的特征函數(shù)提取，當(dāng)特征函數(shù)達(dá)到設(shè)定閾值時拾取相應(yīng)的震相到時。多個臺站的特征函數(shù)和震相信息成為后續(xù)地震識別和定位的基本數(shù)據(jù)。目前，震相自動拾取方法種類繁多，前人對這些方法按照時間域、頻率域、時頻域、綜合方法等進行過分類綜述（Withers el al，1998；周彥文，劉希強，2007； Nippress el al，2010；王彩霞等，2013：Ross，Ben-Zion，2014；劉翰林，吳慶舉，2017）。本文選取常用的長短窗、赤池信息準(zhǔn)則、峰態(tài)和偏態(tài)、偏振分析等單特征方法和基于整體形態(tài)的互相關(guān)方法、基于多特征的機器學(xué)習(xí)方法，分別闡述基本原理和特征函數(shù)演化，并比較分析其特點。

1.1 長短窗均值比方法

長短窗均值比方法（Short-Term-Average/Long-Term-Average，簡稱STA/LTA）基于震相到達(dá)時特征值的平均值在短窗口比長窗口變化快的特點拾取震相。STA/LTA方法的發(fā)展歷史久遠(yuǎn)，應(yīng)用廣泛，形成了反映波形不同變化特點的眾多特征函數(shù)。Stevenson（1976）首次提出長短窗方法時，使用振幅的絕對值|Y（i）|作為特征函數(shù)，其中Y（i）為第i個采樣點的振幅值。但是使用|Y（i）|或者Y（i）2作為特征函數(shù)時，長短窗比值只對振幅的變化敏感，因此Allen（1978，1982）引入了特征函數(shù)CF=Y（i）2+K（Y（i）-Y（i-1））2，其中K為根據(jù)采樣率和臺站噪聲特征設(shè)置的權(quán)重常數(shù)，該特征函數(shù)的引入使得長短窗均值比同時隨振幅和頻率變化而變化（圖1）。Bear和Kradolfer（1987）使用基于包絡(luò)函數(shù)平方的特征函數(shù)，同時引入動態(tài)閾值，既避免了瞬變噪聲引起的誤觸發(fā)，也適用于拾取地方震到遠(yuǎn)震的不同震相，同時對低信噪比微震信號也變得更敏感。為了更好地拾取各種微弱的次生震相，Earle和Shearer（1994）定義了基于Hilbert變換的包絡(luò)函數(shù)作為特征函數(shù)。高淑芳等（2008）也提出了反映振幅和頻率變化的新特征函數(shù)。余建華等（2011）引入加權(quán)系數(shù)法減少震相漏撿率。

1.2 赤池信息準(zhǔn)則方法

赤池信息準(zhǔn)則方法（Akaike Information Criteri-on，簡稱AIC）沿每個采樣點將信號分為2段，當(dāng)分割方法與2個獨立過程對應(yīng)時，對應(yīng)的AIC值最小，即全局最小值對應(yīng)震相到時。AIC方法也有非常多的變種，以下介紹幾種較為常見的AIC函數(shù)。AR-AIC（Sleeman，Van Eek，1999；Leonard，Kennett，1999）使用自回歸方法計算AIC值，即分割點對應(yīng)2個自回歸過程，自回歸階數(shù)需要反復(fù)試驗得到。Maeda（1985）提出的基于二階統(tǒng)計量的VAR-AIC方法，可以直接從地震信號計算AIC值（圖2）.0劉希強等（2009）提出TOC-AIC方法，使用三階累積量代替方差來壓制高斯有色噪聲，在低信噪比情況下有更高的精度。趙大鵬等（2012）提出峰度一赤池信息準(zhǔn)則方法（Kurtosis-AIC），使用峰度作為特征函數(shù)替代VAR一AIC中的方差，有效減少了誤報率和漏報率。

1.3 偏態(tài)/峰態(tài)方法

偏態(tài)/峰態(tài)方法（P Arrival Identification-Skewness/Kurtosis，簡稱PAI-S/K）由Saragiotis等（2002）提出，當(dāng)震相到達(dá)觀測臺站時，記錄數(shù)據(jù)的分布由具有高斯性的隨機噪聲向非高斯性轉(zhuǎn)變。峰態(tài)（Kurtosis）和偏態(tài)（Skewness）是反映數(shù)據(jù)分布形態(tài)的高級統(tǒng)計量（Higher-Order Statis-tics，HOS），數(shù)據(jù)服從高斯分布時峰態(tài)等于3，分布變寬時小于3，變窄時大于3，同理偏態(tài)在數(shù)據(jù)對稱分布時等于0，左偏時小于0，右偏時大于0（圖3）。劉勁松等（2013）引入了新的到時判定算法，改進了PAI-SIK方法的拾取精度。Baillard等（2014）使用滑動窗峰態(tài)函數(shù)和極化分析識別P、S波。峰態(tài)和偏態(tài)值方法有較好的低信噪比拾取能力和精度，屬于常用的震相拾取方法之一（Panagiotakis et al，2008；Kiiperkoch et al，2010；Nippress et al，2010；Hibert et al，2014；Langet etal，2014；Ross，2016）。

1.4 偏振分析方法

常用的2種偏振分析方法是協(xié)方差矩陣方法和奇異值分解方法，在多數(shù)情況下二者效果接近。偏振分析的思想最早由Flinn（1965）提出，他使用協(xié)方差矩陣的本征值分解來識別S波偏振，進而約束震源機制解。Montalbetti和Kanasewieh（1970）首次詳細(xì)闡述如何用對稱本征問題來進行偏振分析。分解有限長度三分量數(shù)據(jù)的協(xié)方差矩陣，得到本征值和本征向量矩陣，可以用來計算質(zhì)點運動的多種偏振特征，沿著時間軸滑動窗口，重復(fù)計算就可以得到偏振特性隨時間的變化（圖4），進而完成P或S波能量壓制、波形分解等操作，是提升S震相拾取率的常用方法（Jurkevies，1988；Robertset al，1989；Bai，Kennett，2000：Diallo et al，2006：Amoroso et al，2012；Ross et al，2014，2016）。Jackson等（1991）首次引入奇異值分解方法（Singular Value Decomposition，簡稱SVD）進行偏振分析，一步即可完成數(shù)據(jù)的分解操作。在沿時間軸滑動時，傳統(tǒng)的偏振分析方法要使用窗口內(nèi)數(shù)據(jù)重復(fù)求解協(xié)方差矩陣，Rosenberger（2010）引入迭代奇異值分解法，通過追加數(shù)據(jù)矢量來更新解，提升了奇異值分解方法的計算效率，使其適用于實時數(shù)據(jù)處理。Kurzon等（2014）基于迭代奇異值分解方法將地震波形實時分解成P分量和S分量，然后使用視出射角、奇異值、信噪比拾取器拾取P、S震相，并且嘗試多種指示變量來增強P、S波的拾取。

1.5 互相關(guān)方法

互相關(guān)系數(shù)是描述2段信號相似程度的量，地震波形的相似性主要由震源機制和經(jīng)過的路徑?jīng)Q定，因此互相關(guān)方法可以將事件波形作為模板，從連續(xù)波形中尋找相似地震（Poupinet et al，1984；Rubin，2002；Schaff，Richards，2004，2011），該方法對于拾取低信噪比震相和震顫等特殊信號非常有用（Shelly et al，2007；Peng，Gomberg，2010），能檢測到比正式目錄震級低1級的地震事件（Schaff，2010；Schaff，Waldhauser，2010；Gib-blon，Ringdal，2006）。Yang等（2009）提出聯(lián)合三分量的滑動窗互相關(guān)方法（Sliding-WindowCross-correlation，簡稱SCC），比單獨的三分量互相關(guān)方法拾取到更多低信噪比震相。

互相關(guān)方法必須先制作好模板，并且靈敏度高度依賴模板和待檢測信號的相似性，Brown等（2008）提出自相關(guān)方法來，能夠從連續(xù)波形中提取模板，但是計算十分耗時。Yoon等（2015）基于數(shù)據(jù)挖掘算法提出（Fingerprint And SimilarityThresholding，簡稱FAST），用關(guān)鍵特征代替波形，速度是自相關(guān)方法的140倍，能夠在2h內(nèi)處理單臺一周的連續(xù)波形數(shù)據(jù)。FAST方法有效地提升了模板的制作速度，但是運算量依然會隨數(shù)據(jù)長度的增加而快速增長，處理長達(dá)數(shù)月甚至數(shù)年的余震活動數(shù)據(jù)時，運算效率依然是需要解決的問題。Shelly等（2007）為了拾取更多的低頻地震，嘗試降低時間精度來延拓互相關(guān)系數(shù)的時間范圍，提高對空間位置變化的適應(yīng)性。

1.6 機器學(xué)習(xí)方法

機器學(xué)習(xí)（Machine Learning，簡稱ML）是一類從數(shù)據(jù)中尋找潛在模式，進而提升性能的龐大算法集合的總稱，大致分為經(jīng)典機器學(xué)習(xí)算法和人工神經(jīng)網(wǎng)絡(luò)算法（Artificial Neural Networks，簡稱ANN），這2種算法都包括樣本數(shù)據(jù)訓(xùn)練模型和模型預(yù)測新數(shù)據(jù)的2個階段。ANN是ML的子集，用簡單數(shù)學(xué)函數(shù)構(gòu)建神經(jīng)元網(wǎng)絡(luò)，模擬高度復(fù)雜的非線性關(guān)系，孔令文（2014）按照網(wǎng)絡(luò)結(jié)構(gòu)進行了分類。ML在地震數(shù)據(jù)處理的很多領(lǐng)域均有應(yīng)用，但是主要集中在處理分類問題。Rawles和Thurber（2015）在微震定位中使用鄰域法測定潛在震相的類別和到時。Kortstrom等（2016）使用支持向量機方法，基于多頻帶不同震相段的特征函數(shù)值，從監(jiān)測系統(tǒng)自動識別的地震中分離出爆破和誤觸發(fā)事件。Maggi等（2017）使用隨機森林方法實現(xiàn)1套自動地震事件分類器，使用多個臺站的波形和頻譜特征組合，對留尼旺火山島的地震事件進行分類。曲均浩等（2012）使用反向傳播神經(jīng)網(wǎng)絡(luò)分類近震和遠(yuǎn)震。Esposito等（2013）使用多層感知機分類意大利南部斯特龍博利火山的滑坡和其他事件。Horstmann等（2013）使用自組織映射神經(jīng)網(wǎng)絡(luò)檢測和分類地震。Giudicepietro等（2017）使用多層感知機分類地方震、區(qū)域震和遠(yuǎn)震。

震相自動拾取主要采用反向傳播神經(jīng)網(wǎng)絡(luò)（Back Propagation Neural Network，簡稱BPNN），Wang和Teng（1995，1997）設(shè)計使用長短窗序列和滑動窗頻率譜作為輸入的神經(jīng)網(wǎng)絡(luò)算法，與STA/LTA相比，在拾取精度和低信噪比適應(yīng)性上取得一定優(yōu)勢。Dai和Mocbeth（1995）對ANN在震相到時拾取方面的性能進行了測試。張范民等（1998）使用三分量地震數(shù)據(jù)的矢量模作為ANN輸入，確定震相類型和拾取到時。Zhao和Takano（1999）使用BPNN從IRIS數(shù)據(jù)中拾取初至，用于地幔結(jié)構(gòu)層析成像研究。Gentili和Michelini（2006）用具有高泛化能力的樹形神經(jīng)網(wǎng)絡(luò)IU-ANT2拾取P和S震相。

1.7 總結(jié)與分析

STA/LTA、AIC、PAI-S/K、偏振分析屬于單特征方法，所使用的特征函數(shù)反映波形記錄的某一項或幾項特征，有相對明確的物理含義，都發(fā)展出了大量的衍生變種和組合算法。其中STA/LTA方法的突出優(yōu)點是計算簡單快速，不需要先驗信息，只使用震相監(jiān)測點之后少量的采樣點（圖1a），因此適用于對實時性和處理效率要求較高的地震監(jiān)測、地震預(yù)警和大數(shù)據(jù)量處理。STA/LTA方法的缺點之一是低信噪比時，到時拾取精度會受到影響，因此通常配合其他方法進行精確拾取。而Lomax等（2012）基于多頻帶濾波和Bear和Kradolfer（1987）的方法提出FilterPicker方法，較好地解決了余震震相受其他地震尾波干擾的問題。STA/LTA方法的另一個缺點是需要根據(jù)不同應(yīng)用調(diào)整參數(shù)和選擇特征函數(shù)，Vassallo等（2012）對自動拾取方法的優(yōu)化策略進行了研究，劉晗和張建中（2014）對各種特征函數(shù)和改進算法進行了分析。

AIC方法的突出優(yōu)點是到時拾取精度高。其計算效率相對較低，并且假設(shè)背景噪聲和震相信號是2個獨立的過程，例如，當(dāng)數(shù)據(jù)窗口內(nèi)的震相數(shù)量大于一時便不符合這一假設(shè)，拾取結(jié)果就會受到影響。因此，AIC方祛的窗口選擇十分重要，通常使用STA/LTA提供的先驗信息進行精確拾?。R強等，2013；Bastin et al，2013；牟磊育，吳建平，2015）

PAI-S/K方法與STA/LTA一樣，是常用的初至拾取方法，使用統(tǒng)計信息代替信噪比，在背景噪聲接近高斯分布時有很好的拾取精度，計算簡單高效，低信噪比環(huán)境下表現(xiàn)更穩(wěn)定。缺點是在背景噪聲不符合高斯分布時，拾取精度會受到影響。

偏振分析方法充分利用三分量信號的運動特征，可以得到多種偏振特征。震相的類型是前述方法無法或者很難直接獲取的。由于與P波尾波重疊，出現(xiàn)轉(zhuǎn)換震相，或者可能出現(xiàn)剪切波分裂，拾取S震相變得很困難。通過偏振分析方法能夠完成震相類型識別、P波能量壓制，進而提高S震相拾取率，這對于識別余震很有幫助。

余震與主震在構(gòu)造上有密切聯(lián)系，發(fā)震機理相似，空間上高度密集，大多數(shù)情況下波形相似度非常高，前述單特征方法都沒有利用這些特點。由于噪聲分布有隨機性，互相關(guān)操作時會被有效的壓制。互相關(guān)方法比對事件波形，是目前低信噪比拾取能力最強，拾取精度最高的方法。當(dāng)然，事件波形形態(tài)受震源機制和傳播路徑等與具體事件相關(guān)的因素影響，這也導(dǎo)致波形互相關(guān)方法缺乏通用性。

機器學(xué)習(xí)算法基于多特征值和概率預(yù)測，最終的效果主要由具體算法和參數(shù)、特征選擇和訓(xùn)練數(shù)據(jù)決定，震相拾取最常用的是經(jīng)典機器學(xué)習(xí)算法和BPNN。經(jīng)典機器學(xué)習(xí)算法需要人工提取和選擇特征值。BPNN實際上比SVM、隨機森林等經(jīng)典機器學(xué)習(xí)方法更早出現(xiàn)，理論上通過訓(xùn)練可以自行提取特征值，但是層數(shù)增加時需要訓(xùn)練的權(quán)重系數(shù)很多，容易產(chǎn)生過擬合，受限于當(dāng)時的硬件條件通常只采用三層網(wǎng)絡(luò)，依然依賴人為選擇特征，并且反復(fù)調(diào)參試驗，實際使用需要一定的理論基礎(chǔ)和經(jīng)驗。經(jīng)典機器學(xué)習(xí)和全連接神經(jīng)網(wǎng)絡(luò)取得了優(yōu)于傳統(tǒng)算法的低信噪比適應(yīng)性和拾取精度，通用性也要優(yōu)于互相關(guān)方法，但是在實際監(jiān)測工作中并沒有得到廣泛應(yīng)用。

余震震相識別最大的難度來源j二震相信號受到主震和其他余震尾波的強烈干擾，從單臺拾取震相到時的角度，F(xiàn)ilterPicker和互相關(guān)方法的適用性和穩(wěn)定性最好。但是隨著偏移疊加技術(shù)的興起，傳統(tǒng)單特征方法所發(fā)展出的眾多特征函數(shù)也可以用于低信噪比環(huán)境下的余震識別，本文將在第三章中進行介紹。

2 走時定位法研究現(xiàn)狀

基于走時的常規(guī)方法被地震監(jiān)測和預(yù)警所普遍采用，可以分為2部分：從連續(xù)觀測波形中拾取震相；使用震相信息確定震源位置和發(fā)震時刻（圖5）。一次余震有唯一的初始破裂位置和發(fā)震時刻，如果已知臺站的空間位置，使用合適的理論模型可以預(yù)估震相到時。因此，利用多個臺站拾取到的震相，通過走時定位方法，能夠?qū)ふ业胶线m的震源位置和發(fā)震時刻，使得走時殘差小于設(shè)定閾值，即可識別出余震。

走時定位方法通過最小化觀測與理論走時的殘差確定震源位置和發(fā)震時刻，其發(fā)展已經(jīng)較為成熟，目前的難點在于深度定位（鄭勇，謝祖軍，2017），

常用的走時定位方法依據(jù)求解方式可分為傳統(tǒng)幾何定位方法、線性定位法和非線性方法，依據(jù)殘差計算方式可分為絕對定位方法和相對定位方法。

幾何定位法使用均勻介質(zhì)走時方程和幾何知識，在平面上作圖確定震源點，是早期地震定位中的主要方法，常見的有方位角法、石川法、交切法。雖然精度有限，幾何方法在余震識別中仍然是有效的約束手段。

線性定位方法最早由Geiger（1912）提出，其基本原理是將關(guān)于震源時空參數(shù)的非線性方程組作泰勒展開，簡化為線性方程組，用最小二乘法求解，通過取初值迭代得到震源參數(shù)。由于理論速度模型與真實速度模型存在偏差，在Geiger方法基礎(chǔ)上衍生出眾多改進型方法。Douglas（1967）提出的JED方法和Dewey（1972）提出的JHD方法，認(rèn)為殘差由臺站下方速度結(jié)構(gòu)的差異引起，引入臺站校正系數(shù)校正；Crosson（1976）提出的SSH方法，認(rèn)為走時殘差由速度結(jié)構(gòu)偏差引起，同時進行地震定位和地下速度結(jié)構(gòu)反演。

當(dāng)2次地震發(fā)生在同一個震源區(qū)，并且震源區(qū)的范圍遠(yuǎn)小于到某個臺站的距離，相同震相到達(dá)這個臺站所經(jīng)過的路徑絕大部分是相同的，也就是說速度模型引起的走時殘差近似相等，因此可以假設(shè)2個震相間的走時差是由震源間的距離引起的，利用該思想的定位方法稱之為相對定位方法。主事件方法以一個定位比較準(zhǔn)確的地震作為參考事件，對周圍的其他地震進行相對定位（Evernden，1969a，b； Spence，1980），由于被定位事件必須距離參考事件足夠近，主事件方法的定位范圍受到一定限制。Waldhauser（2000）提出雙差定位方法（DD），將震源區(qū)距離較近的地震構(gòu)成地震對，利用事件對的走時差約束相對位置，由于不需要主事件，適用范圍比主事件更廣。

地震定位是一個非線性反演問題，可以使用非線性最優(yōu)化方法求解，國內(nèi)也有很多研究，常見的有Powell方法（唐國興，1979）、遺傳算法（萬永革，李鴻吉，1995；周民都等，1999；王文萍，王慶良，1999）、模擬退火法（高星等，2002；孫麗，2016）、網(wǎng)格搜索法（趙仲和，牟磊育，2005；宋維琪，楊曉東，2011）、單純形法（趙珠等，1994）。非線性方法的優(yōu)勢是可以較好地避免結(jié)果陷人局部極小值，但是運算量通常要高于線性方法。

研究發(fā)震構(gòu)造的空間展布是余震自動識別的重要應(yīng)用，對余震的定位精度提出較高的要求，由于絕對定位會引入理論模型和真實地下結(jié)構(gòu)在整個射線路徑上的偏差，因此無法滿足精度需求。余震序列在空間上高度密集，這一特征是采用主事件和DD等相對定位方法的有利條件。DD方法目前是震后重定位余震，確定發(fā)震構(gòu)造最常用的方法（黃媛等，2008；房立華等，2011；蘇金蓉等，2013：張廣偉等，2014）。

3 偏移疊加法研究現(xiàn)狀

基于走時定位識別余震只使用震相到時和類型，識別能力受制于單臺震相拾取算法，如果震相引起的特征變化無法達(dá)到觸發(fā)閾值，顯然無法通過聯(lián)合走時來定位和識別余震。一次余震會引發(fā)多個臺站的波形記錄變化，在不同臺站記錄中獨立進行震相的識別，沒有利用這種多臺站震相信息間的一致性。偏移疊加方法不需要拾取到時或者識別震相類型，直接對齊各個臺站的波形或者特征值，加權(quán)平均，能有效壓制噪聲，提高地震識別準(zhǔn)確率和定位精度。由于石油、頁巖氣、水庫監(jiān)測等領(lǐng)域?qū)φT發(fā)微震監(jiān)測的投入，偏移疊加方法發(fā)展十分迅速，可以將其分為2類：一類是基于波動方程的波場反向傳播聚焦方法，這類方法完全通過數(shù)值計算，因而計算量很大，也容易受強噪聲和地下結(jié)構(gòu)不均勻的影響（Mcmechan et al，1982，1985；Gajewski，Tessmer，2005）；另一類是直接疊加特征函數(shù)的震源掃描方法（Source Scanning Algorithm，簡稱SSA）及其衍生方法。

SSA由Kao和Shan（2004，2007）提出，在大致確定震源參數(shù)的情況下，搜索潛在震源位置和發(fā)震時間，疊加各個臺站歸一化的記錄振幅絕對值，得到的亮度函數(shù)最大值點即認(rèn)為是最佳解。Baker等（2005）用P波包絡(luò)進行疊加，并引入基爾霍夫重建法（Kirchhoff reconstruction method）。Gharti等（2010）根據(jù)理論預(yù)測到時，截取數(shù)據(jù)并旋轉(zhuǎn)至射線路徑，最后同時使用P波和S波包絡(luò)進行疊加。Grigoli等（2013，2014，2016）在SSA方法基礎(chǔ)的上使用STA/LTA作為特征函數(shù)進行疊加，并引入主事件方法和震源臺站修正項，進一步提高定位精度。Langet等（2014）則提出基于連續(xù)峰度函數(shù)的偏移疊加識別微震。

模板匹配技術(shù)（Matched Filter Technique，簡稱MFT）結(jié)合波形互相關(guān)和臺陣疊加技術(shù)，將多通道的互相關(guān)波形偏移到模板事件的發(fā)震時刻，疊加取平均值得到平均互相關(guān)波形，最大值超過設(shè)定閾值就認(rèn)為檢測到一次相似地震（Gibbons，Ringdal，2006；Shelly et al，2007；Peng，Zhao，2009），模板的位置視為被匹配地震的位置，模板發(fā)震時刻偏移后確定發(fā)震時刻。MFT被廣泛運用于完備目錄中遺漏的早期余震（Peng，Zhao，2009；Lenglinee etal，2012：Meng et al，2013：Rastin，2013：Kato，2012，2014a，2014b； Wu et al，2014）和研究遠(yuǎn)場動態(tài)觸發(fā)地震（Yukutake，2013； Meng，Peng，2014）0在MFT和SSA基礎(chǔ)上，Zhang和Wen（2015a）提出匹配定位方法（Match And Locate，簡稱M&L;），通過搜索模板事件附近的三維空間網(wǎng)格，將各個臺站的互相關(guān)波形對齊疊加，互相關(guān)系數(shù)最高的點即為被檢測事件的空間位置。M&L;方法不僅進一步提高對低信噪比事件的檢測能力，同時也給出精確的空間位置，被用于精確定位朝鮮核試驗（Zhang，Wen，2015b）和監(jiān)測火山噴發(fā)前的微震活動（Zhang，Wen，2015c）。

偏移疊加方法直接疊加特征函數(shù)，比傳統(tǒng)走時定位方法有更高的地震識別能力，其函數(shù)既可以選擇傳統(tǒng)震相拾取方法所用特征函數(shù)，也可以是互相關(guān)波形，不同偏移疊加方法的識別能力和特點與所選疊加函數(shù)有關(guān)。余震震源機制和震中位置重復(fù)性好，這對基于互相關(guān)方法的MFT和M&L;十分有利，MFT是目前檢測能力最強的余震自動識別方法之一，可以識別出10倍于正式目錄的地震，常被應(yīng)用于完備余震目錄。M&L;采用相對定位方式，非常適合余震的高精度定位，是已知事件大致位置和時間的情況下，目前識別能力和定位精度最高的方法。但是二者對模板有非常高的依賴性，需要大量的先驗信息，通用性較差，在識別余震時需要準(zhǔn)備大量的模板。因此，選用其他單特征震相拾取方法所用的特征函數(shù)，更利于提升性能和通用性，但是檢測能力會有所降低。余震在空間上密集分布，能將搜索區(qū)域限定在相對小的范圍內(nèi)，這對于計算量龐大的偏移疊加方法也是十分有利的條件。

4 討論與展望

如前所述，余震自動識別仍然有許多技術(shù)難點亟待解決，震相自動識別技術(shù)是余震自動識別技術(shù)的基礎(chǔ)，也是最大的難點所在，多臺震相聯(lián)合方法是增強余震識別能力的有效手段。本文以傳統(tǒng)震相拾取方法和互相關(guān)方法為出發(fā)點，首先總結(jié)現(xiàn)有震相識別技術(shù)及其特征函數(shù)，然后對基于走時定位的識別方法進行綜述，同時總結(jié)了近些年興起的偏移疊加方法的發(fā)展，并分別分析了各項技術(shù)在余震自動識別中的優(yōu)缺點和適用性。下面將從偏移疊加和震相識別方法2方面，分析存在的問題和原因，并簡述對余震自動識別技術(shù)未來工作的設(shè)想。

4.1 優(yōu)化偏移疊加方法

偏移疊加方法通過疊加特征函數(shù)，不需要拾取震相，能識別單臺震相信號無法達(dá)到觸發(fā)閾值的余震，識別能力明顯強于基于走時定位的常規(guī)方法，各種衍生方法都可以用于余震的識別。常規(guī)走時定位方法的優(yōu)點是計算簡單快速，容易并行執(zhí)行，不需要任何先驗信息，適用于日常監(jiān)測和地震預(yù)警等對時效要求比較高的應(yīng)用。偏移疊加方法有非常高的識別能力，對潛在震源時間和空間進行搜索需要耗費大量的計算資源，雖然余震和人工誘發(fā)地震類似，震源的密集空間分布縮小了搜索范圍，但是時效上仍然需要提升。因此，從實際應(yīng)用的角度出發(fā)，偏移疊加方法的性能優(yōu)化將成為未來余震識別的重點發(fā)展方向。

4.2 探索深度學(xué)習(xí)在余震識別領(lǐng)域的應(yīng)用

波形特征函數(shù)是震相到時拾取和偏移疊加方法的基礎(chǔ)，因此單臺震相拾取方法仍然是提升余震自動識別能力的重點。特征值是對觀測數(shù)據(jù)的抽象，傳統(tǒng)方法選取人為設(shè)計的數(shù)種特征值，通常不需要先驗信息，具有計算簡單快速的特點，但是單一的特征值無法描述震相信號的所有變化，參數(shù)設(shè)置也很難適應(yīng)復(fù)雜多變的余震波形。互相關(guān)方法利用波形的整體形態(tài)特征，壓制噪聲的能力非常強，震相拾取能力遠(yuǎn)超傳統(tǒng)方法。波形形態(tài)中包含與具體事件相關(guān)的信息，如震源機制和空間位置，余震的震源機制和空間位置雖然多數(shù)非常接近，但也存在偏差較大的情況。因此，互相關(guān)方法通常需要準(zhǔn)備大量模板，如Peng和Zhao（2009）在識別余震時使用3647個地震事件的波形作為模板，這導(dǎo)致其通用性和時效性較差，自相關(guān)方法和FAST方法的出現(xiàn)也沒有從根本上解決這一問題。

筆者認(rèn)為介于傳統(tǒng)方法和互相關(guān)方法之間的綜合特征方法是未來發(fā)展的方向，即從波形數(shù)據(jù)中盡可能多地抽象出與具體余震事件無關(guān)的特征，綜合這些特征進行震相拾取，應(yīng)當(dāng)注意的是，這里的綜合特征并不等同于多特征。多特征方法中最典型的屬經(jīng)典機器學(xué)習(xí)方法，這些方法依賴于人工從數(shù)據(jù)中提取出的數(shù)種特征值。BPNN等全連接神經(jīng)網(wǎng)絡(luò)方法理論上可以自動提取特征值，但是層數(shù)增加時需要訓(xùn)練的權(quán)重過多，實際應(yīng)用中以3層網(wǎng)絡(luò)為主，一般通過增加神經(jīng)元數(shù)量來提升性能，與經(jīng)典機器學(xué)習(xí)方法一樣，比較依賴人為選擇。

深度學(xué)習(xí)方法是近年興起的ANN子集，優(yōu)點是擁有更多隱含層，通過多層次的抽象比簡單的增加神經(jīng)元數(shù)量具有更強的描述能力，可以從原始數(shù)據(jù)中抽象出所需特征，這些特征也許我們從未注意到過，人工智能現(xiàn)今的優(yōu)異成績幾乎都是由深度學(xué)習(xí)方法所取得。Perol等（2017）使用CNN方法識別和分類美國中部的誘發(fā)地震，在取得比匹配濾波方法更高識別能力的同時，計算耗時和存儲使用量比自相關(guān)方法和FAST方法明顯降低。卷積神經(jīng)元網(wǎng)絡(luò)（Convolutional Neural Network，簡稱CNN）是深度神經(jīng)網(wǎng)絡(luò)的一種，使用卷積核作為中介，識別數(shù)據(jù)中存在的固有局部結(jié)構(gòu)，降低權(quán)重參數(shù)數(shù)量，可以有效避免過擬合，十分適合處理圖像和聲學(xué)數(shù)據(jù)。地震連續(xù)波形數(shù)據(jù)與聲學(xué)數(shù)據(jù)有很多共同點，都是記錄介質(zhì)的震動信息，并且存在固有的局部結(jié)構(gòu)，卷積神經(jīng)網(wǎng)絡(luò)方法識別震相可能會取得與處理聲學(xué)數(shù)據(jù)相近的效果。

深度學(xué)習(xí)方法的運算量集中在模型訓(xùn)練階段，使用模型進行預(yù)測的速度則非?？欤Ｐ偷姆夯梢蕴嵘m應(yīng)性，因此具備同時兼顧運行效率和通用性的能力，在余震識別領(lǐng)域尚處于起步階段。深度學(xué)習(xí)方法，屬于黑盒算法，即無法從其訓(xùn)練模型中分析震相特征的抽象過程，以及如何組合特征識別震相，因此無法提升對地震信號的認(rèn)識，但是優(yōu)異的性能對于解決震相識別這類實際問題很有意義，未來將成為震相識別，乃至地震監(jiān)測技術(shù)的重要發(fā)展方向。

如果深度學(xué)習(xí)能夠準(zhǔn)確拾取震相到時，那么在識別余震時，常規(guī)走時定位法所面臨的低信噪比震相無法拾取和精確性不足問題就迎刃而解，而偏移疊加方法也能夠更好地約束搜索范圍，提高定位精度。當(dāng)然，立足更長遠(yuǎn)的未來，深度學(xué)習(xí)技術(shù)也有可能替代偏移疊加，以更復(fù)雜和更抽象的方式實現(xiàn)多臺震相聯(lián)合，甚至進一步由余震識別這個小的領(lǐng)域向地震監(jiān)測和地震數(shù)據(jù)處理發(fā)展。

本人關(guān)于地震定位、震相拾取的理論學(xué)習(xí)和算法實現(xiàn)都是在中國地震臺網(wǎng)中心馬延路博士的指導(dǎo)下完成，審稿專家對本文也提出了寶貴意見，在此一并感謝。

參考文獻(xiàn)：

房立華，吳建平，張?zhí)熘?，?2011.2011年云南盈江M55.8地震及其余震序列重定位[J].地震學(xué)報，33（2）：262-267.

高淑芳，李山有，武東坡，等.2008.一種改進的STA/LTA震相自動識別方法[J].世界地震工程，24（2）：37-41.

高星，王衛(wèi)民，姚振興.2002.用于地震定位的SAMS方法[J].地球物理學(xué)報，45（4）：525-532.

黃媛，吳建平，張?zhí)熘?，等?008.汶川8.0級大地震及其余震序列重定位研究[J].中國科學(xué)：地球科學(xué)，38（10）：1242-1249.

孔令文.2014.神經(jīng)網(wǎng)絡(luò)的分類及其應(yīng)用[J].電子技術(shù)與軟件工程，（17）：25.

劉晗，張建中.2014.微震信號自動檢測的STA/LTA算法及其改進分析[J].地球物理學(xué)進展，29（4）：1708-1714.

劉翰林，吳慶舉.2017.地震自動識別及震相自動拾取方法研究進展[J].地球物理學(xué)進展，32（3）：1000-1007.

劉勁松，王赟，姚振興.2013.微地震信號到時自動拾取方法[J].地球物理學(xué)報，56（5）：1660-1666.

劉希強，周彥文，曲均浩，等.2009.應(yīng)用單臺垂向記錄進行區(qū)域地震事件實時檢測和直達(dá)P波初動自動識別[J].地震學(xué)報，31（3）：260-271.

馬強，金星，李山有，等.2013.用于地震預(yù)警的P波震相到時自動拾取[J].地球物理學(xué)報，56（7）；2313-2321.

牟磊育，吳建平.2015.利用單臺觸發(fā)檢測臺陣數(shù)據(jù)庫事件[J].地球物理學(xué)進展，30（6）；2511-2516.

曲均浩，劉希強，吳丹彤，等.2012.基于神經(jīng)網(wǎng)絡(luò)的近震與遠(yuǎn)震識別[J].地震研究，30（3）：360-366.

宋金，蔣海昆.2009.序列衰減與余震激發(fā)研究進展及應(yīng)用成果[J].地震地質(zhì)，31（3）：559-571.

宋維琪，楊曉東.2011.解域約束下的微地震事件網(wǎng)格搜索法、遺傳算法聯(lián)合反演[J].石油地球物理勘探，46（2）：259-266，160.

蘇金蓉，鄭鈺，楊建思，等.2013.2013年4月20日四川蘆山7.0級地震與余震精確定位及發(fā)震構(gòu)造初探[J].地球物理學(xué)報，56（8）：2636-2644.

孫麗.2016.國家地震臺網(wǎng)地震定位方法的改進[J].國際地震動態(tài)，（11）：12-16.

唐國興.1979.用計算機確定地震參數(shù)的一個通用方法[dl.地震學(xué)報，1（2）：186-196.

萬永革，李鴻吉.1995.遺傳算法在確定震源位置中的應(yīng)用[J].地震地磁觀測與研究，16（6）：1-7.

王彩霞，白超英，王馨.2013.地震震相初至自動檢測技術(shù)綜述[J].地球物理學(xué)進展，28（5）：2363-2375.

王文萍，王慶良.1999.利用遺傳算法和最小二乘聯(lián)合反演共和地震位錯參數(shù)[J].地震學(xué)報，21（3）；62-67.

徐錫偉，陳桂華，于貴華，等.2013a.蘆山地震發(fā)震構(gòu)造及其與汶川地震關(guān)系詞論[J].地學(xué)前緣，20（3）：11-20.

徐錫偉，聞學(xué)澤，韓竹軍，等.2013b.四川蘆山7.0級強震：一次典型的盲逆斷層型地震[J].科學(xué)通報，58（20）：1887-1893.

徐錫偉，聞學(xué)澤，葉建青，等.2008.汶川MS8.0地震地表破裂帶及其發(fā)震構(gòu)造[J].地震地質(zhì)，30（3）；597-629.

余建華，李丹丹，韓國棟.2011.特征函數(shù)響應(yīng)特性分析及STA/LTA方法的改進[J].國外電子測量技術(shù)，30（7）：17-19，23，27.

余娜，張曉清，袁伏全.2016.2016年青海門源MS6.4地震余震序列參數(shù)穩(wěn)定性分析[J].地震研究，39（S1）；62-68，133.

張范民，李清河，張元生，等.1998.利用人工神經(jīng)網(wǎng)絡(luò)理論對地震信號及地震震相進行識別[J].西北地震學(xué)報，20（4）；43-49.

張廣偉，雷建設(shè)，梁姍姍，等.2014.2014年8月3日云南魯?shù)镸S6.5級地震序列重定位與震源機制研究[J].地球物理學(xué)報，57（9）：3018-3027.

張?zhí)K平，孫艷萍，陳文凱.2013.余震信息在岷縣漳縣6.6級地震震后應(yīng)急救援重災(zāi)區(qū)快速判定時作用探討[J].地震工程學(xué)報，35（3）：465-470.

趙大鵬，劉希強，李紅，等.2012.峰度和AIC方法在區(qū)域地震事件和直達(dá)P波初動自動識別方面的應(yīng)用[J].地震研究，35（2）；220-225，295.

趙仲和，牟磊育.2005.單臺地震自動定位網(wǎng)格搜索法及其MATLAB試驗[J]地震地磁觀測與研究，26（4）：1-12.

趙珠，丁志峰，易桂喜，等.1994.西藏地震定位——一種使用單純形優(yōu)化的非線性方法[J].地震學(xué)報，16（2）；212-219.

鄭勇，謝祖軍.2017.地震震源深度定位研究的現(xiàn)狀與展望[J].地震研究，40（2）：167-175，333.

鄭韻，姜立新，楊天青，等.2015：利用余震頻度分布進行宏觀震中快速判定[J].地震，35（2）：121-132.

周民都，張元生，張樹勛.1999.遺傳算法在地震定位中的應(yīng)用[J].西北地震學(xué)報，21（2）：167-171.

周彥文，劉希強.2007.初至震相自動識別方法研究與發(fā)展趨勢[J].華北地震科學(xué)，25（4）：18-22.

莊建倉，馬麗.2000.主震和余震——從大森公式到ETAS模型[J].國際地震動態(tài)，（5）：12-18.

Allen R V.1978.Automatic earthquake recognition and timing from singletraces[J].Bulletin of the Seismological Society of America，68（5）：1521-1532.

Allen R.1982.Automatic phase pickers；Their present use and futureprospects[J].Bulletin of the Seismological Society of America，72（6）：225-242.

Amoroso O，Maercklin N，Zollo A.2012.S-wave identification by polari-zation filtering and waveform coherence analyses[J].Bulletin of theSeismological Society of America，102（2）：854-861.

Baer M，Kradolfer U.1987.An automatic phase picker for local andteleseismic events[J].Bulletin of the Seismological Society of Amer-ica，77（4）：1437-1445.

Bai C Y，Kennett B L N.2000.Automatic phase-detection and identifi-cation by full use of a single three-component broadband seismo-gram[J].Bulletin of the Seismological Society of America，90（1）：187-198.

Baillard C，Crawford W C，Ballu V，et al.2014.An automatic kurtosis-based P-and S-phase picker designed for local seismic networks[J].Bulletin of the Seismological Society of America，104（1）：394-409.

Baker T，Granat R，Clayton R W.2005.Real-time earthquake locationusing Kirchhoff reconstruction[J].Bulletin of the Seismological So-ciety of America，95（2）：699-707.

Brown J R，Beroza G C，Shelly D R.2008.An autocorrelation method todetect low frequency earthquakes within tremor[J].Geophysical Re-search Letters，35（16）：1-5.

Crosson R S.1976.Crustal structure modeling of earthquake data，1，si-multameous least squares estimation of hypocenter and velocity pa-rameters[J].Geophys Res，81（17）：3036-3046.

Dai H，Macbeth C.1995.Automatic picking of seismic arrivals in localearthquake data using an artificial neural network[J].GeophysicalJournal International，120（3）：758-774.

Dewey James W.1972.Seismicity and tectonics of western venezuela[J].Bulletin of the Seismological Society of America，62（6）：1711-1751.

Diallo M S，Kulesh M，Holschneider M，et al.2006.Instantaneous polari-zation attributes based on an adaptive approximate covariance method[J].Geophysics，71（5）：V99.

Douglas A.1967.Joint epicentre determination[J].Nature，215（5096）：47-48.

Earle P S，Shearer P M.1994.Characterization of global seismograms u-sing an automatic-picking algorithm[J].Bulletin of the Seismologi-cal Society of America，84（2）：366-376.

Esposito A M，D Auria L，Giudicepietro F，et al.2013.Automatic recogni-tion of landslides based on neural network analysis of seismic sig-nals：an application to the monitoring of stromboli volcano（southernItaly）[J].Pure and Applied Geophysics，170（11）：1821-1832.

Evernden J F.1969a.Identification of earthquakes and explosions by useof teleseismic data[J].Journal of Geophysical Research，74（15）：3828-3856.

Evernden J F.1969b.Precision of epicenters obtained by small numbersof world-wide stations[J].Bulletin of the Seismological Society ofAmerica，59（3）：1365-1398.

Flinn E A.1965.Signal analysis using rectilinearity and direction of parti-cle motion[J].Proceedings of the IEEE，53（12）：1874-1876.

Gajewski D，Tessmer E.2005.Reverse modelling for seismic event charac-terization[J].Geophysical Journal International，163（1）：276-284.

Geiger L.1912.Probability method for the determination of earthquakeepicenters from arrival time only[J].Bull St Louis Univ，8；60-71.

Gentili S，Michelini A.2006.Automatic picking of P and S phases using aneural tree[J].Journal of Seismology，10（1）：39-63.

Gharti H N，Oye V，roth M，Kiihn D.2010.Automated microearthquakelocation using envelope stacking and robust global optimization[J].Geophysics，75（4）：27-46.

Gibbons S J，Ringdal F.2006.The detection of low magnitude seismic e-vents using array-based waveform correlation[J].GeophysicalJournal International，165（1）：149-166.

Giudicepietro F，Esposito A M，Ricciolino P.2017.Fast discrimination oflocal earthquakes using a neural approach[J].Seismological Re-search Letters，88（4）：1089-1096.

Grigoli F，Cesca S，Amoroso O，et al.2014.Automated seismic event loca-tion by waveform coherence analysis[J].Geophysical Journal Inter-national，196（3）：1742-1753.

Grigoli f，Cesca S，Krieger L，et al.2016.Automated microseismic eventlocation using Master-Event Waveform Stacking[J].Scientific Re-ports，6：1-13.

Grigoli F，Cesca S，Vassallo M，et al.2013.Automated Seismic Event Lo-cation by Travel-Time Stacking；An Application to Mining InducedSeismicity[J].Seismological Research Letters，84（4）：666-677.

Hibert C，Mangeney A，Grandjean G，et al.2014.Automated identifica-tion，location，and volume estimation of rockfalls at Piton de la Four-naise volcano[J].Journal of Geophysical Research：Earth Surface，119（5）：1082-1105.

Horstmann T，Harrington R M，Cochran E S.2013.Semiautomated tremordetection using a combined cross-correlation and neural networkapproach[J].Journal of Geophysical Research；Solid Earth，118（9）：4827-4846.

Jackson G M，Mason I M，Greenhalgh S A.1991.Principal componenttransforms of triaxial recordings by singular value decomposition[J].GEOPHYSICS，56（4）：528-533.

Jurkevics A.1988.Polarization analysis of three-component array data[J].Bulletin of the Seismological Society of America，78（5）：1725-1743.

Kao H，Shan S J.2004.The Source-Scanning Algorithm；Mapping thedistribution of seismic sources in time and space[J].GeophysicalJournal International，157（2）：589-594.

Kao H，Shan S J.2007.Rapid identification of earthquake rupture planeusing Source-Scanning Algorithm[J].Geophysical Journal Interna-tional，168（3）：1011-1020.

Kato A，Nakagawa S.2014a Multiple slow-slip events during a fore-shock sequence of the 2014 Iquique，Chile Mw 8.1 earthquake[J].Geophysical Research Letters，41（15）：5420-5427.

Kato A，Obara K，Igarashi T，et al.2012.Propagation of slow slip leadingup to the 2011 MW9.0 Tohoku-Oki earthquake[J].Science，335（6069）：705-708.

Kato A，Obara K.2014b.Step-like migration of early aftershocks follow-ing the 2007 Mw 6.7 Noto-Hanto earthquake，Japan[J].Geophysi-cal Research Letters，41（11）：3864-3869.

Kortstrom J，Uski M，Tiira T.2016.Automatic classification of seismic e-vents within a regional seismograph network[J].Computers and Ge-osciences，87：22-30.

Kiiperkoch L，Meier T，Lee]，et al.2010.Automated determination of P-phase arrival times at regional and local distances using higher orderstatistics[J].Geophysical Journal International，181：1159-1170.

Kurzon I，Vernon F L，Rosenberger A，et al.2014.Real-time automaticdetectors of P and S waves using singular value decomposition[J].Bulletin of the Seismological Society of America，104（4）：1696-1708.

Langet N，Maggi A，Michelini A，et al.2014.Continuous kurtosis-basedmigration for seismic event detection and location，with application toPiton de la Fournaise volcano，La Reunion[J].Bulletin of the Seis-mological Society of America，104（1）：229-246.

Lenglinee O，Enescu B，Peng Z，et al.2012.Decay and expansion of theearly aftershock activity following the 2011，MW9.0 Tohoku earth-quake[J].Geophysical Research Letters，39（18）：L18309.

Leonard M，Kennett B L N.1999.Multi-component autoregressive tech-niques for the analysis of seismograms[J].Physics of the Earth andPlanetary Interiors，113（1-4）：247-263.

Lomax A，Satriano C，Vassallo M.2012.Automatic picker developmentsand optimization：FilterPicker-a robust，broadband picker for real-time seismic monitoring and earthquake early warning[J].Seismo-logical Research Letters，83（3）：531-540.

Maeda N.1985.A method for reading and checking phase time in auto-processing system of seismic wave data[J].Zisin，38（3）：365-379.

Maggi A，F(xiàn)errazzini V，Hibert C，et al.2017.Implementation of a multista-tion approach for automated event classification at piton de la Four-naise volcano[J].Seismological Research Letters，88（3）：878-891.

Mcmechan G A，Luetgert J H，Mooney W D.1985.Imaging of earthquakesources in Long Valley Caldera，alifornia，1983 [J].Bulletin of theSeismological Society of America，75（4）：1005-1020.

Mcmechan G A.1982.Determination of source parameters by wavefieldextrapolation[J].Geophysical Journal of the Royal Astronomical So-ciety，71（3）：613-628.