云南地区GNSS应变率场异常识别方法及地震预测效能评估

摘要：利用1999—2020年云南及邻区近300个GNSS测站的观测数据解算获取的速度场为约束，采用克里金插值方法分时段估计了1999—2007年，2009—2014年，2015—2020年三期区域应变率场；通过回溯各个观测时段之后3年内MS≥5.0地震事件，分析区域地壳形变特征与地震事件发生地点之间的相关性，结果表明，绝大部分地震都发生在面应变高梯度带的张压转换区和最大剪应变率沿断层方向的高值区或其边缘。基于上述应变率场异常特征，提出格网地震危险因子算法，建立地震危险区识别模型，通过估计格网最大剪应变率和面应变率风险区划因子，定量提取异常区地震危险指标，结果显示采用数值模型识别出的异常区与地震事件具有较好地对应关系；进一步采用R值评分的方式对应变率场异常区模型识别方法进行地震预测效能量化评估与分析，结果显示3期应变率场预测结果均通过R值评分检测。

关键词：GNSS；应变率场；地震预测；R值评分；异常指标；云南地区

中图分类号：P315.725"" 文献标识码：A"" 文章编号：1000-0666（2025）01-0049-12

doi：10.20015/j.cnki.ISSN1000-0666.2025.0006

0 引言

地震孕育的能量主要来源于地壳相对运动与变形造成的弹性应变能积累，因此，自从中国大陆GNSS观测站获取地壳形变信息产出后，便被作为认知地震孕育过程最直接有效的途径之一，被广泛应用于实际地震预测工作中。研究人员基于中国大陆地壳应变率场与构造动力背景关联的基本格局、分区特征等（江在森等，2003；朱守彪等，2005；邓起东等，2014），分析研究GNSS应变率场动态特征，获取的动态信息用于强震活动主体区带的分析。Wu等（2022）统计了中国大陆不同时间窗口中GNSS应变率与地震活动的动态相关性，发现中国大陆高应变率区域与1900年以来的中强地震震中位置相当吻合。部分学者（武艳强等，2013；徐克科，李伟，2017；魏文薪等，2018；魏斌，陈长云，2024）对一些典型地震震前或者同震的应变率场影响进行研究。

近年来，随着GNSS观测站点的大量布设和震例资料的不断累积，地震前兆异常区的相关研究获得越来越多的关注，不少学者利用GNSS资料监测川滇地区地壳形变动态变化特征（党亚民等，2019），基于GNSS应变时序对云南地区构造运动与地震事件孕震模式进行分析，探讨GNSS应变场时空演化特征与地震事件的关联（洪敏等，2014；高涵等，2018；王伶俐等，2020），为探索地壳运动规律和提取与地震相关的前兆异常信息提供了思路。但以往的研究过程中，多以震例与现象之间关系分析为主，危险区的判定多依赖个人经验，对应变率场异常区的判定缺乏定量化识别方法，更缺乏对识别效果的系统评价。在会商工作逐步规范化的进程中，由于缺乏异常判定标准及定量化的预测指标提取规则，导致应变场分析结果在地震会商中应用受限。

本文通过研究云南地区地壳形变时空演化特征，探讨区域地壳形变特征与地震事件发生地点之间的相关性，并尝试通过数学模型自动分析应变率场异常特征，对风险区域进行自动识别，避免人为主观判断的影响。在此基础上，利用R值评分，提取其具有地震前兆指示意义的地震预测预报指标，以期为后续使用该方法进行地震危险地点的判定提供参考。

1 数据解算

1.1 数据资料收集

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

1.2 数据解算

本文采用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 云南地区GNSS应变场与地震事件关联分析

2.1 区域应变率场计算方法

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

假设某个测点的位移为u、v，其应变状态分量为εx、εy、γxy，其中εxy=12γxy，dx、dy为两点间距离分量的变化量，ω为旋转量，α 为两格网点间的坐标方位角。则与其无限接近的一点的位移分量可表示为：

u′=u+εxdx+εxydy-ωdyv′=v+εxydx+εydy+ωdx（1）

式（1）两边同时除两点间距离，可转变为线应变与方位角之间的关系为：

Δud=εxcosα+εxysinα-ωsinαΔvd=εxycosα+εysinα+ωcosα（2）

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

2.2 应变率场异常与地震关联特征分析

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

根据区域地壳面应变率和最大剪应变率的空间变化以及相应时段之后的3年内MS≥5.0地震事件分布特征，绘制不同时期云南地区GNSS应变率场及相应时段之后3年内MS≥5.0地震事件分布图（图3），其中面应变率场图中约定面压缩为负、面膨胀为正。图4给出了不同时期云南地区应变率场残差标量分布，从图中可以看到，除少量网格外，绝大多数网格的应变率场的拟合残差优势分布在1.5×10-8/a范围内。

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

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

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

3 应变率场异常地震风险区划指标提取

3.1 危险区识别模型

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

（1）最大剪应变风险区划因子的处理

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

Rgi=maxk（Ri-Rk） （k=i-1，i+1）（3）

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

Rdi=0（Rgilt;Avgg）

Rdi=Rgi×2 （Rilt;Avg-Std且Rgigt;Avgg）

Rdi=Rgi ""（Rigt;Avg-Std且Rgigt;Avgg）（4）

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

（2）面应变风险区划因子的处理

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

EDif=Ag-As（5）

WDif=1-p-qp+q（6）

E′Dif=EDif×WDif（7）

式中：EDif用于表征扫描区内应变的差异程度；WDif为权重，反映应变差异的显著程度，取值0～1，如果正负应变的网格数量差异不大时，该变量趋近于1，表明该网格处于张压转换带；E′Dif最终反映了该扫描区的应变差异程度。完成所有网格的扫描后，计算E′Dif系列的均值，表示为Ad，则任意网格面应变的危险性（Edi）为：

Edi=E′Dif （E′Dif gt;Ad）

Edi=0 （E′Dif≤Ad）（8）

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

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

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

Mdi=（EDi+RDi）×Wimaxi=1…n（Mdi）（9）

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

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

3.2 应变率场异常风险区划指标预报效能评估

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

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

R=c-b=报准地震数实际地震数-预报占用时间和空间预报研究总时间和空间（10）

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

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

R=报准地震数实际地震数-预测有地震的网格数网格总数（11）

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

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

4 结论

本文利用云南地区GNSS速度场为约束，采用克里金插值方法分时段估计了1999—2007年、2009—2014年、2015—2020年区域应变率场；通过对各个观测时段之后3年内MS≥5.0地震事件进行回溯分析，归纳出识别强震危险地点的应变率场异常判据，并建立地震危险区识别数学模型对应变率场异常区域进行自动识别，最后采用R值评分对模型识别结果进行地震预测效能评估，得出以下结论：

（1）云南地区MS≥5.0地震大多发生于前期GNSS观测到的面应变高梯度带的张压转换区及最大剪应变率沿断层方向的高值区域及边缘地带，可见研究区各个观测时段GNSS应变率场对其之后3年内的中强震发生区域有一定的指示意义。

（2）从地震危险区识别模型计算结果提取的异常危险区来看，划定的危险区与预测区间内大多数地震具有较好的对应关系，模型结果匹配较好。该模型一定程度上丰富了中长期地震预测工作中发震位置的判定方法。

（3）3个时期应变率场预测结果均通过R值评分检测，预测指标有效，证明地震危险区模型识别方法信度较高。但本文方法获得的预测指标均基于概率统计，该类指标大多会受地震样本数多少的影响。后续还要围绕干扰排除、震例时空相关性等进行更深入研究，对地震事件进行不断的检验和优化，对模型内阈值参数设置不断测试和调整，才能为后期震情跟踪工作提供更合理的定量化指标。

本文在GPS区域网观测资料收集处理过程中得到湖北省地震局赵斌研究员、王东振工程师的帮助，在此表示衷心感谢。

参考文献：

陈小斌.2007.中国陆地现今水平形变状况及其驱动机制［J］.中国科学：地球科学，37（8）：1056-1064.Chen X B.2007.Present-day horizontal deformation status of continental China and its driving mechanism［J］.Scientia Sinica Terra，37（8）：1056-1064.（in Chinese）

党亚民，杨强，薛树强，等.2019.GNSS监测的川滇地区地壳形变动态变化特征［J］.大地测量与地球动力学，39（2）：111-116.Dang Y M，Yang Q，Xue S Q，et al.2019.GNSS monitoring dynamic variation characteristics of crustal deformation in the Sichuan-Yunnan region［J］.Journal of Geodesy and Geodynamics，39（2）：111-116.（in Chinese）

邓起东，程绍平，马冀，等.2014.青藏高原地震活动特征及当前地震活动形势［J］.地球物理学报，57（7）：2025-2042.Deng Q D，Cheng S P，Ma J，et al.2014.Seismic activities and earthquake potential in the Tibetan plateau［J］.Chinese Journal of Geophysics，57（7）：2025-2042.（in Chinese）

高涵，张明，秦姗兰，等.2018.基于多尺度球面小波分析云南地区应变特征［J］.大地测量与地球动力学，38（12）：1232-1237.Gao H，Zhang M，Qing S L，et al.2018.Analysis the strain characteristics of Yunnan：based on multi-scale spherical wavelet［J］.Journal of Geodesy and Geodynamics，38（12）：1232-1237.（in Chinese）

国家地震局科技监测司.1990.地震学分析预报方法程式指南［M］.北京：地震出版社.Department of Science and Technology Monitoring，National Earthquake Administration.1990.Program guide for seismological analysis and prediction methods［M］.Beijing：Seismological Press.（in Chinese）

洪敏，张勇，邵德盛，等.2014.云南地区近期地壳活动特征［J］.地震研究，37（3）：367-372.Hong M，Zhang Y，Shao D S，et al.2014.Recent tectonic activity features of Yunnan region［J］.Journal of Seismological Research，37（3）：367-372.（in Chinese）

江在森，方颖，武艳强，等.2009.汶川8.0级地震前区域地壳运动与变形动态过程［J］.地球物理学报，52（2）：505-518.Jiang Z S，Fang Y，Wu Y Q，et al.2009.The dynamic process of regional crustal movement and deformation before Wenchuan MS8.0 earthquake［J］.Chinese Journal of Geophysics，52（2）：505-518.（in Chinese）

江在森，马宗晋，张希，等.2003.GPS初步结果揭示的中国大陆水平应变场与构造变形［J］.地球物理学报，46（3）：352-358.Jiang Z S，Ma Z J，Zhang X，et al.2003.Horizontal strain field and tectonic deformation of China mainland revealed by preliminary GPS result［J］.Chinese Journal of Geophysics，46（3）：352-358.（in Chinese）

江在森，武艳强，邹镇宇，等.2020.GNSS在中国大陆的地震预测应用研究进展与展望［J］.中国地震，36（4）：693-707.Jiang Z S，Wu Y Q，Zou Z Y，et al.2020.Progress and prospect of earthquake forecasting based on GNSS data in China［J］.Earthquake Research in China，36（4）：693-707.（in Chinese）

江在森，武艳强.2012.地壳形变与强震地点预测问题与认识［J］.地震，32（2）：8-21.Jiang Z S，Wu Y Q.2012.Crustal deformation and location forecast of strong earthquakes：understandings and questions［J］.Earthquake，32（2）：8-21.（in Chinese）

马宏生，刘杰，吴昊，等.2004.基于R值评分的年度地震预报能力评价［J］.地震，4（2）：31-37.Ma H S，Liu J，Wu H，et al.2004.Scientific evaluation of annual earthquake prediction efficiency based on R-value［J］.Earthquake，4（2）：31-37.（in Chinese）

田晓，占伟，郑洪艳，等.2021.川滇地区现今三维地壳运动特征［J］.大地测量与地球动力学，41（7）：739-746.Tian X，Zhan W，Zheng H Y，et al.2021.Characteristics of present-day 3D crustal movement of Sichuan-Yunnan region［J］.Journal of Geodesy and Geodynamics，41（7）：739-746.（in Chinese）

王伶俐，洪敏，张勇，等.2020.云南地区GNSS应变场时空演化与地震事件关联分析［J］.中国地震，36（1）：91-104.Wang L L，Hong M，Zhang Y，et al.2020.Analysis of the relationship between the spatiotemporal evolution of GNSS strain field and earthquake events in Yunnan Region［J］.Earthquake Research in China，36（1）：91-104.（in Chinese）

魏斌，陈长云.2024新疆及邻区现今GNSS变形特征与地震关系研究［J］.地震研究，47（3）：419-429.Wei B，Chen C Y.2024.Study on the Relationship between GNSS deformation characteristics and earthquakes in Xinjiang and its adjacent regions［J］.Journal of Seismological Research，47（3）：419-429.（in Chinese）

魏文薪，江在森，邵德胜，等.2018.2014年鲁甸6.5级地震GPS同震位移及反演分析［J］.地球物理学报，61（4）：1258-1265.Wei W X，Jiang Z S，Shao D S，et al.2018.Coseismic displacements from GPS and inversion analysis for the 2014 Ludian 6.5 earthquake［J］.Chinese Journal of Geophysics，61（4）：1258-1265.（in Chinese）

吴中海，赵根模.2013.地震预报现状及相关问题综述［J］.地质通报，32（10）：1493-1512.Wu Z H，Zhao G M.2013.The earthquake prediction status and related problems：A review.［J］.Geological Bulletin of China，32（10）：1493-1512.（in Chinese）

武艳强，江在森，王敏，等.2013.GNSS监测的芦山7.0级地震前应变积累及同震位移场初步结果［J］.科学通报，58（20）：1910-1916.Wu Y Q，Jiang Z S，Wang M，et al.2013.Preliminary results of the co-seismic displacement and pre-seismic strain accumulation of the Lushan MS7.0 earthquake reflected by the GPS surveying.［J］.China Sci Bull，58（20）：1910-1916.（in Chinese）

谢张迪，韩竹军.2019.2014年云南景谷MS6.6地震发震断层及其动力学参数［J］.地震地质，41（4）：887-912.Xie Z D，Han Z J.2019.Study on the seismogenic fault and dynamics parameters of the 2014 MS6.6 Jinggu Earthquake in Yunnan［J］.Seismology and Geology，41（4）：887-912.（in Chinese）

徐克科，李伟.2017.利用GNSS基线分析芦山MS7.0级地震前后应变演变特征［J］.武汉大学学报·信息科学版，42（8）：1054-1060. Xu K K，Li W.2017.Strain evolution characteristics before and after Lushan MS7.0 earthquake using GNSS baseline［J］.Geomatics and Information Science of Wuhan University，42（8）：1054-1060.（in Chinese）

许绍燮.1989.地震预报能力评分［C］//地震预报方法实用化研究文集：地震学专辑.北京：学术书刊出版社，586-590.Xu S X.1989.Earthquake prediction ability score.Seismic prediction method practical research collection：Seismological Album［C］.Beijing：Academic Press，586-589.（in Chinese）

赵静，刘杰，任金卫，等.2019.汶川、芦山地震前龙门山断裂带地壳形变特征对比分析［J］.地球学报，40（1）：186-198.Zhao J，Liu J，Ren J W，et al.2019.A contrastive analysis of crustal deformation characteristics along the Longmen Shan fault zone before the MS8.0 Wenchuan earthquake and the MS7.0 Lushan earthquake［J］.Acta Geoscientica Sinica，40（1）：186-198.（in Chinese）

朱守彪，蔡永恩，石耀霖.2005.青藏高原及邻区现今地应变率场的计算及其结果的地球动力学意义［J］.地球物理学报，48（5）：1053-1061.Zhu S B，Cai Y E，Shi Y L，et al.2005.Computation of the present-day strain rate field of the Qinghai Tibetan plateau and its geodynamic implications［J］.Chinese Journal of Geophysics，48（5）：1053-1061.（in Chinese）

朱治国，秦姗兰，艾力夏提·玉山，等.2018.基于GPS的精河MS6.6地震前地壳变形动态特征研究［J］.中国地震，34（3）：428-434. Zhu Z G，Qin S L，Ailixiati Y S，et al.2018.Study on the dynamic characteristics of crustal deformation before the Jinghe MS6.6 earthquake based on the GPS［J］.Earthquake Research in China，34（3）：428-434.（in Chinese）

Hong M，Shao D S，Wu T F，et al.2018.Short-impending earthquake anomaly index extraction of GNSS continuous observation data in Yunnan，Southwestern China［J］.Journal of Earth Science，29（1）：230-236.

Wu Y Q，Jiang Z S，Pang Y，et al.2022.Statistical correlation of seismicity and geodetic strain rate in the Chinese mainland［J］.Seismological Research Letters，93（1）：268-276.

Identification Method of GNSS Strain Rate Field Anomaly and Evaluationof Earthquake Prediction Efficiency in Yunnan

WANG Lingli1，3，HONG Min1，NIU Tian1，LI Qin2，3，YANG Xinjun1，Yu Shixian1

（1.Yunnan Earthquake Agency，Kunming 650201，Yunnan，China）

（2.Yunnan Institute of Geo-environment Monitoring，Kunming 650216，Yunnan，China）

（3.Yunnan Key Laboratory of Geohazard Forecast and Geo-ecological Restorationin Plateau Mountainous Area，Kunming 650216，Yunnan，China）

Abstract

In this paper，the velocity field is calculated by using the observation data from nearly 300 GNSS stations in Yunnan and its adjacent areas in the period from 1999 to 2020.Based on the velocity field，three stages（1999-2007，2009-2014，and 2015-2020）of the regional strain rate fields are estimated by Kriging interpolation method.By backtracking the earthquake events（MS≥5.0）in the following 3 years of each observation stage，the correlation between the characteristics of the regional crustal deformation and the location of the earthquake events is analyzed，and some abnormal criteria of strain rate field identifying the risk locations of the strong earthquakes are summarized.That is，most of the earthquake cases occur in the tension-pressure transition zone of the high gradient zone of surface strain，and the high value zone or the edge of the maximum shear strain along the fault.Based on the abnormal characteristics of strain rate field，this paper proposes a grid algorithm of the seismic risk factor，and establishes a model of the seismic-risk area identification.By estimating the recognition factors of the risk zones of the maximum shear strain and the surface strain with the grid algorithm，the seismic risk indicators of the abnormal areas are quantitatively extracted.The abnormal areas identified by the numerical model corresponds well with the seismic events.In order to provide an objective evaluation basis for the subsequent use of this model to track earthquake conditions and to determine earthquake risk locations，this paper adopts the R-value scoring method to evaluate the efficiency of this proposed model.The results show that the prediction of the strain rate field in the three stages passes the R-value scoring.

Keywords：GNSS；strain rate field；earthquake prediction；R-value scoring；abnormal indicators；Yunnan area