王鵬 侯金欣 吳朋

1）山東省地震局，濟南市歷城區(qū)港西路2066號 250102

2）中國地震局地球物理研究所，北京市海淀區(qū)民族大學南路5號 100081

3）四川省地震局，成都 610041

0 引言

2017年8月8日21時19分，四川省阿壩州九寨溝縣發(fā)生MS7.0地震，震源深度20km。此次地震是繼2008年汶川MS8.0地震和2013年蘆山MS7.0地震之后在青藏高原東部邊緣發(fā)生的又一次破壞性地震，地震波及到四川、甘肅、青海、寧夏、陜西等省份，共造成24人死亡，近500人受傷（郭安寧等，2017）。地震發(fā)生在東昆侖斷裂帶東端、構造轉換地區(qū)、呈帚狀分叉的次級走滑樹正斷裂帶上，向南東為雄黃山隆起區(qū)，其邊界斷裂為岷江斷裂、塔藏斷裂、雪山斷裂和虎牙斷裂等（袁道陽等，2016）。此次地震為走滑型地震（USGS），發(fā)震斷層——樹正斷裂位于岷江斷裂帶與塔藏斷裂中間，地震時斷層向兩端破裂擴展，北段向塔藏斷裂方向擴展，南段與虎牙斷裂貫通。該區(qū)歷史上地震活動較為活躍，震中附近200km范圍內(nèi)共發(fā)生26次5級以上地震（高曙德等，2017）。

當一次強震發(fā)生后，人們最關心的問題莫過于震后是否還會有強余震或者更大的破壞性地震發(fā)生。對余震序列衰減特征以及早期序列參數(shù)特征的研究一直是地震學和地球物理學研究的熱點問題（蔣海昆等，2007）。余震序列衰減特征研究是震后趨勢判定、強余震預測和發(fā)震構造分析等研究的重要組成部分，對于理解主震發(fā)震過程、分析地震危險性和抗震救災等工作具有重要作用（宋金等，2009）。余震序列的衰減特征和激發(fā)余震的能力可通過序列的統(tǒng)計參數(shù)來表征，由于震源區(qū)的孕震條件和物理狀態(tài)的不同，不同的構造區(qū)域、序列類型、主震的震源機制類型、震源區(qū)局部構造應力水平及大地熱流等地球物理特征都會有不同的序列參數(shù)特征（Kagan et al，2010）。因此，深入分析強震的余震序列活動特征具有重要的意義。

余震序列活動特征主要包括序列時間演化、空間展布及強度分布特征等3個方面。其中，序列時間演化，一方面可基于b值、h值、k值等序列活動參數(shù)判定早期發(fā)展趨勢（劉正榮等，1986；劉正榮，1984；朱傳鎮(zhèn)等，1989；吳開統(tǒng)等，1996；馬玉虎等，2012；宋春燕等，2014）；另一方面可采用具有一定物理背景的統(tǒng)計模型，如修正的大森公式（Utzu，1961、1965；譚毅培等，2015a、2015b）、ETAS模型（Ogata，1988、1989、2001；蔣海昆等，2007；蔣長勝等，2013）以及 BASS模型（Sheherbakov et al，2004；Turcotte et al，2007）等來判定。

地震目錄是研究余震序列衰減的基礎資料，目錄的完整性直接關系著研究結果的科學性和可信性。強震發(fā)生后短時間內(nèi)往往伴隨著大量的余震事件，其地震波形受主震和強余震的尾波干擾，波形相互疊加，一些地震事件的震相難以清晰識別，因此，測震臺網(wǎng)的地震目錄往往會遺漏大量的地震事件。研究顯示，主震后短時間內(nèi)遺漏的余震可達到地震目錄給出余震數(shù)量的數(shù)倍甚至十幾倍（Peng et al，2007、2009；Lenglinéet al，2012；譚毅培等，2015c），如此大量的遺漏地震勢必會影響余震序列地震目錄的完整性，進而影響以地震目錄為基礎的地震活動性、發(fā)震構造和發(fā)震機理等研究的可靠性。

基于波形互相關的模板匹配技術近期發(fā)展迅速且應用廣泛。Gibbons等（2006）證實了模板匹配方法在微震檢測方面更有優(yōu)勢；Shelly等（2007）研究了非火山地脈動與低頻地震間的關系；Schaff（2008）及Schaff等（2010）進一步證實了在地震監(jiān)測中模板匹配方法的識別能力；Peng等（2009）完備了2004年 Parkfield地震的余震目錄并分析了余震的遷移特征；Yang等（2009）確認了2008年 Illinois Monut CarmelMW5.8地震的斷層面；Meng等（2013）發(fā)現(xiàn)了2003年San SimeonM6.5地震之后Parkfield地區(qū)的地震活動性與靜剪切應力間的對應關系；譚毅培等（2016）研究了灤縣震群的發(fā)震構造；馬騰飛（2015）研究了2008年汶川地震序列中的“重復地震”；Zhang等（2015）實現(xiàn)了對朝鮮地區(qū)多次核爆位置的精確檢測。雖然模板匹配技術已成功應用在地震學研究的多個領域，但由于其計算量過于龐大，傳統(tǒng)的CPU計算會花費大量的時間，Meng等（2012）利用GPU并行計算實現(xiàn)了模板匹配技術算法加速，并完備了2010年Mayor-CucapahMW7.2地震之后Salton Sea地區(qū)的地震目錄。

本文利用基于GPU加速的模板匹配技術來檢測九寨溝MS7.0地震前后的遺漏地震，再基于補充遺漏地震后的完備地震目錄，利用b值和Mc值對地震序列的時間演化特征進行分析。

1 遺漏地震檢測

強震序列震后短時間內(nèi)波形相互交疊，不易標識，容易造成大量余震的缺失。研究表明，西南地區(qū)的強震序列在震后5天以內(nèi)衰減的能量最為顯著，此后，激發(fā)次級余震的能力明顯減弱（郭路杰等，2014）。2017年8月8日九寨溝MS7.0地震發(fā)生之前，僅在8月7日發(fā)生1次1.5級地震，前震序列偏少。因此，為了考察完整的地震序列的活動特征，本研究選取了主震前7天和主震后5天的連續(xù)波形進行分析計算，所用數(shù)據(jù)為地震序列周邊的九寨溝臺（JZG）、若爾蓋臺（REG）、松潘臺（SPA）、平武臺（PWU）、舟曲臺（ZHQ）和文縣臺（WXT）等6個固定臺站記錄的三分量連續(xù)波形（圖1），數(shù)據(jù)采樣率為100Hz。

圖1 臺網(wǎng)目錄地震和臺站分布

模板匹配技術是在滑動窗互相關（sliding-window cross-correlation，SCC）檢測技術的基礎上發(fā)展起來的（Yang et al，2009），是在低信噪比和稀疏臺網(wǎng)情況下提取信號的一種有效方法。模板匹配技術是以已知的震相清晰的地震事件作為模板，利用波形相關性，掃描連續(xù)波形數(shù)據(jù)，尋找與模板事件波形相似度高的一段連續(xù)波形作為新檢測事件，以實現(xiàn)對遺漏地震的自動識別和檢測。

利用研究時段內(nèi)臺網(wǎng)記錄到的1336個ML≥0地震事件作為候選模板，其中，最大震級為MS7.0。對數(shù)據(jù)進行去均值、去線性趨勢等預處理后，采用2～12Hz的Butterworth帶通濾波器濾波；利用震相報告中的P、S波到時計算模板事件的信噪比，舍棄信噪比小于5的模板事件，最終選取了滿足條件的1033個地震事件作為模板。利用P波、S波到時的前1～5s作為掃描窗長，分別與被檢測連續(xù)波形的垂直分量和水平分量作互相關計算

其中，、y分別為x（t）、y（t）的均值；t0、t1為掃描時窗的開始、結束時間；CC為互相關系數(shù)（cross-correlation coefficient，CC）。

將各臺站垂直、水平分量的CC值所對應的時間分別減去模板對應的P、S波到時與連續(xù)波形起始時間之差，將各臺站、各分量的互相關系數(shù)序列進行偏移后的疊加處理，最后，除以所用分量的有效個數(shù)來計算平均CC值，得到平均互相關系數(shù)波列。同時，計算CC值波列的絕對離差中位數(shù)（median absolute deviation，MAD），若平均CC值大于其閾值，則判定為1次地震事件。本文采用12倍MAD作為閾值（Peng et al，2009；譚毅培等，2014；侯金欣等，2017），考慮到采樣率為100Hz時，0.01s的移動掃描間距會造成大量的誤檢測事件，因此，取2s內(nèi)相關系數(shù)最大的地震為檢測事件。

基于模板匹配方法的假設，檢測事件與模板事件具有相同的走時，由此得到檢測事件的時間。檢測事件的震級可以通過檢測事件與其對應的模板事件在各臺站的水平分量上S波到時前2s至后2s內(nèi)的最大振幅比的平均值求得（圖2）。

圖2 模板匹配方法檢測遺漏地震示例

通過計算，共檢測到九寨溝地震序列2017年8月1～12日間的5887個地震事件，其中，包括1033個模板事件的自檢測，因此，共4854個檢測事件。中國地震臺網(wǎng)中心的地震目錄共記錄到地震事件3057個，其中，可定位地震事件1460個，單臺事件1597個。由此可見，本文的檢測事件個數(shù)為臺網(wǎng)可定位地震個數(shù)的3倍多，共檢測到1797個遺漏地震事件，該數(shù)量遠大于可定位地震和單臺事件的數(shù)量。單臺事件為震后僅在震中距最近的1個臺站能檢測到P、S波震相的地震事件，多出現(xiàn)在主震發(fā)生后較短的時間內(nèi)，是受強震后續(xù)震相干擾、余震波形相互重疊等因素的影響，因記錄的余震波形信噪比較低而造成的。此次地震序列的單臺事件數(shù)略大于可定位事件數(shù)，約占臺網(wǎng)目錄的52.2%，與譚毅培對岷縣漳縣MS6.6地震序列主震后1h內(nèi)單臺事件所占的比例一致（譚毅培等，2015b）。單臺事件震中缺失且震級可靠性較差，而檢測事件是基于多個臺站檢測到的地震事件，對震級的估計更為準確，因此，基于補充了遺漏地震的地震目錄對地震序列的活動特征進行分析具有更強的科學性和可靠性。

將檢測到的遺漏地震補充到臺網(wǎng)目錄中（臺網(wǎng)目錄截至2017年9月30日），得到更加完整的地震序列目錄，這使得8月8日主震之前的前震序列更加豐富完整，而不僅是臺網(wǎng)目錄中記錄的1次地震（圖3）。從前震序列的演化可見，8月1日以來存在著地震活動的平靜—增強過程，特別是主震前地震活動短時間內(nèi)的增強，可能是應力的積累導致震源區(qū)發(fā)生微小破裂所致，這一階段可能是深部更大的斷層發(fā)生破裂之前加速蠕動的臨震階段（牛志仁，1978）。同時，主震后5天內(nèi)的余震序列目錄也更加完整，主要體現(xiàn)在0.0～1.0級震級檔的補充，遺漏地震事件大多集中在這一震級區(qū)間內(nèi)（圖3），而從震后5天以后的臺網(wǎng)目錄記錄地震的分布可見，0.0～1.0級震級區(qū)間的地震較震后5天內(nèi)的數(shù)量明顯增多，臺網(wǎng)目錄的完整性較好，因此，僅對震后5天作了遺漏地震檢測。從補充后的整個地震序列來看，震級主要集中在1.0～2.0級區(qū)間內(nèi)，地震頻次的統(tǒng)計也呈現(xiàn)較好的正態(tài)分布。

圖3 地震序列震級-時間（a）和震級-頻度（b）分布

2 地震目錄完備性分析

地震目錄是進行地震活動性分析、活動構造研究、地震預測與地震危險性評估的重要基礎資料（徐偉進等，2014），而地震目錄的質量則依賴于研究區(qū)域地震臺網(wǎng)的檢測能力。完備震級Mc是國際上普遍采用的評估地震臺網(wǎng)檢測能力的定量標準，通常被定義為在一個時空范圍內(nèi)，地震能被臺網(wǎng)100%監(jiān)測到的最小震級（Rydelek et al，1989）。假設地震的發(fā)生是自相似過程，則觀測到的地震累積發(fā)生頻率應滿足G-R定律。在震級較大處，地震發(fā)生頻率偏離線性G-R定律，其原因可能是因大地震數(shù)目太少而引起的隨機波動或特征地震現(xiàn)象（黃亦磊等，2016）。而小震級處的偏離，則被解釋為臺網(wǎng)監(jiān)測能力的不足造成震級小于完備震級的地震事件可能未被臺網(wǎng)監(jiān)測而在地震臺網(wǎng)目錄中被遺漏，導致地震目錄的不完備，從而影響到基于地震目錄的相關研究結果的可靠性（李智超等，2014）。而通過檢測遺漏地震來補充地震臺網(wǎng)目錄，則會使得地震目錄更加完備（Schaff，2008），更有利于提高分析結果的科學性和可靠性。

本文利用最大曲率法 MAXC（The Maximum Curvature technique）（Wiemer et al，2000）和擬合優(yōu)度測試法 GFT（The Goodness of Fit Test）（Wiemer et al，2000）分析地震序列的完備震級，其中，最大曲率法選取震級-頻率曲線中斜率最大值所對應的震級作為Mc。在實際中，這個震級往往對應非累積震級-頻率分布中擁有最多地震數(shù)目的震級。而擬合優(yōu)度法則是通過計算不同假定起始震級與 G-R關系間的擬合情況來確定起始震級的方法。

圖4 2017年8月1～12日2種不同目錄頻次及b值對比

通過最大曲率法對8月1～12日的檢測目錄和臺網(wǎng)目錄進行分析（圖4）。由圖4可見，檢測目錄的頻次明顯大于臺網(wǎng)目錄，但非累計頻次與震級的關系均顯示2種目錄在1.2級處記錄的地震頻次最高，因此，完備震級都是1.2級。利用最大似然法對累計頻度與震級的擬合發(fā)現(xiàn)，完備震級也相差不大，檢測目錄為1.2級，臺網(wǎng)目錄為1.3級。檢測目錄擬合得到的b值為0.8558，略高于臺網(wǎng)目錄的0.7682（圖4）。由于擬合優(yōu)度法對于每一個震級下限都采用最似然估計計算，因此，在處理臺網(wǎng)探測地震能力隨震級變化較快和有不均勻性的地震目錄時更加有優(yōu)勢。在擬合優(yōu)度測試法的結果中，檢測目錄的完備震級為1.4級，小于臺網(wǎng)目錄的1.6級（圖5）。研究表明，主震后短時間內(nèi)記錄的單臺事件給出的震級明顯偏低（譚毅培等，2005b）。2種目錄的完備震級相差不大，可能是因為大量單臺事件的加入使得臺網(wǎng)目錄在小震級區(qū)間內(nèi)的地震數(shù)增多，使得完備震級被低估。

圖5 擬合優(yōu)度測試法計算的2種目錄的完備震級

3 序列活動特征分析

基于補充了遺漏地震的完整地震目錄（截至2017年9月30日），以天為單位，逐漸增大統(tǒng)計時窗，分別利用擬合優(yōu)度法和最大似然法計算統(tǒng)計時段的完備震級和b值，分析余震序列的完備震級和b值隨時間的演化特征。主震發(fā)生后短時間內(nèi)的完備震級較高，出現(xiàn)3.5級的突跳，這可能與主震尾波的嚴重干擾有關，隨著時間的推移，完備震級緩慢降低并趨于穩(wěn)定，在1.2級左右波動，且呈現(xiàn)準周期的波動（圖6），這可能與晚間噪聲水平較低有關（侯金欣等，2017）。

圖6 地震序列完備震級和b值隨時間的演化

研究表明，b值與研究區(qū)域的應力狀態(tài)（Chan et al，2012；易桂喜等，2011）等因素有關，序列的低b值反映了震源區(qū)的高應力狀態(tài)。主震前的低b值反映了震源區(qū)應力積累到了較高的水平，伴隨著主震前地震的密集增強后而發(fā)震，此時可能正處于主震前應力加速釋放的臨震階段。主震發(fā)生后，隨著震源區(qū)應力的釋放，b值快速升高，8月15日以后逐漸平緩，在0.5左右波動，并仍有緩慢升高的趨勢（圖6），表明余震序列的衰減過程減緩，而該區(qū)的背景b值為1.0左右（易桂喜等，2011），說明余震序列仍會在未來較長一段時間內(nèi)處于衰減的狀態(tài)。

4 討論與結論

本文利用模板匹配技術，結合 GPU并行運算加速，識別了 2017年 8月 8日九寨溝MS7.0地震序列在主震前7天和震后5天的遺漏地震，檢測地震事件數(shù)量是臺網(wǎng)可定位地震事件的3倍多，將余震的完備震級從1.6下降到1.4級，并基于補充了遺漏地震的完整地震目錄，分析了余震序列的活動特征。

以九寨溝MS7.0地震為中心，對半徑為100km的范圍進行統(tǒng)計，2008年汶川MS8.0地震后地震活動開始明顯增強，尤其是2013年蘆山MS7.0地震以來地震活動進一步增強，一直處于較活躍狀態(tài)。但2017年以來，該區(qū)域地震活動明顯減弱，特別是3月之后，地震活動處于低活動背景，這與傳統(tǒng)意義上強震前地震活動由增強到平靜的過程相符。

從完備目錄的震級-時間分布來看，此次地震序列為前震-主震-余震型。前震序列存在平靜—增強的過程，特別是主震前的地震活動短時間內(nèi)增強，明顯高于該區(qū)域2017年以來每天1次的地震日次均值，這種高于背景地震活動水平的增強可能是應力的積累導致震源區(qū)發(fā)生微小破裂所致，這一階段可能是深部更大的斷層發(fā)生破裂之前加速蠕動的臨震階段。與之對應，主震前序列的b值處于低值水平，表明了震源區(qū)處于較高的應力狀態(tài)，而主震發(fā)生后b值隨時間的推移快速回升，維持在0.5左右，仍存在緩慢升高的趨勢。與之相反，序列的完備震級在主震發(fā)生后短時間內(nèi)達到最高值，這可能與主震尾波的嚴重干擾有關，隨后，完備震級緩慢降低并逐漸穩(wěn)定在1.2級左右。與該區(qū)域1.0的背景b值相比，目前序列的b值仍處于低值水平，從b值緩慢升高的趨勢來看，余震序列仍會在未來較長一段時間處于衰減的狀態(tài)。

本文僅對九寨溝MS7.0地震序列震前、震后幾天的連續(xù)波形進行掃描，并基于發(fā)震前后幾天的完備目錄結合b值進行初步分析。對于序列衰減系數(shù)和其他地震參數(shù)的分析，以及對于前震序列的相關推斷的論證，將在后續(xù)工作中進一步研究。