王 璇
(上海對外經貿大學統計與信息學院,上海201620)
地震發生時間間隔的分布特征一直以來都是地震學家感興趣的話題,了解區域地震發生時間間隔對地震區劃、地震預報、地震危險性分析和地震災害預測具有重要意義。中國位于環太平洋地震帶和歐亞地震帶的交匯部位,是一個震災嚴重的國家,地震活動頻度高、強度大、震源淺、分布廣。基于地質的地震研究的作用不可否認,通過運用沉積地質學、地貌學和構造地質學等常用方法有助于評估未來地震發生概率(劉靜等,2021),這些方法的局限性在于許多地震不會在確定的斷層出現,因此,將統計學應用到地震發生和預測上很有意義。本文將地震發生看作一個統計過程,假設地震發生時間間隔是與一些概率分布模型相關聯的隨機變量,找出地震發生時間間隔的最佳擬合概率分布函數,它的分布特征對研究地震活動的復發和預測未來的地震有重要作用。
自從概率統計方法與地震預報研究相結合以來,許多學者對地震發生時間間隔的統計規律進行了不同程度的研究。于利民等(1996)對中國大陸有史以來的強地震發生時間間隔進行分析,得出強震時間間隔的經驗概率和模擬函數。范琦(2001)以甘肅天祝地區為例,利用最大熵原理進行分析檢驗,結果認為地震震級和地震間隔時間的變化服從威布爾分布。李強和徐桂明(2002)利用威布爾分布建立了江蘇及鄰區中強地震的時間間隔分布模型,并對該地區未來幾年發生中強地震的概率進行預測。王煒和劉震華(1987)探討了華北地區小震時間間隔的統計分布,經過統計檢驗認為它服從威布爾分布。張亞奎(1996)對吉林省伊通—舒蘭斷裂帶以及我國南北地震帶南段安寧河—滇東地震帶發生的地震進行了統計分析,發現在大的地震帶上發震時間間隔服從正態分布。Utsu(1984)使用威布爾分布、伽馬分布、對數正態分布和指數分布估計日本地區的地震發生時間間隔。Sadeghian(2018)研究適用于伊朗地區的地震時間間隔數據的概率密度函數,結果為伽馬分布、貝塔分布、指數分布、對數正態分布、正態分布和威布爾分布。Khan等(2021)根據擬合優度檢驗,得出威布爾分布和對數邏輯分布是興都庫什地區最佳擬合模型,對數正態分布和伽馬分布分別為該地區第二和第三合適的概率模型。Somette和Knopoff(1997)基于威布爾分布、對數正態分布、指數分布、伽馬分布和貝塔分布等5種概率分布研究了地震時間間隔并預測地震復發時間。
2008年5月12日,四川汶川發生8.0級大地震,破壞性極強、波及范圍廣,造成了非常嚴重的人員傷亡和經濟損失。在四川發生過的中小地震更是數不勝數,因此本文選取中國四川地區發生的中強地震作為研究對象,從指數分布、伽馬分布、對數正態分布和威布爾分布4種概率分布出發,分別研究1970年至2007年和1970年至今發生的M5.0及以上地震,通過擬合優度檢驗找出地震發生時間間隔的最佳擬合分布,并評估未來發生M5.0及以上地震的概率。
根據國家地震科學數據中心記載(1970年1月1日至2021年7月31日的地震數據),四川自1970年至2021年7月曾發生過164次M5.0及以上的地震。由于2008年汶川地震余震過于頻繁,若使用2008年全部地震數據預測未來發生強震則會導致概率值偏高,因此剔除2008年發生的余震數據。本文使用的數據見表1,共121條數據。

表1 地震數據
本文主要探討兩部分內容:(1)使用1970年至2007年的地震數據,計算兩次地震時間間隔(以天為單位),根據計算結果做模型擬合和預測,通過預測2008年發生M5.0及以上地震的可能性從而驗證模型選擇方法的可靠性;(2)使用表1全部數據用同樣的方法預測未來發生M5.0及以上地震的可能性。
本文分析4種概率分布,首先介紹這4種概率分布的概率密度函數、模型參數(見表2)。

表2 4種概率分布的概率密度函數
(1)指數分布:指數分布是伽馬分布和威布爾分布的特殊形式,具有無記憶性,即對隨機變量T,有P(T>s+t|T>t)=P(T>s),且(t,s>0),常用來表示獨立隨機事件發生的時間間隔,其參數λ表示單位時間內事件發生的平均次數。其在可靠性研究具有廣泛應用。
(2)伽馬分布:伽馬分布的參數α為形狀參數,β為尺度參數,當α=1時,伽馬分布為指數分布。龔平等(2001)指出當地震的發生為相互獨立事件且僅依賴于時間間隔,與時間起點無關時,可得出地震發生時間服從伽馬分布。
(3)對數正態分布:對數正態分布的參數μ為位置參數,σ為形狀參數,與正態分布相似,在可靠性研究中,數據若不符合正態分布,則常取其對數使之符合正態分布,因此稱作對數正態分布。
(4)威布爾分布:威布爾分布的參數η為比例參數,m為形狀參數,當m=1時,威布爾分布為指數分布。該分布在地震研究中也具有廣泛的應用。
本文使用常用的模型選擇方法——赤池信息準則(AIC準則)、貝葉斯信息準則(BIC準則)和K-S檢驗。
AIC可以衡量統計模型擬合的優良性,它建立在熵的概念上,提供了權衡估計模型復雜度和擬合數據優良性的標準。AIC定義為式(1):

式中k為模型參數個數,L為似然函數。從一組可供選擇的模型中選擇最佳模型時,通常選擇AIC最小的模型。AIC為模型選擇提供了有效的規則,但也有不足之處。當樣本容量很大時,在AIC準則中擬合誤差提供的信息就要受到樣本容量的放大,而參數個數的懲罰因子卻和樣本容量沒關系,因此當樣本容量很大時,使用AIC準則選擇的模型不收斂于真實模型,它通常比真實模型所含的未知參數個數要多。
BIC貝葉斯信息準則與AIC相似,可用于模型選擇,并且改進了AIC的不足之處。BIC定義為式(2):

式中k為模型參數個數,n為樣本數量,L為似然函數。BIC值越小,說明模型擬合程度越高。
單樣本K-S檢驗是一種擬合優度的非參數檢驗方法,利用樣本數據推斷總體是否服從某一理論分布。K-S檢驗一般返回兩個值:D和P值。其中D值表示兩條累計分布曲線之間的最大垂直距離,所以D值越小,這兩個分布的差距越小,分布越一致。P值是假設檢驗里面的P值,如果P值大于0.05,那么就不能拒絕原假設,所以P值越大,分布越一致。
本文首先擬合四川省2008年以前發生M5.0及以上地震時間間隔的概率分布模型,計算參數值,并根據AIC、BIC、D值和P值來檢驗指數分布、伽馬分布、對數正態分布和威布爾分布對地震發生時間間隔的擬合程度,運算結果如表3所示。根據數據的累積概率對應4種分布累積概率繪制p-p plot(見圖1),可以直觀的看出樣本數據是否服從某一分布。由于指數分布的AIC和BIC值相較于其他3種分布偏大很多,說明擬合效果最差。比較伽馬分布、對數正態分布和威布爾分布的K-S距離,可以看到在K-S檢驗中伽馬分布返回的D值最小,P值最大。由圖1可見,樣本數據與伽馬分布擬合的較好。因此,選用伽馬分布作為四川省2008年以前發生M5.0以上地震時間間隔的最佳擬合分布。

圖1 1970年至2007年地震數據4種概率分布P-P plot圖

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

表4 1970年至今地震數據(剔除2008年余震數據)建模的擬合檢驗結果

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

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

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