孫印 潘素珍 劉明軍

中國地震局地球物理勘探中心，鄭州市文化路75號 450002

0 引言

地震事件的識別與拾取是地震數(shù)據(jù)處理過程中非常重要的一步，拾取方法的優(yōu)劣關(guān)系著事件拾取的精度、拾取的效率，好的拾取方法不僅可以提高地震學(xué)研究效率，而且還能促進地震速報、地震預(yù)警等防震減災(zāi)工作的進步。起初，地震工作者在進行地震數(shù)據(jù)處理時，主要通過人工肉眼分析震相來識別地震事件，這種方法效率低下，并且容易出現(xiàn)人為性疏漏。隨著地震臺網(wǎng)越來越密集，數(shù)據(jù)量越來越龐大，人工拾取地震事件的質(zhì)量和效率已遠遠不能滿足地震數(shù)據(jù)處理的需求。另外，隨著地震學(xué)的發(fā)展，其研究對象已不僅局限于較大震級地震，對于微小震級地震事件的研究也不斷深入，這也對能有效識別微小地震的地震事件拾取技術(shù)提出了更高的要求，即需要效率更高、精度更高、錯誤更少的自動識別地震事件技術(shù)。圖1為Tamaribuchi等（2016）在研究熊本地震過程中自動拾取震相與人工拾取震相的對比圖。由圖1可見，自動拾取事件數(shù)量與人工拾取數(shù)量之比高達14，并且對于微地震有更強大的識別能力，而對于這些微震，人工拾取往往難以分辨，Tamaribuchi等（2016）的研究表明，發(fā)展自動震相拾取技術(shù)是地震學(xué)發(fā)展的必經(jīng)之路。

圖1 自動拾取震相與人工拾取震相對比（據(jù)Tamaribuchi等（2016））

20世紀70年代以來，自動識別地震事件和自動拾取震相的方法越來越豐富，也越來越成熟。無論是對于流動臺網(wǎng)還是固定臺網(wǎng)的地震數(shù)據(jù)，其應(yīng)用都越來越廣泛。最初的震相拾取是通過長短時窗的比值檢測來進行的（Stevenson，1976），其后，Allen（1978）用其他特征函數(shù)代替原始序列提高了精度；之后又有學(xué)者提出不同形式的特征函數(shù)，這進一步提高了震相識別的精度（高淑芳等，2008；余建華等，2011）。赤池準則法最早由Sleeman等（1999）應(yīng)用于地震信號的識別，后來，劉希強等（2009）、趙大鵬等（2013）分別用偏斜度和峰度代替方差，取得了較好的結(jié)果。Schaff等（2004）對重復(fù)地震給出了實際操作層面的定義，之后模板匹配技術(shù)蓬勃發(fā)展，有研究者將其應(yīng)用到2008年汶川8.0級地震和2015年灤縣震群（最大震級ML3.2，在1976年唐山地震余震區(qū)內(nèi)）中，均取得了較好結(jié)果（呂鵬等，2011；譚毅培等，2016）。Brown等（2008）為了避免震相拾取對模板信息的依賴，提出了基于自相關(guān)的徹底盲搜索算法。Yoon等（2015）又將LSH（locality sensitive hashing）算法與震相拾取相結(jié)合，提出了FAST（fingerprint and similarity thresholding）算法，該算法避免了計算量大的問題。也有學(xué)者試圖將神經(jīng)網(wǎng)絡(luò)引入到地震學(xué)研究中來（Dai et al，1997；裴正林等，1999；王金峰等，2006）。馬強等（2013）已在福建省地震預(yù)警實驗系統(tǒng)中應(yīng)用Delaunay三角剖分法來剔除非幾何相關(guān)的干擾信號。近年來，Grigoli等（2018）的全波形疊加、何先龍等（2016）的二次方自回歸以及Lomax等（2012）的多頻率震相識別等方法均在拾取精度或運算速度上有所突破。

本文總結(jié)了目前常用的幾種震相識別技術(shù)，如長短時窗法、赤池準則法、基于互相關(guān)的模板匹配濾波法、基于自相關(guān)盲搜索的FAST法、對于S波震相拾取的偏振分析法、神經(jīng)網(wǎng)絡(luò)等，以及近幾年發(fā)展起來的全波形疊加、多頻率震相識別等在地震學(xué)中的應(yīng)用，并分析了每種方法的優(yōu)勢和局限性，以期為震相自動拾取方面的研究提供相關(guān)信息。

1 震相自動拾取技術(shù)簡介

震相拾取方法種類很多，我們主要通過2個指標來評價某個方法的優(yōu)劣：一是精度，這要求拾取震相精確可靠，也包括在低信噪比下的適用性；二是效率，要有足夠的效率以滿足需求。以下在介紹具體方法時也從這2個方面進行評價。事實上，沒有哪個方法能完美地同時滿足精度和效率，而將不同方法的優(yōu)勢結(jié)合、劣勢互補就顯得尤為重要。下面分別介紹震相的拾取過程及方法。

1.1 長短時窗法

長短時窗法最早由 Stevenso（1976）提出?？s寫為 STA／LTA法，STA為短時窗平均值，LTA為長時窗平均值，兩者之間比值的變化可以反映信號能量的變化。STA、LTA分別表示如下

圖2 4種特征函數(shù)效果（據(jù)蔣策（2017））

M、N為采樣點數(shù)；R為長短時窗比。當STA比LTA變化大時，R顯著增大，這表明有信號突變，此時通過設(shè)定好的閾值便可判斷是否為地震事件。早期的地震事件拾取工作直接利用原始的地震序列，結(jié)果并不是很好，后來選擇其他特征值代替進行計算，特征值選取可有以下幾種長短時窗法是一種操作簡單、計算速度較高的方法，它不依賴于模板，擁有廣泛的適用空間。然而，其對信噪比較為敏感，在低信噪比及高頻脈沖的影響下，長短時窗法只能拾取地震事件的大致到時，甚至?xí)板e和誤判。所以，實際工作中先使用長短時窗法得到粗略的到時來縮小范圍，再采用其他方法使結(jié)果更加精確。

1.2 赤池準則法

赤池信息準則（Akaike information criterion，AIC）最初用于統(tǒng)計學(xué)，是統(tǒng)計模型擬合度的衡量標準，由日本統(tǒng)計學(xué)家赤池弘次創(chuàng)立和發(fā)展。赤池函數(shù)表示為

Sleeman等（1999）將赤池函數(shù)引入地震學(xué)，他認為以突變時刻n作為分界點的話，地震信號可視為2個完全不同的趨于穩(wěn)態(tài)的模型，然后分別計算突變前后的AIC函數(shù)值，不斷調(diào)整n，使AIC值最小的模型則被認為是最佳模型，而此時的n即被視為初動到時

其中，N為信號總長度；M為自回歸階數(shù)；n為分界點；σ1、σ2分別為前后2個模型的方差。

引入數(shù)學(xué)中高階統(tǒng)計量與赤池準則進行結(jié)合判定，數(shù)學(xué)期望為一階統(tǒng)計量，方差為二階，偏斜度為三階，峰度為四階，表達式分別如下

劉希強等（2009）、趙大鵬等（2013）將式（10）～（13）中方差替代為偏斜度和峰度，得出在信噪比不高的情況下由偏斜度和峰度得到的結(jié)果較由方差得到的結(jié)果更精確的結(jié)論。沈統(tǒng)等（2017）利用三分量微震數(shù)據(jù)協(xié)方差矩陣對AIC方法進行約束，提高了計算精度。其中，趙大鵬等（2013）用2個時間段數(shù)據(jù)的峰度作為特征函數(shù)的方法，被稱為Kurtosis-AIC方法，表示為

其中，L為特征函數(shù)的長度；k取值為1～L。圖3為效果圖。

1.3 模板匹配法

模板匹配技術(shù)（matched filter technique）是基于滑動互相關(guān)技術(shù)發(fā)展起來的，Gibbons等（2006）在研究1997年Kara Sea地震的1次mb3.5余震的過程中發(fā)現(xiàn)，傳統(tǒng)的長短時窗方法無法檢測到的微震可以通過模板匹配方法識別出來。Kato等（2014）利用模板匹配技術(shù)對2014年智利8.1級地震的研究得到了10倍于傳統(tǒng)方法拾取的地震事件數(shù)。

模板匹配技術(shù)通過計算地震波形與模板事件的互相關(guān)來衡量二者的相似程度，從而判斷地震事件，因此對模板事件的選取很重要。以下為模板匹配法具體的流程：

（1）模板的選取。強震過后發(fā)生的余震與主震具有相近的震源機制和傳播路徑，所以波形互相關(guān)很好，選取震相清晰且信噪比較高的事件作為模板事件。

圖3 峰度作為特征值的Kurtosis-AIC方法在2009年梁山臺記錄地震中的應(yīng)用（據(jù)趙大鵬等（2013））（a）垂直向記錄；（b）圖（a）的虛框放大圖；（c）對圖（b）的濾波結(jié)果；（d）峰度曲線；（e）Kurtosis-AIC變化曲線

（3）地震事件判別。將各個臺站不同分量的互相關(guān)進行疊加平均得到平均互相關(guān)波形，將其與設(shè)定好的閾值進行比對，超過閾值的則認為是地震事件。

通過掃描在連續(xù)波形上識別出與模板事件波形相似的地震事件（圖4）。

圖4 模板匹配法識別震相示例（據(jù)李璐等（2017））藍色波形為模板事件；數(shù)值為互相關(guān)系數(shù)值

隨著模板匹配濾波方法的發(fā)展，研究者對傳統(tǒng)方法進行了不斷的改進。Meng等（2012）將連續(xù)波形分成N個部分，分別在不同線程下進行掃描，極大地提高了計算速度。Peng等（2009）同時選取P波、S波附近垂直、水平分量的時間窗，在互相關(guān)疊加時，放大模板事件波形，壓制了不相關(guān)信號，提高了信噪比。Harris（2006）、Barrett等（2014）都試圖將需要匹配檢測的模板波形多樣化，因而分別提出了子空間和經(jīng)驗子空間檢測方法，這使得模板匹配濾波方法更具有一般性。Zhang等（2015）提出了一種改進思路，其被稱為空間掃描算法，即首先在疊加之前在模板事件附近劃分三維網(wǎng)格，每個網(wǎng)格都被認為是潛在事件位置，基于潛在事件與模板事件的走時差進行校正，校正方案如下

1.4 基于自相關(guān)盲搜索的FAST法

圖5 M＆L方法與傳統(tǒng)模板匹配法對比（據(jù)Zhang等（2015））

基于互相關(guān)的模板匹配技術(shù)計算效率較高，對低信噪比的信號也有一定的適用性，但其過于依賴模板的選取，而模板的選取通常是人工下載拾取，這使其在適用性上受到限制。Brown等（2008）在研究地顫動低頻搜索問題中提出了一種基于波形自相關(guān)的檢測技術(shù)，該技術(shù)避免了選取先驗?zāi)０逍畔ⅰＦ渌悸啡缦拢菏紫葘⒚總€臺站每個分量地震序列按一定長度劃分成一定數(shù)量的時間窗，然后對這些時間窗作兩兩互相關(guān)，并每個臺站對應(yīng)互相關(guān)波形疊加平均，與設(shè)定好的閾值比對，超過閾值的則視為候選事件；之后對同一臺站的波形和這些候選事件作互相關(guān)，得到各自的互相關(guān)系數(shù)，再按照高低分組疊加來提高信噪比，最后用該高信噪比模板進行模板匹配濾波技術(shù)補充搜索（Zhang et al，2015）

其中，n為臺站號；CC為互相關(guān)值；u為時窗ti、tj下的地動位移。圖6為波形自相關(guān)方法原理圖。

基于波形自相關(guān)的檢測方法不僅與互相關(guān)運算一樣對噪聲有良好的壓制作用，而且擺脫了對模板的依賴性，其不足之處在于計算量較大，對于大數(shù)據(jù)并不適用，因而出現(xiàn)了改進算法 FAST（fingerprint and similarity thresholding）。

FAST算法思想：2個信號在哈希函數(shù)轉(zhuǎn)換前后的數(shù)據(jù)空間中都有很高的相似性，或者變換前后相似性都不高。大體思路是對原始波形作頻譜分析，提取關(guān)鍵特征，并將其壓縮成二進制“指紋”，然后采用LSH（locality-sensitive hashing）算法（即局部敏感哈希算法）將該指紋與數(shù)據(jù)庫中的二進制“指紋”作快速相似性搜索成對，最后生成地震事件檢測列表（劉翰林等，2017）。

圖6 波形自相關(guān)方法原理圖（據(jù)Brown等（2008））

1.5 對于S波拾取的偏振分析法

對于S波的震相拾取，由于受P波尾波信號的干擾，長短時窗法識別并不可靠。Roberts等（1989）提出利用P波、S波偏振方向的不同來識別震相，Amoroso等（2011）對此作了深入研究，具體方法如下：

第1步，計算三分量地震記錄的協(xié)方差矩陣，最大特征值可以認為是P波偏振方向，協(xié)方差矩陣表示如下

其中，N、E、Z分別為北向、東向、垂向分量；cov為2個分量協(xié)方差。

第2步，將方位坐標系旋轉(zhuǎn)至偏振坐標系，偏振三坐標軸分別記作 L、Q、T，分別表示 P波震動方向、Sv波振動方向、Sh波振動方向，用ψ表示射線入射角，β表示直達P波反方位角，則有

第3步，計算偏振函數(shù)，取L、Q、T三分量協(xié)方差矩陣3個特征值λ1、λ2、λ3，偏振函數(shù)表示如下

當偏振函數(shù)取最大值時，即為S波震相到時。Cichowicz（1993）提出聯(lián)合偏斜度、偏振函數(shù)、能量比3個因素進行S波信號檢測。Amoroso等（2011）用加權(quán)因子對實際值和觀測值進行應(yīng)用對比，效果顯著。Wang等（2017）對汶川8.0級地震P波、S波的拾取效果作了對比，研究發(fā)現(xiàn)P波、S波的檢出率分別為91%、85%，P波優(yōu)于S波，認為這是因為P波是具有比S波路徑特征更簡單的壓縮波，而S波的產(chǎn)生和傳播相對復(fù)雜，會在一定程度上受傳播路徑的影響。

1.6 人工智能方法

神經(jīng)網(wǎng)絡(luò)最早由心理學(xué)家和神經(jīng)學(xué)家創(chuàng)立，發(fā)展至今已成為一種非線性建模識別方法，相比于其他方法，它有很多顯著的優(yōu)點，如識別精度極高、非線性品質(zhì)良好、學(xué)習(xí)方式靈活而有效、模型結(jié)構(gòu)層次分明、存儲結(jié)構(gòu)呈完全分布式等，基于這些特點，神經(jīng)網(wǎng)絡(luò)被廣泛應(yīng)用于信號處理、模式識別及非線性優(yōu)化等領(lǐng)域，并顯示出了巨大潛力（王金峰，2006）。1986年Rumelhart等（1986）提出了經(jīng)典的BP（error back propagation）算法，將神經(jīng)網(wǎng)絡(luò)分為輸入層、隱含層、輸出層等3部分，隱含層通常不止1層。相鄰2層節(jié)點完備連接（即以上一層輸出作為輸入），同一層內(nèi)無連接，跨層無連接，每個連接都有1個權(quán)重。工作過程分為2個階段，第1個階段為學(xué)習(xí)過程，在這個過程中調(diào)整與連接相關(guān)聯(lián)的權(quán)重，使輸出滿足精度要求以達到學(xué)習(xí)的目的；然后進入第2個階段，利用學(xué)習(xí)好的神經(jīng)網(wǎng)絡(luò)進行工作，輸入地震數(shù)據(jù)，輸出識別結(jié)果。

經(jīng)典的BP算法中使用最速下降法達到收斂結(jié)果，但其具有一定的局限性，比如容易局部收斂，無法達到全局最優(yōu)。另外，經(jīng)典BP算法的學(xué)習(xí)率為固定值，始終保持不變，這就顯得不夠靈活。如果學(xué)習(xí)率選的較高，雖然收斂速度較快，但同時也會發(fā)生振蕩，導(dǎo)致不能收斂到期望值；而如果學(xué)習(xí)率選的太低，雖然可以保證收斂，但收斂速度又太慢。王金峰等（2006）提出了基于梯度的學(xué)習(xí)效率自適應(yīng)、基于誤差的學(xué)習(xí)速率自適應(yīng)等算法，提高了學(xué)習(xí)效率。宋建國等（2018）提出了級聯(lián)（cascade-correlation）相關(guān)算法，該方法可以擴展網(wǎng)絡(luò)拓撲結(jié)構(gòu)以學(xué)習(xí)新樣本，并自行確定網(wǎng)絡(luò)結(jié)構(gòu)，提高了收斂速度。

Perol等（2018）闡述了Conv Net Quake算法應(yīng)用于地震事件檢測定位的理論，其算法的核心是由8個卷積層（z1～z8）組成的前饋系統(tǒng)和1個完全連接層z來輸出，所有層包含多個通道，由二維張量表示。8個卷積層的每個通道由前一層與1個線性一維濾波器卷積得到，其表示為

其中，σ為ReLU非線性激勵函數(shù)；c、c′為輸出與輸入通道；t、t′為時間；Ci為第i層的通道數(shù)。與標準卷積的不同在于，步幅s＝2即以2個樣本的增量滑動（而非1）。因而，每層之后沿時間軸向下采樣2個單位的數(shù)據(jù)，8層之后將4×32的張量Z8矢量化為128個元素的矢量之后通過1個線性完全連接層處理得到分類值zc

W為權(quán)重系數(shù)；b為修正值。通過此線性完全連接層把各種地震信號與分類值聯(lián)系起來，特定分類值輸出代表地震事件，由此便可拾取地震事件。表1為Conv Net Quake算法的人工智能方法與 FAST、Autocorrelation算法的對比。由表1可見，前者的運算速度比 Autocorrelation快13500倍，比 FAST快 48倍；而其可靠度可與Autocorrelation相媲美，都要優(yōu)于FAST算法，說明其有很強的實用性。

表1 Conv Net Quake算法與 FAST、Autocorrelation算法結(jié)果對比（據(jù)Perol等（2018））

1.7 多頻率震相識別方法

FP（filter picking）算法是Lomax等（2012）開發(fā)的一種用于實時地震監(jiān)測和地震早期預(yù)警的算法。該算法可連續(xù)、實時、穩(wěn)定地處理寬帶信號，并避免在大地震中的過度拾取。下面介紹其原理。

首先，對原始的離散信號進行一階差分

然后進行濾波，先是2個單極高通濾波器濾波，HP（high pass）代表高通濾波

隨后單極低通濾波器濾波，LP（low pass）代表低通濾波

特征函數(shù)定義為

FP算法是對離散時間序列進行操作的，其可以運用少量沒有進行預(yù)處理的數(shù)據(jù)，甚至無需濾波和去平均，而且，F(xiàn)P算法中只有基本的算術(shù)運算（即平方根等而沒有指數(shù)、對數(shù)或變換算法等），因此運算效率較高。

1.8 全波形疊加方法

Grigoli等（2018）提出的全波形疊加方法以整個波形作為輸入數(shù)據(jù)，利用不同臺站記錄地震波形的相關(guān)性進行事件的識別和定位，類似于地震學(xué)中的偏移技術(shù)。值得一提的是，此方法在識別之前對每個臺站的數(shù)據(jù)進行區(qū)域歸一化，具體如下

其中，Tα、Tβ為時窗起止點；Tω作為緩沖窗口以彌補處理波形之間的間隙；i為臺站名；為 P／S波疊加函數(shù)；為P／S波相對整個系統(tǒng)的初至P波的相對走時。區(qū)域歸一化可視為對每個臺站的貢獻的一種加權(quán)，在壓制大值的同時定性保留其相對其他值的大小，它可以減弱尖脈沖和由于臺站與震源太近對疊加波形的影響；然后，將區(qū)域歸一化后的P／S波疊加函數(shù)作為輸入，在一定時間窗內(nèi)計算相關(guān)函數(shù)

其中，δ為狄拉克函數(shù)，此時式（31）為空間和時間的函數(shù)，下面通過提取每個時間t上的空間最大值獲得時間相關(guān)函數(shù)

在搜索空間x下計算體積V下的平均相關(guān)函數(shù)

預(yù)設(shè)閾值H1、H2，如果在某個時刻te滿足以下條件，則判別為地震事件

雙閾值的操作用以減少錯誤判定事件的發(fā)生。這種基于波形相關(guān)方法的初衷是同步實現(xiàn)事件檢測和定位，事實上，僅在拾取事件上就有不錯的效果，尤其在噪聲較大的環(huán)境下，效果更為顯著。然而，其局限性在于當臺站布設(shè)較廣而稀疏時，需要掃描大片的區(qū)域，計算量爆炸增長，不適合實時監(jiān)測與定位。

1.9 二次方自回歸模型

何先龍等（2016）采用一種二次方自回歸模型對初至附近能量變化率曲線進行二次方自回歸處理，精確拾取了P、S波初至。首先，他們將三分量離散的地震信號表示為能量變化率序列

其中，xi為地震信號；Δt為樣本的時間間隔，對式（35）的能量變化率序列應(yīng)用傳統(tǒng)的長短時窗方法進行粗略拾取，得到大致的到時區(qū)間，然后，基于二次方自回歸模型法精確拾取震相。假設(shè)理想的二次方自回歸模型方程為

其中，未知參數(shù)為ω、g、f，通過下式近似求解這3個未知數(shù)，得出其近似序列c′i（3個未知數(shù)為其系數(shù)）

其中，2s為序列長度；i為序數(shù)。由最小二乘法求解系數(shù) ω、g、f，將其代入 ci表達式（式（36）），使ci為最大值時對應(yīng)的i時刻即為精確到時。表2為何先龍等（2016）在汶川8.0級地震研究中將二次方自回歸模型方法與傳統(tǒng)的長短窗法、赤池準則法進行比較得到的結(jié)果。由表2可見，雖然在計算效率上稍遜色，但二次方自回歸模型方法對地震事件檢測的精度與準確率較高，不失為一種優(yōu)秀的方法。

表2 二次方自回歸模型拾取與經(jīng)典方法對比表（據(jù)何先龍等（2016））

2 討論和結(jié)論

相對于早期人工拾取地震事件的低效率、低精度，20世紀以來的各種自動拾取方法為地震學(xué)研究的發(fā)展做出了巨大貢獻。Tamaribuchi等（2016）在對2016年7.1級熊本地震的研究中對自動拾取和人工拾取作了對比（圖1），其對微地震拾取數(shù)量是人工拾取的14倍，不僅在速度上有強大優(yōu)勢，而且數(shù)量上更遠優(yōu)越于人工拾取。

震相自動拾取的每種方法各有優(yōu)劣，基于長短時窗一類的方法計算簡單，速度快，不依賴先驗的模板事件信息，但是在面對低信噪比的數(shù)據(jù)記錄時效果不佳，閾值的設(shè)定很難，閾值過大可能導(dǎo)致漏檢，過小可能導(dǎo)致誤檢。基于波形互相關(guān)模板匹配的方法能有效壓制背景噪聲，在處理低信噪比的信號時有優(yōu)勢，但對模板的選取過于依賴，在地下結(jié)構(gòu)十分復(fù)雜的情況下模板選取不具代表性，選取困難。基于波形自相關(guān)的搜索方法較好地解決了上述2個問題，但其計算量過大，易導(dǎo)致實際應(yīng)用效果不佳，即使引入了數(shù)據(jù)挖掘中的優(yōu)化算法，依然有著龐大計算量。神經(jīng)網(wǎng)絡(luò)應(yīng)用于地震震相拾取對噪聲具有較高的承受能力，這使其占有比較重要的地位，在Perol等（2018）的研究中拾取準確率近100%，準確率得以保障。但是，該方法需要較長的訓(xùn)練時間，而且需要大量依靠經(jīng)驗的參數(shù)，這使得其在可解釋性方面較差。近幾年發(fā)展起來的多頻率震相識別、全波形疊加、二次方自回歸模型等方法都從現(xiàn)代的視角分析問題，盡管有一些問題（如運算速度等）有待解決，但在實際應(yīng)用中都可以發(fā)揮一定的作用，如二次方自回歸模型（何先龍等，2016）用于拾取汶川8.0級地震震相時已達到P波誤差大于0.1s的事件數(shù)僅為1個，S波誤差大于0.1s的僅為3個，比傳統(tǒng)的長短時窗法等已有了很大的突破。

一種震相自動拾取方法很難同時兼顧精度和效率，因此，合理運用每種方法的優(yōu)勢，規(guī)避劣勢可以在提高震相識別效率的同時又提高精度。比如，可以將長短時窗法與AIC方法結(jié)合起來，先用前者粗略拾取到時，然后設(shè)置時間窗用后者使到時更精確?？傊?，可以通過聯(lián)合運用各種方法來綜合識別震相，從而得到滿意的效果。