基于Gutenberg-Richter定律快速估算最大余震震級：以2017年九寨溝MS 7.0地震為例

解孟雨 孟令媛 申文豪 史保平

1）中國科學院大學地球科學學院，北京市玉泉路19號（甲） 100049

2）中國地震臺網(wǎng)中心，北京 100045

3）中國地震局地殼應力研究所，北京 100085

0 引言

據(jù)中國地震臺網(wǎng)中心測定，2017年8月8日21時19分四川省阿壩州九寨溝縣發(fā)生MS7.0地震（以下簡稱九寨溝地震），震中 33.22°N、103.83°E，震源深度 20km，矩震級MW6.5。地震震中位于岷江斷裂、塔藏斷裂和虎牙斷裂附近，推測其發(fā)震構(gòu)造為塔藏斷裂南側(cè)分支或虎牙斷裂北段。此次地震是繼2008年汶川MS8.0地震和2013年蘆山MS7.0地震后在青藏高原東緣發(fā)生的又一次破壞性地震，據(jù)統(tǒng)計，截至8月13日20時，地震共造成25人死亡，525人受傷，6人失聯(lián)。九寨溝縣17個鄉(xiāng)鎮(zhèn)的1680戶房屋受損，52個風景名勝中的30個受到影響，造成較為嚴重的人員經(jīng)濟損失。

中強地震發(fā)生之后，常伴隨著一系列震級小于主震震級的地震，即余震。主震發(fā)生后為了避免較大余震產(chǎn)生的2次危害，有必要快速準確地預測主震震中周邊余震的時空分布以及可能發(fā)生的最大余震。對于最大余震的震級，B?th（1965）指出，主震震級與最大余震震級之差ΔM是不受主震震級影響的常數(shù)（即貝特定律）

其中，Mmain為主震震級；MRLA為最大余震震級；且ΔM的平均值約為1.2。研究該定律對于理解地震破裂過程和余震發(fā)生規(guī)模有著重要意義（Shcherbakov et al，2005；Shearer，2012；Rodríguez-Pérez et al，2016）。眾多研究者對 ΔM的大小、影響因素及物理機制等進行了研究，結(jié)果顯示，貝特定律的存在表明地震觸發(fā)機制中存在一定程度的自相似過程（Console et al，2003；Shearer，2012；?alohar，2014）；貝特定律并不嚴格適合所有余震序列，常數(shù) ΔM僅在平均意義上接近 1.2（Kisslinger et al，1991；Console et al，2003；Shcherbakov et al，2005、2013；Shearer，2012；Tahir et al，2012），其值為 0～3（Kisslinger et al，1991；Helmstetter et al，2003；Shcherbakov et al，2004、2013；蔣海昆，2010），且可能受震源機制（蘇有錦等，2008；Tahir et al，2012；Rodríguez-Pérez et al，2016）、斷層之間的相互作用情況（?alohar，2014）、震源深度（B?th，1965；蘇有錦等，2014）、地震序列類型（蘇有錦等，2014）、區(qū)域構(gòu)造特征（Rodríguez-Pérez et al，2016）等因素的影響。為了更好地估計最大余震震級，我國學者依據(jù) Gutenberg-Richter定律（以下簡稱G-R定律）首先提出了b值截距法（吳開統(tǒng)等，1984；毛春長，1989；張智等，1989；國家地震局科技監(jiān)測司，1990）。該方法是將震級-頻度圖上 G-R定律的擬合曲線與震級軸交點的震級值作為最大余震震級的估計值。之后，研究者一方面不斷驗證該方法；另一方面，利用該方法估計最大余震震級，同時也計算了在震后不同時間內(nèi)估計結(jié)果的變化情況（劉正榮，1995；錢曉東等，2008；蘇有錦等，2008、2014；付虹等，2008）。研究結(jié)果顯示，雖然使用該方法得到的最大余震震級存在一定的偏差，但考慮到誤差的影響，該方法仍可以很好地估計最大余震震級。在國外，Shcherbakov等（2004）依據(jù)G-R定律提出了改進貝特定律，并用于對美國加州包括 Landers7.3級地震在內(nèi)的10個中強地震余震序列的研究。他們提出了推定最大余震震級MILA的概念，并給出了MILA的計算公式，其實質(zhì)與b值截距法相同，不同的僅是相應的描述方式。與b值截距法一樣，對于不同余震序列，計算得到的推定最大余震震級MILA與真實最大余震震級MRLA也基本一致（Shcherbakov et al，2004、2013）。鑒于推定最大余震震級MILA的含義更為明確，故本文采用推定最大余震震級MILA進行方法說明。本文利用推定最大余震震級MILA的計算公式，針對九寨溝地震余震序列，探討了在不同時間尺度下提取最大余震震級估算值的可行方案，并與實際測量值進行比較。依據(jù)對九寨溝地震余震序列的計算和結(jié)果分析，本文給出了快速評估最大余震震級的方法。

1 數(shù)據(jù)及完備性分析

本文選取九寨溝地震發(fā)震后至10月5日3時18分（主震發(fā)生后約57天）的余震數(shù)據(jù)用于建模和參數(shù)確定，該目錄來源于中國地震臺網(wǎng)中心，其中包含里氏震級ML和面波震級MS兩類數(shù)據(jù)。若僅考慮里氏震級目錄，則共有余震事件5864個，ML≥3.0事件共92個，其中，3.0≤ML＜4.0地震 78次，4.0≤ML＜5.0地震 13次，5.0≤ML＜6.0地震 1次。特別需要強調(diào)的是，在所選取的地震目錄中，最大震級的余震為8月9日10時17分（距主震發(fā)生時刻約0.54天）發(fā)生的ML5.2（或MS4.8）地震。圖1為九寨溝地震余震的空間分布情況。由圖1可見，九寨溝地震震中位于汶川MS8.0地震、蘆山MS7.0地震震中的北側(cè)，且九寨溝地震余震主要集中于主震附近，近似條帶狀分布。在使用地震目錄的數(shù)據(jù)時，首先應確定目錄的完備震級Mc，即在給定的時間、空間范圍內(nèi)可全部測定地震的最小震級。如果完備震級Mc的估算值太高，則會導致地震數(shù)據(jù)使用不充足，不能獲得完整的地震信息；而太低的完備震級Mc值又會引入記錄不完整的數(shù)據(jù)，導致錯誤的模型參數(shù)估計，進而得出錯誤的結(jié)論（Mignan et al，2012；蔣長勝等，2013）。因此，對于九寨溝地震的余震序列，我們需要選擇和確定快速且準確的完備震級Mc估算方法。完備震級的估算方法主要有最大曲率、GFT、EMR、MBS等方法（Hamdache et al，2017），其中，利用最大曲率（以下簡稱 MAXC）方法進行計算，能較為適合快速評估的需要。首先，該方法在實際計算時將非累積頻率-震級分布中取得最大值的震級檔作為完備震級Mc，因而它不需要任何的參數(shù)擬合，是一種快速直接且對數(shù)據(jù)強健的方法（Mignan et al，2012）。此外，Mignan等（2012）還指出，MAXC方法也可在目錄數(shù)據(jù)較少的情況下得出合理的完備震級，適合在主震發(fā)生后短時間內(nèi)估算完備震級Mc的大小。不過MAXC方法會低估地震目錄的完備震級，因此，為更接近實際真值，可以添加相應的修正系數(shù)，即Mc=Mc（MAXC）＋0.2（Woessner et al，2005）。

圖1 九寨溝地震余震分布與周邊近期發(fā)生的強震

因此，本文采用修正的MAXC方法對九寨溝地震的余震目錄進行完備震級的估算。對于估算完備震級的余震目錄，我們選擇主震后2h（約為0.083天）內(nèi)的余震序列目錄。之所以采用這樣的選擇，一方面是為了能夠在主震發(fā)生后很短的時間尺度下快速得到完備震級，另一方面也是由于采用余震序列早期的目錄才可得到較為準確的完備震級（Enescu et al，2002；米琦等，2015）。經(jīng)過對九寨溝地震2類目錄（里氏震級ML、面波震級MS目錄）的計算得到，對于里氏震級ML的余震目錄，其完備震級Mc約為2.4；而震級類型為面波震級MS的余震目錄，其完備震級Mc=1.6。2類目錄完備震級估算值的不確定度均約為0.007，可忽略不計。

2 最大余震震級的估算

無論是全球尺度的地震還是區(qū)域尺度的地震，除去少數(shù)個例，基本都滿足Gutenberg-Richter定律（簡稱 G-R定律）（Marzocchi et al，2003；Shcherbakov et al，2004），即在給定區(qū)域和時間段內(nèi)，其震級大于等于M的地震總個數(shù)N（≥M），滿足

其中，M≥Mc；a、b為正常數(shù)，a給出了選定區(qū)域和時間段內(nèi)0級以上地震個數(shù)的對數(shù)，它依賴于選定區(qū)域和觀測數(shù)據(jù)時間尺度的大?。℉amdache et al，2017），而b值則刻畫了區(qū)域內(nèi)不同震級地震的相對分布情況（Hamdache et al，2017），其值為 0.6～1.4（Marzocchi et al，2003；Shcherbakov et al，2004；Wang et al，2015；Hamdache et al，2017）。b值的大小與區(qū)域應力狀態(tài)和斷層的破裂過程有關(guān)，可用來指示區(qū)域地震危險性水平（吳開統(tǒng)等，1984；Hamdache et al，2017）。

圖2 九寨溝地震發(fā)生后1天內(nèi)余震目錄的震級-頻度圖

對于余震序列而言，G-R定律依舊成立。在這種情況下，N（≥M）表示余震序列中滿足震級大于等于M（M≥Mc）的地震的個數(shù)，而且對于余震序列，相應b的取值范圍也與區(qū)域或全球的統(tǒng)計結(jié)果一致（Shcherbakov et al，2004）。圖2為九寨溝地震發(fā)生后1天內(nèi)余震目錄的震級-頻度圖。由圖2可見，對于實際余震目錄而言，N=1的最大震級檔（M=5.2）即為主震之后最大余震震級MRLA，即N（≥MRLA）=1。據(jù)此，Shcherbakov等（2004）定義推定最大余震震級MILA為用G-R定律描述的余震震級-頻度分布滿足N（≥M）=1時的地震震級檔M（圖 2），即

于是，貝特定律可改寫為

式（4）即為Shcherbakov等（2004）提出的改進貝特定律。Shcherbakov等（2004）利用改進貝特定律分析了加州1987～2003年間發(fā)生的10個震級大于5.5級地震的最大余震震級，結(jié)果顯示，由該定律得到的推定最大余震震級MILA或震級差（ΔM′）與實際值一致。因此，本文中我們使用推定最大余震震級MILA（式（3））估計實際最大余震震級MRLA。

為了計算推定最大余震震級MILA，首先要估算出b、a值。對于b值的估計，可采用最小二乘法和最大似然法。前人的研究指出，利用最小二乘法計算的b值存在系統(tǒng)性偏差（Sandri et al，2007；Hamdache et al，2017）。因此，本文選擇 Marzocchi等（2003）描述的最大似然法估計b值。具體而言，b值及不確定度δb（Uncertainty）的計算公式為

而相應的不確定度δa為（Taylor，1997）（完備震級Mc的不確定度約為 0.007；相對于b值的不確定度可忽略不計）

于是，推定最大余震震級MILA的不確定度δMILA可表示為（Taylor，1997）

3 結(jié)果分析與討論

針對九寨溝地震的余震序列，我們分別對里氏震級ML目錄和面波震級MS目錄計算了不同時間尺度下推定最大余震震級MILA。圖3、圖4、表1和表2為相應的計算結(jié)果。由圖3（a）、3（c）、4及表1、2可見，隨著選取余震目錄時間范圍（即從主震發(fā)生后到截止時刻經(jīng)過的時間）的不斷擴大，推定最大余震震級MILA逐漸增大并穩(wěn)定于固定值。對于里氏震級ML目錄而言，該值為MILA=5.3±0.1；而對于面波震級MS，則為MILA=5.0±0.1，這與實際測得的最大余震震級基本一致（MRLA=ML5.2或MRLA=MS4.8）。這也說明利用推定最大余震震級MILA可對最大余震震級進行有效估算，并且不受地震目錄中震級類型的限制，即可獲得與實際觀測較一致的數(shù)值。而MILA的不確定度隨時間的流逝不斷減少，這是由于隨時間的增加，地震目錄不斷增加，有效信息不斷增多。這也提醒我們，如果能加強臺網(wǎng)對余震的識別能力，增加余震數(shù)據(jù)，就可以給出更為準確的推定最大余震震級。

圖3 不同時間尺度下由九寨溝地震序列得到的推定最大余震震級MILA隨時間變化的誤差棒圖

對于較小的截止時間，即較短的余震目錄時間，圖3（b）、3（d）顯示，在0.6天（14.4h）之前推定最大余震震級MILA并不穩(wěn)定，并且整體上不確定性較大（余震個數(shù)較少，信息不足）。不過，隨著截止時間的增加，MILA的值趨于穩(wěn)定。截止時間為1、2天的MILA值分別為5.0±0.3（ML）、4.6±0.3（MS）和 5.2±0.3（ML）、4.9±0.2（MS），這些值與實際測得的最大余震震級MRLA（ML5.2和MS4.8）基本一致。因此，在短時間內(nèi)，仍可以利用推定最大余震震級MILA估算最大余震震級MRLA，而且也依然不受震級類型的影響。此外，在實際最大震級余震發(fā)生前（約為0.53天前），由于地震目錄中余震個數(shù)較少，利用0.5天的數(shù)據(jù)得到的推定最大余震震級ML（4.9±0.4）和MS（4.5±0.3）均小于實際最大余震震級MRLA（MRLA=ML5.2或MRLA=MS4.8），不過它們卻可以作為MRLA值的合理下界，這也說明利用推定最大余震震級確實可以有效地預測最大余震震級。據(jù)此，在主震發(fā)生后利用推定最大余震震級，我們可以在1天或2天時得到較為合理的最大余震里氏震級或面波震級的估計值；另一方面，在更短的時間內(nèi)（0.5天或12h），對于現(xiàn)有臺網(wǎng)識別得到的余震數(shù)據(jù)，利用推定最大余震震級可以給出MRLA合適的下界。隨著臺網(wǎng)識別能力的提升，一方面可以增加各個震級余震的個數(shù)N，降低不同參數(shù)的不確定度，進而得到更為準確的MRLA估計值；另一方面也會降低完備震級Mc，在主震后短時間內(nèi)就可有大量可使用的余震數(shù)據(jù)，實現(xiàn)在更短時間內(nèi)得出MRLA的有效估計值。

圖4 九寨溝地震發(fā)生后不同時間段內(nèi)余震目錄的震級-頻度圖

表1 不同時間尺度下計算得到的九寨溝地震推定最大余震里氏震級MILA（ML）

表2 不同時間尺度下計算得到的九寨溝地震推定最大余震面波震級MILA（MS）

4 結(jié)論

本文以2017年九寨溝MS7.0（或ML7.1）地震震后57天的余震目錄（包含ML、MS目錄）為資料，通過對余震目錄完備震級、a值和b值的計算，并利用推定最大余震震級，在不同時間尺度上估算了九寨溝地震的最大余震震級。結(jié)果顯示，推定最大余震震級可有效地預測最大余震震級，并給出合理的估計值。具體而言，在長時間尺度上推定最大余震震級可以給出與實際觀測最大余震震級一致的估計值；而在短時間尺度上，1、2天的結(jié)果也可作為最大余震震級的合適估計，且0.5天（12h）的結(jié)果就可作為最大余震震級的合理下界。

此外，利用推定最大余震震級估計實際最大余震震級的方法，不受震級類型的限制，對于MS、ML目錄均可得到較好的結(jié)果。因此，我們建議可采用以下步驟進行最大余震震級的快速評估：①根據(jù)震后2h的目錄，利用修正的MAXC方法計算完備震級；②利用推定最大余震震級分別在主震后12h、1天、2天時估算最大余震震級，將其分別作為最大余震震級的下界和估計值。另外，計算結(jié)果也顯示，如果可以提高余震的識別能力，利用推定最大余震震級則可在短時間內(nèi)給出更準確的最大余震震級的估計值。