2018, Vol. 38
2. 山东科技大学矿山灾害预防控制省部共建国家重点实验室，青岛市前湾港路579号，266590;
3. 河海大学地球科学与工程学院, 南京市佛城西路8号，211100
2 016-12-25智利发生MW7.6大地震(43.4°S, 73.9°W)，是典型的发生在2个板块俯冲带之间的逆冲型地震，震源深度38.0 km。随着电磁观测技术的快速发展，地震电磁学已经成为地震预报监测的重要手段之一。许多学者已经证实地震会影响电磁场，进而造成电离层的异常扰动[1-2]。美国阿拉斯加大地震时首次探测到地震震前的电离层异常扰动现象[3]。之后，更多的学者开始研究震前电离层异常效应，探测到在大地震发生的前几天甚至前几小时，震中及其附近的电离层中总电子含量(TEC)大多存在扰动异常[4-6]。随着GNSS技术的发展，其价格低廉、容易获取的观测数据被用于探测震前电离层异常[7-8]。
传统的探测手段包括静态探测和动态探测，静态探测包括中位数法、平均值法等，动态探测以四分位距法为主。无论是哪种探测方法，都存在参考背景值精度较低的问题，且容易受到多种外界环境因素的干扰，这也是造成目前传统方法研究震前电离层扰动异常结果不一的重要原因[9]。本文利用奇异谱分析(SSA)对智利MW7.6地震进行震前电离层异常探测，意在提取TEC的主成分，获得高精度的背景参考值，然后将SSA探测结果和IQR探测结果对比分析，进一步研究电离层异常的时空变化规律，更加有效地提取地震引起的电离层异常信息。
1 数据来源 1.1 TEC数据采用IGS提供的GIM数据(ftp://cddis.gsfc.nasa.gov/gnss/products/ionex/)对智利地震进行震前电离层异常探测。GIM小范围的电离层异常扰动可以通过平滑或者内插过滤掉，因此其非常适合进行大尺度的震前电离层异常探测[10]。
在TEC时间序列选择上，一方面，太阳自转周期为27 d；另一方面，X射线、紫外线等太阳活动影响地球的地磁活动，从而导致TEC产生27 d的周期变化。本文首先选取距离智利地震震中最近的4个格网点的GPS TEC值，通过双线性插值获取震中2016-11-29~12-25共27 d的TEC时间序列，利用SSA和IQR进行异常探测，找到其主要的异常时间段，然后提取相应时段的全球格网TEC数据，分析电离层空间异常特征，并判断其和地震的关系。
1.2 太阳和地磁数据电离层容易受到太阳和地磁活动的干扰，需在电离层异常探测前对太阳和地磁活动进行分析。本文采用太阳射电通量F10.7反映太阳活动，采用赤道地磁活动指数Dst和全球地磁活动指数Kp反映地磁活动。F10.7数据由中国科学院空间中心提供(http://www.sepc.ac.cn)，时间分辨率为1 d，太阳平静状态下小于100 SFU[11]。Dst指数和Kp指数分别由日本京都地磁数据中心(http://wdc.kugi.kyoto-u.ac.jp/dstae/index.html)和美国国家海洋和大气局(ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/KP_AP/2016)提供。Dst指数分辨率为1 h，地磁平静状态下-30~30 nT[12]；Kp指数时间分辨率为3 h，平静状态下小于4[13]。
2 电离层异常探测方法 2.1 奇异谱分析奇异谱分析可以较好地从含噪声的有限尺度时间序列中提取信息，目前已应用于多种时间序列的分析。
奇异谱分析大致分为嵌入、SVD分解、分组和重建4步。对给定兴趣点的27 d的TEC数据，组成x1，x2，…，xN的时间序列{x}，建立时滞矩阵X：
| $ \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {{x_1}}&{{x_2}}& \cdots &{{x_{N - M + 1}}}\\ {{x_2}}&{{x_3}}& \cdots &{{x_{N - M + 2}}}\\ \cdots&\cdots&\cdots&\cdots \\ {{x_4}}&{{x_{M + 1}}}& \cdots &{{x_N}} \end{array}} \right] $ | (1) |
式中，N为时间序列长度，M为嵌套维数(时间窗口)，1≤M≤
| $ x \approx y = \sum\limits_k^p {{x^k}} {\rm{, }}\;{\rm{ }}1 \le p \le l $ | (2) |
SSA提取的TEC时间序列主成分TECmain包含电离层短周期信号，可以获得高精度的背景参考值。利用这27 d的TEC数据计算背景噪声ε，取两倍ε作为上下限值，由此得到电离层TEC的正常变化范围:
| $ \left. \begin{array}{l} {\rm{up = TE}}{{\rm{C}}_{{\rm{main}}}}{\rm{ + }}2\varepsilon \\ {\rm{low = TE}}{{\rm{C}}_{{\rm{main}}}}{\rm{ - }}2\varepsilon \end{array} \right\} $ | (3) |
以up和low作为TEC序列变化的上下界，超过该范围则认为是异常值。
2.2 四分位距法四分位距是一种稳健统计技术中用于表示数据离散度的一个量，常用来检查数据的异常情况。本文选取27 d作为时间窗口，在分析每天的数据时取前27 d同一时刻的TEC组成一组时间序列，并按照从小到大的顺序对其进行排列为x1, x2, …, x27，计算上四分位数Q1、中四分位数Q2、下四分位数Q3和四分位距IQR[15]：
| $ {Q_1} = {x_7}, {Q_2} = {x_{14}}, {Q_3} = {x_{21}} $ | (4) |
| $ {\rm{IQR}} = {Q_3} - {Q_1} $ | (5) |
采用Q2±1.5IQR作为TEC异常检验上下边界的判断阈值，如果TEC值没有超过上下界，则ΔTEC=0；超过上界ΔTEC＞0，表示正异常；低于下界值ΔTEC＜0，表示负异常。
2.3 可行性分析采用2种方法对近年来的2个大地震(东日本地震(2011-03-11, MW9.0)和尼泊尔地震(2015-04-25, MW8.1))震中时间序列进行异常探测。
图 1和2分别给出东日本和尼泊尔地震的探测结果。可以看出，无论是哪种方法，最主要的TEC异常时段都一样，均发生在2个地震震前3 d，这和马一方等[16]及姚宜斌等[17]的探测结果一致，说明2种方法均具有可行性。SSA探测TEC异常时段较少，振幅和临震时间呈正比，但TEC异常幅度小，将2种方法结合有利于地震短期预报。
|
图 1 东日本地震震中时间序列探测结果 Fig. 1 Analysis of TEC series of east Japan earthquake epicenter |
|
图 2 尼泊尔地震震中时间序列探测结果 Fig. 2 Analysis of TEC series of Nepal earthquake epicenter |
太阳和地磁活动会对电离层造成很大的扰动，在分析电离层扰动前应该考虑这2个因素的变化情况。图 3分别给出了智利地震震前近1个月的F10.7、Dst和Kp指数的变化情况。由图 3(a)可以看出，震前1个月太阳活动比较平静，对电离层的影响比较小。图 3(b)、3(c)显示，地磁主要在2个阶段比较活跃，一个是震前16 d，Kp指数达到了5，Dst指数超过了30 nT，该天共有6 h达到了小型磁暴的水平(UTC 18：00~24：00);另一个阶段是震前4~3 d，Kp指数达到了6，Dst指数同样超过了30 nT，地磁活动共有3 h达到了中型磁暴水平(震前4 d UTC15：00~18：00)，3 h达到小型磁暴水平(震前3 d UTC 00:00~03：00)。
|
图 3 太阳地磁环境 Fig. 3 Solar-terrestrial environment variations |
利用SSA和IQR对智利震中震前26~0 d的TEC时间序列分别进行异常探测，图 4、5分别为IQR、SSA探测的震中TEC异常结果。
|
图 4 智利地震IQR异常探测结果 Fig. 4 TEC anomalies over Chile earthquake epicenter based on IQR |
|
图 5 智利地震SSA异常探测结果 Fig. 5 TEC anomalies over Chile earthquake epicenter based on SSA |
从图 4可以看出，智利震中的TEC值日变化比较明显，通常情况下，TEC日时间序列具有单峰现象，但这27 d观测数据中多天出现了2个峰值。IQR探测出2个TEC主要异常时间段，一个是震前4 d(2016-12-21)出现明显的正异常，异常值超过7 TECu；另外一个是震前6 d(2016-12-19)出现明显的负异常，异常值达4.1 TECu。其余天数也出现振幅较小的TEC正异常，比如震前16 d TEC异常达3.2 TECu，震前19 d和震前23~22 d TEC异常也超过2 TECu。由图 4可知，震前16 d出现小型磁暴，且震前23 d秘鲁地区发生MW6.2地震，可能会对智利地震震中上方的TEC造成一定的影响。
从图 5可以看出，由SSA提取的TEC主成分包含TEC时间序列的日周期信息，SSA探测出的2个时间段和IQR异常探测的主要时间段一致。其中，震前4 d(2016-12-21)出现TEC正异常，异常值为6.0 TECu；震前6 d(2016-12-19)出现TEC负异常，异常值为2.7 TECu。相比IQR探测，SSA探测的TEC异常值较小，但出现在震前较短的时间段内，异常的时间段数也明显减少，对于地震的短期预报具有一定的意义。为了判断是否和地震有关，本文利用2种方法从空间上探测这2个时段电离层的异常特征。
5 空间电离层异常探测以2016-11-29~12-25的全球电离层格网数据为背景值，2倍中误差作为上下限，利用SSA进行探测电离层异常。如果TEC值与背景值之差小于2ε，表明没有明显异常；如果大于2ε，且ΔTEC＞0，表明有正异常；如果大于2ε，且ΔTEC＜0，表明有负异常。
此外，分别提取这2 d前27 d的全球电离层格网数据，按照四分位距方法确定这2 d的TEC上下限，利用滑动四分位法判断相应时段的TEC数据是否超限，对四分位距和SSA的探测结果进行分析。
5.1 2016-10-19空间异常探测图 6为IQR探测的2016-10-18(震前7 d)22:00~10-19(震前6 d)04:00(UTC)区域TEC异常分布，图 7为SSA探测的区域TEC异常。可以看出，SSA和IQR探测的TEC异常分布基本一致，TEC异常区域的运动趋势(从东向西)、形成和消散过程基本一样。2016-10-18 22:00(UTC)大西洋南部出现小幅度的TEC负异常。2016-10-19 00:00 TEC负异常区域到达智利地震震中上空，此时TEC异常幅度值最大；02:00 TEC异常区域继续向西移动，在地震震中西南部达到该时间段的最大范围；04:00 TEC负异常开始消散。
|
图 6 IQR探测2016-10-18~10-19区域电离层异常 Fig. 6 TEC anomalies prior to earthquake based on IQR from October 18th to 19th, 2016 |
|
图 7 SSA探测2016-10-18~10-19区域电离层异常 Fig. 7 TEC anomalies prior to earthquake based on SSA from October 18th to 19th, 2016 |
进一步对比发现，IQR异常时间较长，达到8 h，最大异常幅度较大，为5.5 TECu；而SSA探测TEC异常区域比较小，异常峰值点更加靠近震区中心，这有利于对地震震中位置的预测。该天震中区域外电离层比较平静，地磁活动和太阳活动干扰小，结合2种探测方法认为，该时间段的TEC异常是地震前兆信息。
5.2 2016-10-21空间异常探测图 8给出了IQR探测的2016-10-21(震前4 d)14:00~22:00(UTC)区域TEC异常分布。14:00大西洋南部出现小幅度的电离层异常；16:00电离层异常区域突然扩大，达到该时段的最大范围，除大西洋南部，整个南极洲东部地区也发生小幅度的电离层异常扰动，并出现向西运动的趋势；18:00电离层异常区域到达了震区南部，共轭区出现小范围的电离层异常，但振幅很小，不超过2 TECu；20:00电离层异常幅度均达到最大，TEC异常最大值位置靠近震中，TEC异常值达到16 TECu，共轭区域电离层的正扰动同样比较明显；22:00起异常区域开始消散，电离层恢复平静。从图 3可以看出，有一个时段Kp指数达到6(15：00~18:00)，Dst指数主相变化不强烈(Dst＜50)，表明该时段发生缓始型中型磁暴[18]。这类磁暴造成的TEC最大扰动通常出现在Dst峰值附近，且幅度较小，一般为0.2~2.5 TECu[19-20]，16:00南极洲上空的TEC异常很可能由中型磁暴引起。
|
图 8 IQR探测2016-10-21区域电离层异常 Fig. 8 TEC anomalies prior to earthquake based on IQR on October 21th, 2016 |
图 9为SSA探测的区域TEC异常分布。2种方法的不同点主要在异常范围、异常幅度和共轭区域3个方面。SSA探测结果显示，电离层异常范围更加集中，异常值峰值点更加靠近地震震中，16:00 SSA探测南极洲上空异常范围比IQR小得多，在一定程度上削减了地磁活动造成的TEC扰动。此外，SSA探测共轭区域电离层异常更加明显，共轭区域的电离层异常是判断地震发生的一个重要标准，对于预测地震有很大的帮助。但相比于IQR探测，其探测到的异常振幅较小，最大值为12 TECu。结合全球背景，该天的TEC异常可以认为是地震发生的前兆。
|
图 9 SSA探测2016-10-21区域电离层异常 Fig. 9 TEC anomalies prior to earthquake based on SSA on October 21th, 2016 |
2016-10-19和2016-10-21的异常探测结果显示，电离层异常可以作为地震预测的有效方法。随着地震发生的临近，震区电离层异常幅度、范围和持续时间都更加明显。对比2种方法发现，SSA可以削弱地磁活动对电离层的影响，相比IQR探测，其TEC异常峰值更加靠近地震震中，且共轭区域出现明显的扰动，但异常范围较小，2种方法可以相互补充以达到更好的探测效果。值得一提的是，通常情况下震前电离层异常区域出现在赤道异常双峰线附近[21]，但智利地震震中位于中高纬地区，结合全球背景仍然出现了比较明显的地震电离层异常，这对于发生在中高纬地区的地震震前电离层异常探测有着积极的意义。
6 结语本文基于奇异谱分析法和四分位距法，利用GIM数据对2016-12-25智利MW7.6地震震中TEC时间序列和空间电离层进行异常探测及分析。结果表明，震前7 d震中南部出现明显的电离层负异常扰动，震前4 d震中南部出现明显的电离层正异常扰动，同时共轭区有明显的电离层正异常变化。2个时间段的电离层异常区域具有同样的运动趋势，均自东向西移动，异常幅度和范围均在震中附近达到最大。将2种方法互为补充可有效地突出主要异常时间段及孕震区和共轭区的异常特征。此次智利地震震中属于中高纬地区，仍可探测到明显的地震电离层异常效应，这对于该区域的地震短期预报有一定的意义。
| [1] |
Hayakama M. Atmospheric and Ionospheric Electromagnetic Phenomena Associated with Earthquakes[J]. J Atmos Sol-Terr Phys, 2000, 62(3): 225-227 DOI:10.1016/S1364-6826(99)00110-8
(  0) |
| [2] |
Le H, Liu J Y, Liu L. A Statistical Analysis of Ionospheric Anomalies before 736 M6.0+ Earthquakes during 2002-2010[J]. J Geophys Res, 2011, 116(A2)
( 0) |
| [3] |
Leonard R S, Barnes R A. Observations of Ionospheric Disturbances Following the Alaska Earthquake[J]. J Geophys Res, 1965, 70(5): 1 250-1 253 DOI:10.1029/JZ070i005p01250
( 0) |
| [4] |
Liu J Y, Chen Y I, Pulinets S A, et al. Seismo-Ionospheric Signatures Prior to M≥6.0 Taiwan Earthquakes[J]. Geophys Res Lett, 2000, 27(19): 3 113-3 116 DOI:10.1029/2000GL011395
( 0) |
| [5] |
Liu J Y, Chen Y I, Jhuang H K, et al. Ionospheric foF2 and TEC Anomalous Days Associated with M≥5.0 Earthquakes in Taiwan during 1997-1999[J]. Terr Atmos Ocean Sci, 2004, 15(3): 371-384 DOI:10.3319/TAO.2004.15.3.371(EP)
( 0) |
| [6] |
Zakharenkova I E, Shagimuratov I I, Tepenitzina N Y, et al. Anomalous Modification of the Ionospheric Total Electron Content Prior to the 26 September 2005 Peru Earthquake[J]. J Atmos Sol-Terr Phys, 2008, 70(15): 1 919-1 928 DOI:10.1016/j.jastp.2008.06.003
( 0) |
| [7] |
Calais E, Minster J B. GPS Detection of Ionospheric Perturbations Following the January 17, 1994, Northridge Earthquake[J]. Geophys Res Lett, 1995, 22(9): 1 045-1 048 DOI:10.1029/95GL00168
( 0) |
| [8] |
曾中超, 张蓓, 方广有, 等. 利用DEMETER卫星数据分析汶川地震前的电离层异常[J]. 地球物理学报, 2009, 52(1): 11-19 (Zeng Zhongchao, Zhang Bei, Fang Guangyou, et al. The Analysis of Ionospheric Variations before Wenchuan Earthquake with DEMETER Data[J]. Chinese Journal of Geophysics, 2009, 52(1): 11-19)
( 0) |
| [9] |
祝芙英, 吴云, 林剑, 等. 汶川MS8.0地震前电离层TEC异常分析[J]. 大地测量与地球动力学, 2008, 28(6): 16-21 (Zhu Fuying, Wu Yun, Lin Jian, et al. Study on Ionospheric TEC Anomaly Prior to Wenchuan MS8.0 Earthquake[J]. Journal of Geodesy and Geodynamics, 2008, 28(6): 16-21)
( 0) |
| [10] |
Mannucci A J, Wilson B D, Yuan D N, et al. Global Mapping Technique for GPS-Derived Ionospheric Total Electron Content Measurements[J]. Radio Sci, 1998, 33(3): 565-582 DOI:10.1029/97RS02707
( 0) |
| [11] |
Fejer B G, Paula E R D, Batista I S, et al. Equatorial F Region Vertical Plasma Drifts during Solar Maxima[J]. Journal of Geophysical Research, 1989, 94(A9): 12 049-12 054 DOI:10.1029/JA094iA09p12049
( 0) |
| [12] |
沈晓飞, 倪彬彬, 顾旭东, 等. 第23太阳活动周期太阳风参数及地磁指数的统计分析[J]. 地球物理学报, 2015, 58(2): 362-370 (Shen Xiaofei, Ni Binbin, Gu Xudong, et al. A Statistical Analysis of Solar Wind Parameters and Geomagnetic Indices for the Solar Cycle 23[J]. Chinese Journal of Geophysics, 2015, 58(2): 362-370)
( 0) |
| [13] |
Gordon R. Geomagnetic Indices[J]. Rev Geophys Space Phys, 1972, 10(4): 935-950 DOI:10.1029/RG010i004p00935
( 0) |
| [14] |
Hassani H. Singular Spectrum Analysis: Methodology and Comparison[J]. Journal of Data Science, 2007, 5(2): 239-257
( 0) |
| [15] |
Li W, Guo J Y, Yu X M, et al. Analysis of Ionospheric Anomaly Preceding the MW7.3 Yutian earthquake[J]. Geodesy and Geodynamics, 2014, 5(2): 54-60 DOI:10.3724/SP.J.1246.2014.02054
( 0) |
| [16] |
马一方, 匡翠林, 周晓慧. 基于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) |
| [17] |
姚宜斌, 翟长治, 孔建, 等. 2015年尼泊尔地震的震前电离层异常探测[J]. 测绘学报, 2016, 45(4): 385-395 (Yao Yibin, Zhai Changzhi, Kong Jian, et al. The Pre-Earthquake Ionosphere Anomaly of the 2015 Nepal Earthquake[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(4): 385-395)
( 0) |
| [18] |
李义红.缓始磁暴与电离层TEC扰动特点分析[C].第9届国家安全地球物理专题研讨会, 武汉, 2013 (Li Yihong. Research on Characteristic Analysis of Gradual Commencement Geomagnetic Storms and the Corresponding Ionosphere TEC Disturbance[C].The 9th National Security Geophysical Symposium, Wuhan, 2013) http://cpfd.cnki.com.cn/Article/CPFDTOTAL-ZGDW201310031072.htm
( 0) |
| [19] |
詹贤鋆. 磁暴的形态分类与选定[J]. 地震地磁观测与研究, 1981(1): 12-18 (Zhan Xianyun. The Morphology and Selection of Geomagnetic Storm[J]. Seismological and Geomagnetic Observation and Research, 1981(1): 12-18)
( 0) |
| [20] |
杨越.电离层受中小型磁暴影响特征的GNSS数据分析[D].沈阳: 东北大学, 2014 (Yang Yue. Analyzing of the Characteristics of Ionospheric Disturbances Derived by Medium and Small Geomagnetic Storms Based on GNSS Data[D]. Shenyang: Northeastern University, 2014) http://cdmd.cnki.com.cn/Article/CDMD-10145-1016011229.htm
( 0) |
| [21] |
Pulinets S A. Seismic Activity as a Source of the Ionospheric Variability[J]. Advances in Space Research, 1998, 22(6): 903-906 DOI:10.1016/S0273-1177(98)00121-5
( 0) |
2. State Key Laboratory of Mining Disaster Prevention and Control Cofounded by Shandong Province and Ministry of Science & Technology, Shandong University of Science and Technology, 579 Qianwangang Road, Qingdao 266590, China;
3. College of Earth Science and Engineering, Hohai University, 8 West-Focheng Road, Nanjing 211100, China