田 平，蔣 策，陳 瑛，梁 明

（1.廣東省地震局，廣州 510070；2.中國地震局地震監(jiān)測與減災技術(shù)重點實驗室，廣州 510070）

0 引言

廣東省新豐江水庫自1959 年截流蓄水以來，水庫區(qū)內(nèi)的地震活動明顯增加，并在1962 年觸發(fā)了6.1 級地震，成為華南地區(qū)地震活動最為活躍的地區(qū)之一。為研究庫區(qū)的地震活動特征及發(fā)震機制，廣東省地震局在新豐江地區(qū)開展了大量的研究工作，不斷加密測震臺網(wǎng)，目前該地區(qū)的水平定位誤差小于1千米，達到百米級別的定位能力[1]。水庫區(qū)密集的測震臺站可以高質(zhì)量記錄到區(qū)域內(nèi)發(fā)生的大量小震微震事件，對這些地震活動及其時空分布特征進行精確分析，完善地震目錄，有利于對新豐江地區(qū)的發(fā)震斷層結(jié)構(gòu)的精細探測、水庫蓄水對地震的觸發(fā)作用以及地震預測等方面研究工作的開展。

波形模板匹配（matched filter）方法是一種有效的微震檢測方法，該方法將模板地震的各分量波形與連續(xù)波形進行互相關(guān)計算，通過對互相關(guān)波形進行疊加來判斷地震事件[2-4]。波形模板匹配方法對微弱的地震信號敏感，即使在低信噪比的情況下也可以檢測到震級極小的微震事件，廣泛應用于地震前震、余震檢測[5-8]、觸發(fā)地震[9-10]以及非天然地震檢測[11-12]等方面。Zhang 等[13]基于波形模板匹配方法發(fā)展了匹配定位方法（Match and Locate，簡稱M&L 方法），在相關(guān)波形疊加前考慮待檢測事件與模板事件之間的相對位置，在檢測微震的同時，獲取高精度的發(fā)震位置信息。Zhang等[14]利用該方法檢測了日本Ontake 火山2007 和2014 年兩次噴發(fā)前的密集微震活動，得到30 余倍日本氣象廳的地震目錄結(jié)果；馬曉靜等[15]將該方法應用于2014 年6 月河源地區(qū)的小震群事件研究中，結(jié)果顯示模板檢測到的地震事件數(shù)目相較于地震目錄結(jié)果增加了一倍；王志偉等[16]利用匹配定位方法對2008 年10 月至2011 年7 月間重慶榮昌地區(qū)的微震進行檢測與定位，將該地區(qū)地震目錄的震級下限從ML1.0 降低至ML0.3；Liu 等[17]在匹配定位方法的基礎上，開發(fā)了GPU 版本的匹配定位計算程序，顯著地提高了模板匹配計算的速度；杜瑤等[18]在水庫區(qū)爆破事件的檢測研究中，發(fā)現(xiàn)匹配定位方法可以有效區(qū)分不同類型的地震活動，進而提高對疑爆事件的識別率。

在本研究中，我們利用M&L 方法檢測新豐江河源地區(qū)的微震事件，通過數(shù)字地震臺網(wǎng)中心數(shù)據(jù)處理系統(tǒng)建立適合新豐江地區(qū)的微震模板事件庫，實現(xiàn)對地震臺網(wǎng)記錄的連續(xù)波形的自動檢測，將其應用于臺網(wǎng)的地震目錄完善工作。

1 原理與方法

M&L 方法不需要假定待檢測地震與模板地震發(fā)生在同一位置。該方法將各臺站分量上記錄的模板波形與包含潛在地震信號的連續(xù)波形進行互相關(guān)計算，在對相關(guān)波形進行疊加前，考慮待檢測地震與模板地震發(fā)生位置之間的差異，在模板地震周圍的三維空間進行網(wǎng)格搜索，計算模板地震與每個假定待檢測地震位置之間的走時差，在此基礎上對互相關(guān)波形進行走時矯正和疊加，計算疊加后的互相關(guān)波形中的平均互相關(guān)（Cross-Correlation，CC）系數(shù)和信噪比（Signal-Noise ratio，SNR）。

兩個距離很近的地震波形之間的歸一化互相關(guān)可表示為：

互相關(guān)系數(shù)與地震波形的相似度成正比，可以作為事件檢測的閾值標準。但僅使用較低互相關(guān)系數(shù)作為閾值時，盡管可以增加檢測事件的數(shù)量，同時也會極大提高了事件的誤檢概率；信噪比的引用可以對其進行約束，在較低閾值標準設定的情況下，檢測到更多地震，同時避免增加誤檢概率。

M&L 方法將平均互相關(guān)系數(shù)和信噪比二者聯(lián)立作為地震事件的判斷標準。當其值超過設定的相關(guān)系數(shù)閾值和信噪比閾值時，即判斷為檢測到地震事件。之后，根據(jù)所有臺站三分量波形與參考事件波形的振幅比的中位數(shù)確定地震事件的震級。

2 研究內(nèi)容

2.1 河源地區(qū)微震事件檢測

2018—2019 年期間，廣東測震臺網(wǎng)記錄到新豐江地區(qū)發(fā)生的天然地震共2967件（圖1），其中震級ML≤2.0 的地震事件共2935 件，震級ML≤1.0 的地震事件共2625 件，微震小震事件占新豐江地區(qū)地震事件總數(shù)的88%以上。

圖1 2018—2019年間新豐江地區(qū)地震目錄Fig.1 Earthquake catalog in Xinfengjiang area from 2018 to 2019

上述地震目錄全部為人工瀏覽連續(xù)地震波形識別、編目得到，由于人工識別誤差、震級未達到編目要求（ML＜0.0）等原因，還存在大量微震事件未被記錄在內(nèi)。本研究將地震目錄中的部分地震事件作為地震模板，檢測2018—2019 年間新豐江河源地區(qū)的地震事件，將其與人工地震目錄結(jié)果比較，測試M&L方法在新豐江河源地區(qū)的適用性。

M&L 方法的檢測結(jié)果依賴于地震模板的豐富程度，無法利用有限的模板識別出所有的微震事件。地震模板的種類和數(shù)量越豐富，檢測的結(jié)果越豐富越接近真實情況。相應地，模板的大量使用也極大增加計算時間，進而降低檢測效率。在使用模板匹配檢測方法時，一般將研究時間范圍內(nèi)的所有事件進行篩選作為地震模板，對連續(xù)波形進行掃描。新豐江庫區(qū)地震數(shù)量較多，為測試在該地區(qū)模板檢測的能力及效率，對使用的地震模板進行了人為震級篩選。最終確定選定2019 年間新豐江河源地區(qū)（23.700°～23.766°N，114.57°～114.75°E），震級范圍0.5≤ML≤1.0，共125條地震事件作為地震模板（圖2中紅點所示）。新豐江水庫區(qū)測震臺站較為密集，經(jīng)過對臺站數(shù)據(jù)的記錄質(zhì)量以及相對研究區(qū)域的位置分布，選取其中的10 個臺站（圖4）記錄的數(shù)據(jù)進行分析處理。M&L 方法受地殼速度模型影響較小，本文采用PREM 速度模型[20]，利用TauP 工具包[23]計算理論震相到時。微震事件的判斷閾值參考采用9 倍背景平均相關(guān)系數(shù)。另外，新豐江水庫區(qū)的地震深度普遍較淺，深度分布比較集中，為提高檢測效率，在地震檢測中沒有對震源深度進行搜索，固定深度為模板地震的深度。

圖2 2018—2019年間利用模板檢測得到的地震事件M-T分布結(jié)果Fig.2 M-T distribution results of earthquake events obtained by template detection from 2018 to 2019

通過對連續(xù)波形的模板匹配掃描，共檢測到2018—2019 年間新豐江河源地區(qū)的地震事件共2160 條（圖2），最大震級ML3.08（圖3），最小震級ML0.8。對檢測結(jié)果的震級分布進行統(tǒng)計分析，震級0.0≤ML≤1.0 的微震事件共625 條，占地震事件總數(shù)的28.9 %；震級ML＜0.0 的地震事件共1438 條，占地震事件的66.5%。檢測到的地震事件在空間上分布在地震模板的發(fā)震位置附近。2018—2019 年地震目錄結(jié)果（圖4）顯示新豐江河源區(qū)域記錄地震共1764 件，震級0.0≤ML＜0.5 的微震事件共1163條，震級0.5≤ML≤1.0 的微震事件共425 條，震級1.0≤ML≤2.0 的地震事件共159 條，ML＞2.0 的地震事件17條。

圖3 模板檢測結(jié)果示例（最大震級為ML 3.08，發(fā)震時刻2018-07-25 08：53：45）Fig.3 Example of template detection results（the maximum magnitude is ML 3.08，which occurred at 08：53：45 2018-07-25）

上述僅利用少量震級0.5≤ML≤1.0 的地震模板，對2018-2019年間新豐江河源地區(qū)的地震事件的檢測結(jié)果表明，地震模板事件全部被自檢出，相關(guān)系數(shù)接近于1.0。將這些事件進行比較，發(fā)現(xiàn)檢測結(jié)果的發(fā)震時刻和震級與編目結(jié)果具有較小的差異，具有較高的準確度。但是，由于缺少2018 年0.5≤ML≤1.0 的地震模板，對2018 年期間檢測出的震級0.5≤ML≤1.0 的地震事件僅92 件，少于廣東臺網(wǎng)編目結(jié)果300條，M&L方法的檢測結(jié)果受地震模板的豐富程度影響較大，僅利用部分模板無法對所有地震事件進行檢測。

2.2 檢測程序改進

上述的模板檢測結(jié)果表明，M&L 方法適用于新豐江水庫區(qū)的微震檢測，而且可以對檢測到的微震事件進行定位以及震相標記，可將其應用于測震臺網(wǎng)的地震編目工作中。在數(shù)字臺網(wǎng)中心的服務器上進行軟件部署，利用數(shù)字臺網(wǎng)中心數(shù)據(jù)處理軟件系統(tǒng)JOPENS[21-22]獲取地震事件以及連續(xù)波形，實時擴充地震模板庫，并且對其調(diào)用在連續(xù)波形上自動進行微震檢測，達到完善地震目錄的目的。

根據(jù)M&L 模板匹配計算程序的處理流程（圖5），為達到對地震波形的自動檢測的要求，需要改進的主要部分為流程中最主要的輸入部分：即“連續(xù)波形”和“地震模板”。

圖5 模板檢測程序的流程Fig.5 Process of template detection program

對于“連續(xù)波形”部分，通過數(shù)據(jù)處理系統(tǒng)獲取到連續(xù)波形數(shù)據(jù)后，需要對其進行格式轉(zhuǎn)換，完成波形的拼接、濾波以及降采樣等數(shù)據(jù)處理工作。我們利用JOPENS 的數(shù)據(jù)導出功能，將其設置為定時將新入庫的連續(xù)波形數(shù)據(jù)導出為seed 文件，保存在臺網(wǎng)中心服務器端。檢測程序部署在信息中心機群，利用腳本從臺網(wǎng)中心服務器下載波形數(shù)據(jù)seed 文件，并自動將其轉(zhuǎn)換為sac 格式數(shù)據(jù)文件；對連續(xù)波形進行1～20 Hz 的濾波以及降采樣處理工作，統(tǒng)一命名放置于M&L 檢測程序內(nèi)部的Trace目錄，每日對其進行更新。

對于“檢測模板”部分，臺網(wǎng)編目報告的地震事件seed 文件以及地震目錄catalog 文件經(jīng)數(shù)據(jù)處理系統(tǒng)導出上傳至中心機群，放置于目標目錄下。系統(tǒng)定時在該目錄下進行搜索，發(fā)現(xiàn)地震事件seed 文件并將其轉(zhuǎn)換為sac 格式數(shù)據(jù)文件，對事件波形進行1-20 Hz 的濾波以及降采樣處理，保證與連續(xù)波形數(shù)據(jù)格式相同。經(jīng)上得到的新的地震模板文件，統(tǒng)一命名放置于M&L 檢測程序內(nèi)部的Template目錄，實時對地震模板庫進行更新。

上述利用震級范圍0.5≤ML≤1.0 的地震模板進行檢測的結(jié)果表明，僅利用部分地震模板不能檢測出潛在的所有地震，地震模板的豐富程度對檢測結(jié)果產(chǎn)生非常大的影響。在保證檢測效率的同時，為了檢測盡量多的微震事件，在部署自動檢測軟件中，將待檢測長度為24 小時的連續(xù)波形包括其前24小時以及后24小時，共72小時內(nèi)臺網(wǎng)編目的所有地震事件作為模板，對連續(xù)波形進行檢測。若此時間段內(nèi)不存在地震事件，自動將調(diào)取時間范圍進行延長；另一方面，系統(tǒng)也可隨時根據(jù)需要調(diào)用數(shù)據(jù)庫中的其他地震模板，以保證充足的地震模板進行微震檢測。自相關(guān)系數(shù)閾值等參數(shù)設置采用與上文相同設置。

2.3 河源M 4.3地震余震檢測

基于對程序的上述配置，重新檢驗其在新豐江地區(qū)地震檢測中的適用性。2023 年2 月11 日10時41分廣東省河源市源城區(qū)發(fā)生了M4.3地震，震后發(fā)生一系列余震，廣東數(shù)字地震臺網(wǎng)中心的編目結(jié)果顯示，截至2月12日0時，共記錄到該地區(qū)余震171 件，其中，0≤ML＜1.0，143 件；1.0≤ML＜2.0，23件；ML≥2.0，5件；最大余震ML3.4。

利用上述改進的程序配置，對其在該時間段內(nèi)對連續(xù)波形的自動檢測結(jié)果進行分析。結(jié)果顯示，截至2023 年2 月12 日0 時，自動檢測配置共檢測到河源M4.3 地震的余震265 件，與人工余震目錄相比，新檢測到的地震事件中，有15 件事件的震級在ML0.0 左右，最大的檢測震級為ML0.64；其余震級均在ML＜0.0（圖6），最小的震級為ML-0.65。

圖6 檢測出的地震事件（ML 0.37）Fig.6 One detected earthquake event with magnitude of ML0.37

地震目錄記錄的171件地震事件，共檢測出其中的167件，經(jīng)人工對波形校對發(fā)現(xiàn)，未檢測出的4 件地震事件均為雙震事件，震級在0≤ML＜1.0 范圍。通過對模板的檢查發(fā)現(xiàn)，雙震中的兩起事件發(fā)生間隔時間過短（圖7），生成的模板同時包含了兩起事件，在連續(xù)波形計算過程中對另一事件無法檢測，導致對雙震事件的檢測失效。

圖7 未檢測出的雙震事件Fig.7 Undetected double-earthquake event

將地震模板中被重新檢測出中被地震目錄收錄的167件地震事件的發(fā)震時刻、震級與人工結(jié)果進行比較（圖8），發(fā)現(xiàn)大部分震級與人工編目的震級結(jié)果相近，差值一般小于ML0.5；發(fā)震時刻與人工編目的發(fā)震結(jié)果相近，差值小于1.0 s。

圖8 編目報告與檢測結(jié)果中同事件的震級（a）與發(fā)震時刻（b）的差值統(tǒng)計比較Fig.8 Statistical comparison of the difference between the magnitude（a）and the occurrence time（b）of the same event in the catalog report and the detection results

3 討論

3.1 雙震事件檢測

利用上述配置對河源M4.3 的余震自動檢測結(jié)果顯示，多個雙震事件沒能被掃描識別出來。經(jīng)過對地震目錄的分析，發(fā)現(xiàn)這幾次雙震事件發(fā)生的時間間隔很短，普遍在2 s 以內(nèi)。由于系統(tǒng)的配置未對地震事件進行篩選，全部處理為地震模板并應用于模板匹配，導致部分單個地震模板內(nèi)包含了多次地震事件。在對波形進行掃描時，直接將其與連續(xù)波形進行滑動互相關(guān)并輸出匹配結(jié)果，當雙震的震級差比較大時，輸出為掃描時間窗內(nèi)震級較大的地震事件。

在利用M&L 方法進行匹配定位時，一般需要對地震模板進行篩選，尤其是多個地震事件的發(fā)震時間間隔較近，前一事件的s 震相淹沒在后一事件更大的p 震相中的情況。盡管在處理地震模板時，可以控制地震模板的波形長度，在一定程度上消除另一事件的影響，但這也可能導致其他地震模板的長度過短，不能包含所有必需的震相信息；另外，雙震事件作為模板也將導致檢測出的事件的定位結(jié)果和震級產(chǎn)生較大的誤差，雙震以及多地震事件不適合作為地震模板進行匹配定位檢測。在對地震模板庫進行自動更新，需結(jié)合地震目錄，剔除使用發(fā)震間隔較近、且發(fā)震位置也較近的地震事件。

3.2 臺網(wǎng)編目工作中的適用性分析

本文與前人[15]在新豐江河源地區(qū)的微震檢測中，對于平均互相關(guān)系數(shù)的參數(shù)設定存在差異，這種差異對微震的檢測結(jié)果數(shù)量產(chǎn)生影響。在不同的研究區(qū)域，需要對研究區(qū)域內(nèi)擬使用的地震模板進行自檢來估計平均背景相關(guān)系數(shù)，以此確定檢測標準。數(shù)字臺網(wǎng)中心利用M&L 方法對不同區(qū)域的微震進行檢測時，需要重新確定檢測標準以保證檢測質(zhì)量。

M&L 方法在地震模板處理中，利用TauP 工具包[23]根據(jù)給定地殼速度模型標記理論震相信息截取地震模板。通過計算模板事件與待檢測事件的可能位置在同一臺站上的走時差來進行定位計算，對模型的依賴較小。本文利用PREM 模型生成的地震模板，發(fā)現(xiàn)其震相標注與實際震相的到時差異較大，這對于檢測微震以及對其進行定位的結(jié)果沒有影響。但是，地震目錄包含精確的震相標記信息，可以利用這些震相信息替代理論到時標記，對地震模板進行截取、標記處理，在匹配定位檢測的同時，輸出比較準確的震相信息，更加適合地震編目的完善工作。

另外，M&L 程序已經(jīng)發(fā)布GPU（Graphics Processing Unit-Based）版本，該程序的計算速度提高到M&L 程序的4.5倍，同時在地震模板的各分量波形上引入權(quán)重系數(shù)顯著提高檢測結(jié)果的穩(wěn)定性[13]。這可以在一定程度上緩解由于地震模板過多導致的檢測效率過低的問題，適合臺網(wǎng)處理大量微震檢測工作。

4 結(jié)論

本文利用匹配定位M&L 方法，通過對匹配定位程序與數(shù)據(jù)處理系統(tǒng)軟件進行配置和參數(shù)設置，在新豐江河源地區(qū)建立地震模板數(shù)據(jù)庫，進行微震事件自動檢測，并對河源M4.3 地震震后余震的自動檢測結(jié)果進行分析，得到以下結(jié)論：

（1）在新豐江地區(qū)，存在有大量未被記錄的微震事件，匹配定位方法可以準確地對這些事件進行檢測以及定位；

（2）利用數(shù)據(jù)處理系統(tǒng)軟件將數(shù)據(jù)庫內(nèi)的地震事件自動轉(zhuǎn)換為地震模板，建立新豐江地區(qū)的地震模板庫并對其進行實時擴充，方便在微震檢測中進行調(diào)用；

（3）根據(jù)對匹配定位程序的配置和參數(shù)設置，檢測出河源M4.3 地震震后大量ML＜0.0 的地震事件，可以達到完善地震目錄的要求。

致謝：感謝張淼博士提供地震模板匹配程序（Match and Locate），本文部分圖件使用GMT繪制。