基于贝叶斯定理的地震危险性概率预测研究

1966年邢台地震以来，我国有计划地开展了以地震预测为目的的地震科学研究(张国民等，2001)，逐步形成中国特色的长、中、短、临渐进式地震预测科学思路(梅世蓉等，1993；张国民等，2002；刘桂萍，2006)。在地震预测的实践研究中，探索了多种观测技术，包括地震活动性、地壳形变、重力、地电、地磁、地下水流和地球化学等。朱令人(2002)、张国民等(2002)等认为地震综合预测需要积累多手段的前兆连续观测，探测正常背景下的异常变化，分析异常群体的时空强综合特征及其演化过程；需要应用从大量震例经验和理论、实验研究取得的对孕震过程阶段性发展的认识，以及各阶段中异常群体特征的综合判据与指标，对孕震过程进行追踪分析，并对地震发生的时间、地点和强度进行以物理为基础的概率性预测(张晓东等，2003)。许多学者探索了多种地震综合预测手段，分别利用多种预测方法对地震进行综合预测并取得了一定效果(平建军等，2003；李均之等2007；武安绪等，2008；Yu et al，2013；余怀忠等，2013)。由于缺乏对地震孕育发生的物理机理的足够理解，我国的地震预测主要是经验性的(陈运泰，2015)，对于地震危险性的判定主要是以圈定地震危险区的方式(即圈内有震、圈外无震的确定性预测)进行。例如我国年度地震重点危险区相应的预测是圈定未来一年有可能发生“目标”震级以上地震的空间区域，并给出震级范围的确定性预测(庄建仓等，2012)。每一个危险区都是综合考虑了多个学科、单位的预测意见得到的。这种地震预测方式在我国防震减灾事业中发挥了重要作用，但由于地震孕育时间长、孕育过程复杂，地震观测干扰源多且很难识别，地震预测方法的预测效率差异较大，因此以概率的形式来表达地震预测意见的不确定性更加科学合理。事实上，地震概率预测也是目前地震科学中一个迅速发展的领域(陈运泰，2015)。

美国联邦地质调查局(USGS)的加州地震发生率工作组(working group on california earthquake probabilities，简称WGCEP)自1988年以来一直对美国加州的圣安德斯断裂的地震危险性进行研究(WGCEP，1988；Field et al，2014、2015、2017)。在30多年的研究中，该工作组对加州主要潜在震源区的地质构造、应力环境、断层运动状态、历史地震等方面进行了详尽的调查，取得了大量资料(李昌珑等，2016)。2014年WGCEP发布了第3版加州地震破裂预测模型(UCERF 3)，对加州可能发生破坏性地震的震级、震中位置和平均发生率进行了估计(Field et al，2014、2015)。为了服务公众，USGS综合吸收全国多个工作组的成果，在2014年发布了美国未来50年的长期地震风险预测(Petersen et al，2014)，其中包括WGCEP发布的UCERF 3模型。在长期预测模型的基础上，USGS每年会根据上一年的地震活动情况并参考历史地震来进一步计算未来1年的中期地震风险预测(Petersen et al，2018)。

我国学者也采用了多种概率模型对我国多个地震构造区进行地震危险性研究，探讨了不同构造区域的中长期地震概率。董瑞树等(1996)将分段泊松模型用于汾渭地震带的地震危险性分析，闻学泽(1998)、杨明等(2000)分别在四川地区和祁连山中东段地区开展了时间相依的概率地震危险性分析。刘杰(2000)在Vere-Jones提出的应力释放模型的基础上，考虑各部分之间的相互作用，提出耦合应力释放模型，并将该模型应用于华北地区的历史地震目录中。郭星(2015)提出了强震复发的随机特征滑动模型，该模型对弹性回跳理论做了改进，根据应力积累速率恒定的假设，提出了震级和地震发生率都是时间相依的模型，并用于中国大陆部分地震构造区。石耀霖等(2018)从物理的角度给出中国地震数值预测系统框架的设想，其中包括了概率预测。

对背景地震概率的估计是长期预测的主要内容之一，对于中短期地震预测而言，背景地震概率只是一个起点(张天中，1998、1999)。地震的前兆异常能够对中短期的地震危险性判定提供信息。前兆异常概率增益的概念在地震概率预测中有重要意义。贝叶斯定理(Bayes et al，1763；Lee，2012)表明，多个独立的前兆的总概率增益近似为各前兆概率增益之积，而某一前兆异常的概率增益可以通过其预测成功率与前兆时间的比值来表示(Utsu，1979)。Aki(1981)基于贝叶斯公式推导得到前兆异常的概率增益为有前兆条件下的地震概率与无条件背景地震概率的比值。Ogata(2017)基于ETAS模型获取累积地震活动性资料，探讨了地震活动性方面的异常概率增益。

本文结合我国长期以来形成的地震预测工作思路和实践经验，基于测震学、地壳形变、地磁、地下流体、重力等多种单项预测方法给出的短期或年度预测意见，考虑基于泊松分布的预测区域背景地震概率和预测方法的历史预测效能，采用贝叶斯定理计算得到单项预测方法的短期或年度地震危险概率预测结果，进一步采用综合概率方法(Field，2007)，计算得到基于多种单项预测方法的短期或年度地震危险概率预测结果。为保证本研究工作的顺利开展，中国地震台网中心各学科组依据本学科的观测资料及分析结果，每月对我国大陆地震危险性提出短期或年度的预测判定意见。该方法被应用于2018年2~9月中国大陆地震危险概率预测实践并取得了较好的效果，中国大陆约72%的5级以上地震都发生于相对高概率预测区域。

为了用概率形式表达不同预测方法给出的预测区域，首先基于泊松分布假设计算预测区域的背景地震概率；然后，考虑背景地震概率和预测方法的历史预测效能，采用贝叶斯定理计算得到单个方法给出的预测区域地震危险概率预测结果，同时考虑危险区边缘的空间概率衰减；最后，采用综合概率方法给出基于多个预测方法的地震危险概率预测分布。具体计算步骤如图 1所示。

在特定时间范围内，地震的发生受到确定性和随机性双重因素的制约，很难确定性地预测在给定时段内是否一定会发生地震，采用概率表示发生地震的可能性是一种更加合理的选择(闻学泽，1998)。假定某地区(断裂)的地震复发间隔概率密度函数为f(T)(图 2)，自上次地震发生以来的离逝时间为Te，若从上次地震至Te时刻没有目标震级以上的地震发生，那么在Te至Te+ΔT之间发生地震的概率为(WGCEP，1990；闻学泽，1998)

传统的地震危险性概率分析模型假设地震的发生服从泊松分布，从上次地震发生以来的离逝时间Te内发生k次地震的概率为：

式中，k为地震发生次数；λ为单位时间(或单位面积)内随机事件的平均发生率，对于地震事件而言，使用归一化值$\frac{T_{e}}{\overline{T}}$作为λ，其中T为研究区内地震复发间隔的中位数，这样能够消除不同断裂段或不同震级范围地震复发间隔的绝对差别(闻学泽，1998)。在地震离逝时间Te内发生地震(k~1)等价于不发生地震(k=0)的反事件。离逝时间内不发生地震(k=0)的概率计算相对简单

因此，在离逝时间Te范围内发生地震(k≥1)的概率为

实际工作中面临的情况是，已知某区域离逝时间内尚未发生目标震级以上的地震，希望计算得到未来一定时间范围内发生目标震级以上地震的概率，因此，可以基于泊松分布利用公式(1)来计算未来一定时间内的地震发生概率。需要注意的是，式(4)为累积分布函数，即

因此，式(1)可以改写为

根据式(5)可见，基于泊松分布的假设，在未来一定时间ΔT内的地震概率与历史地震复发间隔中位数 有关，与地震离逝时间Te无关，即基于泊松分布的地震发生率不随时间变化，而是时间独立模型。为了保证地震样本数足够多，我们根据历史地震分布以及地质构造将中国大陆进行了分区(图 3)，利用式(6)即可计算各个分区在一定时间范围内的背景地震概率，表 1给出了不同分区的背景地震概率。

为了将单个预测区域的地震危险以概率形式表示，综合考虑预测区域的背景地震概率以及单方法的历史预测效能，采用贝叶斯定理计算得到单方法给出的预测区域地震危险概率预测结果。首先，定义发生地震为事件A，地震发生的概率可表示为P(A)；定义预测地震为事件B，那么预测地震发生的概率为P(B)；发生地震且预测准确的概率可表述为在A发生条件下B发生的概率P(B|A)，即地震有报率(地震发生之前有预测的概率)。实际工作中，我们希望得到P(A|B)，即预测了地震且发生了地震的概率，这里可以利用贝叶斯定理(Bayes et al，1763；Lee，2012)得到

式中，P(B|A)可以根据单项预测方法以往的预测结果来确定，地震发生的背景概率P(A)可以通过1.1节给出的方法计算得到。此外，根据历史预测经验可以计算在预测时间内异常出现的概率P(B)。

另外，为了使概率表达更合理，我们假设地震发生概率在危险区边界外随距离衰减。地震危险概率空间衰减可有多种表达方式，本文暂以余弦平方衰减公式对危险区外围d km内的地震概率进行衰减处理

式中，x为某位置与危险区边界的距离；P为危险区内的背景地震概率；d为概率衰减范围；c为最小衰减值；d和c均为预先给定的常数。例如，危险区内的地震概率为20%，我们希望将危险区外100km处的地震概率衰减到1%，那么P=0.2，d=100，c=0.01，在100km范围内概率以余弦平方的形式衰减(图 4)，100km范围外的地震概率赋予空值。

在实际地震预测中，往往要将不同预测方法提出的多个危险区地震危险概率预测结果进行综合计算。对于同一地点，考虑多种方法预测得到的地震危险概率，计算得到一个综合的地震危险概率(Field，2007)

式中，P为综合概率；P1、P2、P3…为通过用不同方法计算得到的对同一区域的地震概率。

本文研究中所用到的基础资料包括地震目录和各学科的单项预测判定意见。地震目录用于计算背景地震概率，在此基础上进一步结合各学科预测意见计算得到我国大陆中短期的地震危险概率预测结果。另外，我国东、西部地震活动性存在显著差异，西部地震活动明显强于东部，难以用统一的标准展示东西部地震概率，同时由于经济发展、人口密集程度等因素的差异，东部地区对地震震级的敏感度更高。因此，在资料选取与分析计算过程中，西部地区(108°E以西)的预测目标震级为5.5级，东部地区(108°E以东)的预测目标震级为4.5级。

选取1970年以来预测目标震级以上的历史地震目录，并对所选取的地震目录剔除了余震、震群。基于选取的地震目录，利用本文1.1节给出的方法计算各区域的背景地震概率。

为保证本研究工作的顺利开展，中国地震台网中心各学科组依据本学科观测资料及分析结果，每月提出未来3个月、每季度提出未来1年对我国大陆地震危险性的预测判定意见。依据测震学、地壳形变、地磁、地下流体、重力等预测方法提出预测危险区域。2018年2月以来，中国地震台网中心每月产出未来3个月、每季度产出未来1年我国大陆的地震危险概率预测图，实现了动态的年尺度或短期地震危险性概率预测(王海涛，2005)。另外，对分析结果进行跟踪检验，为算法改进提供依据。

选取2018年2~9月产出的8期未来3个月我国大陆的地震危险概率预测图进行检验分析。在计算地震危险概率时，西部地区的目标震级为5.5级，东部地区的目标震级为4.5级，在实际检验中我国大陆在预测期间发生目标震级的地震较少，难以验证所提方法的有效性。综合考虑东西部地震活动情况，将目标震级调整为5级进行检验，同时展示东部4级以上地震。由于我国大陆2018年1~4月未发生5级以上地震，具体检验分析5月份以来产出的地震危险概率预测结果：

(1) 图 5为4月底产出的5~7月地震危险概率预测结果以及大陆西部5级以上、东部4级以上地震分布，预测期内发生5级以上地震2次，分别为2018年5月6日称多5.3级、5月28日松原5.7级地震，均发生在地震危险相对高概率区域。

(2) 图 6为5月底产出的6~8月地震危险概率预测结果以及大陆西部5级以上、东部4级以上地震分布，预测期内发生5级以上地震4次，其中，8月13、14日通海2次5.0级地震均发生在地震危险相对高概率区域，8月3日治多5.1级、8月4日日土5.2级地震未落入地震危险概率预测区域。

(3) 图 7为6月底产出的7~9月地震危险概率预测结果以及大陆西部5级以上、东部4级以上地震分布，预测期内发生5级以上地震8次，其中，8月13、14日通海2次5.0级、9月4日伽师5.5级、9月8日墨江5.9级、9月12日宁强5.3级地震均发生在地震危险相对高概率区域，8月3日治多5.1级、8月4日日土5.2级、9月28日日土5.1级地震未落入地震危险概率预测区域。

(4) 图 8为7月底产出的8~10月地震危险概率预测结果以及大陆西部5级以上、东部4级以上地震分布，预测期内发生5级以上地震10次，其中，8月13、14日通海2次5.0级、9月4日伽师5.5级、9月8日墨江5.9级、9月12日宁强5.3级、10月16日精河5.4级、10月31日西昌5.1级地震均发生在地震危险相对高概率区域，8月3日治多5.1级、8月4日日土5.2级、9月28日日土5.1级地震未落入地震危险概率预测区域。

(5) 图 9为8月底产出的9~11月地震危险概率预测结果以及大陆西部5级以上、东部4级以上地震分布，预测期内发生5级以上地震8次，其中，9月4日伽师5.5级、9月8日墨江5.9级、9月12日宁强5.3级、10月16日精河5.4级、10月31日西昌5.1级、11月4日阿图什5.1级地震均发生在地震危险相对较高概率区域，9月28日日土5.1级、11月26日台湾海峡6.2级地震未落入地震危险概率预测区域。

总的来看，2018年2~9月中国地震台网中心产出的8期未来3个月我国大陆的地震危险概率预测结果检验如图 10和表 2所示，每期MS≥5.0地震命中数平均为72.5%。综合分析上述结果表明，该方法在地震预测实践中具有一定的预测效能，但存在一定的区域差异性，新疆、川滇等多震区域的预测检验效果好于其他地区，表明背景概率在其中发挥着重要作用。

(1) 2018年2~11月中国大陆72%的5级以上地震都位于相对高概率预测区域，表明本文基于贝叶斯定理的地震危险概率预测是一种有效的地震短临预测方法。这种基于概率的地震预测结果相比于确定性的预测更具科学性，能够更直观地展示各地区地震发生的危险性，有利于防震减灾工作的决策。同时，通过不断的实践检验积累，能够对不同学科预测方法进行更加完备的检验，为预测方法进入预测业务工作提供了好的检验平台。

(2) 本文在计算危险区所在地区的背景地震概率时，采用了大的分区平均，不能体现具体的断层构造和多震或少震特点，后续可考虑断层构造特征和地震分布，将全国网格化计算地震背景发生概率。

(3) 本文在计算危险区所在地区的背景地震概率时假设地震的发生服从泊松分布，这是一种时间独立的地震危险性概率分析，但是许多地区(断层)的地震活动性随时间有较明显的变化(杨明等, 2000, 李昌珑等，2016)，后续可考虑采用时间相依模型计算背景地震概率(Field et al，2015)。

(4) 目前对预测效果的检验只考虑了落入相对高概率区域的地震比率，后续应进一步考虑风险高概率区面积占有率。