“十五”和“十一五”期间,由中国地震局牵头的国家重大科学工程项目“中国大陆构造环境监测网络工程”(简称陆态网络)启动,项目在全国范围内布设了高密度、大范围的GNSS观测站,为GNSS观测资料应用于地壳形变分析、地球动力学研究及地震预测预报提供了基础观测数据支撑。近年来,随着“云南省政府十项重点工程”“川滇地震预报实验场”以及“北斗地基增强”等项目的相继实施,观测站点数量大幅提升,本文围绕云南地区的65 个GNSS 连续观测站点、200多个流动观测站点多年积累的观测资料为基础开展研究(图1)。

图1 云南及邻区构造背景和GNSS观测站分布
Fig.1 Active faults and GNSS stations in Yunnan and its adjacent areas

本文采用GAMIT/GLOBK软件对数据进行处理。选取BJFS、SUWN、USUD、WUHN、LHAZ、PIMO等十几个中国及周边地区的IGS跟踪站与研究区域内GNSS站点进行联合解算。以每24 h作为一个分析时段,利用GAMIT软件获取包含测站位置、大气参数、地球定向参数以及相应的方差-协方差矩阵的单日松弛解; 利用GLOBK卡尔曼滤波程序将GAMIT软件处理得到的单日松弛解和IGS全球单天无约束解(igs1-igs9)的联合平差,选取全球范围内均匀分布、时间序列稳定的IGS连续观测站作为框架约束,获取在该框架下待解算测站的点位结果。综合几年来的单日解,可得到ITRF2014参考框架下点位时间序列,通过对不同时间段时间序列线性拟合获得在ITRF2014参考框架下站点的速度场。为排除后期应变计算时误将同震位移产生的应变场扰动当成异常,本文对数据做了同震阶跃的处理,具体方法是:在GAMIT数据处理的过程中,在eq文件中定义地震,根据震级决定震中距离,对指定震中距中所有的测站采用tsfit程序进行参数估计。考虑到地震后调整所引起的应变场变化可能对其他地区产生应力加载而进一步诱发其他地震,所以对于震后的黏弹性松弛影响并未进行处理。此外,在处理过程中对数据观测质量较差的历元(或者测站)以及最终速率中误差大于2 mm/a的测站进行剔除。最终得到1999—2007、2009—2014、2015—2020年3个时段云南地区相对于稳定欧亚板块的运动速度场(表2),后期的应变分析均以此为约束。

图2 云南地区相对于欧亚板块的运动速度场
Fig.2 Velocity field in Yunnan relative to the Eurasian plate

选择适当的数学模型通过速度场求取区域应变特征参数是准确理解区域地壳形变特征的关键,为避免点位密度差异带来的影响,引入克里金插值方法进行应变率场估计,该方法理论基础严密,是一种线性、无偏、最优的内插估计算法(陈小斌,2007; 朱守彪等,2005)。本文在获得了各期扣除了大震同震影响速度场的基础上,将速度场转化为区间位移场,由于云南地区站点分布密度平均在30 km,为了客观反映每个点位变化对应变场的影响,采用0.2°×0.2°为单元进行网格化,每个网格分辨率大约为20 km,基本和点位密度相匹配。实际计算过程中,通过GLOBK平差得到的坐标N、E分量计算网格中心点经纬度信息,通过大地反算法计算任意两个网格之间的大地线长度以及大地线长度的变化量,并计算网格中心点之间的大地方位角。

假设某个测点的位移为u、v,其应变状态分量为ε_x、ε_y、γ_xy,其中ε_xy= 1/2 γ_xy,d_x、d_y为两点间距离分量的变化量,ω为旋转量,α 为两格网点间的坐标方位角。则与其无限接近的一点的位移分量可表示为:

式(1)两边同时除两点间距离,可转变为线应变与方位角之间的关系为:

每个格网与其他各个相邻格网点通过(2)联立方程组,再通过最小二乘法求解得到应变状态分量,进一步可以计算得到其他的应变参数,包括最大剪应变、面应变等(Hong et al,2018)。

中强地震的发生需经历长期应变积累并逐渐趋于极限状态,因此,探讨应力应变积累背景与动态分布特征,可为地震危险地点的判定提供一定依据。 单期GNSS 应变率场结果可辨识不同区域变形强弱,多期GNSS应变率场结果则能进一步反映同一区域构造变形增强或减弱的时空演化特征,进而跟踪应力应变增强区与弱化区,可了解地震发生前震中附近的应变状态变化情况。因此,本文利用上述方法获取1999—2007、2009—2014、2015—2020年云南地区的面应变率和最大剪应变率场,为分析地震事件与区域内地壳活动的内在联系,统计了研究区域内(21°～30°N,97°～106°E)3个观测时段之后3年内M_S≥5.0地震事件信息(表1)。

根据区域地壳面应变率和最大剪应变率的空间变化以及相应时段之后的3年内M_S≥5.0地震事件分布特征,绘制不同时期云南地区GNSS应变

率场及相应时段之后3年内M_S≥5.0地震事件分布图(图3),其中面应变率场图中约定面压缩为负、面膨胀为正。图4给出了不同时期云南地区应变率场残差标量分布,从图中可以看到,除少量网格外,绝大多数网格的应变率场的拟合残差优势分布在1.5×10^-8/a范围内。

通过回溯性研究显著地震前后应变率场的空间的动态变化,探讨应变积累背景异常特征对未来强震地点的指示意义。从3个不同时期的应变率场来看,云南地区地质构造和地震孕育活动错综复杂,面应变率交替出现压缩和膨胀的特征,最大剪应变率高值区位于川滇菱形块体的整个东边界,沿着鲜水河—安宁河—则木河—小江断裂带展布,与该断裂带的走向基本一致,西边界与丽江—小金河断裂带交会,滇西南腾冲、普洱等高值区域。上述结果与相关研究结果总体特征较为一致(田晓等,2021),如小江断裂带和滇西北东条带整体为最大剪应变率高值区,以及滇西北区域整体为拉张等特征均表现出了宏观上的一致性。对比不同时期面应变率及最大剪应变率不难看出空间分布的相似性,说明该区域的构造及相对运动的背景还是保持一定的稳定性和继承性,但也存在一些高值区的细微差异,异于长期背景的应变能高值异常区域,理论上后期更容易触发地震。

表1 云南地区 MS≥5.0目标地震统计(2011—2022年)
Tab.1 Statistics of MS≥5.0 earthquakes in Yunnan from 2011 to 2022

图3 不同时期云南地区GNSS应变率场及之后3年内MS≥5.0地震分布
Fig.3 Strain rate of GNSS in different stages from 2009 to 2020 and the MS≥5.0 earthquakes in the following three years

图4 不同时期云南地区GNSS应变率残差标量分布
Fig.4 Distribution of the residuals scalar of the GNSS strain rate in Yunnan area in different stages from 1999 to 2020

通过统计不同时期的应变率场与其之后3年内云南地区M_S≥5.0地震事件发震地点发现:绝大多数中强地震都发生在面应变高梯度带的张压转换区和最大剪应变高梯度区或其边缘。2008年8月30日攀枝花M_S6.1地震位于1999—2007年最大剪应变高梯度区和面应变张压转换区,随着时间的推移,2009—2014年整个小江断裂带最大剪应变强度逐渐减弱,高值区域范围逐渐缩小,该发震区域剪应变梯度不太明显(图3b)。2009—2014年在景谷—思茅一带形成了面拉张与压缩高值区交替变化异常(图3e),之后发生了2014年景谷M_S6.6地震,震中位置正好位于面应变正负交替异常区。从2015—2020年应变率场来看(图3c、f),2022年9月5日沪定M_S6.8地震也位于最大剪应变高梯度区和面应变正负交替异常区。

王伶俐等(2020)计算了1999—2004、2004—2007、2009—2013、2013—2015、2015—2018年5个时段的区域应变率场,发现相应时段之后的3年内M_S≥5.0地震事件震中位置也有类似特征,说明GNSS应变率场对其之后1～3 a内的中强震发生区域有一定的指示意义。该现象也与江在森等(2012,2009,2020)多位学者对于强震危险区的认识相符。可见应变率异常区的判定对于中长期地震危险区识别有一定指示意义。

应变率异常区与历史地震事件有较好的对应关系,研究区域各个观测时段GNSS应变率场异常区对其之后3年内的中强震发生地点有一定的指示意义。因此,以面应变高梯度带的张压转换区,以及最大剪应变率高值区域及边缘地带来判定未来地震危险区是可行的。但是如何确定危险区范围具有较强的主观因素,因此,有必要进一步以数学模型替代人为的主观判断来识别异常区的危险性。本文提出了数值计算模型:通过估计格网最大剪应变和面应变风险区划因子,定量提取异常区危险指标,该模型思路如下:

(1)最大剪应变风险区划因子的处理

最大剪应变最为活跃的区域为断层区域。如果对最大剪应变进行整体分析,由于断层区域和非断层区域最大剪应变量差异较大,平行断层走向会形成最大剪应变高梯度带。当识别最大剪应变梯度带时,得到的结果将可预见的平行于走滑断层走向展布,而这些区域并非所期望的最大剪应变梯度区。相关震例研究结果显示,沿断层走向的张压转换形变梯度带和最大剪应变梯度带,是需要识别的潜在的地震危险区(朱治国等,2018),因此,本文以断层为单元进行危险性分析,首先判定网格与断层的关系,假设某一条断层穿过了n个网格,任意网格最大剪应变率为R_i(i=1,…n),其中i为断层穿过网格的顺序号,非网格编号。定义每个网格的最大剪应变率沿断层方向的梯度极大值为：

最大剪应变沿着断层方向变化越快(梯度越大)的地方越危险,因此,将该值作为表征区域危险性的基础值。在此基础上,计算穿过同一断层所有网格的最大剪应变率的均值(Avgg)、标准差(Std)以及最大剪应变梯度的均值(Avgg),设断层穿过的网格地震危险性为Rd_i(i=1,…,n),令:

如果某个网格的应变梯度超出了所有网格梯度均值时,认为该网格应变梯度偏高,存在危险性,反之认为无危险。在梯度偏高的情况下,[HJ2.3mm]如果其最大剪应变率小于断层上所有网格最大剪应变率的一倍标准差时,认为该段内断层剪切速率衰减快,且断层存在弱剪切现象,有闭锁特征,其危险倍率加倍。

(2)面应变风险区划因子的处理

云南地区以走滑断层分布为主,面应变的梯度区并不像最大剪应变一样沿断层有显著的相关性,也难以像逆冲断层一样以单一断层为研究目标,可通过分析断层不同段的应变差异来识别危险区(赵静等,2019)。因此,首先识别所有的面应变张压转换带,作为潜在的地震风险源。为了更好地扫描出面应变的张压转换带,以一定的空间尺度为单元,对每个网格的面应变率场开展扫描分析。设计算得到的网格的面应变率结果为E_i(i=1,…,n),i为网格编号。对于任意网格i,寻找其周边100 km(尺度可自定义,定义半径50 km,对单网格来说扫描空间区域为100 km)范围内的网格有m个,将这些网格按照应变率的正负进行分类,假设应变率为张性的网格有p个,其应变率表示为E_ig(g=1,…,p),应变率的均值表示为Ag; 应变率为压性的网格有q个,其应变率表示为E_is(s=1,…,q),应变率的均值表示为As。定义

式中:E_Dif用于表征扫描区内应变的差异程度; W_Dif为权重,反映应变差异的显著程度,取值0～1,如果正负应变的网格数量差异不大时,该变量趋近于1,表明该网格处于张压转换带; E' _Dif 最终反映了该扫描区的应变差异程度。完成所有网格的扫描后,计算E' _Dif系列的均值,表示为Ad,则任意网格面应变的危险性(Ed_i)为:

从式(8)可见,当网格的应变梯度小于均值时,表示该区域应变梯度不显著,危险性不大; 当应变梯度大于均值时,将应变梯度作为危险因子赋值给该网格。

基于上述模型,计算得到2009—2014年最大剪应变率与面应变率场异常特征区域(图5),并从GlobalCMT下载了对应统计时间范围内的有关地震(表1)的震源机制解。依据面应变率和最大剪应变率异常区提出的地区危险区识别模型更清楚地显示了发震风险区域。

从图5可以看到,通过地震危险区模型可以基本提取出面应变高梯度带的张压转换区及最大剪应变率高值区域及边缘地带。图5a显示的剪应变率异常区在未来3年内较大概率对应发生走滑地震,如2016年云龙M_S5.0地震,2017年漾濞M_S5.1地震,但也存在非异常区域发生了地震的情况,如2014年景谷M_S6.6地震。有研究表明景谷地震震中无明显与已知断层相关的地表破裂,初步认为其发震构造为一条新生的兼具正断性质的右旋走滑断层(谢张迪,韩竹军,2019)。剪应变率异常区域对该地震对应不是很理想,究其原因,可能是数值模拟过程中,所采用的断层模型不够全面,无法对新生发震构造起到约束,因此,在后期可能还需要进一步完善断层模型。但是漾濞地震在面应变异常区得到了较好的对应效果。图5b显示面应变异常区域也较大概率发生了具有正断或逆冲分量的地震。因此,在实际震情跟踪应用过程中,不仅要关注面应变高梯度张压转化帯,还要同时关注剪切变形弱化区,可以将每个格网的面应变率与最大剪应变危险因子值进行求和计算,综合判识可能的异常区。即对于任意网格,其地震风险定义为:

图5 2009—2014年应变率场异常区模型识别结果
Fig.5 The result from the automatic identification of the abnormal region of the strain rate field from 2009 to 2014 by the model

式中:W_i为权重,由应变计算误差来计算。具体计算方法是:假设δ为第i个网格计算得到的应变中误差,其应变量为θ,则将1-(δ/θ)作为该网格危险因子的权重,从而实现将相对中误差作为权重参与模型计算,当相对中误差越小时,权重越大,可得到更可靠的地震风险性量化因子。鉴于风险因子是无单位的量值,对其结果进行归一化处理,得到3个时段的区域地震风险性综合评估结果(图6)。

图6 1999—2020年区域地震风险性综合评估结果
Fig.6 Results from the synthetic estimation of the earthquake risk in Yunnan and its adjacent areas from 1999 to 2020

从图6可看出,模型结果匹配较好。划定区域与预测区间内大多数地震具有较好的对应关系,该模型一定程度上丰富了中长期地震预测工作中发震位置的判定方法,当然,实验过程中也存在部分异常区域未发生较大震级地震或未异常区域发生了地震的情况,这可能与模型内阈值参数设置有关,因此模型参数优化问题后续值得进行更深入研究。

前期的试验结果显示采用地震危险区模型识别出的异常区与地震事件具有较好地对应关系,表明该方法判定依据应用于地震预测预报具有可行性。不过在实际的震情跟踪分析和会商研判过程中对方法的好坏的评估还缺乏客观的评价依据,异常区分析结果以及相关指标的信度高低无法区分。因此本文在危险区模型识别的基础上,对该方法以多期速度场结果进行地震预测效能量化评估与分析。

目前地震预测研究中最常用的一种评价预报效能的方法是由许绍燮(1989)提出的R值评分法,马宏生等(2004)、吴中海和赵根模(2013)对地震预报的 R值评分方法进行了阐述和发展。近些年来最常用的时空 R值评分公式为:

式中:c是报准率; b是预报时空占有率; R值的计算中,减去了预报的时空占有率,即扣除了随机概率的预报成功率,当给定一个置信水平(一般取97.5%)(国家地震局科技监测司,1990),按照二项式分布原则编制最低R值对应表(R₀值),当 R>R₀时,通常认为预报效能R已经达到97.5%置信度,即有预报意义。

针对异常风险区划指标的效能评估,只需考虑预报的空间占有率。基于前期计算得到的危险性区划预测因子,与地震结合起来,通过R值评分方法进行预报效能的评价,为了获得最优的预报效果, 求取所有网格地震危险因子(M_di)的均值()和标准差(σ_M),以均值为最低阈 值线,阈值动态增加,每次增加步长为0.1×σ_M,阈值上限为,动态扫描30个不同阈值,求取地震危险因子的均值和标准差,并按照标准差的倍率进行阈值动态扫描,危险因子超过阈值的区域认为存在危险,这样获得的危险区域和地震进行对应,通过空间占有率进行R值计算,计算得到的R值取最大值为最优阈值结果。

在指标提取中,地震样本的选取见表1。因此,本文提取的指标也主要对云南地区M_S≥5.0地震区的划分有指示意义。各时期应变率场异常地震风险预测指标预报效能评估结果见表2。表中准确率和漏报率是以参与统计的地震样本数为分母来统计,虚报率是以数值模型中实际参与评估的1 437个格网(排除了境外无实际测点的网格)为统计样本。

表2 预测指标预报效能评估结果
Tab.2 Results from the evaluation of the mid-term prediction performance

表2中针对3期应变率场异常风险区域地震预测效能量化评估结果R值均高于R₀,说明该预测结果通过R值评分检测,预测指标有效,证明危险区模型识别方法是有效和可行的。虽然该模型丰富了地震发震位置的判定方法,但也存在部分异常区域漏报率相对较高的形况,究其原因为:一方面是与点位的分布疏密程度及点位反映构造活动能力差异有关。时间较长的1999—2007年应变率场映震结果相对较好,可能也表明区域速度场的稳定性对于映震效果有着显著影响; 另一方面可能与模型内阈值参数设置有关,后期阈值参数的扩展调整及评估因素的增加可能使获取指标相对更优、评估更加合理。