黃亦磊， 周仕勇*， 莊建倉

1 北京大學(xué)理論與應(yīng)用地球物理研究所， 北京　100871 2 日本統(tǒng)計(jì)數(shù)理研究所， 東京　190-8562

基于地震目錄估計(jì)完備震級方法的數(shù)值實(shí)驗(yàn)

黃亦磊1， 周仕勇1*， 莊建倉2

1 北京大學(xué)理論與應(yīng)用地球物理研究所， 北京100871 2 日本統(tǒng)計(jì)數(shù)理研究所， 東京190-8562

摘要本文將5種估計(jì)完備震級(magnitude of completeness,簡記為Mc)的方法運(yùn)用在三個(gè)不同模型產(chǎn)生的理論地震目錄上，進(jìn)而對比它們的優(yōu)缺點(diǎn).我們發(fā)現(xiàn)分段斜率中值分析法MBASS(The Median-based analysis of the segment slope)適用于不完備部分臺網(wǎng)探測地震能力隨震級變化快及監(jiān)測能力在時(shí)間上存在不均勻性(heterogeneity)的目錄中，但是要求目錄包含大量的地震事件，而b值穩(wěn)定法MBS(The Mc by b-value stability approach)則適合用于臺網(wǎng)探測地震能力隨震級減小衰減慢地震目錄，但是比較費(fèi)時(shí).最大曲率法MAXC(The Maximum Curvature technique)和擬合優(yōu)度測試法GFT(The Goodness-of-Fit Test)在使用時(shí)都低估Mc,需要加調(diào)整量.完整性震級范圍法EMR(Mc from Entire Magnitude Range)則一般給出比較穩(wěn)定、適中的Mc估計(jì)值.這種方法適用在地震數(shù)目少，且對地震丟失容忍度比較高的情況.在實(shí)踐中針對不同性質(zhì)的地震目錄，我們希望這項(xiàng)研究能幫助研究者選擇最合適估計(jì)完備震級Mc的方法，并指出了一些估計(jì)完備震級中應(yīng)當(dāng)避免的問題.

關(guān)鍵詞完備震級Mc; b值; 地震目錄完備性

1引言

地震目錄是地震活動(dòng)性(seismicity)分析、地震預(yù)測與地震危險(xiǎn)性評估的重要基礎(chǔ)資料(黃瑋瓊等, 1994; 劉杰等, 1996；徐偉進(jìn)和高孟潭, 2014).例如，利用地震目錄資料計(jì)算的b值可以作為反映區(qū)域構(gòu)造應(yīng)力狀態(tài)的一個(gè)指標(biāo)(Schorlemmer et al., 2005)，其他應(yīng)用例子有b值余震序列分析(e.g. Woessner et al., 2004; 譚毅陪等，2015)和地震活動(dòng)性中動(dòng)態(tài)觸發(fā)現(xiàn)象研究等(dynamic triggering)(Stein, 1999; Jia et al., 2012, 2014).研究表明(Ishimoto and Iida,1939; Gutenberg and Richter,1944)地震震級與地震發(fā)生頻率(Frequency Magnitude Distribution, 簡記為FMD)間滿足

(1)

公式中N(≥M)是大于等于震級M的地震事件累積頻率，a和b分別描述背景地震頻率和小地震與大地震的相對分布關(guān)系.從地震風(fēng)險(xiǎn)評估的角度來看，準(zhǔn)確和穩(wěn)健地估算b值具有重要意義.Aki(1965)推導(dǎo)了b值最似然估計(jì)：

(2)

這里〈M〉是地震目錄中所有震級大于等于Mc事件的平均震級，ΔM是震級劃分的最小分度值(magnitude bin, 一般取為0.1).從公式(2)可以看出，若欲從一個(gè)地震目錄中估計(jì)合理的b值，則依賴于Mc的正確選取.

地震臺在空間上的分布不均勻，臺網(wǎng)監(jiān)測能力在時(shí)間上也隨著地震臺數(shù)目和技術(shù)的改進(jìn)有所變化，另外各個(gè)臺站處理地震信號的方式也可能各不相同，所以在使用地震目錄前有必要評定地震目錄的質(zhì)量和一致性(consistency).在小震級范圍內(nèi)，實(shí)際記錄的地震發(fā)生累積頻率與震級關(guān)系對G-R定律的偏離可能由以下4個(gè)因素造成：(1) 地震事件震級太小，所產(chǎn)生的信號被埋沒在背景噪聲中；(2) 地震事件震級太小不足以被足夠多的臺站所記錄；(3) 臺網(wǎng)的工作人員認(rèn)為在一定震級以下的地震事件不被處理；(4) 在大地震之后，一些小地震可能無法從地震尾波中探測出來(Mignan and Woessner, 2012).由以上四個(gè)因素，我們可以看到，地震的震級越大，臺網(wǎng)對其記錄、定位并列入地震目錄的概率越大，這個(gè)探測概率會隨地震的震級增大逐漸收斂到100%.因而，地震學(xué)家們(Ogata and Katsura，1993, Woessner and Wiemer，2005，Iwata，2008)常用累積概率密度函數(shù)(cumulative probability function)刻畫逐漸隨震級增大而增加的地震探測概率.

Mc通常被定義成在一個(gè)時(shí)空范圍內(nèi)，地震能被臺網(wǎng)100%監(jiān)測到的最小震級(Rydelek and Sacks, 1989).而實(shí)際中，Mc還被用于從真實(shí)記錄的地震目錄中選取完備子目錄用于地震相關(guān)研究.在這篇文章中，我們關(guān)注一些常用估計(jì)Mc方法能否給出一個(gè)有效的Mc估計(jì)值，作為選取完備子目錄的依據(jù).假設(shè)地震發(fā)生是自相似過程，觀測到的地震累積發(fā)生頻率應(yīng)滿足G-R定律(公式(1)).在大震級處，地震發(fā)生頻率偏離線性G-R定律，可能原因是大地震數(shù)目太少引起的隨機(jī)波動(dòng)或是因?yàn)樘卣鞯卣瓞F(xiàn)象(characteristic earthquake phenomenon)(Schwartz and Coppersmith, 1984).而小震級處的偏離，如前文所述，則被解釋為臺網(wǎng)監(jiān)測能力不足.

公式(2)表明估計(jì)的b值依賴于選定的Mc(Aki, 1965).為了正確地估計(jì)b值，必須從地震目錄中挑選出一個(gè)完備的子目錄.由于地震探測概率在小震級處比較低，選取小的Mc值截取的子目錄可能不完備，導(dǎo)致錯(cuò)誤的參數(shù)估計(jì)，從而使分析有所偏離.一個(gè)比較安全的做法是選取一個(gè)足夠大的Mc值，但是這又將導(dǎo)致可用的數(shù)據(jù)變少.所以這項(xiàng)研究將重點(diǎn)權(quán)衡地震事件接近100%探測和可用于研究的數(shù)據(jù)量要求.

本項(xiàng)研究旨在通過數(shù)值實(shí)驗(yàn)來評估不同估計(jì)Mc方法的表現(xiàn).通過將這些方法運(yùn)用在人工合成的地震目錄(synthetic catalogs)上，我們希望得到這些方法在處理不同性質(zhì)地震目錄上的特點(diǎn)，從而為處理實(shí)際資料時(shí)提供一些理論指導(dǎo).目前，有兩類估計(jì)Mc的方法：第一類基于地震目錄(catalog-based method)，只使用地震目錄的數(shù)據(jù).第二類方法則是基于波形數(shù)據(jù)(waveform-based method)，這類方法(Gomberg, 1991; Sereno and Bratt, 1989)利用波形數(shù)據(jù)計(jì)算信噪比(signal-to-noise ratio)或者用震相識別的數(shù)據(jù)(phase-pick data)來確定Mc.絕大多數(shù)的估計(jì)完備震級的方法估計(jì)的是一個(gè)時(shí)空范圍的總體Mc，而由Schorlemmer和Woessner(2008)提出的基于概率的完備震級方法(Probability-based Magnitude of Completeness, PMC)和由Mignan等(2011)提出的貝葉斯完整性震級法(Bayesian Magnitude of Completeness, BMC)則可以給出Mc的時(shí)空分布特征，例如李智超和黃清華(2014)用PMC方法評估了首都圈的臺網(wǎng)監(jiān)測能力.由于基于地震目錄的方法相比基于波形的方法更加省時(shí)和容易操作，在實(shí)際操作中較多使用基于地震目錄的方法，例如在我國很多地區(qū)的完備震級估計(jì)中有運(yùn)用了此類方法(李志海等，2011；馮建剛等，2012).而且在大部分地震活動(dòng)性分析中，只需要用到一個(gè)時(shí)空范圍整體的Mc值，所以我們在這項(xiàng)研究中只研究基于地震目錄方法中5種比較流行的算法：

(1) 最大曲率法(The Maximum Curvature technique，簡記為MAXC)(Wiemer and Wyss, 2000)

這種方法選取震級頻率曲線(Frequency-magnitude Curve)中斜率最大值所對應(yīng)的震級作為Mc.在實(shí)際當(dāng)中，這個(gè)震級往往對應(yīng)非累積震級頻率分布(non-cumulative frequency magnitude distribution)中擁有最多地震數(shù)目的震級.

(2)擬合優(yōu)度測試法(The Goodness-of-Fit Test，簡記為GFT)(Wiemer and Wyss, 2000)

(3)

在公式(3)中，Bi和Si分別是每個(gè)震級區(qū)間(magnitude bin)實(shí)際觀測和理論計(jì)算的累積地震數(shù)目.Mc將取第一個(gè)使得R達(dá)到一定置信度(confidence level)的最小Mc0.置信度一般取成90%或95%.如果95%置信度可以達(dá)到，則不使用90%置信度.

(3)B值穩(wěn)定法(TheMcbyb-value stability approach,簡記為MBS)(Cao and Gao, 2002)

B值穩(wěn)定法是由Cao和Gao(2002)提出，這種方法將b值的穩(wěn)定性視為Mc0的函數(shù)，并假設(shè)b值將會隨著Mc0接近Mc而增大，當(dāng)Mc0≥Mc時(shí)，b值將保持不變.在Cao和Gao(2002)中，他們將b值穩(wěn)定性的標(biāo)準(zhǔn)設(shè)為0.03,然而這個(gè)值并不是在所有情況下都穩(wěn)定.所以Woessner和Wiemer(2005)使用b值的不確定度δb(Shi and Bolt, 1982)代替0.03，從而改善了這種方法.δb的表達(dá)式如下：

(4)

其中〈M〉是震級在Mc0以上事件震級的平均值，N是事件數(shù)目.Mc取使Δb=|bave-b|≤δb成立的最小Mc0.其中bave是Mc0相鄰的dM范圍內(nèi)每一個(gè)震級對應(yīng)b值的平均值，即

其中ΔM是震級分度值0.1，dM取0.5.

(4) 分段斜率中值分析法(The Median-based analysis of the segment slope,簡記為MBASS)(Amorese, 2007)

分段斜率中值分析由Amorese (2007)提出，是一種迭代尋找累積FMD中斜率序列多次改變點(diǎn)(multiple changes)的方法.此方法用Wilcoxon-Mann-Whitney (WMW) Test (Mann and Whitney, 1947; Wilcoxon, 1945)從迭代斜率序列中尋找FMD中的斜率不連續(xù)點(diǎn)，其中最主要的不連續(xù)點(diǎn)就對應(yīng)Mc，具體算法可以參見Amorese(2007).

(5) 完整性震級范圍法(Mcfrom Entire Magnitude Range，簡記為EMR)(Woessner and Wiemer, 2005)

完整性震級范圍法由Wossner和Wiemer(2005)提出，他們設(shè)立一個(gè)包含兩部分的模型.完備記錄部分用G-R定律描述，地震事件100%記錄.不完備部分用累積正態(tài)分布函數(shù)q(M|μ,σ)來表示臺網(wǎng)探測地震概率.在此模型中，震級M處的臺網(wǎng)探測地震概率q(M|μ, σ)可以表示為

(5)

在公式(5)中，μ是有50%概率探測到地震對應(yīng)的震級，用于描述一個(gè)臺網(wǎng)監(jiān)測地震的能力，σ是標(biāo)準(zhǔn)差，刻畫了不完備部分中臺網(wǎng)探測地震能力隨震級變化快慢程度，σ越大，探測地震的概率隨震級變化的速度越慢.由于Mc在探測密度函數(shù)中顯式表示，所以可以用最似然估計(jì)來估計(jì)Mc.

本項(xiàng)研究主要測試這5種基于地震目錄來估計(jì)Mc的方法在不同性質(zhì)地震目錄下的表現(xiàn).接下來我們就用3個(gè)模型來測試臺網(wǎng)探測地震概率隨震級變化快慢程度σ、地震數(shù)目和時(shí)空的不均勻性如何影響這些方法在Mc估計(jì)上的表現(xiàn).以上的5種方法將被用于由3個(gè)模型產(chǎn)生的人工地震目錄的檢測中.

2人工地震目錄的產(chǎn)生

人工地震目錄的產(chǎn)生依據(jù)的是Ogata和Katsura(1993)的模型(之后簡記為OK1993).模型用累積正態(tài)分布函數(shù)來描述臺網(wǎng)探測地震的能力，如下式所示：

(6)

其中所有參數(shù)的含義都和公式(5)中一樣.Woessner和Wiemmer(2005)與OK1993都使用了累積正態(tài)分布函數(shù)，但是前者只用在不完備部分，而后者用到了所有的震級范圍.在Woessner和Wiemmer(2005)中，他們使用了很多累積分布函數(shù)來描述臺網(wǎng)探測地震的概率函數(shù)，最后發(fā)現(xiàn)累積正態(tài)分布對實(shí)際數(shù)據(jù)擬合得最好.我們決定使用OK1993模型而不是Woessner和Wiemer(2005)的模型來產(chǎn)生人工地震目錄，因?yàn)榈卣鹛綔y概率隨著震級逐漸變化，概率不連續(xù)點(diǎn)在探測過程中是不合理的.在模擬人工地震目錄時(shí)，我們使用了拒絕法(Rejection Method)(Zhuang and Touati, 2015).在模擬的過程中，我們首先用理論地震分布概率密度和地震探測概率函數(shù)相乘推出觀測的地震概率密度函數(shù).基于由公式(1)給出的G-R定律，理論震級分布密度可以表示為

(7)

其中β=bln10.經(jīng)過歸一化的觀測地震概率密度可以寫成

(8)

我們將使用3個(gè)模型來模擬人工地震目錄.第一個(gè)模型(之后記為Model 1)，地震的概率密度服從公式(8)，q(M)中參數(shù)b=0.9,μ=1.5,σ=0.2.我們將b取成0.9，因?yàn)閷?shí)際中b值一般為0.8-1.1.σ取成0.2，與前人工作一致，其中μ描述的是臺網(wǎng)監(jiān)測地震的能力，取為1.5，表達(dá)的是若研究區(qū)域發(fā)生了一個(gè)1.5級的地震，臺網(wǎng)有50%的可能性對定位并放入地震目錄中.由于Mc是相對μ值大小的概念，因此μ的選取不影響基本結(jié)論.第二個(gè)模型(之后記為Model 2)，除σ=0.4之外，其他參數(shù)和Model 1一致.在Model 3中，一個(gè)從Model 1中產(chǎn)生的地震目錄和一個(gè)等數(shù)目震級在1.5以上被完全記錄的目錄混合在一起.Model 1和Model 2屬于同分布均勻模型，模型中σ越大，對應(yīng)著臺網(wǎng)對地震的監(jiān)測能力隨地震的震級變化越慢(當(dāng)σ趨于無窮時(shí)，表示臺網(wǎng)對各個(gè)震級的地震監(jiān)測能力一樣)，因此我們改變σ取值以檢測不同方法對完備震級估計(jì)在不同虛擬臺網(wǎng)(不同σ值)下的表現(xiàn).考慮到臺網(wǎng)探測地震能力不斷隨時(shí)間改變，使得地震目錄具有不均勻性，將Model 1和震級在1.5級以上完全記錄的目錄混合在一起得到的Model 3，一定程度上可以表示這種時(shí)間上的不均勻性，因此可以用Model 3 生成的地震目錄測試臺網(wǎng)探測能力隨時(shí)間變化對不同估計(jì)完備震級方法的影響.Models 1、2和3的概率密度分布圖如圖1(a—c)所示.

我們計(jì)算了3個(gè)模型在每100,500,1000,5000和10000個(gè)事件中期望缺震數(shù)為1所對應(yīng)的5個(gè)完備震級并列于表1.這些震級可以通過解下面的方程得出：

(9)

其中N是總地震事件數(shù)目，q(m)是公式(6)中的臺網(wǎng)探測地震的概率函數(shù).這些數(shù)值如表1所示，并在圖1中被畫出.在這項(xiàng)研究中，我們選擇每500個(gè)地震事件中期望缺震數(shù)為1對應(yīng)的震級為Mc的標(biāo)準(zhǔn).由于實(shí)際記錄的地震震級只保留一位小數(shù)，所以對于Model 1、Model 2和Model 3的Mc標(biāo)準(zhǔn)四舍五入之后得到1.9, 2.4和1.8.為了研究不同Mc截取地震子目錄對于b值和擬合度的影響，我們另外從每個(gè)模型中產(chǎn)生了100個(gè)包含200000個(gè)地震事件的目錄來計(jì)算Mc從1.5到3對應(yīng)的平均b值和由公式(3)所定義的擬合度.計(jì)算出的b值和擬合度，以及我們所選定的Mc標(biāo)準(zhǔn)均展現(xiàn)在圖2中.在之后的數(shù)值實(shí)驗(yàn)中，所有Mc估計(jì)結(jié)果都將與我們所選定的標(biāo)準(zhǔn)相比，讀者也可以根據(jù)自己對地震缺失的容忍度來選擇Mc的標(biāo)準(zhǔn)，并且對比這些實(shí)驗(yàn)結(jié)果.

為了研究地震數(shù)目對方法的影響，我們從每個(gè)模型產(chǎn)生3組數(shù)據(jù)，每組數(shù)據(jù)都包含1000個(gè)人工地震目錄，3組數(shù)據(jù)包含的地震數(shù)目不同，從第一組到第三組，地震數(shù)目分別為10000，50000和100000.接下來，我們將估計(jì)Mc的方法用于這些人工產(chǎn)生的地震目錄中，并且對比它們的表現(xiàn).

表1　Mc的可能理論標(biāo)準(zhǔn)值

注：*我們在這篇文章選擇的Mc的理論標(biāo)準(zhǔn)值.

圖1　三個(gè)模型的概率密度函數(shù)圖(a), (b), (c)分別展示了模型1到3的自然對數(shù)下的地震分布概率密度.在圖(a)、(b)中50%地震被探測的概率用帶星號的豎線表示.理論的Mc參考值在圖(a), (b)和(c)中用其他標(biāo)識5條豎線表示.在每100，500，1000，5000和100000個(gè)地震中期望缺震數(shù)為1的理論Mc值分別被從左到右的豎線表示.在這些豎線中，不帶點(diǎn)的虛線代表我們選擇的Mc參考值.Fig.1　Probability density function of three models for simulationNatural logarithm probability density (%) of Models 1 to 3 are shown respectively in (a), (b) and (c). The theoretical Mc are represented by five vertical lines and the 50% detection rate in (a) and (b) are represented by a vertical line marked with star. The theoretical Mc at which expectation of missing 1 event in every 100,500,1000,5000 and 100000 events are represented by vertical lines from left to right. Among the vertical lines, the dashed one without dot represent the criterion for Mc we choose in this study.

圖2　不同Mc對應(yīng)的平均b值和擬合度100次計(jì)算得到3個(gè)模型的平均b值和平均的擬合度分別如圖(a)、(b)所示.在每次運(yùn)算中使用了200000個(gè)地震事件，三條豎線是我們給三個(gè)模型選定的Mc的標(biāo)準(zhǔn).根據(jù)Mc選定地震子目錄對計(jì)算的b值和擬合度的影響可由此圖看出.Fig.2　Mean b value and goodness of fit corresponding to different McMean b value and goodness of fit due to different Mc by 100 times calculation of 3 models are respectively shown in Fig.2a and Fig.2b. In each calculation, 200000 events are used. The vertical lines are the Mc criterion we choose for this study. The influence of selecting a subset according to Mc on the estimation of b value and goodness of fit can be seen from this figure.

3計(jì)算Mc的實(shí)驗(yàn)結(jié)果

以上介紹的5種方法被分別用到從3個(gè)模型產(chǎn)生的9組人工地震目錄中.結(jié)果討論如下.

(1) MAXC.圖3帶圈的線表示MAXC的結(jié)果.從圖3(a—c)中可以看出對于Model 1的3組人工地震目錄，估計(jì)的Mc基本平均地分布在1.6和1.7，低于所選的標(biāo)準(zhǔn)1.9.由于地震數(shù)目的改變對估計(jì)Mc的分布影響很小，所以這種方法只要較少的地震數(shù)目就能達(dá)到穩(wěn)定的Mc估計(jì).當(dāng)標(biāo)準(zhǔn)差σ增加到0.4時(shí)，估計(jì)的Mc方差也相應(yīng)增大并且主要探測的Mc值下降至1.5或1.4(圖3(d—f))，遠(yuǎn)遠(yuǎn)低于2.4.從這里，我們可以看出，當(dāng)臺網(wǎng)探測地震概率隨震級變化變緩慢時(shí)，MAXC往往表現(xiàn)得更差.我們認(rèn)為這是由于地震事件震級更加分散分布所致，最多事件對應(yīng)的震級往往會更加偏離Mc.當(dāng)這種方法運(yùn)用到混合目錄時(shí)，估計(jì)的Mc一致為1.6，低于1.8(圖3(g—i)).我們總結(jié)在同分布的地震目錄中，MAXC會低估Mc，且低估的程度與臺網(wǎng)探測地震的概率隨震級改變的速率相關(guān)，但是這個(gè)方法能用很少的數(shù)據(jù)得出一個(gè)穩(wěn)定的結(jié)果.另外，這種方法也不適用于像Model 3產(chǎn)生的異質(zhì)性的目錄中.

我們所得出的結(jié)論與Mignan等(2011)認(rèn)為MAXC的低估是由于監(jiān)測臺網(wǎng)時(shí)空的不均勻性所造成的觀點(diǎn)有所不一致.因?yàn)樵谖覀兙鶆蛲植嫉哪Ｐ椭?例如Model 1和Model 2)，Mc同樣被低估.這種低估我們認(rèn)為是尋找最多地震事件對應(yīng)震級的算法所造成的.隨著震級的增大，地震被探測的概率增加，但是G-R定律則反映出隨震級增加，地震發(fā)生頻率下降，所以兩者綜合的結(jié)果可能導(dǎo)致Mc對應(yīng)的震級并不是擁有最多事件的震級.Model 2中Mc的低估程度和地震被探測的概率隨震級變化的速度有負(fù)相關(guān)關(guān)系也印證了這一點(diǎn).

(2) GFT.由于地震目錄的擬合度可以達(dá)到95%，所以我們選擇95%作為Mc置信度的標(biāo)準(zhǔn).Mc的估計(jì)值在圖3中用帶有三角形的線所表示，從圖3(a—e)我們可以看出，Model 1和Model 2的Mc估計(jì)值隨著地震數(shù)目的增加分別歸一到1.6和1.7.而從圖3(g—i)可知，Mc的估計(jì)值主要為1.6，少數(shù)分布在1.5.當(dāng)?shù)卣鹗录?shù)目增加時(shí)，結(jié)果變得更加穩(wěn)定，相比于MAXC，GFT在處理臺網(wǎng)探測地震能力隨震級變化比較快和有不均勻性的地震目錄時(shí)，更加有優(yōu)勢，但是同樣低估了每一個(gè)模型的Mc值.

低估的一個(gè)可能原因，我們認(rèn)為是這種方法取滿足置信度的最小震級作為Mc，因此忽略了所有滿足置信度的其他震級.高的置信度可能會給出一個(gè)更加接近真實(shí)Mc的結(jié)果.與MAXC相比，由于GFT運(yùn)用所有大于Mc0以上的事件，相比MAXC,GFT對臺網(wǎng)探測地震的概率隨震級變化緩慢和不均勻性的目錄更加具有抵抗力.

(3) MBS.我們用改進(jìn)后的MBS(Woessner and Wiemer, 2005)測試人工地震目錄，改進(jìn)的MBS使用Δb=|bave-b|≤δb作為標(biāo)準(zhǔn).標(biāo)準(zhǔn)差δb是用100次的有放回的重采樣計(jì)算得到，因此這個(gè)算法相對費(fèi)時(shí).相比于其他的方法，MBS擁有最高的Mc估計(jì)值，這種方法典型的特征是Mc估計(jì)值有個(gè)長尾巴.圖3(a—c)中，估計(jì)的Mc值隨著地震數(shù)目的增加，估計(jì)值眾數(shù)從1.8變化到1.9，而1.9恰好是我們所選的Mc標(biāo)準(zhǔn).從圖3(d—f)展示了用MBS估計(jì)的Mc對于臺網(wǎng)探測地震能力隨著震級改變速度是最敏感的.隨著地震數(shù)目的增加，Mc的主要估計(jì)值從2變化到了2.1，并有一定的概率取到2.2或2.3，這些值僅比2.4低一些.當(dāng)將MBS運(yùn)用到Model 3時(shí)，Mc的主要估計(jì)值為所選的Mc標(biāo)準(zhǔn)1.8，有時(shí)也到1.9(圖3(g—i)).

我們可以總結(jié)出MBS雖然是一種很費(fèi)時(shí)間的方法，但是也是相對保守的，對于Model 1和Model 3的Mc估計(jì)值基本上就是我們選定的Mc標(biāo)準(zhǔn)，而對于Model 2，Mc的估計(jì)值也是最高的，雖然比標(biāo)準(zhǔn)Mc要低了一些.這個(gè)方法在處理混合地震目錄時(shí)也表現(xiàn)很不錯(cuò)，但是這個(gè)方法需要相對大量地震數(shù)目才能實(shí)現(xiàn)穩(wěn)定的估計(jì).

(4) MBASS.MBASS也是一個(gè)依賴于地震事件數(shù)目的方法，需要相對大量的地震數(shù)目達(dá)到一個(gè)穩(wěn)定的Mc的估計(jì)值.圖3(a—c)展示了主要估計(jì)值隨著地震數(shù)目的增加從1.8變成1.9.當(dāng)臺網(wǎng)探測地震能力隨震級變化變緩時(shí)，可以很明顯地從圖3(d—f)看到Mc的估計(jì)結(jié)果強(qiáng)烈地依賴于地震數(shù)量，從1.6逐漸變成1.9.至于混合目錄，MBASS給出了和MBS類似的結(jié)果(圖3(g—i))，Mc估計(jì)值主要為1.8.

MBASS在臺網(wǎng)探測地震能力隨震級變化比較快和不均勻性存在時(shí)，比較適用.這個(gè)方法與MBS對比非常省時(shí)，而且在處理臺網(wǎng)探測地震能力隨震級變化較緩慢時(shí)，這個(gè)方法給出了除MBS以外最大的Mc估計(jì)值，雖然Mc也被低估.

圖3　5種方法的測試結(jié)果Model 1到Model 3所有人工合成地震目錄估計(jì)Mc的分布.(a)、(b)、(c)分別展示了Model 1中組1到3人工地震目錄的Mc.(d)、(e)、(f)展示了Model 2中組1到3的估計(jì)Mc分布，(g)、(h)、(i)則是Model 3中組1到3的Mc分布.豎直虛線則是我們選擇的Mc的標(biāo)準(zhǔn).MAXC, GFT, MBS, MBASS和EMR的結(jié)果分別用帶圈，三角形，菱形，加號和星號所代表.Fig.3　Estimation results of 5 methodsThe Mc distribution for syntheticcatalogs generated from Models 1 to 3 are shown in this figure. (a), (b) and (c) are estimated Mc distributions of Group 1 to 3 generated from Model 1. (d) to (f) shows the estimated Mc distribution of Group 1 to 3 for Model 2 and (g) to (i) for Model 3. The dashed vertical lines represent the criterion for Mc we choose. Results of MAXC, GFT, MBS, MBASS and EMR are respectively represented by lines marked with star, triangle, diamond, plus and star.

(5) EMR.對于Model 1，Mc的主要估計(jì)值一致是1.7(圖3(a—c)).當(dāng)臺網(wǎng)探測地震能力隨震級變化變緩時(shí)，Mc估計(jì)值分布也相應(yīng)變廣.Model 2的結(jié)果顯示這個(gè)方法對于臺網(wǎng)探測地震能力隨震級變化的速度不敏感，主要Mc估計(jì)值仍是1.7(圖3(d—f)).當(dāng)對混合目錄使用EMR時(shí)，估計(jì)的Mc主要為1.5或1.6(圖3(g—i)).這個(gè)方法也不需要大量的事件就能達(dá)到一個(gè)穩(wěn)定的估計(jì)，但是傾向于低估Mc.與其他4種方法相比，EMR給出了一個(gè)介于MAXC、GFT和MBASS、MBS的一個(gè)Mc估計(jì)值.我們沒有考察不完備部分不滿足累積正態(tài)分布的地震目錄，但是我們懷疑這個(gè)方法可能對這些目錄效果不是很好.

4討論

從圖3可以看出，當(dāng)模型中標(biāo)準(zhǔn)差σ變成0.4時(shí)，所有方法估計(jì)的Mc都偏低，所以我們建議當(dāng)臺網(wǎng)探測地震的概率隨著震級緩慢增加時(shí)，這樣記錄的地震目錄要小心處理.如果臺網(wǎng)探測地震的概率隨震級變化比較快，或者是具有不均勻性的地震目錄，當(dāng)?shù)卣饠?shù)目足夠多時(shí)，我們推薦MBASS.當(dāng)臺網(wǎng)探測地震概率隨震級變化緩慢且不考慮計(jì)算時(shí)間時(shí)，我們推薦MBS.

這項(xiàng)研究所用的人工地震目錄均是由OK1993產(chǎn)生,如前所述，這個(gè)模型有它的優(yōu)點(diǎn).但是在不完備部分由于和EMR模型有同樣的分布，所以可能在測試EMR時(shí)，結(jié)果說服力有待考慮.可以考慮用其他的模型來模擬人工地震目錄測試EMR，但是我們懷疑EMR的表現(xiàn)可能不那么好.由于選擇累積正態(tài)分布來描述地震目錄不完備部分并不是基于物理模型，可能存在其他更適合用于測試的模型.

我們關(guān)于EMR的結(jié)論和Woessner和Wiemer(2005)的有所不同.他們將EMR和MAXC、GFT和MBS進(jìn)行比較，發(fā)現(xiàn)EMR在估計(jì)Mc上表現(xiàn)得更好，但是我們在這項(xiàng)研究中發(fā)現(xiàn)EMR只是給出了一個(gè)一般介于MAXC & GFT和MBASS & MBS之間的結(jié)果.在他們合成的地震目錄中，他們設(shè)定震級大于或等于1.5的地震探測概率函數(shù)為1，我們嘗試了這種地震目錄，EMR確實(shí)表現(xiàn)很好，但是探測概率在一定震級以上為1會導(dǎo)致不合理的間斷點(diǎn).另外一個(gè)考慮是Woessner和Wiemer(2005)衡量了Mc的不確定度(uncertainty)與地震事件數(shù)量的關(guān)系，發(fā)現(xiàn)EMR需要最少事件就能達(dá)到一個(gè)很小的Mc的方差(variance)，在這一點(diǎn)上我們贊同他們，因?yàn)镋MR是一個(gè)基于4個(gè)參數(shù)估計(jì)的算法，所以結(jié)果的穩(wěn)定性相對高.我們認(rèn)為EMR方法低估Mc的可能原因之一是這個(gè)方法運(yùn)用Kolmogorov-Smirnov測試，設(shè)定0.05的顯著水平來接受或拒絕原假設(shè)——真實(shí)數(shù)據(jù)和計(jì)算的模型數(shù)據(jù)具有相同的分布.所以EMR方法可以看成是兩部分?jǐn)M合的GFT法，用分段函數(shù)擬合地震目錄的完備和不完備部分.所以我們認(rèn)為EMR低估Mc的原因和GFT可能一樣，與設(shè)定的顯著性水平相關(guān).另外一個(gè)可能原因是如果完備部分的擬合度很高，非完備部分的擬合度就可以相對低一些，從而導(dǎo)致Mc估計(jì)的偏低.EMR由于運(yùn)用了全部的地震事件，所以會給出一個(gè)相對GFT保守的Mc估計(jì).

有些方法需要一定數(shù)量的地震事件數(shù)量來達(dá)到一個(gè)穩(wěn)定的估計(jì)，而剩下的其他方法，例如MAXC、GFT和EMR，則相對不依賴地震事件數(shù)目.在數(shù)值實(shí)驗(yàn)中，地震事件數(shù)量可以人為操控，但是實(shí)際當(dāng)中并非如此.當(dāng)只有有限的地震數(shù)目時(shí)，有放回地對地震目錄重采樣(bootstrap)常常被用來增加Mc估計(jì)的穩(wěn)定性.但是有放回的重采樣和從一個(gè)模型中直接產(chǎn)生地震事件的區(qū)別我們沒有研究.正如文獻(xiàn)Woessner和Wiemer(2005)中表示，累積正態(tài)分布函數(shù)對實(shí)際地震目錄的不完備部分?jǐn)M合最好，那么重采樣的結(jié)果和直接從模型產(chǎn)生事件的區(qū)別就很小，則本文結(jié)論適用于重采樣得到的數(shù)據(jù)，這些結(jié)論則同樣適用于實(shí)踐當(dāng)中，可以作為選擇Mc估計(jì)方法的參考依據(jù).

5結(jié)論

Mc不僅指示一個(gè)臺網(wǎng)的監(jiān)測地震能力，而且常常被用來選擇一個(gè)真實(shí)地震目錄的子目錄，用于估計(jì)b值和a值等地震活動(dòng)性分析所用的參數(shù).用Mc截取地震子目錄對b值和擬合度的影響可以從圖2a和2b中看出.從圖2可以看出，如果從b與擬合度的角度來看，我們選擇的Mc標(biāo)準(zhǔn)是比較合理的.

我們將結(jié)果和表1我們所選擇的Mc標(biāo)準(zhǔn)進(jìn)行比較，結(jié)合討論，得出每種方法的優(yōu)缺點(diǎn)如下:

(1) MAXC是一種省時(shí)、簡單易行的估計(jì)Mc的方法，但是低估Mc.這種方法不需要大量的事件數(shù)目就能得到一個(gè)相對穩(wěn)定的結(jié)果，但是這個(gè)方法處理不了有不均勻性的地震目錄.在我們的數(shù)值實(shí)驗(yàn)中，當(dāng)臺網(wǎng)探測地震能力隨震級變化減緩時(shí)，Mc被更加低估，我們建議使用這種方法時(shí)，應(yīng)該根據(jù)臺網(wǎng)探測地震能力隨震級的變化快慢程度合理加上一個(gè)Mc調(diào)整量.

(2) GFT和MAXC一樣低估Mc，但是比MAXC對不均勻性和臺網(wǎng)探測地震能力隨震級的變化速率更有抵抗力.這種方法也很省時(shí)，并且對目錄所包含地震事件數(shù)要求不高.但對使用這種方法所估計(jì)的Mc值，為保險(xiǎn)起見，我們建議根據(jù)地震探測能力隨震級變化的快慢程度，適當(dāng)加上調(diào)整，例如在本文的實(shí)驗(yàn)中，當(dāng)?shù)卣鹛綔y能力隨震級變化快時(shí)，可以加上0.3的調(diào)整.

(3) 雖然MBS在臺網(wǎng)探測地震能力隨震級變化緩慢時(shí)略有低估Mc，但是在3個(gè)模型的人工地震目錄測試中表現(xiàn)得都很好.這個(gè)方法需要大量的地震事件來達(dá)到一個(gè)穩(wěn)定的Mc估計(jì)，耗時(shí)量也大于其他方法.在使用此方法時(shí)，若地震數(shù)目足夠多且不限制時(shí)間，則多次重復(fù)MBS計(jì)算的估計(jì)結(jié)果，即可用做Mc.

(4) MBASS 能準(zhǔn)確地估計(jì)地震數(shù)目多而且臺網(wǎng)探測地震能力隨震級變化較快或有不均勻性目錄的Mc.當(dāng)臺網(wǎng)探測地震能力隨震級變化緩慢時(shí)，這個(gè)方法表現(xiàn)得要比MBS差.所以，在處理臺網(wǎng)探測地震能力隨震級不是那么緩慢和有不均勻性的目錄時(shí)，我們推薦用MBASS對Mc進(jìn)行估計(jì).

(5) EMR方法低估了所有的情況下的Mc，但是一般給出比GFT和MAXC更大的但是小于MBASS和MBS給出的Mc估計(jì)值，并且也對臺網(wǎng)探測地震能力隨震級變化速率不敏感.由于EMR方法是基于4個(gè)參數(shù)的估計(jì)，這種方法也相對穩(wěn)定，并且比MBS以外方法耗時(shí).所以，當(dāng)臺網(wǎng)探測地震能力隨震級變化比較快，地震數(shù)目不是很足夠但要求結(jié)果穩(wěn)定，且對地震丟失容忍度比較高的話，我們建議用此方法，也可以加上一個(gè)小的Mc調(diào)整.

致謝我們在此感謝兩位匿名審稿人的中肯建議，感謝蘇黎世聯(lián)邦理工大學(xué)Woessner教授與我們討論計(jì)算Mc的程序.最后，我們感謝范文淵博士對這篇文章提出建設(shè)性的意見.

References

Aki K. 1965. Maximum likelihood estimate ofbin the formula and its confidence limits.Bull.Earthq.Res.Inst.,43: 237-239.

Amorese D. 2007. Applying a change-point detection method on frequency-magnitude distributions.Bull.Seismol.Soc.Am., 97(5): 1742-1749.Cao A M, Gao SS. 2002. Temporal variation of seismicb-values beneath northeastern Japan island arc.Geophys.Res.Lett., 29(9): 48-1-48-3.

Gomberg J. 1991. Seismicity and detection/location threshold in the southern Great Basin seismic network.JournalofGeophysicalResearch:SolidEarth(1978—2012), 96(B10): 16401-16414.

Gutenberg B, Richter C F. 1944. Frequency of earthquakes in California.Bull.Seismol.Soc.Am., 34(4): 185-188.

Ishimoto M, Iida K. 1939. Observations of earthquakes registered with the microseismograph constructed recently.Bull.Earthq.Res.Inst., 17: 443-478.

Iwata T. 2008. Low detection capability of global earthquakes after the occurrence of large earthquakes: Investigation of the Harvard CMT catalogue.Geophys.J.Int., 174(3): 849-856.

Jia K, Zhou S, Wang R. 2012. Stress interactions within the strong earthquake sequence from 2001 to 2010 in the Bayankala block of eastern Tibet.Bull.Seismol.Soc.Am., 102(5): 2157-2164.Jia K, Zhou S Y, Zhuang J C, et al. 2014. Possibility of the independence between the 2013 Lushan earthquake and the 2008 Wenchuan earthquake on Longmen Shan fault, Sichuan, China.Seismol.Res.Lett., 85(1): 60-67.

Mann H B, Whitney D R. 1947. On a test of whether one of two random variables is stochastically larger than the other.TheAnnalsofMathematicalStatistics, 18(1): 50-60.

Mignan A, Werner M J, Wiemer S, et al. 2011. Bayesian estimation of the spatially varying completeness magnitude of earthquake catalogs.Bull.Seismol.Soc.Am., 101(3): 1371-1385.Mignan A, Woessner J. 2012. Estimating the magnitude of completeness for earthquake catalogs.CommunityOnlineResourceforStatisticalSeismicityAnalysis, doi: 10.5078/corssa-00180805. Ogata Y, Katsura K. 1993. Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues.Geophys.J.Int., 113(3): 727-738.

Rydelek P A, Sacks I S. 1989. Testing the completeness of earthquake catalogues and the hypothesis of self-similarity.Nature, 337(6204): 251-253.Schorlemmer D, Wiemer S, Wyss M. 2005. Variations in earthquake-size distribution across different stress regimes.Nature, 437(7058): 539-542.SchorlemmerD,Woessner J. 2008. Probability of detecting an earthquake.Bull.Seismol.Soc.Am., 98(5):2103-2117.

Schwartz D P, Coppersmith K J. 1984. Fault behavior and characteristic earthquakes: Examples from the Wasatch and San Andreas fault zones.JournalofGeophysicalResearch:SolidEarth(1978—2012), 89(B7):5681-5698.

Sereno T J Jr, Bratt S R. 1989. Seismic detection capability at NORESS and implications for the detection threshold of a hypothetical network in the Soviet Union.JournalofGeophysicalResearch:SolidEarth(1978—2012), 94(B8):10397-10414.

Shi Y L, Bolt B A. 1982. The standard error of the magnitude-frequencybvalue.Bull.Seismol.Soc.Am., 72(5):1677-1687.Stein R S. 1999. The role of stress transfer in earthquake occurrence.Nature, 402(6762):605-609.

Wiemer S, Wyss M. 2000. Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the western United States, and Japan.Bull.Seismol.Soc.Am., 90(4):859-869.

Wilcoxon F. 1945. Individual comparisons by ranking methods.BiometricsBulletin, 1(6): 80-83.

Woessner J, Hauksson E, Wiemer S, et al. 2004.The 1997 Kagoshima (Japan) earthquake doublet: A quantitative analysis of aftershock rate changes.Geophys.Res.Lett., 31, doi: 10.1029/2003GL018858.

Woessner J,Wiemer S. 2005. Assessing the Quality of Earthquake Catalogues: Estimating the Magnitude of Completeness and Its Uncertainty.BSeismolSocAm, 95, 684-698.Zhuang, J C. and Touati S. 2015. Stochastic simulation of earthquake catalogs.CommunityOnlineResourceforStatisticalSeismicityAnalysis,doi:10.5078/corssa-43806322.

附中文參考文獻(xiàn)

馮建剛，蔣長勝，韓立波等. 2012. 甘肅測震臺監(jiān)測能力及地震目錄完整性分析. 地震學(xué)報(bào). 34(5): 646-658.

李智超，黃清華. 2014. 基于概率完備震級評估首都圈地震臺網(wǎng)檢測能力.地球物理學(xué)報(bào)，57(8): 2584-2593, doi: 10.6038/cjg20140818.

李志海，蔣長勝，黃瑜等. 2011. 新疆地區(qū)地震目錄最小完整性震級和臺網(wǎng)科學(xué)布局研究.地震學(xué)報(bào).33(6): 763-775.

劉杰，陳棋福，陳颙. 1996. 華北地區(qū)地震目錄完全性分析.地震, 16(1): 59-67.

譚毅培，曹井泉，陳繼鋒等. 2015. 2013年甘肅岷縣漳縣Ms6.6地震余震序列時(shí)域衰減特征分析.地球物理學(xué)報(bào)，58(9): 3222-3231, doi:10.6038/cjg20150917.

黃瑋瓊，李文香，曹學(xué)鋒. 1994. 中國大陸地震資料完整性研究之二. 地震學(xué)報(bào). 16(4): 423-432.

徐偉進(jìn)，高孟潭. 2014. 中國大陸及周緣地震目錄完整性統(tǒng)計(jì)分析.地球物理學(xué)報(bào)，57(9): 2802-2812，doi:10.6038/cjg20140907.

(本文編輯胡素芳)

Numerical tests on catalog-based methods to estimate magnitude of completeness

HUANG Yi-Lei1, ZHOU Shi-Yong1*, ZHUANG Jian-Cang2

1InstituteofTheoreticalandAppliedGeophysics,PekingUniversity,Beijing100871,China2InstituteofStatisticalMathematics, 10-3Midori-cho,Tachikawa,Tokyo190-8562,Japan

AbstractThis study compares five methods for estimating the completeness magnitude threshold of earthquake catalogs through applying them to synthetic catalogs generated from 3 different models. We have found that the Median-based analysis of the segment slope (MBASS) method is suitable for catalogs recorded by networks whose detection ability improves rapidly with magnitude and for those with temporal heterogeneity if the amount of earthquakes is large enough. The Mc by b-value stability approach (MBS) is optimal in dealing with catalogs recorded by networks whose detection ability improves slowly with magnitude, but it is time-consuming. The Maximum Curvature technique (MAXC)& The Goodness-of-Fit Test (GFT) method underestimate Mc and need an Mc criterion. The Mc from Entire Magnitude Range (EMR) method gives a moderate and stable Mc estimation. It is recommended when the amount of events is not large and the tolerance of missing events is relatively high.This study helps us to choose the optimal Mc estimation method in practice to cope with different earthquake catalogs and points out some potential problems caused by these methods.

KeywordsMagnitude threshold of completeness; b value; Completeness of catalogues

基金項(xiàng)目中國地震研究專項(xiàng)(201508009)和中國國家自然科學(xué)基金(41474033)聯(lián)合資助.

作者簡介黃亦磊，男，1993年出生，北京大學(xué)研究生在讀，主要從事統(tǒng)計(jì)地震學(xué)方面的研究. E-mail: huangyilei@pku.edu.cn *通訊作者周仕勇，男，1962年出生，教授，主要從事地震物理學(xué)方面的研究. E-mail: zsy@pku.edu.cn

doi:10.6038/cjg20160416 中圖分類號P315

收稿日期2015-10-15，2015-11-28收修定稿

黃亦磊， 周仕勇， 莊建倉. 2016. 基于地震目錄估計(jì)完備震級方法的數(shù)值實(shí)驗(yàn). 地球物理學(xué)報(bào)，59(4):1350-1358,doi:10.6038/cjg20160416.

Huang Y L, Zhou S Y, Zhuang J C. 2016. Numerical tests on catalog-based methods to estimate magnitude of completeness.ChineseJ.Geophys. (in Chinese),59(4):1350-1358,doi:10.6038/cjg20160416.