0 引言

中国地震预报长期以经验为主，缺乏严格统一的预报效能评估体系，会制约由观测现象、经验认识到理论模型构建和孕震机理解释的进程(张国民等，2002)。基于合理的统计检验模型，对不同时空尺度的地震危险性预测结果做出定量、科学、客观的评价，成为近年来地震预测预报领域的重要研究内容。结合地震发生情况科学系统地评价预测结果，一方面可以验证预测方法的可行性、模型参数配置的合理性，另一方面也是优化预报策略、完善预测模型的客观需求(蒋长胜等，2011；高朝军等，2013；闫伟等，2019)。

本研究系统梳理了常用地震预测模型产出的结果，大致可分为：时间序列类预测结果、线状要素预测结果(如地磁低点位移、重力固体潮零值线等)、边界确定的面状要素预测结果(如年度地震危险区)、边界模糊的面状要素预测结果(如综合概率预测结果等)、包含时间演化特征的面状要素预测结果等。据此国内外学者开展了大量相关研究，提出针对不同预测结果的评价指标，并应用于地震预测效能评价的研究实践。例如，美国南加州地震中心(SCEC)发起的全球“地震可预测性合作研究”计划(Collaboratory for the Study of Earthquake Predictability，CSEP)，应用多种预测模型开展向前预测实践，针对不同模型的预测结果，选用地震数检验N-test(Kagan et al，1995；Schorlemmer et al，2007)、似然检验L-test(Helmstetter et al，2006)、比率检验Ratio-test、Molchan图表法(Molchan，1997；Molchan et al，1992、2008)、区域技能评分ASS-test和接受者操作特性检验ROC(Holliday et al，2005)6种方法进行预测效能的科学评价和模型的对比研究，通过系统总结进而辅助改进预测模型的构建(Schorlemmer et al，2010；Werner et al，2010)。我国自20世纪60年代开展大规模地震预报研究以来，逐步形成了长、中、短、临的预报思路和工作流程(梅世蓉等，1993)，中期预报指对未来一、两年内可能发生破坏性地震的地区和震级的预测，并以年度预报的形式来实施，在我国地震预报实践中应用最广、积累资料最为完善。其主要思路是利用测震和前兆观测资料，结合专家经验和震例研究成果圈定未来一年可能发生破坏性地震的地区并给出震级强度，即年度地震危险区。国内研究人员在年度危险区预测效果的定量检验方面进行了一系列探索，例如，许绍燮(1989)在对1972年提出的震兆鉴定方法总结的基础上，提出了用于评价多边形区域地震预测效能的R值评分方法；其后，石耀霖等(2000)改进了R值计算方法，用虚报面积占全部无震面积代替预报时间窗口总和与总时长比例，并将该方法用于评价1990—1998年间国家地震局年度地震危险区的预报效能，认为中国大陆地区年度危险区预测效果总体上好于随机预报；张国民等(2002)为了讨论地震能否预报和评价当时中国地震预报能力等问题，将中国大陆分为931个1°×1°的网格，以报有震网格数占总网格数比例与虚报网格数占总无震网格数比例为评分值，评价1990—2000年年度危险区预报效果，认为虽然预报水平还相当低且不稳定，但年度危险区预测结果好于随机预测；马宏生等(2004)将中国大陆以107°E为界分为东部和西部，针对1990—2002年的年度地震危险区，计算东西部不同震级下限时的R值评分，结果表明，虽然年度危险区预测水平还很低且不稳定，但总体呈上升趋势。

综上，在众多地震预测结果评价方法中，R值评分方法已被广泛应用于年度危险区预报效果的评价和总结，该方法在优化改进年度危险区预测技术方法上发挥了重要的作用。由预测意见到结果评分，是一种“后验”的统计评价模型，然而对从R值评分方法和地震预测水平、地震活动性统计结果到预测危险区个数、空间尺度参数的“先验”性研究较少，本研究在前人工作的基础上，考虑历史地震的年度发生个数、地震“丛集”分布的统计学特征，利用随机分布和R值评分反推模型，从统计学的角度研究了年度危险区个数的最佳估值，对地震预测预报实践具有一定的参考意义。

1 年度地震危险区最大空间占有率计算公式推导 1.1 R值评分原理

R值评分方法是当前普遍认可和使用的地震预测效能评价方法(许绍燮，1989；张国民等，2002；闫伟等，2019)，20世纪80年代以来不断完善并广泛应用于年度地震危险区的预测效能检验，本文采用许绍燮(1989)的计算方法，用地震报准率和预测的时间或空间占有率之差计算R值评分，具体公式如下

$R=h-\tau=\frac{N_p}{N_g}-\frac{S_p}{S_g} \geqslant R_0$

(1)

其中，h表示地震报准率，为报准地震的个数N_p与年度发生地震的总数N_g的比值；τ表示预测区域的空间占有率，为年度地震危险区的面积S_p与国土总面积S_g的比值，R值评分即为地震报准率与预报区域的空间占有率之差。最佳的预测效果对应于在最高的报准率(h→1)下付出最小的代价(τ→0)，即R→1；反之，最差的预测效果对应于在最低的报准率(h→0)下付出最大的代价(τ→1)，即R→- 1；通常认为R＞0时，预报的成功率高于随机概率的成功率，该预报结果具有一定的预测效能(张国民等，2002)。

假设地震的发生是相互独立的，一定时段内发生的地震个数符合n重伯努利分布(泊松分布的离散化表达)，该时段内单个地震的发生概率用预测区域的空间占有率τ表示，相应的显著性水平满足

$\alpha=\sum\limits_{i=H}^N\left[\left(\begin{array}{c} N \\ i \end{array}\right) \tau^i(1-\tau)^{N-i}\right], \left(\begin{array}{c} N \\ i \end{array}\right)=\frac{N !}{i !(N-i) !}$

(2)

其中，N表示地震总数，H表示报准地震的个数，α表示置信度。取α=5%，根据τ值反算伯努利分布下的报准率h，进而可以求得对应于置信度α的R值评分，并以R₀表示，当R>R₀时，认为该预测结果具有较高的统计显著性(Zechar et al，2008；Jiang et al，2010)。相应的概率增益Gain的计算方法为

$\text { Gain }=\frac{1-m}{\tau}=\frac{h}{\tau}$

(3)

其中，m表示地震的漏报率。R₀的计算方法为

$R_0=\frac{H}{N}-\tau_H$

(4)

其中，τ_H表示根据报准地震的个数H和5%的显著性水平反推出来的预测区域空间占有率。

1.2 R值评分与Molchan图表法的关系

Molchan图表法(Molchan Error Diagram)是CSEP检验中心常规使用的预测效能检验方法，目前，已被广泛应用于确定性和概率性预测的统计检验和预测效能评估(蒋长胜等，2011；孙丽娜等，2012；张琳琳等，2019)。如图 1(a)所示，该方法主要使用预测区域的时空占有率τ和漏报率m两个变量进行评价统计，对于某一空间展布的预测变量，随着异常判定阈值的降低，预测区域的空间占有率逐渐增大，同时漏报率逐渐降低，二者的变化构成Molchan检验曲线。对于Molchan图表中的某点P₀(τ₀，m₀)，由于其与直线τ+m-1=0(即直线R=0)的距离为$D=\left|\tau_0+m_0-1\right| / \sqrt{2}$，结合R值计算公式可推导出$D=R / \sqrt{2}$，表明Molchan图表中每个点与直线τ+m-1=0的距离D与R值评分成正比，二者可以相互转换(图 1)。以直线τ+m-1=0为分界线，其左下方区域的R值评分大于0，右上方区域的R值评分小于0。二者的不同之处在于：①Molchan图表法适用于边界模糊的概率性预测结果的效能检验；②相比于R值评分方法，Molchan图表法具有更多、更加直观的预测效能对比指标，如式(3)定义的概率增益和显著性曲线，便于对多种预测模型进行预测效能的对比研究。Molchan曲线上各点距离R=0的距离与R值评分的定量转换关系表明，目前国内外常用的预测效能评价方法具有一定的通用性，鉴于年度危险区是边界确定的预测结果，本研究中采用R值评分的方法对其进行效能评价和最佳个数推估。

1.3 基于R值评分的年度地震危险区最大空间占有率计算公式推导

由中国地震台网中心历年的年度危险区预测实效的统计结果可知，目前年度地震危险区报准率的一般水平约为30%(Yu et al，2022)，即h≈0.3；将年度地震总个数N_g、报准率和置信水平α，代入式(1)、(2)，可计算出对应于置信水平α的R₀；令R>R₀，即可求得年度危险区空间占有率的上限，具体的推导公式为

$\tau \leqslant h-R_0$

(5)

$S_p=S_g \cdot \tau=\frac{n {\rm{ \mathsf{ π} }} a b}{4}$

(6)

式(6)中，年度地震危险区通常为椭圆形封闭区域，a、b分别表示其长、短轴，本研究中均取其平均值约为300km；S_g表示国土总面积，约等于9.6×10⁶km²，进而可以反推出年度地震危险区个数的最佳估值n，即

$n \leqslant \frac{4 \cdot S_g \cdot\left(h-R_0\right)}{{\rm{ \mathsf{ π} }} a b}$

(7)

最后的关键问题在于如何科学合理地确定年度地震的发生个数N_g。

2 年度M_S5.0以上地震发生个数 2.1 基于统计分布检验的年度地震发生个数计算方法

年度地震危险区一般用于判定中国大陆东部(107°E以东) M_S5.5以上、西部(107°E以西) M_S6.0以上地震的可能发生区域，因此在年度地震发生个数的第一种计算方法中，以上述起始震级作为目标地震的筛选标准，分别统计了1900—2019年间各年度目标地震的发生个数以及其对应出现的频次，结果详见表 1。进一步计算得到该组数据的均值约为4.5个/年，方差约为8.96，表明该组数据的离散度较高，若利用均值作为年度地震的发生个数，会引入较大的不确定性，并影响年度地震危险区个数评估结果的合理性。

针对上述问题，拟引入统计分布检验方法，结合置信区间确定年度地震发生个数的估值，尽可能减小数据的不确定性，保证漏报率最低；同时考虑到地震目录的完备性，将1900—2019年间的统计结果截断为1950—2019年。依次对该数据源进行正态分布、Gamma分布、泊松分布、指数分布以及rayleigh分布最小二乘拟合，进而利用卡方检验(χ²)筛选出最优分布模型。具体操作方法为：将目标地震的年度发生个数代入各分布模型的拟合概率密度曲线，求出对应的观察概率值X并代入χ²检验计算公式，即可求出各分布模型的拟合概率密度曲线χ²值；对比置信水平为95%、自由度为17的χ²分布临界值，当χ²值小于该临界值时，表明观测值与实际值无显著差异，数据源可以通过该概率分布检验。统计分布检验的计算结果表明，地震年发生个数的最佳拟合分布为Gamma分布，其概率密度函数为

$f(X)=\frac{X^{\alpha-1} \lambda^\alpha e^{-\lambda X}}{\varGamma(\alpha)}$

(8)

其中，α为形状参数，β为尺度参数，λ=1/β。拟合得到的数据分布频次及拟合概率密度曲线如图 2所示。

2.2 改进聚类算法计算年度地震“丛集区”和“孤立发生地震”的个数

由于区域构造应力的集中和地震间相互触发等因素的影响，结合历年来发生地震的空间分布图像，可以得出年度M_S5.0以上地震具有一定的“丛集性”。年度地震危险区一般为具有确定边界的空间区域，丛集性的存在导致多个地震可能同时落入一个地震危险区，当年度地震危险区的范围与地震发生丛集区域的重合率最高时，将取得最佳的预测效能。因此本研究中确定年度地震发生个数的第二种方法为：拟利用适当算法识别多年来年度地震“丛集区”和“孤立发生地震”的个数之和，进而结合多种统计分布模型和χ²检验确定出最佳概率密度曲线，最后代入式(7)计算年度地震危险区个数的最佳估值。

聚类算法是按照某个特定的标准(如距离)将一个数据集分割成不同的类或簇，使得同一个簇内的数据对象的相似性尽可能大，同时不在同一簇中的数据对象的差异性也尽可能的大。本研究拟采用基于距离的聚类算法自动识别地震丛集区，传统聚类算法如K-means等首先需要设置初始中心，进而通过迭代计算各点到初始中心的距离并求解新的中心坐标循环分类，直到达到一定的阈值或迭代次数为止。上述算法严重依赖于初始簇中心的选择，若直接将其应用到丛集地震的识别，难以确定“迭代结束阈值”而影响最终的识别效果；另外，每次计算都是在全数据域中，因此识别的效率较低。鉴于其局限性，本研究提出了一种简洁准确且不需要设置“迭代结束阈值”的地震“丛集区”和“孤立发生地震”的识别算法，具体操作方法如下：

(1) 假设一年中发生了n次地震，首先构建由这些地震间距离组成的初始矩阵D(d₁₁，d₁₂，…，d_1n；…；d_n1，d_n2，…，d_nn)，d_ij表示第i和j个地震之间的距离，且d_ij=d_ji；对该矩阵中的各元素进行判断，若其小于指定的“分离距离” (丛集条件，本研究中取年度地震危险区长、短轴的平均值约为300km)，则标记为1，否则标记为0；对上述标记数组逐列求和，筛选出小于1的列号，其对应的地震为“孤立发生的地震”(图 3)。

(2) 构建由[i，j，d_ij]组成的两两地震间的“离数组”，并根据第三列d_ij的取值判断，若大于“分离距离”，则该行置空，进而按照第三列d_ij取值的大小对“距离数组”进行升序排列。

(3) 如图 4(a)所示，距离最小的为1、2连线，其次为4、5连线，对1、2、3、4号点组成的数组进行组合运算，若新组合中的每一对点间连线的距离均小于“分离距离”，则将上述4个点划分为一个“丛集区”；反之，继续下一组连线重新进行组合运算和距离判断。

(4) 遍历完成第一轮以后，得到一系列地震编号组成的“丛集区一”，重新检索距离数组并将与“丛集区一”中点号相关的所有行置空，重复步骤(3)，直到遍历完整个“距离数组”。

以2018年1月1日—12月31日期间中国大陆地区发生的M_S5.0以上地震为例，利用上述算法扫描得到的丛集区结果如图 4(b)所示，其中，蓝色实心圆点代表应用该算法识别出的“孤立发生地震”，红色实心圆点代表识别出来的“丛集区”内的地震，红色空心圆代表直径为300km的“丛集区”，从中可以看出该改进算法具有较好的识别分类效果。当2018年的危险区个数等于本年度“孤立发生地震”的个数加上“丛集区”个数，即10个时，报准率保持在往年平均值的情况下(≈30%)，理论上的空间占有率达到最低，能够取得最佳的年度预测效果，因此本研究中确定年度地震发生个数的第二种方法为：采用年度地震“丛集区”和“孤立发生地震”的个数之和作为统计结果，结合多种概率分布模型确定其最佳估值。

由于地震的发生在空间上存在一定“丛集性”，采用上述方法和前文中目标震级统计出的年度地震发生个数偏小，本方法中将目标震级下调为中国大陆地区M_S5.0以上。利用上述算法扫描并统计1950—2019年间满足目标震级的各年度地震“丛集区”和“孤立发生地震”的个数之和，分别代入正态分布、Gamma分布、泊松分布、指数分布以及rayleigh分布并进行最小二乘拟合，结合卡方检验(χ²)筛选出的最优分布模型仍然为Gamma分布，拟合得到的数据分布频次及拟合概率密度曲线如图 5所示。

3 结果与讨论

利用2.1中的方法，分别取α=0.15、α=0.10、α=0.05，通过计算拟合Gamma概率密度函数的分位数，得到85%、90%、95%三个置信水平下的年度地震发生个数的上限，以此作为背景地震个数。通过式(6)，进一步计算得出了各置信水平对应的最大危险区个数及其空间占有率，详见表 2。3个置信水平下的最大空间占有率τ均满足年度地震危险区划定规则中要求的小于10%的限制，取h-τ最大的置信度90%这一组结果。在确保最低漏报率的情况下，得出具有统计学意义的年度地震危险区个数上限为9个(数据不含中国台湾)，此为当前预测水平条件下，为获取最高预测评分所建议的每年危险区个数。

利用2.2中的方法，同样取α=0.15、α=0.10、α=0.05，通过计算拟合Gamma概率密度函数的分位数，得到85%、90%、95%三个置信水平下的年度地震“丛集区”与“孤立发生地震”个数之和的上限。通过式(6)，进一步计算得出了各置信水平对应的最大危险区个数及其空间占有率，详见表 3所示。只有95%的置信水平下对应的最大空间占有率τ略大于10%的限制，三组结果基本均满足条件，因此取h-τ最大的置信度95%这一组结果，即13个(数据不含中国台湾)作为年度地震危险区的最佳估值，这与中国地震台网中心年度会商所划定的危险区个数最为接近。

4 结论

本文针对年度地震危险区划分的问题开展了一些尝试性的研究，主要应用目前国内常用的地震预测效能检验方法——R值评分，推导了年度地震危险区个数的评估公式；结合多种统计分布和检验模型，给出了危险区个数的最佳估值，为年度会商中危险区个数的确定提供具有统计意义的参考。由于统计检验环节的输入数据不同，最终得到的年度危险区个数的参考值也略有差别。

(1) 对于第一种方法，即利用年度地震发生个数的统计值估算危险区个数，由于目标震级设定稍高(107°E以东为M_S5.5以上，107°E以西为M_S6.0以上)，得出的年度地震危险区个数的参考值为9个，低于中国地震台网中心目前划定年度地震危险区个数的平均值(约13个)，且偏离较大，原因在于该方法与实际检验方式不符合，其更适合于稍大震级的判定。

(2) 对于第二种方法，鉴于年度地震危险区预测效能检验过程中，通常选取M_S5.0作为起始震级，筛选年度发生地震时将中国大陆东西部地区的目标震级统一设定为M_S5.0以上，同时利用改进聚类算法对空间距离较近的地震进行了“丛集”化扫描，得出的年度地震危险区个数的参考值为13个，与实际工作中使用的个数一致，证明了目前的年度危险区划分方案具有一定的统计意义，且对判定年度M_S5.0以上地震具有较好的应用实效。