徐偉進(jìn) 高孟潭 李宗超
(中國地震局地球物理研究所,北京 100081)
地震目錄是指按照時間順序,對地震的時、空、強(qiáng)等參數(shù)進(jìn)行收錄,編輯成冊的目錄資料。地震目錄是地震學(xué)研究的基礎(chǔ)資料,對地震目錄的分析研究是地震學(xué)研究的基本內(nèi)容之一。一般來說,地震目錄的分析處理主要包括:余震刪除、完整性分析、地震的時空特征分析等。
本文通過中蒙兩國地震工作者的共同努力,搜集整理得到了蒙古國及其周邊地區(qū)的地震目錄。目前關(guān)于對蒙古國地震目錄分析處理的工作還比較少見。為了使廣大地震科研工作者能夠進(jìn)一步認(rèn)識蒙古國的地震活動特征,以及使該目錄能夠方便地應(yīng)用于蒙古國的地震危險性分析、地震預(yù)測等地震學(xué)研究中,作者對該目錄做了余震刪除、時間完整性分析以及時間統(tǒng)計(jì)分布特征分析等工作。
當(dāng)前的地震危險性分析及地震預(yù)測模型中,大都將地震時間分布看成是泊松分布(Cornell,1968;高孟潭,1996;胡聿賢等,2001;Kagan等,1994;Frankel,1995)。然而,小區(qū)域、短時間內(nèi)的余震叢集是不符合泊松分布的,這就需要把余震從地震序列中刪除。目前的余震刪除方法主要有:時空窗法(Utsu,1969;Gardner等,1974)、基于地震事件之間時空距離的鏈接法(Resenberg,1985)、隨機(jī)除叢法(Zhuang等,2002)。Wang等(2010)系統(tǒng)比較了幾種余震刪除方法的差異,結(jié)果顯示,沒有明顯的證據(jù)表明哪一種方法更有優(yōu)勢。考慮到計(jì)算的簡便和使用的廣泛性,文中作者使用Gardner等(1974)的方法進(jìn)行余震刪除。
地震目錄時間完整性分析是地震目錄分析的重要工作。在地震預(yù)測、地震危險性分析及其它地震活動性研究中,不完整的地震目錄會給地震活動性參數(shù)的估計(jì)帶來偏差,因此,在研究一個地區(qū)的地震活動性時,分析這一地區(qū)地震目錄的可靠性是非常有必要的。本文中作者將采用Albarello等(2001)提出的新的統(tǒng)計(jì)方法對蒙古國的地震目錄進(jìn)行完整性分析。
地震的時間獨(dú)立性,即地震發(fā)生在時間上是一個泊松過程,是被地震學(xué)界廣泛接受的假說之一,這一假說是地震危險性分析及其它地震學(xué)研究中的最基本理論之一。文中作者采用χ2檢驗(yàn)、K-S檢驗(yàn)、游程檢驗(yàn)(Run-test)對蒙古國地震目錄的時間統(tǒng)計(jì)分布特征做了初步分析。關(guān)于檢驗(yàn)方法的詳細(xì)介紹可參考Gibbons等(2003)的文章。
文中采用的地震目錄由兩套目錄合并而成。一套目錄是中蒙邊境中國一側(cè)的區(qū)域臺網(wǎng)記錄的目錄,另一套為蒙古國天文與地球物理研究所記錄的目錄。由于兩套目錄中有相同的地震記錄,因此在合并前進(jìn)行了甄選和剔除,并對部分地震進(jìn)行了重新定位。此外目錄中還包含蒙古國根據(jù)歷史文獻(xiàn)記載整理的歷史地震目錄。目錄的起始震級為MS3.5級,記錄周期為公元1119—2010年,共計(jì)7500條。
圖1為蒙古國及周邊地區(qū)地震目錄空間分布圖。從圖中可以看出,蒙古國地震活動主要集中在中部和西部地區(qū),東部地區(qū)地震活動較弱且較為分散。圖中灰色粗邊框以內(nèi)的地區(qū)為文中的研究區(qū)域,作者將選取這一區(qū)域內(nèi)的地震目錄進(jìn)行分析。
由于地震臺網(wǎng)監(jiān)測能力的差異,地震目錄的完整性在空間上往往是有差異的,因此在進(jìn)行地震目錄的完整性分析時,可采用空間網(wǎng)格劃分的方法進(jìn)行計(jì)算。然而蒙古國地震目錄樣本較少,不適合采用空間網(wǎng)格劃分的方法進(jìn)行計(jì)算。因此,文中為了表示地震完整性的空間差異,根據(jù)蒙古國的地質(zhì)構(gòu)造及地震活動性分布劃分了10個研究區(qū)域進(jìn)行地震目錄的完整性計(jì)算(見圖1灰色區(qū)域)。圖中用編號表示研究區(qū)域的名稱,如An表示第n個研究區(qū)域。其中A0表示整個研究區(qū)域。
圖1 蒙古國地震目錄空間分布及研究區(qū)示意圖Fig. 1 Epicentral distribution of earthquakes(M≥3.5)reported in the Mongolian seismic catalog and study areas(gray areas)
在前言中已經(jīng)簡單介紹了地震目錄的余震刪除方法。在目前的地震活動性研究中,特別是概率地震危險性分析及概率地震預(yù)測中,人們普遍采用Gardner等(1974)提出的時空窗余震刪除方法,其主要原因是該方法簡單易操作,并且具有扎實(shí)的數(shù)學(xué)、物理基礎(chǔ)。因此作者也采用這一方法進(jìn)行余震刪除。在該方法中,空間窗口是根據(jù)下式確定的(Knopoff等,1972):
式中,a、b為固定參數(shù)值,a、b的取值分別為0.5和-1.78(劉杰等,1996)。時間窗口根據(jù)表1給出,震級不在表中的根據(jù)線性插值得出。
表1 不同震級地震的余震時間窗(Gardner等,1974)Table 1 Time window algorithm for aftershocks with different magnitudes(Gardner et al.,1974)
圖2為刪除余震前后地震活動率的對比曲線,可以看出在未刪除余震前地震活動率出現(xiàn)幾個峰值,這是由于在這一時間有震級較大的地震發(fā)生,從而在短時間內(nèi)產(chǎn)生大量余震事件。而在刪除余震后,地震活動率不再有明顯的峰值,背景地震活動隨時間變化基本平穩(wěn)。
圖2 余震刪除前后地震活動率比較(M≥3.5)Fig. 2 Seismic rate as a function of time for cluster(red line)and decluster(blue line)catalogs(M≥3.5)
首先采用Albarello等(2001)提出的新的統(tǒng)計(jì)方法,分析蒙古國及其周邊地區(qū)地震目錄的完整起始時間。其方法為:
假設(shè)一套地震目錄是完整的(假設(shè)C)、具有代表性的(假設(shè)R),則可將計(jì)算每年的平均年發(fā)生頻次作為地震活動率,將其稱為樣本估計(jì)。一般來說,這里計(jì)算的地震活動率是符合泊松分布的,但是根據(jù)中心極限定理,可以近似將其看作是符合高斯分布的。由于高斯分布的對稱性,在兩個不同時間段itΔ和jtΔ內(nèi)觀測到的地震活動率in和jn應(yīng)滿足下式:
式中,p是相關(guān)概率。
如果接受假設(shè)R,即地震的發(fā)生是一個獨(dú)立的過程,那么就認(rèn)為在ΔT時間段內(nèi)地震目錄是完整的。
為了實(shí)現(xiàn)上述假設(shè),首先將地震目錄的時間段TΔ以tδ為間隔分成2N個不重合的子段,對于每一子段itδ計(jì)算該段內(nèi)的地震活動率in,基于這些地震活動率定義一個二進(jìn)制的隨機(jī)變量D:
實(shí)際上,變量D是地震目錄前一部分地震活動率和后一部分地震活動率的比較。若D=0的次數(shù)較多,就意味著時間較早的地震有缺失,地震目錄完整程度下降。
對于N′個可用的地震目錄子段,m為D=0的次數(shù)。如果目錄是完整的,那么實(shí)際觀測中計(jì)算得到的次數(shù)大于等于m的概率可用二項(xiàng)式計(jì)算:
根據(jù)式(1),p值取0.5。這一概率主要依賴于假設(shè)R的置信程度。概率值越小,意味著目錄完整程度越低。
上述過程是地震子目錄在整個目錄中具有代表性的條件下其完整性的概率。顯然希望得到的是某一具體實(shí)際段iTΔ內(nèi)地震完整的非條件概率:表示地震目錄在時間段iTΔ內(nèi)具有統(tǒng)計(jì)代表性,即假設(shè)R的置信程度。
為了估計(jì)假設(shè)R的概率,應(yīng)當(dāng)考慮地震孕育過程的長期統(tǒng)計(jì)特征,當(dāng)然這一概率并非是一個先驗(yàn)概率,可將假設(shè)R的概率簡單寫為:式中,
式中,ΔTmax是采用的整個地震目錄的記錄周期;k是一個大于1的定值。
因此,隨著ΔTi的增大,是線性遞增的。
聯(lián)合式(3)、(4)、(5),則:
然而一般來說這是很難得到的,因此需要將 Pi值重新標(biāo)度,用以得到新的完整性函數(shù),即:
這一估計(jì)的不確定性可以用四分位來計(jì)算:
其中
文中主要分析了圖1中10個研究區(qū)域和蒙古全境MS≥3.5級地震的完整起始時間。此外,還分析了蒙古全境M≥5級和M≥6級地震的完整起始時間。
表2為各研究區(qū)地震目錄的完整起始時間及其上下限。其中,cT為完整起始時間;uT為不確定性的上限值;lT為不確定性的下限值。從中可以看出,蒙古國MS≥3.5級地震的完整起始時間大都在1970年以后,誤差約為±10年左右。對于儀器記錄的地震而言,這一完整起始時間相對較晚且離散度較大,這可能與蒙古國臺網(wǎng)架設(shè)密度及經(jīng)濟(jì)、科學(xué)技術(shù)發(fā)展水平有關(guān)。蒙古國M≥5.0級地震的完整起始時間約為1945年前后,誤差約為±15年;M≥6.0級地震的完整起始時間約為1900年前后,誤差約為±30年。可以看出,隨著起始震級的增大,完整起始年的誤差也越大,這是由于隨著震級的增大,地震事件樣本長度急劇減小,從而導(dǎo)致了不確定性的增加。其次,對于大地震來說其復(fù)發(fā)周期很長,從而也導(dǎo)致了其完整起始年不確定較大。
從表2第3列各研究區(qū)的地震目錄完整起始記錄時間可以看出,蒙古國及其周邊地區(qū)地震目錄的完整起始時間具有一定的空間差異,這可能與每個地區(qū)地震儀器的布設(shè)密度有關(guān)。在地震儀器布設(shè)密度較大的地區(qū),地震監(jiān)測能力強(qiáng),地震的完整起始時間就稍早一些,反之則稍晚一些。此外,與儀器性能、人員素質(zhì)以及該地區(qū)地震地質(zhì)研究工作開展的早晚等也有一定關(guān)系。
表2 各研究區(qū)不同起始震級地震目錄完整起始時間 cT及不確定性 uT和 lTTable 2 Completeness timecTvalues and uncertaintyuTandlTvalues for seismic zones
地震的時間獨(dú)立性,即地震發(fā)生在時間上是一個泊松過程,是被地震學(xué)界廣泛接受的假說之一,這一假說是地震危險性分析及其它地震學(xué)研究中的最基本理論之一。本文以地震時間分布符合泊松分布為零假設(shè)H0,采用χ2檢驗(yàn)、K-S檢驗(yàn)、游程檢驗(yàn)(Run-test)對蒙古國地震目錄在時間上是否符合泊松分布做了檢驗(yàn)。并采用Ripley’s K-function(Ripley,1976;1977)對地震時間的叢集特征做了分析。
(1)2χ檢驗(yàn)
2χ檢驗(yàn)法是在總體X的分布未知時,根據(jù)來自總體的樣本,檢驗(yàn)關(guān)于總體分布假設(shè)的一種檢驗(yàn)方法。使用2χ檢驗(yàn)對總體分布進(jìn)行檢驗(yàn)時,首先提出的原假設(shè)為:
H0:總體X的分布函數(shù)為F0(x)
然后,根據(jù)樣本的經(jīng)驗(yàn)分布和所假設(shè)的理論分布之間的吻合程度,來決定是否接受原假設(shè)時間相關(guān)性檢驗(yàn)。2χ檢驗(yàn)的基本步驟為:
②將落入第k個小區(qū)間 Ak的樣本值的個數(shù)記作 nk,稱為實(shí)際頻數(shù)。
③根據(jù)所假設(shè)的理論分布,可以算出總體X值落入每個kA的概率kp,于是,knp就是落入kA的樣本值的理論頻數(shù):
如何表示實(shí)際頻數(shù)與理論頻數(shù)的差異,Pearson(1900)引入了如下統(tǒng)計(jì)量來表示經(jīng)驗(yàn)分布與理論分布的差異:
2χ檢驗(yàn)比較的是經(jīng)驗(yàn)分布與理論分布概率密度函數(shù)的差異,并且在實(shí)際應(yīng)用中第一步小區(qū)間間隔的選取是十分重要的。具體2χ檢驗(yàn)方法可參考《應(yīng)用統(tǒng)計(jì)學(xué)》(倪加勛等,1993)。
(2)Kolmogorov-Smirnov(K-S)檢驗(yàn)
Kolmogorov-Smirnov檢驗(yàn)是用來檢驗(yàn)一組具有未知分布XF的隨機(jī)數(shù)是否來自某一個特定連續(xù)的分布函數(shù)0F,其原假設(shè)為:
Glivenko-Cantelli theorem (Gibbons等,2003) 認(rèn)為,樣本的經(jīng)驗(yàn)分布函數(shù)(empirical distribution function) Sn(x)在所有的x值處為)提供了點(diǎn)估計(jì)。因此在樣本容量較大時,經(jīng)驗(yàn)分布函數(shù)值與真實(shí)分布函數(shù)值的差別是很小的。因此Kolmogorov-Smirnov檢驗(yàn)的統(tǒng)計(jì)量可寫為:
同樣在給定顯著性水平α后,可通過查表獲得統(tǒng)計(jì)量在該顯著性水平下統(tǒng)計(jì)量的臨界值。關(guān)于K-S檢驗(yàn)的方法可參考Gibbons等(2003)的文章。
(3)游程檢驗(yàn)(Runs-test)
游程檢驗(yàn)是一個以隨機(jī)過程為平穩(wěn)的且不具有自相關(guān)性為原假設(shè)的檢驗(yàn)方法。在任何一個由實(shí)數(shù)組成的時間序列里,一個游程是指一個連續(xù)不間斷的且具有相同符號的數(shù)值組成的子序列,并且時正、時負(fù)交替進(jìn)行。這里需要特別指出的是,可以考慮關(guān)于樣本均值的游程隨機(jī)性,對于大于均值的數(shù)據(jù)將其定為正號,小于均值的數(shù)據(jù)定為負(fù)號。假設(shè)有N個數(shù)據(jù),其中p個定為正號數(shù)據(jù),n個定為負(fù)號數(shù)據(jù),那么,得到的游程個數(shù)小于等于觀測值R的概率為:
若 n, p都大于10,那么R的分布可以近似看成符合以均值為,方差為的高斯分布,顯然,可以將該正態(tài)分布寫成統(tǒng)計(jì)量為
(4)K-function叢集檢驗(yàn)
一種廣泛應(yīng)用于估計(jì)單一點(diǎn)過程在時間和空間上是否叢集的方法為 Ripley’s K-function(Ripley,1976;1977)。其定義為:
式中,λ為單位時間或面積內(nèi)的事件數(shù)目,即事件密度(intensity)。對于一維的時間點(diǎn)過程K-function,可用下式進(jìn)行計(jì)算:
式中,L為總的時間段;n為總的觀測事件數(shù);t為時間;ijd為時間i和時間j之間的間隔;I等于1,若ijd<t,否則為0;ijw為邊界修正項(xiàng)。
為了更加直觀地表現(xiàn)出點(diǎn)過程在不同時間尺度上的分布特征,Hwang(2005)提出了關(guān)于時間點(diǎn)過程的L函數(shù):
表3為采用三種檢驗(yàn)方法進(jìn)行檢驗(yàn)得到的結(jié)果,可以看出采用游程檢驗(yàn)和卡方檢驗(yàn),在顯著性水平α=0.05和α=0.01的情況下,均接受地震目錄在時間上符合泊松分布這一假設(shè)。只是采用卡方檢驗(yàn)時p值小于0.1,接受H0的理由不是很強(qiáng)。采用K-S檢驗(yàn),在顯著性水平α=0.05時不接受H0,在α=0.01時接受H0,由此可以看出接受H0的理由是比較弱的。
上述檢驗(yàn)結(jié)果表明,蒙古國及其周邊地區(qū)地震目錄符合泊松分布的理由是不夠強(qiáng)的。這可能是由于背景地震活動的時間非平穩(wěn)性決定的。為此,作者采用時間K-function分析方法研究了蒙古國及其周邊地區(qū)地震目錄的時間叢集特征。為了更加直觀地表現(xiàn)出地震事件的時間叢集特征,往往把K-function換算成L-function。關(guān)于K-function分析方法的具體介紹可參考徐偉進(jìn)(2012)的博士論文。
圖3為蒙古國及其周邊地區(qū)地震目錄L-function隨時間的變化曲線。圖中黑線表示泊松分布L(t)-t的理論值。實(shí)測數(shù)值若在黑線之上表示事件具有時間叢集性;實(shí)測數(shù)值若在黑線之下則表示事件分布為規(guī)則分布,實(shí)際計(jì)算中這一現(xiàn)象往往是由邊界效應(yīng)造成的。圖中紅線為采用未刪除余震的地震目錄計(jì)算得到的L(t)-t值,可以看出未刪除余震地震目錄的L(t)-t值基本在黑線之上,說明地震目錄在時間上是叢集的。一般認(rèn)為,L(t)-t值的峰值所對應(yīng)的時間為叢集周期,因此,從紅線可以看出未刪除余震的目錄具有6年和15年左右的叢集周期。圖3中藍(lán)線為刪除余震后地震目錄的L(t)-t值,可以看出藍(lán)線在6年左右的時間點(diǎn)對應(yīng)一個峰值,而后L(t)-t逐漸減小,趨近于黑線。這意味著,即使在刪除余震后,地震目錄仍具有6年左右的叢集周期,說明蒙古國及其周邊地區(qū)地震活動具有時間非平穩(wěn)性。
表3 刪除余震后的地震目錄卡方、K-S及游程檢驗(yàn)結(jié)果Table 3 Result of χ2-test,K-S test,Run-test for decluster catalog
圖3 蒙古國地震目錄時間L-function變化曲線,紅色曲線為未刪除余震的地震目錄的L-function變化曲線,藍(lán)色曲線為刪除余震后的地震目錄的L-function變化曲線,黑線為泊松分布L-function理論值Fig. 3 L(t)-t value as a function of time for cluster (red curve) and decluster(blue curve) catalog. Black curve shows theoretical result of Poisson Model
通過上述對蒙古國及其周邊地區(qū)地震目錄的初步分析,作者認(rèn)為,蒙古國 M≥3.5級地震完整記錄起始時間為1977年;M≥5.0級地震完整記錄起始時間為1945年;M≥6.0級地震完整記錄起始時間為1900年。地震目錄完整記錄起始時間具有一定的空間差異,這可能與每個地區(qū)地震儀器的布設(shè)密度有關(guān),在地震儀器布設(shè)密度較大的地區(qū),地震監(jiān)測能力強(qiáng),地震的完整起始時間就稍早一些,反之,則稍晚一些。此外,與儀器性能、人員素質(zhì)以及該地區(qū)地震地質(zhì)研究工作開展的早晚等也有關(guān)系。
研究還表明,刪除余震后蒙古國地震活動隨時間變化基本平穩(wěn),時間統(tǒng)計(jì)分布特征符合泊松分布。但是在刪除余震后,蒙古國地震活動仍有6年左右的叢集周期,這說明蒙古國及其周邊地區(qū)地震活動具有時間非平穩(wěn)性。地震活動在時間上表現(xiàn)出的非平穩(wěn)性,可能是由于一下幾個方面的原因造成的:①在強(qiáng)震發(fā)生后,其形成的破裂帶會長期控制該區(qū)域的地震活動特征,從而地震活動呈長期叢集;②地震構(gòu)造的活動速率(如:滑移、蠕動等)是隨時間變化的(Lombardi等,2007);③區(qū)域或全球范圍內(nèi)長期應(yīng)力變化的相互影響(Ogata等,1991)。因此,在以后的地震危險性估計(jì)中為了更加科學(xué)合理地估計(jì)地震活動率,應(yīng)當(dāng)充分利用地球物理、大地測量、地球動力學(xué)等方面的資料對地震的活動趨勢進(jìn)行研究,以減小地震活動的非平穩(wěn)性對地震活動率估計(jì)值的影響。
本文的研究結(jié)果可以供其他學(xué)者進(jìn)行對比研究使用,得到的完整的地震目錄可以作為地震學(xué)研究的基礎(chǔ)資料。
高孟潭,1996. 基于泊松分布的地震烈度發(fā)生概率模型. 中國地震,12(2):195—201.
胡聿賢,高孟潭,2001.《中國地震動參數(shù)區(qū)劃圖(GB 18306-2001)》宣貫教材. 北京:中國標(biāo)準(zhǔn)出版社.
劉杰,陳棋福,陳颙,1996. 華北地區(qū)地震目錄完全性分析. 地震,16(1):59—67.
倪加勛,袁衛(wèi),易丹輝,蔡志洲,1993. 應(yīng)用統(tǒng)計(jì)學(xué). 北京:中國人民大學(xué)出版社,262.
徐偉進(jìn),2012. 地震危險性分析中地震時空統(tǒng)計(jì)分布模型研究 [博士學(xué)位論文]. 中國地震局地球物理研究所.
Albarello D.,Romano Camassi and Alessandro Rebez,2001. Detection of Space and Time Heterogeneity in the Completeness of a Seismic Catalog by a Statistical Approach:An Application to the Italian Area. Bull. Seism. Soc.Am.,91(6):1694—1703.
Cornell C.A.,1968. Engineeing seismic risk analysis. Bulletin of the Seismological Society of America,58:1583—1606.
Frankel A.,1995. Mapping Seismic Hazard in the Central and Eastern United States. Seismol. Res. Lett.,66:8—21.
Gardner J.K. and Knopoff L.,1974. Is the Sequence of Earthquakes in Southern California,with Aftershocks Removed,Poissonian? Bull. Seism. Soc. Am.,64(5):1363—1367.
Gibbons J.D. and Chakraborti S.,2003. Non-parametric Statistical Inference.Boca Raton,F(xiàn)la.:CRC Press.
Hwang Sungsoon,2005. Extending Spatial Hot Spot Detection Techniques to Temporal Dimensions. Environment and Planning B Planning and Design,36(5):865—882.
Kagan Y.Y. and Jackson D.D.,1994. Long-term Probabilistic Forecasting of Earthquakes. J. Geophys. Res.,99:13685—13700.
Knopoff L. and Gardner J.K.,1972. Higher Seismic Activity during Local Night on the Raw Worldwide Earthquake Catalog. Geophys. J.,28:311—313.
Lombardi A.M. and Marzocchi W.,2007. Evidence of Clustering and Nonstationarity in the Time Distribution of Large Worldwide Earthquakes. J. Geophys. Res.,112:B02303,doi:10.1029/2006JB004568.
Ogata Y. and Abe K.,1991. Some Statistical Features of the Long-term Variation of the Global and Regional Seismic Activity. Int. Stat. Rev.,59:139—161.
Pearson K.,1900. On the criterion that a system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling. Phil. Mag. Ser. 5,50:157—175.
Ripley B.,1976. The Second-order Analysis of Stationary Point Processes. Journal of Applied Probability,13(2):255—266.
Ripley B.,1977. Modelling Spatial Patterns (with discussion). Journal of the Royal Statistical Society,Series B 39:172—212.
Resenberg P.,1985. Second-Order Moment of Central California Seismicity,1969-1982. Journal of Geophysical Research,90(B7):5479—5495.
Utsu T.,1969. Aftershock and Earthquake Statistics(I):Some Parameters Which Characterize an Aftershock Sequence and Their Interrelations. Journal of the Faculty of Science,129—195.
Wang Q.,Jackson D.D. and Zhuang J.,2010. Are Spontaneous Earthquakes Stationary in California? J. Geophys.Res.,115:B08310,doi:10.1029/2009JB007031.
Zhuang J.,Ogata Y. and Vere-Jones D.,2002. Stochastic Declustering of Space-time Earthquake Occurrences. J. Am.Stat. Assoc.,97:369—380.