引用本文
宋冬梅, 向亮, 单新建, 尹京苑, 王斌, 崔建勇. 基于SVR模型的电离层TEC背景场构建方法[J]. 地震地质, 2019,41(6): 1511-1528.
SONG Dong-mei, XIANG Liang, SHAN Xin-jian, YIN Jing-yuan, WANG Bin, CUI Jian-yong. THE METHOD OF CONSTRUCTING IONOSPHERIC TEC BACKGROUND FIELD BASED ON SVR MODEL[J]. SEISMOLOGY AND GEOLOGY, 2019,41(6): 1511-1528.  
DOI:10.3969/j.issn.0253-4967.2019.06.013
Permissions
Copyright©2019, 《地震地质》编辑部
《地震地质》编辑部 所有
基于SVR模型的电离层TEC背景场构建方法
宋冬梅1,2), 向亮3), 单新建4), 尹京苑5), 王斌1), 崔建勇1)
1)中国石油大学(华东), 海洋与空间信息学院, 青岛 266580
2)青岛海洋科学与技术试点国家实验室, 海洋矿产资源评价与探测技术功能实验室, 青岛 266071
3)中国科学院海洋研究所, 青岛 266071
4)中国地震局地质研究所, 地震动力学国家重点实验室, 北京 100029
5)上海市地震局, 上海 200062

〔作者简介〕 宋冬梅, 女, 1973年生, 教授, 主要从事灾害遥感研究, E-mail: songdongmei@upc.edu.com。

收稿日期: 2018-12-05
基金项目: 国家自然科学基金(41772350)、 上海市科学与技术委员会科研计划项目(16dz1206000)和地震动力学国家重点实验室开放基金(LED2012B02)共同资助
摘要

引起电离层电子浓度总含量(Total Election Content, TEC)变化的因素有很多, 而地震TEC扰动只是很小的一部分。 文中尝试构建一个考虑了太阳活动与地磁活动影响的TEC非震动态背景场, 对比分析TEC非震动态背景场和传统的滑动时窗背景相对于原始TEC值的残差情况, 滑动时窗背景的残差存在明显的月周期与半年周期, 该周期性残差将对后续的地震电离层异常探测造成重要影响。 同时, 利用非震动态背景场法与滑动时窗法探测了汶川附近研究点(30°N, 100°E)2008年3月1日—9月26日长时间序列的TEC异常情况, 发现在太阳活动、 地磁活动等非震干扰的情况下, 相较于传统的滑动时窗法, TEC非震动态背景场法很少检测出TEC异常扰动; 而在地震发生前, 基于文中的方法提取的TEC异常较滑动时窗法提取的异常强度更大、 异常次数更多。 最后, 文中分析了汶川地区2008年5月12日、 8月21日、 8月30日3次地震前的TEC异常表现, 均主要为负异常扰动, 且异常主要分布在震中靠近赤道一侧。

关键词: 电离层; 电离层电子总量; SVR模型; 小波分析; 汶川地震
中图分类号:P315.72+1 文献标志码:A 文章编号:0253-4967(2019)06-1511-18
THE METHOD OF CONSTRUCTING IONOSPHERIC TEC BACKGROUND FIELD BASED ON SVR MODEL
SONG Dong-mei1,2), XIANG Liang3), SHAN Xin-jian4), YIN Jing-yuan5), WANG Bin1), CUI Jian-yong1)
1)College of Ocean and Space Information, China University of Petroleum(East China), Qingdao 266580, China
2)The Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China
3)Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China
4)State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China
5)Shanghai Earthquake Agency, Shanghai 200062, China
Abstract

There are many factors related to the variations of TEC, and the changes of TEC caused by earthquake only occupy a small portion. Therefore, it is vital how to exclude the ionospheric interference of non-seismic factors accurately in the process of seismic ionospheric anomaly extraction. This study constructed a TEC non-seismic dynamic background field considering the influence of solar and geomagnetic activities. Firstly, the TEC components of half-year cycle and annual cycle are extracted by wavelet decomposition. Then, it establishes a regression model between TEC in which periodic factors are removed and solar activity index, geomagnetic activity index with SVR method(support vector regression)in non-seismic period. Finally, based on the constructed model, the solar activity index and geomagnetic activity index is used to reconstruct aperiodic components of TEC in earthquake's period. From the reconstructed aperiodic components of TEC plus the half-year periodic components and annual periodic components of TEC in the same period, the non-seismic dynamic background field is obtained. Comparing the residuals relative to original TEC values in non-seismic dynamic background field and traditional sliding window background, there are apparent monthly periodic change and semi-annual periodic change in the residuals of sliding window background, which can have obvious impacts on the subsequent seismic ionospheric anomaly detection. In order to test the validity of seismic TEC anomaly detection based on the background field construction method, this paper investigated the long time series TEC anomalies near Wenchuan city(30°N, 100°E)from March 1 to September 26 in 2008. It is found that under the condition of non-seismic disturbance such as solar activity and geomagnetic activity, TEC abnormal disturbance is rarely detected by non-seismic dynamic background field method, when compared with the traditional sliding time-window method. And before the earthquake, more TEC anomalies were detected based on the proposed method, also, they were more intense than those extracted by sliding window method. Therefore, the TEC background field construction method based on SVR(support vector regression)has superiorities in both system errors elimination, which are caused by solar, geomagnetism, the non-seismic ionospheric disturbance events and periodic fluctuations of TEC, and in reducing the false alarm rate of seismic TEC anomaly. Moreover, it can also improve the seismic TEC anomaly detection ability. In addition, this paper analyzed the time-spatial distribution of TEC anomaly before three earthquakes on May 12, August 21 and August 30, 2008. They were mainly negative abnormal perturbations and often distributed on the equatorial side of epicenter.

Keyword: ionosphere; TEC(Total Electron Content); SVR(Support Vector Regression); wavelet analysis; Wenchuan earthquake
0 引言

自1964年Leonard等(1965)利用离子探空仪在美国阿拉斯加大地震前发现电离层行进式扰动现象后, 地震电离层前兆研究逐渐被人关注, 吸引了大量研究人员投入地震电离层前兆研究工作中, 并在许多大地震前发现了电离层异常前兆现象(Pulinets et al., 2005; Zhao et al., 2008)。 近20a来, 随着电离层遥感技术的发展, 电离层监测能力逐渐提升。 地震电离层前兆探测, 尤其是针对电离层电子浓度总量(TEC)的研究成为地震电离层前兆研究的热点之一。 Liu等(2001)对1999年7.7级集集地震的TEC进行时空分析, 发现在震前4d内震中附近多次出现TEC值明显降低的现象, 且异常主要分布在震中靠近赤道一侧。 Aggarwal(2015)对2010年青海6.9级地震前多个IGS站点的数据开展调查, 发现在震前第4天出现低值异常, 并且异常从距离震中较远的站点开始逐渐向震中附近移动。 Ho等(2013)通过对比TEC和DEMETER卫星的电子密度数据, 发现在2010年智利8.8级地震前9~19d夜间震中附近均出现高值异常, 但2006— 2009年同时期数据并未出现类似的正异常现象。 Le等(2011)针对TEC数据对2002— 2010年736个6级以上地震进行统计分析, 结果显示TEC异常的发生率随着地震临近而升高, 并且大震级、 浅源地震更加明显。

尽管在很多地震前均发现了电离层的TEC异常现象, 但关于地震电离层TEC的研究仍然存在一些争议。 Afraimovich(2004)对1999年加利福尼亚7.1级地震震前的TEC展开调查, 分析认为观测到的TEC变化似乎受到当地时间和地磁活动的控制, 而不与即将到来的地震有任何关联。 Thomas等(2017)对2000— 2014年1 279个6级以上地震前的TEC展开统计分析, 结果显示TEC异常与即将到来的地震在统计意义上并没有关联性。

地球电离层作为近地空间环境的重要组成部分, 其状态主要由不同波长范围内的太阳辐射通量所决定(Balan et al., 1996; Afraimovich et al., 2008; Liu et al., 2011), 同样地, TEC值与太阳活动高度相关。 此外, TEC波动还受到极端天气(张雅雯等, 2013)、 天文现象(Reinisch et al., 2018)、 地磁环境(Belehaki et al., 2017)等诸多因素的干扰, 而地震电离层TEC异常探测的关键在于如何从众多形式的TEC波动中确定出与地震有关的扰动。 目前的地震电离层异常研究主要是利用滑动时窗法对震前、 震后十多d的TEC进行异常检测, 虽然有一定的合理性, 但研究时段较短, 且存在以下不足: 1)由于背景场的构建过程中未考虑TEC自身周期性波动的影响, 导致提取的异常信号里仍然存在周期性信息; 2)由于背景场中未考虑太阳活动等非震因素对电离层TEC扰动的影响, 使得基于此背景场方法的异常检测在平静期时也能提取出大量TEC异常, 即存在较高的虚警率; 3)在地震事件和其它非震扰动事件同时发生时, 无法判别TEC异常的来源。

鉴于已有方法的不足, 本文提出了一种新的TEC非震动态背景场构建方法, 该方法考虑了太阳活动、 地磁活动等空间环境扰动以及TEC自身的周期性波动。 同时, 还分别利用本文所述方法与传统的滑动时窗法提取2008年3月1日— 9月26日期间汶川附近研究点(30° N, 100° E)的TEC异常信息并进行比较分析。

1 数据情况
1.1 全球TEC地图数据

本文主要选用欧洲定轨中心(Center of Orbit Determination in Europe, CODE)发布的全球TEC地图数据①(①http://cddis.gsfc.nasa.gov/pub/gps/produvts/ionex。) (Global Ionosphere Map, GIM)进行研究, 其基于全球200多个GPS地面接收站的数据插值而成, 时间分辨率为2h, 每天从0时到24时共生成12幅地图。 其经度覆盖范围为180° W— 180° E, 经向空间分辨率为5° ; 纬度覆盖范围为87.5° N— 87.5° S, 纬向空间分辨率为2.5° , 因此, 每张地图包括5i183(71× 73)个格网点。 TEC的单位为TECU, 1TECU表示1m2有1016 个电子。 由于TEC格网地图具有空间上固定、 时间上连续的全球观测能力, 目前已被广泛应用于电离层特性分析(Guo et al., 2015; Klimenko et al., 2016)以及地震电离层前兆研究中(Tao et al., 2017), 特别是电离层TEC异常与地震的统计分析中(Shah et al., 2015; Zhu et al., 2016)。

1.2 太阳和地磁活动数据

理想情况下, EUV辐射通量数据为衡量不同波长范围太阳活动的最佳指标(Bilitza, 2001), 但该数据无法直接在地面观测, 只能利用卫星紫外线仪获得, 且无法获得稳定连续的观测数据。 因此, 为了分析太阳活动与TEC长周期特性的关系并构建两者之间的模型, 本文选择太阳黑子数(SSN)和太阳10.7cm波长的辐射通量①(① https ://omniweb.gsfc.nasa.gov/form/dx1.html。) (F10.7指数)作为衡量太阳活动的指标, 两者的时间分辨率为1d, 均可在地面观测到, 并且可以获得长时期的数据资料。

TEC波动除了受太阳活动的影响外, 还受到地磁活动与宇宙环境的影响, 大量的报道记录了这些因素造成的TEC扰动现象(Orus et al., 2007; 闫相相等, 2014)。 本文在非震动态背景场的构建过程中充分考虑了这些因素, 主要选择KP指数(每日3h内的地磁扰动总强度)、 AP指数(全球全日地磁扰动指数)、 DST指数(磁赤道附近环型电流扰动)、 AE指数(极光带每h正变化与最大负变化的绝对值)和IMF-Bz数据(电离层平均磁场强度) ②(② https://omniweb.gsfc.nasa.gov/form/dx1.htm。)来衡量这些因素的影响。

1.3 地震数据

本文主要针对汶川及其周边区域进行研究, 研究范围为60° ~150° E, 35° S— 62.5° N。 从USGS(United States Geological Survey)下载该区域2004— 2011年5级及以上地震目录③(③ https://earthquake.usgs.gov/earthquakes/search/。), 用于剔除非震动态背景场构建过程中训练集内受地震影响的数据。

2 TEC非震动态背景场建立的原理与方法
2.1 TEC非震动态背景场的概念

虽然利用全球电子总量地图(GIM)可以获得空间电离层TEC全球性、 长时间序列的连续观测数据, 但引起TEC波动的因素很多, 包括太阳活动、 地磁活动、 空间环境以及极端天气等, 而地震引起的TEC扰动只是众多扰动中很弱小的一部分, 直接从GIM中提取与地震相关的信息依然很困难。 因此, 如何排除非震因素的强干扰, 提取与地震有关的电离层TEC弱扰动是地震电离层前兆研究的关键。

在地震电离层TEC扰动的研究过程中, 首先应找到地震造成的电离层TEC扰动与其它因素(太阳活动、 地磁活动、 空间环境与极端天气等)造成的TEC扰动之间的区别, 否则, 就难以从电离层TEC波动中准确提取与地震有关的扰动。 然而, 目前的研究尚不能明确把握地震造成的TEC扰动的特性。 在此情况下, 先找出主要由非震因素引起的扰动, 作为正常的背景场, 这样既可以排除非震强背景的干扰, 又可放大地震造成的TEC扰动的弱信号, 然后再从剩余的TEC信息中寻找地震造成的TEC异常扰动。

建立非震动态背景场, 需要对非震干扰因素有充分的理解。 图1a给出了汶川地区附近点(30° N, 100° E)2000— 2016年的归一化TEC日均值与归一化F10.7值的情况。 从图中可以直观地看出, TEC值与F10.7在长期趋势上存在很明显的相关性。 2007— 2009年为太阳活动低年, 归一化F10.7值与归一化TEC值存在更大偏差, 即两者并非简单的线性相关。 Afraimovich等(2006)对全球TEC和F10.7值进行了分析, 发现两者在太阳活动低年存在较好的线性相关性, 但在太阳活动高年, 随着F10.7值增加到一定程度, TEC值将不再随之改变。

图 1 归一化TEC值与归一化F10.7值
a 原始曲线; b 频谱分布Fig. 1 Normalized TEC and F10.7(located at 30° N, 100° E).

为了进一步分析太阳活动与TEC之间的关系, 对归一化后的F10.7值、 TEC值进行傅里叶变换, 分析两者的频谱图。 从图1b中可见, TEC与F10.7在各周期的幅度基本一致, 且信号能量集中在大于半年的长周期成分上, 而地震信号主要集中在短周期成分。 红色曲线显示TEC在约180d周期和360d周期存在相对较强的信号, 而F10.7在这2个周期成分上没有出现较强的信号。 180d周期与360d周期正好对应了TEC的半年周期与年周期成分, 说明TEC存在半年周期与年周期特性。 据实验观测, 去除半年周期与年周期成分后, TEC与F10.7值的吻合度更高。 因此, 基于此观测结果, 本文的方法首先剔除TEC数据的半年周期和年周期性成分, 再利用SVR模型进行TEC与非震因素的回归建模。

某一地点的电子总量T0受到太阳活动、 地磁活动等多因素的干扰, 在时间序列上表现出不断波动, 则有:

T0=T1+T2+T3+T4+Te(1)

式中, T1为太阳周期性活动的影响, 主要表现在长周期变化上, T2为TEC自身稳定的年周期及半年周期成分波动, T3为地磁和空间宇宙环境活动的影响, T4为极端天气、海啸等其它因素的TEC扰动, Te为地震引起的TEC扰动。其中, T1、T2和T3这3个部分为TEC波动的主要来源, 如果能准确模拟出这3个因素的扰动情况, 即认为找到了非震动态TEC背景场Tback。本文首先利用小波多尺度分解提取T2部分的信息, 然后根据太阳活动与地磁活动和TEC信号的相关性, 利用支持向量机回归算法(SVR)构建太阳、地磁活动指数与T1、T3 2部分信息之间的回归预测模型。最后整合3个部分的信息, 得到TEC非震动态背景场Tback:

Tback=T1+T2+T3(2)

在剩余的TEC成分Trest中, 地震TEC扰动信号得到放大。通过分析Trest, 可以更为准确地提取地震电离层TEC扰动。

Trest=T4+Te(3)

2.2 小波分析

小波变换分析方法是一种将信号转换到时间-频率域上进行分析的方法, 其克服了傅里叶变换时频不可分的局限性, 能够同时对信号的时间域与频率域进行分析, 因此被广泛地应用于地震信号去噪、地震信号时频分析、异常检测等工作中(韩冰, 2014), 其公式为

Ψa, b(t)=1a　Ψt-ba(4)

式中, , b均为常数, Ψ (a, b) (t)由Ψ (t)函数经伸缩和平移变换得到, a为伸缩因子, b为平移因子, 因此a、b的组合决定了Ψ (a, b) (t) a的形式。

本文主要利用小波分析中的一维小波多尺度分解结果, 其可以看作是一组高、低通滤波器对信号进行层叠滤波的结果, 其中高通滤波器的输出反映各尺度下的高频细节成分, 低通滤波器的输出反映信号的低频近似成分(熊攀, 2009)。图2所示的流程通过对TEC原始信号T0进行多层分解可得到不同周期的高频细节成分, 其中如图2b所示, TEC的原始信号中7阶高频细节成分正好对应着半年周期成分, TEC原始信号中8阶高频细节成分正好对应着年周期成分, 提取7阶与8阶高频细节成分, 即为公式(1)中TEC的自身周期成分T2, 而TEC原始信号减去TEC自身周期成分T2得到TEC的残差Tr。

图 2 TEC小波多尺度分解
a 小波分解及Tr的计算过程流程图; b 小波分解的高频成分曲线。T0为原始的TEC时间序列, H1、H2、H3、H4、H5、H6、H7和H8为小波分解的18阶高频成分, L1、L2、L3、L4、L5、L6、L7和L8为小波分解的18阶低频成分, T2为TEC自身年周期及半年周期成分波动, Tr为剔除稳定周期成分后的TEC残差Fig. 2 Wavelet multiscale decomposition of TEC.

2.3 支持向量机回归模型

支持向量机(Support Vector Machine, SVM)最早是基于统计学习理论与结果风险最小化理论下发展起来的一种处理二分类问题的方法。Vapnik(1999)引入了ε 不敏感损失函数, 将SVM推广到非线性回归估计中, 得到支持向量回归算法(Support Vector Regression, SVR)。SVR的主要思想是通过核函数将线性不可分的样本数据T= x1, y1, ……, xl, yl∈ x×yl映射到线性高维特征空间, 然后通过解决线性最优化问题得到回归函数f(x)=ω · φ (x)+b, 其中φ (x)为非线性映射函数。本文采用的ε 不敏感损失函数和核函数的支持向量回归数学形式为

minω, ε=12ω2+Cl∑i=1lξi+ξi* (5)

其中, C为惩罚参数, ξi、ξi* 皆为模型的松弛变量。本文选用不敏感损失函数的目的是当实际值和预测值的误差不超过提前设定的ε 时, 则认为在该时间点的预测不存在损失。

SVR对于解决小样本、非线性、多维数据的回归预测问题具有许多独有的优势。由前文分析可知, 地磁活动、太阳活动与TEC存在较大的非线性相关性, 且衡量太阳活动与电磁活动的相关指标较多, SVR正好适用于这种情况。本文主要利用SVR构建太阳活动、地磁活动指数与Tr(剔除年变周期、半年周期成分的TEC)之间的回归模型。模型输入量为表征太阳活动和地磁活动的9个变量, 见表1, 模型输出量为Tr。

表1 SVR模型输入参量 Table1 Introduction to SVR model input parameters

不同地理位置、 不同时刻的TEC变化存在差异, 对太阳活动等因素的响应也不同。 为了减少建模过程中地理位置与时刻带来的误差, 本文将研究区(35° S— 62.5° N, 60° ~150° E)内2004— 2011年的数据按不同地理位置(共19× 40=760个网格)、 不同时刻(共12个时刻)分成不同的数据集(即将同一位置、 同一时刻的数据放到同一个数据集), 则研究区总共有12× 760=9 120个数据集, 每个数据集有8(a)× 365(d)=2 920个数据。 然后针对每个数据集分别建立TEC观测值Tr与太阳活动、 地磁活动等空间环境变量之间的SVR回归模型。

SVR回归建模的具体流程如图3所示。在回归建模之前, 首先需要对每个样本集的原始TEC残差Tr与不同时间分辨率的太阳、地磁活动指数进行匹配。然后, 为了保证模型训练的输入、输出均为未受地震干扰的数据, 需要对每个数据集进行非震样本筛选, 即选择不受地震干扰的数据输入模型进行训练。筛选原则为:将时间上前后15d、空间上R=100.43M(R为孕震区范围, M为地震的震级)范围内出现5级以上地震的数据剔除(Dobrovolsky et al., 1979)。在进行SVR回归建模之前, 需要对回归建模数据指标进行标准化处理, 这里选择的是最大最小0-1归一化方法(汤荣志等, 2016)。最后, 对归一化后的指标数据, 随机选取数据集中的80%作为SVR模型的训练集, 剩下的20%作为SVR模型的测试集。SVR模型选用RBF核函数, 输出为模型参数、模型的平均误差σ (即TEC的预测值T'r与实际Tr之间的差值)及模型精度。SVR建模过程中先将训练集数据输入SVR回归模型, 其中输入为表征太阳活动和地磁活动的9个变量, 见表1, 模型的输出为Tr。通过不断调整模型参数进行训练, 找到能刻画输入、 输出数据回归关系的最优模型, 并保存最优SVR模型的参数。 按照每个时刻(每天12个时刻)、 每个地理空间位置(共760个网格)分别进行训练与建模, 因此总的模型数量为12× 760=9 120个。然后将测试集数据中表征太阳活动和地磁活动的9个变量输入最优SVR模型得到预测的TEC残差T'r, T'r即表示太阳活动等非震因素对TEC值影响的部分。最后计算出模型预测的T'r与测试集原始Tr之间的误差和精度, 即为模型的精度与模型平均误差σ 。

图 3 SVR建模流程图Fig. 3 The flow chart of SVR modeling.

2.4 TEC非震动态背景场的建立与分析

如图 4 所示, 对于每一个时刻(每2h)、 每一地理位置(经纬度网格为5° × 2.5° ), 将其对应的表征太阳活动和地磁活动的9个变量输入到已经训练好的SVR模型, 输出该时间点、 该位置空间环境中TEC的预测值, 然后加上该时间点、 该位置的原始信号T0经小波分解得到的小波7阶高频值H7和8阶高频值H8, 即得到该时刻、 该空间位置的TEC非震动态背景值Tback。 根据上述过程, 分别计算出研究区所有时刻(以2h为间隔, 每天共12个时刻)、 所有地理位置(经纬度网格为5° × 2.5° )的非震动态背景值, 得到整个研究区的TEC非震动态背景场。

图 4 TEC非震动态背景场的构建流程图Fig. 4 The establishment process of TEC dynamic non-seismic background field.

图 5 为研究区2004— 2011年TEC非震动态背景场的平均精度。 由图 5 可知, 研究区TEC非震动态背景场的平均精度 > 80% , 其中低纬度地区的精度高于高纬度地区, 且在20° N附近均存在1条相较于两侧区域精度更低的条带。 此外, 在(60° N, 100° E)附近出现1片相对于整个研究区而言精度极低的区域, 汶川地区(30° N, 100° E)的精度约87%。

图 5 2004— 2011年TEC非震背景场的平均精度图Fig. 5 Average accuracy of TEC non-seismic background field from 2004 to 2011.

为了分析本文方法与传统的滑动时窗法背景场的差异, 分别对汶川地区研究点(30° N, 100° E)的TEC非震动态背景场时间序列残差和传统的滑动时窗法背景(测试集每个数据选择其前15d对应时刻数据的均值作为该数据的背景)时间序列残差进行傅里叶变换, 结果如图 6 所示。 在30d周期与180d周期上, 蓝色曲线出现较高的峰值, 幅值高于其它周期2倍以上, 分别达到1.2TECU与1.6TECU, 说明滑动时窗法的残差中还存在月周期与半年周期的变化成分, 而这些非震因素引起的波动将对地震电离层TEC异常的探测造成巨大影响。 红色曲线各周期成分的幅值分布相对较为平均, 基本在0.8TECU以下, 并没有出现某个周期信号特别强烈的情况。 因此, 本文提出的非震动态背景场能更有效地去除TEC残差中的周期性系统误差。

图 6 TEC背景残差傅里叶频谱
红色曲线表示TEC非震动态背景残差, 蓝色曲线表示滑动时窗背景残差Fig. 6 The Fourier spectrum of TEC background residuals.

3 汶川地区TEC长时间序列异常探测的实例分析
3.1 数据选取及异常探测策略

本文选取汶川附近研究点(30° N, 100° E)作为异常探测对象, 分别使用非震动态背景场与滑动时窗背景来提取2008年第61~270天近7个月的TEC异常情况。 其中基于非震动态背景场的TEC异常提取方法以TEC非震动态背景场为背景值, 以1.75倍数据集平均残差作为阈值偏差, 当TEC值与非震动态背景场的偏差超过1.75倍数据集平均残差时, 即认为该数据为异常数据。 滑动时窗法主要是对于每个时刻、 每个位置, 取前15d的数据均值作为背景值, 以前15d数据的1.75倍标准差作为阈值偏差来探测异常。 在探测时段中, 研究区附近发生了5次地震, 如表2 所示, 其中2008年5月12日发生的汶川地震距离研究点最近, 且震级最大。

表2 2008年3— 9月汶川附近5次MW6.0地震参数 Table2 The introduction to 5 earthquakes with MW≥ 6.0 happening near Wenchuan from March to September in 2008

此外, 本文还选取研究时段内太阳活动与地磁活动指数数据, 来分析基于不同方法探测得到的异常与太阳活动、 地磁活动的关系。 图 7 展示了2008年3月1日— 9月28日(DOY61— DOY270)共210d的太阳与地磁活动指数的变化情况, 其中横坐标DOY表示一年中的第几天。 从F10.7指数变化(图6b)可以看出, 在第84~97天出现了1个较大的峰值, 反映了该时段太阳活动的增强极可能导致TEC值上升。 此外, 为了分析地磁与空间环境对电离层的干扰影响, 当DST指数< -40nT、 KP指数> 40、 IMF指数> 10nT、 AP指数> 40、 AE指数> 600时, 即认为该时刻电离层发生了非震扰动事件。 从图8a、 c、 d、 e和f可看出, 在第61、 68、 69、 86、 87、 88、 96、 97、 122、 115、 142、 167、 168、 194、 205、 223、 231、 248和259天共19d出现非震电离层扰动事件。

图 7 2008年3月1日— 9月26日(DOY61— DOY270)太阳与地磁活动指数的变化情况
a— f依次为IMF值、 F10.7指数、 KP指数、 AP指数、 DST指数和AE指数; 红色竖线表示电离层扰动事件Fig. 7 Variations of the sun and geomagnetic activity index(DOY61— DOY270)from March 1 to September 26, 2008.

图 8 2008年3月1日— 9月26日基于TEC非震动态背景场法的TEC异常探测结果Fig. 8 Detection results of TEC anomaly based on non-seismic background field method from March 1 to September 26 in 2008.

3.2 异常探测结果分析

本文分别用基于非震动态背景场法与滑动时窗法对汶川地区研究点(30° N, 100° E)2008年3— 9月的TEC时间序列进行异常探测, 结果如图 8 与图 9 所示。 其中红色曲线为TEC原始曲线, 绿色曲线为异常探测的上边界, 蓝色曲线为异常探测的下边界, 当原始值超出异常探测的上边界即认为该时刻出现正异常扰动(以红色柱体进行标记), 当原始值低于异常探测的下边界即认为该时刻出现负异常扰动(以蓝色柱体进行标记)。 为了更为明显地观测到异常扰动情况, 图 8 与图 9 中的柱状图显示的正负异常扰动为放大3倍后的结果。 此外, 图中黑色虚线表示地震事件, 黑色圆点表示电离层非震扰动事件。

图 9 2008年3月1日— 9月26日滑动时窗法的TEC异常探测结果
图中红色曲线、 绿色曲线及蓝色曲线分别代表原始TEC值、 TEC异常探测的上界限和下界限; 红色柱体为正异常扰动, 蓝色柱体为负异常扰动, 黑色虚线代表地震事件, 黑色圆点为电离层非震干扰事件Fig. 9 Detection results of TEC anomaly based on sliding time window method from March 1 to September 26 in 2008.

由图 8 与图 9 可知, 在整个研究时段内(2008年3月1日— 2008年9月26日)基于本文方法提取的异常次数明显少于滑动时窗法提取到的异常次数。 结合电离层非震扰动情况, 在第86~88天、 第167~168天、 第223天与第248天4次较强的电离层非震干扰事件中, 利用滑动时窗法均探测到较强的正异常扰动, 且整个探测时段正异常扰动出现的时间基本相隔约27d, 这可能是图 6 中背景残差27d周期成分的影响。 而基于本文方法提取的异常较弱或几乎没有, 说明本文方法能够很好地降低由非震异常影响造成的虚警率。 闫相相等(2014)通过滑动四分位法探测中国西南区域的TEC异常并分析发现异常出现的时间基本相隔约27d, 与太阳自转周期27d一致。

结合5次震例事件(如图 8 和图 9 中的虚线所示), 分析图 8 与图 9 2种方法提取的TEC异常情况。 A地震前, 利用本文方法在地震前第2~6天(DOY127~131, 即5月6日— 5月10日)探测到多次负异常扰动, 其中地震前第3天(DOY130, 即5月9日)探测到正异常扰动。 而利用滑动时窗法只在地震前第6天(DOY127, 即5月6日)和第3天(DOY130, 即5月9日)探测到负异常扰动, 在震前第3天(DOY130, 即5月9日)同样还探测到较强的正异常扰动。 2种方法在汶川地震前均探测到不同程度的负异常扰动与正异常扰动, 与前人的异常探测结果基本一致(Liu et al., 2009; 闫相相等, 2014), 但相比于滑动时窗法, 本文的方法检测到的异常强度相对较高, 且负异常探测次数相对较多。

B地震前, 利用本文方法在地震前6d内没有探测到较为明显的异常现象, 而利用滑动时窗法在地震前第4天(DOY142, 即5月21日)发现明显的正异常, 但当天同时也发生了电离层扰动事件, 且B地震距离研究对象点(30° N, 100° E)的距离为659km, 大于其地震孕震范围的半径420km(根据R=100.43M计算), 因此并不能肯定5月21日的正异常扰动与其后的地震有关。

C地震前, 利用本文方法在地震前6d内没有探测到较为明显的异常现象, 而利用滑动时窗法在地震前第6天(DOY212, 即7月30日)和第2天(DOY216, 即8月3日)发现较弱的正异常扰动。 C地震距离研究对象点(30° N, 100° E)为675km, 大于其地震孕震范围380km(根据R=100.43M计算), 因此该次地震前的2次正异常扰动是否与地震有关, 尚不能完全肯定。

D地震前, 利用本文方法与滑动时窗法在震前第2天(DOY232, 即8月19日)均探测到较弱的负异常扰动。 E地震的异常探测结果与D地震基本相似, 即在震前第2天(DOY241, 即8月28日)2种方法均探测到较弱的负异常扰动。

综上所述, 2种方法探测到的震前异常扰动基本上均为负异常扰动, 而地磁活动的影响主要为正异常扰动。 为了进一步分析负异常扰动与地震的关系, 对2008年3月1日— 9月26日共271d探测到的负异常采用10d时间窗滑动统计, 结果如图 10 所示。 图10a显示了2008年3月1日— 9月26日基于本文方法得到的负异常次数的滑动统计结果, 可以明显地看到在5月16日(DOY136)出现1个明显的峰值, 即说明该天的前10d(DOY127~136)负异常活动出现较频繁, 达到13次, 远大于其它天数的负异常滑动统计结果, 此时段频繁的负异常扰动可能与5月12日(DOY133)的汶川7.9级地震(A地震)有关。 图10b中滑动时窗法负异常叠加统计的结果显示负异常次数呈现出月周期的波动特征, 且整体上累加统计的负异常次数高于图10a中的负异常次数, 这说明滑动时窗法的负异常探测结果中可能还存在其它非震因素的干扰和系统误差的扰动。 为了对比非震年份与地震年份异常扰动的差异, 本文还对2007年同时段的负异常次数进行了10d窗方法的滑动统计, 结果如图10c所示。 2007年负异常10d窗方法的滑动次数基本不超过5次, 与2008年非震时期的结果相似, 说明利用本文方法在非震年份或非震时期检测到的负异常次数较少, 即虚警率低。

图 10 负异常次数10d窗滑动统计结果
a 2008年3月1日— 9月26日基于本文方法提取的TEC异常结果; b 2008年3月1日— 9月26日基于滑动时窗法提取的TEC异常结果; c 2007年3月1日— 9月26日基于本文方法提取的TEC异常结果Fig. 10 The negative abnormal statistical results.

图 11— 13显示了基于本文方法提取出的研究区范围内(62.5° N— 35° S, 60° ~150° E)震前TEC异常的空间分布。 图 11 显示了2008年5月6— 10日TEC扰动的空间分布情况, 其中A为2008年5月12日汶川地震的震中位置、 F为2008年5月7日日本地震的震中位置。 由图可见, 从5月6日06UT(UT表示世界时)开始, A地震附近逐渐出现负异常扰动, 并不断增大, 到08UT时异常达到最大值, 且相对震中向赤道一侧偏移, 异常整体呈椭圆状, 经向范围为90° ~120° E, 纬向范围为20° ~35° N。 同时, 在南半球磁共轭区域也出现了负异常, 不过相对于震中附近的异常较弱。 接下来几天, 5月7日08— 10UT、 8日06— 08UT、 10日04— 06UT在震中附近均出现不同程度的负异常扰动, 其空间分布形态与5月6日的异常分布基本类似, 同时在南半球磁共轭区也多次出现对应的扰动情况。 此外, 在5月6日, F地震发生的前1d出现了较弱的负异常扰动, 可能与接下来的F地震有关, 这与Liu等(2009)对汶川地震前的异常检测结果基本一致。 值得一提的是, 本文方法在5月9日同样检测出较强的正异常现象, 这一次的电离层正异常现象也一度成为众多学者讨论的热点(Zhao et al., 2008; Liu et al., 2009)。

图 11 2008年5月6日— 5月10日基于本文方法提取的TEC异常空间分布Fig. 11 Spatial distribution of TEC anomaly extracted by the proposed methods from May 6 to May 10 in 2008.

图 12 基于本文方法提取的2008年8月19日TEC扰动的空间分布Fig. 12 Spatial distribution of TEC anomaly extracted by the proposed methods on August 19, 2008.

图 13 基于本文方法提取的2008年8月28日TEC扰动的空间分布Fig. 13 Spatial distribution of TEC anomaly extracted by the proposed methods on August 28, 2008.

图 12 显示了2008年8月19日基于本文方法提取的TEC异常的空间分布, 其中D为2008年8月21日的中国与缅甸交界的6.0级地震的震中位置, G为2008年8月27日贝加尔湖区域6.3级地震的震中位置。 由图可见, 在08UT, D地震与G地震附近均出现了负异常, 且异常强度与空间展布范围远小于图 11 中探测到的负异常扰动。

图 13 显示了8月28日基于本文方法提取的TEC扰动的空间分布情况, E为2008年8月30日四川与云南交界的6.0级地震的震中位置。 在06— 08UT时段内, E地震附近出现了负异常扰动, 且负异常分布向赤道一侧偏移。 而在G区域附近出现正异常扰动, 这可能与极地电离层活动有关, 而与本地震无太大关联。

综合分析A、 D和E 3次中国区域的地震可知, 在3次地震前均检测到一定程度的负异常扰动, 而且负异常均表现为从震中偏向赤道一侧。 其中, A地震(汶川地震)前异常强度最大、 范围最广, 且出现多次负异常扰动, 这可能与A地震震级较大有关。 另外, 在A地震前负异常存在明显的关于磁力线共轭的现象, 这可能与电离层异常沿着磁力线移动到震中共轭区域有关(Pulinets et al., 2004)。

4 结论

引起电离层TEC扰动的因素有很多, 而地震引起的TEC扰动只是影响因素之一。 本文尝试构建一个考虑了太阳活动与地磁活动影响的TEC非震动态背景场, 详细介绍了TEC非震动态背景场的构建思路与流程, 并对比分析了本文方法与传统的滑动时窗背景法的TEC残差情况。 结果显示, 滑动时窗背景法的TEC残差存在明显的月周期与半年周期, 这对后续电离层异常的探测将造成重要影响。 为了检验基于本文背景场构建方法的地震TEC异常探测的效果, 以汶川地区为研究对象, 分别用本文方法与滑动时窗法探测了研究点(30° N, 100° E)2008年3月1日— 9月26日长时间序列的TEC异常情况, 得到如下结论: 在太阳活动、 地磁活动等非震因素干扰的情况下, 相较于传统的滑动时窗法, 本文的方法几乎很少检测出TEC异常扰动, 说明本文方法有效降低了虚警率; 而在地震发生前, 基于本文方法提取的异常较滑动时窗法强度更大, 异常次数更明显, 且异常通常表现为负异常扰动。 因此, 本文提出的背景场构建新方法能够很好地消除太阳、 地磁非震电离层干扰和TEC自身周期波动带来的系统误差, 同时还能提高地震TEC异常探测的能力。

参考文献
[1] 韩冰. 2014. 小波变换及其在地震电磁信号分析中的应用研究 [D]. 北京: 中国地震局地质研究所.
HAN Bing. 2014. Wavelet transformation and its application in analyzing seismic electromagnetic signal [D]. Institute of Geology, China Earthquake Administration, Beijing(in Chinese). [本文引用:1]
[2] 汤荣志, 段会川, 孙海涛. 2016. SVM训练数据归一化研究[J]. 山东师范大学学报(自然科学版), 31(4): 60—65.
TANG Rong-zhi, DUAN Hui-chuan, SUN Hai-tao. 2016. Research on data normalization for SVM training[J]. Journal of Shand ong Normal University(Natural Science), 31(4): 60—65(in Chinese). [本文引用:1]
[3] 熊攀. 2009. 小波方法在地震遥感信息提取中的应用 [D]. 北京: 中国地震局地震预测研究所.
XIONG Pan. 2009. Wavelet-based methods for detecting earthquake information in remote sensing data [D]. Institute of Earthquake Forecasting, China Earthquake Administration, Beijing(in Chinese). [本文引用:1]
[4] 闫相相, 单新建, 曹晋滨, 等. 2014. 中国西南区域孕震区电离层TEC变化长时间序列分析[J]. 地震地质, 36(1): 253—265. doi: 103969/j.issn.0253-4967.2014.01.02.
YAN Xiang-xiang, SHAN Xin-jian, CAO Jin-bin, et al. 2014. Long time series analysis of ionospheric TEC disturbance over seismically active region in southwest China during low solar activity[J]. Seismology and Geology, 36(1): 253—265(in Chinese). [本文引用:3]
[5] 张雅雯, 陈美红, 李传起. 2013. 寒潮期间电离层TEC与地面天气参数的相关性[J]. 南京信息工程大学学报, 5(3): 278—283.
ZHANG Ya-wen, CHEN Mei-hong, LI Chuan-qi. 2013. Correlation analysis between ionospheric TEC and ground weather parameters during cold wave[J]. Journal of Nanjing University of Information Science & Technology, 5(3): 278—283(in Chinese). [本文引用:1]
[6] Afraimovich E L. 2004. Variations of the total electron content in the ionosphere from GPS data recorded during the Hector Mine earthquake of October 16, 1999, California[J]. Russian Journal of Earth Sciences, 6(5): 339—354. [本文引用:1]
[7] Afraimovich E L, Astafyeva E I, Oinats A V, et al. 2008. Global electron content: A new conception to track solar activity[J]. Annales Geophysicae, 26(2): 763—769. [本文引用:1]
[8] Afraimovich E L, Astafyeva E I, Zhivetiev I V. 2006. Solar activity and global electron content[J]. Doklady Earth Sciences, 409(2): 921—924. [本文引用:1]
[9] Aggarwal M. 2015. Anomalous changes in ionospheric TEC during an earthquake event of 13—14 April 2010 in the Chinese Sector[J]. Advances in Space Research, 56(7): 1400—1412. [本文引用:1]
[10] Balan N, Bailey G J, Su Y Z. 1996. Variations of the ionosphere and related solar fluxes during solar cycles 21 and 22[J]. Advances in Space Research, 18(3): 11—14. [本文引用:1]
[11] Belehaki A, Kutiev I, Marinov P, et al. 2017. Ionospheric electron density perturbations during the 7—10 March 2012 geomagnetic storm period[J]. Advances in Space Research, 59(4): 1041—1056. [本文引用:1]
[12] Bilitza D. 2001. International Reference Ionosphere 2000[J]. Radio Science, 36(2): 261—275. [本文引用:1]
[13] Dobrovolsky I P, Zubkov S I, Miachkin V I. 1979. Estimation of the size of earthquake preparation zones[J]. Pure and Applied Geophysics, 117(5): 1025—1044. [本文引用:1]
[14] Guo J, Wang L, Liu X, et al. 2015. Temporal-spatial variation of global GPS-derived total electron content, 1999—2013[J]. PLoS One, 10(7): e0133378. [本文引用:1]
[15] Ho Y Y, Jhuang H K, Su Y C, et al. 2013. Seismo-ionospheric anomalies in total electron content of the GIM and electron density of DEMETER before the 27 February 2010 M8. 8 Chile earthquake[J]. Advances in Space Research, 51(12): 2309—2315. [本文引用:1]
[16] Klimenko M V, Klimenko V V, Zakharenkova I E, et al. 2016. Longitudinal variation in the ionosphere-plasmasphere system at the minimum of solar and geomagnetic activity: Investigation of temporal and latitudinal dependences[J]. Radio Science, 51(12): 1864—1875. [本文引用:1]
[17] Le H, Liu J Y, Liu L. 2011. A statistical analysis of ionospheric anomalies before 736 M6. 0+earthquakes during 2002—2010[J]. Journal of Geophysical Research: Space Physics, 116(A2): A02303. doi: 10.1029/2010JA015781. [本文引用:1]
[18] Leonard R S, Jr Barnes R A. 1965. Observation of ionospheric disturbances following the Alaska earthquake[J]. Journal of Geophysical Research, 70(5): 1250—1253. [本文引用:1]
[19] Liu J Y, Chen Y I, Chen C H, et al. 2009. Seismoionospheric GPS total electron content anomalies observed before the 12 May 2008 MW7. 9 Wenchuan earthquake[J]. Journal of Geophysical Research: Space Physics, 114(A4): 231—261. [本文引用:3]
[20] Liu J Y, Chen Y I, Chuo Y J, et al. 2001. Variations of ionospheric total electron content during the Chi-Chi earthquake[J]. Geophysical Research Letters, 28(1): 1383—1386. [本文引用:1]
[21] Liu L B, Wan W X, Chen Y D, et al. 2011. Solar activity effects of the ionosphere: A brief review[J]. Chinese Science Bulletin, 56(12): 1202—1211. [本文引用:1]
[22] Orus R, Cand er L R, Hernand ez-Pajares M. 2007. Testing regional vertical total electron content maps over Europe during the 17—21 January 2005 sudden space weather event[J]. Radio Science, 42(3): RS3004. doi: 10.1029/2006Rs003515. [本文引用:1]
[23] Pulinets S A, Contreras A L, Bisiacchigiraldi G, et al. 2005. Total electron content variations in the ionosphere before the Colima, Mexico, earthquake of 21 January 2003[J]. Geofísica Internacional, 44(4): 369—377. [本文引用:1]
[24] Pulinets S A, Gaivoronska T B, Contreras A L, et al. 2004. Correlation analysis technique revealing ionospheric precursors of earthquakes[J]. Natural Hazards and Earth System Sciences, 4(5): 697—702. [本文引用:1]
[25] Reinisch B W, Dand enault P B, Galkin I A, et al. 2018. Investigation of the electron density variation during the 21 August 2017 solar eclipse[J]. Geophysical Research Letters, 45(3): 1253—1261. [本文引用:1]
[26] Shah M, Jin S. 2015. Statistical characteristics of seismo-ionospheric GPS TEC disturbances prior to global MW≥5. 0 earthquakes(1998—2014)[J]. Journal of Geodynamics, 92: 42—49. [本文引用:1]
[27] Tao D, Cao J, Battiston R, et al. 2017. Seismo-ionospheric anomalies in ionospheric TEC and plasma density before the 17 July 2006 M7. 7 south of Java earthquake[J]. Annales Geophysicae, 35(3): 589—598. [本文引用:1]
[28] Thomas J N, Huard J, Masci F. 2017. A statistical study of global ionospheric map total electron content changes prior to occurrences of M≥6. 0 earthquakes during 2000—2014[J]. Journal of Geophysical Research: Space Physics, 122(2): 2151—2161. [本文引用:1]
[29] Vapnik V N. 1999. An overview of statistical learning theory[J]. IEEE Transactions on Neural Networks, 10(5): 988—999. [本文引用:1]
[30] Zhao B, Wang M, Yu T, et al. 2008. Is an unusual large enhancement of ionospheric electron density linked with the 2008 great Wenchuan earthquake?[J]. Journal of Geophysical Research: Space Physics, 113(A11): A11304. doi: 10.1029/2008JA013613. [本文引用:2]
[31] Zhu F, Lin J, Su F, et al. 2016. A spatial analysis on the ionospheric TEC anomalies prior to M7. 0+earthquakes during 2003—2014[J]. Advances in Space Research, 58(9): 1732—1738. [本文引用:1]