余娜 張曉清 袁伏全
青海省地震局,西寧市興海路1號 810001
地震序列是地震活動的主要表現形式之一,是統計地震學的重要研究對象。不同的構造區域、地震序列類型、震源區構造應力水平和區域大地熱流值等均可能表現為地震序列參數的差異(Kagan et al,2010),且余震序列的衰減特征和激發余震的能力可通過序列的統計參數來表征,因此,獲得準確可靠的地震序列參數對快速判定地震序列類型、研究震源區特征和評估后續地震危險性等均具有重要的參考價值(蔣海昆等,2007;Guo et al,1997)。目前國際上對地震序列參數的計算主要采用時間序列的“傳染型余震序列”(epidemic type aftershock sequence,簡稱 ETAS)模型(Ogata,1988、1989;Zhuang,2011)。蔣海昆等(2006)系統研究了中國大陸5級以上地震序列的統計特征,并從物理上討論了不同條件下的序列衰減和余震激發問題。蔣長勝等(2013a、2013b、2013c)對1976年唐山地震序列、2013年岷縣漳縣MS6.6地震序列和蘆山 MS7.0地震序列利用ETAS模型進行了研究。余娜等(2016a、2016b)采用時間序列的ETAS模型計算了2010年玉樹MS7.1、2016年門源MS6.4地震序列。
青海地區是西北地震活動最為活躍的地區之一,是地震學家高度關注的地震危險區域,歷史上多次強震的發生,積累該區域內更多、更可靠的地震序列參數,有利于震后地震趨勢估計和余震危險性研究。
用于擬合地震序列的 ETAS模型(Ogata,1988、1989;Zhuang,2011),其假設所有的余震均可按照大森-宇津公式(Omori,1894;Utsu,1961)激發自己的余震

式中,f(t)為t時刻對應的地震序列發生率;t為主震發生后的離逝時間;μ為背景地震發生率;p表示余震序列衰減的快慢;c為主震后余震頻次達到峰值時所對應的時間;K為常數,用于描述余震的活躍程度。ETAS模型假定主震的發生為初始零時刻,在其后的一個觀測時間段[0,T]內的地震序列{(ti,Mi);i=1,2,…,N}的強度函數可表示為(Ogata Y,1988)

其中,Mi和ti分別為第i個事件的震級和發生時間;M0為參考震級,一般可取截止震級。式(2)中的參數p表示地震序列衰減的快慢,p越大衰減越快,反之則越慢;參數α表示觸發次級余震的能力(Ogata,1989),α越大表示觸發次級余震的能力越弱,反之亦然。對于震群型序列,一般情況下 α<1,而當地震序列中無明顯的被激發的次級余震時,α一般大于 1(Ogata,2001);參數c為數值極小的正常數,一般用于調節當地震序列在主震發生后極短時間內余震記錄不完整,以及避免式(2)中分母為0的情況。為確保參數擬合的穩定性,當區域內背景地震發生率不高時,可通過設置背景地震的發生率μ=0來降低擬合難度。
本文選取的“青海地區”的空間范圍為 32.0°~39.0°N,89.5°~103.0°E。該研究區域位于青藏高原中北部及東北緣,地質構造規模大、活動性強,主要有阿爾金斷裂帶、祁連山-河西走廊活動斷裂帶、昆侖山斷裂帶、巴顏喀拉斷裂帶和唐古拉斷裂帶等,以NNW、NW-NWW向斷裂帶最為發育。地震具有頻次高、強度大、分布廣的特點(馬玉虎等,2006;曾秋生,1999)。2009年以來研究區主震震級MS≥5.0地震有12個,其中,5.0≤MS<6.0的有8個,6.0≤MS<7.0的有3個,MS≥7.0的有1個。受ETAS模型計算對地震序列質量的要求,本文使用了記錄較好的5級以上地震序列共9個(圖1),分別為2009年8月28日大柴旦MS6.4、2010年3月24日西藏聶榮MS5.8、2010年4月14日玉樹MS7.1、2011年6月26日囊謙MS5.2、2013年1月30日雜多MS5.2、2013年9月20日門源MS5.1、2015年11月23日祁連MS5.0、2016年1月21日門源MS6.4和2016年10月17日雜多MS6.2地震。
本研究所采用的地震目錄為中國地震臺網中心提供的《全國統一正式目錄》。青海地區地震目錄存在著明顯的空間不均勻現象,說明不同時期、不同地區地震記錄水平存在明顯的差異。在實際的地震序列選取中,采用了蔣長勝等(2014)區分地震序列與背景地震的時-空分布界限的“自然邊界法”,通過緯度-時間圖、經度-時間圖及震中分布圖相結合的方式選取地震序列。圖2給出了2009年大柴旦MS6.4地震序列的選取結果。依照上面的選取方法,選取了其余8個地震序列。此外,由于在一些強震發生后的短期內,主震的波形振幅較大,面波等波列持續時間較長,可出現將隨后大量發生的震級較小的余震“淹沒”而導致余震區甚至更大范圍內的地震監測能力顯著降低(Akaike,1974)的現象。為確保地震序列的完整性,采用“震級-序號”法(Huang,2006;蔣長勝等,2014;Zhuang et al,2012)確定了每個地震序列完整性震級Mc(圖3)。“震級-序號”法按地震發生時間的先后順序排序,地震密度較大的區域的連線大致為截止震級Mc的時序變化。

圖1 研究區9個地震序列主震的空間分布

圖2 2009年大柴旦MS 6.4地震序列的選取結果
ETAS模型參數一般使用最大似然法進行估計。在擬合時間[S,T]范圍內(0<S<T),似然函數L的形式為

將式(2)帶入式(3),即可對序列參數[μ,K,c,α,p]進行最大似然估計。作為示例,圖4給出了2009年海西MS6.4地震序列在Mc=ML1.9情況下的ETAS模型參數結果,獲得ETAS模型參數為 μ=0.0000,K=0.0280,c=0.0035,α=1.2411和 p=1.2650。為進一步考察 ETAS模型的擬合效果,常使用“轉換時間”域的“殘差分析”法(Ogata,1988)。“殘差分析”方法可將復雜的地震序列{ti}轉換為齊次泊松過程,并在“轉換時間”域{τ}中進行擬合效果分析,相應的轉換過程為

這種經過轉化的{τ}域的數據被稱為“殘差點過程”(RPP)。若在殘差分析中,轉換時間域的RPP在{τ}上的累計頻次近似直線,并與理論直線較為接近(Zhuang et al,2012),則ETAS模型擬合效果較好。以2009年大柴旦MS6.4地震序列為例,圖4(a)給出了主震之后0~91天和Mc=ML1.9時的“殘差分析”結果,可以看出,序列累計地震數與理論ETAS擬合曲線較為接近,ETAS模型對此地震序列的擬合較為理想;圖4(b)為相應的{τ}域的 M-τ分布。研究區9個地震序列參數的估算結果見表1。
ETAS模型參數與截止震級的選取有一定的關系(蔣長勝等,2013a、2013b)。本文以2009年8月28日大柴旦MS6.4、2010年4月14日玉樹MS7.1地震為例,研究了不同截止震級對地震序列參數的影響。2009年大柴旦MS6.4地震選取截止震級Mc=ML1.8、1.9、2.1,發現地震序列參數α值和p值在Mc=ML2.1時變化較小,誤差也小(圖5)。2010年玉樹MS7.1地震選取截止震級Mc=ML1.5、1.6、…、2.1,當截止震級 Mc從 ML1.5到2.1逐漸增大時,α值和p值隨截止震級的變化均有較大變化,參數k值逐漸減小,α值總體上有增加的趨勢,p值隨截止震級的增大而減小,但是變化幅度不大。因此,在考察地震序列參數的早期特征時,須考慮在何種截止震級條件下進行討論。
由于板內地震所處的構造運動速率不同,故地震序列持續時間存在明顯差異。因此,所選取的地震序列的擬合截止時間不同,獲得的ETAS模型參數可能不同,擬合參數的標準差也將存在差異(Iwata,2008)。以 2010年4月14日青海玉樹 7.1級地震為例,設定 tend=[1.0,2.0,3.0,…,42.0],當 Mc=ML1.5時,ETAS模型參數 α值和 p值均在震后 14天出現突跳變化,標準差幅值較大,尤其是激發次級余震能力的α值變化最為劇烈,p值在震后18天內有一定幅度波動。b值變化相對穩定,但在震后14天內有逐漸增加的變化,隨著序列持續時間的增加α值減小。地震早期ETAS模型參數的α值、p值變化幅度較大,其后,趨于相對穩定。

圖3 研究區9個地震序列的“震級-序號”法處理結果

圖4 2009年8月28日大柴旦MS 6.4地震序列擬合情況的“殘差分析”

表1 研究區9個地震序列主震參數和地震序列ETAS模型參數值

圖5 2009年8月28日大柴旦MS 6.4地震 ETAS模型計算結果
地震序列參數的早期變化比較劇烈,反映了地震序列主震發生后應力快速調整的過程,余震也主要集中發生在地震早期階段。為了研究地震序列參數在時間上的變化特征,本文將地震序列不同截止條件下的ETAS模型參數穩定時間和震級進行統計(表2),可以看出,地震序列以不同的截止震級作為研究對象時,地震序列的穩定時間 Tα和Tp是不同的,也就是說,在一個地震序列中,地震序列參數α值和p值的穩定時間Tα和 Tp與截止震級沒有明顯的正相關關系,當 Mc=ML1.8時,Tα=19,Tp=19;當 Mc=ML1.9時,Tα=17,Tp=21;當 Mc=ML2.1時,Tα=23,Tp=19。地震序列參數 α值和 p值的穩定時間Tα和 Tp與主震震級間沒有明顯的相依關系,主震震級較大的其穩定時間也可能較短,震級較小的地震穩定時間反而較長。比如玉樹MS7.1地震,地震序列參數α值和p值的穩定時間為:當Mc=ML2.0時,Tα=15,Tp=15;雜多 MS5.2地震,當 Mc=ML2.0時,Tα=18,Tp=13;反映了玉樹 7.1地震震后應力迅速調整的特征。
研究區不同震級檔地震序列參數的最長穩定時間和最短穩定時間分別為:地震序列主震震級 MS5.0~6.0,Tα,min=9,Tα,max=47,Tp,min=13,Tp,max=45;主震震級 MS6.0~7.0,Tα,min=17,Tα,max=57,Tp,min=11,Tp,max=26。

表2 地震序列參數穩定時間統計表
為系統考察青海地區地震序列的ETAS模型,并對地震序列的震后早期特征進行分析,本文運用中國地震臺網中心提供的《全國統一正式目錄》,對青海地區的地震序列進行了篩選,發現青海地區地震目錄存在著明顯的空間不均勻現象,不同時期、不同地區地震記錄水平存在明顯的差異。在實際的地震序列選取中采用“自然邊界法”選取了研究區2009年以來的9個地震序列,利用“震級-序號”法確定了每個地震序列的完整性震級Mc,使用最大似然法對每個地震序列進行參數估算,獲得了如下認識:
(1)為檢測地震序列參數的穩定性,考察了截止震級和擬合截止時間的選取對ETAS模型參數影響。研究發現,截止震級Mc對α、k和p值有一定的影響。在考察地震序列參數的早期特征時,須考慮在何種截止震級下進行討論;地震早期ETAS模型參數的α值、p值變化幅度較大,反映了地震主震發生后震源區應力的快速調整過程,其后,地震參數的α值和p值趨于相對穩定。
(2)青海地區地震序列的ETAS模型參數α值、p值的穩定時間Tα、Tp與截止震級間沒有明顯的正相關關系,與主震震級間也無明顯的相依關系,主震震級較大的地震事件穩定時間可能較短,震級較小的地震事件其穩定時間也有可能較長。
(3)青海地區不同震級檔地震序列參數的最長穩定時間和最短穩定時間分別為:地震序列主震震級 MS5.0~6.0,Tα,min=9,Tα,max=47,Tp,min=13,Tp,max=45;主震震級 MS6.0~7.0,Tα,min=17,Tα,max=57,Tp,min=11,Tp,max=26。
由于研究區域強震數目和ETAS模型計算條件所限,僅對9個地震序列進行了估算,并對其時間特征進行了討論,且所得出的結論僅反映部分統計特征,今后還需要繼續積累更多的震例進行更深入研究。
致謝:研究中使用了中國地震臺網中心“全國地震編目系統”提供的“統一正式目錄”,日本統計數理研究所莊建倉教授編制、提供了本文所使用的ETAS模型和余震短期概率預測程序,在此一并表示感謝。