2021, Vol. 41
表1(Table 1)
2. 云南大理滇西北地壳构造活动野外科学观测研究站，云南省大理市滨海大道，671000
现阶段，GNSS站点主要分为连续站和区域站，区域站网点分布密集，空间分辨较高，可通过定期监测获取板块间的相对运动，来描述地壳趋势运动，解释区域地球动力学的背景等[1-2]；连续站观测从区域站的定期变为准实时，短周期信息较为丰富，具有较高的时间分辨率，能实时监测地壳运动的动态变化特征。GNSS连续观测的特点决定了其应该具有捕获地震短临异常前兆的能力。
不少学者在GNSS用于短临预测方面进行尝试，取得了一些进展。张晓亮等[3]通过数学形变滤波的方式处理单条基线的变化，从中分离出与地震时间有关的形变分量，用于提取地震前的短期形变异常。吴云等[4]利用GNSS基准站的时间序列，采用小波变换分析部分震例认为，GNSS基准站连续观测序列的低频段中，地震前6个月或稍长时段有一定的前兆异常出现，可作为地震中短期预测的依据。张燕等[5]利用小波变换对GNSS时间序列进行分析，捕获了大震前1~2 a各站点运动趋势的改变，并认为该变化可能是一种地震前兆信号。邵德盛等[2]基于云南地区28个GNSS连续站2011-01~2016-04数据，以21°~29°N、97°~106°E范围内的20个5级以上地震为例，提出面应变综合预测指标法，获取了适用于云南地区的地震短临异常识别指标。
利用GNSS观测资料，可进行应变参数的求解，获取区域应变场。应变参数主要包括面应变和最大剪应变等，各应变参数实际上就是对不同特性形变信息的一种客观分离[6]。面应变参数能直接反映区域挤压(收缩)或拉张(膨胀)强弱的特性，而最大剪应变是反映剪切应变强度的参数[6]。震源机制解研究表明，云南地区的M≥5.0地震多为走滑型地震，地震的破裂类型多含有剪切破裂，震前应变积累包含一定的剪切应变积累(走滑型地震以剪应变为主，而非走滑型地震也带有一定的走滑成分)[6]，因而，地震发生前最大剪应变会表现出一定的异常特征，利用最大剪应变资料提取适用于云南地区的短临异常指标并进行短临预测应该是可行的。
鉴于此，本文在前人研究的基础上，基于云南地区的43个GNSS连续站2011-01~2018-12数据，分别提取适用于云南地区的面应变及最大剪应变短临异常指标，并以其间发生的27个M≥5.0地震为样本，采用R值评分法[7]对异常指标进行预测效能评价，以期将GNSS连续观测资料更好地应用于云南地区的震情跟踪和预测预报工作。
1 数据收集整理本文的研究范围(21°~29°N，97°~107°E)内共包含43个GNSS连续站。GNSS数据为中国地震局GNSS数据产品服务平台(http://www.cgps.ac.cn)提供的站点位移时间序列原始数据，数据起止时间为2011-01~2018-12。该平台采用美国麻省理工学院(MIT)和加州大学圣地亚哥分校Scripps海洋研究所(SIO)研制的GAMIT/GLOBK软件(10.4版)进行数据处理。研究区2011-01~2018-12期间共发生27个M≥5.0地震。GNSS站点(点位信息来源于中国地震局GNSS数据产品服务平台(http://www.cgps.ac.cn))及地震震中分布(地震信息来源于中国地震台网(http://news.ceic.ac.cn))如图 1所示，相关地震参数见表 1。
|
图 1 GNSS连续站、地震震中分布及格网划分结果 Fig. 1 The grid division and the distribution of GNSS continuous stations and earthquakes |
|
表 1 云南地区27个M≥5.0地震事件列表(2011-01~2018-12) Tab. 1 List of 27 earthquakes with M≥5.0 in Yunnan area (2011-01-2018-12) |
使用云南地区43个GNSS连续站的站点位移时间序列原始数据进行应变场的求解。考虑到GNSS连续站分布的不均匀性，在获取了各个站点位移时间序列的基础上，采用克里金插值方法对点位数据进行1°×1°格网化插值(格网划分结果及编号见图 1)，获取均匀分布于80个格网的位移时间序列，再通过位移时间序列求取区域应变场[2, 8]，具体解算方法见文献[8]。
通过计算可获取面应变、最大剪应变、主应变等应变参数[6]。其中，面应变θ和最大剪应变γmax计算公式为[6, 8]：
| $ \left\{\begin{array}{l} \theta=\varepsilon_{x}+\varepsilon_{y} \\ \gamma_{\max }=\sqrt{\left(\varepsilon_{x}-\varepsilon_{y}\right)^{2}+\gamma_{x y}^{2}} \end{array}\right. $ | (1) |
本文的短临异常指标提取采用邵德盛等[2]提出的应变异常识别方法。在§2.1应变场计算的基础上，分别提取80个格网的应变(面应变、最大剪应变)时间序列，并进行线性去趋势处理，以获取长趋势变化背景下的短期异常波动。对去除趋势项后的时间序列，划定±2倍标准差作为单条曲线的异常指标线[2]。
图 2为24号格网的面应变和最大剪应变时间序列图。图 2(a)显示，该格网存在一定的趋势压缩，年压缩速率为-0.40×10-8，rms为±1.07。通过线性去趋势处理，并设定±2倍标准差作为异常阈值，可以明显看到(图 2(b))，该格网的面应变时间序列在腾冲M5.2、宁蒗M5.7地震前存在较为明显的超指标异常，即面应变时间序列值超过±2倍标准差线。图 2(c)显示，该格网的最大剪应变年速率为3.41×10-8，rms为±0.94。通过去除趋势项，并设定±2倍标准差作为异常阈值，可以明显看到(图 2(d))，该格网的最大剪应变时间序列在腾冲M5.2、彝良M5.7、洱源M5.0地震前存在较为明显的超指标异常，即最大剪应变时间序列值超过±2倍标准差线。
|
图 2 24号格网的面应变和最大剪应变时间序列 Fig. 2 Time series diagram of surface strain and maximum shear strain of grid No.24 |
在对应变格网时序进行线性去趋势处理后，需要对单个格网时间序列的映震能力进行评价。评价方式以±2倍标准差线作为阈值，超出该指标线的均为异常，并根据异常出现时间与地震事件的对应关系，对面应变、最大剪应变曲线的映震能力进行评分[2]。
由于需获取短临异常信息，因此地震事件间对应的时间间隔设定为不超过90 d(3个月)。假设在某时段内第i个格网面应变及最大剪应变时序曲线报准次数为Ni，漏报次数为Li，虚报次数为Xi，那么该曲线的最终评价得分为[2]：
| $ S_{i}=\frac{N_{i}}{N_{i}+L_{i}+X_{i}}, i=0, 1, \cdots, 79 $ | (2) |
经过以上处理，可分别获取相应时段80个格网的预测能力评分。若在该时段内第i、j、l、m格网出现异常，则该时段的异常总评分为[2]：
| $ S=S_{i}+S_{j}+S_{l}+S_{m} $ | (3) |
在此基础上通过逐个时间点对所有格网进行时空扫描，可分别获取云南地区的面应变及最大剪应变异常评分时间序列，称之为面应变及最大剪应变短临异常指标。当该异常评分值超出一定阈值后(此处设定2倍均值线作为异常阈值)，表示异常指标出现，认为云南地区存在发生M≥5.0地震的危险性。图 3为分别提取的2011-01~2018-12云南地区的面应变及最大剪应变短临异常指标。
|
图 3 云南地区的面应变及最大剪应变短临异常指标(2011-01~2018-12) Fig. 3 Short-term anomaly indicators of surface strain and maximum shear strain in Yunnan area (2011-01-2018-12) |
在分别提取了云南地区的面应变及最大剪应变短临异常指标后(图 3)，需要结合相应的地震事件对异常指标进行预测效能评价。现阶段应用最为广泛的预测效能评价方法是1973年由许绍燮院
士提出的R值评分法[9]。R值评分法主要依据预测指标实际预测地震的有震报准率c和预测占时率(预测时间占有率)b[7]，反映的是预测方法或指标与地震的相关程度。R值的计算公式为：
| $ \begin{gathered} \quad R=c-b= \\ \frac{\text { 预测有震命中次数 }}{\text { 地震的总次数 }}-\frac{\text { 预测有震的时间 }}{\text { 预测研究总时间 }} \end{gathered} $ | (4) |
式中，R=1表示全报对，R=0表示预测没有起作用。R值越大，预测效果越好。当预测效能R＞R0(R0为保证97.5%置信水平的最低R值)时[10]，表明所评估的预测方法通过统计检验，高于随机预测效能。
采用R值评分法对云南地区的面应变及最大剪应变短临异常指标进行效能评价，两种指标不同预测窗长的R值遍历曲线如图 4所示。由图可知，最大剪应变短临异常指标的R值分布较面应变短临异常指标更为集中，R值主要分布于预测窗长250 d以内，最优R值为0.360，所对应的预测窗长为30 d。而面应变短临异常指标的R值主要分布于预测窗长450 d以内，最优R值为0.280，相对应的预测窗长为40 d和70 d。
|
图 4 不同预测窗长的R值遍历曲线 Fig. 4 R value traversal curves of different prediction window lengths |
由于本文获取的是云南地区的短临异常指标，因此，在对异常指标进行预测效能评价时，预测窗长设定为不超过90 d(3个月)。在图 4的基础上，我们分别对面应变短临异常指标进行预测窗长为30 d、40 d、60 d、70 d、90 d的预测效能评价，对最大剪应变短临异常指标进行预测窗长为30 d、60 d、90 d的预测效能评价，具体效能评价结果见表 2和表 3。
|
|
表 2 面应变短临异常指标预测效能评价结果 Tab. 2 Evaluation results of prediction efficiency of surface strain short-term anomaly indicators |
|
|
表 3 最大剪应变短临异常指标预测效能评价结果 Tab. 3 Evaluation results of the prediction efficiency of the short-term anomaly index of the maximum shear strain |
由表 2可知，面应变短临异常指标在预测窗长为30 d、40 d、60 d、70 d、90 d时都通过了效能检验。其中，当预测窗长设置为30 d时，预测效能评分R为0.23，R0为0.21，R＞R0，对于实际发生的27个地震，该指标报对地震15次，漏报12次，准确率为55.56%，漏报率为44.44%；当预测窗长为60 d时，预测效能评分R为0.21，R0为0.20，R＞R0，该指标报对地震19次，漏报8次，准确率为70.37%，漏报率为29.63%；当预测窗长为90 d时，预测效能评分R为0.24，R0为0.19，R＞R0，该指标报对地震23次，漏报4次，准确率为85.19%，漏报率为14.81%；当根据最优R值设置预测窗长为40 d和70 d时，该指标预测准确率分别为66.67%和81.48%，漏报率分别为33.33%和18.52%。通过进一步对比发现，从30 d到90 d，该指标预测准确率逐渐升高，到90 d时预测准确率达到85.19%。
由表 3可知，当最大剪应变短临异常指标的预测窗长设置为30 d时(最优R值对应的预测窗长)，预测效能评分R为0.36，R0为0.20，R＞R0，该预测指标通过了效能检验，对于实际发生的27个地震，该指标报对地震21次，漏报6次，准确率为77.78%，漏报率为22.22%；当预测窗长为60 d时，预测效能评分R为0.25，R0为0.19，R＞R0，预测指标通过了效能检验，该指标报对地震23次，漏报4次，准确率为85.19%，漏报率为14.81%；当预测窗长为90 d时，预测效能评分R为0.12，R0为0.19，R＜R0，该预测指标未通过效能检验，其预测结果与60 d窗长相同。
对比表 2和表 3可知，当预测窗长设置为30 d时，面应变短临异常指标的预测效能评分R为0.23，R/R0为1.12，预测准确率为55.56%；最大剪应变短临异常指标的预测效能评分R为0.36，R/R0为1.82，预测准确率为77.78%，前者的R值、R/R0值和准确率均低于后者。当预测窗长设置为60 d时，面应变短临异常指标的预测效能评分R为0.21，R/R0为1.03，预测准确率为70.37%；最大剪应变短临异常指标的预测效能评分R为0.25，R/R0为1.30，预测准确率为85.19%，前者的R值、R/R0值和准确率均低于后者。当预测窗长为90 d时，二者的预测准确率相同，但是最大剪应变短临异常指标未通过效能检验。
4 讨论根据R值评分规则，在97.5%置信度下，当R＞R0值时，预测效能检验通过，因此，当异常指标超过控制线时可认为，云南地区未来30~90 d内发生M≥5.0地震的可能性较高(预测窗长为90 d时，最大剪应变短临异常指标除外)。纵观不同预测窗长的R值大小发现，两种异常指标的R值都不太高，主要原因是M≥5.0地震在云南地区的自然发生概率较高。
从两种异常指标的预测效能评价结果来看，当预测窗长设置为30 d或60 d时，无论是R值、R/R0值还是预测准确率，最大剪应变短临异常指标均优于面应变短临异常指标。究其原因，可能与本文所涉及到的27个M≥5.0地震多为走滑型地震，最大剪应变对走滑型地震前的应力变化过程更为敏感有关。
结合R值、R/R0值、预测准确率，对两种异常指标进行不同预测窗长的对比发现，面应变短临异常指标更适合于预测窗长为90 d的短临预测，而最大剪应变短临异常指标更适合预测窗长为60 d的短临预测。
结合图 3中的短临异常指标线来看，两种异常指标与地震事件之间均存在较好的对应关系。2014年云南地区相继发生05-30盈江6.1级、08-03鲁甸6.5级、10-07景谷6.6级地震。3次M≥6.0地震发生前，面应变短临异常指标分别于05-09、07-18、09-16出现明显的超指标异常，最大剪应变短临异常指标分别于05-07、07-11、09-22出现明显的超指标异常。两种异常指标都对3次地震作出了准确的预测，应该说3次M≥6.0地震的发生分别对两种异常指标进行了3次很好的检验。
5 结语本文通过对云南地区43个GNSS连续站2011-01~2018-12站点位移时间序列原始数据进行深加工处理，在计算应变参数格网时间序列的基础上，分别提取适合云南地区M≥5.0地震短临预测的面应变及最大剪应变短临异常指标；并基于其间发生的27个M≥5.0地震，采用R值评分法对异常指标进行预测效能评价，形成以下认识：
1) 面应变、最大剪应变都是重要的应变参数，基于二者分别提取的面应变及最大剪应变短临异常指标物理意义明确，两种指标对云南地区的M≥5.0地震的发震时间均具有较好的指示意义。
2) 通过对两种指标进行不同预测窗长的效能检验发现，当预测窗长设置为30 d或60 d时，在相同预测窗长下，最大剪应变短临异常指标的预测准确率高于面应变短临异常指标；通过对同一异常指标进行不同预测窗长的对比发现，面应变短临异常指标的最佳预测窗长为90 d，最大剪应变短临异常指标的最佳预测窗长为60 d。
3) 本文异常指标提取过程中设定±2倍标准差作为单个格网时间序列映震能力的异常阈值，设定2倍均值作为异常指标评分的异常阈值，主要是基于提高信噪比和异常识别能力而考虑的。对于选取更高或更低的异常阈值标准下，两种指标的短临预测和效能评价结果如何，有待进一步研究。
4) 现阶段，云南地区的GNSS连续站点密度较低，空间分辨率亦较低。随着中国地震科学实验场更多项目的开展，云南地区的GNSS连续站点数量将大幅提高，到时可基于更高空间分辨率的GNSS连续观测数据提取具有更强短临异常识别能力的综合指标，这对于云南地区的震情跟踪及预测预报工作具有重要的意义。
致谢: 本文GNSS数据来源于中国地震局GNSS数据产品服务平台(http://www.cgps.ac.cn)，数据分析软件由云南省地震局信息中心提供，在此一并表示感谢。
| [1] |
江在森, 马宗晋, 张希, 等. GPS初步结果揭示的中国大陆水平应变场与构造变形[J]. 地球物理学报, 2003, 46(3): 352-358 (Jiang Zaisen, Ma Zongjin, Zhang Xi, et al. Horizontal Strain Field and Tectonic Deformation of China Mainland Revealed by Preliminary GPS Result[J]. Chinese Journal of Geophysics, 2003, 46(3): 352-358 DOI:10.3321/j.issn:0001-5733.2003.03.012)
(  1) |
| [2] |
邵德盛, 洪敏, 张勇, 等. 云南地区形变观测资料短临异常指标提取[J]. 武汉大学学报: 信息科学版, 2017, 42(9): 1223-1228 (Shao Desheng, Hong Min, Zhang Yong, et al. Extract Short Impending Abnormal Indicator from Deformation Data in Yunnan Province[J]. Geomatics and Information Science of Wuhan University, 2017, 42(9): 1223-1228)
( 8) |
| [3] |
张晓亮, 江在森, 王敏, 等. 利用GPS连续站资料研究地壳运动与地震的关系[J]. 大地测量与地球动力学, 2006, 26(4): 63-68 (Zhang Xiaoliang, Jiang Zaisen, Wang Min, et al. On Relation between Crustal Movement and Earthquake by Use of Observations of Continuous GPS Sites[J]. Journal of Geodesy and Geodynamics, 2006, 26(4): 63-68)
( 1) |
| [4] |
吴云, 孙建中, 乔学军, 等. GPS揭示的现今地壳运动与地震前兆特征[J]. 大地测量与地球动力学, 2003, 23(3): 14-20 (Wu Yun, Sun Jianzhong, Qiao Xuejun, et al. Characteristics of Current Crust Movement and Seismic Precursors Revealed by GPS Survey[J]. Journal of Geodesy and Geodynamics, 2003, 23(3): 14-20)
( 1) |
| [5] |
张燕, 吴云, 施顺英, 等. GPS时间序列揭示地震前兆的初步探索[J]. 大地测量与地球动力学, 2005, 25(3): 96-99 (Zhang Yan, Wu Yun, Shi Shunying, et al. Preliminary Discussion on GPS Time Series Manifesting Earthquake Precursor[J]. Journal of Geodesy and Geodynamics, 2005, 25(3): 96-99)
( 1) |
| [6] |
江在森, 张希, 张晶, 等. 地壳形变动态图像提取与强震预测技术研究[M]. 北京: 地震出版社, 2013 (Jiang Zaisen, Zhang Xi, Zhang Jing, et al. Crustal Deformation Dynamic Image Extraction and Strong Earthquake Prediction Technology Research[M]. Beijing: Seismological Press, 2013)
( 5) |
| [7] |
许绍燮. 地震预报能力评分——震兆预报方法小结[R]. 中国地震局地球物理研究所, 北京, 1989 (Xu Shaoxie. Earthquake Prediction Ability Score-Summary of Earthquake Prediction Methods[R]. Institute of Geophysics, CEA, Beijing, 1989)
( 2) |
| [8] |
洪敏, 张勇, 邵德盛, 等. 云南地区近期地壳活动特征[J]. 地震研究, 2014, 37(3): 483-488 (Hong Min, Zhang Yong, Shao Desheng, et al. Recent Tectonic Activity Features of Yunnan Region[J]. Journal of Seismological Research, 2014, 37(3): 483-488)
( 3) |
| [9] |
陆远忠, 吴云, 王炜, 等. 地震中短期预报的动态图像方法[M]. 北京: 地震出版社, 2001 (Lu Yuanzhong, Wu Yun, Wang Wei, et al. Dynamic Image Processing Method of Mid Short-Term Earthquake Prediction[M]. Beijing: Seismological Press, 2001)
( 1) |
| [10] |
国家地震局科技监测司. 地震学分析预报方法程式指南[M]. 北京: 地震出版社, 1990 (Department of Science and Technology Monitoring, National Earthquake Administration. Program Guide for Seismological Analysis and Prediction Methods[M]. Beijing: Seismological Press, 1990)
( 1) |
2. Field Scientific Observation and Research Station on Crustal Tectonic Activities in Northwest Yunnan, Dali, Yunnan, Binhai Road, Dali 671000, China