天然地震識別與震相自動拾取技術進展

孫印 潘素珍 劉明軍

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

0 引言

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

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

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

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

1 震相自動拾取技術簡介

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

1.1 長短時窗法

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

圖2 4種特征函數效果（據蔣策（2017））

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

1.2 赤池準則法

赤池信息準則（Akaike information criterion，AIC）最初用于統計學，是統計模型擬合度的衡量標準，由日本統計學家赤池弘次創立和發展。赤池函數表示為

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

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

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

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

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

1.3 模板匹配法

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

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

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

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

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

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

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

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

1.4 基于自相關盲搜索的FAST法

圖5 M＆L方法與傳統模板匹配法對比（據Zhang等（2015））

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

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

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

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

圖6 波形自相關方法原理圖（據Brown等（2008））

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

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

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

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

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

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

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

1.6 人工智能方法

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

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

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

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

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

表1 Conv Net Quake算法與 FAST、Autocorrelation算法結果對比（據Perol等（2018））

1.7 多頻率震相識別方法

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

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

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

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

特征函數定義為

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

1.8 全波形疊加方法

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

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

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

在搜索空間x下計算體積V下的平均相關函數

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

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

1.9 二次方自回歸模型

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

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

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

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

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

2 討論和結論

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

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

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