Bayesian estimation of completeness magnitude in the Capital Circle Region
-
摘要:
地震目录的完整性评估是地震活动性分析的基础性工作，常用的基于地震目录的评估方法未能考虑台站信息，对于少震地区无法给出评估结果，并且受到主观选取的计算参数的影响。本文采用基于贝叶斯统计的完整性震级（BMC）评估方法，对2010年至2023年期间中国地震台网记录的首都圈地区的地震目录进行分析，通过迭代优化获得计算最小完整性震级MC的最优扫描半径和先验的MC模型，根据数据误差的高斯分布特征推导出MC的先验和似然的概率分布，最后得到MC的后验估计。BMC将地震台站分布的先验信息与局部观测值相结合，权重由各自的不确定性决定，给出了少震区域的MC估计值，并且降低了结果的不确定性。评估结果显示，2010年至今首都圈地区地震监测能力较强，但MC的空间分布不均匀，监测能力最强的地区其MC可达到0左右，较弱地区仅能达到2.7左右。采用最大曲率法评估了研究区域1966年至2023年整体的完整性震级的变化，结果显示首都圈的地震监测能力逐渐提升，2010年后提升较显著。此外，本文还对比了最大曲率法（MAXC）、拟合优度法（GFT）和中位数分段斜率方法（MBASS）在研究区域的结果，认为方法的选择和计算参数对评估结果有不同程度的影响。
Abstract:Earthquake catalog completeness, crucial for seismicity analyses, is defined by the completeness magnitude （MC）: The lowest magnitude at which all earthquakes are reliably detected. Accurate MC estimation is essential for seismic hazard assessments and earthquake studies. Traditional methods often solely rely on the frequency-magnitude distribution （FMD） and Gutenberg-Richter （G-R） law, neglecting station coverage, making them unsuitable for regions with low seismicity and susceptible to the subjective selection of calculation parameters. This study utilizes the Bayesian magnitude of completeness （BMC） method to analyze the Capital Circle Region of China’s （37°N—42°N, 114°E—120°E） earthquake catalog from 2010 to 2023, a period marked by significant network upgrades. Initially, we assessed the overall catalog completeness from 1966 to 2023 using the maximum curvature （MAXC） method, and the results revealed improved monitoring capabilities, especially after 2010, with MC consistently between 0.5 and 1.5. Focusing on 2010—2023 （
23546 events, 153 stations）, we employed the BMC method with a two-step process: ① Optimizing spatial resolution and prior model parameters based on the relationship between MC and station density （distance to the k-th nearest station）; ② Integrating prior information with observed MC values using Bayesian inference. Iterative optimization yielded the optimal scanning radius and prior MC model, from which assuming Gaussian-distributed errors, prior and likelihood distributions were derived, leading to a posterior MC estimate; the BMC method integrates station distribution priors with local MC observations, weighted by their uncertainties, enabling MC estimation in low-seismicity regions and reducing overall uncertainty. The optimized scanning radius R varied spatially, smaller in densely instrumented areas like Beijing. The prior model of Capital Circle Region differed significantly from that of Taiwan, highlighting the need for region-specific models. Posterior MC estimates from BMC showed reduced uncertainty compared to observed MC from MAXC, demonstrating the value of integrating station data. Results revealed spatially heterogeneous monitoring capabilities, with MC reaching 2.7 in regions with relatively weak monitoring capabilities. We compared BMC with three FMD-based methods $ ［ $MAXC, goodness-of-fit test （GFT）, and median-based analysis of segment slope （MBASS）$ ］ $ at varying radii （5—75 km）. GFT yielded the most conservative estimates $ ［ $MC（GFT）＞MC（MBASS）＞MC（MAXC）$ ］ $. Larger radii smoothed spatial variations and potentially overestimated MC in well-monitored areas, emphasizing the importance of careful radius selection. BMC, unlike probability-based magnitude of completeness （PMC）, optimizes R and incorporates station, but does not account for variations among individual stations. PMC, while not reliant on the G-R model, has limitations due to its preset starting magnitude. Our findings show a significant improvement in the seismic network’s monitoring capabilities after 2010 compared to previous levels. The spatial MC variability highlights the importance of localized assessments for hazard analysis. This study demonstrates BMC’s efficacy for robust MC estimation, crucial for accurate seismic hazard characterization and mitigation strategies.-
Keywords:
- completeness magnitude /
- Capital Circle Region /
- BMC method /
- Bayesian estimation
-
引言
地震目录的完整性几乎对于所有地震活动分析都至关重要（Habermann，1987），完整性震级（completeness magnitude，缩写为MC）被定义为在一个特定的空间-时间范围内，可以可靠地检测到所有地震事件的最低震级（Rydelek，Sacks，1989；Woessner，Wiemer，2005）。研究地震参数的分布、地震预测以及概率地震危险性评估等通常需要考虑到地震目录的完整性。如果一组地震事件的震级分布符合古登堡—里克特（Gutenberg-Richter）模型（Gutenberg，Richter，1944），也称为G-R模型，基于对震级频率分布（frequency-magnitude distribution，缩写为FMD）的分析，MC可以被定义为累积频率震级分布和G-R模型偏离的最小震级（Zúñiga，Wyss，1995）。大多数现有的MC估计方法都基于G-R模型的假设，即假设地震事件具有自相似性（Wiemer，Katsumata，1999；Wiemer，Wyss，2000；Cao，Gao，2002；Woessner，Wiemer，2005；Amorèse，2007；Marsan，Daniel，2007）。Rydelek和Sacks （1989）通过比较白天和夜晚事件在每个震级的增量来估计MC。还有一些方法则是直接基于地震台站的噪声谱来确定地震检测的概率（Gomberg，1991；Kværna，Ringdal，1999）。仅基于地震目录的统计方法受到样本数的限制，在地震活动较少的地区可能无法给出评估结果。一些结合了地震目录和地震台站信息的评估方法可以将评估结果扩展到地震活动较少或无地震活动的地区，从而克服了仅基于地震目录的方法对样本数的依赖。这类方法包括基于概率的完整性震级（probability-based magnitude of completeness，缩写为PMC）方法（Schorlemmer，Woessner，2008；李智超，黄清华，2014；刘芳等，2014；蒋长胜等，2015）和贝叶斯完整性震级（Bayesian magnitude of completeness，缩写为BMC）方法（Mignan et al，2011，2013）。
BMC方法通过对空间分辨率进行优化以最小化MC估计中的空间不均匀性和不确定性。然后基于贝叶斯方法，将关于MC与地震台站密度之间关系的先验信息与本地观测到的MC相结合，并根据不确定性分配权重。这一方法引入了灵活的空间分辨率和贝叶斯估计，并且与b值等地震活动参数的计算具有兼容性，为地震研究人员提供了一个更精细的、更准确的完整性震级估计方法（Mignan et al，2011），可以适用于不同地质和地震活动背景，在BMC方法基础上发展出的基于层次贝叶斯模型的完整性震级评估方法（H-BMC）进一步提高了精度和准确性（Feng et al，2022）。
首都圈作为一个人口密集、经济繁荣的地区，地震活动和地震地质背景对该地区的安全和可持续发展具有重要影响，这个地区也是中国地震台网台站分布较密集的地区，已有的研究显示，随着地震台网的建设升级，研究区域的地震监测能力显著提升，自2006年，各主要城市都能够检测到ML1.6以上的地震事件，而在北京及其周边地区，完整性震级可以达到ML1.0 （李智超，黄清华，2014）。经过“十五”数字化改造，地震台网升级并稳定运行后，首都圈地区已经积累了超过十年的观测资料。为了客观评估这段时间地震目录的完整性，本文采用贝叶斯完整性震级方法对该地区的地震数据进行分析。
1. 数据
本文研究的区域为首都圈地区（37°—42°N，114°—120°E），包括北京及其周边地区、河北、天津等地，地处张家口—渤海地震带、华北平原地震带和山西地震带的交会部位，新构造活动频繁，历史上曾发生过多次破坏性地震（周依等，2022）。本文使用了中国地震台网产出的1966年至2023年的5万3732条地震目录，震级单位为ML（如无特殊说明，下文震级单位默认为ML）。为了考察研究区域地震目录完整性在时间上的变化，采用最大曲率法（maximum curvature，缩写为MAXC）（Wyss et al，1999；Wiemer，Wyss，2000）以500个地震事件为窗长评估研究区域整体的MC。时序特征（图1）显示，研究区域的MC在1966年至2023年逐渐降低，这与监测能力的提升和地震台网的加密等因素相关。1966年至1976年间，MC维持在2.0左右；自1976年唐山MS7.8地震发生后，研究区域的监测能力有短期的下降，MC仅达到4.0左右；在1977年恢复，1977年至2009年MC在1.0—2.0之间波动。1975年辽宁海城地震后，北京地区子台数量从最初的8个增加到20个，台站的分布延伸到京津唐地区。1998年至2001年期间，首都圈数字地震台网实现实时传输，台站数量增加到107个。1996年开始，中国地震局进行了“中国数字地震监测系统”的建设，地震监测能力进一步提升（刘瑞丰等，2003）。2003年中国地震局启动了“中国数字地震观测网络”项目，到2007年完成（刘瑞丰等，2008），2010以后观测网络稳定运行，监测能力有明显提升，如图1所示，2010年后研究区域的MC保持在0.5—1.5之间，有较大程度的下降。
![]()
图 1 首都圈地区域（37°N—42°N，114°E—120°E） MC随时间的变化δMC为200次重采样给出的误差估计Figure 1. Temporal variation of MC in the Capital Circle Region of China （37°N—42°N，114°E—120°E）δMC is the bootstrap error estimate （200 resamples）本文着重研究2010年至2023年的中国地震台网完成升级并均稳定运行的阶段首都圈地区MC的空间分布，研究时段内地震目录包含了2万
3546 个地震事件，有效使用的台站共153个。图2a给出了地震事件和台站的分布，可以看出，地震事件的空间分布不均匀，主要集中在唐山地区和张渤地震带，而台站的分布也存在不均匀性，北京和天津地区的台站密度相对较高。图2b给出了地震观测报告中的震级-频次分布情况，地震目录的震级范围为−0.9—5.3。图2c给出了研究区域地震的深度-频次分布，地震深度主要分布在30 km以内，优势深度区间为5—15 km。图2d给出了研究区域地震的月频次时序图，2010年以后研究区域的地震月频次相对平稳。![]()
2 2010年至2023年7月首都圈地区地震目录（a） 地震和台站分布图2. Earthquake catalog data for the Capital Circle Region from January 2010 to July 2023（a） Distribution of earthquakes and stations![]()
2 2010年至2023年7月首都圈地区地震目录（b） 震级-频次分布/震级累积分布函数；（c） 深度-频次分布图；（d） 时间-频次分布图2. Earthquake catalog data for the Capital Circle Region from January 2010 to July 2023（b） Magnitude-frequency distribution/cumulative magnitude distribution frequency；（c） Depth-frequency distribution；（d） Time-frequency distribution2. BMC方法原理
BMC （Bayesian magnitude of completeness）方法包括先验模型拟合和后验结果求解两个步骤：第一个步骤是同时优化先验模型参数和扫描半径R；第二个步骤是根据先验模型和观测结果，通过贝叶斯方法获得MC的后验估计。
2.1 空间分辨率和先验模型优化
在第一个步骤中，BMC方法通过优化空间分辨率减小MC估计中的空间不均匀性和不确定性，并获得一个先验模型。先验模型反映了MC与台站密度的经验关系，其中台站的密度使用研究位置与按照由近及远的顺序排列的第k个台站的距离$ {d}_{k} $来表征，模型包含3个参数$ ［ $式（1）$ ］ $。优化过程通过迭代实现：首先使用网格扫描方法计算不同位置的MC值，然后根据式（1）拟合$ {d}_{k} $与MC的经验关系。一般情况下，4个以上台站能够较好地约束震中位置和震源深度，本文根据中国地震台网的实际情况选择k＝4。之后根据式（2）优化空间扫描半径R，采用优化后的R重新计算MC，并根据式（1）再次修正先验模型。经过迭代优化先验模型参数c1，c2，c3和扫描半径R，使得先验模型稳定。参数c1，c2和c3是经验性的。
$$ {M}_{{\mathrm{C}}} （ {d}_{k} ） {\text{＝}} {c}_{1}{d}_{k}^{{c}_{2}} {\text{＋}}{c}_{3} \text{，} $$ (1) $$ R {\text{＝}} \frac{1}{2}\left[\left(\frac{{c}_{1}{d}_{k}^{{c}_{2}} {\text{＋}} {\sigma }_{{\mathrm{prior}}}}{{c}_{1}}\right)^{\tfrac{1}{{c}_{2}}} {\text{－}}\left(\frac{{c}_{1}{d}_{k}^{{c}_{2}}-{\sigma }_{{\mathrm{prior}}}}{{c}_{1}}\right)^{\tfrac{1}{{c}_{2}}}\right] \text{，} $$ (2) 式中，σ为先验MC值与观测到的MC值之间残差的标准偏差。
2.2 MC后验估计
在BMC方法计算的第二个步骤中，采用贝叶斯方法整合两方面信息：MC与台站密度关系的先验信息，以及局部MC观测值。考虑到两者各自的不确定性，进行加权计算。MC与地震台站距离dk之间的观测经验关系被视为先验信息，用先验概率分布P（MC）描述$ ［ $式（3）$ ］ $，考虑了模型的不确定性σ。先验模型可以不依赖观测数据，基于台站分布来获得MC的预测模型。由实际观测资料估计得到的MC可以通过似然概率分布P（${M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}} $|MC）表示$ ［ $式（4）$ ］ $，这个似然分布考虑了观测数据的不确定性σ0，表示对于给定的MC值观测到${M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}} $的可能性。在贝叶斯框架中，根据观测数据和似然函数，可以使用贝叶斯定理计算MC的后验分布$ ［ $式（5）$ ］ $，后验分布P（MC|${M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}} $）表示在考虑观测数据后对MC的认识的更新，通过最大后验估计方法给出MC的估计值$ ［ $式（6）$ ］ $。
根据先验模型的拟合情况，预测的MC值（${M}_{\mathrm{C}}^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}} $）与观测的MC值（${M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}}$）之间的残差近似为一个高斯分布，因此将先验概率P （${M}_{\mathrm{C}}$）分布定义为
$$ P （ {M}_{{\mathrm{C}}} ） {\text{＝}} \frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{ {\text{－}}{ （ {M}_{{\mathrm{C}}} {\text{－}}{M}_{{\mathrm{C}}}^{{\mathrm{pred}}} ） }^{2}}{2{\sigma }^{2}}\right) \text{，} $$ (3) 同时将似然概率分布表示为
$$ P （ {M}_{{\mathrm{C}}}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mid {M}_{{\mathrm{C}}} ） {\text{＝}} \frac{1}{\sqrt{2\pi }{\sigma }_{0}}{\mathrm{exp}}\left( \frac{ {\text{－}}{ （ {M}_{{\mathrm{C}}}^{{\mathrm{obs}}} {\text{－}}{M}_{{\mathrm{C}}} ） }^{2}}{2{\sigma }_{0}^{2}}\right) \text{，} $$ (4) 式中，${M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}} $是基于最优局部半径（通过空间分辨率优化过程获得）的局部观测值，而σ0是通过自助法（bootstrap）（Woessner，Wiemer，2005；Amorèse，2007）计算得出的局部标准误差。
根据贝叶斯原理，后验概率为
$$ P （ {M}_{\mathrm{C}}\mid {M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}} ） {\text{＝}}\frac{P （ {M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mid {M}_{\mathrm{C}} ） P （ {M}_{\mathrm{C}} ） }{P （ {M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}} ） }\text{，} $$ (5) 将式（3）和式（4）带入式（5）后，后验$ {M}_{\mathrm{C}}^{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}} $可表示为两个高斯分布的乘积，根据最大后验估计可以导出后验的MC估计（$ {M}_{\mathrm{C}}^{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}} $）为
$$ {M}_{\mathrm{C}}^{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}{\text{＝}}\frac{{M}_{\mathrm{C}}^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}{\sigma }_{0}^{2}{\text{＋}}{M}_{\mathrm{C}}^{\mathrm{o}\mathrm{b}\mathrm{s}}{\sigma }^{2}}{{\sigma }^{2}{\text{＋}}{\sigma }_{0}^{2}} \text{，} $$ (6) $ {M}_{\mathrm{C}}^{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}} $可以理解为预测值和观测值根据各自的不确定性的加权平均。后验标准差σpost可以表示为
$$ {\sigma }_{\text{post}}{\text{＝}}\sqrt{\frac{{\sigma }^{2}{\sigma }_{0}^{2}}{{\sigma }^{2}{\text{＋}}{\sigma }_{0}^{2}}} \text{，} $$ (7) 式中，σ0表示观测数据的标准差，σ表示先验模型拟合的标准差。
3. 基于BMC的完整性震级分析
3.1 分辨率优化和先验模型拟合
对研究区域按0.1°×0.1°的间隔划分网格，使用MAXC方法计算每个格点的MC值，使用蒙特卡洛近似自助法（bootstrap）对用于计算的每个格点的原始样本进行有放回的采样，获得用于计算MC的等量样本集，参考Woessner和Wiemer （2005）重复200次，将自助法得到的结果的算术平均值作为MC的估计，将标准差σ0作为MC误差的估计。然后根据式（1）拟合获得初步的先验模型参数，模型误差由先验MC值与观测到的MC值之间残差的标准偏差σ来表示。之后根据式（2）优化MC的扫描半径，再使用经过优化后的扫描半径重新计算每个格点的MC值以修正式（1）中的模型参数。进行迭代直到模型参数稳定，从而获得最优扫描半径R和BMC先验模型。
图3a给出了优化后的扫描半径R的空间分布图，在台站密度高的区域，最优的R较小，如北京的部分地区R在20 km以下，R反映的是MC计算的最优分辨率。图3b中红色曲线给出了优化后的BMC先验模型，横坐标为研究点到与之最近的第4个台站的距离$ {d}_{k} $，纵坐标为MC观测值，蓝色圆点给出了优化后MC观测值与$ {d}_{k} $的对应关系。图3b右下方的子图给出了拟合残差的概率密度分布图，符合高斯分布（子图中蓝色曲线为拟合的高斯曲线），标准差σ为0.44，作为先验模型的误差估计。先验模型参数c1为5.96，c2为0.08，c3为−7。图3b还展示了台湾地区BMC先验模型（黑色曲线，Mignan et al，2011），首都圈地区的实际观测数据拟合得到的先验模型与台湾地区的模型存在较大的差异，这种差异可能反映了震级标度的系统性差异，在相同的台站密度下，台湾地区模型的MC高于首都圈地区的模型。
![]()
图 3 优化后的半径和BMC先验模型（a） 扫描半径R的空间分布图；（b） BMC先验模型，右下角为拟合残差概率密度分布图Figure 3. Optimized radius and BMC prior model（a） Spatial distribution map of scanning radius R；（b） BMC prior model，the bottom right panel shows the probability density distribution of the fitting residuals3.2 后验结果计算
采用优化后半径和MAXC方法得到了观测的MC分布（图4a）、先验的MC分布（图4c）以及后验的MC分布（图4d），并给出评估结果和误差估计分布图，图4b展示了观测的MC分布的不确定性估计，图4e提供了后验BMC结果的不确定性估计。在MAXC方法中采用了200次自助法来估计结果的误差范围，符合计算条件格点数为
1816 （总格点数3331 ，可计算格点占比58%），图4b显示大部分误差分布在0.7以下，个别地区的误差较高。根据先验模型生成了研究区域完整的MC分布（图4c），由拟合残差的标准差估计先验模型的误差。根据贝叶斯定理$ ［ $式（6）$ ］ $将先验信息与观测值进行了融合，生成了预测和观测完整性震级的加权平均值，作为MC的后验估计$ ［ $式（6），图4d$ ］ $，其中权重与它们各自的不确定性成正比（Mignan et al，2011），同时估计了其误差范围$ ［ $式（7），图4e$ ］ $。图5b给出了观测值MC （MAXC）、先验模型MC （BMC先验）和后验结果MC （BMC后验）的公共可计算格点的MC累积分布函数（cumulative distribution function，缩写为CDF），MC （BMC后验）的CDF介于MC （BMC先验）和观测值MC （MAXC）之间。误差的CDF比较（图5a）显示，由于后验模型融合了台站分布和实际观测数据，误差水平有所降低。图4f展示了后验模型与先验模型的差值分布，可以看出研究区域北侧为负值，研究区域南侧为正值，这表明地震较丰富且台站较密集的地区后验结果的MC值低于先验模型的预测值，而地震较少且台站密度较低的地区后验结果的MC值高于先验模型的预测值，先验模型根据台站密度的分布和观测值的经验关系给出研究区域每个格点的预测值，对于没有观测值的区域，先验模型给出合理的估计值，对于已有观测值的区域，BMC方法获得的后验结果根据贝叶斯原理融合了观测结果和台站密度信息。在先验模型中，台站密度与MC的分布关系仍然存在一定的离散性（图3b），这表明除了密度因素外，还有其它主要的未知因素影响MC，因此MC （BMC先验）与MC （BMC后验）有一定差别。![]()
图 4 BMC方法结果分布图（a） MC （MAXC）；（b） δMC （MAXC）；（c） MC （BMC先验）；（d） MC （BMC后验）；（e） δMC （BMC后验）；（f） MC （BMC后验）−MC （BMC先验）Figure 4. Distribution map of BMC method results（a） MC （MAXC）；（b） δMC （MAXC）；（c） MC （BMC prior）；（d） MC （BMC posterior）；（e） δMC （BMC posterior）；（f） MC （BMC posterior）−MC （BMC prior）![]()
图 5 误差（a）和BMC结果（b）的累积分布CDF对比Figure 5. Comparison of CDF of error （a） and BMC results （b）4. 不同方法和计算参数的比较
为了分析方法和参数的选择对MC评估结果的影响，选择三种方法和不同的计算参数对研究区域的MC进行评估，并对比结果。
4.1 方法介绍
在绘制MC的空间分布图时，有三种常用的数据选取方法，每种方法都有其优点和限制。第一种方法是网格划分方法，这种方法将研究区域划分成网格，使用每个网格内的地震目录数据来计算MC。统计方法评估MC需要地震目录中包含一定数量的地震，划分网格后，一些网格内的地震数量可能不足以满足统计需求。第二种方法是固定半径法，这种方法是通过对每个格点采用一定半径范围内的地震目录数据来计算MC，需要考虑扫描半径R和最低样本数N。第三种方法是固定样本数量法，即固定样本数量N，选择最近的N个事件来计算MC。这种方法的优势在于MC的不确定性被约束，缺点是高地震密度梯度区域的结果可能不真实或变化较大。本文采用固定半径方法对比3种方法的评估结果。MC的不确定性通常随着样本数量的增加而减小，因此大多数评估方法在处理更多事件时会具有更低的不确定性。选择合适的数据选取方法、扫描半径R和样本数阈值N对于准确估计地震目录的完整性震级MC至关重要。这些参数的选择应该基于具体研究区域、地震活动密度、地质背景和可用数据等因素进行权衡和调整，以确保MC估计的可靠性。
我们采用了三种不同的基于FMD的MC估计方法，包括MAXC方法，拟合优度方法（goodness-of-fit test，缩写为GFT）（Wiemer，Wyss，2000）和中位数分段斜率方法（median-based analysis of the segment slope，缩写为MBASS）（Amorèse，2007），并分别设置了不同的扫描半径R来评估研究区域的MC值分布。MAXC方法通过在FMD图中找到具有最高频率的震级区间来估算MC。GFT方法是根据实际和理论震级-频度分布下的拟合度来搜索得出完整性震级，将观察到的FMD与合成的G-R模型分布进行比较，并将观测FMD与合成FMD之间残差低于某个阈值（本文选取5%，即95%拟合度）的最小截止震级作为MC的估计。MBASS方法通过迭代方法在非累积FMD中寻找斜率变化，通过Wilcoxon秩和检验来接受或拒绝零假设（斜率无变化），检测出主要不连续性（概率最大的斜率变化处）作为MC的估计。
图6展示了三种方法在点（40.5°N，116.5°E）的评估结果（采用半径20 km范围内的地震目录），图6a给出MAXC方法的结果，根据MAXC方法，在非累积频次分布的最高值对应着MC为0.0。图6b给出了GFT方法的评估结果，图中展示了观测FMD与合成的G-R分布的残差随着截止震级的变化曲线，由残差低于5%的最小截止震级得到MC估计值为0.2。图6c给出了MBASS方法的评估结果，图中展示了
1000 个自助法样本中检测到主要不连续性的分布，频次最高值对应的MC估计值为0.0。结果显示，MAXC和MBASS结果相同，GFT结果较保守。![]()
图 6 MAXC （a），GFT （b）和MBASS （c）三种方法在点（40.5°N，116.5°E）的评估结果Figure 6. The results of MAXC （a），GFT （b） and MBASS （c） at point （40.5°N，116.5°E）4.2 计算结果和对比分析
图7给出由MAXC，GFT和MBASS在不同扫描半径（5，10，20，50，75 km）下的MC分布图，N设为50。结果表明：随着扫描半径R的增加，每种方法可以计算的格点数量也增加，这是由于较大的扫描半径可以提供更多的地震事件样本，从而有更多格点的样本数达到N，同时降低MC估计的不确定性。在较小的扫描半径（5 km和10 km）下，三种方法仅能计算少量的格点，其中MAXC和GFT方法的可计算格点比MBASS多。这表明MBASS方法对样本数的要求更加严格。当扫描半径R在20 km以下时，三种方法的结果接近，当扫描半径R在50 km以上时，不同方法出现了差异。
![]()
图 7 MAXC （a），GFT （b）和MBASS （c）三种方法在不同扫描半径R下的MC分布图Figure 7. MC map of MAXC （a），GFT （b） and MBASS （c） at different scanning radii R比较MAXC方法在不同的扫描半径的公共可计算格点的结果，经验累积分布函数CDF对比显示（图8a），随着R值的增大，MC的值更加集中，这表明较大的R值可以抑制空间差异，使结果更加均匀，此外扫描半径增大时，MC评估结果整体上逐渐降低。三种方法结果的CDF对比（图8b，8c）显示，MC （GFT）＞MC （MBASS）＞MC （MAXC），MAXC和MBASS结果接近，GFT结果较保守（MC偏高），与图5的结论一致。
![]()
图 8 不同方法和扫描半径R下MC值的CDF对比（a） MAXC方法在三种扫描半径下的结果；（b）−（c） MAXC，GFT和MBASS方法在R＝50 km和R＝25 km时的结果Figure 8. Comparison of CDF of MC at different methods and R（a） The MAXC results for different radii；（b）−（c） Results of three methods with R＝50 km and R＝25 kmBMC方法和三种基于FMD的方法给出的MC评估结果均显示，北京地区的监测能力较强，监测能力最强的区域达到0左右。为了直观地检验计算结果和考察扫描半径对结果的影响，采用MAXC方法，以10，20和50 km三种扫描半径计算点（40.5°N，116.5°E）的MC，并给出相应的FMD分布图（图9），在三种半径下累积频次分布均符合G-R模型。当R＝10 km时，MC为−0.1，FMD的顶部较尖锐；当R＝20 km时，MC为0.0，FMD的顶部比R＝10 km时平缓；在R＝30 km时，MC为0.3，FMD的顶部更加平滑，呈现圆弧状。随着R的增加，MC升高，FMD的分布顶端趋于平缓。Mignan （2012）采用模拟数据评估了FMD分布形态对结果的影响，将FMD分布分为顶部尖锐的（Angular）和顶部平缓的（Gradually Curved），顶部平缓的FMD模型会导致MC高估，图9支持了这个结论。究其原因，较大的扫描半径可能包含了不同性质的数据集，会使频率-震级分布（FMD）趋向平缓，而观测的累积频次分布就会在较高的值处偏离G-R模型，从而导致MC被高估。这种高估的情况仅出现在监测能力较好的区域。图8a显示，当R较大时，MAXC方法的结果更加集中，而整体上是降低的，较大的R模糊了空间细节特征。因此，在选择扫描参数时，需要仔细权衡这些因素，以确保获得可靠的MC估计结果。
![]()
图 9 扫描半径R为10 km （a），20 km （b）和50 km （c）时点（40.5°N，116.5°E）的MAXC结果Figure 9. The MAXC results with scanning radius R＝10 km （a），20 km （b） and 50 km （c） at point （40.5°N，116.5°E）5. 讨论与结论
本文采用MAXC方法分析了1966年以来首都圈地区完整性震级随时间的变化，采用BMC方法和三种基于FMD的方法着重分析了2010年以后的地震目录完整性。为了考察参数对评估结果的影响，我们对比了三种基于FMD的方法，即MAXC，GFT和MBASS在不同参数下的评估结果。随着扫描半径的增大，不同方法差异增大，局部差异被抑制，且监测能力较好的地区的MC可能被高估，这表明扫描半径的选取会影响结果的分辨率和可靠性。我们采用BMC方法对2010年以后的资料进行分析，获得了先验、观测和后验的评估结果，通过迭代优化确定扫描半径先验模型，避免了固定扫描半径的主观性。我们发现首都圈的先验模型与台湾地区有系统性的差异，这可能是震级标度差异造成的。通过BMC方法，我们利用台站密度和观测值拟合得到先验模型，将评估结果延拓到少震地区，得到研究区域完整的评估结果，克服了样本数量的限制，同时基于贝叶斯原理通过结合先验知识和观测数据降低不确定性。
和常用的基于FMD的评估方法比较，BMC方法有三方面优点：优化了扫描半径R以提高MC的分辨率；给出了样本不足区域的MC的合理估计；降低了结果的不确定性。与PMC方法比较，BMC方法和PMC方法都可以给出少震地区的MC评估，可以用于研究和规划现有地震台网的改进，但是这两种方法基于不同的假设：PMC方法假设台网能够记录到更完整的目录，在这个假设下估计单台的检测概率时考虑了不同台站的差异性，不需依赖G-R模型的假设即可给出完整震级的估计，但预设的起算震级限制了可评估的震级下限，台网目录完整的假设也可能不符合实际情况；BMC方法则是基于台网密度和监测能力的经验关系，其先验模型比较简单，没有考虑不同台站的差异性，也没有考虑可能影响监测能力的其它因素，其预测结果与实际评估结果可能存在一定偏差。
综合计算结果显示，首都圈地区及其周边地区的地震检测能力较高，北京市境内部分区域的MC值约为ML0左右，整个研究区域的MC值低于2.7。武敏捷等（2013）对截至2011年以前首都圈地区地震监测能力的分区和分时段研究结果显示，2000以来MC最低达到1.2，高于本文评估结果，本文采用的数据时段为2010年至2023年，结果的差异反映了2010年以来监测能力的提高。李智超和黄清华（2014）采用PMC方法估计首都圈的地震监测能力，数据截至2009年，结果显示各主要城市都达到1.6，在北京及其周边地区可以达到1.0以内，本文评估结果在部分区域低于上述结果，除了反映出2010年以后监测能力的提升，也体现了方法的差异，PMC方法的计算过程预设了起算震级，可能限制了MC的评估下限。
采用MAXC方法分析了研究区域整体的MC随时间的变化，我们发现首都圈地区的监测能力逐步提升，2010年后MC值保持在0.5至1.5之间。
我们对比了MAXC，GFT和MBASS三种方法在五种扫描半径下的MC评估结果，在相同的扫描版半径下，MC （GFT）＞MC （MBASS）＞MC （MAXC），MAXC和MBASS结果接近，GFT结果较保守。当扫描半径增大时，MAXC方法在监测能力较好的区域会高估MC值，MC值趋于集中，整体上是降低的，这表明较大的R模糊了空间细节特征。
采用首都圈地区的数据，我们同时优化了BMC的先验模型和最佳扫描半径，首都圈地区的先验模型与台湾地区的模型存在较大差异。
综合先验模型和观测结果，我们得到了首都圈地区的BMC后验评估结果。与观测结果相比，后验结果的不确定性有所降低。计算结果显示首都圈地区及其周边地区的地震检测能力较高，北京境内的部分区域的MC值达到0左右，研究区域整体的MC值低于2.7。
本研究采用了zmap，rseismNet和GRTo软件包进行计算和分析，两位评审专家提供的建议对论文质量的提升起到了至关重要的作用，作者在此一并表示感谢。
-
![]()
图 1 首都圈地区域（37°N—42°N，114°E—120°E） MC随时间的变化
δMC为200次重采样给出的误差估计
Figure 1. Temporal variation of MC in the Capital Circle Region of China （37°N—42°N，114°E—120°E）
δMC is the bootstrap error estimate （200 resamples）
![]()
2 2010年至2023年7月首都圈地区地震目录
（a） 地震和台站分布图
2. Earthquake catalog data for the Capital Circle Region from January 2010 to July 2023
（a） Distribution of earthquakes and stations
![]()
2 2010年至2023年7月首都圈地区地震目录
（b） 震级-频次分布/震级累积分布函数；（c） 深度-频次分布图；（d） 时间-频次分布图
2. Earthquake catalog data for the Capital Circle Region from January 2010 to July 2023
（b） Magnitude-frequency distribution/cumulative magnitude distribution frequency；（c） Depth-frequency distribution；（d） Time-frequency distribution
![]()
图 3 优化后的半径和BMC先验模型
（a） 扫描半径R的空间分布图；（b） BMC先验模型，右下角为拟合残差概率密度分布图
Figure 3. Optimized radius and BMC prior model
（a） Spatial distribution map of scanning radius R；（b） BMC prior model，the bottom right panel shows the probability density distribution of the fitting residuals
![]()
图 4 BMC方法结果分布图
（a） MC （MAXC）；（b） δMC （MAXC）；（c） MC （BMC先验）；（d） MC （BMC后验）；（e） δMC （BMC后验）；（f） MC （BMC后验）−MC （BMC先验）
Figure 4. Distribution map of BMC method results
（a） MC （MAXC）；（b） δMC （MAXC）；（c） MC （BMC prior）；（d） MC （BMC posterior）；（e） δMC （BMC posterior）；（f） MC （BMC posterior）−MC （BMC prior）
![]()
图 5 误差（a）和BMC结果（b）的累积分布CDF对比
Figure 5. Comparison of CDF of error （a） and BMC results （b）
![]()
图 6 MAXC （a），GFT （b）和MBASS （c）三种方法在点（40.5°N，116.5°E）的评估结果
Figure 6. The results of MAXC （a），GFT （b） and MBASS （c） at point （40.5°N，116.5°E）
![]()
图 7 MAXC （a），GFT （b）和MBASS （c）三种方法在不同扫描半径R下的MC分布图
Figure 7. MC map of MAXC （a），GFT （b） and MBASS （c） at different scanning radii R
![]()
图 8 不同方法和扫描半径R下MC值的CDF对比
（a） MAXC方法在三种扫描半径下的结果；（b）−（c） MAXC，GFT和MBASS方法在R＝50 km和R＝25 km时的结果
Figure 8. Comparison of CDF of MC at different methods and R
（a） The MAXC results for different radii；（b）−（c） Results of three methods with R＝50 km and R＝25 km
-
蒋长胜,房立华,韩立波,王未来,郭路杰. 2015. 利用PMC方法评估地震台阵的地震检测能力:以西昌流动地震台阵为例[J]. 地球物理学报,58(3):832–843. doi: 10.6038/cjg20150313 Jiang C S,Fang L H,Han L B,Wang W L,Guo L J. 2015. Assessment of earthquake detection capability for the seismic array:A case study of the Xichang seismic array[J]. Chinese Journal of Geophysics,58(3):832–843 (in Chinese).
李智超,黄清华. 2014. 基于概率完备震级评估首都圈地震台网检测能力[J]. 地球物理学报,57(8):2584–2593. doi: 10.6038/cjg20140818 Li Z C,Huang Q H. 2014. Assessment of detectability of the Capital-Circle Seismic Network by using the probability-based magnitude of completeness (PMC) method[J]. Chinese Journal of Geophysics,57(8):2584–2593 (in Chinese).
刘芳,蒋长胜,张帆,杨彦明,梁莹,王磊,苗春兰. 2014. 内蒙古区域地震台网监测能力研究[J]. 地震学报,36(5):919–929. Liu F,Jiang C S,Zhang F,Yang Y M,Liang Y,Wang L,Miao C L. 2014. A study on detection capability of the Inner Mongolia regional seismic network[J]. Acta Seismologica Sinica,36(5):919–929 (in Chinese).
刘瑞丰,吴忠良,阴朝民,陈运泰,庄灿涛. 2003. 中国地震台网数字化改造的进展[J]. 地震学报,25(5):535–540. doi: 10.3321/j.issn:0253-3782.2003.05.009 Liu R F,Wu Z L,Yin C M,Chen Y T,Zhuang C T. 2003. Development of China digital seismological observational systems[J]. Acta Seismologica Sinica,25(5):535–540 (in Chinese).
刘瑞丰,高景春,陈运泰,吴忠良,黄志斌,徐志国,孙丽. 2008. 中国数字地震台网的建设与发展[J]. 地震学报,30(5):533–539. doi: 10.3321/j.issn:0253-3782.2008.05.012 Liu R F,Gao J C,Chen Y T,Wu Z L,Huang Z B,Xu Z G,Sun L. 2008. Construction and development of digital seismograph networks in China[J]. Acta Seismologica Sinica,30(5):533–539 (in Chinese).
武敏捷,武安绪,龙锋,林向东,赵雯佳. 2013. 首都圈地震活动完整性震级时空分布特征研究[J]. 地震,33(3):98–106. doi: 10.3969/j.issn.1000-3274.2013.03.011 Wu M J,Wu A X,Long F,Lin X D,Zhao W J. 2013. Temporal-spatial distribution characteristics of minimum magnitude of completeness in the capital region of China[J]. Earthquake,33(3):98–106 (in Chinese).
周依,郭蕾,王亚玲. 2022. 首都圈地区中强地震前固体潮调制比特征分析[J]. 大地测量与地球动力学,42(11):1166–1170. Zhou Y,Guo L,Wang Y L. 2022. The characteristics of earth tide modulation ratio before moderate-strong earthquakes in the capital circle area of China[J]. Journal of Geodesy and Geodynamics,42(11):1166–1170 (in Chinese).
Amorèse D. 2007. Applying a change-point detection method on frequency-magnitude distributions[J]. Bull Seismol Soc Am,97(5):1742–1749. doi: 10.1785/0120060181
Cao A M,Gao S S. 2002. Temporal variation of seismic b‐values beneath northeastern Japan island arc[J]. Geophys Res Lett,29(9):1334.
Feng Y,Mignan A,Sornette D,Li J W. 2022. Hierarchical Bayesian modeling for improved high‐resolution mapping of the completeness magnitude of earthquake catalogs[J]. Seismol Soc Am,93(4):2126–2137.
Gomberg J. 1991. Seismicity and detection/location threshold in the southern Great Basin seismic network[J]. J Geophys Res:Solid Earth,96(B10):16401–16414. doi: 10.1029/91JB01593
Gutenberg B,Richter C F. 1944. Frequency of earthquakes in California[J]. Bull Seismol Soc Am,34(4):185–188. doi: 10.1785/BSSA0340040185
Habermann R E. 1987. Man-made changes of seismicity rates[J]. Bull Seismol Soc Am,77(1):141–159.
Kværna T,Ringdal F. 1999. Seismic threshold monitoring for continuous assessment of global detection capability[J]. Bull Seismol Soc Am,89(4):946–959. doi: 10.1785/BSSA0890040946
Marsan D,Daniel G. 2007. Measuring the heterogeneity of the coseismic stress change following the 1999 MW7.6 Chi‐Chi earthquake[J]. J Geophys Res:Solid Earth,112(B7):B07305.
Mignan A,Werner M J,Wiemer S,Chen C C,Wu Y M. 2011. Bayesian estimation of the spatially varying completeness magnitude of earthquake catalogs[J]. Bull Seismol Soc Am,101(3):1371–1385. doi: 10.1785/0120100223
Mignan A. 2012. Functional shape of the earthquake frequency‐magnitude distribution and completeness magnitude[J]. J Geophys Res:Solid Earth,117(B8):B08302.
Mignan A,Jiang C,Zechar J D,Wiemer S,Wu Z,Huang Z. 2013. Completeness of the Mainland China earthquake catalog and implications for the setup of the China Earthquake Forecast Testing Center[J]. Bull Seismol Soc Am,103(2A):845–859. doi: 10.1785/0120120052
Rydelek P A,Sacks I S. 1989. Testing the completeness of earthquake catalogues and the hypothesis of self-similarity[J]. Nature,337(6204):251–253. doi: 10.1038/337251a0
Schorlemmer D,Woessner J. 2008. Probability of detecting an earthquake[J]. Bull Seismol Soc Am,98(5):2103–2117. doi: 10.1785/0120070105
Wiemer S,Katsumata K. 1999. Spatial variability of seismicity parameters in aftershock zones[J]. J Geophys Res:Solid Earth,104(B6):13135–13151. doi: 10.1029/1999JB900032
Wiemer S,Wyss M. 2000. Minimum magnitude of completeness in earthquake catalogs:Examples from Alaska,the western United States,and Japan[J]. Bull Seismol Soc Am,90(4):859–869. doi: 10.1785/0119990114
Woessner J,Wiemer S. 2005. Assessing the quality of earthquake catalogues:Estimating the magnitude of completeness and its uncertainty[J]. Bull Seismol Soc Am,95(2):684–698. doi: 10.1785/0120040007
Wyss M,Shimazaki K,Ito A. 1999. Introduction:Seismicity Patterns,Their Statistical Significance and Physical Meaning[M]. Basel:Birkhäuser Basel:203−207.
Zúñiga F R,Wyss M. 1995. Inadvertent changes in magnitude reported in earthquake catalogs:Their evaluation through b-value estimates[J]. Bull Seismol Soc Am,85(6):1858–1866. doi: 10.1785/BSSA0850061858
下载:
计量
- 文章访问数: 90
- HTML全文浏览量: 21
- PDF下载量: 33
