基于時-空ETAS模型的新疆伽師地區背景及觸發地震活動的探討

張盛峰 張永仙 范曉易

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

0 引言

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

通過傳統方法判斷某個地震事件屬于背景地震活動還是被觸發的余震活動往往較為困難(Helmstetteretal., 2003), 通常需要對其及主震所處的時空范圍進行對比加以確定(Keilis-Boroketal., 1980)。然而, 地震活動的叢集特征在不同的研究區域存在差異, 這為有效識別余震活動帶來了困難。因此, 若干識別方法被相繼提出, 如基于窗口(windows-based)(Utsu, 1970；Gardneretal., 1974)和基于關聯(link-based)(Reasenberg, 1985)的去叢方法等。在使用這些方法進行判定時, 優化時空窗參數或關聯距離的選取較為依賴于研究人員的主觀經驗, 同時, 一些不應該被判定為余震的事件常被確定性地刪除, 從而導致有用信息受到損失(Zhuangetal., 2005)。以傳染型余震序列(Epidemic Type Aftershock Sequence model, ETAS模型)為代表的統計模型主要是以分支點過程的形式研究地震叢集的時空發生行為(Kagan, 1991；Musmecietal., 1992；Rathbun, 1996；Ogata, 1998, 2004；Consoleetal., 2001, 2003；Zhuangetal., 2002；Ogataetal., 2003)。通常, 這些模型將地震活動性分為2部分, 即背景地震活動和叢集活動。其中每個地震事件, 無論是來自背景成分(通常假定為時空泊松過程, 穩定或非平穩, 均質或非均質), 還是由另一個事件觸發生成, 都會根據一些分支規則產生(觸發)自己的后代(余震)。為了獲得準確去叢的目錄, Zhuang等(2002, 2004)提出了一種隨機去叢方法作為傳統方法的替代方法。在這種方法中, 不再劃定地震是背景事件還是由其他事件觸發, 而認為每個事件都有可能是背景事件或其他事件觸發的直接后代, 根據用于描述地震叢集特征的一些模型估計每個事件的概率。由于這一方法所基于的時-空ETAS模型可以較好地描述地震活動的行為, 并且可在概率分布上定量分析背景地震活動和觸發地震活動, 因此在研究背景地震或觸發地震方面得到了廣泛應用(蔣長勝等, 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), 這些模型一般用條件強度函數的形式(Daleyetal., 2007)表示地震發生率。本工作使用的隨機除叢法主要基于Ogata(1998)給出的時-空ETAS模型, 其涉及的基本原理已在相關文獻中有系統描述(Omori, 1894；Utsu, 1970；Ogata, 1978, 1998；Rathbun, 1993, 1996；Consoleetal., 2003；Daleyetal., 2003；Ogataetal., 2006；Zhuangetal., 2006), 這里不再贅述。

1.2 瘦化算法和隨機除叢法

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

(1)

其中,

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

(2)

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

(3)

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

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

圖1 1970年1月1日—2020年6月21日研究區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級及以上地震事件, 左上角給出了本研究區所在的位置

2 研究區域及數據

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

表 1 研究區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

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

圖2 研究區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中垂直虛線表示研究區發生的MS6.5以上地震事件的發震時刻

3 模型參數估計與地震活動強度

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

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

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

圖4 時-空ETAS模型計算得到的空間總地震發生率(a)、背景地震發生率(b)、叢集地震發生率(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 單位為事件個數/((°)2·d),每個網格點的叢集率計算誤差可通過總地震發生率或地震事件空間分布密度粗略估計

圖5 隨機除叢方法對背景地震活動及地震事件觸發作用的分析結果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地震發生前地震活動的空間分布, 其中色標表示這些事件對此次MS6.5地震(黑色圓圈)的觸發貢獻, 紅色圓點表示對其觸發貢獻最大的2020年1月18日MS5.3地震

4 背景地震事件與觸發地震事件分析

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

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

圖6 基于隨機除叢方法計算得到的地震事件背景地震概率隨地震序號(a)、時間(b)變化的分布和其他事件對2020年6月26日MS6.5地震觸發貢獻隨地震序號(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 研究區自1970年1月1日以來ML3.5以上地震事件對其他個體的觸發能力的空間分布Fig.7 Spatial distribution of triggering contribution of each event above ML3.5 to others in the study region since Jan.1, 1970.

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

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

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

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

6 討論和結論

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

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

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