近年来，随着全球定位系统（Global Positioning System，GPS）等卫星星座的不断完善，基于卫星的导航、定位、通讯以及各种科学研究越来越普遍。而电离层对电磁波的折射、散射、反射和吸收作用对卫星信号的传播影响巨大，其中电离层垂直总电子含量（vertical total electron content，VTEC）是决定电离层延迟的重要参数，对其进行预报研究具有重要的现实意义^[1-8]。传统经验模型对电离层VTEC的估计精度仅有50%~60%，且无法实时提供延迟改正^[9]，因此有学者尝试将时间序列作为一种电离层短期预报方法，并取得了一些有益结论。其中最早被引入电离层预报领域的自回归模型（auto regressive model，AR）虽然待估参数较少、建模简单，但未能较好考虑电离层VTEC的季节性变化，预报精度较低^[10]；在其基础上发展的自回归积分滑动平均模型（auto regressive integrated moving average model，ARIMA）较为全面地考虑了VTEC的短期趋势与周期变化，在平静期的不同纬度地区相对精度可达80%以上，有效提高了预报值的可靠性^[11-12]。然而，电离层受日地相对距离、太阳的周期性活动等因素影响较大^[13]，在太阳活动的活跃期，受太阳风高能粒子影响使得电离层总电子含量的周日变化较平静期更为显著^[14-15]，研究表明，采用同一方法对VTEC进行预报时，活跃期预报值的残差较平静期明显增大^[16-17]。当磁暴发生时，伴随地磁水平经历初相、主相和恢复相的过程，电离层VTEC同样会产生异常扰动^[18]，其非线性与非平稳性特征增强，大幅降低了预报精度。

为此，本文将小波分解与时间序列分析法相结合改进了ARIMA模型，并分别采用16个GPS观测站的实际VTEC解算值和全球电离层模型（global ionosphere model，GIM）格网数据，利用长度为10 d的样本序列进行预报，验证改进模型的预报精度与适用性。

1 小波分解改进ARIMA模型 1.1 小波分解基本原理

非平稳时间序列随时间的变化可以分为慢变和快变两部分，慢变对应时间序列的低频部分，描绘了序列的主体信息；快变对应时间序列的高频部分，表示序列的细节纹理。小波分解在局部范围内具有伸缩、平移和放大功能，对信号可以进行多尺度分析，在处理非平稳时间序列时具有一定的优势。小波分解与重构采用的是Mallat塔式算法，分解过程的表达式如下：

$ \left\{ {\begin{array}{*{20}{l}} {a{'_{i + 1}} = H{a_i}}\\ {d{'_{i + 1}} = G{d_i}} \end{array}} \right. $

(1)

式中，H为低通滤波器；G为高通滤波器；a_i、d_i分别为高频分量与低频分量；$ a{'_{i + 1}}$和$ d{'_{i + 1}}$分别是1/2减采样后得到的低频与高频系数。通过对分解后得到的低频系数与高频系数进行重构，得到与原始序列长度相等的低频序列a_i+1与高频序列d_i+1，重构过程的Mallat算法表达式为：

$ \left\{ {\begin{array}{*{20}{l}} {{a_{i + 1}} = {H^{\rm{*}}}a{'_{i + 1}} + {G^{\rm{*}}}d{'_{i + 1}}}\\ {{d_{i + 1}} = {x_i} - {a_{i + 1}}} \end{array}} \right. $

(2)

式中，H^*为重构低通滤波器；G^*为重构高通滤波器，分别与H、G互为对偶算子。按上述算法对原始序列{a}进行i级分解后得到的各级分量与原始序列的关系为：a=a_i+d_i+d_i-1+…+d₁。

1.2 ARIMA模型及改进

利用ARIMA模型预报电离层VTEC时需要先判断样本序列的平稳性，若序列不平稳，则需要先对其进行平稳性处理，得到平稳序列后根据样本序列表现出的自相关与偏自相关特性选取合适的模型参数完成建模。VTEC序列是季节性与短期相关性显著的时间序列，故本文采用季节乘积模型ARIMA(p, d, q)×(P, D, Q)进行预报实验，该模型的完整表达为：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\Phi }\left( B \right){\mathit{\Phi }_S}\left( B \right){\nabla ^d}\nabla _S^D{x_t} = \mathit{\Theta }\left( B \right){\mathit{\Theta }_S}\left( B \right){\varepsilon _t}}\\ {{\mathit{\Phi }_S}\left( B \right) = 1 - {\varphi _1}{B^S} - \ldots - {\varphi _p}{B^{PS}}}\\ {{\mathit{\Theta }_S}\left( B \right) = 1 - {\theta _1}{B^S} - \ldots - {\theta _q}{B^{QS}}} \end{array}} \right. $

(3)

式中，S为季节性周期；$ \nabla _S^D = {(1 - {B^S})^D}$是以周期S为步长的季节差分过程；D为差分阶数；${\mathit{\Phi }_S}\left( B \right) $是季节性自回归算子；P为季节性自回归阶数；${\varphi _1}, {\varphi _2} \ldots {\varphi _p} $为季节性自回归参数；${\mathit{\Theta }_S}\left( B \right) $为季节性滑动平均算子；Q为季节性滑动平均阶数；$ {\theta _1}, {\theta _2} \ldots {\theta _q}$为季节性滑动平均参数。

利用小波分解能够快速高效地将VTEC序列的细节成分进行分离，使得样本序列的短期趋势与周期变化更为显著，从而提高预报精度。但由于预报结果是通过对小波分解后各分量的预报值重构而获得，过多的分解层数会导致预报误差的积累。以往的研究通常采用1~3级分解对样本序列进行处理^{[7, 15]}，本文通过实验证明一级分解后的预报结果精度最高，同时操作简单，避免了反复建模过程中误差的累积。通过小波分解改进的ARIMA模型建模流程如下：

1）利用db4正交小波对样本序列进行一级分解，得到长度减半的低频系数与高频系数，分别对其重构得到与原始序列等长的低频与高频分量；

2）判断各分量是否平稳，若不平稳，则进行差分处理，直到平稳后根据其表现出的自相关与偏自相关特性确定模型的其他参数；

3）分别对低频序列和高频序列进行预报并重构预报值，得到最终预报结果。

2 基于多测站观测值的VTEC预报 2.1 数据来源与地磁活动分析

本文选取中国境内的16个GPS观测站于2011年年积日66~78天（3月7日至3月19日）的观测数据，采用载波相位平滑伪距法（参见文献[1]）以15 min为间隔分别对单站VTEC进行解算，为方便实验，以1 h为间隔提取每天24个VTEC值构成样本序列。为更好地研究活跃期电离层的变化特性以及预报值均方差的分布情况，本文主要研究范围为75°E~125°E，17.5°N~45°N，使得测站能够尽可能均匀分布于实验区域内。

为考察预报期实验区域内的地磁扰动强度，图 1给出了预报期内扰动暴实时指数（disturbance storm time index，Dst）与K_p指数。Dst指数表征环电流强度，时间分辨率为1 h，而K_p指数用以描述每3 h内的平均地磁强度。由图 1可知，在太阳活动剧烈时，受太阳风高能粒子吹拂，地球自东向西的环电流强度会大幅增强，其电磁场与地磁场反向，故Dst指数会下降；当Dst指数低于-30 nT即可能发生小磁暴，小于-50 nT时有可能发生中等磁暴。实验期间内3月11日，多个时刻的Dst指数小于-50 nT，当天与之相对应的K_p指数多个时刻达到5或6（发生中小磁暴），当日K_p指数之和大于30，故可以确定当日地磁活动剧烈且有磁暴发生。

为分析磁暴对电离层总电子含量的影响，选取实验区域内不同经纬度的3个测站（BJSH（116.20°E，40.20°N）、WUHN（114.50°E，30.50°N）、QION（109.84°E，19.03°N））解算VTEC值，结果如图 2所示。从图 2中可以看出，3个测站的VTEC的总体变化趋势随Dst指数的减小和K_p指数的增大而增大，在磁暴发生当日VTEC达到最大值，受磁暴影响，VTEC序列周日变化幅度较大且总体平稳性较差。

可见，VTEC值随纬度的降低总体增大，其残差也呈现增大趋势。ARIMA及其改进模型在每日VTEC的极大值处预报误差较大，在极小值处误差较小，改进模型的预报值残差明显减小，尤其在每日VTEC的极大值点附近能够较大程度地削弱残差峰值，较原模型具有更高的预报精度。

2.2 模型精度及适用性分析

单站解算得到的VTEC经过筛选提取后构成的时间序列时间分辨率为1 h，为使实验过程中原模型与改进模型的样本序列最大程度地受到磁暴的影响，以年积日66~75天（前10 d）的VTEC作为样本序列，对年积日76~78天（后3 d）的VTEC值进行预报，并将预报值与实际解算值进行比较分析。由于测站数量较多，本文选取6个测站（DXIN、BJSH、ZHNZ、WUHN、LHAS、QION）的预报结果与预报值残差（Δ）进行分析，结果如图 3所示，其他测站表现了相同的特性。

进一步分析上述6个测站的预报值残差并进行统计，表 1给出了以1、2、3 TECu为节点的不同区间内的残差比例。可以看出，原模型与改进模型在纬度较低测站的残差要大于高纬度的测站，原模型对纬度较低的QION测站的预报结果中残差小于1 TECu的仅有8.33%，而大于3 TECu的达到了70.83%；即使在预报效果较好的DXIN测站，小于1 TECu的残差值也仅有45.83%；改进模型较ARIMA模型的预报精度有明显提升，其中小于1 TECu的残差比例最高可达63.89%，而预报结果不理想的QION测站改进后的模型也较大程度地提高了小于1 TECu的残差比例，降低了大于3 TECu的残差比例。

图 4给出了原模型和改进模型在实验区域内16个测站3天预报值各1 152个残差的分布情况，可以直观看出改进模型的残差较原模型明显减小。通过统计分析得到ARIMA模型与改进模型置信度为95%的残差区间分别为（2.77±8.56）TECu和（1.60±6.89）TECu，就残差分布的统计结果而言，改进模型的预报精度要高于原模型。

均方根误差（root mean square error，RMSE）能够较好地反映预报值的可靠性^{[7, 14-15]}，本文采用实验区域内RMSE的大小来评定模型的预报精度。同时，为论证改进模型在实验区域内是否具有良好的适用性，分别利用ARIMA模型和改进模型对区域内的16个测站进行VTEC预报，通过预报值分别求取年积日为76~78天的每日RMSE，统计值见表 2，图 5描绘了区域内3 d的RMSE分布情况。从图 5中可以直观看出，预报过程中第2日（年积日77天）的RMSE最大，预报精度最低，原模型在纬度较低的局部地区峰值可达12 TECu左右，而改进模型能将该区域内的RMSE峰值有效降低约4 TECu；此外，预报的第1日和第3日原模型预报值在大部分区域的RMSE均在5 TECu以内，最大峰值不超过9 TECu，与之对应的改进模型预报值在大部分区域的RMSE均在3 TECu左右，且最大峰值不超过6 TECu。可见，在实验区域内，改进模型预报结果的RMSE较原模型总体降低，具有更高的预报精度和良好的区域适用性。

3 基于GIM的VTEC预报

GIM是欧洲定轨中心（Center for Orbit Determination in Europe, CODE）通过选取IGS分布在全球范围内的150个GPS观测站的观测数据，采用15阶球谐函数拟合的适用于全球范围的电离层模型，该模型能够有效反映全球范围的电离层特性。GIM数据的空间分辨率为2.5°×5°（纬度×经度），时间分辨率为2 h，一天共提供13幅全球电离层图像。为方便建模，本文选取单个格网点上00：00~22：00间隔为2 h，每天共12个VTEC数据构成样本序列（年积日的选取与§2实验一致）。为使实验区域保持一致，选取75°E~125°E(间隔10°)、17.5°N~45°N（间隔7.5°），共30个格网点，分别利用ARIMA与改进模型进行3 d预报。

利用§2.2中方法对30个格网点上1 080个预报值的残差进行统计分析，得到原模型和改进模型置信度为95%的残差区间分别为(3.01±9.23) TECu和(1.81±8.36) TECu。表 3给出了实验区域内30个格网点上3 d预报值的RMSE均值，除个别格网点外，改进模型的RMSE均小于原模型。实验区域内各格网点上每天的RMSE分布见图 6，总体反映了实验区域内ARIMA及其改进模型预报值RMSE的分布情况。与以上实验结果相似，预报过程中RMSE最大峰值出现在预报第2日的较低纬度区域，原模型在该区域的误差峰值可达14 TECu左右，而改进模型降低了约3 TECu；在预报的第1日与第3日，两种模型的预报效果均比第2日理想，单日的RMSE峰值也并不明显；但从图 6上可以直观反映出改进模型在这两日中减小了区域范围内的RMSE。因此，改进模型较大程度地提高了活跃期VTEC单点预报的精度，且在实验区域内具有良好的适用性。

4 结语

本文分析了太阳活跃期磁暴发生前后地磁指数与电离层总电子含量变化的关系，利用db4正交小波分解样本序列改进ARIMA模型，并通过测站实际解算的VTEC数据与IGS提供的GIM分别验证了改进模型的预报精度，得到以下结论。

1) 在磁暴发生过程的前后，Dst指数经历了先增大后减小再增大，最终趋于平稳的过程，其中磁暴发生时Dst指数达到最小值；而K_p指数表现为先减小后增大再减小，最终趋于平稳，其中磁暴发生时K_p指数达到最大值；对应时刻的总电子含量同样在磁暴发生时达到最大值，该影响对低纬度地区最为显著，随纬度升高，影响逐渐减弱。

2) 利用db4小波分解改进的ARIMA模型在活跃期的预报结果残差较原模型明显减小；通过不同数据源分析实验区域内预报值RMSE的分布情况，结果表明，改进模型的预报精度要高于原模型，且在实验区域内具有良好的适用性。

另外，电离层除受到太阳活动影响外，还会受到日地相对距离、磁层、热层、中层大气的物理化学等变化的影响，对流层的雷暴等天气亦会对电离层底部电子密度产生影响，单层模型无法细致反映这些因素引起的电离层变化。因此，在后续研究中，可以考虑利用电离层的层析方法，同时结合各种电离层影响因素的观测资料，进行更加详细、可靠的相关性分析，进一步提高VTEC预报精度。