引言

中强地震具有较强的破坏性，给人们生命财产带来巨大的损失。然而，地震产生的机理相当复杂，世界各国对地震的预测预报研究一直处于探索阶段。如何提升地震监测预报的有效性一直是各国研究人员面对的难题，众多的地震工作者不断努力的利用各种分析方法在不同领域进行探索。

自1964年美国阿拉斯加大地震时，Leonard等^[1]发现地震期间电离层有异常扰动现象后，各国学者便开始从事电离层的变化与地震之间关系的研究工作，并逐渐成为当前研究的热点之一。Antsilevich^[2]发现1966年塔什干地震期间电子含量有明显增加的现象；Weaver等^[3]发现1969年千岛群岛地震期间电离层也出现了异常扰动。大量研究发现，大于5级的地震发生前几天到几个小时，震中附近上空电离层都会出现异常现象^[4-6]。

近年来，我国有众多相关领域的研究人员在关于电离层TEC变化与地震的相关性方面做了大量的研究工作。马一方等^[7]研究发现芦山地震前在震中附近上空东向和北向梯度VTEC均出现了负异常。姚宜斌等^[8]利用震中附近GNSS观测数据，采用滑动窗口法发现2011年3月11日本地震前电离层有扰动现象。张小红等^[9]利用IGS提供的震中附近4个格网点TEC数据，采用时间序列法得到2012年1月10日苏门答腊岛7.2级地震前13天电离层TEC时间序列，并分析了时间序列法、传统滑动时窗法和四分位法在预测电离层TEC参考背景值和精度方面的差异。姚璐等^[10]和张明敏等^[11]先后研究了玉树地震和九寨沟地震前后电离层TEC的变化异常。

本文采用迭代法和滑动四分位距法对同一震例进行电离层TEC数据分析。说明不同方法在具体震例中的应用及各自的特点，为通过分析电离层TEC异常预测地震发生积累经验。

1.   数据震例及分析方法

1.1   数据及震例

本文采用经NASA处理后的IGS全球电离层格网数据进行电离层异常分析，其空间分辨率为5°（经度）×2.5°（纬度），时间分辨率为2小时。

选取震例为北京时间2019年10月12日22时55分发生在广西壮族自治区玉林市北流市（22.18°N，110.51°E）的5.2级地震（以下简称北流地震），震源深度10 km，图1为震中及格网点示意图。

图  1  北流地震震中及格网点示意图

Figure  1.  Schematic diagram of epicentral and grid points of Beiliu earthquake

下载: 全尺寸图片 幻灯片

选取该震例主要是为了便于在迭代法研究预测发震区域的基础上，采用不同方法对于其预测区域发生地震震中周边电离层TEC数据进行数据分析，在较短时间序列内发现电离层TEC变化异常，从而可以组成中长期预测指标和短临预测指标相结合的前兆异常分析方法。

1.2   分析方法

1.2.1   迭代法

迭代法是一种不断用变量的旧值递推新值的过程。使用迭代法进行电离层TEC数据分析过程中，空间分辨率为5°（经度）×2.5°（纬度），时间分辨率为2 h，连续对20年来全球经度（−180°—180°），纬度（−87.5°—87.5°）范围内的5183（73×71）个格点TEC数据进行系统处理分析。作者研究中为了更好的在较大区域范围内找到显著的电离层TEC异常，经过多次尝试后最终选择7×7组成的大格网（经度跨度35°，纬度跨度17.5°）区域进行分析应用。在数据分析中，将各格点TEC值数据转化成数据变化率，在大格网区域内，以本格网为中心，将本格网变化率数据与周边48个网格单元对比，若本格网数据全部大于或小于周边网格，则判断当前格网出现了一次异常，然后以365日为滑动窗口计算格网的异常频次年均线，最后将年均值与历史最低年均值比较，当年均线值比历史最低年均线值增大1倍时视为均线异常。

由文献[12]可知，在TEC数据分析中，以年均值线的突变现象作为判断标准，突升是TEC活跃性增强异常，反之则是TEC活跃性减弱异常，突变指数越大异常越明显。经对比研究发现：日本3•11大地震、汶川地震、玉树地震、雅安地震、九寨沟地震等地震震前震中区或周邻地区均出现了明显的TEC活跃性变化特征。以2008年5月12日发生的汶川8.0地震（31°N，103°E）为例，地震前后震中所在区域TEC活跃性发生明显变化，从2005年开始缓慢上升，2007年上升加剧，短时间内就达到分析方法设定的异常判断标准，即增大1倍，2007年底达到峰值，随后呈下降趋势。之所以出现这种情况，根据地震发生的机理判断可能是因为2007年地下板块运动达到顶峰，可运动空间变小，导致活动性减弱所致。TEC活跃性至地震发生时仍处于异常判断标准之上，地震发生后其活跃性快速下降至异常判断标准之下（图2a）。分析以汶川为中心的大区域TEC活跃性，发现汶川周围地区震前出现高值异常，且这种现象在地震前持续时间超过1年（图2b）。

图  2  汶川地震TEC活跃性指数变化图（a）与TEC活跃指数时段平面图（2008年5月） （b）

Figure  2.  Changes of TEC activity index in Wenchuan earthquake （a） and plan diagram of TEC activity index period （May 2008） （b）

下载: 全尺寸图片 幻灯片

1.2.2   滑动四分位距法

作者在进行迭代法研究电离层变化异常工作的同时，提出采用滑动四分位距法对预测区域发生地震前后电离层数据进行分析的思路，力求寻找电离层TEC变化短临异常与中强地震的相关性。

四分位距是一种稳健统计技术中用于表示数据离散度的一个量，常用来检查数据的异常情况。四分位距法在数据分析中具有能较为准确地获得预测参考背景值，并可以比较准确的计算出背景参考值上下限的特点。所谓的四分位数就是将数列分为4部分，一个数列可以有3个四分位数，即下四分位数、中位数和上四分位数。以一个16个数的数列为例，将其从小到大排列为x₁，x₂，……，x₁₆，则

上四分位数：Q₁=$\dfrac{{{x_4} + {x_5}}}{2}$

中位数：Q₂=$\dfrac{{{x_8} + {x_9}}}{2}$

下四分位数：Q₃=$\dfrac{{{x_{12}} + {x_{13}}}}{2}$

四分位距：IQR=Q₃−Q₁

其中，Q₁表示在该数值以下的数据占总数的25%；Q₂表示在该数值以下的数据占总数的50%；Q₃表示在该数值以下的数据占总数的75%。

在统计学上IQR=1.34σ，即四分位距的期望值是标准差的1.34倍。在本文数据分析中将Q₂±1.5IQR作为判定TEC是否为异常值的标准控制值，即上限=Q₂+1.5IQR，下限=Q₂−1.5IQR。其阈值约为标准差的2倍，异常检验置信度约为95%。

滑动四分位距法可以看做滑动窗口法与四分位距法的结合，就是对一段时间序列内的数据采用四分位距法进行数据处理，从而得到较为合理的控制限值用来判读变化异常。由于电离层变化受季节影响较大，选择时间窗口不宜过长。因此，在具体数据分析中以15天作为时间窗口，即分析每天数据时取当天及前14天数据组成数列，进行数据分析，实际观测值超出上下限值即为异常值。

2.   TEC数据分析处理

北京时间2019年10月12日22时55分在广西玉林市北流市（22.18°N，110.51°E）发生5.2级地震，震源深度10 km。

迭代法分析结果显示，自2018年中广西北部湾出现TEC异常且持续，直至2019年7月异常达到近10年来最大值，随后异常出现降低，但仍处于分析方法设定值以上，直至地震发生，地震发生后TEC异常呈快速降低形态，2020年初基本回归原活动范围（图3a）。分析以震中区域为中心的大区域TEC活跃性，发现震中周围地区震前出现高值异常，且这种现象在地震前持续时间超过1年（图3b）。

图  3  北流地震TEC活跃性指数变化图（a）与TEC活跃指数时段平面图（2019年7月） （b）

Figure  3.  Variation chart of TEC activity index of Beiliu earthquake （a） and time interval plan of TEC activity index （July 2019） （b）

下载: 全尺寸图片 幻灯片

选取2019年9月1日—10月16日期间，地震发生时间前15天，发生后3天组成的时间序列，以15日为时间窗口应用滑动四分位法进行数据分析处理。下面图中依次显示了震中周围各格网点的TEC异常变化情况。通过数据分析可以得出，各格网点TEC异常变化具有较好的一致性，震前10—14天、7天、3天、2天、当日均出现正值异常，震前1天出现负值异常（图4）。

图  4  北流地震震中周围格网点TEC时间序列（黑色柱线为发震时间）

Figure  4.  TEC time series of grid points around the epicenter of the Beiliu earthquake （the black column is the time of seismogenesis）

下载: 全尺寸图片 幻灯片

根据相关研究文献[13]可知，电离层异常变化易受到如太阳活动变化、地磁活动异常、天气变化因素等多种因素的影响。因此，研究地震前后电离层异常变化需要考虑时间窗口内影响因素的变化情况。在下面的异常分析中将结合F_10.7指数、A_p指数、K_p指数、D_st指数等主要影响因素进行分析。相关太阳活动和地磁活动数据均来源于国家科学院空间环境预报中心（图5）。

图  5  9月27日—10月15日北流地震前后太阳活动和地磁变化情况

Figure  5.  Changes of solar activity and geomagnetism before and after the Beiliu earthquake from September 27 to October 15

下载: 全尺寸图片 幻灯片

根据国家科学院空间环境预报中心对2019年9—10月份的太阳活动和地磁观测数据报告可知，9月份太阳活动水平较低，无C级以上级别耀斑产生；地磁活动9月27—30日出现小磁暴现象，对TEC变化有所影响。10月份太阳活动水平依然较低，无C级以上级别耀斑产生；地磁活动10月9日K_p指数、A_p指数处于高值，10月10—11日，太阳风速上升至510 m/s左右，地磁活动水平较高，会对TEC变化有所影响。

因此，可以判断震前11—14天，震前1—3天TEC变化异常受到地磁影响因素较大，不宜作为地震前兆异常。

通过剔除对电离层TEC异常影响因素后，获得电离层TEC异常后还有震前10天、7天及地震当日3组异常，即10月2日、10月5日、10月12日。

3.   TEC空间分布异常分析

为了进一步分析震前3组TEC异常是否与此次地震有关，本文分别根据数据分析结果给出了3组异常高值前后（5°—40°N，75°—125°E）范围内电离层TEC异常分布图。

根据数据分析结果给出的10月5日异常高值出现在14:00UT，据此给出10月5日10:00UT—20:00UT之间的TEC空间分布图，发现震中附近区域仅在18:00UT出现较短时间异常（图6）。因此，可以认为该数据分析异常不宜作为地震前兆异常。

图  6  10月5日数据异常高值时间TEC空间分布图（黄色标记为震中位置）

Figure  6.  Spatial distribution of TEC with abnormally high data on October 5 （epicenter location is marked in yellow）

下载: 全尺寸图片 幻灯片

根据数据分析10月2日异常高值时间给出10月2日2:00UT—12:00UT间隔两小时震中周边区域TEC空间分布图。可以明显看出，TEC异常区域自2:00UT出现在震中东南方向，随时间推移逐渐沿22°N由东向西移动，8:00UT异常区域出现在震中附近，直至12:00UT远离震中区域（图7a）。

图  7  10月2日、12日TEC空间分布图（黄色标记为震中位置）

Figure  7.  The spatial distribution of TEC on October 2 and 12 （the epicenter is marked in yellow）

下载: 全尺寸图片 幻灯片

根据数据分析10月12日异常高值给出10月12日2:00UT—12:00UT间隔两小时震中周边区域TEC空间分布图。可以明显看出，TEC异常区域变化同10月2日情况基本一致，除异常区域自东向西移动外，6:00UT异常区域出现在震中附近，直至8:00UT异常区域远离震中区域（图7b）。

综合考虑可以认为，10月2日与10月12日的电离层TEC异常与广西北流地震有关。

因此，对中强地震发生前后数天内电离层TEC数据采用滑动四分位距法进行分析，尽管易受到太阳活动和地磁变化的影响，但将去除干扰因素后的异常作为地震前兆短临判断指标，仍具有一定的有效性。这一结论与其他研究人员分析地震前电离层TEC变化异常的研究结果基本一致^[14-17]。

4.   结语

本文采用迭代法和滑动四分位距法对同一震例中电离层TEC异常进行了分析。迭代法具有分析较长时间内大区域电离层TEC异常的特点，滑动四分位距法可以对单个格网较短时间序列内电离层TEC异常进行分析。在地震监测预报日常工作中可以尝试将迭代法和滑动四分位距法分析电离层TEC异常的方法综合使用来进行地震前兆判断研究，以达到中长期预报和短临预报相结合的目的。但由于地震发生的机理复杂多变，利用电离层扰动对地震进行预测存在明显的局限性，还有待于参考多种不同分析方法进行研究和震例印证，来提高地震预测的有效性。

致谢

感谢陕西省地震局方伟高工给予的指导和支持，以及西安市地震监测中心同事在研究过程中提供的帮助。