空间科学与试验学报

空间科学与试验学报 >

2026 , Vol. 3 >Issue 1: 26 - 35

DOI: https://doi.org/10.19963/j.cnki.2097-4302.2026.01.004

空间科学专栏

基于地球自转动力学与加权最小二乘法的极移预报及其在定轨中的应用

  • 黄俊海 , 1 ,
  • 朱柏洋 1 ,
  • 刘晓东 , 1, 2, * ,
  • 吴志刚 1 ,
  • 王巍 3
展开
  • 1. 中山大学航空航天学院, 深圳 518107
  • 2. 中山大学深圳市智能微小卫星星座技术与应用重点实验室, 深圳 518107
  • 3. 中国航天科技集团有限公司, 北京 100048
（1985—），男，教授。主要研究方向为深空探测、轨道动力学、地球自转、太阳系尘埃/空间碎片/星座动力学、行星形成。（本文通信作者）通信地址：广东省深圳市光明区公常路66号中山大学深圳校区（518107）电子邮箱：

（2001—），男，硕士研究生。主要研究方向为机器学习、地球自转、极移预测。通信地址：广东省深圳市光明区公常路66号中山大学深圳校区（518107）电子邮箱：

网络出版日期: 2026-04-10

基金资助

国家自然科学基金(62388101)；空间碎片与近地小行星防御科研专项(KJSP2023020202)资助

收起

Dynamics-Driven Polar Motion Prediction Using Weighted Least Squares: Application to Precision Orbit Determination

  • Junhai HUANG , 1 ,
  • Boyang ZHU 1 ,
  • Xiaodong LIU , 1, 2, * ,
  • Zhigang WU 1 ,
  • Wei WANG 3
Expand
  • 1. School of Aeronautics and Astronautics, Sun Yat-sen University, Shenzhen 518107, China
  • 2. Shenzhen Key Laboratory of Intelligent Microsatellite Constellation, Sun Yat-sen University, Shenzhen 518107, China
  • 3. China Aerospace Science and Technology Corporation, Beijing 100048, China

Online published: 2026-04-10

Fold

摘要

航天器精密定轨依赖于高精度的参考坐标系转换，而极移是实现该转换的必要参数之一。针对实测数据存在约13天时滞、传统预报缺乏物理机制的局限，提出一种基于地球自转动力学与加权最小二乘法的融合预报方法。该方法首先利用刘维尔方程将极移反演为激发函数，随后引入时间衰减因子进行加权最小二乘法外推，最终通过动力学积分重构预报序列。基于2019—2023年数据的验证实验表明，在1~365天内，该方法的精度优于最小二乘+自回归方法及国际地球自转和参考系服务发布的Bulletin A预报产品：较最小二乘+自回归方法提升精度约43.61%～46.50%，较Bulletin A提升精度约16.26%～21.08%。通过90天定轨仿真显示，采用该方法后，全球定位系统和“北斗”卫星导航系统等的定轨精度较Bulletin A提升23.97%～35.93%。结果表明，引入动力学约束能有效克服中长周期预报发散，减少定轨中的坐标系转换误差，为高精度自主定轨及深空探测提供可靠的技术支撑。

关键词： 地球自转预报; 极移预报; 加权最小二乘法; 卫星定轨

本文引用格式

黄俊海 , 朱柏洋 , 刘晓东 , 吴志刚 , 王巍 . 基于地球自转动力学与加权最小二乘法的极移预报及其在定轨中的应用[J]. 空间科学与试验学报, 2026 , 3(1) : 26 -35 . DOI: 10.19963/j.cnki.2097-4302.2026.01.004

Abstract

High-accuracy reference frame transformations for spacecraft precision orbit determination rely on polar motion. To address the 13-day data latency in observed data and the lack of physical mechanisms in traditional prediction models, we propose a fusion prediction method combining Earth rotation dynamics and weighted least squares. This method inverts polar motion series into excitation functions via the Liouville equation, extrapolates them using a weighted least squares model with a time decay factor, and reconstructs predictions via dynamic integration. Experiments using 2019—2023 data demonstrate that this approach outperforms least squares+autoregressive model and International Earth Rotation and Reference Systems Service Bulletin A over 1~365 days, improving accuracy by 43.61%~46.50% and 16.26%~21.08%, respectively. Furthermore, 90-day simulations for navigation satellites such as Global Positioning System and Beidou Navigation Satellite System reveal a 23.97%~35.93% accuracy enhancement over Bulletin A. Results indicate that dynamic constraints effectively mitigate medium-to-long-term divergence, reducing frame transformation errors and supporting high-precision autonomous precision orbit determination and deep space exploration.

Key words： Earth rotation prediction; polar motion prediction; weighted least squares; satellite orbit determination

0 引 言

地球极移是指地球瞬时自转轴在地球本体坐标系内的运动，表现为地极点在地球表面位置的周期性变化。作为地球方向参数(Earth Orientation Parameters， EOP)的核心分量，极移是连接国际天球参考系与国际地球参考系的关键转换参数[1-2]。
随着航天强国战略的深入推进，我国相继部署了“北斗”全球卫星导航系统(Beidou Navigation Satellite System，BDS)[3]、“天问”二号探测器[4]等一系列重大任务。为确保这些任务的顺利实施，必须实现航天器的高精度定轨，而这一过程严格依赖于国际天球参考系与国际地球参考系之间的精确转换。因此，这也对极移参数的高精度测量与预报提出了更为严苛的要求[5]。
目前，空间大地测量技术已发展出多种极移观测手段，包括月球激光测距、甚长基线干涉测量、全球定位系统(Global Positioning System，GPS)和多普勒轨道与卫星无线电定位集成系统[6]。近年来，随着大型高精度光学陀螺仪的精度大幅提升，激光陀螺仪、光纤干涉仪也可以实现地球自转和极移的高精度测量[7]。然而，受限于传统空间大地测量技术观测数据的物理传输延迟以及复杂参数解算流程的计算耗时，高精度极移观测产品普遍存在显著的时间滞后性。根据国际全球导航卫星系统服务组织发布的EOP最终产品通常存在约13天的延迟[8]。为了弥补这个缺口，需要发展高精度的极移预测技术。
针对极移预报问题，国内外学者开展了广泛深入的研究，提出了多种预报模型与方法。早期研究中，Chao[9]率先利用最小二乘法(Least Squares，LS)提取极移序列中的长期趋势与周期项，构建了外推模型进行极移的预报。随后，Kosek等[10]在LS模型基础上引入自回归(AutoRegressive，AR)技术，对LS拟合后的残差序列进行随机建模。凭借结构简单且在短期预报中表现稳健的特点，LS+AR组合模型逐渐成为极移预报领域的经典范式[11]。此后，许多学者对该范式进行了优化与改进，如张昊等[12]通过引入加权策略构建加权最小二乘(Weighted Least Squares，WLS)模型，有效增强了中长期预报的稳定性。随着人工智能技术的快速发展，其处理复杂非线性映射的强大能力为极移预报注入了新活力。Schuh等[13]最早探索了人工神经网络在极移预报中的应用。近年来，长短期记忆网络等深度学习模型也被使用在极移的预报任务中，进一步提升了极移预报的精度[14-15]。
然而，上述纯数据驱动方法通常侧重于挖掘极移序列本身的统计规律，在一定程度上未能充分考量地球自转背后的物理激发机制。为了弥补这一缺陷，学界将目光转向极移的动力学激发的研究[16-18]。Dill等[19]在极移的物理激发域的基础上，融合了大气与海洋的有效角动量数据，提升了极移的短期预报精度。在第二届地球定向参数预报比较活动中[11]，Dill等[20]凭借融合大气、海洋及陆地水圈物理信息，结合最小二乘外推方法进行建模，在极移的预报项目中表现卓越。这表明，将地球自转动力学机制与数学模型深度融合，是突破当前极移预报精度瓶颈的重要途径。
本文旨在研究融合地球自转动力学机制与WLS的极移预报方法。该方法首先依据地球自转动力学方程，将极移观测序列转换为对应的激发函数；继而使用WLS提取激发函数中的主要周期项，并引入时间衰减因子对历史建模数据实施动态加权，增强近期观测数据对模型的贡献度；最后，通过动力学积分将预测的激发函数转换回极移数据，完成极移序列的外推预报。该方法通过结合动力学约束与数据驱动，在确保短期预报精度的同时，提高了中长期预报的稳定性。此外，本文还进一步分析本文预报方法在航天器定轨中的应用效果，量化不同预报时长对轨道转换精度的影响，为我国航天任务中的高精度参考坐标系转换提供理论依据与技术支持。

1 极移预测的动力学机制与建模

1.1 极移的激发函数

地球的自转运动不是简单的刚体旋转，而是受到地球外部力矩和地球内部多圈层影响的复杂运动。考虑到地球的非刚体性和各圈层导致的质量与角动量变化的方程，被称为欧拉-刘维尔方程。通过对欧拉-刘维尔方程进一步简化，地球的动力学方程可简化为一组更易于分析的线性微分方程，即线性刘维尔方程。其中，极移的动力学方程可以表示为[6,21]：
$ \dfrac{i}{\sigma_{\mathrm{cw}}}\dfrac{\mathrm{d}{p}(t)}{\mathrm{d}t}+{p}(t)={\chi}(t) $
式中，$ {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]：
$ {{\chi}}_{\mathrm{geo}}(t)=\dfrac{{{p}}(t)-{{p}}(t-\Delta t)\mathrm{e}^{i{{\sigma}}_{\mathrm{cw}}\Delta t}}{1-\mathrm{e}^{i{{\sigma}}_{\mathrm{cw}}\Delta t}} $
式中，$ \Delta t $是离散极移数据的时间步长。本文选用IERS发布的EOP 20 C04数据集，该数据集代表了当前极移观测数据的最高精度水平。如图1所示，基于该数据集提供的日分辨率极移序列，利用式(2)反演计算出大地测量激发函数。
图 1 大地测量激发

Fig.1 The geodetic excitation

1.2 激发函数频率分析

利用快速傅里叶变换方法对大地测量激发函数$ {{\chi }}_{\text{geo}}(t) $进行频谱分析。1993年以来，随着全球定位系统等现代空间大地测量技术的引入，极移观测精度得到了显著提升，从而有效捕捉到了更为真实和丰富的高频信号[23]。本文选取1993—2024年观测数据进行分析，以避免早期数据中潜在的高频插值误差干扰。如图2的频率谱图所示，在周期小于400天的范围内，极移的激发函数呈现显著的周期性特征，其能量主要周期集中在365天、180天及120天等季节性周期附近。这种周期性特征的物理根源，可追溯至极移的主要激发源——大气质量和海洋质量的周期性迁移，以及大气表层与海底压力变化引发的地球形变。这些动力学过程通过角动量守恒机制，将自身周期性变化传递至地球自转，最终在极移的大地测量激发函数中显示出可观测的周期性信号[24]。
图 2 大地测量激发的频率谱

Fig.2 Frequency spectrum diagram of the geodetic excitation

同时，在周期小于30天的范围内，极移还拥有丰富的高频信号，主要是由潮汐效应和大气、洋流在短期内导致的气压和海底压力变化驱动[25-27]。由于极移的低频分量特征明显且其能量稳定，而高频分量有着极强的波动性，难以进行准确的建模外推。本文选择仅对主要的低频特征进行建模，对于高频分量使用低通滤波器进行滤波处理，以防止高频能量泄漏到低频信号中。

1.3 激发函数的建模与外推

在获取历史时段的激发函数$ {\boldsymbol{\chi }}_{\text{geo}}(t) $后，根据第1.2节所述的频谱分析结果，确定了其主要周期成分(T1=365.250天，T2=182.625天，T3=121.750天)。为了从激发函数的原始序列中重建出这些主要成分，本文构建了包含线性趋势和周期项的谐波模型$ S(t) $：
$ S(t)={a}_{0}+{a}_{1}t+\sum\limits_{k=1}^{n}[{C}_{k}\cos ({\omega }_{k}t)+{D}_{k}\sin ({\omega }_{k}t)] $
式中，t为观测历史窗口时间，$ {\omega }_{k}=2{\text π} /{T}_{k} $为主要周期的角频率，a0为常数项，a1、Ck、Dk分别为线性趋势项和周期项待定系数。为了求解上述待定系数，将离散的观测方程表示为矩阵形式：
$ \boldsymbol{Y}=\boldsymbol{A}\boldsymbol{X}+\boldsymbol{V} $
式中，Y为观测向量，A是由常数项、线性趋势项和周期项构成的设计矩阵，X为待定系数向量，V是残差向量。使用最小二乘法对原始序列进行拟合，通过最小化残差平方和，可得待定系数的最优解为：
$ {{\overset{\frown }{\boldsymbol{X}}}}_{\text{OLS}}={({{\boldsymbol{A}}^{\text{T}}}\boldsymbol{A})}^{-1}{\boldsymbol{A}}^{\text{T}}\boldsymbol{Y} $
由于极移激发函数受大气、海洋等多源驱动机制的影响，具有明显的时变性，所以近期观测数据对模型构建具有更强的信息价值。为此，本文构建了基于时间的权重函数：通过动态调整历史观测数据的权重系数，设计权重随时间衰减，赋予近期数据高权重，从而将模型的注意力转移到近期的数据中。时间权重wi为：
$ {w}_{i}=\max ({\lambda }^{({{t}_{n}}-{{t}_{i}})},{w}_{\min }） $
式中，i=1, 2, ··· , n为历史观测数据窗口中的采样编号，ti为对应的时间点。设置遗忘因子λ=0.999，其对应半衰期时长约为两年；同时设定权重下限wmin=0.5，防止历史数据被过度衰减导致序列有效长度缩短。将权重组成权重矩阵W=diag{w1, w2, ···,wn}，则WLS的参数估计为：
$ {{\overset{\frown }{\boldsymbol{X}}}}_{\text{WLS}}={({{\boldsymbol{A}}^{\text{T}}}\boldsymbol{W}\boldsymbol{A})}^{-1}{\boldsymbol{A}}^{\text{T}}\boldsymbol{W}\boldsymbol{Y} $

1.4 极移积分预测

为实现极移的动力学预测，本文构建了基于刘维尔方程的积分框架。首先推导了极移刘维尔方程的积分形式表达式，然后采用四阶龙格-库塔(Runge-Kutta，RK4)数值方法求解该动力学方程，从而建立激发函数与极移轨迹的映射关系。基于极移刘维尔方程的解析形式，将极移位移量、激发函数及钱德勒频率的复数域表达式代入动力学模型，可推导得到极移动力学微分方程组：
$ \begin{cases} {\dot{p}}_{x}=\dfrac{2{\text π} }{{T}_{\text{cw}}}\bigg(({p}_{y}+{\chi }_{y})-\dfrac{1}{2Q}({p}_{x}-{\chi }_{x})\bigg)\\{\dot{p}}_{y}=-\dfrac{2{\text π} }{{T}_{\text{cw}}}\bigg(({p}_{x}-{\chi }_{x})+\dfrac{1}{2Q}({p}_{y}+{\chi }_{y})\bigg)\end{cases} $
式中，$ {\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} $的更新如下：
$ \begin{cases} {k}_{1}=f({t}_{n},{\overset{\frown }{p}}_{n},{\overset{\frown }{\chi }}_{n})\\{k}_{2}=f({t}_{n}+\dfrac{h}{2},{\overset{\frown }{p}}_{n}+\dfrac{h}{2}{k}_{1},{\overset{\frown }{\chi }}_{n+1/2})\\{k}_{3}=f({t}_{n}+\dfrac{h}{2},{\overset{\frown }{p}}_{n}+\dfrac{h}{2}{k}_{2},{\overset{\frown }{\chi }}_{n+1/2})\\{k}_{4}=f({t}_{n}+h,{\overset{\frown }{p}}_{n}+h{k}_{3},{\overset{\frown }{\chi }}_{n+1})\\{\overset{\frown }{p}}_{n+1}={\overset{\frown }{p}}_{n}+\dfrac{h}{6}({k}_{1}+2{k}_{2}+2{k}_{3}+{k}_{4})\end{cases} $
式中，$ 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 基准选择

选取IERS Bulletin A发布的官方预报产品，与经典的LS+AR方法共同作为对比基准。其中，Bulletin A代表当前极移预报业务化应用的主流水平，LS+AR方法则是极移预测领域长期认可的传统统计学基准算法[11]。该方法通过6年长度的历史极移序列对钱德勒摆动(周期433.000天)、周年摆动(周期365.250天)及半周年摆动(周期182.625天)进行参数建模；随后，对最小二乘拟合后的残差序列使用AR建模，并通过赤池信息量准则动态优选最优模型阶数，从而实现对随机残差分量的精确表征与预测[10]。

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)。衡量单次预测值与真实值之间的偏差大小，直接反映个别预测的准确程度：
$ {E}_{\text{AE}}=|{P}_{i}-{O}_{i}| $
(2)平均绝对误差(Mean Absolute Error，MAE)。用于衡量所有预测误差的平均绝对水平，因其对异常值不敏感、更具鲁棒性，能够更稳健地评估模型整体预测性能：
$ {E}_{\text{MAE}}=\sum\limits_{i=1}^{n}|{P}_{i}-{O}_{i}| $
(3)均方根(Root Mean Square，RMS)误差。衡量预测误差的偏差程度，能更敏锐地捕捉大偏差数据，从而确保模型在捕捉极端信号波动时的可靠性：
$ {E}_{\text{RMS}}=\sqrt{\dfrac{1}{n}\sum\limits_{i=1}^{n}{({{P}_{i}}-{{O}_{i}})}^{2}} $
式中，Pi为第i时刻的极移预测值，即本文方法或Bulletin A预测值；Oi为第i时刻的IERS C04极移观测值。
(4)误差降低率(Error Rate Reduction，ERR)。用于定量评估本文方法相对其他方法在预报结果上的误差降低程度：
$ {R}_{\text{err}}=\dfrac{{E}_{\text{base}}-{E}_{\text{new}}}{{E}_{\text{base}}}\times 100\mathrm{\% } $
式中，Ebase为基准模型误差，本文中为Bulletin A的RMS误差；Enew为本文提出新方法的RMS误差；Rerr为正值代表误差降低，负值代表误差提高。

2.2 结果分析

图3展示了2019—2023年期间，基于本文方法、在7天滑动预报间隔下，预报跨度为365天的极移预测误差统计结果。由图3可以看出，在365天长周期预报中，X分量的最大误差为66.00 mas，略高于Y分量的63.00 mas；在90天中周期预报中，X与Y分量的最大误差分别为32.00 mas和27.00 mas；在10天短周期预报中，两者分别为11.00 mas和6.00 mas。随着预报长度的增加，由于积分过程的误差积累效应，预测的误差明显增大。
图 3 本文算法的绝对误差

Fig.3 Absolute error of the proposed method

为剔除极端异常值的影响并评估模型的稳健性，进一步计算了全预报时段的综合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分量，最终导致其统计误差偏大。
图4呈现了2019—2023年期间，本文方法与Bulletin A及传统LS+AR方法的综合RMS误差对比。三种方法的预报误差均随预测跨度延长呈现非线性递增趋势：其中，在短期预报阶段(0~10天)，本文方法与Bulletin A的误差水平保持相近，二者误差幅度差异较小。从表1的ERR数据可以看出，预报长度为1天时，本文方法对比Bulletin A的X分量为−1.19%、Y分量为7.12%；预报长度为10天时，对应X、Y分量分别为−5.00%、3.42%，二者误差幅度差异较小。
图 4 不同方法的RMS误差对比

Fig.4 Comparison of RMS error of different methods

表 1 本文方法对比Bulletin A和LS+AR方法的ERR

Table 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 极移预报算法对定轨影响研究

极移对航天器精密定轨的影响机制主要包含两个方面[28]：其一，在坐标系转换过程中，极移预报误差直接导致地心天球参考系(Geocentric Celestial Reference System，GCRS)与ITRS之间转换的旋转矩阵偏差，从而产生几何位置误差；其二，在航天器轨道积分过程中，极移误差在动力学模型引入了附加摄动项。这种极移摄动项源于地球非球形重力场的定义方式——重力场系数通常定义在随地球旋转的ITRS中，而轨道积分则在惯性系GCRS中进行，两者之间的转换使得极移误差与动力学方程耦合，形成时变的极移摄动加速度[29]。
然而，现有研究表明这两类影响的量级存在显著差异。极移误差在轨道积分过程中产生的累积影响通常在亚毫米量级，相比之下，GCRS到ITRS的坐标系转换过程是极移误差影响定轨精度的主导因素[28]。鉴于此，本文将重点聚焦于坐标系转换过程，对比研究不同极移预报方法对航天器定轨产生的几何误差影响。

3.1 极移在坐标系转换中误差的传播机理

在航天器精密定轨任务中，实现GCRS与ITRS之间的高精度转换至关重要。根据IERS规范，惯性系下的位置矢量$ {\boldsymbol{r}}_{\text{GCRS}} $与地固系下的位置矢量$ {\boldsymbol{r}}_{\text{ITRS}} $之间的转换关系可表示为[1]：
$ {\boldsymbol{r}}_{\text{ITRS}}=\boldsymbol{W}(t)\boldsymbol{R}(t)\boldsymbol{Q}(t){\boldsymbol{r}}_{\text{GCRS}} $
式中，R(t)、Q(t)分别是地球自转和岁差章动引起的旋转矩阵；W(t)则是极移px和py决定的旋转矩阵。由于极移是微小量，其旋转矩阵可以近似展开为：
$ \boldsymbol{W}(t)={\boldsymbol{R}}_{1}({p}_{x}){\boldsymbol{R}}_{2}({p}_{y})\approx \left[\begin{matrix}1 & 0 & -{p}_{x}\\0 & 1 & -{p}_{y}\\{p}_{x} & {p}_{y} & 1\end{matrix}\right] $
当进行极移预报时，会引入极移误差Δpx和Δpy。该误差会通过旋转矩阵传递给航天器，从而造成航天器的定轨偏差。对于在GCRS下位置为$ {\boldsymbol{r}}_{\text{GCRS}} $的航天器，则由极移误差影响造成的偏差Δr可表示为[28]：
$ \Delta \boldsymbol{r}=\left[\begin{matrix}0 & 0 & \Delta {p}_{x}\\0 & 0 & -\Delta {p}_{y}\\-\Delta {p}_{x} & \Delta {p}_{y} & 0\end{matrix}\right]{\boldsymbol{r}}_{\text{GCRS}} $

3.2 数值仿真设置

为了量化极移预报方法对航天器轨道确定的影响，并验证本文预报方法的有效性，设计了涵盖多种卫星轨道的对比仿真实验。
(1)仿真对象。为了全面评估极移预报误差对航天器定轨的影响，本文从全球导航卫星系统(GPS、GLONASS、Galileo及BDS)中选取6颗具有代表性的卫星的参数作为仿真对象。如表2所示，这些卫星覆盖了中地球轨道(Medium Earth Orbit，MEO)、地球静止轨道(Geostationary Earth Orbit，GEO)及倾斜地球同步轨道(Inclined Geosynchronous Orbit，IGSO)多种轨道类型，且轨道形状均设置为圆轨道。
表 2 仿真卫星参数[30-33]

Table 2 Parameters of the simulation satellite[30-33]

卫星星座 轨道类型 轨道高度/km 倾角/(°)
GPS
GLONASS
Galileo
BDS
BDS
BDS		 MEO
MEO
MEO
MEO
IGSO
GEO		 20 189
19 100
23 189
21 528
35 786
35 786		 55.0
64.8
56.0
55.0
55.0
0
(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 结果分析

图5展示了2023年52组仿真的平均结果，反映了本文预报方法在不同轨道卫星的坐标系转换中，由此引入的在径向、切向和法向上的误差随时间演化特性。结果表明，由极移预报引起的定轨误差具有明显的方向性与轨道依赖性。
图 5 本文方法在不同轨道卫星上的坐标系转换误差演化

Fig.5 Temporal evolution of frame transformation errors for different satellites using the proposed method

从三轴方向分解来看，所有卫星的径向误差均处于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%。
图 6 90天预报期Bulletin A与本文方法的坐标系转换均方根误差对比及误差降低率

Fig.6 Comparison of 90-day frame transformation RMSE and ERR between Bulletin A and the proposed method

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%的精度提升。综上，本文方法在极移的中长期预报中带来了明显的精度提升，可有效降低定轨中的参考坐标系转换误差，为高精度自主定轨及深空探测任务提供可靠支撑。
尽管本文方法在中长周期预报中取得了良好效果，但仍存在进一步优化的空间。未来研究可尝试引入大气角动量、海洋角动量等高精度地球物理流体预报产品，进一步捕捉极移的高频变化特征，提高短期预报精度；在技术方法方面，可结合人工智能技术与地球动力学机制，构建融合多源地球物理信息与物理机制的混合预报模型，探索地球自转变化中潜在的深层机制。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
1
PETIT G, LUZUM B. IERS conventions (2010)[R]. Frankfurt: Verlag des Bundesamts für Kartographie und Geodäsie, 2010.

2
王巍, 冯文帅, 张首刚. 基于光纤干涉仪的世界时高精度测量技术[M]. 北京: 科学出版社, 2024.

3
项宇, 蒋孝卿, 杨建华, 等. 地球定向参数预报误差及其对北斗三号卫星定轨精度的影响[J]. 天文学进展, 2023, 41 (2): 269- 278.

XIANG Y, JIANG X Q, YANG J H, et al. Influence of Earth orientation parameter forecast error on the orbit determination accuracy of Beidou-3 satellite[J]. Progress in Astronomy, 2023, 41 (2): 269- 278.

4
王波, 鄢建国, 高梧桐, 等. EOP预报误差对深空探测器精密定轨结果影响分析[J]. 武汉大学学报（信息科学版）, 2024, 49 (9): 1538- 1545.

DOI

WANG B, YAN J G, GAO W T, et al. Impact analysis of EOP prediction errors on orbit determination of deep-space spacecraft[J]. Geomatics and Information Science of Wuhan University, 2024, 49 (9): 1538- 1545.

DOI

5
张卫星, 刘万科, 龚晓颖. EOP预报误差对自主定轨结果影响分析[J]. 大地测量与地球动力学, 2011, 31 (5): 106- 110.

ZHANG W X, LIU W K, GONG X Y. Influence of EOP prediction error on autonomous orbit determination of navigation satellite[J]. Journal of Geodesy and Geodynamics, 2011, 31 (5): 106- 110.

6
GROSS R S. Earth rotation variations-long period[J]. Treatise on Geophysics, 2007, 3, 239- 294.

DOI

7
王巍, 冯文帅, 张首刚, 等. 基于高精度光纤干涉仪的世界时测量技术研究[J]. 导航与控制, 2023, 22 (5): 1- 11.

DOI

WANG W, FENG W S, ZHANG S G, et al. Research on universal time measurement technology based on high-precision fiber optic interferometer[J]. Navigation and Control, 2023, 22 (5): 1- 11.

DOI

8
DOW J M, NEILAN R E, RIZOS C. The international GNSS service in a changing landscape of Global Navigation Satellite Systems[J]. Journal of Geodesy, 2009, 83 (3): 191- 198.

9
CHAO B F. Predictability of the Earth’s polar motion[J]. Journal of Geodesy, 1985, 59 (1): 81- 93.

DOI

10
KOSEK W, MCCARTHY D D, LUZUM B J. Possible improvement of Earth orientation forecast using autocovariance prediction procedures[J]. Journal of Geodesy, 1998, 72 (4): 189- 199.

DOI

11
ŚLIWIŃSKA J, KUR T, WIŃSKA M, et al. Second earth orientation parameters prediction comparison campaign （2nd EOP PCC）: Overview[J]. Artificial Satellites, 2022, 57 (s1): 1- 14.

12
张昊, 王琪洁, 朱建军, 等. 加权最小二乘法与AR组合模型在极移预测中的应用研究[J]. 天文学进展, 2011, 29 (3): 343- 352.

DOI

ZHANG H, WANG Q J, ZHU J J, et al. Joint model of weighted least-squares and AR in prediction of polar motion[J]. Progress in Astronomy, 2011, 29 (3): 343- 352.

DOI

13
SCHUH H, ULRICH M, EGGER D, et al. Prediction of Earth orientation parameters by artificial neural networks[J]. Journal of Geodesy, 2002, 76 (5): 247- 258.

DOI

14
WANG C, ZHANG P. Improving the accuracy of polar motion prediction using a hybrid least squares and long short-term memory model[J]. Earth, Planets and Space, 2023, 75 (1): 153.

DOI

15
YU K, YANG K, SHEN T, et al. Estimation of earth rotation parameters and prediction of polar motion using hybrid CNN-LSTM model[J]. Remote Sensing, 2023, 15 (2): 427.

DOI

16
BARNES R T H, HIDE R, WHITE A A, et al. Atmospheric angular momentum fluctuations, length-of-day changes and polar motion[J]. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 1983, 387 (1792): 31- 73.

DOI

17
BRZEZIŃSKI A. Polar motion excitation by variations of the effective angular momentum function: Considerations concerning deconvolution problem[J]. Manuscripta Geodaetica, 1992, 17 (1): 3- 20.

DOI

18
GROSS R S, FUKUMORI I, MENEMENLIS D. Atmospheric and oceanic excitation of the Earth’s wobbles during 1980—2000[J]. Journal of Geophysical Research: Solid Earth, 2003, 108 (B8): 2370.

DOI

19
DILL R, DOBSLAW H. Short-term polar motion forecasts from earth system modeling data[J]. Journal of Geodesy, 2010, 84 (9): 529- 536.

DOI

20
DILL R, DOBSLAW H, THOMAS M. Improved 90-day Earth orientation predictions from angular momentum forecasts of atmosphere, ocean, and terrestrial hydrosphere[J]. Journal of Geodesy, 2019, 93 (3): 287- 295.

DOI

21
MUNK W H, MACDONALD G J F. The rotation of the Earth: A geophysical discussion[M]. Cambridge: Cambridge University Press, 1960.

22
WILSON C R. Discrete polar motion equations[J]. Geophysical Journal of the Royal Astronomical Society, 1985, 80 (2): 551- 554.

DOI

23
BIZOUARD C, GAMBIS D. The combined solution C04 for Earth orientation parameters consistent with international terrestrial reference frame 2005[C]. IAG Symposium Munich, Berlin, Germany, October 9—14, 2006.

24
GROSS R S. The excitation of the Chandler wobble[J]. Geophysical Research Letters, 2000, 27 (15): 2329- 2332.

DOI

25
DICKMAN S R. Dynamic ocean-tide effects on Earth’s rotation[J]. Geophysical Journal International, 1993, 112 (3): 448- 470.

DOI

26
SALSTEIN D A, ROSEN R D. Regional contributions to the atmospheric excitation of rapid polar motions[J]. Journal of Geophysical Research: Atmospheres, 1989, 94 (D7): 9971- 9978.

DOI

27
PONTE R M. Oceanic excitation of daily to seasonal signals in Earth rotation: Results from a constant-density numerical model[J]. Geophysical Journal International, 1997, 130 (2): 469- 474.

DOI

28
WANG Q X, DU L, ZHOU X H, et al. Impacts of Earth rotation parameters on GNSS ultra-rapid orbit prediction: Derivation and real-time correction[J]. Advances in Space Research, 2017, 60 (12): 2855- 2870.

DOI

29
MMONTENBRUCK O, GILL E, LUTZE F H. Satellite orbits: Models, methods, and applications[M]. Berlin: Springer, 2000.

30
ENGE P K. The global positioning system: Signals, measurements, and performance[J]. International Journal of Wireless Information Networks, 1994, 1 (2): 83- 105.

31
RUSSIAN INSTITUTE OF SPACE DEVICE ENGINEERING. GLONASS interface control document[R]. Moscow: RNIIKP, 2008.

32
EUROPEAN UNION. European GNSS (Galileo) open service definition document[R]. Brussels: European Union, 2021.

33
CHINA SATELLITE NAVIGATION OFFICE. BeiDou navigation satellite system open service performance standard (Version 3.0)[S]. Beijing: China Satellite Navigation Office, 2021.

Options
文章导航

/

〈 〉