2021, Vol. 41
震前电离层扰动与电离层耦合机制是地震前兆研究的热点之一[1-3]。研究表明，在5级以上地震发生前，震中及其附近地区的电离层TEC会存在异常扰动的现象[2-6]。传统的探测震前TEC异常的方法主要包括中位数法、四分位距法和滑动时窗法等[6], 虽然这些方法已成功分析了许多震前TEC异常现象，但其计算TEC参考背景值的精度较低，同时又缺乏统一的限差判决门限，因此会影响电离层异常扰动事件的准确判断[5]。基于上述原因，许多新的震前TEC异常探测方法[6-8]被提出，其中基于时间序列分析理论的方法受到大量关注。但该方法主要用于7级以上的地震震例分析，对5、6级地震震例的研究很少。本文利用ARMA模型对TEC预测数据进行残差修正，获得TEC参考背景值，提出相对误差和绝对误差相互制约的限差判决门限，并以2016-01-21青海门源MS6.4地震为例，研究震前14 d和地震当天震中及周围区域的TEC异常。
1 残差修正的ARMA模型(ARMA-RC)ARMA模型利用动态数据本身的结构与规律，应用数理统计方法进行预报[7]。如果基于历史数据Xt, Xt-1, …在t时刻预测未来值Xt+l(l > 0), 则可证明
理论上，如果ARMA模型成立，则残差序列
| $\hat{X}_{t+l}(R)=\hat{X}_{t+l}+\hat{\varepsilon}_{t+l}$ | (1) |
式中，
基于IGS提供的2014-11-01~12-29网格点(0°N, 90°E)TEC数据，分别采用ARMA-RC模型、滑动时窗法(时窗长度为10 d)和四分位距法预测2014-12-25~29共120 h的TEC小时数据。引入相对误差(Ur)和绝对误差(Δ)对结果进行分析：
| $\left\{\begin{array}{l}U_{r}=\frac{\left|\mathrm{TEC}_{\mathrm{pre}}-\mathrm{TEC}_{\mathrm{IGS}}\right|}{\mathrm{TEC}_{\mathrm{IGS}}} \times 100 \% \\ \Delta=\mathrm{TEC}_{\mathrm{pre}}-\mathrm{TEC}_{\mathrm{IGS}}\end{array}\right.$ | (2) |
式中，TECpre、TECIGS分别为TEC预报值和IGS发布值。
由表 1可知，在120 h预报时间内，无论是从|Δ| < 2 TECu的数据量占全部数据的比例，还是|Δ| < 1 TECu所表示的较高精度来看，ARMA-RC模型的预测精度都最高。图 1横轴为从预报开始时刻起算的小时数，可以看出，3种方法中，ARMA-RC模型预测TEC值的绝对误差和相对误差波动范围都最小，说明其预测TEC的精度高于传统方法。
|
表 1 3种方法预测TEC值的误差统计 Tab. 1 Error statistics of predicted TEC of three methods |
|
图 1 ARMA-RC模型与传统方法预测TEC的误差对比 Fig. 1 Comparison of errors of predicted TEC between ARMA-RC model and traditional methods |
门源MS6.4地震发生于2016-01-21 UT 01:13, 震中位置为37.7°N、101.6°E。基于IGS发布的时间分辨率为1 h的TEC数据，通过双线性插值求出震中位置TEC值，数据时段为2015-11-14~2016-01-21。
3.1 ARMA-RC模型计算TEC参考背景值将震前数据分成3组：2015-11-14~2016-01-06、2015-11-19~2016-01-11、2015-11-24~2016-01-16, 分别作为训练数据用于参数估计及模型建立。将震前14 d和地震当天分为3个时段：2016-01-07~11、2015-11-12~2016-01-16、2015-11-17~2016-01-21, 其相应数据作为观测值。
基于训练数据，参照文献[7]的计算流程，利用ARMA-RC模型，分别对未来5 d的TEC小时数据进行滑动递推预测，得到2016-01-07~11、2015-11-12~2016-01-16、2015-11-17~2016-01-21时段TEC参考背景值。对新加入的训练数据进行异常剔除，处理方法为利用相邻内插法计算正常值并将其替换[6]。
3.2 TEC异常限差判定策略若数据点满足以下条件，则视为TEC异常：1)预测残差|ΔTEC| > 2σ1; 2)相对误差|ΔTEC|/TECARMA > 2σ2; 3)探测时段未发生地磁扰动或太阳异常活动。其中，σ1、σ2分别为预测残差和相对误差数据的标准差。设εn为数据序列，预测时长为T(本文T=120), 利用式(3)可计算预测数据的均值εt和标准差σt:
| $\left\{\begin{array}{l}\bar{\varepsilon}_{t}=\frac{1}{T} \sum\limits_{n=t+1}^{t+T} \varepsilon_{n}, t=0, 1, 2, \cdots \\ \sigma_{t}=\sqrt{\frac{1}{T-1} \sum\limits_{n=1}^{T}\left(\varepsilon_{n}-\bar{\varepsilon}_{t}\right)^{2}}, t=0, 1, 2, \cdots\end{array}\right.$ | (3) |
基于门源县TEC数据，初步发现10 d共17个时刻出现TEC异常。为排除太阳、地磁异常活动等因素对TEC数据的影响，对地震前14 d各时段地磁活动的Kp指数、Dst指数和太阳辐射流量F10.7进行分析。从图 2可以看出，在探测时段内，F10.7远小于150 SFU, 指数变化相对稳定，可以排除太阳活动引起TEC异常的可能性；01-07、01-11~13、01-19~21地磁活动Kp指数均超过3, 表明存在一定的地磁干扰；除01-21外，其他日期Dst指数绝对值均在30 nT以内。因此，剔除上述日期内的TEC异常后，其他异常极可能由此次地震引起。
|
图 2 震前(01-07~01-21)Dst、Kp和F10.7指数变化 Fig. 2 Variation of Dst, Kp and F10.7 before Menyuan earthquake (from January 7 to 21) |
表 2为2016-01-07~21时段TEC异常的相关数据，其中，ΔTEC=TECIGS-TECARMA, ΔTEC > 0表示正异常，ΔTEC < 0表示负异常；正负异常持续时间(ΔT)为ΔTEC连续同号的小时数。图 3(a)、3(b)分别为利用滑动时窗法和ARMA-RC模型探测的TEC异常结果，数字1~7表示异常数据点，EQ为地震发生时刻位置，绿色区域表示TEC参考背景值上下限值区间，黑色曲线为实际探测的TEC值。
|
|
表 2 门源地震震中上空TEC异常扰动数据 Tab. 2 TEC anomalies data over epicenter of Menyuan earthquake |
|
图 3 门源地震前(2016-01-07~21) TEC异常扰动结果 Fig. 3 The results of TEC anomalies before Menyuan earthquake (from January 7 to 21, 2016) |
由图 3可知，2种方法探测的结果存在明显不同。滑动时窗法由于预报的参考背景值精度较低，其发生异常的天数明显多于ARMA-RC模型，且异常一旦出现，随后几天均会出现异常，同时正负异常无明显分布规律，该结论与前人的讨论结果相符[1, 6]。由于滑动时窗法探测结果的准确性存在一定问题，因此本文只对ARMA-RC模型的探测结果作进一步分析。
由表 2及图 3(b)可知，震前14 d内共发生4次正异常和3次负异常，正负异常次数基本相同；正异常最长持续时间为3 h, 在震前第5天UT 05:00达到此次正异常峰值(位置7), 该点也为震前正异常的最大值时刻；负异常最长持续时间为3 h, 在震前第11天UT 09:00达到此次负异常峰值(位置4), 但负异常最大值时刻并未与该点重合，而是发生在震前第5天UT 00:00(位置6)。对于此次地震而言，正负异常的最大值发生处均为最临近地震的2次异常，并且负异常发生时间早于正异常。
3.4.2 考虑IGS数据解算精度的结果分析IGS发布的电离层TEC数据的解算精度为2~8 TECu, 为进一步验证表 2中TEC异常结果的可靠性，对上述TEC异常作进一步筛选。综合考虑IGS数据的解算精度，制定如下异常筛选方案：1)将表 2中|ΔTEC|与IGS发布的RMS值进行对比，剔除|ΔTEC|远小于RMS的数据点；2)通过区域二维电离层地图进行异常分析。根据方案1, 表 2中序号6、7的数据点可作为进一步待分析的异常候选点，其他数据点受IGS数据解算精度所限而无法有效判断，暂不进行分析。
为进一步分析震前TEC异常的二维空间分布情况，利用ARMA-RC模型对0°~60°N、50°~150°E范围内的GIM电离层地图进行异常探测。探测数据时段与震中相同，时间分辨率为1 h, 空间分辨率为经度5°、纬度2.5°, 共计525个探测点。图 4为2016-01-15 UT 23:00~2016-01-16 UT 07:00的TEC异常分布图，图中五角星为震中位置，图 4(b)、4(g)分别对应表 2中序号6、7的震中异常时刻。
|
图 4 ARMA-RC模型探测区域电离层异常(2016-01-15 UT 23:00~2016-01-16 UT 07:00) Fig. 4 TEC anomalies based on ARMA-RC model(from UT 23:00, January 15, 2016 to UT 07:00, January 16, 2016) |
从图 4可以看出，UT 23:00震中区域电离层非常平静，震中西南方向(15°N, 70°E)附近出现明显负异常，最大异常峰值为-2.9 TECu; 随着时间的推移，异常区域向震中方向移动并逐渐增大，且呈现集中现象，在UT 00:00形成覆盖震中的负异常区域，震中位于该区域的偏东方向，最大异常峰值为-2.6 TECu, 位于震中以南；随后异常区域远离震中并向东南方向扩散，在UT 01:00形成较大范围、分散的负异常区域，最大异常峰值为-3.7 TECu, 并在10°N、150°E附近出现明显正异常区域，最大异常峰值为4.4 TECu; UT 02:00震中东南方向出现异常峰值为3.8 TECu和-4.6 TECu的2个异常区域，查询2016年全球地震动态可知，2016-01-27加罗林群岛西部(8.3°N, 137.8°E)、2016-01-19中国台湾台东县海域(22.9°N, 121.3°E)分别发生5.2级和5.8级地震，2次地震震中恰好位于图 4(d)中正、负异常区域，2次异常可能与这2次地震有关，具体关联还需进一步研究；UT 03:00震中附近电离层又趋于平静；UT 04:00震中西部出现峰值为-2.7 TECu的负异常；UT 05:00形成覆盖震中的正异常区域，震中位于该区域的偏南方向，正异常峰值同样位于震中以南，最大异常峰值为2.8 TECu; 随后异常区域向东北方向移动并逐渐减小，震中附近区域异常再次消失。
综上所述，2016-01-16 UT 00:00和UT 05:00在震中及其附近形成范围非常集中的电离层异常区域，异常峰值点均位于震中偏南方向，同时震中位置并不是异常区域的中心点：UT 00:00发生负异常时，震中位于负异常区域的偏东方向；UT 06:00发生正异常时，震中位于正异常区域的偏南方向。上述分析特征与前人的研究结论基本一致[4-6, 11], 表明2016-01-16 UT 00:00和UT 05:00的电离层异常现象与本次地震密切相关。
4 结语本文利用ARMA模型对TEC预测残差进行修正，得到TEC参考背景值，并设定相对误差和绝对误差相互制约的限差判定策略，探测门源MS6.4地震前TEC异常。结果表明：1)震前第5天存在明显的TEC异常现象，且异常与太阳活动、地磁活动等影响无关，从而验证了本方法对震前TEC异常探测的有效性；2)震前TEC正负异常发生次数基本相同，并且负异常发生在正异常之前；3)震中及附近更大区域内形成范围非常集中的电离层异常区域，异常峰值点位于震中以南，震中位置与异常区域的中心点并不重合，位于负异常区域偏东方向和正异常区域偏南方向。
本方法可为震前TEC异常探测中参考背景值的确定和限差判决门限的设定提供重要参考。为了更好地验证本文方法的有效性，下一步将选取更多、更经典的震例进行地震前后异常情况探测。
致谢: 感谢IGS、日本京都地磁中心、中国科学院空间中心为本文提供电离层、Kp、Dst和F10.7数据。
| [1] |
汤俊, 毛文飞. 多尺度ARMA残差修正模型震前电离层TEC异常探测[J]. 武汉大学学报: 信息科学版, 2019, 44(6): 791-798 (Tang Jun, Mao Wenfei. Multi-Scale ARMA Residual Correction Model for Detecting Pre-Earthquake Ionospheric Anomalies[J]. Geomatics and Information Science of Wuhan University, 2019, 44(6): 791-798)
(  0) |
| [2] |
Şentürk E, Çepni M S. A Statistical Analysis of Seismo-Ionospheric TEC Anomalies before 63 MW≥5.0 Earthquakes in Turkey during 2003-2016[J]. Acta Geophysica, 2018, 66(6): 1 495-1 507 DOI:10.1007/s11600-018-0214-2
( 0) |
| [3] |
Shah M, Tariq M A, Ahmad J, et al. Seismo Ionospheric Anomalies before the 2007 M7.7 Chile Earthquake from GPS TEC and DEMETER[J]. Journal of Geodynamics, 2019, 127: 42-51 DOI:10.1016/j.jog.2019.05.004
( 0) |
| [4] |
马一方, 匡翠林, 周晓慧. 基于GPS数据的震前电离层异常分析[J]. 测绘通报, 2015(4): 1-4 (Ma Yifang, Kuang Cuilin, Zhou Xiaohui. Investigation of Ionospheric VTEC Anomalies before Earthquakes Based on GPS Data[J]. Bulletin of Surveying and Mapping, 2015(4): 1-4)
( 0) |
| [5] |
聂兆生, 祝芙英, 付宁波. Kalman滤波在地震电离层TEC异常探测中的应用[J]. 大地测量与地球动力学, 2011, 31(3): 47-50 (Nie Zhaosheng, Zhu Fuying, Fu Ningbo. Application of Kalman Filtering in Detecting Ionospheric TEC Anomaly Prior to Earthquake[J]. Journal of Geodesy and Geodynamics, 2011, 31(3): 47-50)
( 0) |
| [6] |
张小红, 任晓东, 吴风波, 等. 震前电离层TEC异常探测新方法[J]. 地球物理学报, 2013, 56(2): 441-449 (Zhang Xiaohong, Ren Xiaodong, Wu Fengbo, et al. A New Method for Detection of Pre-Earthquake Ionospheric Anomalies[J]. Chinese Journal of Geophysics, 2013, 56(2): 441-449)
( 0) |
| [7] |
李磊. 基于GNSS的电离层总电子含量的预测与应用研究[D]. 大连: 大连海事大学, 2015 (Li Lei. Research on Forecasting and Application of the Ionospheric Total Electron Content Based on GNSS[D]. Dalian: Dalian Maritime University, 2015)
( 0) |
| [8] |
杨力, 赵海山, 董明, 等. 日本九州岛地震震前电离层TEC异常[J]. 测绘学报, 2016, 45(增2): 139-146 (Yang Li, Zhao Haishan, Dong Ming, et al. Ionospheric Anomaly before Kyushu, Japan Earthquake[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(S2): 139-146)
( 0) |
| [9] |
杨叔子, 吴雅, 轩建平, 等. 时间序列分析的工程应用[M]. 武汉: 华中科技大学出版社, 2007 (Yang Shuzi, Wu Ya, Xuan Jianping, et al. Time Series Analysis in Engineering Application[M]. Wuhan: Huazhong University of Science and Technology Press, 2007)
( 0) |
| [10] |
刘军, 柴洪洲, 常宜峰, 等. 改进的修正预测法预报电离层[J]. 测绘科学技术学报, 2011, 28(1): 19-22 (Liu Jun, Chai Hongzhou, Chang Yifeng, et al. Predicting Ionospheric Using Improved Correcting Prediction Method[J]. Journal of Geomatics Science and Technology, 2011, 28(1): 19-22)
( 0) |
| [11] |
史坤朋, 郭金运, 刘智敏, 等. 2016-12-25智利MW7.6地震震前电离层TEC异常探测[J]. 大地测量与地球动力学, 2018, 38(9): 979-984 (Shi Kunpeng, Guo Jinyun, Liu Zhimin, et al. Detection of Ionospheric TEC Anomaly before Chile MW7.6 Earthquake on December 25, 2016[J]. Journal of Geodesy and Geodynamics, 2018, 38(9): 979-984)
( 0) |