基于PI算法的青藏块体强震回溯性研究

余 娜 张晓清 袁伏全 赵燕杰

（中国西宁 810000 青海省地震局）

0 引言

近年来，基于统计物理学与地震学的图像信息学（Pattern Informatics，简称PI）算法在地震较活跃地区的中长期预测研究中得到了广泛应用。Rundle 等（2002）将PI 算法应用到美国南加州地区地震的预测研究中，认为该方法对加州地区地震活动有较好的预测效果；国内一些研究者也将PI 算法应用到中国台湾（Chen et al，2005）、中国大陆西部（蒋长胜等，2008；Zhang et al，2010，2013；Jiang et al，2010，2011，2013；张小涛等，2014；袁伏全等，2016）等地区地震的研究中，认为PI 算法的预测效果显著优于随机概率法。Zhang 等（2013）和张小涛等（2014）利用PI 算法分别对2008 年汶川8.0 级、2014 年于田7.3 级、2013 年芦山7.0 级地震进行地震危险性的回溯性研究，结果表明模型计算参数的选取对PI 算法的预测结果有一定影响，对于7 级以上大地震，选择较大的网格尺度和较长的预测时间窗可能会取得更好的预测效果；袁伏全等（2016）利用PI 算法对青海地区中强地震进行地震危险性的回溯性研究，通过采用不同的网格尺度进行回溯性预测检验，认为对于青海地区采用0.2°×0.2°网格尺度时预测效果更好。

印度板块对欧亚板块的挤压作用使青藏块体区域构造运动强烈，地震频发，中国大陆7 级以上强震主要发生在青藏块体及邻近区域（吴哲，2018）。本文研究区（30.0°—41.0°N，88.0°—105.0°E）包括羌塘活动地块、巴颜喀拉活动地块、柴达木活动地块和祁连活动地块等4 个Ⅱ级活动地块。活动块体边界带是我国大陆强震的主体带和集中区，因此，研究活动地块边界带的强震趋势显得尤为重要（张浪平等，2010）。本文在前人工作的基础上，利用PI 算法对研究区MS≥6.0 地震进行回溯性研究，对不同时间尺度下的预测结果进行ROC 检验，并给出适合该地区的模型计算参数。

1 计算方法

1.1 PI 算法

PI 算法的实质是通过对研究区地震活动增强和平静进行分析，计算并预测在中长期时间尺度上可能发生中强地震的概率。主要步骤：首先，将研究区按一定尺度划分空间网格，对落入网格的不小于截止震级的地震事件构建时间序列；然后，对地震活动强度变化进行归一化处理，计算每个网格发震概率；最后，发震概率减去背景概率得到发震概率较高的区域，即“PI 热点分布图”（Tiampo et al，2002；Nanjo et al，2006；蒋长胜等，2008）。具体计算步骤如下。

（1）对研究区进行空间网格划分，每一网格为xi。

（2）对落入网格且不小于截止震级的地震事件的网格构建1 个时间序列Ni(t)，其中，Ni(t)为中心坐标和相邻的8 个格点（Moor 近邻）（Moore，1962）在时间t 内单位时间的地震数目；t0为研究区地震目录的起始时刻；tb为滑动变化各时间序列的起始时刻；t1为地震活动中异常学习的起始时刻；t2和t3分别为预测时段的起始时刻和终止时刻（图1）。

图1 空间网格划分及时间序列构建Fig.1 Division of spatial grid and construction of time series

（3）地震活动强度函数Ii（tb，t）为从时刻tb到t 的单位时间内发生在网格i 中的不小于截止震级的平均地震数

（4）对不同时间段的地震活动强度求平均，再除以标准偏差，进行标准化处理

（5）计算地震活动强度函数平均变化量，以较少随机扰动的影响

（6）未来强震发生概率Pi（t0，t1，t2）为地震强度函数平均变化量的平方。第i 网格的概率值减去所有网格概率的平均值得到强震发生在第i 网格的概率

将ΔPi(t0,t1,t2) ＞0的格点视为研究区的地震热点。在PI 图像中将地震活动性较强或地震频次较高的30%格点画出（陈建志等，2013），由其反映预测时间窗内空间上发生目标地震概率较高的区域，用lg(ΔP/ΔPmax)的值表示其发震概率。

1.2 ROC 检验

ROC 检验方法是在国际合作项目（CSEP）中进行科学预测试验的方法之一，其目的是检验预测结果的优劣程度。ROC 检验既要考虑“命中率”，同时也要考虑“虚报率”。有意义的预测必须是ROC 值大于0.5。ROC 值为ROC 曲线与随机预测曲线所围成的有效面积，ROC 值越大，预测效能越好。

2 研究区域及模型参数

2.1 研究区域及资料选取

选取青藏块体为（30.0°—41.0°N，88.0°—105.0°E）研究对象（图2）。历史上该区域曾发生过多次7 级以上强震。本文使用中国地震台网中心提供的1970—2019 年8 月全国ML≥2.0 地震目录，共计20 896 次地震。其中，ML2.0—2.9 地震3 286 次；ML3.0—3.9 地震14 611 次；ML4.0—4.9 地震2 307 次；ML5.0—5.9 地震605 次；ML6.0—6.9地震74 次；ML≥7.0 地震13 次。Jiang 等（2011）基于川滇地区强震的回溯性研究讨论了强余震对PI 算法的影响，认为余震对PI 算法计算结果的影响时间不超过1 年；Tiampo 等（2002）认为震后余震反映了区域高应力的释放。因此，本文没有对地震目录进行剔除余震处理。

图2 研究区空间位置矩形框为研究区域Fig.2 The spatial location of the study area

截止震级Mc是PI 算法中重要的参数之一，而截止震级Mc的选取与该地区的最小完备震级有关（Woessner et al，2005）。在计算中，若截止震级Mc选取过高，则参与计算的地震数目会减少；若截止震级Mc选取过低，则有些区域因监测能力较弱而记录不到地震。因此，本文将截止震级Mc设定为ML3.0，这既能保证充足的地震数据量，又能满足截止震级应至少小于目标震级2 个震级单位的要求（Holliday et al，2005，2006）。

2.2 模型计算参数设置

一 些 研 究 者（Rundle et al，2002；Chen et al，2005；Nanjo et al，2006；Jiang et al，2010；袁伏全等，2016）利用PI 算法对中国大陆西部地震进行回溯性检验，认为网格参数选取为0.2°×0.2°，预测时间窗选为8 a 时，PI 算法的预测效果较好。本文参照前人的研究结果（Rundle et al，2002；Chen et al，2005；Nanjo et al，2006；Jiang et al，2010；袁伏全等，2016），选取0.2°×0.2°空间网格尺度，以30 d 为步长连续向前滑动，目标震级≥MS6.0。研究区域中强地震频发，2009 年8 月至2019 年8 月研究区共发生9 次MS≥6.0 地震（图3，表1），为了考察研究区MS≥6.0 地震在不同模型参数下的回溯性预测效果，在研究区范围、网格尺度、截止震级以及滑动步长不变的情况下，选取预测时间窗为10 a、5 a、3 a 进行回溯性研究，模型参数设置结果见表2。

图3 2009 年8 月至2019 年8 月研究区MS≥6.0 地震震中分布Fig.3 Distribution of MS≥6.0 earthquakes since August 2009

表1 2009 年8 月至2019 年8 月研究区MS≥6.0 地震参数Table 1 MS≥6.0 earthquakes parameters in the study area since 2009

表2 模型计算参数设置Table 2 Selection of calculation parameters

3 计算结果

根据表2 的模型计算参数，得到了3 个回溯性预测检验时段的PI 地震热点分布图（图4）。

图4 研究区回溯性PI 预测图像及ROC 检验结果（a）、（b）2009 年8 月1 日—2019 年8 月1 日；（c）、（d）2014 年8 月1 日—2019 年8 月1 日；（e）、（f）2016 年8 月1 日—2019 年8 月1 日Fig.4 PI retrospective forecast and ROC test for the research area

选定研究区地震目录的起始时刻 t0为1970 年1 月1 日，分别选取地震活动中异常学习的起始时刻t1为1999 年8 月1 日、2009 年8 月1 日和2013 年8 月1 日，预测时段的起始时刻t2分别为2009 年8 月1 日、2014 年8 月1 日 和2016 年8 月1 日，预测时段的终止时刻t3为2019 年8 月1 日，因此共3 个回溯性预测时间段。图4（a）为10 a 预测时间窗（2009 年8 月1 日至2019 年8 月1 日）内的地震热点分布。由图4（a）可见，异常出现在大柴旦—宗务隆山断裂、鄂拉山断裂、达布逊湖—霍布逊湖断裂、祁连山北缘断裂、岷江断裂和汶川—北川断裂附近。根据预测规则可知，这些区域是预测时间窗内空间上相对危险的地区，发生目标地震的概率较高。实际上，该预测时间窗内研究区共发生了9 次MS≥6.0地震，其中有4 次地震震中落在由PI 算法计算得出的地震热点区域，分别为2009 年海西6.4级、2013 年芦山7.0 级、2016 年门源6.4 级、2017 年九寨沟7.0 级地震。

从5 a 预测时间窗（2014 年8 月1 日至2019 年8 月1 日）内的PI 预测图像[图4（c）]可见，2016 年门源6.4 级、2017 年九寨沟7.0 级地震震中均位于PI 预测图像的地震热点区域；从3 a 预测时间窗（2016 年8 月1 日至2019 年8 月1 日）内的PI 预测图像[图4（e）]可见，2017 年九寨沟7.0 级地震震中位于PI 预测图像的地震热点区域。

采用ROC 方法回溯性检验不同预测时间窗的预测效果，发现10 a 时间尺度的ROC 值为0.610，5 a 时间尺度的为0.725，3 年时间尺度的为0.793。可以看出，预测时间窗内的强震震中基本都在PI 预测图像的地震热点区域，而且PI 算法预测结果明显优于随机概率法[图4（b），4（d），4（f）]。虽然可供参考的强震事件数较少，但在对PI 算法的回溯性检验中发现，2009 年海西6.4 级、2016 年门源6.4 级、2013 年芦山7.0 级、2017年九寨沟7.0 级地震震中都落在地震热点的丛集区，这表明PI 算法对青藏块体的强震具有一定的预测效能，并且将预测时间窗设定为3 a 时间尺度时其预测效能较好。

4 结论与讨论

2009 年8 月至2019 年8 月研究区共发生9 次MS≥6.0 地震，其中，6 次发生在活动块体边界带上（图3），分别为2009 年海西6.4 级、2010 年玉树7.1 级、2013 年芦山7.0级、2014 年康定6.3 级、2016 年门源6.4 级、2017 年九寨沟7.0 级地震。考虑到研究区强震大部分发生在活动块体边界带，且活动块体边界带也是我国大陆强震的主体带和集中区，因此定性分析了活动块体边界带“目标”地震与地震热点丛集间的关系，发现青藏块体MS≥6.0 地震与活动块体边界带上地震热点间对应关系较好，活动块体边界带的“目标”地震都发生在地震热点附近区域，而在活动块体内部的MS≥6.0 地震震前无地震热点情况。仅就PI 算法而言，很难对该现象给出具体原因。但在地震中长期预测中应用PI 算法时，将其计算结果与活动块体边界带结合起来考虑强震相对危险区域，在地震预测研究中可能有一定参考价值。

本文利用中国地震台网中心提供的1970 年以来全国地震目录，在分析研究区最小完备震级的基础上，应用PI 算法，选取空间网格尺度0.2°×0.2°，截止震级为ML3.0，目标震级≥ML6.0，将预测时间窗设定为10 a、5 a、3 a 不同时间尺度进行回溯性研究，并对预测结果进行ROC 检验，发现PI 算法对青藏块体强震有一定预测效能，预测效果明显优于随机概率法，并且将预测时间窗设定为3 a 时间尺度时预测效能较好。