前人的研究中对PI 方法已有详细介绍,本文只简述其计算方法。主要分为4个步骤:①首先对选定的研究区域进行时空网格化,参与计算的目录数据其截止震级M_C在不低于最小完备性震级的同时至少要小于“目标地震震级M_T”两个震级单位(Holliday et al,2005; Zhang et al,2013)。②对落入相应网格内的地震活动频度构建其时间序列,并定义为3个时段(图1)。通常设定研究资料起始时刻为t₀; 滑动变化的时间序列起始时刻为t_b(t_b≥t₀), t_b从t₀时刻向[t₁-(t₂-t₁)]时刻滑动(t_b<t₁); 背景参考时段为[t_b,t₁], 学习(考察)时段为[t₁,t₂], 预测时段为[t₂,t₃]。背景参考时段通常要偏长一些,后两个时段通常取t₂-t₁=t₃-t₂。③对每个网格的地震频度时间序列进行时空归一化处理并计算其相对于背景值的偏离程度。④用某个网格的发震概率减去所有网格发震概率的平均值即为对应网格的概率增益ΔP_i。利用lg(ΔP_i/ΔP_max)对概率增益进行归一化,获得“PI热点图像”。最终以热点的形式表示出研究区域内预测窗时段内目标地震(M≥M_T)可能发生的位置。根据lg(ΔP_i/ΔP_max)的下限值来约束热点数量(Holliday et al,2005; Nanjo et al,2006a; Kawamura et al,2013; Zhang et al,2013)。在研究过程中所使用的地震数据资料,其背景参考时段和考察时段的地震目录一定是已实际监测到的,即t₂之前属于“过去式”。而[t₂-t₃]既可能是过去时段,也可能是未来时段。当预测窗[t₂-t₃]处于过去时段,则为回溯性预测研究(可检验); 当[t₂-t₃]处于未来时段,则只能预测,并不能检验。

图1 PI方法时间序列节点示意图
Fig.1 Schematic diagram of the nodes of time series of PI method

针对2021年青海玛多M_S7.4地震,选取(30°～40°N,94°～106°E)为研究区,并从中国地震台网中心获取1980年1月1日至2021年5月22日玛多地震前研究区的地震目录。近20年来中国大陆7级以上地震主要发生在巴颜喀拉地块周缘,但玛多M_S7.4地震发生于巴颜喀拉地块内部,本文主要针对研究区内4个二级块体(图2)进行分析讨论。由于PI方法计算中的各项参数对于热点效果都存在一定的影响,需要根据震例实际情况对参与计算的震级下限、震源深度、网格尺度、时间窗长等进行选择设置(Zhang et al,2013; 张小涛等,2014; 宋程等,2018; 余娜等; 2021)。

在袁伏全等(2016)对青海地区M_S≥5.0目

图2 研究区内二级块体及2008年以来7级以上地震分布
Fig.2 Distribution of the secondary blocks and the MS≥7 earthquakes in the study area since 2008

标地震和余娜等(2021)对青藏块体(30.0°～41.0°N,88.0°～105.0°E)M_S≥6.0目标地震的回溯研究中,认为在网格尺度0.2°×0.2°、预测时间窗长3 a、截止震级M_C≥M_L2.8(1980年以来,满足不低于最小完备性震级)的参数组合下热点效能较好。最小完备性震级可以有效地反映地震台网的监测能力(Schaff,2008; 冯建刚等,2012),而输入目录数据的M_C是PI方法计算中直接影响热点效果的重要参数。M_C偏高会导致输入数据量不足,偏低则会导致输入量过多降低热点对于目标地震的敏锐性。蒋长胜(2010)研究表明,对于板内地震研究的模型参数选取中设定M_C低于“目标地震震级M_T”两个震级单位的常用标准可能会出现问题。而本文针对玛多M_S7.4地震预设目标地震下限为M_T=M_S7.0,明显高于5～6级震例,在保证M_C+2≤M_T的情况下理论上最高可取M_S5.0(Rundle et al,2002; Tiampo et al,2002a; Holliday et al,2005; Nanjo et al,2006a)。为进一步选定合适的M_C,分时段画出1981年至2021年5月21日研究区的G-R关系图(图3),从图3中能够看到4个时段研究区内整体数据拟合较好。但图3b、d中在M_L4.8处局部拟合线与整体拟合线的斜率出现偏转,故综合考虑后本文选择输入震级下限为M_L4.8,亦可保证数据量充足。

前人应用PI方法研究中国西部强震用到的网格尺度参数各有不同,包括青藏地区1.0°×1.0°(夏彩韵等,2015; Zhang et al,2017)和0.2°×0.2°(袁伏全等,2016; 余娜等,2021)、新疆局部地区0.4°×0.4°(李志海,黄瑜,2010; 尼鲁帕尔?买买吐孙等,2018)和0.5°×0.5°(尼鲁帕尔?买买吐孙,2020)、川滇地区1.0°×1.0°(蒋长胜等,2009; Zhang et al,2013)和2.0°×2.0°(张小涛等,2014)。参考前人在应用PI方法研究中国大陆西部及国内外7级以上震例时的参数设置(蒋长胜,2010; Zhang et al,2013,2017; 蒋卉等,2013; 张小涛等,2014; 夏彩韵等,2015; 宋程等,2018),通常网格尺度较大且预测窗长相对较长时,热点效果更为集中和显著(Zhang et al,2013,2017; 张小涛等,2014; 宋程等,2018)。蒋长胜(2010)对于2008年5月12日

图3 1981-01-01—2021-05-21研究区域内分时段G-R关系图
Fig.3 The G-R relation in different periods in the study area from Jan.1,1991 to May 21,2021

汶川M_S8.0地震的PI算法回溯性研究结果表明,震前PI热点会因为网格参数的变化而覆盖到主震破裂区,说明网格尺度要因地、因震制宜。在确定目标地震震级下限和截止震级后,考虑到玛多M_S7.4地震的位置特殊性、网格内计算数据量充足性以及热点迁移性,本文设定网格尺度为0.8°×0.8°。

玛多地震距中国大陆上一次7级以上强震(2017年8月8日四川九寨沟M_S8.0地震)发震时间间隔3.4 a,距青海省内上一次7级以上地震(2010年4月14日青海玉树M_S7.1地震)发震时间间隔11.2 a。7级地震孕震机制较为复杂,结合之前不同学者对于国内中强震例的研究经验(Zhang et al,2013; 张小涛等,2014; 夏彩韵等,2015; Zhang et al,2017; 宋程等,2018),本文设定时间窗长为8 a,阈值选择前人惯用的lg(ΔP_i/ΔP_max)=-0.6(Zhang et al,2013; 张小涛等,2014; 夏彩韵等,2015; 宋程等,2018)。

通常选取与研究的目标地震震源深度相接近的范围内的事件参与计算,更能符合其构造环境、孕震过程、热点追踪的目的和效果。考虑到本次玛多地震为浅缘地震(震源深度17 km)以及保证数据量的充足,本文选择60 km以内的浅源地震事件参与计算。

考虑到余震也是应力释放的一种表现,删除余震可能会造成地震目录一定程度上的畸变(张小涛等,2014),故本文并未将其删除。本次针对的目标地震震级为M_T≥M_S7.0,参与计算地震的最小震级为M_L4.8,删除目录中的较小余震对于计算结果无实质性影响,删除M_L>4.8的所有余震则可能无法恰当的反应局部区域的热点特征。Tiampo 等(2002b)认为余震反映了局部区域应力的高水平释放; 在物理上,较小的地震事件仍能携带较为丰富的强震危险性增强的信息(蒋长胜,2010)。

图4显示2009—2020年研究区内5个连续预测窗口中的热点变化情形,青海玛多县附近持续出现成片密集热点。虽然玛多地震发生在2021年,但在5个预测窗口中,玛多地震震中所在网格及其3个摩尔邻近网格持续存在热点,即这5个回溯性预测窗口中显示了明确的高发震风险覆盖范围,具有强震发震危险性地点指示意义。

图4中分布于巴颜喀拉块体南边界的2010年青海玉树M_S7.1地震、东边界的2013年四川芦山M_S7.0地震和2017年四川九寨沟M_S7.0地震并未出现过有效PI热点。 甚至整个研究区内34°N以南地区几乎无热点,只在2011—2018年预测窗口内青海中南部出现过2个热点。Zhang等(2013)对中国西部大陆(20°～50°N,70°～110°E)进行PI 方法回溯性预测研究时,发现2008 年于田 M_S7.3地震和汶川M_S8.0地震震中存在热点。张小涛等(2014)对川滇地区(21°～35°N,97°～109°E)进行 PI 方法热点回溯性预测研究得出,芦山地震震中落入热点网格内。而图4中芦山地震和九寨沟地震未出现有效热点的原因,除了与所选择的研究区主体范围不同(本文选择研究区相对较小)外,还和PI方法对于模型参数的依赖性有关。张小涛等(2014)对芦山M_S7.0地震前的PI热点特征研究中设置的主要参数为2°×2°/1°×1°网格尺度、9 a预测时间窗长。Zhang等(2013)针对2008年汶川M_S8.0和于田M_S7.3地震的PI热点回溯性研究中,选取的主要参数为2°×2°/1°×1°网格尺度,预测窗长为4～8 a的ROC检验效果较好。以上研究中截止震级皆为M_L4.5,不低于区域最小完备性震级,但低于本文所选的最小输入震级M_L4.8,数据量可能更大。本文所选的预测窗长为8 a,跨度时长接近或重叠于上述时间参数。差异最大的参数为网格尺度,本文设定为0.8°×0.8°。模型参数的差异性,可能是本文研究中其他7级以上地震缺少有效PI 热点的主要原因。

图4 2009—2020年研究区连续热点演化分布图
Fig.4 Retrospective forecast of the hot spots distribution by PI method in the study area from 2009 to 2020

按照相同的时空参数,从2014年开始向未来滑动,则包含玛多地震发震时间(2021年)的窗口有8个。从时间先后顺序来看,8个预测窗口中玛多地震的相关PI热点呈现为“出现—减弱—消失—强密集—收缩”的趋势变化。

在前3个窗口中(图5a～c)热点相对分散,且主要分布于研究区的西北部,但玛多地震震中或摩尔邻近网格出现了热点。相对于图4预测窗口

中部相对稳定连续的密集PI热点,图5中玛多地震附近的热点分布呈不连续、不稳定、丛集性不明显特征。在2017—2024年、2018—2025年2个预测窗口中(图5d、e),预测窗口起始时间距离发震年份为4～5 a,即发震年份处于窗口近乎中间的年份,玛多地震及其摩尔邻近网格未出现相关热点,研究区内热点主要分布在青海最南部,玉树藏族自治州囊谦县及玉树县(青川藏交界附近)。在 2019—2026 年预测窗口中(图5f),玛多地震震中及其摩尔邻近网格出现明显密集深色热

图5 2014—2028年研究区连续热点演化分布图
Fig.5 Retrospective forecast of the hot spot distribution by PI method in the study area from 2014 to 2018

点,且只有2个零散热点,热点整体丛集性较强。至最后2个窗口时段(5g、h),玛多地震摩尔邻近网格密集热点颜色几乎无变化但覆盖面积有所收缩,彰显了发震紧迫性的增强,对玛多地震的发震时间及位置具有一定的指示意义。

在研究区内13个连续滑动的预测窗热点图像中,只有2013—2020年预测窗口出现了19个热点,其余12个预测窗口中热点总数不超过15个。在图5中,只有2个窗口在玛多地震震中或摩尔邻近网格皆无热点,有效热点占比数最高的时段为2019—2026年预测窗。

表1 2014—2028年玛多地震相关热点统计
Tab.1 Statistics of the hot spots related to the Maduo earthquake from 2014 to 2028

从2009年以来研究区内的13个连续滑动预测窗热点图像(图4、图5)中可见,巴颜喀拉地块内部始终存在PI 热点,说明其地震活跃程度相对于其他次级地块更为活跃。其中11个预测窗中,巴颜喀拉地块和柴达木地块交界处存在热点; 在羌塘地块和祁连地块内,各有8个预测窗口出现过至少1个热点。在2009—2020年预测窗口中,5个滑动窗内巴颜喀拉地块和柴达木地块交界处青海玛多县附近持续存在密集热点(图4)。在2012—2019年、2013—2020年预测窗口中,祁连地块出现了明显热点,主要分布在青海海西州德令哈市附近。在2014—2021年、2015—2022年、2016—2023年3个预测窗口时段内(图5a～c),热点主要在南北向上跨巴颜喀拉地块、柴达木地块和祁连地块3个二级块体。

在2014—2021年预测窗口中(图5a)无零散热点,巴颜喀拉地块内玛多地震震中所在网格热点颜色最深,发震危险性最高。从2015—2022年预测窗口开始至2021—2028年窗口(图5b～h),在羌塘地块内始终存在热点。在2016—2023年预测窗口内(图5c),热点总数变少,羌塘地块以北的3个地块内热点数皆有减少。玛多地震东侧巴颜喀拉地块和柴达木地块交界处热点变少但颜色加深,危险性增强。至2017—2024年、2018—2025年预测窗口(图5d、e),玛多地震震中及其摩尔邻近网格不存在有效热点。在2017—2024年预测窗口中,柴达木地块和祁连地块的最北端热点延伸到了甘肃境内。至2018—2025年预测窗口中(图5e),巴颜喀拉地块以北的热点全部消失,浅色热点集中于巴颜喀拉地块和羌塘地块交汇处(青海南部),且在巴颜喀拉地块东北缘(四川北部甘川交界附近)出现了一个深色热点。

在最后3个连续的预测窗口中(图5f～h),玛多地震东侧巴颜喀拉地块(北部)和柴达木地块的交汇处存在密集深色热点,且热点个数逐渐减少。2019—2026年预测窗口(图5f),相对于前2个滑动窗,玛多地震震中及摩尔邻近网格出现深色密集热点,发震危险性显著增强。2020—2027年预测窗口(图5g),相对于上一个窗口热点总个数不变,但零散热点危险性降低。巴颜喀拉地块和柴达木地块的交汇处的密集热点危险性增强,玛多地震震中热点消失。至2021—2028年预测窗口(图5h),巴颜喀拉地块和柴达木地块的交汇处密集热点进一步收缩,体现了后续发生强震的危险性在向玛多地震的东侧迁移,并且可能和震源区及其附近仍然主要受到NEE向近水平挤压应力场的控制有关(吕苗苗等,2022)。

青藏高原具有持续变形的特点(许志琴等,2006),热点的迁移也是相对变化的,本文研究中没有任意2个预测窗口内热点的分布及发震概率是完全相同的。从2014—2021年预测窗口开始至2021—2028年预测窗口,玛多地震的有效热点的个数、分布位置和发震危险性并不具有稳定性、持续性,这体现了7级以上强震孕震物理过程的复杂性。研究区内所有热点在不同次级块体内的迁移也不具有明确的规律性,其内部地震的发生和区域构造应力背景与次级活动地块息息相关。徐志国等(2021)通过应力场反演结果得出,研究区整体的背景应力场具有近EW向挤压应力场特征,与巴颜喀拉块体的水平最大主应力方向总体上为EW向(范桃园等,2013)具有较好的一致性。区域构造应力场能够影响震源区的变形和破裂,本次玛多M_S7.4地震为双侧破裂(李智敏等,2021; 王未来等,2021),但图5不同预测窗口中玛多地震的西侧并没有显著热点分布。中国大陆强震孕育发生的动力学环境整体呈“纵向分层、横向分块”的特征(邵志刚等,2022),块体之间也存在一定的相互作用,而块体运动受深部物质作用驱动。大型走滑型地震发生时,复杂的平移错动引起的浅部破裂面和深部破裂面是存在差异的(徐志国等,2021),其深部物质运移和深部地下结构异常会对后续地震产生影响(许志琴等,2006; 王椿镛等,2008; 嘉世旭等,2017; 宋向辉等,2021)。用PI方法计算得到的热点对局部区域某个时段的发震范围具有参考作用,但不能精准预测发震地点和震源深度。研究区内的热点分布与迁移和浅层破裂展布并不完全对应,可能和强震孕育过程、深部动力、地下深部结构等因素存在一定关系。