徐偉進 李雪婧 謝卓娟 呂悅軍 高戰(zhàn)武

1）中國地震局地球物理研究所，強地震學研究室，北京 100081

2）國家自然災害防治研究院，地震災害研究室，北京 100085

3）中國地震災害防御中心，北京 100029

引言

地震時間分布模型是進行地震預測和地震危險性分析的重要理論基礎，分析海域地震時間分布特征是進行海域地震危險性分析的重要研究內容。地震的發(fā)生在時間上符合泊松模型是當前概率地震危險性分析和地震預測中廣泛接受的基本假設。泊松模型表示地震發(fā)生是相互獨立的，即本次地震發(fā)生在時間上不受上次地震的影響，也不會對下次地震造成影響。地震時間獨立分布模型包含泊松模型、雙泊松模型、分段泊松模型、更新模型以及特征地震模型等（Cornell，1968；Gardner 等，1974；Schwartz 等，1984；胡聿賢，1990；Console 等，2003），其中泊松模型在地震學和工程地震學中應用廣泛。地震學家們針對地震時間分布是否符合泊松模型做了大量統(tǒng)計檢驗工作（Gardner 等，1974；Bufe 等，2005；Michael，2011；Shearer 等，2012；Parsons 等，2012），目前中國、美國等國家在地震區(qū)劃和地震危險性分析中均采用泊松模型（高孟潭，1996；潘華等，2013；Petersen 等，2014）。

相對于地震的時間獨立模型，還有用于描述地震時間分布的時間相依模型，時間相依模型依據Reid（1910）提出的彈性回跳理論。彈性回跳理論認為地震發(fā)生是斷層應變能漸進且連續(xù)的積累過程，當該斷層出現突然的應變釋放時就會發(fā)生構造地震，地震將大大減輕地殼的應變，在穩(wěn)定的構造力作用下，壓力慢慢地重新累積變大，最終發(fā)生下次地震。按照彈性回跳理論，斷層上本次地震的發(fā)生受上次地震的影響，同樣對下次地震產生影響，地震是時間相依的。近幾十年來，科學家們使用不同地區(qū)的地震目錄對地震時間相依性進行了實證研究，發(fā)現在許多情況下地震時間分布并不符合泊松模型，而時間相依模型能夠更好地描述地震的時間分布特征（Utsu，1984；Nishenko 等，1987；Ogata，1991；Tripathi，2006；Sharma 等，2010；Ellsworth 等，2015）?？茖W家們提出了采用伽馬模型（Gamma）、對數正態(tài)模型（Lognormal）、威布爾模型（Weibul）以及布朗過程時間模型（Brownian Passage Time，BPT）描述地震時間分布特征（Utsu，1984；Matthews 等，2002；Tripathi，2006；Pasari 等，2015，2018；Bajaj 等，2019），這些模型為地震預測和地震危險性分析提供了重要的理論支撐。

本研究中，以中國海域統(tǒng)一地震目錄為基礎資料，以泊松模型（指數分布）、伽馬模型、對數正態(tài)模型、威布爾模型以及布朗過程時間模型為目標模型，回歸了各模型參數，并通過AIC、BIC 判定準則和K-S 檢驗選擇能夠描述海域地震活動時間分布特征的最優(yōu)模型。采用擴散熵分析（Diffusion Entropy Analysis，DEA）和標準差分析（Standard Deviation Analysis，SDA）方法對中國海域地震時間叢集特征和時間相關性進行研究。

1 數據、模型與方法

1.1 數據

海域地震事件記錄是進行海域地震活動特征分析、海域地震危險性分析以及海域地震災害預測等的重要基礎數據，近幾十年來，隨著地震觀測技術的提高和更多地震臺網的布設，我國及周邊國家、地區(qū)均記錄了豐富的海域地震事件。地震學家收集了我國海域及周邊國家、地區(qū)的地震臺網記錄，經過分析和整理，編制了我國海域統(tǒng)一地震目錄，共計61 285 條地震事件（Xie 等，2021），該地震目錄為進行我國海域地震危險性分析提供了重要的基礎資料。本研究使用的海域地震目錄是Xie 等（2021）分析整理的中國海域統(tǒng)一地震目錄（圖1（a））。Xie 等（2021）還對海域地震目錄做了震級轉換和完整性分析等工作，本研究中直接采用其海域地震時間完整性分析結果。使用Gardner 等（1974）得出的時空窗法刪除了地震目錄中的余震事件，該方法被廣泛應用于地震區(qū)劃圖編制和地震活動性分析中的余震刪除（Shearer 等，2012；Daub 等，2012；Petersen 等，2014）。刪除余震后，還剩29 450 條主震事件，使用刪除余震后的主震目錄研究地震時間分布特征。

根據新劃分的中國海域及領區(qū)地震帶（圖1（b）），以地震帶為單元研究了海域地震時間分布特征。海域及領區(qū)地震帶是進行海域地震活動性參數計算的基本單元，研究海域地震帶地震時間分布特征對海域地震活動性參數的確定具有重要科學意義。

圖1 中國海域及鄰區(qū)地震分布Fig. 1 Earthquake distribution in China sea areas and adjacent areas

1.2 模型

本研究中選擇泊松模型（指數分布）、伽馬模型、對數正態(tài)模型、威布爾模型以及布朗過程時間模型作為參考分布模型，通過統(tǒng)計檢驗選擇能夠描述海域地震事件分布特征的最優(yōu)模型。

（1）指數模型（指數分布）

指數分布是無記憶分布，可由指數分布在極端情況下推導得出，地震學家們通過檢驗地震時間間隔是否符合指數分布判斷地震是否符合指數分布模型。指數分布概率密度函數為：

（4）威布爾（Weibull）模型由Weibull 等（1951）提出，Kagan（1997）認為該模型可以用于評估地震的復發(fā)概率。Weibull 概率密度函數f(x)表達式如下：

1.3 方法

采用擴散熵分析和標準差分析方法研究中國海域地震的時間叢集特征和時間相關性。DEA 法是Scafetta等（2004）在計算青少年生育現象時首先采用的，該方法在是否為高斯過程的情況下還是在變異數為無窮大的情況下，都能得到正確的標度參數。此后，該方法被廣泛應用于分析不同地區(qū)地震的時間分布特征（Jiménez 等，2006；Zhou 等，2016）。Mega 等（2003）將DEA 應用于解釋1976-2002 年加州地區(qū)不同地震叢集在時間上的分布關系，得出加州地區(qū)大尺度時間上的地震叢集和下個地震叢集具有相關性。Tsai 等（2008）用該方法分析了臺灣地區(qū)地震與地震之間的相關性，發(fā)現臺灣地區(qū)地震時間概率密度函數的標度值為0.83。

DEA 方法簡述如下：設定單位時間Δt，將地震序列分成m個單位時間序列ξ i，設定門檻值k，當第i個單位時間內地震個數大于k時，ξ i=1，否則ξi=0。以t為窗口寬度，移動窗口，每次移動1 個單位時間，由n=0 到n=m-t開始移動，共生成m-t+1 道軌跡，在第n道軌跡，窗口內所有ξi總和為：

2 時間分布模型結果分析

以中國海域地震目錄為基本輸入，采用極大似然估計方法，回歸了上文介紹的各模型參數。為選出最優(yōu)模型，采用赤池信息準則（Akaike Information Criterion，AIC）、貝葉斯信息準則（Bayesian Information Criterion，BIC）以及K-S 檢驗判斷模型擬合優(yōu)度。AIC 和BIC 的定義分別為：

式中，k為模型參數個數；L為模型似然函數；n為樣本量。訓練模型時，增加參數數量相當于增加模型復雜度，會增大似然函數，也會導致過擬合現象，針對該問題，AIC 和BIC 均引入與模型參數個數相關的懲罰項，BIC 的懲罰項較AIC 的懲罰項大，考慮了樣本數量，樣本數量過多時，可有效防止模型精度過高造成的模型復雜度過高。從一組可供選擇的模型中選擇最佳模型時，通常選擇AIC 或BIC 最小的模型。

表1 為根據整個中國海域及鄰區(qū)地震目錄回歸得到的各模型參數值和模型擬合優(yōu)度判定準則參數值。對于5.0 級以上地震，從表1 中可以看出，指數分布、威布爾分布以及伽馬分布的AIC、BIC 以及K-S 檢驗值相近，其中威布爾分布具有最小的AIC 和BIC，指數分布具有最小的K-S 檢驗值（K-S 檢驗值表示地震經驗分布與理論模型的差異，該值越小表示越接近），這意味M≥5 級地震不拒絕時間獨立性的假設，但也具有一定的時間相關性。對于6.0 和7.0 級以上的地震，根據BIC 和K-S 檢驗值，可判定指數分布能夠最優(yōu)地表達海域地震的時間分布特征。由圖2 的累積分布函數曲線可以看出，指數分布、威布爾分布以及伽馬分布與實際數據吻合較好，對數正態(tài)分布和布朗過程時間分布與實際數據吻合較差。由上述分析可知，對于整個中國海域及鄰區(qū)，指數分布能夠最佳描述地震的時間分布特征，威布爾分布和伽馬分布與實際觀測數據吻合良好，但具有2 個參數，模型相對復雜。

圖2 中國海域及鄰區(qū)不同震級地震的時間間隔累積經驗分布函數及其對應的5 個模型累積分布函數Fig. 2 The cumulative empirical distribution function of the time interval of M≥5, 6, 7 earthquakes in the Sea area of China and adjacent areas and the cumulative distribution function of the corresponding five models

表1 中國海域及鄰區(qū)模型參數及AIC、BIC 和K-S 檢驗值Table 1 Model parameters and AIC, BIC and K-S test results for earthquakes in China Sea and adjacent areas

表2 為根據華南沿海地震帶地震目錄回歸得到的各模型參數值和模型擬合優(yōu)度判定準則參數值。從AIC、BIC 和K-S 檢驗值上看，對于M≥4 級的地震，指數分布、威布爾分布和伽馬分布均表現優(yōu)良，模型間差異較小，K-S 檢驗接受地震時間間隔分布符合這3 種分布的假設；對于M≥5 級的地震，指數分布被判定為最優(yōu)模型，威布爾分布和伽馬分布也有良好表現；對于M≥6 的地震，K-S 檢驗顯示所有模型均不拒絕時間獨立性的假設，所有模型的AIC 和BIC 判定準則值均非常接近，特別是對于M≥4 和M≥5 級地震表現較差的對數正態(tài)分布和布朗過程時間分布，對于M≥6 級的地震表現優(yōu)良，擬合優(yōu)度好于威布爾分布和伽馬分布。這意味著在發(fā)震構造特征相似的地震帶內，大地震時間分布模型可用對數正態(tài)分布和布朗過程時間分布描述，這與國際上的研究結果相似。

表2 華南沿海地震帶模型參數及AIC、BIC 和K-S 檢驗值Table 2 Model parameters and AIC, BIC and K-S test results for earthquakes in South China Coastal seismic zone

表3 為根據臺灣西部地震帶的地震目錄回歸得到的各模型參數值和模型擬合優(yōu)度判定準則結果。從AIC、BIC 和K-S 檢驗值上看，對于M≥4 和M≥5 級的地震，指數分布為最優(yōu)分布，威布爾分布和伽馬分布的擬合優(yōu)度較好；對于M≥6.0 級的地震，正態(tài)對數分布和布朗過程時間分布擬合優(yōu)度最好。這再次說明對于大地震，采用對數正態(tài)分布、布朗過程時間分布等模型來描述其時間分布特征更合適。

表3 臺灣西部地震帶模型參數及AIC、BIC 和K-S 檢驗值Table 3 Model parameters and AIC, BIC and K-S test results for earthquakes in the western Taiwan earthquake zone

表4 為根據臺灣-馬尼拉海溝地震帶地震目錄回歸得到的各模型參數值和模型擬合優(yōu)度判定準則參數值。從AIC、BIC 和K-S 檢驗值上看，對于M≥5 級的地震，威布爾分布擬合優(yōu)度最好，伽馬分布與威布爾分布較接近；對于M≥6 級的地震，指數分布為最優(yōu)分布，威布爾分布和伽馬分布表現較好；對于M≥7 級的地震，綜合判定指數分布為最優(yōu)模型，威布爾分布和伽馬分布具有較好的擬合優(yōu)度。從地震事件分布模型上可以看出，在該地震帶內5 級以上地震具有較強的時間相依性，6 級以上地震可接受指數分布假設，但也具有一定的時間相關性。這是由于臺灣-馬尼拉海溝地震帶地震活動非常強烈，經常發(fā)生8 級左右的地震，大震的發(fā)生對稍小震級地震具有一定影響，且震級越小所受影響越大。

表4 臺灣-馬尼拉海溝地震帶模型參數及AIC、BIC 和K-S 檢驗值Table 4 Model parameters and AIC, BIC and K-S test results for earthquakes in Taiwan - Manila trench seismic zone

本文計算了琉球海溝地震帶各模型參數值和模型擬合優(yōu)度判定準則參數值，得到了與上述地震帶相似的結論。特別地，在琉球海溝地震帶，對于M≥7 級地震，正態(tài)對數分布和布朗過程時間分布擬合優(yōu)度最好。這再次說明了對數正態(tài)分布和布朗過程時間分布在描述大地震時間分布特征的方面具有優(yōu)越性，同時意味著琉球海溝地震帶內7 級以上地震具有一定的時間相關性。

根據上述分析可知，對于震級相對較小的地震（M＜6），指數分布、伽馬分布以及威布爾分布均能較好地描述地震時間分布特征。在大的區(qū)域范圍內（如整個海域），震級相對較大的地震（M＞6）可采用指數分布描述其時間分布特征。在較小的區(qū)域范圍內（如地震帶），大地震時間間隔可能更符合對數正態(tài)分布和布朗過程時間分布。

3 地震時間叢集性和相關性分析

地震在時間上呈現叢集現象是相對于地震在時間上是完全隨機分布（指數分布）而言的，是指地震在一段時間內集中成組發(fā)生，而后在一段時間內表現出相對頻度降低的特征，然后又集中成組發(fā)生的現象。地震時間叢集一直是地震學家關心的重要問題，其對地震預測和地震危險性分析具有重要影響。圍繞地震在時間上是否叢集，地震學家們做了許多探索和研究，有的研究表明全球地震在時間上遵從指數分布（Michael，2011；Shearer 等，2012；Daub 等，2012；Parsons 等，2012），有的研究認為地震呈現叢集特征，但每個叢集之間相互獨立（Mega 等，2003），還有的研究表明地震之間存在相互影響且呈現叢集特征（Scafetta 等，2004；Bufe，2005；Kulkarni 等，2013；Salditch，2020）。

本研究從地震時間間隔的變異系數(Coefficient of Variation)入手，采用擴散熵分析和標準差分析方法對中國海域地震時間叢集特征進行研究。

一套地震時間序列中，若存在叢集性，則在統(tǒng)計量上表現為地震之間的時間間隔離散性變大（相對于泊松模型），變異系數是比較2 組數據（本研究為實際觀測數據和泊松模型數據）離散程度大小的理想參數，其定義為標準差與均值的比值。對于實際地震記錄，連續(xù)地震之間的時間間隔是確定的，因此可以計算出1個變異系數。對于指數分布，其變異系數的理論值為1，但由于受樣本量的影響，符合泊松模型的隨機變量的變異系數不可能完全為1，而是呈現以1 為均值的分布。我們采用蒙特卡洛方法模擬1 000 套和實際觀測數據樣本相同的符合指數分布的隨變量，計算每1 套模擬數據的變異系數，將其分布與實際觀測數據的變異系數進行比較。圖3 為整個海域M≥5、M≥6、M≥7 級地震變異系數與泊松模型變異系數的對比，可以看出在1%的置信度下，M≥5 地震被拒絕接受為指數分布的假設，其變異系數即使在1%置信度的極端情況下也大于指數分布的變異系數，這意味著M≥5 地震在時間上具有叢集性。M≥6、M≥7 級地震的變異系數在5%置信度下接受為指數分布的假設，這也意味著海域地震中主要是M5～M6 級地震體現出時間叢集性。

圖3 海域地震時間間隔變異系數與泊松模型變異系數對比Fig. 3 Comparison between the variation coefficient of sea area earthquake time interval and the variation coefficient of Poisson model

采用擴散熵分析和標準差分析方法分析地震的時間相關性。擴散熵分析法的標度參數為δ，標準差分析法的標度參數為H（Hurst exponent）。當隨機變量的δ=H=0.5 時，表示隨機變量符合指數分布，變量之間無相關性，實際情況下若計算的標度值接近0.5±0.05，則認為變量無相關性；當H＞δ＞0.5 時，則認為變量之間具有相關性且具有長期記憶性，即在某一叢集期內，一次地震的發(fā)生促使下一次地震的發(fā)生；當H＜0.5 時，則認為地震之間具有負相互作用，即一次地震的發(fā)生使下次地震延遲了，在實際情況下，由于地震樣本量較少，即使符合指數分布，計算的標度值偏小，這已被蒙特卡洛模擬結果證明，因此若計算的標度值小于0.5，則認為變量符合指數分布。

圖4 為整個海域地震的擴散熵分析法和標準差分析法分析結果，可以看出對于M≥5 級以上地震，標度參數均大于0.5，這意味著海域5 級以上地震之間具有相關性；對于M≥6 和M≥7 地震，標度參數則非常接近0.5，表現出完全隨機性，說明地震之間無相關性。這意味海域地區(qū)M5～M6 地震受更大地震（M≥6、M≥7）活動的影響，大地震發(fā)生后促使了中小地震的發(fā)生。

圖4 整個海域地震的擴散熵和標準差隨時間變化及其標度值Fig. 4 Diffusion entropy and standard deviation of earthquakes as the function of time in the whole sea area and the scale values

針對不同地震帶計算不同震級地震目錄的擴散熵分析法和標準差分析法標度參數（表5），對于華南沿海地震帶，M≥4 級以上地震表現出時間相關性，M≥5、M≥6 級以上地震可認為符合指數分布，這說明該地震帶M4～M5 級地震受M≥5、M≥6 級以上地震的影響。在長江下游-南黃海地震帶，5 級以上地震表現為時間叢集，且受更大地震（M≥6 級以上）的控制。對于臺灣西部地震帶，M≥5、M≥6 級地震均具有時間相關性，這說明該地震帶內M=5、M=6 級地震相互促進或受更大震級地震的控制。在臺灣-馬尼拉海溝地震帶，M≥5 地震呈現出時間叢集性，M≥6、M≥7 級地震則符合指數分布，這說明了M5～M6 級地震受更大地震的影響。在琉球海溝地震帶內，M≥7 級地震表現出相互促進作用，在地震預測和地震危險性分析中需要重點考慮。

表5 各地震帶不同起始震級地震目錄的擴散熵分析法和標準差分析法標度值Table 5 DEA and SDA scale values of earthquake catalogs with different initial magnitudes in each seismic zone

以往的研究中，科學家們關注大地震之間的時間相關性（Bufe，2005；Michael，2011；Shearer等，2012；Daub 等，2012；Kulkarni 等，2013；Salditch，2020），本研究表明，即使在刪除余震和前震后，地震目錄中的大地震與中小地震之間也具有時間相關性，這意味著大地震的發(fā)生促進了中小地震的發(fā)生，使地震活動在統(tǒng)計特征上表現出長期記憶性，中小地震在時間軸上呈現叢集特征。

4 結論與討論

由上述分析可知，在大的區(qū)域范圍內（如整個海域區(qū)域），起始震級稍小的地震（M≤5、M≤6）時間分布多與指數分布吻合更好，威布爾分布和伽馬分布的擬合優(yōu)度也較好；起始震級大的地震（M≥6、M≥7）則完全符合指數分布。在一個較小區(qū)域范圍內（如地震帶），指數分布、威布爾分布和伽馬分布均能較好描述震級較小地震的時間分布特征，對于起始震級大（M≥6、M≥7）的地震，對數正態(tài)分布和布朗過程時間分布有時具有較好的擬合優(yōu)度。這與其他學者的研究結果相似（Matthews 等，2002；Bajaj 等，2019）。由于威布爾分布和伽馬分布在特殊情況下可退化為指數分布，因此，在較大區(qū)域內，在遵從科學性的前提下，本著便于使用的原則，泊松模型可選擇作為進行地震危險性分析的基本模型。

由于大地震發(fā)生率對于地震危險性具有重要的意義，在較小區(qū)域內（如某潛在震源區(qū)或某條斷層），若歷史地震目錄較為豐富，可嘗試采用對數正態(tài)模型或布朗過程時間模型計算高震級地震的發(fā)生概率和發(fā)生頻度，這為地震危險性分析中地震發(fā)生率的計算提供了新途徑，從而使地震危險性計算結果更加科學合理。目前，也有一些學者和機構采用對數正態(tài)分布模型或布朗過程時間模型來計算地震危險性（Hebden 等，2009；Working Group on California Earthquake Probabilities，2013；Petersen 等，2014），并取得了新的認識，為地震災害的預測預防提供了重要支撐。

對海域地震時間叢集性和相關性進行分析，結果表明，在某區(qū)域內稍小震級的地震易受更大地震的影響，在時間上呈現叢集性，而震級大的地震則更多地表現出完全隨機性（符合指數分布）。一方面這是由于目前沒有有效的方法區(qū)分主震和余震，即使刪除余震，地震目錄中仍然包含數量相當的余震（中小震）。另一方面，在大地震發(fā)生后，應力擴散和調整需要較長時間，表現為中小震級地震的發(fā)生（Shearer 等，2012），在統(tǒng)計特征上表現為叢集性和長期記憶性（即H＞0.5），因此，在對某地區(qū)進行地震危險性分析時應充分考慮該地區(qū)最近一次發(fā)生的強震對未來中小地震的影響。本研究還發(fā)現在地震活動非常強烈的琉球海溝地震帶，M≥7 級地震表現出正向的相關性，即某次7 級以上地震的發(fā)生會促使下次7 級以上地震的發(fā)生，需要引起特別注意，這一發(fā)現對于理解地震孕育發(fā)生機理具有一定科學意義。