结合PCA和滑动时窗的电离层异常探测方法研究

谢 婷 陈必焰 匡翠林

1 中南大学地球科学与信息物理学院, 长沙市麓山南路932号, 410083

多年的观测实践与大量的震例研究[1-3]表明，地震发生前电离层会出现异常变化。目前，传统电离层异常探测的方法主要包括中位数法、平均值法、滑动四分位距法和滑动时窗法等，而传统电离层异常探测方法是基于参考背景值与标准差来进行探测的，这就决定了高精度的参考背景值和合理的异常判别限差是获取精准异常探测结果的前提。基于此，研究者们陆续提出了许多异常探测新方法[4-5]，能更有效地提取信息、细化信号，并提高结果的可靠性，对地震异常探测具有积极的影响。随着研究的深入，有学者也尝试了一些无需构造背景场的异常探测方法，如主成分分析(PCA)方法[6-8]。

但是，利用PCA方法探测异常也存在一些不足之处：1)PCA方法无法提取出异常的空间特征，不利于震前异常研究；2)因为主特征值容易出现极端高值(如磁暴)，在进行长时间序列异常探测时，直接利用主特征值进行异常分析会削弱PCA方法的异常探测能力；3)目前的PCA方法直接采用经验值(0.5)为阈值进行异常判断，容易出现疏漏和失误，且检测得到的结果无法区分正负异常。针对这些问题，本文对PCA方法作出了改进，首先利用PCA方法初步过滤空间环境(地磁活动和太阳活动)等非震因素对电离层TEC扰动的影响，再采用滑动时窗方法进行电离层异常探测，同时与传统探测方法(滑动四分位距法和滑动时窗法)进行对比分析，验证结合PCA与滑动时窗方法的有效性，并利用本文方法初探了2020年瓦哈卡地震震前电离层时空变化特性，为后续地震研究工作提供可靠信息。

1 结合PCA与滑动时窗的探测方法

PCA是由Pearson[9]提出的，在数据挖掘和数据分析等领域得到了广泛应用。PCA方法主要是通过特征分解得到数据的主成分(特征向量)与特征值，进而获得数据的主要信息。

本文基于全球电离层网格(GIM)数据，研究PCA提取的TEC主特征值所代表的物理含义，以2014-01-12～02-12主特征值时空图为例进行探讨。从图1可以看出，在纬度±15°处，主特征值与地磁线走势吻合，且主特征值会随时间和位置发生变化，其中02-02于田地震震中上方区域观测到明显的高值异常，该异常关于磁赤道线对称，这与研究发现的震前电离层异常特征[10]是一致的。基于这些特征推测，PCA分离出的主特征值可能是基础背景场与异常的叠加。

图中，五角星为2014-02-12于田地震震中位置，圆圈为地震影响范围，红色曲线为纬度±15°处的地磁线

滑动时窗法是一种在给定的特定窗口大小的数组或字符串上执行要求操作的算法，其操作简单、运行快速，是一种常用的电离层异常探测方法。但该方法探测结果的准确度取决于背景值的精度，易受太阳活动等非震因素的影响，因此其探测精度有限。而PCA提取的主信号能初步过滤掉太阳活动等非震因素对电离层TEC的影响，在PCA方法的基础上采用滑动时窗法进行异常探测能更有效地检测出地震引起的电离层异常，且滑动时窗主要是基于滑动窗口来进行异常探测的，可以很好地解决PCA方法时序探测有限的问题。

结合PCA和滑动时窗(PCA-滑动时窗)方法的具体数据处理流程为：1)基于单点TEC时间序列，利用PCA求解出特征值；2)对特征值进行排序，求得最大特征值(单天1个主特征值)，多天的数据处理可以得到TEC的主特征值序列；3)采用滑动时窗法对主特征值序列进行检测，滑动窗口为16 d，以1.5倍均方差为探测依据，得到异常序列；4)由图1可知，主特征值之间相差较大会导致异常大小无法被很好地对比分析，因此对异常数据进行标准化处理(该处理不影响异常探测结果)。其中，电离层异常的判定标准为：异常出现在震中位置附近独立存在的区域，在位置或数值上与赤道电离层异常区[11](磁纬0°～±15°)有明显分离[2]。

2 分析与讨论

本文分别以2008年汶川地震、2010年玉树地震及2014年于田地震为例，研究PCA-滑动时窗方法的异常探测能力，其中汶川地震和玉树地震震前地磁活动频繁，且玉树地震震前地磁活动的剧烈程度较高，可用来探究PCA方法的抗地磁活动干扰能力；而于田地震震前的太阳活动可用来研究太阳活动对PCA方法探测异常的影响。最后，利用PCA-滑动时窗方法研究2020年瓦哈卡地震的震前异常现象，为后续地震研究工作提供可靠信息。4次地震的基本信息见表1。

表1 4次地震的基本信息

2.1 汶川、玉树及于田地震分析

本文选取了2008-01-01～05-20、2010-01-01～05-20及2014-01-01～05-20的GIM数据，采用PCA-滑动时窗方法进行电离层时空异常探测。由于不同震级对应孕震区的范围不同，因此采用经验公式R=100.43M(R为孕震区半径，M为震级)计算地震的影响范围[12]。此外，本文还选取了太阳活动和地磁活动指数数据，以表征地磁及太阳的活动情况，其中Dst为中低纬地区地磁活动指数，Dst >-30 nT表明地磁处于平静状态，-50 nT < Dst <-30 nT表明可能发生弱磁暴，Dst <-100 nT表明可能发生大磁暴；KP为一种判断地磁扰动总强度的指标，可将KP= 0～2、3、4、5、6、7～9的地磁活动分为平静、不平静、活动、小、大、强烈等不同等级；F10.7为表示太阳活动强弱的常用指标，其值通常在60～300之间，F10.7平稳微弱浮动表示太阳活动平静。

2.1.1 3个地震的时间序列异常分析

滑动四分位距法[13]是目前被认可并广泛使用的异常探测方法，同时滑动时窗法也是传统方法中相对合理的一种方法，因此本文从已有研究中选取这2种方法的震前异常探测结果，并结合PCA方法的探测结果与PCA-滑动时窗方法进行验证对比。由于目前PCA方法仅用于单点探测，因此本文基于GIM数据内插出震中TEC值的时间序列，再采用PCA方法进行异常探测，其中表2为4种方法探测得到的前3个地震的异常结果，滑动四分位距法和滑动时窗法的探测结果在不同研究中有所不同，但总体上是一致的。

表2 基于PCA-滑动时窗方法、PCA和前人研究资料汇总得到的3个地震的探测结果

表2中PCA-滑动时窗方法的探测结果没有排除空间环境(太阳活动和地磁活动)的影响，但得到的异常值与前人研究中排除空间环境影响的结果基本一致，因此可初步推测PCA-滑动时窗方法探测异常不易受地磁活动和太阳活动的影响。分析探测结果可发现，PCA方法在汶川地震震前21 d和19 d探测到了异常，但经过对孕震区内多点的异常进行检测发现该异常不具备区域性，即该异常只存在于相邻两点之间，这也表明仅利用震中上空单点进行异常探测存在局限性，需进一步开展时空异常探测。本文结果有别于汤俊等[6]利用PCA方法探测得到的结果，因为本文选取了2008-01-01～05-20共141 d的数据来进行异常探测。在进行长时间序列异常探测时，主特征值容易出现异常高值，直接利用其进行异常分析会削弱PCA方法的异常探测能力，这也是该方法的缺陷之一。从表2还可以看出，PCA-滑动时窗方法对这一缺陷有很好的改进，其结果与前人研究结果有很好的一致性。由于前人研究地震前兆一般只研究震前16 d左右异常，因此利用滑动四分位距法和滑动时窗法进行震前异常分析的资料汇总时间均不超过16 d。并且通常来说，如果某天有地磁或太阳活动，那么当天的电离层会出现明显异于平静期的扰动，且该扰动一般是全球尺度的，其中赤道区域电离层最易受太阳和地磁影响[17]，其扰动程度会更大。

2.1.2 汶川地震空间异常分析

图2(a)～2(c)给出了2008-04-21～05-14期间各时段的KP指数、Dst指数及F10.7指数，由图可知，汶川地震前后除了04-23、04-24、04-26、04-28、04-30及05-01和05-02时段的KP指数大于3，且04-23、04-24、04-25和04-26时段的Dst指数低于-30 nT外，其余时段的KP指数和Dst指数基本处于较低水平，而F10.7指数在整个时段内都非常稳定。

图2(d)为PCA-滑动时窗方法探测得到的二维异常分布，可以看出，该方法探测到的异常时段为04-29(震前13 d)、05-02(震前10 d)、05-05(震前7 d)、05-06(震前6 d)、05-07(震前5 d)和05-09(震前3 d)，这与前人研究结果具有极高的一致性，其中与滑动四分位距法和滑动时窗法相比，PCA-滑动时窗方法检测到的异常更全面。从图2(e)可以看出，PCA方法检测出的异常时段为04-21和04-23(震前21 d和19 d)。由图2(a)～2(c)可知，04-23～26的Dst指数低于-30 nT，查询空间环境预报中心的报告[18]发现，地磁活动在04-23处于中等扰动水平，对比图2(d)的异常空间分布可知，04-23异常主要出现在磁赤道线附近。从异常情况来看，04-29和05-05～07异常减少，05-02南部异常减弱、东部异常增大，05-09异常增加，其中电离层异常出现在震中偏东南方向，呈共轭结构，这与前人采用其他方法探测得到的结果是一致的[1,15]，证实PCA-滑动时窗方法可准确区分正负异常，弥补了PCA方法的缺陷。

(a)Dst指数；(b)KP指数；(c)F10.7指数(虚线是各指数设定的阈值，Dst为-30 nT，KP为3)；(d)基于PCA-滑动时窗法探测得到的二维异常分布(五角星为2008-05-12汶川地震震中位置，红色曲线为磁赤道线(磁纬0°)和磁纬±15°)；(e)利用PCA分析电离层TEC所得震中位置的主特征值序列(水平红线为0.5阈值线，垂线红线为地震发生日)

2.1.3 玉树地震空间异常分析

通过分析图3的太阳活动和地磁活动数据可以发现，04-05～08和04-11～12全球发生了中等地磁暴，但在图3(d)对应时间段的异常分布中并没有出现全球尺度或大范围的异常，验证了PCA方法探测异常不易受地磁活动干扰这一推测。根据空间环境预报中心发布的报告[18]可知，04-05全球发生了一次强烈的空间天气事件，并从04-06开始，地球同步轨道连续6 d发生高能电子通量增强事件，命名为“清明节事件”。因此，该时段内的TEC异常很可能是地磁活动引起的，本文方法在04-05～10检测到的异常位置主要位于赤道异常区和美洲扇区，与磁暴引起的电离层特征[17]吻合。相比于利用传统方法探测到的全球尺度大范围异常，PCA-滑动时窗方法会削弱地磁活动对电离层的影响。另外，整个时段内太阳活动水平较低，F10.7指数变化幅度小，整体维持在75～80 sfu附近，因此PCA-滑动时窗方法探测到的03-25(震前19 d)异常为正异常，03-30(震前14 d)异常为负异常，04-01(震前12 d)异常为正异常，这3次电离层异常都出现在震中偏东南方向，呈共轭结构。04-13地震当天震中以南部分区域出现TEC异常增加的现象，与前人研究获取的地震电离层扰动特征非常相似[1-2]，很可能是地震引起的电离层异常。

(a)Dst指数；(b)KP指数；(c)F10.7指数(虚线是各指数设定的阈值，Dst为-30 nT，KP为3)；(d)基于PCA-滑动时窗法探测得到的二维异常分布(五角星为2010-04-13玉树地震震中位置，红色曲线为磁赤道线(磁纬0°)和磁纬±15°)；(e)利用PCA分析电离层TEC所得震中位置的主特征值序列(水平红线为0.5阈值线，垂线红线为地震发生日)

2.1.4 于田地震空间异常分析

于田地震前后太阳活动和地磁活动情况如图4(a)～4(c)所示，从图4(c)中可以看出，震前太阳活动并不稳定，在01-22、01-28、02-05和02-06出现了异常上升或下降，太阳活动频繁，其中F10.7为200 sfu的峰值出现在震前6 d、7 d和8 d。

图4(a)和4(b)分别显示了Dst指数和KP指数的时间序列，综合分析2个参数发现，在02-07之前，Dst指数大于-30 nT，KP指数小于3，但在02-07～10期间，Dst指数达到-40 nT，KP指数超过5，查询空间环境预报中心[18]资料发现，该时段内地磁有3 h达到小磁暴水平，21 h达到活跃水平。

综上所述，太阳活动频繁期为01-22(震前21 d)、01-28(震前15 d)、02-04(震前8 d)、02-05(震前7 d)和02-06(震前6 d)，地磁活动频繁期为02-07～10(震前2～5 d)。从图4(d)的异常时空分布中可发现，在太阳活动频繁期和地磁活动频繁期内，异常主要集中在磁赤道附近和美洲扇区，这与电离层对地磁暴和太阳活动的主要响应特征一致，且本文方法在此期间并没有探测到全球大范围尺度的异常。结合前2个震例研究发现，PCA-滑动时窗方法探测异常不易受太阳活动和地磁活动的影响，因此本震例探测的02-02和02-12异常为正异常，且均能在磁共轭区被观测到，该探测结果与前人研究异常分析的结果比较一致[3]。

(a)Dst指数；(b)KP指数；(c)F10.7指数(虚线是各指数设定的阈值，Dst为-30 nT，KP为3)；(d)基于PCA-滑动时窗法探测得到的二维异常分布(五角星为2014-02-12于田地震震中位置，红色曲线为磁赤道线(磁纬0°)和磁纬±15°)；(e)利用PCA分析电离层TEC所得震中位置的主特征值序列(水平红线为0.5阈值线，垂线红线为地震发生日)

这3个震例验证了PCA-滑动时窗方法探测电离层异常的有效性，并具有以下优势：1)在数据处理过程中剔除掉了许多无关因素，计算速度快，探测到的异常情况比较直观、简洁；2)能进行长时间序列的时空异常探测，并能很好地区分正负异常；3)不易受地磁活动和太阳活动的影响。

2.2 2020年瓦哈卡地震研究

墨西哥位于三大构造板块之上，地壳板块的相对运动导致其频繁发生地震和爆发火山。2020-06-23在墨西哥瓦哈卡附近发生了7.4级地震，目前还没有关于该地震震前电离层异常的研究，因此本文利用PCA-滑动时窗方法对该地震的电离层TEC资料进行处理和分析，初探瓦哈卡地震震前电离层TEC的时空变化特性，为后续的地震研究工作提供可靠信息。

图5为瓦哈卡地震震前21 d至地震后2 d共24 d的F10.7、Dst及KP等指数的变化情况，由图可知，这24 d内空间天气比较平静。图5(d)为异常探测结果，可以看出，06-03(震前20 d)TEC异常减少，06-10(震前13 d)TEC异常增加，06-11(震前12 d)、06-12(震前11 d)、06-14(震前9 d)及06-17～19(震前4～6 d)电离层TEC异常减少，但由于06-03在23.274°S、68.468°W发生了6.8级地震，当天异常可能与该地震有关。从空间上来看，地震异常位于震中偏南区域，主要表现为低值异常。

(a)Dst指数；(b)KP指数；(c)F10.7指数(虚线是各指数设定的阈值，Dst为-30 nT，KP为3)；(d)基于PCA-滑动时窗法探测得到的二维异常分布(五角星为2020-06-23瓦哈卡地震震中位置，红色曲线为磁赤道线(磁纬0°)和磁纬±15°)；(e)利用PCA分析电离层TEC所得震中位置的主特征值序列(水平红线为0.5阈值线，垂线红线为地震发生日)

3 结 语

本文基于GIM数据，利用PCA-滑动时窗方法探测了4次地震的震前电离层异常，首先利用PCA-滑动时窗方法探测了3次强震事件(汶川地震、玉树地震及于田地震)，并将结果与滑动时窗法和滑动四分位距法的结果进行了对比验证，认为本文方法具有有效性及探测优势。结果表明：1)PCA-滑动时窗方法能够应用于地震电离层的时空异常探测，且探测结果直观、简洁；2)相较于PCA方法，PCA-滑动时窗方法能够准确探测出正负异常，并能更好地应用于长时间序列的时空异常探测；3)相较于传统方法，本文方法探测电离层异常时不易受太阳活动和地磁活动的影响，探测出的异常结果可信度高。本文在验证了PCA-滑动时窗方法可行性的基础上对2020年瓦哈卡地震进行了异常探测，结果表明，06-10(震前13 d)TEC异常增加，06-11(震前12 d)、06-12(震前11 d)、06-14(震前9 d)及06-17～19(震前4～6 d)电离层TEC异常减少，这些异常很可能是地震的前兆异常。从空间位置来看，异常位于震中偏南区域，TEC主要表现为低值异常。综合发现，电离层TEC异常主要出现在震中偏南方向，且呈共轭结构。由于发震机制复杂，目前电离层与地震之间的耦合关系尚不清楚，因此将TEC异常作为一种前兆异常来进行地震预测非常重要。本文方法可很好地进行多震例异常特征(出现时间、空间方位、正负情况等)的提取及统计研究，从而更好地应用于电离层异常监测。