全球MW≥8.0淺源地震的前震序列研究

薛 艷 解孟雨 劉 杰 莊建倉

1)中國地震臺網中心， 北京 100045 2)日本東京數理統計研究所， 東京190—8562

0 引言

近年來全球發生了多次8級以上大地震， 造成了嚴重的人員傷亡和巨大的財產損失。由于有儀器記錄的時間短， 因此有完整小地震記錄以來的8級地震震例較少。人們對地震序列的研究主要集中于5級以上地震(Wellsetal.， 1994； 蔣海昆等， 2007； Chenetal.， 2013， 2015)。有學者對全球大地震的地震序列進行了一定的研究(蘇有錦等， 2008， 2014； 呂曉健等， 2010)， 但尚未有涉及8級大地震前震序列特征的系統研究。

2011年3月11日日本本州以東海域9.1級地震震前2d發生了7.3級前震序列， 最大前震距主震約43km， 整個前震序列分布在主震NE方向100km范圍內。2012年4月11日蘇門答臘8.6級地震震前90d發生了7.2級前震序列， 最大前震與主震震中相距約20km， 前震序列分布在距主震震中66km的范圍內。2014年4月1日智利伊基克8.2級地震震前15d發生了6.7級前震序列， 最大前震與主震震中相距42km， 前震序列分布在主震周圍75km的范圍內。這3次大地震的前震活動引起了廣泛關注(Meredithetal.， 2011； Xueetal.， 2012； Brodskyetal.， 2014)。

在所有的地震短臨前兆中， 前震是學術界公認的預報強震最有效的指標之一(Jonesetal.， 1979)。陳颙等(2015)總結了前震序列的4個特點： 前震序列是一種高頻率的地震活動， 其活動頻率比背景地震的活動頻率高； 前震序列在主震之前(幾h—幾d)發生； 前震序列發生的位置與主震相同； 主震的震級比前震序列中任何地震的震級都大。由此可見， 前震的確定與其和主震的時間差和距離差的約束有關。Chen等(2015)系統研究了加利福尼亞地區5級以上地震的前震序列特征， 發現48%的5級以上地震在震前30d內， 距離主震5km范圍內至少發生了1次前震。

此外， 陳顒(1978)還發現前震序列的震源機制解類似， 而余震序列的震源機制解差異大， 并據此提出利用震源參數的一致性判別一個序列是前震序列還是正常的主震-余震序列的方法。倪四道等(2010)研究了2010年4月14日青海玉樹7.1級地震前約2h在距離主震震中僅2km處發生的4.7級地震序列的時空分布及波形特征， 認為該4.7級序列符合陳顒(1978)的理論， 是典型的前震序列。

Scholz(1968)在巖石破裂實驗中發現， 微破裂b值與壓應力呈負相關。針對一些特定震例的研究表明，b值會在主震前后發生變化， 與余震相比， 前震序列的b值偏低(Smith， 1981； Xueetal.， 2012)。這一特性與初始的高應力狀態隨著主震的應力降而降低的理論一致， 該結論常被用于識別前震(Scholz， 1968； 李全林等， 1978； 林邦慧等， 1994)。Gulia等(2019)研究認為， 前震b值顯著低于區域背景b值， 并提出了實時判別前震的 “交通燈”準則。本文將系統研究1976年以來全球MW≥8.0淺源地震前震序列的統計特征， 應用 “傳染型余震序列(epidemic-type aftershock sequence， ETAS)”模型計算前震和余震序列參數， 并討論參數的差異性。研究中使用了美國國家地震信息中心(NEIC)提供的全球地震目錄(1)http://earthquake.usgs.gov。， 震源機制解數據來源于美國哈佛大學(2)http://www.globalcmt.org/CMTsearch.html。。

1 地震目錄最小完備性震級的討論

美國國家地震信息中心(NEIC)提供了1973年以來全球地震目錄。該目錄早期使用了多個震級標度， 但是近20a以來震級標度主要為mb和MW， 其中4級及以下地震主要為mb震級， 5級以上地震主要為MW震級。本文對該目錄中的所有地震統一用震級M表示。

采用定量評估的 “最大曲率MAXC”方法(Wiemeretal.， 2002)分析全球地震目錄的最小完備性震級Mc隨時間的變化。最大曲率方法(maximum curvature method)是將震級-頻度曲線一階導數的最大值對應的震級作為Mc。結果顯示(圖 1)： 1973—1994年， 最小完備性震級為4.6～5.0級； 1995—2012年， 最小完備性震級為4.3～4.7級； 2013—2016年， 最小完備性震級為4.2～4.6級。

圖 1 基于MAXC方法計算的全球最小完備性震級隨時間的變化曲線Fig. 1 Curve of global catalogue minimum completeness magnitude with time calculated by “MAXC” method.

2 前震序列的活動特征

1976—2017年全球共發生獨立的MW≥8.0淺源地震29次， 其中逆沖型23次。為了確定前震， 本文先做如下約束： 1)前震序列發生在余震區內， 本文取主震后10d中余震的密集活動區作為余震區； 2)主震前30d， 余震區內的地震活動水平明顯高于背景狀態， 頻度達到背景值的2倍以上； 3)最大前震震級M≥5.0。

依據上述約定， 共挑選出8個前-主-余型地震序列(蔣海昆等， 2007)， 其主震發震斷層滑動角均位于69°～110°范圍內， 即全部為逆沖型地震， 占8級以上逆沖型地震總數(23次)的34.8%。前-主-余型地震主要分布在環太平洋地震帶的西北邊界帶、 西南邊界與東南邊界帶(圖 2)， 這主要受全球8級以上地震空間分布的影響。

圖 2 全球1976—2017年MW≥8.0地震的震源機制解Fig. 2 Global shallow focus mechanism solution with MW≥8.0.紅色震源球為有前震的地震

之后統計了 8 次大地震前震活動開始的時間、 最大前震與主震的震級差及震中距等參數， 結果顯示(表1)， 最大前震與主震的震級差為1.1～2.8級， 震中距為10～53km。87.5%的時間差(最大前震與主震發生的時間差)為2h～15d， 僅2006年11月15日千島群島8.3級地震震前45d發生最大前震。

表1 8次大地震的主震和最大前震的基本參數、 最大前震與主震的時間差、 震級差等相關統計Table1 Basic parameters of the 8 great earthquakes and their largest foreshocks， the difference of occurrence time and magnitude between the main shock and the largest foreshock

圖 3 是8次前-主-余型地震序列的震中分布， 其中余震序列的持續時間為主震后1個月。圖 4 是地震序列沿圖 3 所示AB剖面的震中遷移圖， 即將地震投影到AB線段上， 繪制每個地震投影點與A點的距離D隨時間t的變化情況。可以看出， 前震的空間分布集中， 且集中分布于主震震中附近。圖 5 是8次大地震的前震、 主震和余震震源機制解分布圖。 可見， 前震序列震源機制解與主震震源機制解顯示的斷層錯動類型及2個節面的走向、 傾角都較為一致， 而余震的震源機制解比較復雜， 不僅破裂類型不一致， 且節面的走向、 傾角都有很大差異， 這與陳颙(1978)提出的前震特征相符。

圖 3 8次大地震的前震和余震序列震中分布Fig. 3 Epicenter distribution of foreshock and aftershock sequences of 8 great earthquakes.a 1976-01-14， 克馬德克群島海域， M8.2； b 1985-03-03， 智利中部近海， M8.0； c 1986-05-07， 阿留申群島， M8.0； d2006-11-15， 千島群島， M8.3； e 2007-04-01， 所羅門群島， M8.1； f 2011-03-11， 日本本州近海， M9.1； g 2013-02-06， 圣克魯斯群島， M8.0； h 2014-04-01， 智利北部近海， M8.2。紅色圓形、 藍色圓形和黃色五角星分別表示前震、 余震和主震震中； 圖中虛線框表示表2中背景b值的計算范圍

圖 4 8次大地震的前震和余震序列沿AB剖面(圖 3所示)的震中遷移D-t圖Fig. 4 Epicenter migration map in AB direction.a 1976-01-14， 克馬德克群島海域， M8.2； b 1985-03-03， 智利中部近海， M8.0； c 1986-05-07， 阿留申群島， M8.0； d2006-11-15， 千島群島， M8.3； e 2007-04-01， 所羅門群島， M8.1； f 2011-03-11， 日本本州近海， M9.1； 2013-02-06， 圣克魯斯群島， M8.0； h 2014-04-01， 智利北部近海， M8.2。橫坐標為時間t， 以第1次前震發生的時間為零點， 后續地震距離第1次前震的時間間隔為橫坐標值， 單位為d； 圖中紅色空心圓表示主震

圖 5 8次大地震的前震和余震的震源機制解示意圖Fig. 5 Sketch map of focal mechanism solution for earthquake sequence.a 1976-01-14， 克馬德克群島海域， M8.2； b1985-03-03， 智利中部近海， M8.0； c1986-05-07， 阿留申群島， M8.0； d2006-11-15， 千島群島， M8.3； e2007-04-01， 所羅門群島， M8.1； f 2011-03-11， 日本本州近海， M9.1； g 2013-02-06， 圣克魯斯群島， M8.0； h 2014-04-01， 智利北部近海， M8.2。紅色震源球為前震和主震， 藍色震源球為余震

為了考察前震序列在臨近主震時是否出現加速活動， 統計了主震前半年(180d)余震區內M≥Mc(Mc為最小完備震級)地震累積頻度隨時間的變化。取余震區及其附近主震前3a左右的4級以上地震， 使用G-R關系確定最小完備震級。結果顯示， 1976年克馬德克群島海域8.2級地震余震區附近的最小完備震級為4.9級， 1985年智利中部近海8.0級地震附近的最小完備震級為4.7級， 1986年阿留申群島8.0級地震附近的最小完備震級為4.6級， 2007年4月1日所羅門群島8.1級地震附近的最小完備震級為4.2級。另外4次地震， 即2006年千島群島8.3級、 2011年日本本州以東海域9.1級、 2013年圣克魯斯群島8.0級和2014年智利北部近海8.2級， 最小完備震級分別為4.2級、 4.5級、 4.4級和4.0級。圖6a是8次大地震前半年余震區內M≥Mc地震累積頻度隨時間的變化曲線， 圖6b是累積頻度歸一化的結果。可以看出， 主震發生前半年， 余震區內僅有少量地震活動， 而到主震前約1個月地震活動顯著增強， 累計頻度出現加速增長的特點。具體而言， 在8次大地震中， 1986年8.0級、 2006年8.3級和2007年8.1級地震的前震序列分別自主震前35d、 45d和43d開始頻度顯著增加， 其余5次均在主震前15d內開始加速活動。值得注意的是， 1986年和2007年2個前震序列在主震前1d內、 2006年前震在主震前6d又再次出現加速活動(圖 6)。

圖 6 8次大地震主震前180d余震區內最小完備震級以上(M≥Mc)地震的累積頻度曲線Fig. 6 Cumulative frequency curve of earthquakes with M≥Mc in aftershock area 180 days before 8 main earthquakes.橫坐標以主震發生時刻為0點， 負值表示主震前

3 應用ETAS模型計算地震序列參數

本文使用ETAS模型分別計算前震和余震序列參數。ETAS模型假設所有余震均可以按照大森-宇津公式(Utsu， 1961)激發自身的余震， 且震級的分布是獨立的。假定主震的發生時間為初始零時刻， 在其后某觀測時段[0，T]內地震序列{(ti，Mi)；i=1， 2， …，N}的強度函數可表示為(Ogata， 1988)

(1)

其中，t為主震發生后的離逝時間；M0為計算所用序列的截止震級；Mi、ti分別為第i個事件的震級與距主震的時間長度；μ為背景地震發生率；p表示序列衰減快慢；c為主震后余震頻次達到峰值的時間長度；K表示余震的活躍程度；α表示觸發次級余震的能力(Ogata， 1989， 1992)。

ETAS模型參數可使用最大似然法進行估計。在擬合時段[S，T]內， 似然函數L可表述為

(2)

將式(1)帶入式(2)， 即可對參數[μ，K，c，α，p]進行最大似然估計。

為了考察ETAS模型的擬合效果， 一般使用 “殘差分析”方法(Daleyetal.， 2003)將地震序列轉換為在 “轉換時間”(τ)域的分布， 并考察實際地震序列與理論值的擬合情況。將條件強度函數λ(t)采用時間序列{ti}進行如下轉換：

(3)

由此將{ti}轉化為服從單位速率的穩態泊松分布(3)Zhuang J， Harte D， Werner M J， et al.， 2012， Basic models of seismicity: Temporal models， Community Online Resource for Statistical Analysis， doi: 10.5078/corssa-79905851。{τi}。如果ETAS模型對數據的擬合較好， 則在轉換時間域{τi}的地震累積曲線表現為線性， 接近標準穩態泊松過程的理論直線。

此外， 本文利用最大似然法(Aki， 1965； Utsu， 1966)計算b值：

(4)

(5)

式中，n為參與計算的地震數目。

本文采用震級-序號法確定前震和余震的最小完備震級(Huang， 2006)。圖 7 是4次大地震的前震和余震的震級-序號圖。可以看出， 在主震發生后短時間內， 由于主震波形振幅大、 面波持續時間長， 大量小震級余震被 “淹沒”， 使得地震監測能力降低， 余震最小完備震級上升(Iwata， 2008)。本文通過調節ETAS模型擬合的起始時間C0以確保有足夠多的地震參與計算， 即主震后—C0時刻期間的數據不參與計算， 取C0時刻后最小完備震級以上的地震計算序列參數(蔣長勝等， 2013)。由圖 7 可見， 2006年千島群島8.3級地震的前震序列最小完備震級為4.2級(圖7a)， 主震后0.139d余震的最小完備震級為4.4級(圖7b)； 2011年3月11日日本本州近海9.1級地震前震序列的最小完備震級為4.5級(圖7c)， 余震序列在主震后0.529d的最小完備震級是4.5級(圖7d)； 2013年2月6日圣克魯斯群島8.0級地震前震序列的最小完備震級為4.4級(圖7e)， 余震序列在主震發生后0.164d的最小完備震級為4.3級(圖7f)； 2014年4月1日智利北部海域8.2級地震前震序列的最小完備震級為4.0級(圖7g)， 余震序列在距離主震后0.046d的最小完備震級為4.4級(圖7h)。

圖 7 4次大地震的前震和余震序列的震級-序號圖Fig. 7 The magnitude-number chart of foreshocks and aftershocks of 4 great earthquakes.a和b為2006年11月15日千島群島8.3級地震的前震和余震； c和d為2011年3月11日日本本州近海9.1級地震的前震和余震； e和f為2013年2月6日圣克魯斯群島8.0級地震前震和余震； g和h為2014年4月1日智利北部近海8.2級地震的前震和余震。 圖中2條垂直虛線分別表示最小完備震級Mc和達到最小完備震級的時刻C0， C0為距離主震發生的時間， 單位為d

樣本量N和最低計算震級Mj影響地震序列參數的估算結果(Shietal.， 1982； Bender， 1983； 蔣長勝等， 2013； 吳果等， 2019)。為了確保參數計算的可靠性， 本文約定當前震序列滿足計算樣本量N≥30且最低計算震級Mj≥Mc(Mc為最小完備震級)時， 才同時計算前震和余震序列參數α、p和b值。在8個前震-主震型地震序列中， 有4個可以對比計算其前震和余震序列參數， 其中余震序列的持續時間為主震后1個月(表2)。針對2006年千島群島8.3級和2014年智利8.2級地震的前震序列， 我們分別選了2個最低計算震級， 結果顯示， 不同的最低計算震級可以得到不同的α、p和b值， 但由此造成的差異與計算誤差大體相當， 即ETAS模型的計算結果比較穩定。為了比較， 對同一個地震的前震和余震取相同或相近的最低計算震級。所得結果表明， 前震和余震序列的α值和p值沒有規律性差異， 而b值規律性明顯， 前震b值顯著低于余震(圖 8， 表2)。

為了進一步識別前震， 我們計算余震區及附近較長時間的區域背景b值(確定計算時間主要考慮完備震級且避免期間發生大地震， 計算區域為圖 3 虛線框所示區域)。經對比發現， 4次前震序列的b值均低于區域背景b值， 而3次的余震b值高于區域背景b值， 1次的余震b值與背景b值相當(表2)。計算4次前震序列b值與區域背景b值的差， 再除以區域背景b值， 結果依次為0.243、 0.163、 0.138和0.103， 即前震序列的b值低于區域背景b值的10%～24%。前震b值與區域背景b值之差依次是區域背景b值標準差的4.2倍、 6.1倍、 7.1倍和2.2倍， 具有顯著性。

為了討論前震b值的穩定性， 對樣本量豐富的2個前震序列， 即2006年千島群島8.3級和2014年智利北部8.2級地震的前震序列進行累積滑動計算， 最低計算震級分別為4.2級和4.0級。結果顯示(圖 9)， 在前震發生的最初階段b值較低， 之后逐漸增大， 當計算樣本量N≥70后，b值基本穩定。

圖 8 4次大地震的前震和余震b值(a)、 α值(b)和p值(c)計算結果對比Fig. 8 Parameters b-value(a)， α-value(b)and p-value(c)of foreshocks and aftershocks for 4 great earthquakes.圖件橫坐標為地震編號， 與表2 序號一致

圖 9 2006年千島群島8.3級(a)和2014年智利北部近海8.2級(b)地震前震序列的b值隨計算樣本量的變化Fig. 9 Change of b-value with the calculated sample size of foreshock sequences of the 2006-11-15 Kuril Islands M8.3 and the 2014-04-01 off the coast of northern Chile M8.2 earthquake.圖中直線和虛線表示背景b值及其誤差

4 結論與討論

綜上所述， 得到以下幾點結論。

(1)1976—2017年， 全球共發生29次8級以上淺源地震， 其中8次有前震， 占總數的27.6%， 略低于1900—2011年全球8.5級以上地震的前震比例(Xueetal.， 2012)， 明顯高于中國大陸地區中強以上地震的前震比例(陳颙等， 2015)。

(2)8次前震序列的主震發震斷層的滑動角為69°～110°， 即全部為逆沖型， 占8級以上逆沖型地震總數(23次)的34.8%。

(3)前震序列空間分布集中， 且集中分布在主震震中附近。最大前震與主震的震中距為10～53km， 震級差為1.1～2.8級， 優勢發生時間差為2h～15d， 僅2006年11月15日千島群島8.3級地震震前45d發生最大前震。

(4)主震前半年， 余震區內有少量的地震活動， 而到主震前約1個月， 特別是主震前15d內， 地震活動顯著增強， 累計頻度加速增長， 即前震序列具有高頻度活動的特點。具體而言， 8次大地震中有5次在主震前15d內出現加速活動； 3次的加速活動開始時間超過1個月(35～45d)， 但在主震前1d和6d又再次出現頻次顯著增多的現象。

(5)前震的震源機制解一致性好， 且與主震一致， 而余震的震源機制解比較復雜。

(6)應用ETAS模型計算前震和余震序列的α、p和b值， 結果顯示， 對于反映激發次級余震能力的α值和序列衰減快慢的p值， 前震和余震沒有共性差異； 而b值則共性特征明顯， 前震b值顯著低于余震。與區域背景b值相比， 前震b值明顯偏低； 而余震b值高于背景b值或與其相當。

(7)為了討論前震b值的穩定性， 計算了2個資料豐富的前震序列b值隨樣本量的變化。結果顯示， 在前震序列的開始階段b值較低， 之后逐漸增大， 當計算樣本量N≥70后，b值基本穩定。

雖然前震是短臨地震預報最有效的方法之一， 但前震的識別至今仍較為困難。特別是前震與前兆震群(宋俊高等， 1989)更為相似， 兩者均具有空間分布集中、 高頻次活動和震源機制解一致的特點， 如乳山震群(劉方斌等， 2018)和蓋州震群(王亮等， 2015)。但前兆震群與余震均具有高b值的特點(前兆震群b≥0.65)(宋俊高等， 1989)， 而前震b值則較低(Guliaetal.， 2019)。

b值不僅反映了大、 小地震的比例關系， 更有明確的物理意義。b值的大小主要取決于震源區應力狀態和介質性質， 并與應力呈反比(Scholz， 1968； Wyss， 1973； Urbancicetal.， 1992)， 可用來衡量區域相對應力水平(Wyssetal.， 2000； 易桂喜等， 2013)。本文對比研究了4次大地震的前震和余震序列的b值， 考慮了最低計算震級和樣本量的影響， 并引入區域背景b值， 將前震b值與其進行對比， 結果顯示， 前震b值低于區域背景b值的10%～24%， 前震b值與區域背景b值之差是區域背景b值標準差的2.2～7.1倍， 具有顯著性， 這與Gulia等(2019)的研究結果一致。

此外， 本文研究的震例均為板間地震。有研究顯示(蘇有錦等， 2014)， 對于淺源地震， 板內和板間地震序列具有相似的時、 空、 強統計特征。 因此， 本文的研究結果對板內大地震具有參考意義。