地球物理学报  2016, Vol. 59 Issue (4): 1350-1358   PDF    
引用本文
黄亦磊, 周仕勇, 庄建仓. 2016. 基于地震目录估计完备震级方法的数值实验. 地球物理学报, 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. Chinese J. Geophysics, 59(4): 1350-1358, doi: 10.6038/cjg20160416   
基于地震目录估计完备震级方法的数值实验
黄亦磊1, 周仕勇1, 庄建仓2    
1. 北京大学理论与应用地球物理研究所, 北京 100871;
2. 日本统计数理研究所, 东京 190-8562
收稿日期 2015-10-15, 2015-11-28 收修定稿
基金项目: 中国地震研究专项(201508009)和中国国家自然科学基金(41474033)联合资助.
作者简介: 黄亦磊,男,1993年出生,北京大学研究生在读,主要从事统计地震学方面的研究.E-mail:huangyilei@pku.edu.cn
通讯作者: 周仕勇,男,1962年出生,教授,主要从事地震物理学方面的研究.E-mail:zsy@pku.edu.cn
摘要: 本文将5种估计完备震级(magnitude of completeness,简记为Mc)的方法运用在三个不同模型产生的理论地震目录上,进而对比它们的优缺点.我们发现分段斜率中值分析法MBASS(The Median-based analysis of the segment slope)适用于不完备部分台网探测地震能力随震级变化快及监测能力在时间上存在不均匀性(heterogeneity)的目录中,但是要求目录包含大量的地震事件,而b值稳定法MBS(The Mc by b-value stability approach)则适合用于台网探测地震能力随震级减小衰减慢地震目录,但是比较费时.最大曲率法MAXC(The Maximum Curvature technique)和拟合优度测试法GFT(The Goodness-of-Fit Test)在使用时都低估Mc,需要加调整量.完整性震级范围法EMR(Mc from Entire Magnitude Range)则一般给出比较稳定、适中的Mc估计值.这种方法适用在地震数目少,且对地震丢失容忍度比较高的情况.在实践中针对不同性质的地震目录,我们希望这项研究能帮助研究者选择最合适估计完备震级Mc的方法,并指出了一些估计完备震级中应当避免的问题.
关键词: 完备震级Mc     b值     地震目录完备性    
Numerical tests on catalog-based methods to estimate magnitude of completeness
HUANG Yi-Lei1, ZHOU Shi-Yong1, ZHUANG Jian-Cang2    
1. Institute of Theoretical and Applied Geophysics, Peking University, Beijing 100871, China;
2. Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan
Abstract: This 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.
Key words: Magnitude threshold of completeness     b value     Completeness of catalogues    
1 引言

地震目录是地震活动性(seismicity)分析、地震预测与地震危险性评估的重要基础资料(黄玮琼等，1994; 刘杰等，1996；徐伟进和高孟潭，2014).例如，利用地震目录资料计算的b值可以作为反映区域构造应力状态的一个指标(Schorlemmer et al.，2005)，其他应用例子有b值余震序列分析(e.g. Woessner et al.，2004; 谭毅陪等，2015)和地震活动性中动态触发现象研究等(dynamic triggering)(Stein，1999; Jia et al.，2012，2014).研究表明(Ishimoto and Iida，1939; Gutenberg and Richter，1944)地震震级与地震发生频率(Frequency Magnitude Distribution，简记为FMD)间满足

公式中N(≥M)是大于等于震级M的地震事件累积频率，a和b分别描述背景地震频率和小地震与大地震的相对分布关系.从地震风险评估的角度来看，准确和稳健地估算b值具有重要意义.Aki(1965)推导了b值最似然估计：
这里〈M〉是地震目录中所有震级大于等于Mc事件的平均震级，ΔM是震级划分的最小分度值(magnitude bin，一般取为0.1).从公式(2)可以看出，若欲从一个地震目录中估计合理的b值，则依赖于Mc的正确选取.

地震台在空间上的分布不均匀，台网监测能力在时间上也随着地震台数目和技术的改进有所变化，另外各个台站处理地震信号的方式也可能各不相同，所以在使用地震目录前有必要评定地震目录的质量和一致性(consistency).在小震级范围内，实际记录的地震发生累积频率与震级关系对G-R定律的偏离可能由以下4个因素造成：(1)地震事件震级太小，所产生的信号被埋没在背景噪声中；(2)地震事件震级太小不足以被足够多的台站所记录；(3)台网的工作人员认为在一定震级以下的地震事件不被处理；(4)在大地震之后，一些小地震可能无法从地震尾波中探测出来(Mignan and Woessner，2012).由以上四个因素，我们可以看到，地震的震级越大，台网对其记录、定位并列入地震目录的概率越大，这个探测概率会随地震的震级增大逐渐收敛到100%.因而，地震学家们(Ogata and Katsura，1993，Woessner and Wiemer，2005，Iwata，2008)常用累积概率密度函数(cumulative probability function)刻画逐渐随震级增大而增加的地震探测概率.

Mc通常被定义成在一个时空范围内，地震能被台网100%监测到的最小震级(Rydelek and Sacks，1989).而实际中，Mc还被用于从真实记录的地震目录中选取完备子目录用于地震相关研究.在这篇文章中，我们关注一些常用估计Mc方法能否给出一个有效的Mc估计值，作为选取完备子目录的依据.假设地震发生是自相似过程，观测到的地震累积发生频率应满足G-R定律(公式(1)).在大震级处，地震发生频率偏离线性G-R定律，可能原因是大地震数目太少引起的随机波动或是因为特征地震现象(characteristic earthquake phenomenon)(Schwartz and Coppersmith，1984).而小震级处的偏离，如前文所述，则被解释为台网监测能力不足.

公式(2)表明估计的b值依赖于选定的Mc(Aki，1965).为了正确地估计b值，必须从地震目录中挑选出一个完备的子目录.由于地震探测概率在小震级处比较低，选取小的Mc值截取的子目录可能不完备，导致错误的参数估计，从而使分析有所偏离.一个比较安全的做法是选取一个足够大的Mc值，但是这又将导致可用的数据变少.所以这项研究将重点权衡地震事件接近100%探测和可用于研究的数据量要求.

本项研究旨在通过数值实验来评估不同估计Mc方法的表现.通过将这些方法运用在人工合成的地震目录(synthetic catalogs)上，我们希望得到这些方法在处理不同性质地震目录上的特点，从而为处理实际资料时提供一些理论指导.目前，有两类估计Mc的方法：第一类基于地震目录(catalog-based method)，只使用地震目录的数据.第二类方法则是基于波形数据(waveform-based method)，这类方法(Gomberg，1991; Sereno and Bratt，1989)利用波形数据计算信噪比(signal-to-noise ratio)或者用震相识别的数据(phase-pick data)来确定Mc.绝大多数的估计完备震级的方法估计的是一个时空范围的总体Mc，而由Schorlemmer和Woessner(2008)提出的基于概率的完备震级方法(Probability-based Magnitude of Completeness，PMC)和由Mignan等(2011)提出的贝叶斯完整性震级法(Bayesian Magnitude of Completeness，BMC)则可以给出Mc的时空分布特征，例如李智超和黄清华(2014)用PMC方法评估了首都圈的台网监测能力.由于基于地震目录的方法相比基于波形的方法更加省时和容易操作，在实际操作中较多使用基于地震目录的方法，例如在我国很多地区的完备震级估计中有运用了此类方法(李志海等，2011；冯建刚等，2012).而且在大部分地震活动性分析中，只需要用到一个时空范围整体的Mc值，所以我们在这项研究中只研究基于地震目录方法中5种比较流行的算法：

(1)最大曲率法(The Maximum Curvature technique，简记为MAXC)(Wiemer and Wyss，2000)

这种方法选取震级频率曲线(Frequency-magnitude Curve)中斜率最大值所对应的震级作为Mc.在实际当中，这个震级往往对应非累积震级频率分布(non-cumulative frequency magnitude distribution)中拥有最多地震数目的震级.

(2)拟合优度测试法(The Goodness-of-Fit Test，简记为GFT)(Wiemer and Wyss，2000)

拟合优度测试法测度的是拟合的FMD与实际观测的FMD之间的差异，当目录的不完备部分被包含时，拟合的FMD与实际观测FMD差异会增大，也就是拟合度(goodness of fit)低.在操作过程中，我们对于每一个震级下限，都采用最似然估计的方法，计算在震级下限以上事件对应的a与b值，从而构造出一个理论上严格服从G-R定律的FMD，定义相对误差R:

在公式(3)中，Bi和Si分别是每个震级区间(magnitude bin)实际观测和理论计算的累积地震数目.Mc将取第一个使得R达到一定置信度(confidence level)的最小Mc0.置信度一般取成90%或95%.如果95%置信度可以达到，则不使用90%置信度.

(3)B值稳定法(The Mc by b-value stability approach，简记为MBS)(Cao and Gao，2002)

B值稳定法是由Cao和Gao(2002)提出，这种方法将b值的稳定性视为Mc0的函数，并假设b值将会随着Mc0接近Mc而增大，当Mc0≥Mc时，b值将保持不变.在Cao和Gao(2002)中，他们将b值稳定性的标准设为0.03，然而这个值并不是在所有情况下都稳定.所以Woessner和Wiemer(2005)使用b值的不确定度δb(Shi and Bolt，1982)代替0.03，从而改善了这种方法.δb的表达式如下：

其中〈M〉是震级在Mc0以上事件震级的平均值，N是事件数目.Mc取使Δb=｜bave-b｜≤δb成立的最小Mc0.其中bave是Mc0相邻的dM范围内每一个震级对应b值的平均值，即
其中ΔM是震级分度值0.1，dM取0.5.

(4)分段斜率中值分析法(The Median-based analysis of the segment slope，简记为MBASS)(Amorese，2007)

分段斜率中值分析由Amorese(2007)提出，是一种迭代寻找累积FMD中斜率序列多次改变点(multiple changes)的方法.此方法用Wilcoxon-Mann-Whitney(WMW)Test(Mann and Whitney，1947; Wilcoxon，1945)从迭代斜率序列中寻找FMD中的斜率不连续点，其中最主要的不连续点就对应Mc，具体算法可以参见Amorese(2007).

(5)完整性震级范围法(Mc from Entire Magnitude Range，简记为EMR)(Woessner and Wiemer，2005)

完整性震级范围法由Wossner和Wiemer(2005)提出，他们设立一个包含两部分的模型.完备记录部分用G-R定律描述，地震事件100%记录.不完备部分用累积正态分布函数q(M｜μ，σ)来表示台网探测地震概率.在此模型中，震级M处的台网探测地震概率q(M｜μ，σ)可以表示为

在公式(5)中，μ是有50%概率探测到地震对应的震级，用于描述一个台网监测地震的能力，σ是标准差，刻画了不完备部分中台网探测地震能力随震级变化快慢程度，σ越大，探测地震的概率随震级变化的速度越慢.由于Mc在探测密度函数中显式表示，所以可以用最似然估计来估计Mc.

本项研究主要测试这5种基于地震目录来估计Mc的方法在不同性质地震目录下的表现.接下来我们就用3个模型来测试台网探测地震概率随震级变化快慢程度σ、地震数目和时空的不均匀性如何影响这些方法在Mc估计上的表现.以上的5种方法将被用于由3个模型产生的人工地震目录的检测中.

2 人工地震目录的产生

人工地震目录的产生依据的是Ogata和Katsura(1993)的模型(之后简记为OK1993).模型用累积正态分布函数来描述台网探测地震的能力，如下式所示：

其中所有参数的含义都和公式(5)中一样.Woessner和Wiemmer(2005)与OK1993都使用了累积正态分布函数，但是前者只用在不完备部分，而后者用到了所有的震级范围.在Woessner和Wiemmer(2005)中，他们使用了很多累积分布函数来描述台网探测地震的概率函数，最后发现累积正态分布对实际数据拟合得最好.我们决定使用OK1993模型而不是Woessner和Wiemer(2005)的模型来产生人工地震目录，因为地震探测概率随着震级逐渐变化，概率不连续点在探测过程中是不合理的.在模拟人工地震目录时，我们使用了拒绝法(Rejection Method)(Zhuang and Touati，2015).在模拟的过程中，我们首先用理论地震分布概率密度和地震探测概率函数相乘推出观测的地震概率密度函数.基于由公式(1)给出的G-R定律，理论震级分布密度可以表示为
其中β=bln10.经过归一化的观测地震概率密度可以写成

我们将使用3个模型来模拟人工地震目录.第一个模型(之后记为Model 1)，地震的概率密度服从公式(8)，q(M)中参数b=0.9，μ=1.5，σ=0.2.我们将b取成0.9，因为实际中b值一般为0.8-1.1.σ取成0.2，与前人工作一致，其中μ描述的是台网监测地震的能力，取为1.5，表达的是若研究区域发生了一个1.5级的地震，台网有50%的可能性对定位并放入地震目录中.由于Mc是相对μ值大小的概念，因此μ的选取不影响基本结论.第二个模型(之后记为Model 2)，除σ=0.4之外，其他参数和Model 1一致.在Model 3中，一个从Model 1中产生的地震目录和一个等数目震级在1.5以上被完全记录的目录混合在一起.Model 1和Model 2属于同分布均匀模型，模型中σ越大，对应着台网对地震的监测能力随地震的震级变化越慢(当σ趋于无穷时，表示台网对各个震级的地震监测能力一样)，因此我们改变σ取值以检测不同方法对完备震级估计在不同虚拟台网(不同σ值)下的表现.考虑到台网探测地震能力不断随时间改变，使得地震目录具有不均匀性，将Model 1和震级在1.5级以上完全记录的目录混合在一起得到的Model 3，一定程度上可以表示这种时间上的不均匀性，因此可以用Model 3 生成的地震目录测试台网探测能力随时间变化对不同估计完备震级方法的影响.Models 1、2和3的概率密度分布图如图 1(a—c)所示.

图 1 三个模型的概率密度函数
图(a), (b), (c)分别展示了模型1到3的自然对数下的地震分布概率密度.在图(a)、(b)中50%地震被探测的概率用带星号的竖线表示.理论的Mc参考值在图(a), (b)和(c)中用其他标识5条竖线表示.在每100，500，1000，5000和100000个地震中期望缺震数为1的理论Mc值分别被从左到右的竖线表示.在这些竖线中，不带点的虚线代表我们选择的Mc参考值. Fig. 1 Probability density function of three models for simulation
Natural 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.

我们计算了3个模型在每100，500，1000，5000和10000个事件中期望缺震数为1所对应的5个完备震级并列于表 1.这些震级可以通过解下面的方程得出：

其中N是总地震事件数目，q(m)是公式(6)中的台网探测地震的概率函数.这些数值如表 1所示，并在图 1中被画出.在这项研究中，我们选择每500个地震事件中期望缺震数为1对应的震级为Mc的标准.由于实际记录的地震震级只保留一位小数，所以对于Model 1、Model 2和Model 3的Mc标准四舍五入之后得到1.9，2.4和1.8.为了研究不同Mc截取地震子目录对于b值和拟合度的影响，我们另外从每个模型中产生了100个包含200000个地震事件的目录来计算Mc从1.5到3对应的平均b值和由公式(3)所定义的拟合度.计算出的b值和拟合度，以及我们所选定的Mc标准均展现在图 2中.在之后的数值实验中，所有Mc估计结果都将与我们所选定的标准相比，读者也可以根据自己对地震缺失的容忍度来选择Mc的标准，并且对比这些实验结果.
表 1 Mc的可能理论标准值 Table 1 Possible theoretical criterion for Mc

图 2 不同Mc对应的平均b值和拟合度
100次计算得到3个模型的平均b值和平均的拟合度分别如图(a)、(b)所示.在每次运算中使用了200000个地震事件，三条竖线是 我们给三个模型选定的Mc的标准.根据Mc选定地震子目录对计算的b值和拟合度的影响可由此图看出. Fig. 2 Mean b value and goodness of fit corresponding to different Mc
Mean 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组数据，每组数据都包含1000个人工地震目录，3组数据包含的地震数目不同，从第一组到第三组，地震数目分别为10000，50000和100000.接下来，我们将估计Mc的方法用于这些人工产生的地震目录中，并且对比它们的表现.

3 计算Mc的实验结果

以上介绍的5种方法被分别用到从3个模型产生的9组人工地震目录中.结果讨论如下.

(1)MAXC.图 3带圈的线表示MAXC的结果.从图 3(a—c)中可以看出对于Model 1的3组人工地震目录，估计的Mc基本平均地分布在1.6和1.7，低于所选的标准1.9.由于地震数目的改变对估计Mc的分布影响很小，所以这种方法只要较少的地震数目就能达到稳定的Mc估计.当标准差σ增加到0.4时，估计的Mc方差也相应增大并且主要探测的Mc值下降至1.5或1.4(图 3(d—f))，远远低于2.4.从这里，我们可以看出，当台网探测地震概率随震级变化变缓慢时，MAXC往往表现得更差.我们认为这是由于地震事件震级更加分散分布所致，最多事件对应的震级往往会更加偏离Mc.当这种方法运用到混合目录时，估计的Mc一致为1.6，低于1.8(图 3(g—i)).我们总结在同分布的地震目录中，MAXC会低估Mc，且低估的程度与台网探测地震的概率随震级改变的速率相关，但是这个方法能用很少的数据得出一个稳定的结果.另外，这种方法也不适用于像Model 3产生的异质性的目录中.

图 3 5种方法的测试结果
Model 1到Model 3所有人工合成地震目录估计Mc的分布.(a)、(b)、(c)分别展示了Model 1中组1到3人工地震目录的Mc.(d)、(e)、(f)展示了Model 2中组1到3的估计Mc分布，(g)、(h)、(i)则是Model 3中组1到3的Mc分布.竖直虚线则是我们选择的Mc的标准.MAXC, GFT, MBS, MBASS和EMR的结果分别用带圈，三角形，菱形，加号和星号所代表. Fig. 3 Estimation results of 5 methods
The 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.

我们所得出的结论与Mignan等(2011)认为MAXC的低估是由于监测台网时空的不均匀性所造成的观点有所不一致.因为在我们均匀同分布的模型中(例如Model 1和Model 2)，Mc同样被低估.这种低估我们认为是寻找最多地震事件对应震级的算法所造成的.随着震级的增大，地震被探测的概率增加，但是G-R定律则反映出随震级增加，地震发生频率下降，所以两者综合的结果可能导致Mc对应的震级并不是拥有最多事件的震级.Model 2中Mc的低估程度和地震被探测的概率随震级变化的速度有负相关关系也印证了这一点.

(2)GFT.由于地震目录的拟合度可以达到95%，所以我们选择95%作为Mc置信度的标准.Mc的估计值在图 3中用带有三角形的线所表示，从图 3(a—e)我们可以看出，Model 1和Model 2的Mc估计值随着地震数目的增加分别归一到1.6和1.7.而从图 3(g—i)可知，Mc的估计值主要为1.6，少数分布在1.5.当地震事件数目增加时，结果变得更加稳定，相比于MAXC，GFT在处理台网探测地震能力随震级变化比较快和有不均匀性的地震目录时，更加有优势，但是同样低估了每一个模型的Mc值.

低估的一个可能原因，我们认为是这种方法取满足置信度的最小震级作为Mc，因此忽略了所有满足置信度的其他震级.高的置信度可能会给出一个更加接近真实Mc的结果.与MAXC相比，由于GFT运用所有大于Mc0以上的事件，相比MAXC，GFT对台网探测地震的概率随震级变化缓慢和不均匀性的目录更加具有抵抗力.

(3)MBS.我们用改进后的MBS(Woessner and Wiemer，2005)测试人工地震目录，改进的MBS使用Δb=｜bave-b｜≤δb作为标准.标准差δb是用100次的有放回的重采样计算得到，因此这个算法相对费时.相比于其他的方法，MBS拥有最高的Mc估计值，这种方法典型的特征是Mc估计值有个长尾巴.图 3(a—c)中，估计的Mc值随着地震数目的增加，估计值众数从1.8变化到1.9，而1.9恰好是我们所选的Mc标准.从图 3(d—f)展示了用MBS估计的Mc对于台网探测地震能力随着震级改变速度是最敏感的.随着地震数目的增加，Mc的主要估计值从2变化到了2.1，并有一定的概率取到2.2或2.3，这些值仅比2.4低一些.当将MBS运用到Model 3时，Mc的主要估计值为所选的Mc标准1.8，有时也到1.9(图 3(g—i)).

我们可以总结出MBS虽然是一种很费时间的方法，但是也是相对保守的，对于Model 1和Model 3的Mc估计值基本上就是我们选定的Mc标准，而对于Model 2，Mc的估计值也是最高的，虽然比标准Mc要低了一些.这个方法在处理混合地震目录时也表现很不错，但是这个方法需要相对大量地震数目才能实现稳定的估计.

(4)MBASS.MBASS也是一个依赖于地震事件数目的方法，需要相对大量的地震数目达到一个稳定的Mc的估计值.图 3(a—c)展示了主要估计值随着地震数目的增加从1.8变成1.9.当台网探测地震能力随震级变化变缓时，可以很明显地从图 3(d—f)看到Mc的估计结果强烈地依赖于地震数量，从1.6逐渐变成1.9.至于混合目录，MBASS给出了和MBS类似的结果(图 3(g—i))，Mc估计值主要为1.8.

MBASS在台网探测地震能力随震级变化比较快和不均匀性存在时，比较适用.这个方法与MBS对比非常省时，而且在处理台网探测地震能力随震级变化较缓慢时，这个方法给出了除MBS以外最大的Mc估计值，虽然Mc也被低估.

(5)EMR.对于Model 1，Mc的主要估计值一致是1.7(图 3(a—c)).当台网探测地震能力随震级变化变缓时，Mc估计值分布也相应变广.Model 2的结果显示这个方法对于台网探测地震能力随震级变化的速度不敏感，主要Mc估计值仍是1.7(图 3(d—f)).当对混合目录使用EMR时，估计的Mc主要为1.5或1.6(图 3(g—i)).这个方法也不需要大量的事件就能达到一个稳定的估计，但是倾向于低估Mc.与其他4种方法相比，EMR给出了一个介于MAXC、GFT和MBASS、MBS的一个Mc估计值.我们没有考察不完备部分不满足累积正态分布的地震目录，但是我们怀疑这个方法可能对这些目录效果不是很好.

4 讨论

从图 3可以看出，当模型中标准差σ变成0.4时，所有方法估计的Mc都偏低，所以我们建议当台网探测地震的概率随着震级缓慢增加时，这样记录的地震目录要小心处理.如果台网探测地震的概率随震级变化比较快，或者是具有不均匀性的地震目录，当地震数目足够多时，我们推荐MBASS.当台网探测地震概率随震级变化缓慢且不考虑计算时间时，我们推荐MBS.

这项研究所用的人工地震目录均是由OK1993产生，如前所述，这个模型有它的优点.但是在不完备部分由于和EMR模型有同样的分布，所以可能在测试EMR时，结果说服力有待考虑.可以考虑用其他的模型来模拟人工地震目录测试EMR，但是我们怀疑EMR的表现可能不那么好.由于选择累积正态分布来描述地震目录不完备部分并不是基于物理模型，可能存在其他更适合用于测试的模型.

我们关于EMR的结论和Woessner和Wiemer(2005)的有所不同.他们将EMR和MAXC、GFT和MBS进行比较，发现EMR在估计Mc上表现得更好，但是我们在这项研究中发现EMR只是给出了一个一般介于MAXC & GFT和MBASS & MBS之间的结果.在他们合成的地震目录中，他们设定震级大于或等于1.5的地震探测概率函数为1，我们尝试了这种地震目录，EMR确实表现很好，但是探测概率在一定震级以上为1会导致不合理的间断点.另外一个考虑是Woessner和Wiemer(2005)衡量了Mc的不确定度(uncertainty)与地震事件数量的关系，发现EMR需要最少事件就能达到一个很小的Mc的方差(variance)，在这一点上我们赞同他们，因为EMR是一个基于4个参数估计的算法，所以结果的稳定性相对高.我们认为EMR方法低估Mc的可能原因之一是这个方法运用Kolmogorov-Smirnov测试，设定0.05的显著水平来接受或拒绝原假设——真实数据和计算的模型数据具有相同的分布.所以EMR方法可以看成是两部分拟合的GFT法，用分段函数拟合地震目录的完备和不完备部分.所以我们认为EMR低估Mc的原因和GFT可能一样，与设定的显著性水平相关.另外一个可能原因是如果完备部分的拟合度很高，非完备部分的拟合度就可以相对低一些，从而导致Mc估计的偏低.EMR由于运用了全部的地震事件，所以会给出一个相对GFT保守的Mc估计.

有些方法需要一定数量的地震事件数量来达到一个稳定的估计，而剩下的其他方法，例如MAXC、GFT和EMR，则相对不依赖地震事件数目.在数值实验中，地震事件数量可以人为操控，但是实际当中并非如此.当只有有限的地震数目时，有放回地对地震目录重采样(bootstrap)常常被用来增加Mc估计的稳定性.但是有放回的重采样和从一个模型中直接产生地震事件的区别我们没有研究.正如文献Woessner和Wiemer(2005)中表示，累积正态分布函数对实际地震目录的不完备部分拟合最好，那么重采样的结果和直接从模型产生事件的区别就很小，则本文结论适用于重采样得到的数据，这些结论则同样适用于实践当中，可以作为选择Mc估计方法的参考依据.

5 结论

Mc不仅指示一个台网的监测地震能力，而且常常被用来选择一个真实地震目录的子目录，用于估计b值和a值等地震活动性分析所用的参数.用Mc截取地震子目录对b值和拟合度的影响可以从图 2a和2b中看出.从图 2可以看出，如果从b与拟合度的角度来看，我们选择的Mc标准是比较合理的.

我们将结果和表 1我们所选择的Mc标准进行比较，结合讨论，得出每种方法的优缺点如下:

(1)MAXC是一种省时、简单易行的估计Mc的方法，但是低估Mc.这种方法不需要大量的事件数目就能得到一个相对稳定的结果，但是这个方法处理不了有不均匀性的地震目录.在我们的数值实验中，当台网探测地震能力随震级变化减缓时，Mc被更加低估，我们建议使用这种方法时，应该根据台网探测地震能力随震级的变化快慢程度合理加上一个Mc调整量.

(2)GFT和MAXC一样低估Mc，但是比MAXC对不均匀性和台网探测地震能力随震级的变化速率更有抵抗力.这种方法也很省时，并且对目录所包含地震事件数要求不高.但对使用这种方法所估计的Mc值，为保险起见，我们建议根据地震探测能力随震级变化的快慢程度，适当加上调整，例如在本文的实验中，当地震探测能力随震级变化快时，可以加上0.3的调整.

(3)虽然MBS在台网探测地震能力随震级变化缓慢时略有低估Mc，但是在3个模型的人工地震目录测试中表现得都很好.这个方法需要大量的地震事件来达到一个稳定的Mc估计，耗时量也大于其他方法.在使用此方法时，若地震数目足够多且不限制时间，则多次重复MBS计算的估计结果，即可用做Mc.

(4)MBASS 能准确地估计地震数目多而且台网探测地震能力随震级变化较快或有不均匀性目录的Mc.当台网探测地震能力随震级变化缓慢时，这个方法表现得要比MBS差.所以，在处理台网探测地震能力随震级不是那么缓慢和有不均匀性的目录时，我们推荐用MBASS对Mc进行估计.

(5)EMR方法低估了所有的情况下的Mc，但是一般给出比GFT和MAXC更大的但是小于MBASS和MBS给出的Mc估计值，并且也对台网探测地震能力随震级变化速率不敏感.由于EMR方法是基于4个参数的估计，这种方法也相对稳定，并且比MBS以外方法耗时.所以，当台网探测地震能力随震级变化比较快，地震数目不是很足够但要求结果稳定，且对地震丢失容忍度比较高的话，我们建议用此方法，也可以加上一个小的Mc调整.

致谢 我们在此感谢两位匿名审稿人的中肯建议，感谢苏黎世联邦理工大学Woessner教授与我们讨论计算Mc的程序.最后，我们感谢范文渊博士对这篇文章提出建设性的意见.
参考文献
[1] Aki K. 1965. Maximum likelihood estimate of b in the formula and its confidence limits.Bull. Earthq.Res.Inst.,43: 237-239.
[2] Amorese D. 2007. Applying a change-point detection method on frequency-magnitude distributions. Bull. Seismol. Soc. Am., 97(5): 1742-1749.
[3] Cao A M, Gao SS. 2002. Temporal variation of seismic b-values beneath northeastern Japan island arc. Geophys. Res. Lett., 29(9): 48-1-48-3.
[4] Gomberg J. 1991. Seismicity and detection/location threshold in the southern Great Basin seismic network. Journal of Geophysical Research: Solid Earth (1978—2012), 96(B10): 16401-16414.
[5] Gutenberg B, Richter C F. 1944. Frequency of earthquakes in California. Bull. Seismol. Soc. Am., 34(4): 185-188.
[6] Ishimoto M, Iida K. 1939. Observations of earthquakes registered with the microseismograph constructed recently. Bull. Earthq. Res. Inst., 17: 443-478.
[7] 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.
[8] 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.
[9] 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.
[10] Mann H B, Whitney D R. 1947. On a test of whether one of two random variables is stochastically larger than the other. The Annals of Mathematical Statistics, 18(1): 50-60.
[11] 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.
[12] Mignan A, Woessner J. 2012. Estimating the magnitude of completeness for earthquake catalogs. Community Online Resource for Statistical Seismicity Analysis, doi: 10.5078/corssa-00180805.
[13] 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.
[14] Rydelek P A, Sacks I S. 1989. Testing the completeness of earthquake catalogues and the hypothesis of self-similarity. Nature, 337(6204): 251-253.
[15] Schorlemmer D, Wiemer S, Wyss M. 2005. Variations in earthquake-size distribution across different stress regimes. Nature, 437(7058): 539-542.
[16] Schorlemmer D,Woessner J. 2008. Probability of detecting an earthquake. Bull. Seismol. Soc. Am., 98(5):2103-2117.
[17] Schwartz D P, Coppersmith K J. 1984. Fault behavior and characteristic earthquakes: Examples from the Wasatch and San Andreas fault zones. Journal of Geophysical Research: Solid Earth (1978—2012), 89(B7):5681-5698.
[18] 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. Journal of Geophysical Research: Solid Earth (1978—2012), 94(B8):10397-10414.
[19] Shi Y L, Bolt B A. 1982. The standard error of the magnitude-frequency b value. Bull. Seismol. Soc. Am., 72(5):1677-1687.
[20] Stein R S. 1999. The role of stress transfer in earthquake occurrence. Nature, 402(6762):605-609.
[21] 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.
[22] Wilcoxon F. 1945. Individual comparisons by ranking methods. Biometrics Bulletin, 1(6): 80-83.
[23] 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.
[24] Woessner J,Wiemer S. 2005. Assessing the Quality of Earthquake Catalogues: Estimating the Magnitude of Completeness and Its Uncertainty. B SeismolSoc Am, 95, 684-698.
[25] Zhuang, J C. and Touati S. 2015. Stochastic simulation of earthquake catalogs. Community Online Resource for Statistical Seismicity Analysis,doi:10.5078/corssa-43806322.
[26] 冯建刚,蒋长胜,韩立波等. 2012. 甘肃测震台监测能力及地震目录完整性分析. 地震学报. 34(5): 646-658.
[27] 李智超,黄清华. 2014. 基于概率完备震级评估首都圈地震台网检测能力.地球物理学报,57(8): 2584-2593, doi: 10.6038/cjg20140818.
[28] 李志海,蒋长胜,黄瑜等. 2011. 新疆地区地震目录最小完整性震级和台网科学布局研究.地震学报.33(6): 763-775.
[29] 刘杰,陈棋福,陈颙. 1996. 华北地区地震目录完全性分析.地震, 16(1): 59-67.
[30] 谭毅培,曹井泉,陈继锋等. 2015. 2013年甘肃岷县漳县Ms6.6地震余震序列时域衰减特征分析.地球物理学报,58(9): 3222-3231, doi:10.6038/cjg20150917.
[31] 黄玮琼,李文香,曹学锋. 1994. 中国大陆地震资料完整性研究之二. 地震学报. 16(4): 423-432.
[32] 徐伟进,高孟潭. 2014. 中国大陆及周缘地震目录完整性统计分析.地球物理学报,57(9): 2802-2812,doi:10.6038/cjg20140907.