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

付裕 黃暉 徐鳴潔

南京大學(xué)地球科學(xué)與工程學(xué)院，南京市棲霞區(qū)仙林大道163號 210046

0 引言

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

檢測地震事件的常用方法有 STA／LTA（長短時間窗平均值比）、MFT方法。STA／LTA（Stevenson，1976）方法主要基于長時間窗信號平均值與短時間窗信號平均值的比值來判斷地震事件。但這種方法只能檢測出高信噪比事件，仍會遺漏很多早期余震。MFT（the matched filter technique，模板匹配濾波）（Gibbons et al，2006）方法通過模板地震事件在連續(xù)波形上滑動求取互相關(guān)值來檢測地震事件。此方法對信噪比要求較低，更依賴于事件波形之間的相關(guān)性，相比STA／LTA方法，能夠行之有效地檢測出微小的地震事件，模板匹配濾波技術(shù)在微地震檢測方面具有明顯優(yōu)勢（李璐等，2017）。基于該性質(zhì)，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）、油氣勘探中水壓致裂誘發(fā)地震（Bao et al，2016）、爆炸監(jiān)測（Zhao et al，2016）等研究領(lǐng)域。

2016年8月24日1時36分（國際標(biāo)準(zhǔn)時），意大利中部亞平寧山脈附近發(fā)生MW6.0地震，震源深度約8km，地震對當(dāng)?shù)鼐用窈徒ㄖ斐闪溯^大的傷害。2個月之后，2016年10月26日在MW6.0地震震中的北部，2h內(nèi)相繼發(fā)生了MW5.5、MW5.9地震。緊接著4天以后，在MW6.0、MW5.9地震震中的中間區(qū)域又發(fā)生了 MW6.5地震，形成了序列地震。在過去20年里，該區(qū)域曾發(fā)生2次MW≥6.0地震（即1997年Colfiorito MW6.0、2009年L'Aquila MW6.3地震）（圖1（a））。

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

1 地質(zhì)背景

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

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

圖1 意大利中部地質(zhì)構(gòu)造及地震臺站、模板事件分布

此外，研究區(qū)域內(nèi)的地震活動常常以序列地震的形式發(fā)生。1997年、2009年序列地震都包含了數(shù)個 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月底的地震活動主要發(fā)生在1997年、2009年地震空區(qū)處（Liu et al，2017；Chiaraluce et al，2017），余震主要分布在BF、VF與NF兩組平行的斷裂帶附近。

2 數(shù)據(jù)與方法

2.1 地震數(shù)據(jù)來源

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

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

2.2 地震模板截取

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

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

2.3 互相關(guān)計算

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

2.4 檢測事件的確認(rèn)

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

3 結(jié)果與分析

3.1 地震目錄完備性分析

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

圖2 MFT檢測結(jié)果實例

圖3 地震事件統(tǒng)計

3.2 余震分布的總體特征

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

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

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

圖5 余震總體特征分布

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

3.3 早期余震的空間分布

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

圖6 余震空間分布

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

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

3.4 余震活動性與遷移速率

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

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

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

3.5 對余震成因機制的討論

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

4 結(jié)論

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

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

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