2017年九寨沟MS7.0级地震震前电离层异常

朱军桃,赵苗兴,龚朝飞,王 雷

(桂林理工大学 测绘地理信息学院,广西 桂林 541006)

地震具有突发性强、破坏性大、防御难度大且社会影响深远等特点,严重危害人类的生命财产安全。据中国地震台网(http://www.ceic.ac.cn)统计,2016年9月至今,在世界范围内超过6.0级的地震有198起,造成的人员伤亡和经济损失让人们难以承受。由于地震产生原因复杂,影响因素众多,地震预报特别是短临预报目前仍然处于探索阶段[1]。1964年美国阿拉斯加大地震时,Leonard等[2]首次发现电离层扰动与地震之间存在某些联系,从此拉开了震前电离层异常扰动研究的序幕;Pulinets等[3-4]提出了利用卫星数据进行地震短临预报,统计分析了10年的地震电离层数据,确认电离层在地震前后有显著的异常现象;林剑等[5]利用GPS数据研究了汶川地震前后电离层的异常扰动,发现地震前后一周,孕震区上空连续出现电离层异常扰动,异常形态具有共轭结构,且呈现向赤道漂移趋势;童杨津等[6]采用滑动四分位法发现墨西哥M7.2地震前电离层TEC出现正异常现象;杨力等[7]采用滑动时窗法和相关性分析法研究了日本九州岛地震震前电离层异常变化。随着电离层探测技术的不断发展与完善,孕震区电离层异常扰动成了当前地震预报的热点之一。除了传统的电离层异常分析方法外,姚宜斌等[8]提出了一种基于奇异谱分析的电离层异常探测方法,结合三维层析发现2015年尼泊尔地震前出现大范围电离层正异常,异常随高度变化规律与电离层电子密度垂直分布规律相一致;张小红等[9]提出一种利用时间序列法(ARIMA模型)进行震前电离层异常探测的新方法,发现时间序列法预测背景值的精度要明显高于传统方法,且预报背景值的平均偏差要比传统方法小2倍左右;黄海莎等[10]介绍了震前电离层TEC异常探测原理的研究进展和主要的异常探测方法,总结了时间序列法、Kalman滤波和小波变换在地震电离层异常分析中的应用, 为研究震前电离层异常扰动提供了新的思路; 杨可可等[11]利用IGS数据分析中心提供的全球电离层格网数据,用双线性插值法和滑动四分位法对日本九州岛地震上空电离层变化特性进行分析,探讨地震对电离层变化情况的影响。

2017-08-08T21:19:46(年积日220 d, 协调世界时UT 13:19:46)九寨沟(33.20°N, 103.82°E)发生MS7.0级地震,震源深度20 km,造成重大生命财产损失。本文以此次地震为背景,利用GIM数据分析地震前后电离层异常扰动,充分考虑太阳活动、地磁变化等因素，分析电离层异常与地震的关联性。

1 数据和方法

本文采用CODE(The Center for Orbit Determination in Europe)提供的GIM(global ionosphere maps)数据(ftp://cddis.gsfc.nasa.gov/)进行电离层异常研究, 该数据时间分辨率为1 h, 空间分辨率为2.5°×5°, 单位为TECU(1 TECU=1016个电子/m2)。CODE位于瑞士伯尼尔大学, 是IGS(International GNSS Service)的7个分析中心之一, 其发布的GIM数据经过一定的插值及平滑处理,过滤掉了小范围的电离层扰动,非常适合研究全球大尺度的震前电离层异常扰动。本次九寨沟地震震级高、强度大,根据以往的研究成果分析，其很可能引起了大尺度的电离层扰动。在CODE制作GIM的过程中,使用了武汉、拉萨、乌鲁木齐、昆明、上海等GPS观测站的数据,因此可以采用CODE的GIM数据分析九寨沟地震期间的电离层异常扰动现象。在日地环境方面,本文采用中科院空间环境研究预报中心(space environment prediction center, SEPC)发布的赤道地区地磁活动指数(Dst)、 全球地磁活动指数(Kp)和10.7 cm射电流量(F10.7)( http://www.sepc.ac.cn/)反映太阳活动和地磁变化情况。

在异常检测之前,采用“db4”正交小波和软阈值函数对GIM数据进行小波阈值去噪处理,去噪的原理及具体过程见文献[12-13]。小波阈值去噪可以很好地保护有用的信号尖峰和突变信号,并且有效去除暂态信号和瞬态信号引起的噪声。经小波阈值去噪后,有效地去除了地震异常检测中的噪声干扰。

本文采用滑动四分位法[14-15]进行震前电离层的异常检测。在时窗长度的选取方面,考虑到TEC变化具有季节效应,时窗长度不宜过长,且本文主要研究震前的电离层TEC异常现象,因此时窗长度选择15天。选取年积日第205～223天共19 d的TEC数据作为待检测数据,数据包含震前15天和震后3天,滑动四分位法检测TEC异常的具体过程如下:

将某一时刻待测数据及其前15天同一时刻数据按升序排列,得到I1、I2、 …、I16, 计算其四分位数Q1、Q2、Q3[11]

(1)

四分位(IQR)为

IQR=Q3-Q1,

(2)

该时刻TEC上限UP和下限LOW可表示为

(3)

其中k为倍常数, 决定筛选异常的阈值范围。本文k值取1.5,即取1.5倍的IQR作为异常筛选的上下限。 假设待检测数据服从正态分布, UP和LOW用均值(μ)和标准差(δ)可表示为μ±2δ,即与平均值的偏差超过2倍标准差的值被认定为异常值,置信度为95%[16]。若某一时刻待检测数据高于上限或低于下限,则表明该时刻TEC出现正异常或负异常。

2 电离层异常扰动分析

2.1 日地环境分析

太阳活动、地磁变化等诸多因素都会引起电离层的异常扰动,因此,在研究地震与电离层异常的关系时首先要分析孕震区地震前后日地环境变化。本文采用滑动四分位法对年积日第205～223天的Dst、Kp和F10.7进行了异常检测(倍常数k=1.5),Dst、Kp、F10.7变化情况及异常检测结果如图1所示。

图1 Dst、Kp、F10.7变化情况(a)及异常探测结果(b)Fig.1 Variation(a) and abnormal detection results(b) of Dst,Kp,F10.7 index

结果显示,年积日第216天Dst和Kp均出现正异常现象,说明该日地磁活动强烈且存在磁暴现象,若该日出现电离层TEC异常扰动现象需要特别关注。除年积日第216天以外,此次地震前后太阳及地磁活动平静,日地环境平稳。

2.2 TEC时间序列异常分析

本文选取震中(33.20°N, 10.3.82°E)作为检测点, 震中TEC采用双线性插值法内插获得。 同时, 选取同一纬度的点A(32.5°N, 100°E)和点B(32.5°N, 110°E)、 同一经度的点C(35°N, 105°E)和点D(30°N, 105°E)作为参考点, 如图2所示,

图2 地震信息分布Fig.2 Seismic information distribution

孕震区半径由公式R=100.43M计算得出(R为半径,M为震级)。利用滑动四分位法对震中及各参考点TEC时间序列进行检测, 并绘制异常分布图, 结果如图3、 图4所示, 其中黑色纵线为震发时刻。

图3 震中TEC时间序列及异常分布异常Fig.3 Epicenter TEC time series and anomaly distribution

震中在年积日第216和218天检测到TEC正异常现象, 异常持续时间分别为7和5 h, 异常时刻集中在(UT)00:00—08:00, 即当地时间(LT)08:00—16:00, 其中, 第216天异常峰值为2 TECU, 第218天异常峰值较大, 达6 TECU。 在年积日第208天, 震中检测到TEC负异常现象, 异常值为-2 TECU, 异常时刻为(UT)23:00, 即当地时间(LT)次日的07:00。 整体来看, 震中在震前有2天出现TEC正异常, 有1天出现TEC负异常, 正异常总时长为12 h, 负异常总时长为1 h; 震后未发现TEC异常现象。

分析图4,参考点的检测结果与震中存在高度的一致性。年积日第218天,参考点A和参考点D发现的TEC正异常规模大于其他检测点; 年积日第216天, 震中和参考点C发现的TEC正异常峰值较大; 年积日第208天, 参考点C发现的TEC负异常达-4 TECU,参考点A和参考点D负异常只有-0.1 TECU左右。因此,震前出现的TEC异常区域不垂直对应震中地区。

图4 参考点TEC时间序列及异常分布Fig.4 TEC time series and anomaly distribution of reference point

综上所述, 孕震区在当地时间2017年7月28日(年积日第208天)、 8月4日(年积日第216天)和8月6日(年积日第218天)检测到TEC异常现象,这些异常现象是否由孕震引起还需结合日地环境和异常空间分布进一步分析。

2.3 全球TEC异常分析

利用滑动四分位法只能进行单点的TEC异常检测,不能确定异常的形态与分布。为了进一步判断电离层异常与九寨沟地震的相关性,本文对GIM数据的每个格网点进行异常检测,并制得当地时间2017年7月28日、8月4日和8月6日的全球TEC异常分布图,如图5—7所示,其中“★”表示震中,圆圈代表孕震区,横线代表赤道。

图5给出了2017年7月28日间隔2 h的全球TEC异常分布,可以看出,该日全球TEC异常主要为负异常且集中在赤道附近(15°S—15°N),并沿赤道自东向西移动,该现象属于赤道异常,与孕震关系不大。孕震区周围在LT06:00开始出现TEC负异常, 异常位于孕震区西南方向, 随后异常逐渐扩大、 向西移动并靠近赤道, 至LT10:00异常消失, 负异常最大约-2 TECU。 结合该日日地环境平稳,且其他区域未出现类似现象,因此,该日孕震区周围出现的电离层TEC负异常现象可能是孕育地震引起,异常区域没有垂直对应于震中,而是位于震中的西南方向并逐渐向赤道漂移。

图5 2017年7月28日全球TEC异常分布图Fig.5 Global TEC anomaly maps on July 28,2017

图6给出了2017年8月4日间隔2 h的全球TEC异常分布图, 该日全球出现大范围的TEC正异常现象,异常在全球范围内呈带状分布,南半球异常的规模、幅度及持续时间均大于北半球。结合前文对日地环境的分析,该日地磁活动强烈且存在磁暴现象,因此,该日电离层TEC正异常现象可能是地磁扰动引起。虽然孕震区域也出现了长达8 h的TEC正异常,但是该异常现象是全球性的,在其他区域也有出现,并不能作为地震前兆。

图6 2017年8月4日全球TEC异常分布图Fig.6 Global TEC anomaly maps on August 4,2017

图7给出了2017年8月6日间隔2 h的全球TEC异常分布图。除了赤道异常以外,北半球只有孕震区周围在该日出现了长时间的TEC正异常现象,异常从LT04:00持续到LT14:00,最大异常值达10 TECU,最大异常出现在LT14:00,位于震中西南方,异常区域由震中的东北方逐渐向西南方移动,表现出向赤道漂移的趋势。同时,孕震区对应的赤道共轭区内也发现明显的TEC正异常现象,赤道共轭区的异常范围大于孕震区的异常范围,异常区域也表现出向赤道漂移的趋势。考虑到该日日地环境平稳,所以孕震区周围出现的电离层TEC正异常现象可能是孕育地震引起,异常区域及异常峰值没有垂直对应于震中,异常区域自东北方向穿过震中逐渐向赤道漂移且具有共轭结构。

图7 2017年8月6日全球TEC异常分布图Fig.7 Global TEC anomaly maps on August 6,2017

3 结 论

本文基于CODE提供的GIM数据, 采用滑动四分位法作为异常检测方法, 分析了2017年8月8日九寨沟MS7.0级地震前后电离层的异常情况。 通过对震中及周围4个参考点TEC时间序列的异常检测, 发现地震孕育区内, 2017年7月28日出现电离层TEC负异常现象, 2017年8月4日和8月6日出现电离层TEC正异常现象。 结合日地环境和全球TEC异常分布进一步分析后得出: 2017年7月28日出现的电离层TEC负异常现象和8月6日出现的电离层TEC正异常现象可能是本次地震的前兆,其中7月28日异常出现在LT06:00—10:00, 异常区域没有垂直对应于震中, 而是位于震中的西南方向并逐渐向赤道漂移, 8月6日异常出现在LT04:00—14:00,异常区域没有垂直对应于震中,异常区域自东北方向穿过震中逐渐向赤道漂移且具有共轭结构;而2017年8月4日出现的电离层TEC正异常现象可能是由地磁扰动引起。

地震电离层异常的耦合机制是个非常复杂的物理和化学过程,目前,还没有很好的物理模型对其进行解释。现阶段地震电离层异常扰动的研究仍处于探索阶段,要将震前电离层异常扰动做为地震预报的一种手段还有许多问题亟待解决。本文的研究成果再次证明了强震前确实存在电离层异常扰动,为研究地震电离层异常的耦合机制提供了参考和依据。