基于贝叶斯概率的强震综合概率预测指标建立及其在四川地区的应用
2022-09-25何畅,路茜,龙锋,宫悦,李丽
何 畅,路 茜,龙 锋,宫 悦,李 丽
(1.四川省地震局,四川 成都 610041;2.山西省地震局,山西 太原 030021)
长期的地球物理场观测与地震预报实践表明,由于地球的不可入性以及地震孕育过程的复杂性(陈运泰,2008),单一手段的地震预测方法往往很少有较高的信度以支撑起对某次强震的准确预测,采用算法联合多种监测预报手段进行综合概率预测是一条可行探索的途径,已有众多基于不同算法的实验性工作得以开展并取得成效,如支持向量机回归方法(蒋淳等,2006;卫定军等,2014)、概率增益建模(邱玉荣,2011)、遗传神经网络(陈秀琼等,2005)、模糊聚类法(狄莉莎等,2000)等。利用多手段综合预报,可以“抹平”各手段的随机预测噪声,同时可在一定程度上增强前兆信号的强度,以提高预测水平。
在设计综合概率预测方法的过程中,如何保证各单一指标以不同的权重(信度、预测效能等)参与进来,并产出符合先验的、相对的强震发震期望(概率或某些可度量的数值)是这一方法的核心部分。本文在对四川地区开展地震预测实践中,从贝叶斯概率出发,推导出了单个指标的发生强震的后验概率,采用联合概率计算所有异常的综合概率,并对2013~2019年四川区域的3次6级以上强震进行了回溯性研究,并尝试提取不同震级的综合概率值。
1 方法
不同于古典频率学派,贝叶斯概率公式侧重于基于已有的先验知识预测某件事情发生的概率(四川大学数学系高等数学教研室,1990),已经在参数回归(Hoff,2009)、风险决策(王学军等,2003)、机器学习(茅伟强,2008)等方面得到广泛应用。针对单个地震预测指标的贝叶斯概率的公式为:
式中:Y表示异常,D表示地震。P(Y)和P(D)分别是该指标异常出现以及地震发生的先验概率,P(Y|D)是地震发生前必然存在异常的条件概率,是似然值。而所求的P(D|Y)是异常出现后必然出现地震的条件概率。由于地震预测是二元分类问题——即有异常、无异常、有地震、无地震,因此依据全概率公式,异常出现的总概率P(Y)共有两项,可以写为:
式中:N表示无地震。为了方便计算,同时考虑到异常和地震发生的低概率性,将它们离散到时间上(以天为单位),用时间的占有率来代替各自的概率。设总的研究时段为Tall,地震发生次数(或天数)为Neq,异常出现次数(或天数)为Na,异常出现后发生预测中的地震次数(即“报准”的次数或天数)为Nb,则:
由(3)、(4)两式可得:
另有:
P(Y|N)为未发生地震时出现异常的概率,可表示为未发生地震时未出现异常的概率的补集,即:
式中:X为未出现异常的时段,t为预报指标中单次预报时长。联合(5)~(7)式,最终可得:
上式表明,某个指标异常出现后对应地震的后验概率可从其具体的预报指标和历史验证效能中获得。进一步将式(8)上下同时除以Na,则可得到:
式中:Rc为报准率,从上式可以看出高的报准率和短的预报时长会明显提升预报概率,这与实际相符。另外也可看出式(9)存在分母为0的函数奇点。由于Rc∈[0,1],而t往往大于1,因此分母是否为0取决于t与Neq/Na的大小关系,也就是说,当预报时长与研究区总地震数和异常出现次数的比例接近时,概率值可能为无穷;且当这一比例明显大于预报时长t时,概率值为负数。
2 数据
2.1 单测项月尺度概率
针对四川地区目前仍在正常运作的各观测项进行了6级以上历史震例的回溯性研究,提取了能通过R值检验(许绍燮,1973①许绍燮.1973.地震技术资料汇编:震兆分析一例.北京:《地震战线》.;许绍燮等,1989)的预测指标,最终获得87个测项指标的贝叶斯后验概率计算结果(表1)。在常规的地震预测实践中,我们更期望获得月这一时间尺度上的强震发生概率,因此在通过式⑻计算概率时,往往将指标的预测时间t归算到月。在表1中,地震活动性的部分异常如空区、条带、震群、b值、震源机制解一致性空间扫描、以及宏观异常或其它临时出现但具备一定预报意义等共计9项指标由于无法溯及过往总数,无法计算P(D|Y)值,依相应学科专家依经验给定数值;另有9项指标由于P(D|Y)值计算为负数或无穷大,人为给定该值为0.010(表1中带②的指标项)。
通过对P(D|Y)进行统计,所有指标(图1a)与不考虑人为给定概率的指标(图1b)在统计分布上并无明显差异,后验概率在0.1~0.2是优势概率区间,表明四川地区大部分预测指标在月尺度上对6级以上地震的预测并不具备强约束。具有较高预测概率的3个指标分别是巴塘川-58泉(305K)水温预测概率为0.33、普格川-65泉水温预测概率为0.6和康定二道桥川-55泉水温预测概率为0.75,均属于地下流体测项。从表1中震例回溯的统计数据来看,这3个测项的报准率并不高(占50%或更低),显著低于部分跨断层形变测量测项,如尔乌水准3-A测边(85%以上)、龙灯坝基线A-D测边(100%)等,且有着偏高的虚报率,但由于其较短的预报期限(3个月),使其在月尺度上有着较高预测概率值。
图1 P(D|Y)值所有指标(a)和不考虑认为给定概率的指标(b)统计柱状图
表1 四川地区强震指标及相关震例
续表1
续表1
续表1
2.2 典型震例回溯研究
考虑到历史强震震例存在部分地球物理场观测项失效的问题,这里选择针对2013年以来四川地区的3次6级以上震例进行了概率分析(表1)。统计结果显示,2013年芦山7.0级地震前存在异常40项,其中预报期限在6个月以内的中短期异常12项。为了综合考量所有异常指标对未来强震的预测概率,采用了联合概率的方法,设由式(8)计算得到的第i个指标的预测概率为P ci,则n个异常指标在月尺度上6级以上强震的发震概率可表示为:
根据式(10),若假定P ci=0.01,需229项异常可使P值达到0.9,而当P ci=0.1时,则需22项异常。
以2013年芦山7.0级地震而言,震前40项异常计算出来的发震概率为0.992。绘制发震概率随异常项数变化的曲线绘制出来(图2中蓝色实线),可以看出第7项异常(例如康定二道桥川-55泉水温,P c=0.75;普格川-65泉水温,P c=0.60等)对提升总体概率贡献极大,从0.45提升至0.86。当异常项数达到15项而概率超过0.90后,概率曲线平缓,这是因为式(10)为极限为1的幂函数,这意味着尽管异常项数增加,但最终发震概率值变化幅度较小,通常在小数点后两位变化,这无益于异常信号的识别和预报指标的提取。为此,本文试图从两个方面解决这个问题:(a)仅使用中短期异常。芦山7.0级地震前12项中短期异常的发震概率同样达到了0.92(图2中红色实线),这样既保证了概率值可达到较高的水平,同时相对较少的异常数量使得相邻概率值之间的区分度较高。这种方法并不能保证未来随着监测能力的提升,更多的中短期异常项数依旧会使得概率曲线变得平坦。(b)考虑到概率值在0~1之间变化,通过反正弦函数将概率值转换为角度,即:
图2 芦山地震前概率值和角度值随异常项数变化曲线(实线为概率,虚线为角度)
可在一定程度上消弭概率曲线的平坦效应,增加概率的区分度。从图2可以看出,对于芦山7.0级地震前的40项异常,相较于概率值在第15~40项异常之间仅有0.07的变化幅度,相应区间的角度变化值有15°(蓝色虚线),方便于指标提取时阈值的确立。芦山7.0级地震前40项异常的角度值为82.79°(蓝色虚线),如果只考虑12项中短期异常,则角度值为68.43°(红色虚线)。从实际应用的效果来看,四川地区背景性异常长期维持在30项左右,这一数量使得发震概率P值始终在0.9左右,不利于区分当前强震的紧迫性。因此,仅采用6个月以内的中短期异常来计算概率符合当前按月来追踪震情形势的需求,且概率变幅较大有助于强震紧迫性的区分。
本文收集整理了2014年康定6.3级和2019年长宁6.0级地震前的中短期异常,并计算了它们的发震概率值和概率角度值,其中康定6.3级地震前中短期异常5项,概率和概率角度值分别为0.52和31.10°,而长宁6.0级地震前中短期异常4项,概率和概率角度值分别为0.48和26.08°。
3 结论
为了建立四川地区强震综合概率预测指标,从贝叶斯概率出发,推导出了单个指标出现异常后强震发生的后验概率;基于联合概率,引入多项指标出现异常后的综合概率,并对芦山7.0级、康定6.3级和长宁6.0级等强震开展了回溯性研究。获得的结论如下:(1)单个指标的贝叶斯后验概率的大小与报准率正相关、与预报时长负相关,这是符合经验直觉的。但该概率同时也被预报时长和总地震数及总预报数的比例的相互大小所约束,在计算过程中要规避极小可能出现的无穷或负值。(2)四川地区强震指标的后验概率大多在0.2以内,说明大多数指标短期预测效能一般,对未来强震约束不强。(3)联合概率是一个极限为1的幂函数,当异常数较多时,少数异常的增减在概率上的区分度不高,不便于阈值的选取和指标的提取。为此,建议只采用6个月以内的中短期异常进行联合概率计算,又或者,通过反正弦函数将概率值转换为角度值以获得更好的区分度。(4)采用以上规则,震例回溯研究统计出:芦山7.0级地震前有12项中短期异常,计算得到的联合概率为0.92,概率角度值为68.43°;康定6.3级地震前中短期异常5项,概率和概率角度值分别为0.52和31.10°;而长宁6.0级地震前中短期异常4项,概率和概率角度值分别为0.48和26.08°。(5)通过这三次较少的震例,本文尝试性提出了四川地区6级以上地震的综合概率预报指标为:中短期异常在4项以上,联合概率值大于0.45,概率角度值大于25°;而7级以上地震的综合概率预报指标为:中短期异常在10项以上,联合概率值大于0.90,概率角度值大于65°。(6)四川地区选取6级以上强震震例较少,且地震发生和地球物理场的监测能力均存在时空不均匀性,因此需要开展更多的震例回溯总结和积累工作,以形成更稳健的综合概率预报指标体系。
致谢:本文在数据收集整理过程中,得到了四川地震台张致伟、官致君、苏琴、祁玉萍、马伶俐、芮雪莲、林圣杰、王迪等提供的大力帮助;中国地震局地球物理研究所蒋长胜研究员提供了部分思路,在此一并致谢。