張盛峰 張永仙 范曉易

1)南京大學(xué)地球科學(xué)與工程學(xué)院, 南京 210023 2)中國地震局地震預(yù)測研究所, 北京 100036 3)江蘇省地震局, 南京 210000

0 引言

2020年1月19日21時27分, 新疆喀什地區(qū)伽師縣(39.83°N, 77.21°E)發(fā)生了MS6.5地震, 震源深度為16km。此次地震發(fā)生在青藏高原西構(gòu)造結(jié)的NE側(cè), 與中國地震科學(xué)實驗場所在的青藏高原東構(gòu)造結(jié)遙相呼應(yīng)。為加強對中國地震科學(xué)實驗場區(qū)外發(fā)生的典型地震的研究, 中國地震局開展了針對此次地震的虛擬科考工作。其中, 針對這一地區(qū)的背景地震活動和觸發(fā)地震活動進行時空特征分析, 將為更好地回答與這一地震相關(guān)的科學(xué)問題提供重要參考。為此, 本工作利用時-空ETAS模型分析地震序列特征, 從統(tǒng)計地震學(xué)的角度探索本地區(qū)背景地震活動和觸發(fā)地震活動的分布規(guī)律, 以進一步深入認識此次地震的相關(guān)性質(zhì)。

通過傳統(tǒng)方法判斷某個地震事件屬于背景地震活動還是被觸發(fā)的余震活動往往較為困難(Helmstetteretal., 2003), 通常需要對其及主震所處的時空范圍進行對比加以確定(Keilis-Boroketal., 1980)。然而, 地震活動的叢集特征在不同的研究區(qū)域存在差異, 這為有效識別余震活動帶來了困難。因此, 若干識別方法被相繼提出, 如基于窗口(windows-based)(Utsu, 1970；Gardneretal., 1974)和基于關(guān)聯(lián)(link-based)(Reasenberg, 1985)的去叢方法等。在使用這些方法進行判定時, 優(yōu)化時空窗參數(shù)或關(guān)聯(lián)距離的選取較為依賴于研究人員的主觀經(jīng)驗, 同時, 一些不應(yīng)該被判定為余震的事件常被確定性地刪除, 從而導(dǎo)致有用信息受到損失(Zhuangetal., 2005)。以傳染型余震序列(Epidemic Type Aftershock Sequence model, ETAS模型)為代表的統(tǒng)計模型主要是以分支點過程的形式研究地震叢集的時空發(fā)生行為(Kagan, 1991；Musmecietal., 1992；Rathbun, 1996；Ogata, 1998, 2004；Consoleetal., 2001, 2003；Zhuangetal., 2002；Ogataetal., 2003)。通常, 這些模型將地震活動性分為2部分, 即背景地震活動和叢集活動。其中每個地震事件, 無論是來自背景成分(通常假定為時空泊松過程, 穩(wěn)定或非平穩(wěn), 均質(zhì)或非均質(zhì)), 還是由另一個事件觸發(fā)生成, 都會根據(jù)一些分支規(guī)則產(chǎn)生(觸發(fā))自己的后代(余震)。為了獲得準確去叢的目錄, Zhuang等(2002, 2004)提出了一種隨機去叢方法作為傳統(tǒng)方法的替代方法。在這種方法中, 不再劃定地震是背景事件還是由其他事件觸發(fā), 而認為每個事件都有可能是背景事件或其他事件觸發(fā)的直接后代, 根據(jù)用于描述地震叢集特征的一些模型估計每個事件的概率。由于這一方法所基于的時-空ETAS模型可以較好地描述地震活動的行為, 并且可在概率分布上定量分析背景地震活動和觸發(fā)地震活動, 因此在研究背景地震或觸發(fā)地震方面得到了廣泛應(yīng)用(蔣長勝等, 2010, 2013；龍鋒等, 2010；蔣海昆等, 2012；Kawamuraetal., 2013；Yoderetal., 2014；Kattamanchietal., 2017)。

1 時-空ETAS模型及隨機除叢方法

1.1 時-空ETAS模型

目前, 已有多種分支過程(Branching Processes)模型用于描述地震活動的時-空叢集特征(Ogata, 1988, 1992, 1998；Kagan, 1991；Musmecietal., 1992；Rathbun, 1993；Consoleetal., 2001, 2003), 這些模型一般用條件強度函數(shù)的形式(Daleyetal., 2007)表示地震發(fā)生率。本工作使用的隨機除叢法主要基于Ogata(1998)給出的時-空ETAS模型, 其涉及的基本原理已在相關(guān)文獻中有系統(tǒng)描述(Omori, 1894；Utsu, 1970；Ogata, 1978, 1998；Rathbun, 1993, 1996；Consoleetal., 2003；Daleyetal., 2003；Ogataetal., 2006；Zhuangetal., 2006), 這里不再贅述。

1.2 瘦化算法和隨機除叢法

隨機除叢方法的關(guān)鍵就是對點過程的 “瘦化”運算, 地震i對其后發(fā)生的在地震j處(tj,xj,yj)的總地震發(fā)生率的相對貢獻可表示為

(1)

其中,

ζi(tj,xj,yj)=k(Mi)g(t-ti)f(x-xi,y-yi;Mi)

(2)

表示地震i觸發(fā) “子震”過程的發(fā)生率。因此可以將ρij看作地震j被地震i觸發(fā)的概率, 即地震j作為地震i“子震”的概率。同理, 地震j作為背景地震的概率為

(3)

而地震j被之前地震觸發(fā)的概率, 即叢集地震概率可表示為1-φj。

Zhuang等(2002)的算法分為同時求取背景地震活動強度和模型參數(shù)的迭代算法以及隨機除叢法2部分?？梢钥闯? 與傳統(tǒng)的刪除余震算法所不同的是, 隨機除叢法為每個地震事件的性質(zhì)賦予1個概率值, 利用 “家譜”的形式描述叢集地震, 即每個 “子震”根據(jù)相關(guān)概率隨機地找到自己的 “母震”, 提供了可針對去叢效果進行定量評估的不確定度。

圖1 1970年1月1日—2020年6月21日研究區(qū)ML>3.5地震的空間分布Fig.1 Spatial distribution of events above ML3.5 during Jan.1, 1970 to Jun.21, 2020 in the study region.黑色實線為斷裂，黃色、藍色和紅色實心圓分別表示震級范圍ML3.5～5.9、ML6.0～6.4和6.5級及以上地震事件, 左上角給出了本研究區(qū)所在的位置

2 研究區(qū)域及數(shù)據(jù)

本工作使用了中國地震臺網(wǎng)中心產(chǎn)出的震源區(qū)(38.83°～40.83°N, 76.21°～78.21°E)1970年1月1日—2020年6月21日的正式地震目錄數(shù)據(jù), 區(qū)域地震空間分布見圖1，可見研究區(qū)北部的活動斷裂較多, 多呈NE向分布, 南部地區(qū)多為平坦的沙漠地區(qū), 而本地區(qū)ML3.5以上的地震事件在空間上同樣呈現(xiàn)明顯的南、北分布差異： 北部地區(qū)的地震活動分布相對分散且廣泛, 而南部地區(qū)則相對集中,MS6.0以上地震主要分布于中部地區(qū)。表1 和圖2e給出了本地區(qū)發(fā)生的4個MS6.5以上地震事件的目錄及其在累計頻次曲線上所處的位置。

表 1 研究區(qū)1970年1月1日—2020年6月21日MS6.5以上地震事件Table 1 List of events above MS6.5 from Jan.1, 1970 to Jun.21, 2020 in the study region

科學(xué)評估地震目錄的完整性水平并較為穩(wěn)妥地選擇所需要的截止震級MC是較好地應(yīng)用時-空ETAS模型的前提條件。本工作使用的目錄包含ML和MS2種震級標度, 一般4.5級以上取MS震級, 以下取ML震級。為科學(xué)評估本地區(qū)地震目錄的完整性水平, 這里采用與地震目錄相關(guān)的若干可視性方法對其進行分析。圖2 給出了研究區(qū)所包含地震事件的震級-時間、震級-序號、震級累計頻次等統(tǒng)計結(jié)果。從圖2a、b中可以看出, 自2000年以來研究區(qū)的地震監(jiān)測水平得到大幅提升, 尤其是2010年以后可記錄到研究區(qū)0級以下事件；從圖2c、d中可以看出, 研究區(qū)發(fā)生的1996年MS6.8、1997年MS6.5和2003年MS6.7地震均觸發(fā)了一定數(shù)量的余震事件, 這一現(xiàn)象也在圖2e中的不同震級事件累計頻次隨時間的變化曲線中有所顯示, 但1970年以來的累計頻次曲線同樣說明從ML3.5開始的累計頻次曲線不再明顯地受幾個強震所影響。圖2d、f所示的震級-序號分布顯示，在1997年4月11日MS6.5地震震后短時間內(nèi), 由于背景噪聲或強余震波形的影響, 沒有充分識別ML3.5以下的余震事件, 導(dǎo)致這些事件缺失, 造成震后的 “空白三角”現(xiàn)象。由于時-空ETAS模型進行地震活動擬合時對地震事件的數(shù)量有一定要求, 若綜合考慮整個研究區(qū)域1970年以來的地震序列,MC取ML3.5可較好地保證地震目錄的完整性程度, 如圖2a—d中的紅色橫線所示?？紤]到模型所需地震數(shù)量和完整性水平二者的平衡性, 本文將ML3.5作為模型擬合所使用的參數(shù)。

圖2 研究區(qū)1970年1月1日—2020年6月21日地震目錄的若干分析Fig.2 Analysis for the earthquake catalog during Jan.1, 1970 to Jun.21, 2020.a 全部地震事件的震級-時間分布；b 全部地震事件的震級-序號分布；c ML3.0以上地震事件的震級-時間分布；d ML3.0以上地震事件的震級-序號分布；e 不同震級范圍的地震事件隨時間的累計頻次變化曲線；f ML3.5以上地震事件的震級-序號 分布。a—d中紅色橫線表示ML=3.5標尺的位置；e中垂直虛線表示研究區(qū)發(fā)生的MS6.5以上地震事件的發(fā)震時刻

3 模型參數(shù)估計與地震活動強度

除截止震級MC外, 時-空ETAS模型還需要設(shè)定若干用于計算的參數(shù), 包括目錄起始時間、模型起算時間、模型參數(shù)初始值和最小帶寬等。圖3 為截止震級取ML3.5，擬合起始時刻tc取1970年1月1日時計算得到的時-空ETAS模型參數(shù)隨擬合次數(shù)的變化情況?？梢? 經(jīng)過10次模型擬合, 最終取得了較為穩(wěn)定的模型參數(shù), 而獲取穩(wěn)定可靠的模型參數(shù)又是針對地震活動進行隨機除叢處理的基礎(chǔ)。

圖3 截止震級MC取ML3.5、擬合起始時刻tc取1970年1月1日時計算得到的時-空ETAS模型參數(shù)隨擬合次數(shù)的變化情況Fig.3 Variation of model parameters with fitting times using the input parameters of MC=ML3.5 and tc=1970-01-01.

獲取可靠的模型參數(shù)后, 進而可計算得到研究區(qū)空間總地震發(fā)生率、背景地震發(fā)生率、叢集地震發(fā)生率和叢集率空間分布, 結(jié)果如圖4 所示?？梢钥闯? 研究區(qū)北部的背景地震活動水平較高且分布相對均勻, 南部多為觸發(fā)的叢集活動。對于整個區(qū)域而言, 叢集地震活動對總體地震活動的貢獻最大。圖4d中的叢集率結(jié)果顯示, 南部區(qū)域的叢集率接近1.0, 但范圍更大?？梢? 雖然該區(qū)域的背景地震發(fā)生率不高, 但一旦發(fā)生較強的地震事件, 后期易觸發(fā)一定數(shù)量的余震事件, 同時也容易造成震后短時間內(nèi)的余震記錄缺失現(xiàn)象。這一區(qū)域內(nèi)的地震互相觸發(fā)作用及斷層活動機制也在前人的研究中有所提及(張竹琪等, 2008)。

圖4 時-空ETAS模型計算得到的空間總地震發(fā)生率(a)、背景地震發(fā)生率(b)、叢集地震發(fā)生率(c)和叢集率分布(d)Fig.4 Spatial distribution of total seismicity rate(a), background seismicity rate(b), clustering seismicity rate(c), and cluster ratio(d).a—c 單位為事件個數(shù)/((°)2·d),每個網(wǎng)格點的叢集率計算誤差可通過總地震發(fā)生率或地震事件空間分布密度粗略估計

圖5 隨機除叢方法對背景地震活動及地震事件觸發(fā)作用的分析結(jié)果Fig.5 Analysis results of background events and their triggering contribution to earthquake events in the study region using the stochastic declustering method.a 以色標形式表示背景地震事件概率水平的地震活動空間分布；b 伽師MS6.5地震發(fā)生前地震活動的空間分布, 其中色標表示這些事件對此次MS6.5地震(黑色圓圈)的觸發(fā)貢獻, 紅色圓點表示對其觸發(fā)貢獻最大的2020年1月18日MS5.3地震

4 背景地震事件與觸發(fā)地震事件分析

利用基于時-空ETAS模型的隨機除叢方法可區(qū)分背景地震活動與觸發(fā)叢集地震活動。圖5a為給每個事件賦予背景事件概率后的地震活動空間分布, 若將0.5作為區(qū)分背景地震事件和觸發(fā)地震事件的閾值, 則可看出背景地震事件大多較為均勻地分布于北部區(qū)域, 而受到觸發(fā)的叢集事件主要集中于南部區(qū)域, 與叢集率高值區(qū)域及地震事件空間分布高密度區(qū)域一致。經(jīng)過計算, 此次MS6.5地震被觸發(fā)的概率達99%, 說明其發(fā)生過程存在較強的觸發(fā)作用。圖5b給出了對此次MS6.5地震具有觸發(fā)貢獻的地震活動空間分布, 可以看出大部分事件對此次MS6.5地震的觸發(fā)概率均較低, 僅2020年1月18日發(fā)生的MS5.3地震(圖中紅色圓點所示事件)對此次地震的觸發(fā)貢獻達到了94%(與圖6c，d對應(yīng))。

為研究其他事件對此次地震觸發(fā)貢獻的時序變化, 與圖5 相對應(yīng), 圖6 給出了地震事件的背景地震概率及其他事件對2020年1月19日MS6.5地震觸發(fā)貢獻隨地震序號和時間的變化分布。從圖6a、b中可以看出, 1996年MS6.8、1997年MS6.5和2003年MS6.7地震均明顯觸發(fā)了大量余震事件, 而主要分布于北部地區(qū)的地震事件的背景地震概率結(jié)果在序號和時間上的分布較為均勻；圖6c、d為此次地震被觸發(fā)的概率隨地震序號和時間分布, 顯示結(jié)果與圖5b一致。

圖6 基于隨機除叢方法計算得到的地震事件背景地震概率隨地震序號(a)、時間(b)變化的分布和其他事件對2020年6月26日MS6.5地震觸發(fā)貢獻隨地震序號(c)、時間(d)變化的分布Fig.6 Probability of each event as a background event vs.sequence number(a)and time(b), and contribution of others to the Jun.26, 2020 MS6.5 event vs.sequence number(c)and time(d).

圖7 研究區(qū)自1970年1月1日以來ML3.5以上地震事件對其他個體的觸發(fā)能力的空間分布Fig.7 Spatial distribution of triggering contribution of each event above ML3.5 to others in the study region since Jan.1, 1970.

5 不同地震事件 “干預(yù)”其他事件的能力分析

與分析此次地震作為被觸發(fā)余震的概率類似, 利用隨機除叢方法同樣可以得到研究地區(qū)地震事件個體對其他事件的觸發(fā)貢獻水平, 以便更清楚地看到哪些事件更 “容易”干預(yù)其他事件的發(fā)生過程。圖6 給出了研究區(qū)1970年以來ML3.5以上地震對其他事件觸發(fā)貢獻平均水平的空間分布, 可以看出這些容易觸發(fā)其他個體的事件主要分布于此次MS6.5地震震源區(qū)周圍, 即位于高背景地震活動的北部區(qū)域和多觸發(fā)叢集事件的南部區(qū)域之間。圖8 給出了與圖7 相對應(yīng)的地震觸發(fā)能力隨時間和地震序號變化的圖像, 可見在與圖7 所對應(yīng)的地震事件中, 2020年4月21日ML4.1地震對之后地震事件發(fā)生過程的平均觸發(fā)貢獻達0.505, 表2 給出了 “干預(yù)”其他個體的能力排序靠前的幾個事件的信息。值得一提的是, 這一結(jié)果是基于當前所選區(qū)域和使用數(shù)據(jù)的條件下得到的模型計算結(jié)果, 與實際情況是否吻合, 仍需要開展更加細致的研究和調(diào)查。

表 2 研究區(qū)自1970年1月1日以來的地震事件對其他個體觸發(fā)貢獻平均水平排序靠前的地震事件Table 2 Cases of top average contribution of each event triggering others since Jan.1, 1970 in the study region

圖8 基于隨機除叢方法計算得到的地震事件對其他個體觸發(fā)貢獻隨地震序號(a)和時間(b)的分布Fig.8 Contribution of each event to others vs.sequence number(a)and time(b).

6 討論和結(jié)論

本工作針對新疆伽師地區(qū), 利用時-空ETAS模型分析了震源區(qū)1970年以來的地震活動特征, 基于隨機除叢方法在概率水平上區(qū)分了背景地震活動和觸發(fā)余震地震活動, 結(jié)果顯示研究區(qū)北部與南部分別呈現(xiàn)均勻分布的背景地震活動和觸發(fā)叢集活動, 且叢集地震活動為總體地震活動的主要成分；計算顯示，此次MS6.5地震被觸發(fā)的概率為99%, 2020年1月18日MS5.3地震對其產(chǎn)生的觸發(fā)作用較為明顯, 達94%；在對其他個體具有觸發(fā)貢獻的地震事件中, 2020年4月21日ML4.1地震 “干預(yù)”其他事件發(fā)生過程的平均能力最高, 達0.505, 而這一結(jié)果是否與實際情況吻合, 需要通過其他手段進行詳細研究。值得一提的是, 本工作主要從統(tǒng)計模型的角度對研究區(qū)的群體事件進行了分析, 從概率出發(fā)分析了事件之間可能的觸發(fā)關(guān)系, 而伽師地區(qū)由于受到所處西構(gòu)造結(jié)的影響, 擠壓變形明顯, 強烈的地質(zhì)活動特征與此地區(qū)多個強震的發(fā)生過程有較大關(guān)系(趙翠萍等, 2003；徐錫偉等, 2006)。

以ETAS模型為基礎(chǔ)的隨機除叢方法可從概率水平上判定地震事件屬于背景地震活動還是叢集地震活動, 這種工作思路成為銜接純粹的科學(xué)研究和中國地震會商業(yè)務(wù)體系的有效橋梁。雖然模型的部分結(jié)果受到與地震目錄相關(guān)的因素影響, 如本文計算顯示2020年4月21日ML4.1地震對之后地震事件的觸發(fā)作用相對較高, 這一結(jié)果受到選取目錄的截止時刻影響較大, 目前還無法確定是否與實際情況相符, 但對于一些利用其他手段無法有效解決的問題, 本模型的計算結(jié)果仍會給出參考。在ETAS模型的不同應(yīng)用場景中, 當前應(yīng)用較多的情況仍然是將地震事件當作點源看待, 但隨著這一模型的發(fā)展, 近年來也出現(xiàn)了增加震源深度約束的3D-ETAS模型(Ogataetal., 2019)和斷裂空間形態(tài)的新模型(Guoetal., 2015), 而對于地震發(fā)生后需要快速了解觸發(fā)概率等信息的地震會商工作而言, 建設(shè)與這些新模型相關(guān)的斷裂數(shù)據(jù)庫及提高震源深度的測定精度等是這些新模型能夠得到較好應(yīng)用的前提條件。與模型中背景地震活動水平(μ)設(shè)置為固定值或可變值的操作類似, 判定某個地震事件屬于背景地震活動還是觸發(fā)地震活動的關(guān)鍵因素——概率閾值, 是否也需要根據(jù)特定的環(huán)境而進行調(diào)整？ 例如, 在某些條件下, 當觸發(fā)概率>0.7則認為該事件為被觸發(fā)的余震活動, <0.3則認為其是背景地震活動等。以上問題均需要針對模型和應(yīng)用地區(qū)的地震活動性特點開展進一步研究。

致謝本工作是針對2020年1月19日新疆于田MS6.5地震虛擬科考工作的一部分, 科考隊專家在震后提供的相關(guān)分析結(jié)果為作者更好地認識此次地震提供了條件；作者在日本數(shù)理統(tǒng)計研究所學(xué)習期間與莊建倉準教授、Ogata Y教授、郭一村博士等專家學(xué)習討論了時-空ETAS模型的相關(guān)問題；審稿專家為本文的修改完善付出了辛勤勞動。在此一并表示感謝!