0 引 言
针对极移预报问题，国内外学者开展了广泛深入的研究，提出了多种预报模型与方法。早期研究中，Chao[9]率先利用最小二乘法(Least Squares，LS)提取极移序列中的长期趋势与周期项，构建了外推模型进行极移的预报。随后，Kosek等[10]在LS模型基础上引入自回归(AutoRegressive，AR)技术，对LS拟合后的残差序列进行随机建模。凭借结构简单且在短期预报中表现稳健的特点，LS+AR组合模型逐渐成为极移预报领域的经典范式[11]。此后，许多学者对该范式进行了优化与改进，如张昊等[12]通过引入加权策略构建加权最小二乘(Weighted Least Squares，WLS)模型，有效增强了中长期预报的稳定性。随着人工智能技术的快速发展，其处理复杂非线性映射的强大能力为极移预报注入了新活力。Schuh等[13]最早探索了人工神经网络在极移预报中的应用。近年来，长短期记忆网络等深度学习模型也被使用在极移的预报任务中，进一步提升了极移预报的精度[14-15]。
本文旨在研究融合地球自转动力学机制与WLS的极移预报方法。该方法首先依据地球自转动力学方程，将极移观测序列转换为对应的激发函数；继而使用WLS提取激发函数中的主要周期项，并引入时间衰减因子对历史建模数据实施动态加权，增强近期观测数据对模型的贡献度；最后，通过动力学积分将预测的激发函数转换回极移数据，完成极移序列的外推预报。该方法通过结合动力学约束与数据驱动，在确保短期预报精度的同时，提高了中长期预报的稳定性。此外，本文还进一步分析本文预报方法在航天器定轨中的应用效果，量化不同预报时长对轨道转换精度的影响，为我国航天任务中的高精度参考坐标系转换提供理论依据与技术支持。
1 极移预测的动力学机制与建模
1.1 极移的激发函数
式中，$ {p}(t)={p}_{x}(t)+i{p}_{y}(t) $ 是被定义在复数域的观测极移坐标，其坐标系根据国际地球自转和参考系服务(International Earth Rotation and Reference Systems Service，IERS)规定处于国际地球参考系(International Terrestrial Reference System，ITRS)中：极移X分量$ {p}_{x}(t) $ 正方向沿格林尼治子午线方向，极移Y分量$ {p}_{y}(t) $ 正方向指向西经90°方向；$ {\chi }(t)={\chi }_{x}(t)+i{\chi }_{y}(t) $ 是地球自转的激发函数在赤道平面内的分量，涵盖大气运动、海洋运动、陆地水文等所有可能引起地球自转变化的物理源[21]；$ {\chi }_{x}(t) $ 正方向与$ {p}_{x}(t) $ 相同，而$ {\chi }_{y}(t) $ 正方向与$ {p}_{y}(t) $ 相反，指向东经90°方向，这是由于观测极移和实际极移在转换中存在负号，从而导致两者Y分量的反向[6]；考虑到地球的黏弹性，引入复频率$ {{\sigma }}_{\text{cw}}=2{\text π} [1+i/(2Q)]/{T}_{\text{cw}} $ 表示钱德勒摆动频率与衰减阻尼，本文研究中取钱德勒周期Tcw=433天，品质因子Q=79[6]。
刘维尔方程描述了极移复杂的动力学演化过程，量化了极移的动力学机制，从而将研究视角从极移本身拓展至物理激发域。基于极移观测数据反演大地测量激发函数$ {{\chi }}_{\text{geo}}(t) $ 可以表示为[22]：
式中，$ \Delta t $ 是离散极移数据的时间步长。本文选用IERS发布的EOP 20 C04数据集，该数据集代表了当前极移观测数据的最高精度水平。如图1 所示，基于该数据集提供的日分辨率极移序列，利用式(2)反演计算出大地测量激发函数。
1.2 激发函数频率分析
利用快速傅里叶变换方法对大地测量激发函数$ {{\chi }}_{\text{geo}}(t) $ 进行频谱分析。1993年以来，随着全球定位系统等现代空间大地测量技术的引入，极移观测精度得到了显著提升，从而有效捕捉到了更为真实和丰富的高频信号[23]。本文选取1993—2024年观测数据进行分析，以避免早期数据中潜在的高频插值误差干扰。如图2 的频率谱图所示，在周期小于400天的范围内，极移的激发函数呈现显著的周期性特征，其能量主要周期集中在365天、180天及120天等季节性周期附近。这种周期性特征的物理根源，可追溯至极移的主要激发源——大气质量和海洋质量的周期性迁移，以及大气表层与海底压力变化引发的地球形变。这些动力学过程通过角动量守恒机制，将自身周期性变化传递至地球自转，最终在极移的大地测量激发函数中显示出可观测的周期性信号[24]。
1.3 激发函数的建模与外推
在获取历史时段的激发函数$ {\boldsymbol{\chi }}_{\text{geo}}(t) $ 后，根据第1.2节所述的频谱分析结果，确定了其主要周期成分(T1=365.250天，T2=182.625天，T3=121.750天)。为了从激发函数的原始序列中重建出这些主要成分，本文构建了包含线性趋势和周期项的谐波模型$ S(t) $ ：
式中，t为观测历史窗口时间，$ {\omega }_{k}=2{\text π} /{T}_{k} $ 为主要周期的角频率，a0为常数项，a1、Ck、Dk分别为线性趋势项和周期项待定系数。为了求解上述待定系数，将离散的观测方程表示为矩阵形式：
式中，Y为观测向量，A是由常数项、线性趋势项和周期项构成的设计矩阵，X为待定系数向量，V是残差向量。使用最小二乘法对原始序列进行拟合，通过最小化残差平方和，可得待定系数的最优解为：
由于极移激发函数受大气、海洋等多源驱动机制的影响，具有明显的时变性，所以近期观测数据对模型构建具有更强的信息价值。为此，本文构建了基于时间的权重函数：通过动态调整历史观测数据的权重系数，设计权重随时间衰减，赋予近期数据高权重，从而将模型的注意力转移到近期的数据中。时间权重wi为：
式中，i=1, 2, ··· , n为历史观测数据窗口中的采样编号，ti为对应的时间点。设置遗忘因子λ=0.999，其对应半衰期时长约为两年；同时设定权重下限wmin=0.5，防止历史数据被过度衰减导致序列有效长度缩短。将权重组成权重矩阵W=diag{w1, w2, ···,wn}，则WLS的参数估计为：
1.4 极移积分预测
为实现极移的动力学预测，本文构建了基于刘维尔方程的积分框架。首先推导了极移刘维尔方程的积分形式表达式，然后采用四阶龙格-库塔(Runge-Kutta，RK4)数值方法求解该动力学方程，从而建立激发函数与极移轨迹的映射关系。基于极移刘维尔方程的解析形式，将极移位移量、激发函数及钱德勒频率的复数域表达式代入动力学模型，可推导得到极移动力学微分方程组：
式中，$ {\dot{p}}_{x} $ 和$ {\dot{p}}_{y} $ 分别为极移X分量和Y分量的时间一阶导数；$ {p}_{x} $ 和$ {p}_{y} $ 分别为极移X分量和Y分量；$ {\chi }_{x} $ 和$ {\chi }_{y} $ 分别为对应轴向的激发函数分量。由该方程可知，极移的动力学演化由两项因素共同控制：一是促使极移产生旋转运动的耦合驱动项，其中激发函数提供了持续的外部输入；二是导致能量耗散的阻尼项，它反映了极移克服黏滞阻力，向激发函数所定义的瞬时平衡轴收敛的过程。
基于第1.3节求解的WLS模型对激发函数进行外推预测，获得未来n天的预测激发$ {\overset{\frown }{{\chi }}}(t) $ 。然后将$ {\overset{\frown }{{\chi }}}(t) $ 代入微分方程组进行求解积分，即可基于动力学机制重建得到预测极移$ {\overset{\frown }{{p}}}(t) $ 结果。
为了确保数值解的高精度与稳定性，使用RK4方法，结合初始极移条件与预测激发函数对极移动力学微分方程组进行积分运算。取时间步长h=1天，设$ {t}_{n} $ 时刻的极移状态为$ {\overset{\frown }{p}}_{n} $ ，激发函数为$ {\overset{\frown }{\chi }}_{n} $ ，则第tn+1时刻极移状态$ {\overset{\frown }{p}}_{n+1} $ 的更新如下：
式中，$ f({t}_{n},{\overset{\frown }{p}}_{n},{\overset{\frown }{\chi }}_{n}) $ 为极移动力学微分方程组的右侧函数；激发函数的中间时刻值$ {\overset{\frown }{\chi }}_{n+1/2} $ 通过线性插值获取，即$ {\overset{\frown }{\chi }}_{n+1/2}=（{\overset{\frown }{\chi }}_{n+1}+{\overset{\frown }{\chi }}_{n}）/2 $ 。积分过程的初始条件p0设定为观测数据的最后一个历元值。通过迭代执行上述RK4步骤，可由预测激发函数序列$ {\overset{\frown }{\chi }}_{n} $ 逐步递推出未来极移轨迹$ {\overset{\frown }{p}}_{n} $ 。该方法遵循极移的物理运动规律，有效保障了预测结果的动力学一致性。
2 极移预报方法比较
2.1 评估方案
本文通过以下方案，全面评估本文极移预报方法的有效性。
2.1.1 数据集构建
采用2019—2024年的IERS C04高精度极移时间序列(数据时间间隔为1天)，用于构建大地测量激发函数外推计算所需的历史窗口和未来预报的真实观测结果。以预报时间点之前的6年历史数据作为激发函数WLS模型的输入，并用于确定积分初始条件。
2.1.2 基准选择
2.1.3 预测任务
预测时段选取为2019年至2023年。鉴于Bulletin A每7天发布一期跨度为365天的极移预报，该时段内共包含261个独立预报历元。为确保对比的公平性与数据一致性，本文设定的起报时刻与Bulletin A的发布时间严格对齐。基于前文所述的RK4积分方法，利用WLS模型外推得到的激发函数序列$ {\overset{\frown }{{\chi }}}(t) $ ，递推解算各历元未来365天的极移预测值$ {\overset{\frown }{{p}}}(t) $ ；同时，运行上述构建的LS+AR模型独立生成相同时间窗口的预测序列，并同步提取同期的Bulletin A预报数据，对三者的预测性能进行对比分析。
2.1.4 评价体系
将预测极移序列$ {\overset{\frown }{{p}}}(t) $ 与真实的IERS EOP 20 C04数据集中极移的测量真实值进行比对，本文使用多维度指标来全面量化评估算法的预测性能，使用的指标如下：
(1)绝对误差(Absolute Error，AE)。衡量单次预测值与真实值之间的偏差大小，直接反映个别预测的准确程度：
(2)平均绝对误差(Mean Absolute Error，MAE)。用于衡量所有预测误差的平均绝对水平，因其对异常值不敏感、更具鲁棒性，能够更稳健地评估模型整体预测性能：
(3)均方根(Root Mean Square，RMS)误差。衡量预测误差的偏差程度，能更敏锐地捕捉大偏差数据，从而确保模型在捕捉极端信号波动时的可靠性：
式中，Pi为第i时刻的极移预测值，即本文方法或Bulletin A预测值；Oi为第i时刻的IERS C04极移观测值。
(4)误差降低率(Error Rate Reduction，ERR)。用于定量评估本文方法相对其他方法在预报结果上的误差降低程度：
式中，Ebase为基准模型误差，本文中为Bulletin A的RMS误差；Enew为本文提出新方法的RMS误差；Rerr为正值代表误差降低，负值代表误差提高。
2.2 结果分析
为剔除极端异常值的影响并评估模型的稳健性，进一步计算了全预报时段的综合MAE。统计结果显示，X分量在365天、90天及10天预报跨度下的MAE分别为13.16 mas、6.51 mas和1.83 mas；而Y分量在同期的MAE分别为11.89 mas、4.36 mas和1.06 mas。数据表明，无论在长期还是短中期尺度上，Y分量的预测精度均优于X分量。
这种误差分布的不对称性源于地球自转的动力学耦合机制。根据式(8)，极移X分量的运动主要由激发函数Y分量χy驱动。图2 的频谱分析表明，χy的周年振幅远大于χx，呈现出极强的能量主导性，使得半年和季度等次要周期成分相对微弱。同时，χy巨大的振幅基数导致了建模误差的放大效应——即便相位或幅度存在微小的相对偏差，经过动力学积分后也会转化为极移X分量显著的轨迹偏离。这使得极移X分量的预测难度在物理本质上高于Y分量，最终导致其统计误差偏大。
表 1 本文方法对比Bulletin A和LS+AR方法的ERRTable 1 The ERR of the proposed method compared with Bulletin A and the LS+AR method |
| 预报长度/天 | 对比Bulletin A | 对比LS+AR | |||
| X分量/% | Y分量/% | X分量/% | Y分量/% | ||
| 1 10 90 365 | −1.19 −5.00 13.27 21.08 | 7.12 3.42 23.17 16.26 | 27.01 23.92 52.15 43.61 | 41.21 42.16 63.42 46.50 | |
而在中长期预报区间(>10天)本文方法展现出明显优势，预报性能优于Bulletin A：预报长度为90天时，本文方法对比Bulletin A的X、Y分量分别达到13.27%、23.17%；预报长度为365天时，对应X、Y分量分别为21.08%、16.26%，可见本文方法的预报性能效果更优。
相比之下，LS+AR方法表现最差，其全时间段内的预报精度都弱于本文方法与Bulletin A。从表1 可知，预报长度为1天时，本文方法对比LS+AR的X、Y分量为27.01%、41.21%；预报长度10天时为23.92%、42.16%；预报长度为90天时达52.15%、63.42%；预报长度365天时为43.61%、46.50%，各时长下LS+AR的误差均显著高于本文方法与Bulletin A。
3 极移预报算法对定轨影响研究
然而，现有研究表明这两类影响的量级存在显著差异。极移误差在轨道积分过程中产生的累积影响通常在亚毫米量级，相比之下，GCRS到ITRS的坐标系转换过程是极移误差影响定轨精度的主导因素[28]。鉴于此，本文将重点聚焦于坐标系转换过程，对比研究不同极移预报方法对航天器定轨产生的几何误差影响。
3.1 极移在坐标系转换中误差的传播机理
在航天器精密定轨任务中，实现GCRS与ITRS之间的高精度转换至关重要。根据IERS规范，惯性系下的位置矢量$ {\boldsymbol{r}}_{\text{GCRS}} $ 与地固系下的位置矢量$ {\boldsymbol{r}}_{\text{ITRS}} $ 之间的转换关系可表示为[1]：
式中，R(t)、Q(t)分别是地球自转和岁差章动引起的旋转矩阵；W(t)则是极移px和py决定的旋转矩阵。由于极移是微小量，其旋转矩阵可以近似展开为：
当进行极移预报时，会引入极移误差Δpx和Δpy。该误差会通过旋转矩阵传递给航天器，从而造成航天器的定轨偏差。对于在GCRS下位置为$ {\boldsymbol{r}}_{\text{GCRS}} $ 的航天器，则由极移误差影响造成的偏差Δr可表示为[28]：
3.2 数值仿真设置
为了量化极移预报方法对航天器轨道确定的影响，并验证本文预报方法的有效性，设计了涵盖多种卫星轨道的对比仿真实验。
(1)仿真对象。为了全面评估极移预报误差对航天器定轨的影响，本文从全球导航卫星系统(GPS、GLONASS、Galileo及BDS)中选取6颗具有代表性的卫星的参数作为仿真对象。如表2 所示，这些卫星覆盖了中地球轨道(Medium Earth Orbit，MEO)、地球静止轨道(Geostationary Earth Orbit，GEO)及倾斜地球同步轨道(Inclined Geosynchronous Orbit，IGSO)多种轨道类型，且轨道形状均设置为圆轨道。
(2)仿真方案。时间范围设定在2023年全年。以Bulletin A的发布时间点作为数值仿真初始时间，共设置52组实验。每组从发布日起向后进行90天的极移预报。以IERS EOP 20 C04数据集的极移测量结果作为基准，将Bulletin A预报结果作为对照，与本文方法的预报结果进行对比分析。
(3)定轨误差解算流程。对于每一组仿真，首先基于开普勒动力学方程生成各卫星GCRS下90天的标准轨道位置$ {\boldsymbol{r}}_{\text{GCRS}}(t) $ 和速度$ {\boldsymbol{v}}_{\text{GCRS}}(t) $ ; 随后分别计算Bulletin A与本文预报方法相对于C04的测量真值的极移偏差Δpx、Δpy；利用该偏差构造误差旋转矩阵，解算地固系下的位置偏差矢量Δr；最后将该位置偏差矢量投影至卫星的径向-切向-法向坐标系，并对全年52组仿真结果进行平均处理，从而消除季节性因素带来的偶然影响。
3.3 结果分析
从三轴方向分解来看，所有卫星的径向误差均处于10−16 m量级。误差主要集中在切向与法向，其幅值随极移预报时长的延长呈同步增长趋势。切向误差曲线表现出较好的平滑特征，而法向误差则表现出与轨道特性相关的周期性振荡特征。从卫星轨道差异方面来看，误差表现出明显的放大效应，即随着轨道高度的增加，误差显著增大。如图5 中BDS IGSO卫星的误差曲线要始终高于其他卫星，90天最大误差约为1.5 m。值得注意的是，BDS GEO卫星轨道倾角i=0位于赤道平面，导致极移误差引起的定轨误差几乎全部投影至了法向，呈现出法向误差明显较高且振荡、但径向与切向误差极小的特点。图6 进一步量化了综合90天预报期内的坐标系转换的均方根误差，以及本文方法相较于Bulletin A预报结果的误差降低率。使用Bulletin A预报结果时，MEO卫星的RMS误差为1.1~1.3 m，高轨卫星误差增长至1.5~1.9 m，同样可以看出轨道高度对误差的放大效应。使用本文预报方法后，所有卫星的综合90天RMS误差能够控制在1.0 m以内。其中，对MEO卫星的ERR稳定在23.97%~30.66%；对“北斗”系统卫星的提升幅度为27.93%~35.93%。BDS GEO卫星的改善最为显著，由于其轨道倾角的特殊性，对极移误差十分敏感，其预报效果提升达35.93%。
4 结 语
本文利用刘维尔方程将极移反演到物理激发域，通过引入时间衰减因子的WLS方法对激发函数进行外推，最终通过动力学积分实现极移序列的预报。基于2019—2023年数据的261组极移预报比较表明，本文方法对中长期的极移预报有明显的精度提升。在短期预报中，该方法与现有高精度产品精度相当；而在中长期预报区间内，相较于经典的LS+AR模型，X、Y分量在全时段的平均预报精度提升了约43.61%~46.50%，相较于Bulletin A，精度提升幅度亦达到16.26%~21.08%。
在进一步的定轨数值仿真中，通过对GPS、“北斗”等导航卫星的仿真数据进行参考坐标系转换，验证本文方法在定轨中的影响，并与Bulletin A的结果进行比对。实验结果表明，极移带来的框架转换误差有显著的距离放大效应：卫星轨道的径向误差最小，仅为10−16 m量级，基本可以忽略；切向和法向误差均随预报时长增加而增大，其中法向误差表现出明显的振荡特征。对比Bulletin A的定轨结果，本文方法带来了23.97%~35.93%的精度提升。综上，本文方法在极移的中长期预报中带来了明显的精度提升，可有效降低定轨中的参考坐标系转换误差，为高精度自主定轨及深空探测任务提供可靠支撑。
尽管本文方法在中长周期预报中取得了良好效果，但仍存在进一步优化的空间。未来研究可尝试引入大气角动量、海洋角动量等高精度地球物理流体预报产品，进一步捕捉极移的高频变化特征，提高短期预报精度；在技术方法方面，可结合人工智能技术与地球动力学机制，构建融合多源地球物理信息与物理机制的混合预报模型，探索地球自转变化中潜在的深层机制。