匹配濾波方法研究2016年意大利中部MW 6.0地震的余震遷移

付裕 黃暉 徐鳴潔

南京大學地球科學與工程學院，南京市棲霞區仙林大道163號 210046

0 引言

一般情況下，淺部大地震之后會產生相當多的余震，隨著時間延續，余震的發生率呈現冪次性衰減（Omori，1894；Gutenberg et al，1941；Aki，1965），在空間上，余震由主震破裂區域向外遷移，而余震的時空演化及觸發機制對研究地震斷裂活動性和評估地震災害至關重要。地震的觸發機制可能有多種，如靜態應力的變化（Stein，1999；Meng et al，2013）、慢滑移（Hearn et al，2002；Lohman et al，2007；Peng et al，2009；Borghi et al，2016）、動態應力觸發機制（Brodsky et al，2000）及殼內流體（Shapiro et al，2006）等。大地震之后，會產生數量巨大的余震事件，但大地震的尾波信號往往掩蓋了震級較小的地震信號，而且很多小地震的震相重合，這使得常規方法得到的地震目錄中缺失很多余震事件（Peng et al，2009）。因此，提取出更多的余震事件，獲得較完整的地震目錄，對分析大地震后早期余震分布和遷移規律至關重要。

檢測地震事件的常用方法有 STA／LTA（長短時間窗平均值比）、MFT方法。STA／LTA（Stevenson，1976）方法主要基于長時間窗信號平均值與短時間窗信號平均值的比值來判斷地震事件。但這種方法只能檢測出高信噪比事件，仍會遺漏很多早期余震。MFT（the matched filter technique，模板匹配濾波）（Gibbons et al，2006）方法通過模板地震事件在連續波形上滑動求取互相關值來檢測地震事件。此方法對信噪比要求較低，更依賴于事件波形之間的相關性，相比STA／LTA方法，能夠行之有效地檢測出微小的地震事件，模板匹配濾波技術在微地震檢測方面具有明顯優勢（李璐等，2017）。基于該性質，MFT方法不僅可用來檢測天然地震的前震和余震事件（Peng et al，2009；Kato et al，2012；Sugan et al，2014；譚毅培等，2014；Wu et al，2017；Huang et al，2017；侯金欣等，2017），還可被用于非火山顫動和低頻地震（Shelly et al，2007）、油氣勘探中水壓致裂誘發地震（Bao et al，2016）、爆炸監測（Zhao et al，2016）等研究領域。

2016年8月24日1時36分（國際標準時），意大利中部亞平寧山脈附近發生MW6.0地震，震源深度約8km，地震對當地居民和建筑造成了較大的傷害。2個月之后，2016年10月26日在MW6.0地震震中的北部，2h內相繼發生了MW5.5、MW5.9地震。緊接著4天以后，在MW6.0、MW5.9地震震中的中間區域又發生了 MW6.5地震，形成了序列地震。在過去20年里，該區域曾發生2次MW≥6.0地震（即1997年Colfiorito MW6.0、2009年L'Aquila MW6.3地震）（圖1（a））。

本文利用意大利中部密集臺網的地震記錄，采用模板匹配濾波方法檢測2016年8月24日MW6.0地震發生后80天內的余震事件，研究此地震序列的余震時空分布規律，探索余震遷移和地震觸發機制。

1 地質背景

亞平寧山脈中部地區位于意大利中部，處于阿德里亞微板塊和歐洲板塊的擠壓和俯沖作用帶上（圖1（b）），是意大利境內活躍的地震活動帶之一。前人對該研究區的地質研究顯示，亞平寧中部地區現今主要的地震活動是由一系列NW-SE走向的正斷層構成的斷裂系引起的，該斷裂系形成的地質時期為第四紀（Boncio et al，2004；Galli et al，2008；Pizzi et al，2009；Lavecchia et al，2016），大多數都是傾向SW，以雁列式分布在山脈附近。主要的斷裂有Colforito斷裂（CF）、Bove Faul斷裂（BF）、Vettore斷裂（VF）、Norcia斷裂（NF）和 Gorzano斷裂（GF）等（圖1（a））。該斷裂系又被 2組先存構造帶分割（圖1（a）中的 OAST和 MCT），這2組逆沖斷裂帶形成于第四紀之前的擠壓構造活動（Pizzi et al，2009）。

亞平寧中部地區以中強度的正斷層地震活動為主（Chiarabba et al，2005），圖1（a）中幾次大地震的震源機制解均顯示了拉張為主的活動性質（Chiaraluce et al，2005；http：／／cnt.rm.ingv.it），主張應力軸呈NE-SW向（圖1（a））。研究區的大地測量結果顯示為向山脈兩側擴張，形變速率約3mm／a，方向主要為 NE-SW 向（Serpelloni et al，2005；Petricca et al，2015），與震源機制顯示的拉張主應力軸方向一致。

圖1 意大利中部地質構造及地震臺站、模板事件分布

此外，研究區域內的地震活動常常以序列地震的形式發生。1997年、2009年序列地震都包含了數個 5級以上的地震（Amato et al，1998；Chiarabba et al，2009；Chiaraluce et al，2011），1997年主要的地震分布在2016年地震的NNW方向，2009年L'Aquila地震分布在SSE方向，且震源機制解均顯示以拉張型的斷裂活動為主（Ekstr?m et al，1998；Chiaraluce et al，2011）。2016年8月底至10月底的地震活動主要發生在1997年、2009年地震空區處（Liu et al，2017；Chiaraluce et al，2017），余震主要分布在BF、VF與NF兩組平行的斷裂帶附近。

2 數據與方法

2.1 地震數據來源

本文使用的地震數據是由歐洲臺網中心（http：／／www.seismicportal.eu）獲取的三分量寬頻帶連續波形數據，分別選自意大利CAMP和 AQU等13個臺站（圖1（c）），連續波形采樣率為100Hz，記錄時間從2016年8月24日0時～11月12日0時（UTC時間），總共80天。

地震目錄來源于國際臺網中心 ISC（http：／／www.isc.ac.uk），我們共獲取2016年8月24日～11月1日發生的1426個地震事件作為后續匹配濾波掃描的模板。

2.2 地震模板截取

在模板地震事件截取之前，首先對地震數據進行了預處理，為了削弱遠震事件的影響，采取濾波范圍1～8Hz，連續波形數據重采樣至25Hz，模板數據從預處理后的連續波形數據中截取。

對于微小地震來說，由于S波能量強于P波，更容易從噪聲中檢測出來，因此模板選自S波段。根據地震目錄，選取S波理論到時的前2s、后6s為時間窗截取模板信號，選用P波到時前12～4s作為背景噪聲信號計算模板的信噪比，取信噪比大于5作為初選模板。為保證最終模板的可信度，需足夠多的模板數據用于疊加，因此要求選取的每個模板數據有4個以上臺站的記錄（即至少12個分量的記錄數據）。其中，P、S波理論到時由 TAUP軟件（Crotwell et al，1999）計算得到，速度模型采用當地 CIA模型（Herrmann et al，2011）。最終，選出了1381個模板事件（圖1（c））。

2.3 互相關計算

第1步，從每天記錄的零時開始，計算某一分量模板與相同分量、時長連續波形數據的互相關系數（cross-correlation coefficient，CC），然后將模板在連續波形上后移1個采樣點，每移動1次計算1個互相關系數，對1天的記錄得到該分量的互相關系數序列。第2步，將不同臺站同一模板事件的S波按照到時對齊，對不同臺站所有分量計算平均互相關值（mean CC）。第3步，假定每一天的信號序列是符合高斯分布的隨機過程，信號序列的標準差為1.4826倍絕對中位差（median absolute deviation，MAD），相關值大于12倍 MAD（8.1倍標準偏差）的概率為3.3×10－16。由于1天記錄的數據點為2 160 000個，選取1天的平均互相關值組成序列的12倍MAD作為閾值，每天的錯誤檢測率小于1個。當序列上某個時刻的平均互相關系數高于閾值，則認定是檢測出來的事件（圖2）。

2.4 檢測事件的確認

對于同一時間窗內出現多個模板匹配的情況，保留平均CC值最高的模板。對于新檢測出的事件，由于事件波形和模板事件相似度最高，其震源位置選用模板事件的位置；新事件的發震時刻由CC值最高的時刻相對模板S波到時的差值，再加上模板的發震時刻計算得到；檢測事件的震級，通過與模板在所有分量上最大絕對振幅的比值的中位值計算得到，依據振幅每增加10倍，震級增加1級（Peng et al，2009）。

3 結果與分析

3.1 地震目錄完備性分析

利用選取的1381個模板事件，從2016年8月24日～11月12日總共80天的連續波形中新檢測出71525個余震（12倍MAD），數量是選用模板的數十倍。新檢測出的余震震級范圍為－0.8～5.3級，與模板事件的震級范圍（0.8～6.5級）相比，檢測的余震事件震級更小（圖3）。對新檢測的余震與模板事件作地震頻度與震級大小間關系的分析，采取W iemer等（2000）提出的GOF（goodness of fit）方法，當曲線擬合度達到90%時，得到完備震級Mc，計算結果如圖4所示。由檢測結果和模板數據計算得到的 Mc約為1.0，相比模板目錄的 Mc（2.5），有了大幅度的降低。檢測的余震事件中震級大于1.0的地震事件共有29035個，加上1381個模板事件，組成了新的事件目錄，新目錄共包含30416個地震事件，由此提高了地震目錄的完備性。我們利用新的地震目錄，分析意大利中部地區2016年8月24日～11月12日余震的時空分布規律。

圖2 MFT檢測結果實例

圖3 地震事件統計

3.2 余震分布的總體特征

對于本文研究中2016年意大利中部發生的序列地震，依據其發震時刻和所在位置，將其分為 Event 1（2016年8月24日MW6.0地震）、Event 2（2016年10月26日MW5.5、MW5.9地震）以及Event 3（2016年10月30日MW6.5地震），以進一步分析序列地震中不同事件余震的時空演化過程與差異。

圖4 檢測結果與模板的完備震級Mc（a）及擬合度曲線（b）

圖5為2016年8月24日～11月12日余震分布，縱軸為沿測線A-A′（圖1（c））的距離，圖5一定程度上反映了該序列地震的余震分布規律。

圖5 余震總體特征分布

為了更清晰地分析余震的分布特征，按照4個主事件位置沿測線分成4段進行余震的發震頻度統計（即每小時發生地震的個數）（圖5中彩色曲線）：橙色線表示Event 1（測線0位置）的東南段、投影距離為－20～0km（測線 A′至 Event1之間）的余震頻度曲線；黃色線表示沿測線0～17km（Event1至 Event3之間）的頻度曲線；紅色線為沿測線17～25km（Event3至Event2-2之間）的頻度曲線；綠色線顯示25km西北側的頻度曲線。由圖5可見，Event 1的余震事件主要發生在測線－15～20km的位置，主震東南側較快地趨于平靜（圖5橙色頻度曲線），而西北側的余震遷移至Event 2震中附近，并且持續到Event 2之前仍然較活躍（圖5藍色點分布）；Event 2（包括Event 2-1和Event 2-2）的余震事件主要分布在測線15～45km之間，余震活動以Event 2-1西北段的為主（圖5紅色、綠色頻度曲線）；Event 3的余震事件分布范圍較廣，大致為測線－10～45km的范圍內，其東南側的活動頻率明顯高于西北側（圖5黃色、紅色頻度曲線）。

3.3 早期余震的空間分布

為了進一步分析序列地震中幾次主要地震的早期余震空間分布特征，分別從新余震事件目錄中提取Event 1、Event 2、Event 3發生后3天內的余震，繪制平面分布圖（圖6）。由圖6可見，Event 1發生后104s（～3h）內余震（圖6（a）藍色點）主要集中在 NF斷裂東側，并沿斷裂延伸方向朝主震的兩側擴散，呈現沿斷層走向向南北兩側對稱發生遷移的特征；在104s之后的時間內（圖6（a）紅色點），余震向南擴散的范圍未見明顯變化，而向北的范圍明顯擴大；余震活動主要集中分布在NF、GF、VF斷裂之間（圖6（a））。

圖6 余震空間分布

Event，2的早期余震分布特征（圖6（b））表現為：Event，2-1發生后約 103s（～0.5h）內余震主要集中在震中附近，隨后余震沿斷裂走向逐漸向北遷移至 Event，2-2附近，2h后發生MW5.9地震；在104s（～3h）左右逐漸向東北部轉移到 BF斷裂附近（圖6（b）藍色點），此時向南遷移的余震較少；相反，104s之后（圖6（b）黃、紅色點）余震活動呈現向南遷移的特征，NF、VF斷裂之間逐漸被余震覆蓋，但基本上沒有超過OAST斷裂，表明受到先存斷裂OAST的限制（Pizzi et al，2009）。Event 2的余震活動向南遷移的過程明顯較向北的遷移更晚，呈現“南北不對稱”的余震遷移規律。

Event，3的余震分布范圍更廣。在104s（～3h）以內，余震主要沿著斷裂走向朝南遷移，104s之后，余震在斷裂北端活動加強（圖6（c）紅色點），余震遷移表現出“南北不對稱”的規律。Event，3的大部分余震在平面上分布于NF、BF、VF斷裂附近，僅少部分分布在GF斷裂附近，總體上以分布在OAST斷裂北部為主，西北部基本上沒有超過MCT斷裂，表明先存的MCT逆沖斷裂對VF、BF斷裂的活動起到一定限制作用（Pizzi et al，2009）。

3.4 余震活動性與遷移速率

我們通過沿A-A′方向余震隨對數時間的分布來估算余震向周圍遷移的速率（圖7），以進一步分析余震的遷移活動規律。鑒于MW6.0地震的余震在2個月后已經相當微弱（圖5），我們把 Event 2-1之后時段分離出來分析 Event 2、Event 3余震活動特點與差異（圖7（b））。

圖7 沿A-A′方向余震隨對數時間的分布

在8月24日MW6.0地震之后1天內，余震遷移前端沿斷裂帶向南、北以2～3km／decade速率擴散到15km附近，之后向南部的擴展逐漸減弱，相比之下，北端余震仍不斷擴展，約1周后出現遷移速率明顯加快的現象（圖7（a））。以沿投影距離每2km余震數目達30個以上來計算余震活動曲線（Kato et al，2014）（圖7紅色虛線）可更細致反映余震在遷移過程中的表現，并顯示不同時段、不同位置余震的活躍程度。由圖7可見，在Event 2、Event 3發震前出現了幾處遷移速率明顯加快的現象，余震活動曲線呈階梯狀。Event1發生后的1周內，余震活動較平穩地向北、向南擴散；1周后，北部活動曲線明顯變化，遷移速率加大，顯示向北的活動出現了加快的現象，擴散到Event 2-1震中附近，而向南的余震活動在測線－15 km附近不再向南擴展（圖7（a））。在Event 2-1發生后1～2 h之間，余震活動曲線顯示向北的遷移速率有小幅度加快，向Even 2-2震中遷移，2h后發生了 Even 2-2。Even 2-2發生后，余震活動曲線上向南遷移的余震出現遷移略有加快的現象，達到Event 3震中附近。這可能反映了MW6.0地震的余震遷移在Event 2、Event 3震中附近造成了應力加載，促使了其發生。

3.5 對余震成因機制的討論

余震的觸發、遷移機制可能受多種因素共同影響，如靜態應力、動態應力觸發、慢滑動和流體遷移等。Liu等（2017）、Convertito等（2017）通過計算靜態庫侖應力變化認為，2016年意大利中部MW6.0地震的余震活動主要由同震應力加載所致；Vuan等（2017）擬合余震遷移前端，計算不同深度地殼流體擴散因子D（Shapiro et al，2006）時發現，MW6.0地震的余震活動反映的D值（淺部D值為25～50m2／s，深部D值為160m2／s）過高，排除了粘彈性松弛原因，并認為可能是慢滑動所致。我們通過擬合余震遷移前端發現，MW6.0地震早期余震的兩側遷移較符合lg t的特征，表明可能與慢滑動的遷移有關。除此以外，MW5.5地震（2016年10月26日17時）之后朝南側的遷移也符合lg t特征，說明早期余震遷移活動與慢滑動的遷移有關。我們推測慢滑動遷移可能引起了相鄰區域的應力加載，對后續MW5.9、MW6.5地震的發生起到一定程度的觸發作用。但要嚴格論證慢滑動造成的應力加載，還需進一步的研究工作。

4 結論

我們利用模板匹配濾波方法（MFT）有效地檢測出意大利中部2016年8月24日序列地震的余震，獲得了數十倍于原模板數量的地震檢測結果，得到檢測余震的完備震級Mc由2.5降至1.0，提高了地震目錄的完備性。

對余震時空分布特征的分析顯示，序列地震中不同事件余震的擴散方式和遷移速率存在一定差異。MW6.0地震的早期余震分布呈現向南北兩側對稱擴散的特征，遷移速率2～3km／decade。MW5.5～5.9、MW6.5地震的早期余震分布表現出“南北不對稱”的余震遷移現象：MW5.5～5.9地震的早期余震（3h內）主要向北遷移，然后再向南部擴散；而MW6.5地震早期余震表現出先向南、后向北遷移的現象，且整體的活動范圍更廣。

MW6.0地震向北遷移的余震，在主震發生1周后遷移加快，可能反映了MW6.0地震的余震活動造成后續大地震震中附近的應力加載。MW6.0地震向兩側遷移的早期余震前端和10月26日MW5.5地震向南遷移的早期余震前端符合lg t的特征，說明早期余震遷移活動與慢滑動的遷移有關。