2017年棉兰老岛Mw7.0地震前TEC异常分析

孙 哲，刘智敏，杜自豪，孔晓宇

(山东科技大学 测绘科学与工程学院，山东 青岛 266590)

电离层是一种由大量电子组成的壳层，容易受到各种自然现象的干扰，如火山爆发、太阳耀斑、地震等。研究表明地震孕育过程中，在岩石圈、大气圈以及电离层确实存在电磁异常现象，且电离层前兆的出现有稳定的时间尺度，这一突出特点将使它在短临地震预报中更具实用价值[1]。随着全球定位系统(Global Positioning System，GPS)的发展，基于GPS观测数据反演获取电子总含量(Total Electron Content，TEC)的成本逐渐降低，讨论震前TEC异常变化的规律性成为热点[2-6]。

探测数据异常性的方法有很多，传统静态探测法有中位数法、平均值法等，这些方法不能剔除数据中的异常数据，使得计算的背景参考值存在较大偏差，导致探测结果不准确，为了弥补传统静态探测法背景参考值偏差较大的问题，学者提出动态探测法—滑动时窗法[7]、滑动四分位法[8]。这两种滑动时间窗口方法作为现在主流的动态探测法，虽然较传统静态探测法提高了背景参考值的精度，但对主要异常时段不够突出且探测时窗较长。为提高对异常时段的敏感度和缩短时间窗口长度，本文提出利用平滑Z分法进行震前分析并对2017年4月29日4时23分发生的菲律宾棉兰老岛7.2级地震(震中位于5.51°N，125.08°E，震源深度为57 km)进行震前异常探测，探测结果与滑动时窗法和滑动四分位距法的探测结果基本一致，结合太阳地磁场活动指数、TEC异常空间变化后推测3种方法探测到的异常与此次地震有关，进一步证实了平滑Z分法的可行性、有效性。

1 数据获取及异常检测方法

1.1 数据获取

本文在分析电离层异常时，采用的电离层数据来源于IGS提供的全球电离层格网数据，时间分辨率为 2 h，空间分辨率为 5°×2.5°(经度×纬度)，该数据是通过全球分布的IGS基准站观测资料计算得到的。根据文献[9]提出的孕震区半径公式R=100.43M(M代表震级，R代表孕震区半径(km))，可以得出菲律宾棉兰老岛地震孕震区半径约为1 247 km，全球电离层格网数据的空间分辨率可以满足此次地震电离层异常扰动分析的需求。地震前电离层异常区域往往偏离震中区域一定距离[10]，利用格网数据选取震中附近4个格网点A(125°E，5°N)，B(130°E，5°N)，C(125°E，7.5°N)，D(130°E，7.5°N)震前15 d和震后5 d的TEC时间序列进行分析效果会更好。震中以及4个格网点分布情况如图1所示，其中红色边框五角星表示震中点(红色实心五角星表示震中点在亚洲板块的位置)，红色圆表示格网点。

图1 震中及异常观测点示意图

由于太阳活动和地磁场的变化都会导致电离层的异常变化，在分析电离层异常时需结合该时期的空间环境状况综合分析该段时期的电离层异常变化原因。本文提取了震前15 d及震后5 d的太阳射电量F10.7指数、地磁Dst指数、地磁Kp指数的时间序列，各指数的评判标准已具备[11-12]。

1.2 异常探测方法

1.2.1 滑动时窗法

滑动时窗法是以滑动均值为背景参考值，滑动标准差为异常检测依据的探测方法，对于待探测的时间序列取l天作为时窗窗长，本文选择15 d作为时间窗口来计算均值和均方差[13]。

均值：

(1)

均方差：

(2)

在进行异常探测处理时，通常以倍均方差为探测依据，这里取2倍均方差，置信度超过95%，探测上下限：

(3)

待探测天数同一时刻若超出探测上下限值，即可视该天该时刻的TEC值为异常。若探测第i+1天某一时刻的异常，需计算前一天同一时刻的平均数及均方差，作为背景参考值，按照探测上下限对第i+1天进行判断。

1.2.2 滑动四分位法

结合以往的经验，综合考虑均值法和中位数法的不足，台湾国立中央大学的刘正彦教授提出四分位距法[2]。四分位距( Inter Quartile Rang，IQR) 是一种稳健统计技术中用于表示数据离散度的一个量，常用来检查数据的异常情况。在选择时间窗口时长度不宜过长，因此本文中选择15 d作为时间窗口，取其前15 d同一时刻的TEC组成一组时间序列，并按照从小到大顺序对其进行排列，为x1，x2，…，x15，则：

(4)

(5)

式中：Q1为上四分位，Q2为中位数，Q3为下四分位，k为常倍数，IQR为四分位距，LB为下边界，UB为上边界。本文定义Q2±1.5IQR作为TEC异常检验上下边界的判断阀值，约为2倍均方差，置信度在95%左右。如果TEC值没有超出上下界则ΔTEC=0，超出上界ΔTEC>0，表示正异常，低于下界值ΔTEC<0，表示负异常。

1.2.3 平滑Z分法

针对实时数据中的峰值信号检测，外国学者Jean-Paul构建了一个基于分散原则的稳健算法，称为平滑Z分法(Smoothed Z-score Method)。该算法使用探测时窗扫描一组数据点并计算移动平均值和标准差，高于设定阈值的Z分会被检测出来，Z分的方程：

(6)

平滑Z分法需要设定3个参数：探测时窗(l)，阈值(th)和影响因子(in)。影响因子设置在0～1之间，表示以前观测点对新观测点的影响比。平滑Z分法的组成部分：

(7)

(8)

(9)

式中：xi为检测到的异常值；si为平滑后的值。考虑本文数据特点，参考Mehmet等学者的研究，本文选取影响因子in=0.5，阈值th=2，时间窗口l=4 d，提取震前15 d和震后5 d每个时间节点(全球格网数据2 h为一个时间节点)的TEC值与对应前15 d同一时间节点TEC序列的中值做差得到趋势序列M，将趋势序列M作为观测值进行异常检测。异常判定从观测值第一个数据点开始，将前4 d的趋势序列作为背景参考值，检测到异常值进行平滑后继续向前判定，直至观测值判定完毕。

2 分析与讨论

2.1 时间序列异常分析

本文提取震前15 d到震后5 d的TEC时间序列作为分析对象，采用上述3种方法对A，B，C，D4个点进行TEC异常探测，结果如图2—图5所示。滑动四分位法结果显示4个格网点大部分时间TEC时间序列都在上下边界之内，从震前10 d开始TEC异常出现并逐渐增大，在震前9 d TEC异常达到最大，震前8 d逐渐缩小，震前7 d达到第二个峰值，震前6 d逐渐缩小，震前5 d到震后4 d之间缩小速度增快，震后5 d之后异常消失，4个格网点表现出一致性，整体呈现双峰结构，震前9 d在A点出现最大正异常约为9TECu，震后3 d在C点出现最大负异常约为-4TECu，A，C点异常值比B，D点高；滑动时窗法结果与滑动四分位法结果大致相同，4个格网点在震前10～6 d期间呈现出双峰结构，震前9 d在A点达到最大正异常约5TECu，震后3 d在A点达到最大负异常约-1TECu，相比于滑动四分位法，检测到的TEC异常节点相对较少，最大正异常和最大负异常比滑动四分位法都低，A，C点异常值高出B，D点程度和对异常时段的敏感度比滑动四分位法都低；平滑Z分法结果与前两种方法大致相同，震前9 d和7 d出现两个峰值，其中震前9 d在A点达到最大正异常约11TECu，震后3 d在C点达到最大负异常约-3TECu，相对于滑动四分位法和滑动时窗法异常程度更高，自身A，C点异常值较B，D异常值高，与前两种方法一致，探测到的异常节点较前两种方法相对较少，但震前9 d和震前7 d的异常节点占比提高，更加突出异常时段。观察图1可知A，C点比B，D距震中更近，参考3种方法的异常结果推测格网点异常值与距震中距离相关；图2—图5显示震前9 d和7 d两个异常峰值最为特殊，是否与此次地震相关还需进一步结合地磁活动指数和空间异常进一步分析。

图2 2017-04-29菲律宾棉兰老岛地震A点(125°E，5°N)TEC时间序列及滑动四分位法、滑动时窗法和平滑Z分法法结果(图中0点表示地震当天)

图4 2017-04-29菲律宾棉兰老岛地震C点(125°E，7°N)TEC时间序列及滑动四分位法、滑动时窗法和平滑Z分法法结果(图中0点表示地震当天)

图5 2017-04-29菲律宾棉兰老岛地震D点(130°E，7°N)TEC时间序列及滑动四分位法、滑动时窗法和平滑Z分法结果(图中0点表示地震当天)

结合太阳活动方面因素选取正异常最明显的震前9 d和7 d来进一步分析震前TEC异常原因，本文提取了震前15 d及震后5 d太阳地磁活动时间序列数据如图6所示。从图中可以看出，震前9 d和7 d中只有F10.7指数在设定阈值以内，其余各项指数均有所超标：震前7 d大部分时间的Dst指数超标，最大值达到6，Kp指数也超标，最低值接近-50 nT，说明这几天的太阳地磁活动较强烈，推测图2—图5中用滑动四分位法、滑动时窗法和平滑Z分法探测出震前7 d的TEC异常与太阳地磁活动干扰密切相关。

图6 2017-04-29菲律宾棉兰老岛地震地磁活动时间序列图(图中0点表示地震当天)

震前9 d出现很短时间的小幅度太阳活动异常，在UTC 08:00—09：00时Dst≈-40 nT，Kp=4，其余时间Dst≥-30 nT，Kp<4，F10.7全部低于100SFU，可见这天大部分时间太阳活动很平静，但上述3种方法均检测到这天UTC 00：00—14：00出现大量正异常，是否由孕育区地震引起还需从空间角度进一步研究。

2.2 空间分布异常分析

区分地震扰动与磁暴等活动的一个很重要的评判标准为电离层TEC异常是否具有局地性特征[14]。为了进一步分析震前9 d和震前7 dTEC异常是否与此次地震有关，需要进行空间上大范围的异常检测，为尽可能多的探测异常时段，突出空间表现效果，本文分别给出滑动四分位法在震前9 d和震前7 d，90°E～160°E，10°S～30°N 范围内的电离层异常图。

震前9 d UTC 08：00—18：00间隔2 h的TEC异常情况如图7所示，其中红色五角星为震中位置，由图7可知在UTC 08：00时，震中东南部出现了明显的正异常，异常幅度在6～8TECu；随后在UTC 10：00正异常向震中移动并且异常值逐渐增大，达到峰值，最大异常值约10TECu；UTC 12：00—14：00时正异常开始削减并向西北方移动；UTC 16：00—18：00正异常逐渐削弱直至消失。这一天TEC异常现象在UTC 08：00持续到UTC 16：00，时间较长，异常大部分出现在震中区域附近，而在其他地方没有出现明显的异常，推测震前9 d的电离层异常由此次地震引起。

图7 震前9 d(2017-04-20)TEC异常分布图(红色五角星为震中位置)

震前7 d UTC 02：00—12：00间隔2 h的TEC异常情况如图8所示，其中红色五角星为震中位置，由图8可知在UTC 04：00之前，整个区域中绝大部分的电离层很平静，只有在震中东南方向出现了微弱的电离层异常，异常值为5TECu左右。UTC 04：00在震中东西部开始出现电离层异常，此后，异常区域瞬间向外扩散且异常值也显著增加，UTC 08：00电离层异常达到峰值为12TECu左右，此时异常区域覆盖了整个分析区域的绝大部分。UTC 10：00之后，电离层异常迅速减弱消失，UTC 12：00整个分析区域只有很小一部分微弱异常。此次电离层异常较震前9 d强，覆盖范围广，异常时间缩短，结合图6，可知震前7 d地磁活动剧烈，推测震前7 d电离层异常与地磁活动有关。

图8 震前7 d(2017-04-22)TEC异常分布图(红色五角星为震中位置)

3 结束语

利用IGS提供的全球电离层格网数据对2017年菲律宾棉兰老岛地震分别用滑动四分位法、滑动时窗法和平滑Z分法进行震前TEC异常探测，主要结果如下：

1)3种方法均能探测到震前的TEC值异常信息。3种方法在震前5～10 d 4个点均呈现出明显的正异常，震前9 d和7 d均出现峰值，震前9 d在A点均出现最大正异常，探测结果显示3种方法在4个点A，B，C，D中表现一致，格网点异常强度与距震中距离有关，震前9 d的TEC值异常与地震相关；平滑Z分法较传统的动态测量方法探测到的异常天数较少，但与地震相关异常天数占比更高，对异常数据进行平滑处理，提高对异常时段的敏感度，探测时窗上更短为4 d，有利于短临地震的实时预报。

2)在空间异常分布的图像上来看，震前9 d正异常在震中东南方出现逐渐向震中靠拢并且异常值逐渐增大，随后逐渐向震中西北方散去，异常值也逐渐削减直至消失，期间异常值绝大部分时间都集中在震中附近，具有区域性；震前7 d异常在震中附近突然出现，增加速度极快，随后迅速消失，移动方向无规律，分布范围覆盖整个研究区，具有全域性。

通过对比3种方法的异常探测结果并结合空间异常分布和太阳地磁活动排除干扰因素，推测震前9 d的正异常与此次地震有关，震前7 d的正异常与太阳活动有关。此次探测结果证实了平滑Z分法的可行性和较传统动态测量方法的优越性，将几种方法互为补充可突出主要异常时段，对于短临地震预报有一定价值。