四川地區(qū)M5.0及以上地震時間間隔統(tǒng)計分析與概率預(yù)測

王 璇

（上海對外經(jīng)貿(mào)大學(xué)統(tǒng)計與信息學(xué)院，上海201620）

地震發(fā)生時間間隔的分布特征一直以來都是地震學(xué)家感興趣的話題，了解區(qū)域地震發(fā)生時間間隔對地震區(qū)劃、地震預(yù)報、地震危險性分析和地震災(zāi)害預(yù)測具有重要意義。中國位于環(huán)太平洋地震帶和歐亞地震帶的交匯部位，是一個震災(zāi)嚴(yán)重的國家，地震活動頻度高、強度大、震源淺、分布廣。基于地質(zhì)的地震研究的作用不可否認(rèn)，通過運用沉積地質(zhì)學(xué)、地貌學(xué)和構(gòu)造地質(zhì)學(xué)等常用方法有助于評估未來地震發(fā)生概率（劉靜等，2021），這些方法的局限性在于許多地震不會在確定的斷層出現(xiàn)，因此，將統(tǒng)計學(xué)應(yīng)用到地震發(fā)生和預(yù)測上很有意義。本文將地震發(fā)生看作一個統(tǒng)計過程，假設(shè)地震發(fā)生時間間隔是與一些概率分布模型相關(guān)聯(lián)的隨機變量，找出地震發(fā)生時間間隔的最佳擬合概率分布函數(shù)，它的分布特征對研究地震活動的復(fù)發(fā)和預(yù)測未來的地震有重要作用。

自從概率統(tǒng)計方法與地震預(yù)報研究相結(jié)合以來，許多學(xué)者對地震發(fā)生時間間隔的統(tǒng)計規(guī)律進(jìn)行了不同程度的研究。于利民等（1996）對中國大陸有史以來的強地震發(fā)生時間間隔進(jìn)行分析，得出強震時間間隔的經(jīng)驗概率和模擬函數(shù)。范琦（2001）以甘肅天祝地區(qū)為例，利用最大熵原理進(jìn)行分析檢驗，結(jié)果認(rèn)為地震震級和地震間隔時間的變化服從威布爾分布。李強和徐桂明（2002）利用威布爾分布建立了江蘇及鄰區(qū)中強地震的時間間隔分布模型，并對該地區(qū)未來幾年發(fā)生中強地震的概率進(jìn)行預(yù)測。王煒和劉震華（1987）探討了華北地區(qū)小震時間間隔的統(tǒng)計分布，經(jīng)過統(tǒng)計檢驗認(rèn)為它服從威布爾分布。張亞奎（1996）對吉林省伊通—舒蘭斷裂帶以及我國南北地震帶南段安寧河—滇東地震帶發(fā)生的地震進(jìn)行了統(tǒng)計分析，發(fā)現(xiàn)在大的地震帶上發(fā)震時間間隔服從正態(tài)分布。Utsu（1984）使用威布爾分布、伽馬分布、對數(shù)正態(tài)分布和指數(shù)分布估計日本地區(qū)的地震發(fā)生時間間隔。Sadeghian（2018）研究適用于伊朗地區(qū)的地震時間間隔數(shù)據(jù)的概率密度函數(shù)，結(jié)果為伽馬分布、貝塔分布、指數(shù)分布、對數(shù)正態(tài)分布、正態(tài)分布和威布爾分布。Khan等（2021）根據(jù)擬合優(yōu)度檢驗，得出威布爾分布和對數(shù)邏輯分布是興都庫什地區(qū)最佳擬合模型，對數(shù)正態(tài)分布和伽馬分布分別為該地區(qū)第二和第三合適的概率模型。Somette和Knopoff（1997）基于威布爾分布、對數(shù)正態(tài)分布、指數(shù)分布、伽馬分布和貝塔分布等5種概率分布研究了地震時間間隔并預(yù)測地震復(fù)發(fā)時間。

2008年5月12日，四川汶川發(fā)生8.0級大地震，破壞性極強、波及范圍廣，造成了非常嚴(yán)重的人員傷亡和經(jīng)濟(jì)損失。在四川發(fā)生過的中小地震更是數(shù)不勝數(shù)，因此本文選取中國四川地區(qū)發(fā)生的中強地震作為研究對象，從指數(shù)分布、伽馬分布、對數(shù)正態(tài)分布和威布爾分布4種概率分布出發(fā)，分別研究1970年至2007年和1970年至今發(fā)生的M5.0及以上地震，通過擬合優(yōu)度檢驗找出地震發(fā)生時間間隔的最佳擬合分布，并評估未來發(fā)生M5.0及以上地震的概率。

1 數(shù)據(jù)與方法

1.1 地震數(shù)據(jù)

根據(jù)國家地震科學(xué)數(shù)據(jù)中心記載（1970年1月1日至2021年7月31日的地震數(shù)據(jù)），四川自1970年至2021年7月曾發(fā)生過164次M5.0及以上的地震。由于2008年汶川地震余震過于頻繁，若使用2008年全部地震數(shù)據(jù)預(yù)測未來發(fā)生強震則會導(dǎo)致概率值偏高，因此剔除2008年發(fā)生的余震數(shù)據(jù)。本文使用的數(shù)據(jù)見表1，共121條數(shù)據(jù)。

表1 地震數(shù)據(jù)

本文主要探討兩部分內(nèi)容：（1）使用1970年至2007年的地震數(shù)據(jù)，計算兩次地震時間間隔（以天為單位），根據(jù)計算結(jié)果做模型擬合和預(yù)測，通過預(yù)測2008年發(fā)生M5.0及以上地震的可能性從而驗證模型選擇方法的可靠性；（2）使用表1全部數(shù)據(jù)用同樣的方法預(yù)測未來發(fā)生M5.0及以上地震的可能性。

1.2 模型介紹

本文分析4種概率分布，首先介紹這4種概率分布的概率密度函數(shù)、模型參數(shù)（見表2）。

表2 4種概率分布的概率密度函數(shù)

（1）指數(shù)分布：指數(shù)分布是伽馬分布和威布爾分布的特殊形式，具有無記憶性，即對隨機變量T，有P（T＞s＋t｜T＞t）＝P（T＞s），且（t，s＞0），常用來表示獨立隨機事件發(fā)生的時間間隔，其參數(shù)λ表示單位時間內(nèi)事件發(fā)生的平均次數(shù)。其在可靠性研究具有廣泛應(yīng)用。

（2）伽馬分布：伽馬分布的參數(shù)α為形狀參數(shù)，β為尺度參數(shù)，當(dāng)α＝1時，伽馬分布為指數(shù)分布。龔平等（2001）指出當(dāng)?shù)卣鸬陌l(fā)生為相互獨立事件且僅依賴于時間間隔，與時間起點無關(guān)時，可得出地震發(fā)生時間服從伽馬分布。

（3）對數(shù)正態(tài)分布：對數(shù)正態(tài)分布的參數(shù)μ為位置參數(shù)，σ為形狀參數(shù)，與正態(tài)分布相似，在可靠性研究中，數(shù)據(jù)若不符合正態(tài)分布，則常取其對數(shù)使之符合正態(tài)分布，因此稱作對數(shù)正態(tài)分布。

（4）威布爾分布：威布爾分布的參數(shù)η為比例參數(shù)，m為形狀參數(shù)，當(dāng)m＝1時，威布爾分布為指數(shù)分布。該分布在地震研究中也具有廣泛的應(yīng)用。

1.3 驗證模型

本文使用常用的模型選擇方法——赤池信息準(zhǔn)則（AIC準(zhǔn)則）、貝葉斯信息準(zhǔn)則（BIC準(zhǔn)則）和K-S檢驗。

AIC可以衡量統(tǒng)計模型擬合的優(yōu)良性，它建立在熵的概念上，提供了權(quán)衡估計模型復(fù)雜度和擬合數(shù)據(jù)優(yōu)良性的標(biāo)準(zhǔn)。AIC定義為式（1）：

式中k為模型參數(shù)個數(shù)，L為似然函數(shù)。從一組可供選擇的模型中選擇最佳模型時，通常選擇AIC最小的模型。AIC為模型選擇提供了有效的規(guī)則，但也有不足之處。當(dāng)樣本容量很大時，在AIC準(zhǔn)則中擬合誤差提供的信息就要受到樣本容量的放大，而參數(shù)個數(shù)的懲罰因子卻和樣本容量沒關(guān)系，因此當(dāng)樣本容量很大時，使用AIC準(zhǔn)則選擇的模型不收斂于真實模型，它通常比真實模型所含的未知參數(shù)個數(shù)要多。

BIC貝葉斯信息準(zhǔn)則與AIC相似，可用于模型選擇，并且改進(jìn)了AIC的不足之處。BIC定義為式（2）：

式中k為模型參數(shù)個數(shù)，n為樣本數(shù)量，L為似然函數(shù)。BIC值越小，說明模型擬合程度越高。

單樣本K-S檢驗是一種擬合優(yōu)度的非參數(shù)檢驗方法，利用樣本數(shù)據(jù)推斷總體是否服從某一理論分布。K-S檢驗一般返回兩個值：D和P值。其中D值表示兩條累計分布曲線之間的最大垂直距離，所以D值越小，這兩個分布的差距越小，分布越一致。P值是假設(shè)檢驗里面的P值，如果P值大于0.05，那么就不能拒絕原假設(shè)，所以P值越大，分布越一致。

2 結(jié)果分析

本文首先擬合四川省2008年以前發(fā)生M5.0及以上地震時間間隔的概率分布模型，計算參數(shù)值，并根據(jù)AIC、BIC、D值和P值來檢驗指數(shù)分布、伽馬分布、對數(shù)正態(tài)分布和威布爾分布對地震發(fā)生時間間隔的擬合程度，運算結(jié)果如表3所示。根據(jù)數(shù)據(jù)的累積概率對應(yīng)4種分布累積概率繪制p-p plot（見圖1），可以直觀的看出樣本數(shù)據(jù)是否服從某一分布。由于指數(shù)分布的AIC和BIC值相較于其他3種分布偏大很多，說明擬合效果最差。比較伽馬分布、對數(shù)正態(tài)分布和威布爾分布的K-S距離，可以看到在K-S檢驗中伽馬分布返回的D值最小，P值最大。由圖1可見，樣本數(shù)據(jù)與伽馬分布擬合的較好。因此，選用伽馬分布作為四川省2008年以前發(fā)生M5.0以上地震時間間隔的最佳擬合分布。

圖1 1970年至2007年地震數(shù)據(jù)4種概率分布P-P plot圖

表3 1970年至2007年地震數(shù)據(jù)建模的擬合檢驗結(jié)果

2008年汶川8.0級地震發(fā)生之后，大小余震不斷，使用2008年的全部地震數(shù)據(jù)預(yù)測未來并不合理，因此使用剔除2008年余震數(shù)據(jù)以后的地震數(shù)據(jù)對四川發(fā)生M5.0及以上地震時間間隔的概率分布模型進(jìn)行擬合檢驗，運算結(jié)果如表4所示。根據(jù)數(shù)據(jù)的累積概率對應(yīng)4種分布累積概率繪制p-p plot（見圖2）。指數(shù)分布和對數(shù)正態(tài)分布未通過擬合檢驗，伽馬分布的AIC、BIC和D值均最小，P值最大，并且通過圖2可以直觀的看出樣本數(shù)據(jù)與伽馬分布更加接近，因此，可以認(rèn)為伽馬分布對四川省1973年以來M5.0及以上地震時間間隔數(shù)據(jù)的擬合程度最高。

表4 1970年至今地震數(shù)據(jù)（剔除2008年余震數(shù)據(jù)）建模的擬合檢驗結(jié)果

圖2 1970年至今地震數(shù)據(jù)（不含2008年余震）4種概率分布P-P plot圖

利用建立的伽馬分布對四川的地震危險性進(jìn)行預(yù)測。根據(jù)四川2008年以前發(fā)生M5.0及以上地震數(shù)據(jù)求得的伽馬分布累積函數(shù)（見圖3），其橫坐標(biāo)表示地震發(fā)生的時間間隔，單位為天；縱坐標(biāo)表示某一時間地震發(fā)生的概率。將上一次發(fā)生地震的時間作為坐標(biāo)原點，則從圖中可得出自上一次地震發(fā)生之后再發(fā)生一次地震的累積概率。由表1可知，在2008年5月12日發(fā)生強震的前一次地震在2005年8月5日，從圖3可以看到，在2005年8月5日之后2000天中，發(fā)生地震的概率隨時間的增加而增大，365天（一年）后發(fā)生M5.0及以上地震的概率為0.87左右；739天（兩年）后發(fā)生M5.0及以上地震的概率為0.95左右；預(yù)測到在1011天（約2.77年）時，即2008年5月12日發(fā)生M5.0及以上地震的概率高達(dá)0.98，這與事實基本吻合，說明伽馬分布的預(yù)測效果較好。

圖3 2005年后發(fā)生M5.0及以上地震概率曲線

根據(jù)上述方法，使用全部樣本數(shù)據(jù)預(yù)測未來發(fā)生M5.0及以上地震的概率（見圖4）。由表1可知，截至2021年7月31日，四川發(fā)生上一次強震日期為2020年4月1日，在圖4中，設(shè)2020年4月1日為坐標(biāo)原點，可以看到，639天后（即2021年年底）發(fā)生M5.0及以上地震的概率大約為0.94，1 004天后（即2022年底）發(fā)生M5.0及以上地震的概率大約為0.98，1 369天后（即2023年底）發(fā)生M5.0及以上地震的概率為0.99。

圖4 未來發(fā)生M5.0及以上地震概率曲線

3 結(jié)論

本文討論了1970年至2007年和1970年至今四川地區(qū)發(fā)生M5.0及以上地震時間間隔的分布及概率預(yù)測，得到如下結(jié)論：（1）經(jīng)過擬合檢驗得出，伽馬分布為4種概率分布（指數(shù)分布、伽馬分布、對數(shù)正態(tài)分布和威布爾分布）中對四川地區(qū)M5.0及以上地震發(fā)生時間間隔的最佳擬合。（2）根據(jù)1970年至2007年四川發(fā)生M5.0及以上地震數(shù)據(jù)，通過伽馬分布預(yù)測在2008年會發(fā)生M5.0及以上地震，與事實相符，證明模型選取與預(yù)測具有可靠性。（3）根據(jù)1970年1月1日至2021年7月31日四川地區(qū)發(fā)生M5.0及以上地震數(shù)據(jù)（剔除2008年汶川地震余震），利用伽馬分布預(yù)測未來可能發(fā)生地震的概率，得到四川地區(qū)2021年發(fā)生M5.0及以上地震的概率約0.94，2022年發(fā)生M5.0及以上地震的概率約0.98，2023年以后發(fā)生M5.0及以上地震的概率趨于1，說明發(fā)生M5.0及以上的地震的概率極高。（4）確定最佳的概率模型可以為研究地震提供另一種評估該地區(qū)地震危險性的方法。但本研究存在的不足之處在于僅從4種概率分布出發(fā)研究地震事件，更適用于對地震預(yù)測的初步參考，在后續(xù)的研究中可以研究更多概率分布在地震預(yù)測上的應(yīng)用。