地震给人类带来巨大灾害，是危害最严重的突发性自然灾害之一。地震一旦发生，将可能伴随着人员伤亡、不可逆转的建筑坍塌或严重毁坏等。学者们致力于各种地震预报方法的研究，而震前电离层扰动与地震电离层耦合机制是目前地震预兆研究的热点^[1-13]。Shah等^[4]对1998-2014年间震级Mw≥5.0的地震进行分析，表明震前5 d内电离层会出现异常；Tang等^[5]从时空角度分析了雅安地震震前电离层总电子含量(total electronic content, TEC)扰动，时间上说明震前12~15 d的TEC异常可能是地震引起的，而空间上说明异常点通常出现于震中和震中共轭点；姚宜斌等^[6]分析了日本本州东海岸附近海域发生的Mw9.0级地震前后的电离层异常，确认3月8日的电离层异常可看作是日本“311”地震的前兆；Liu等^[7]对1999-2002年中国台湾地区20次Mw6.0级以上地震震前TEC进行统计分析，结果表明震前5 d电离层TEC扰动主要出现在10:00-14:00 UT，且有80%的大地震震前会出现TEC扰动；Ulukavak和Yalcinkaya^[8]对2010年加利福尼亚州Mw7.2级地震震前TEC变化进行分析，说明震前电离层TEC正负异常均会出现，且与地震有关的电离层扰动出现在震前1~5 d。

震前电离层TEC异常探测的传统方法主要有平均值法、中位数法、四分位距法(inter quartile range, IQR)、滑动时窗法等，其中平均值法和中位数法属于静态探测法，不能剔除样本中的异常数据，使背景值与实际值偏差很大，探测结果不准确^[9]，IQR和滑动时窗法是传统方法中相对合理的两种方法，但背景值预测精度仍不高，有一定的系统偏差；鉴于此，张小红等^[1]利用时间序列差分自回归移动平均(autoregressive integrated moving average, ARIMA)模型进行探测，提高了背景值预测精度，并用该方法对苏门答腊岛Mw7.2级地震电离层异常进行探测，但ARIMA模型极值点附近的预测精度并不高；李磊^[9]运用自回归移动平均(autoregressive moving average, ARMA)残差修正模型对扰动探测背景值进行递推预测，提高了极点附近背景值预测精度，并用其对汶川地震震前TEC扰动进行探测；Akhoondzadeh^[10]利用LSSVM(least squares support vector machine)方法预测尼泊尔Mw7.5级地震发生前周边TEC的参考背景值，并与两种经典方法进行对比，说明震前2 d的电离层TEC异常扰动主要由地震引起，并验证了该方法对震前电离层TEC扰动探测是可行的。

文献[1]、文献[9-10]虽在一定程度上提高了背景值预测精度，但均未考虑电离层TEC受各种因素影响，且在时空上的分布存在差异，运用单一模型对TEC进行预测，精度不高，影响探测结果，因此，本文提出一种将TEC数据分成多个频率尺度，并综合了文献[1]、文献[9]中模型优点的多尺度ARMA残差修正模型(multi-scale ARMA residual correction model，MARCM)。对比分析了MARCM、ARMA模型、IQR和滑动时窗法对电离层TEC参考背景值的预测精度，并探测了2013-04-20四川省雅安市芦山县Mw7.0级地震震前15 d以及地震当天的电离层TEC异常扰动情况。

1 探测原理及方法 1.1 四分位距法和滑动时窗法

四分位距法就是选取探测时段为D天(能被4整除)的电离层TEC数据，分别将D天中对应时刻的TEC数据放一起，然后从小到大排序，得到X₁≤X₂≤…≤X_D，并将其分成4份，分界点分别为Q₁、Q₂、Q₃，其表达式如式(1)所示，而上下界限u、l表示如式(2)所示^{[1, 9]}。

$ \left\{\begin{array}{l}{Q_{1}=\frac{1}{2}\left(X_{D / 4}+X_{D / 4+1}\right)} \\ {Q_{2}=\frac{1}{2}\left(X_{D / 2}+X_{D / 2+1}\right)} \\ {Q_{3}=\frac{1}{2}\left(X_{3 D / 4}+X_{3 D / 4+1}\right)}\end{array}\right. $

(1)

$ \left\{\begin{array}{l}{u=Q_{2}+1.5\left(Q_{3}-Q_{2}\right)} \\ {l=Q_{2}-1.5\left(Q_{2}-Q_{1}\right)}\end{array}\right. $

(2)

滑动时窗法通过选取适当的滑动窗口长度N，计算选定窗口的均值X_K和标准差σ_K，以X_K为参考值，以两倍σ_K(相应检测置信度水平为95.5%)为参考值的上下限值，X_K和σ_K的表达式如式(3)所示，异常值Δ_TEC检测表达式如式(4)所示。

$ \left\{\begin{array}{l}{\overline{X}_{K}=\frac{1}{N} \sum\limits_{L=K+1}^{K+N} X_{L}, K=0, 1, 2 \cdots} \\ {\sigma_{K}=\sqrt{\frac{1}{N} \sum\limits_{n=1}^{N}\left(X_{n}-\overline{X}_{K}\right)^{2}}, K=0, 1, 2 \cdots}\end{array}\right. $

(3)

$ \varDelta_{\mathrm{TEC}}=\left\{\begin{array}{l}{0, \overline{X}_{K}-2 \sigma_{K} \leqslant \varDelta_{\mathrm{TEC}} \leqslant \overline{X}_{K}+2 \sigma_{K}} \\ {Y-\left(\overline{X}_{K}+2 \sigma_{K}\right), \varDelta_{\mathrm{TEC}}>\overline{X}_{K}+2 \sigma_{K}} \\ {Y-\left(\overline{X}_{K}-2 \sigma_{K}\right), \varDelta_{\mathrm{TFC}}<\overline{X}_{K}-2 \sigma_{K}}\end{array}\right. $

(4)

式中，Y为IGS提供的TEC值。

1.2 多尺度ARMA残差修正模型 1.2.1 模型原理

多尺度ARMA残差修正模型主要由改进集总经验模态分解(modified ensemble empirical model decomposition，MEEMD)^[14]及时间序列ARMA模型^[15-16]组合而成。MEEMD模型能有效地抑制模态混叠，减小EEMD方法中不必要的集成平均，提高计算速度，减少无物理意义的IMF(intrinsic mode function)分量，减小由添加白噪声引起的重构误差，保证分解的完备性^[14]，对非平稳数据的处理效果较好，但端点效应仍未得到抑制；ARMA模型简单，短期预测效果较好，但主要针对平稳数据，且在极值点附近的预测精度不高。为了抑制端点效应，提高极值点附近的预测精度，本文引入SVM(support vector machine)对数据两端进行延拓^[17]，并借鉴文献[9]中的残差修正方法。具体步骤如下：

1) 提取原始数据X(t)的线性趋势项，得到趋势项X_q(t)和非趋势项X_d(t)；

2) 利用SVM对X_d(t)进行两端延拓得到X_c(t)；

3) 对X_c(t)进行MEEMD分解^[14]；

4) 将3)中分解得到的分量延拓部分去掉，得:

$ \begin{array}{c} \left\{\begin{array}{l}{\boldsymbol{X}_{t}=\left\{\overline{\boldsymbol{I}}_{j}, j=1, 2 \cdots k\right\}} \\ {\overline{\boldsymbol{I}}_{j}=\left\{x_{t}, t=0, 1, 2 \cdots\right\}}\end{array}\right.\\ {x_{t}-\varphi_{1} x_{t-1}-\varphi_{2} x_{t-2}-\cdots-\varphi_{p} x_{t-p}=} \\ {\varepsilon_{t}-\theta_{1} \varepsilon_{t-1}-\theta_{2} \varepsilon_{t-2}-\cdots-\theta_{q} \varepsilon_{t-q}}\end{array} $

(5)

式中，k为分量个数；X_t为去延拓分量矩阵；I_j为去延拓分量；x_t为分量中t时刻对应的元素；p、q为ARMA模型阶数；ε={ε_t, t=0, ±1…}是方差为σ²的白噪声序列；θ、φ为待定系数^[15]。

采用AIC(Akaike information criterion)准则作为ARMA模型定阶准则，AIC数学式如下：

$ \begin{array}{c}{\operatorname{AIC}(p, q)=\min\limits _{0 \leqslant m, n \leqslant L}(m, n)=} \\ {\min\limits _{0 \leqslant m, n \leqslant L}\left\{\ln \hat{\sigma}^{2}+\frac{2(m+n)}{N}\right\}}\end{array} $

(6)

式中，$\hat{\sigma}^{2}$是对噪声项方差的估计；L为预先给定的最高阶数；N为已知观测数据样本数。利用AIC准则定阶是指在p、q一定变化范围内寻求使得统计量AIC(p, q)达到最小的点(p, q)，用作(p, q)的估计^[15]。理论上L取值越大，得到的预测精度越高，但考虑到计算时间以及L增大时预测精度提高不太明显，所以一般情况下L取值为N/10、ln(N)或$\sqrt{N}$。

5) 利用线性趋势项X_q(t)作为回归模型建模数据，预测趋势项，得到预测值$\tilde{\boldsymbol{X}}_{q}(t)$；

6) 计算最终预测值$\tilde{\mathit{\boldsymbol{{X}}}}(t)$：

$ \tilde{\boldsymbol{X}}(t)=\tilde{\boldsymbol{X}}_{q}(t)+\sum\limits_{j=1}^{k} \tilde{\overline{\boldsymbol{I}}} $

(7)

1.2.2 探测界限的确定以及探测实现

震前电离层TEC异常扰动, 不仅要计算参考背景值，而且要确定比较合适的界限值。通常称2或3倍中误差为极限误差，因此，本文以2倍预测值残差标准差作为异常判断限值。首先，选取震中或者震中位置附近震前第15~90天中不受太阳和地磁活动影响的28 d TEC数据，利用前20 d数据，用MARCM预测后8 d，并计算预测值与真实值的残差，统计百分比在95.5%以上的残差值，作为剔除条件Δ；其次，由于ARMA模型预测长度越，长精度越低，因此预测过程分两次进行；提取震前35 d以及地震当天的TEC数据，选择震前第16~35天这20 d的电离层TEC数据为建模数据，用MARCM预测震前第8~15天的TEC背景值，计算预测值与真实值的残差，判断残差绝对值是否大于剔除条件|Δ|，将残差绝对值大于|Δ|的预测值剔除，并利用插值方法补齐数据，将震前第16~35天这20 d的电离层TEC数据与补齐数据合并成一个新的序列，利用新序列作为建模数据，用MARCM预测震前第1~7天以及地震当天这8 d的TEC值；最后，以这16 d的预测值作为震前电离层TEC扰动探测的参考背景值，并计算预测值与真实值的残差标准差σ，以2σ为参考背景值的上下界限，表达式为：

$ \left\{\begin{array}{l}{u=Y_{\text { multi }}+2 \sigma} \\ {l=Y_{\text { multi }}-2 \sigma}\end{array}\right. $

(8)

式中，Y_multi为MARCM预测的TEC参考背景值；u、l分别为背景值范围的上下界限，若观测值超过界限，则视为存在异常扰动；异常值Δ_TEC计算如式(9)所示：

$ \varDelta_{\mathrm{TEC}}=\left\{\begin{array}{l}{0, l \leqslant Y_{\text { real }} \leqslant u} \\ {Y_{\text { real }}-u, Y_{\text { real }}>u} \\ {Y_{\text { real }}-l, Y_{\text { real }}<u}\end{array}\right. $

(9)

式中，Y_real为TEC观测值。

2 MARCM探测与传统方法探测精度对比

本文提取由IGS中心^[18]提供的2013-02-01-2013-02-28(30.0°N，105.0°E)的TEC数据进行分析，分别采用MARCM、ARMA模型、IQR和滑动时窗法对本时间段电离层TEC背景值进行预测，并利用平均残差Δ、残差标准差σ以及平均相对精度P作为精度衡量指标，表达式如下：

$ \left\{ \begin{align} & \varDelta =\sum\limits_{i=1}^{n}{\left( {{I}_{{\text{pre}}[i]}}-{{I}_{\text{IGS}[i]}} \right)}/n \\ & \sigma =\sqrt{\sum\limits_{i=1}^{n}{\left( {{I}_{\text{pre}[i]}}-{{I}_{\text{IGS}[i]}}-\varDelta \right)}/n} \\ & P=\frac{\sum\limits_{i=1}^{n}{\left( 1-|{{I}_{\text{pre}[i]}}-{{I}_{\text{IGS}[i]}}|/{{I}_{\text{IGS}[i]}} \right)}}{n} \\ & i=1, 2\cdots n \\ \end{align} \right. $

(10)

式中，I_pre为预报值；I_IGS为IGS中心发布的TEC值；n为数据个数。

图 1给出了上述4种方法对2013-02-21-2013-02-28电离层TEC背景值的预测结果，可以看出，IQR与滑动时窗法预测的背景值基本一致，且随着时间变化，这两种方法每天同一时刻的预测值并未发生较大变化，系统偏差较大，不能较好地体现不同时间段电离层TEC的变化特点，因而导致电离层TEC异常扰动探测不准确；而MARCM和ARMA模型的预测值虽然在一些时刻与真实值残差较大，但总体上均能较好地体现不同时间段电离层TEC的变化特点；就MARCM和ARMA模型而言，在极大值附近的预测精度相当，而在极小值附近，MARCM预测精度明显比ARMA模型高。

图 2给出了2013-02-01-2013-02-28地磁和太阳活动平静状态下3种方法的探测结果，可以看出，本文提出的探测方法未发现明显的异常探测结果，该方法总体要优于其他两种方法，进一步说明了新模型与方法探测的有效性。

表 1为4种探测方法每天的平均残差、残差标准差以及相对精度统计。总体上，MARCM、ARMA模型、IQR及滑动时窗法的平均相对精度分别为89.78%、84.60%、88.37%及88.36%，MARCM明显优于其他3种方法；从标准差来看，MARCM比IQR和滑动时窗法略大，而ARMA模型最大，但4者相差不大；由平均残差可以看出，MARCM每天平均残差绝对值均在1 TECU以内，明显优于其他3种方法。

为进一步说明MARCM对参考背景值的预测精度，表 2给出了4种方法残差绝对值范围分布百分比。由表 2可知，MARCM、ARMA模型、IQR以及滑动时窗法预报的残差绝对值小于3.0TECU的分别占91.67%、72.92%、84.38%和85.42%，绝对残差小于5.0 TECU的分别占98.96%、97.92%、97.92%和97.92%，而MARCM和ARMA模型的绝对残差均在6.0 TECU以内，说明本文模型预测精度明显高于其他3种方法。

因此，利用MARCM预测电离层TEC参考背景值来探测震前电离层扰动是可行的。

3 四川雅安市芦山县Mw7.0级地震电离层异常探测

选取北京时间2013-04-20 T 08:02:46中国四川省雅安市芦山县Mw7.0级地震(震中30.3°N，103.0°E，震深13 km)为研究对象，提取2013-03-16-2013-04-20 IGS中心提供的电离层网格数据，并通过双线性插值取得震中位置的TEC值，利用§1.2.2中提到的探测方法分析2013-04-05-2013-04-20这16 d电离层TEC的异常情况。为了方便分析，文中将地方时转换成世界时(UT=LT-8)。

图 3给出了芦山地震前15 d及地震当天表征地磁活动的Kp指数和DST指数以及表征太阳活动的太阳辐射流量F_10.7。一般认为Kp < 3时，无磁暴现象；|DST|>50 nT时，可能发生了中大磁暴；F_10.7指数根据大小分为3个等级，即F_10.7>150 SFU，100 SFU < F_10.7≤150 SFU，F_10.7≤100 SFU，对应强中弱3个等级的太阳活动^{[1, 19]}。由图 3可知，4月5日-20日F_10.7指数除4月19日为99 SFU外，其他15 d均超过100 SFU，其中4月9日和10日接近150 SFU，其他13 d的F_10.7指数变化相对稳定，这说明9日和10日的电离层扰动可能与太阳活动有关；而4月5日-20日电流强度DST的绝对值均在30 nT以内，但14日6:00 UT和9:00 UT的Kp指数超过3，表明该日有一定的地磁干扰，因此，该日电离层TEC扰动可能是地磁扰动引起的；除上述4月9日、10日及14日之外，其他13 d内Kp指数和DST指数均处于较低水平，且F_10.7指数变化相对稳定，表明这13 d内如果存在电离层TEC异常，则极有可能是芦山地震导致的。

图 4给出了MARCM、IQR以及滑动时窗法对震中位置的探测结果。由图 4可看出，在4月5日-14日IQR和滑动时窗法探测到的正异常明显多于MARCM，且异常一旦出现，则随后数天都会出现异常现象，而滑动时窗法尤为明显，与文献[1]、文献[11]中相关论述相符，这主要是由于传统探测方法参考背景值预测精度较低及算法系统偏差导致的；4月15日-20日内，3种探测方法均探测到了明显的负异常。

由于传统探测方法存在一定缺陷，探测不准，因此，本文主要分析MARCM的探测结果。从图 4可以看出，MARCM除了探测到4月9日、10日及14日的TEC正异常外，还探测到4月7日和11日的TEC正异常以及4月16日-20日明显的负异常，结合图 3分析结果，说明4月7日、11日的正异常以及4月16日-19日的负异常极可能是孕育2013-04-20芦山Mw7.0级地震而引起的。

为进一步分析芦山地震震前电离层TEC异常情况，分别用MARCM和IQR探测(32.5°N，105.0°E)、(30.0°N，105.0°E)、(32.5°N，100.0°E)和(30.0°N，100.0°E)4个位置震前15 d及地震当天电离层TEC的异常情况。图 5给出了探测结果，由图中MARCM探测结果可知，4个位置探测的结果与震中位置基本一致，但震中往南((30.0°N，105.0°E)、(30.0°N，100.0°E))震前5 d TEC负异常在时间范围上的分布比其他两位置略大，且与震中位置异常分布更加接近；而震中往东((32.5°N，105.0°E)、(30.0°N，105.0°E))的正异常在数值上略小于其他两位置，且与震中位置正异常在数值上接近；说明芦山地震震中往南的TEC负异常分布更能表征震中位置，而东面的正异常在数值上更能表征震中位置正异常值。

为进一步分析芦山地震震前15 d电离层TEC异常在时间上的分布情况。图 6给出了震中以及震中附近4个方向震前15 d的电离层TEC异常。由图 6可知，震前第6~15天有明显的TEC正异常，主要出现在08:00 UT和10:00 UT，且在10:00 UT时刻，正异常达到最大；而震前第1~4天出现明显的负异常，异常主要出现在00:00-14:00 UT之间，且在12:00 UT时刻，负异常值相对较大；震中往南方向((30.0°N，105.0°E)、(30.0°N，100.0°E))负异常分布与震中位置相似，一定程度上验证了图 5分析的相关结论。

4 结语

本文提出多尺度ARMA残差修正模型，对比分析了该模型与其他3种算法对参考背景值的预测精度，结果表明，MARCM参考背景值预测精度明显优于ARMA模型、IQR和滑动时窗法。还分析了2013-04-20四川雅安市芦山县Mw7.0级地震震前15 d及地震当天的地磁和太阳活动指数，表明了4月9日、10日及14日电离层扰动可能不是地震引起。震中位置附近4个方向的电离层TEC异常，说明芦山地震震中偏南方向的TEC负异常分布更能表征震中位置，而东面的正异常在数值上更能表征震中位置正异常值。

然而，本文采用IGS提供的电离层格网数据作为TEC真值，IGS TEC在全球陆地的精度约5 TECU，由于IGS TEC在建模时利用的中国区域的GNSS测站数量较少，故IGS TEC在中国区域的精度偏低，这在一定程度上影响了TEC异常判定依据构建模型的精度。因此在后续研究中，将利用中国地区密集的地面GNSS区域观测数据重建局部TEC值，从而构建更高精度的预测模型和进行定量分析地震相关研究。同时，由于地震的产生机理复杂，且目前还没有相对完整的理论支撑震前电离层扰动与地震有关，因此，利用电离层扰动与地震相关性还需进一步研究。