地极移动是表征地球自转运动的一个重要参量，简称极移（polar motion，PM）。极移的两个分量x_p、y_p与日长变化统称为地球自转参数（Earth rotation parameters，ERP），极移是实现天球参考系与地球参考系相互转换的必需参数^[1]。甚长基线干涉（very long baseline interferometry，VLBI）和全球卫星导航系统（Global Satellite Navigation System，GNSS）等现代空间测地技术对极移的测量精度可达0.05 mas，分辨率可达1~2 h^[2]。但由于复杂的资料处理过程，现代空间测地技术获取的极移结果存在一定的延迟，而卫星导航和深空探测等领域对极移的实时测量值具有重要需求，因此对极移进行短期预报具有重要的现实意义。此外，由于极移与多种地球物理现象密切相关，这使得极移的中长期预报具有重要的科学意义^[3]。

当前预报极移的方法有多种，如最小二乘（least squares，LS）外推与自回归(autoregressive，AR)模型的组合（简称LS+AR）^[4-6]、LS外推与神经网络（neural network，NN）模型的组合（LS+NN）^[6-8]、小波和模糊推理系统（fuzzy-inference system，FIS）的组合^[9]。全球性的地球定向参数比较竞赛（Earth Orientation Parameters Prediction Comparison Campaign，EOP PCC）的结果表明，LS+AR模型是当前极移预报精度最高的一种方法^[10]。该方法是针对极移序列中的钱德勒摆动（Chandler wobble，CW）与周年摆动（annual wobble，AW）等确定性趋势项进行曲线拟合，建立一个LS外推模型，从而获得极移趋势项的外推值和拟合残差序列，然后利用AR模型对残差序列进行预报，而极移预报值为趋势项外推值和残差项预报值之和。由于极移序列中不仅包含CW与AW项等低频分量，还含有高频分量，因此不宜直接对极移序列进行曲线拟合，即在进行曲线拟合之前，最好将极移序列中的高频分量去除。王小辉等^[7]和Shen等^[11]分别将经验模式分解和奇异谱分析应用于极移高频分量与低频分量分离，有效提高了极移预报精度。本文将小波分解（wavelet decomposition，WD）引入到极移预报中，首先利用WD对极移序列进行多尺度分解，将高频分量和低频分量分离，然后对低频分量建立LS外推模型，获得趋势项外推值和残差序列，最后采用AR模型对高频分量和残差序列之和建模预测。实验表明，基于小波数据预处理方法的极移预报精度得到了显著提高。

1 模型建立 1.1 小波分解

极移包含长期趋势变化、CW、AW、半周年摆动以及更高频尺度上的变化等。本文提及的极移低频分量包含钱德勒摆动、周年摆动和其他长期趋势变化，剩余的季节性、亚季节性以及更高频尺度上的变化称为高频分量。具体而言，振荡周期大于1 a的极移分量称为低频分量，而振荡周期小于1 a的极移分量称为高频分量。由小波分析理论可知，信号可以通过小波分解一层一层分解到不同的频率通道上，分解后的信号在频率成分上比原始信号单一，这些信号按照频率高低依次排列。因此可以通过对分解后的某些信号重构实现原始信号低频分量和高频分量的分离。

对极移x_p、y_p分量的时间序列X、Y进行多尺度一维离散小波分解，其中分解层数j > 1，可以获得各层的近似小波分解系数和细节小波分解系数。对第j层的近似小波分解系数和各层的细节小波分解系数分别进行重构，可以得到一个近似信号A_j和一组频率从高到低排列的细节信号D₁D₂⋯D_j，则可以将极移x_p分量的序列X表示为：

$ \mathit{\boldsymbol{X}} = {\mathit{\boldsymbol{A}}_j} + \sum\limits_{i = 1}^j {{\mathit{\boldsymbol{D}}_i}} $

(1)

由于极移x_p、y_p分量的预报原理相同，本文主要以x_p分量叙述为主，故仅给出序列X的小波分解表达式，下同。

将前k层的细节信号D₁D₂⋯D_k求和，可以获得极移高频分量X'_H，即：

$ \mathit{\boldsymbol{X}}_H^\prime = \sum\limits_{i = 1}^k {{\mathit{\boldsymbol{D}}_i}} $

(2)

易知，极移低频分量X'_L可以表示为：

$ \mathit{\boldsymbol{X}}_L^\prime = \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{X}}_H^\prime = \mathit{\boldsymbol{X}} - \sum\limits_{i = 1}^k {{\mathit{\boldsymbol{D}}_i}} $

(3)

1.2 LS外推模型

对于极移低频分量X'_L、Y'_L，LS外推模型包含3个部分：长期趋势项、AW与CW，模型的数学表达式为:

$ \mathit{\boldsymbol{X}}_L^\prime (t) = \sum\limits_{i = 1}^2 {{c_i}} \cos \left( {2\pi /{p_i}} \right) + {d_i}\sin \left( {2\pi /{p_i}} \right) + a + bt $

(4)

式中，a、b表示低频分量X'_L模型的长期趋势项参数；c₁、d₁表示模型的CW参数；c₂、d₂表示模型的AW参数；p₁、p₂分别表示CW周期和AW周期，取p₁=1.183a，p₂=1a^[5]。Y'_L模型的数学表达式与X'_L类似。

1.3 AR模型

AR模型是对平稳序列$\left\{ {z\left( t \right), t = 1, 2 \cdots n} \right\} $建立的一个概率统计模型，描述了序列z(t)与t时刻以前的规律变化及t时刻的白噪声之间的关系，其数学模型为：

$ z(t) = \sum\limits_{i = 1}^p {{\varphi _i}} z(t - i) + \varepsilon (t) $

(5)

式中，p为模型阶次；ε(t)为白噪声；φ₁, φ₂⋯φ_p为模型参数，可以通过求解Yule-Walker方程来确定^[12]；序列z(t)为极移高频分量与最小二乘拟合残差之和。

模型阶次p直接影响着AR模型性能。本文采用应用相对广泛的AIC（Akaike information criterion）准则来选取阶次p：

$ {\rm{AIC}}(p) = \ln \sigma _p^2 + \frac{{2p}}{n} $

(5)

式中，$ \sigma _p^2$为模型拟合残差的方差，当p时，使AIC(p)取最小值的p就取作AR模型的阶次，一般可将p的最大搜索范围确定为1~ $ \sqrt n $。

综上所述，基于小波预处理的LS+AR模型的极移预报流程如图 1所示, 图 1中“+”表示求和。

2 模型验证 2.1 数据来源

本文采用的极移数据来自国际地球自转与参考系服务（International Earth Rotation and Reference Systems Service，IERS）发布的EOP 08 C04序列，时间跨度为2000-01-01—2016-03-22，采样间隔为1 d。

2.2 精度评定标准

为了衡量极移预报结果的精度，选用平均绝对误差（mean absolute error，MAE）作为精度评定标准，其计算公式为^[9]:

$ {\rm{MA}}{{\rm{E}}_l} = \frac{1}{N}\sum\limits_{i = 1}^N {\left| {P_i^{(l)} - O_i^{(l)}} \right|} $

(6)

式中，l为预报跨度；N为预报期数；$P_i^{(l)} $、$ O_i^{(l)}$分别为i点的预报值和观测值。

2.3 结果分析

首先利用小波分解对极移x_p、y_p分量的时间序列进行分解。小波基函数与小波分解层数和的选取目前尚无明确的理论作为指导，应根据信号变化的情况进行恰当选择。过高的小波分解层数会对信号的确定性的波动现象进行过滤，引起原信号重构失真; 而过低的小波分解层数无法将信号中的不同频率特征的分量分解出来。经过对比分析，本文选取db6小波对2000-01-01—2009-12-31期间的极移x_p、y_p分量的时间序列进行分解，分解层数取7。限于篇幅，这里仅给出x_p分量的序列分解结果，如图 2所示。从图 2可以明显看出，细节信号D₇表现为周期摆动等季节性变化，近似信号A₇表现为年季等长期趋势变化，D₁，D₂⋯D₆信号表现为亚季节性或其他更高频尺度上的变化。将D₁，D₂⋯D₆细节信号重构构成极移高频分量，如图 3(a)所示，将剩余的细节信号D₇和近似信号A₇重构构成极移低频分量，如图 3(b)所示。为了检测极移序列及其低频分量中的主频率成分，利用快速傅里叶变换（fast Fourier transform，FFT）对二者进行频谱分析，结果如图 4所示。从图 4可以看出，极移原始序列及其低频分量均包含CW和AW，说明低频分量中的主频率成分与原始序列一致。

利用LS+AR模型和基于小波预处理的LS+AR模型（WD+LS+AR）对极移分别进行1~180 d跨度的预报，预报期为2010-01-01—2016-03-22，每隔7 d预报一次，总共进行了300期的预报，其中，基础序列的长度为10 a。图 5给出了两种模型在不同跨度的MAE，表 1统计了WD+LS+AR模型相对于LS+AR模型的预报精度改善情况。

从图 5和表 1可以发现，WD+LS+AR模型相对于LS+AR模型的预报精度在各个跨度都有一定程度的提高，这是由于将小波分解应用于极移数据预处理，相当于滤波后进行最小二乘拟合获得趋势项和主周期项，而滤波后的极移低频分量与最小二乘外推模型的匹配性更高，因而能够改善预报精度。从表 1可以看到，对于跨度为1~60 d的预报，WD+LS+AR模型的预报精度提高在10%以内；从第60天开始，WD+LS+AR模型的预报精度提高越来越显著，最大提高达32%，这表明WD+LS+AR模型对极移中长期预报的改善更为显著。从图 5和表 1还可以发现，WD+LS+AR模型的x_p分量的预报精度改善程度要大于y_p分量。究其原因，可以认为小波分解是将时间序列平稳化的一种过程，这将提高模型预报精度，而x_p、y_p分量序列的平稳性存在差异，因而小波分解对x_p、y_p分量预报结果的改善程度也就不同，简言之，小波分解对预报结果的改善程度是与极移观测数据特征相关的。

3 结语

本文将小波分解应用于极移中长期预报。分析结果表明，小波分解能够很好地分离极移序列中的低频分量和高频分量，基于小波预处理的LS+AR模型的极移预报精度在不同跨度均有一定程度改善，特别是在跨度为中长期时提高更为明显，证明了极移预报中小波预处理的必要性和有效性。

在后续的研究工作中，将进一步考虑大气与海洋等极移激发源的影响，在预报模型中加入这些激发因素对极移进行联合预报。