李璐 王寶善 侯金欣

中國地震局地球物理研究所（地震觀測與地球物理成像重點實驗室），北京市民族大學南路5號 100081

0 引言

由于受背景噪聲水平、記錄儀器分布等多種因素的影響，許多微弱的地震信號被淹沒在噪聲中，無法直接觀測。另一方面，大地震發(fā)生后的短時間內，由于受強震尾波的影響以及大量余震事件的發(fā)生，波形相互重疊，難以區(qū)分。如果能夠自動識別出這些微弱的地震信號和遺漏的早期余震，則對于完善地震目錄、全面分析地震活動性、監(jiān)測地震破裂區(qū)域的震后形變等具有重要意義。

目前，在大量地震自動識別方法中，廣泛應用的主要有2種：一種是工業(yè)界普遍使用的長短時窗比（Short-Term-Average／Long-Term-Average，STA／LTA）方法（Stevenson，1976；Earle et al，1994）。該方法利用包絡線上短時窗內信號的平均值與長時窗內信號平均值的比值來表征信號振幅的變化，其中，短時窗用于截取微弱的有用信號，長時窗內信號的平均值則反映背景噪聲水平。雖然算法簡單，計算速度快，但對于信噪比過低的信號，這種方法并不適用。另一種方法是基于互相關技術的模板匹配濾波技術（Matched Filter Technique，MFT）（Gibbons et al，2006；Shelly et al，2007；Peng et al，2009）。該方法在滑動窗互相關（Sliding-Window Cross-Correlation，SCC）檢測技術的基礎上發(fā)展起來（Yang et al，2009），是在信噪比較低的情況下提取微弱信號的一種有效方法。其中，Gibbons等（2006）分析了1997年8月16日Kara Sea地震主震后4hr發(fā)生的1個微地震，發(fā)現(xiàn)傳統(tǒng)的STA／LTA方法無法檢測出此次事件，而模板匹配濾波技術在單通道記錄上即可將其分辨出來。由此可見，模板匹配濾波技術在微地震的自動檢測方面具有明顯優(yōu)勢。

除了地震學，模板匹配濾波技術還在其他領域得到廣泛應用。2015年曾在科學界引起轟動的重大發(fā)現(xiàn)——引力波的探測（Abbott et al，2016）中，科學家們利用干涉疊加提高探測器精度，于2015年9月14日09：50：45（UTC）探測到引力波信號，該信號極其微弱，應變約為1.0×10－21。為了進一步驗證該微弱信號的來源，科學家們以黑洞模型模擬出理論引力波信號作為模板，在實際觀測的連續(xù)波形上掃描，發(fā)現(xiàn)了與模板波形一致的信號，且時間上與實際觀測到的引力波信號一致，從而驗證了引力波信號的可靠性。Zhao等（2016）以爆炸事件波形為模板掃描連續(xù)波形，并未在2015年8月12日天津危險品倉庫爆炸前發(fā)現(xiàn)小爆炸事件。此外，在醫(yī)學、天體物理、電子通信等領域，模板匹配濾波技術還是用于造影、圖像識別、語音識別等的重要手段（Hoover et al，2000；Dong et al，2008；Avadhanulu et al，2013）。

本文主要介紹模板匹配濾波技術的基本原理和數據處理流程。從尋找大地震的前震與余震、遠程動態(tài)觸發(fā)地震、誘發(fā)地震和非火山震顫等方面，詳述該方法在地震數據處理中的應用。此外，對于模板匹配濾波技術的改進與發(fā)展也進行了簡要討論。

1 MFT基本原理

模板匹配濾波技術（MFT）利用現(xiàn)有地震目錄中記錄到的事件作為模板，掃描連續(xù)波形，識別出與模板事件波形具有一定相似度的地震事件。匹配濾波檢測微地震的處理流程如圖1所示。

圖1 模板匹配濾波技術數據處理流程

第1步，選取合適的頻段對模板事件和連續(xù)波形作同樣的濾波處理。微地震檢測一般為區(qū)域地震，因此，選取的頻段應相對高頻。

第2步，標定每個模板事件三分量波形上的P波、S波到時。標定方法可以根據已知速度模型計算理論到時，也可人為手動拾取。再對每一個模板事件的三分量分別計算P波、S波的信噪比。其中，信號項選取P波或S波到時后某一時間窗內的波形計算均方根振幅；噪聲項選取在P波到時前同樣長度的時間窗內的波形計算均方根振幅。分別用P波或S波的信號項除以噪聲項，得到對應的P波或S波信噪比。在后面的波形掃描中，只選取信噪比高于某一閾值的模板事件波形。

第3步，對于某一模板事件，利用下式計算模板事件的波形和對應分量的連續(xù)波形的互相關系數（Cross-Correlation coefficient，CC系數））

其中，t0、t1分別為用于計算互相關系數的窗口的開始與結束時刻；X（t）、Y（t）分別為 t0到 t1時間段內的模板事件波形和連續(xù)波形。計算時每次在連續(xù)波形上移動1個采樣點，最終得到該分量上的一道互相關波形。

第4步，根據模板事件上P波或S波到時與該事件發(fā)震時刻之間的到時差，將每道互相關波形平移至發(fā)震時刻，再對所有臺站、所有分量的互相關波形進行疊加，取平均，得到該模板事件的平均互相關波形。

第5步，設定互相關系數的閾值，當平均互相關波形上某一時刻的互相關系數高于該閾值時，就認為這是一個新檢測事件。選取閾值時，一般選擇絕對中位差（Median Absolute Deviation，MAD）的 9倍（Peng et al，2009）或 12倍（Yao et al，2015）。但是，當臺站分布稀疏、波形信噪比較低時，可以根據需要適當提高閾值（如15倍MAD），以減少模板事件與噪聲信號間的錯誤檢測。圖2為利用MFT以2015年4月26日20：08：33發(fā)生于中國藏南地區(qū)的1次2.40級地震為模板，檢測到發(fā)生于4月25日19：43：39的1次1.81級地震。如圖2所示，通過掃描在連續(xù)波形上可以識別出與模板事件波形相似、沒有在地震目錄上出現(xiàn)的地震事件?；谥貜偷卣鸬奶匦?，認為新檢測事件位置與對應的模板事件位置一致。根據模板事件和新檢測事件在連續(xù)波形上的到時差，可以計算出新檢測事件的發(fā)震時刻。根據某一時間窗內新檢測事件與模板事件的振幅比，可確定新檢測事件的震級。

2 M FT在地震學中的應用

在算法逐漸完善的基礎上，模板匹配濾波技術已經在多種地震事件的檢測中得到應用。

2.1 運用M FT尋找前震和余震

許多大地震發(fā)生前會有1個或多個前震事件發(fā)生，但前震與主震成核過程之間的聯(lián)系還存在爭議。建立更加完整的前震活動序列，有助于加深對地震發(fā)生物理過程的理解，為地震預測提供指導（Jones et al，1979；Ben-zion，2008）。Kato等（2012）通過波形互相關的方法，在2011年日本東北9.0級地震之前識別出2個相互分離的前震地震群，并發(fā)現(xiàn)它們在時間上有向主震初始破裂點緩慢移動的趨勢。通過對滑動速率的分析發(fā)現(xiàn)，第2個前震序列產生主要的應力加載，促進主震的破裂過程。Kato等（2014a）對2014年智利8.1級地震發(fā)生前后的連續(xù)波形進行模板掃描，新檢測到數量10倍于原始地震目錄的地震事件，并發(fā)現(xiàn)前震序列中慢滑移事件同樣有向著主震成核點遷移的趨勢。另一方面，Wu等（2014）對1999年土耳其7.1級地震前65hr內的連續(xù)波形進行了匹配濾波分析，并未發(fā)現(xiàn)前震地震活動對主震破裂的促進作用，這說明前震對主震成核的促進作用以及對主震發(fā)震位置的影響并不是普遍現(xiàn)象。

大地震發(fā)生后的數小時內，在距主震震中幾個破裂長度的范圍內，地震活動性明顯升高，這些被觸發(fā)的地震稱為余震。清楚認識早期余震的活動特性，是分析地震觸發(fā)機制的重要基礎（Helmstetter et al，2009）。模板匹配濾波技術有助于識別出那些淹沒在主震尾波中的早期余震事件，以完善余震序列。Peng等（2009）將該方法運用到2004年帕克菲爾德6.0級地震發(fā)生后的3天內，共檢測到數量11倍于地震目錄的余震事件（圖3）。檢測結果有效補充了早期余震序列，有助于清晰地刻畫出余震的遷移特性。

圖2 模板匹配濾波檢測地震示意圖

同樣運用模板匹配濾波技術，Lengliné等（2012）在2011年日本東北9.0級地震之后的12hr內檢測到數量12倍于原始地震目錄的事件，通過這些新檢測事件在時間和空間上的分布，得到余震衰減速率（p＝1.0）和余震區(qū)沿斷層走向的擴展趨勢；Meng等（2013）通過匹配濾波后更加完善的地震目錄發(fā)現(xiàn)，2003年美國圣西蒙6.5級地震后，近場的觸發(fā)地震與由主震導致的靜態(tài)剪切應力變化間呈正相關，說明在圣安德烈亞斯斷層帕克菲爾德地區(qū)存在可能的低摩擦作用；Meng等（2014）通過對比2個距主震不同震中距的區(qū)域內的地震活動性認為，靜態(tài)、動態(tài)應力的改變均對近場地震的觸發(fā)有影響，其中，靜態(tài)庫侖應力可以作用于長期地震活動性的改變，動態(tài)應力則主要影響更大區(qū)域尺度內的瞬時地震活動性的改變；Kato等（2014b）在2007年日本能登半島6.7級地震發(fā)生后的32天內，同樣觀測到新檢測的余震事件沿破裂斷層的走向擴展，同時發(fā)現(xiàn)該擴展趨勢呈階梯狀，認為可能與主震發(fā)生后的抗震余滑（a seismic after slip）有關。

圖3 對2004年帕克菲爾德6.0級地震后余震序列的模板匹配濾波檢測結果（據Peng等（2009））

2.2 運用M FT尋找遠程動態(tài)觸發(fā)地震

除了在近場帶來應力改變并產生觸發(fā)地震之外，大地震還可以傳播到幾百至上千千米處，在遠場觸發(fā)一些地方震或慢地震（Peng et al，2010；Hill et al，2015）。這種觸發(fā)多與大振幅面波引起的本地動態(tài)應力擾動有關，因此被稱為動態(tài)觸發(fā)。選取合適頻段濾波后，人為手動拾取地方震到時的方法工作量較大，不利于系統(tǒng)觀測。已有研究表明，火山活動區(qū)由于巖漿入侵等流體運動，容易達到應力的臨界狀態(tài)，是尋找動態(tài)觸發(fā)現(xiàn)象的重點區(qū)域（Moran et al，2004；Prejean et al，2004）。運用模板匹配濾波技術，在2011年日本東北9.0級地震的面波持續(xù)期間，Yukutake等（2013）在距主震約450km的箱根火山區(qū)域檢測到由主震觸發(fā)的地方震。另一方面，板塊交界區(qū)域由于需要密切監(jiān)測地震活動，臺站分布密集，有利于動態(tài)觸發(fā)現(xiàn)象的觀測，而對于內陸板塊的動態(tài)觸發(fā)研究則較少。Yao等（2015）通過模板匹配濾波發(fā)現(xiàn)，在西藏中南部由2004年蘇門答臘9.1級地震和2005年尼亞斯8.6級地震動態(tài)觸發(fā)的地方震，并對2次大地震動態(tài)觸發(fā)的持續(xù)時間和觸發(fā)區(qū)域進行了對比。

2.3 運用M FT尋找誘發(fā)地震

在遠場觸發(fā)地震活動中，由人為活動導致的誘發(fā)地震是一個重要組成部分。其中，非常規(guī)天然氣能源——頁巖氣的勘探開發(fā)工作最早起源于美國，并隨著水壓致裂等多項技術的成熟在多個國家和地區(qū)推廣使用。而在該技術迅速發(fā)展的過程中，與生產活動和提高產量有關的廢水注入地下的處理方式是否直接導致了當地地震活動性的增加，一直是科學界和工業(yè)界爭論的話題（Ellsworth，2013）。在實際觀測中，注水帶來的應力擾動十分微弱（約1MPa）（Zoback et al，1997），可觸發(fā)的地震震級較小，不易觀測。運用模板匹配濾波技術可以實現(xiàn)對當地微小地震的更加全面的探測，從而有助于判斷地震活動性變化與注水活動之間的關聯(lián)。van der Hilst等（2013）利用模板匹配濾波技術尋找美國中西部地區(qū)的觸發(fā)地震，發(fā)現(xiàn)有廢水注入的地區(qū)，巖石流體孔隙壓增加，斷層系統(tǒng)更接近臨界狀態(tài)，當有遠場的大地震發(fā)生時，更容易在應力擾動下產生觸發(fā)地震。Walter等（2016）同樣利用模板匹配濾波技術進一步研究了美國得克薩斯州東部和路易斯安那州西北部的地震活動性，新檢測到的40個地震事件在空間上可分為幾個地震群，雖然沒有強烈的證據表明當地地震活動增加與廢水注入間有直接關系，但也并不能完全否定兩者間的相關性，作者認為該地區(qū)需要更多的地震學監(jiān)測。在我國，也有一些研究工作將模板匹配濾波技術應用于微小誘發(fā)地震事件的檢測中。Wang等（2015）在北京房山地區(qū)于2011年日本東北9.0級地震和2012年印度洋8.6級地震發(fā)生后的短時間內，分別新檢測到分布集中的大量淺部地震，并認為可能與當地采礦活動帶來的地下應力變化有關。

2.4 運用M FT尋找非火山震顫

非火山震顫（Non-Volcanic Tremor，NVT）主要發(fā)生在俯沖帶內，由許多相互重疊的低頻事件（Low Frequency Event，LFEs）組成（Shelly et al，2007、2011），具有持續(xù)時間長、頻率較低、體波震相不清晰等特點（Beroza et al，2011），其信噪比較低，不易觀測。Shelly等（2009）在美國帕克菲爾德地區(qū)第1次利用波形互相關觀測到LFEs，并對這一組LFEs組成的震顫進行定位，發(fā)現(xiàn)圣安德烈亞斯斷層可能延伸至地殼底部，比震源深度最大的普通地震還要深約10km。利用模板匹配濾波在震顫中發(fā)現(xiàn)更多重復信號（LFEs），有助于更好地刻畫俯沖帶內斷層的破裂方向與可能的觸發(fā)機制（Aguiar et al，2009；Brown et al，2010）。

3 M FT方法的改進與發(fā)展

在模板匹配濾波技術的基本原理之上，一些研究工作對其算法進行了優(yōu)化。首先，模板匹配濾波是在連續(xù)波形上逐個采樣點掃描，運算量大，耗時長。因此，可以通過降采樣的方式，適當縮短掃描的時間。此外，Meng等（2012）運用圖形加速卡（Graphic Processing Unit，GPU）的并行計算方法，將連續(xù)波形分成幾個區(qū)塊，每個區(qū)塊在GPU的不同線程下同時進行匹配濾波掃描。若同時啟用N個線程，則計算速度可以提高到原來的N倍。其次，前人研究工作中，對模板事件波形的垂直分量和水平分量分別選取P波、S波到時附近的時間窗（Peng et al，2009），當研究區(qū)域內臺站分布稀疏時，可以如圖2所示，對模板事件波形的每一分量同時選取P波和S波附近的時間窗，則模板事件波形的數量翻倍，更有利于互相關波形疊加時對不相關信號的壓制。最后，與重復地震的特征類似，模板匹配濾波的前提假設是，新檢測到的事件與模板事件的發(fā)震位置一致。然而，在實際觀測中，有些遺漏事件可能具有與模板事件高度相似的波形，但位置上略有偏差，這導致它們傳至臺站的用時略有不同，此時，如果將各臺站的互相關波形按模板事件的到時平移至發(fā)震時刻，則會導致該遺漏事件的互相關信號沒有完全對齊，此時的疊加波形上互相關系數減小，降低了該事件被檢測到的可能性。Zhang等（2015a）發(fā)展了一種Match and Locate（M＆L）改進算法，在模板事件附近的1個小三維區(qū)域內進行網格搜索，將新檢測事件定位在每一個網格點上，計算此時該事件和對應的模板事件之間的到時差，將各道互相關波形按對應的到時差平移后再疊加，最后比較所有網格點的疊加互相關系數和信噪比，以確定事件的位置。由于是在小范圍內進行空間搜索，因此，對速度模型的依賴性較小。該方法不僅實現(xiàn)了對新檢測事件的重定位，而且提高了被檢測事件的互相關系數，從而增加了更小震級事件被探測到的可能性。Zhang等（2015c）利用此方法監(jiān)測到朝鮮地區(qū)的多次核爆實驗，因而將全球的探測核實驗能力提高到噸的級別。運用同樣的方法，Zhang等（2015b）在日本Ontake火山2007、2014年2次噴發(fā)前均探測到密集的微震活動，為火山活動的預測提供了依據。

4 結語

通過模板匹配濾波技術可以在背景噪聲水平較高時利用互相關掃描的方法尋找到肉眼難以觀測到的信號，該方法適用于各類微弱信號的檢測。在地震學中，運用該方法監(jiān)測地震活動性對于斷層走向的刻畫、地震破裂過程模擬、地下速度結構成像等方面的研究具有重要意義。基于模板匹配濾波技術的自動掃描、識別信號的特征，今后工作中，可將該技術應用于超密集臺陣的自動化處理，以提高地震監(jiān)測工作的效率。于2015年全面啟動的國家地震烈度速報與預警工程，將在全國范圍內建成由5000余個強震臺站組成的密集觀測臺網用于災情的快速分析，并指導抗震救災工作。利用模板匹配濾波技術，將地震事件的自動識別與定位相結合將是處理這些海量地震數據的有效手段。