微小推力下小卫星的轨道递推与推力标定

张军华，陈建林，孙冲，袁建平

1. 西北工业大学 航天学院航天飞行动力学技术重点实验室，西安 710072

2. 西北工业大学 民航学院，西安 710072

1 引言

随着商业航天的发展，低轨小卫星平台受到越来越广泛地关注。小卫星相对于大卫星来说，由于体积小、质量轻、成本低，具有更灵活的应用潜力，因此近年来小卫星广泛应用于近地航天任务，已经成为航天工程的研究热点之一[1-2]。

通常，小卫星受限于质量约束，无法像大卫星平台一样携带大量的燃料，因此高比冲微小推力推进器小卫星的常用推进方式，主要包括电推进、太阳帆推进、电帆推进、冷气推进等[3-10]。随着航天技术的发展，当前的微小卫星任务对轨道预测和估计的精度要求不断提高，这不仅对卫星轨道动力学建模精度提出了较高的要求，还对轨道预测和估计方法提出了精度要求。用于轨道的预测和估计的数值积分方法，如Dormand-Prince方法、Gauss-Jackson方法、Adams-Bashforth-Moulton方法等往往在误差控制、稳定性等方面存在问题[11-15]。而采用简化的动力学模型往往导致轨道预测的精度无法满足实际需求[16-19]。

进一步，由于微小推力的量级小，难以严格在地面完成高精度的推力标定，需要卫星在轨对其推力进行标定。常见的推力标定包括参数标定、飞轮标定、矢量标定等[20-23]，这些标定方法采用卫星载荷对标定所需参数进行估计。但小卫星搭载载荷有限，上述标定方法无法适用于小卫星的推力在轨标定。另一方面，微小推进器存在推力损失等不确定性因素，如何保证推力在轨标定的精度，是个亟待解决的问题。

本文针对小卫星的高精度轨道预测和微小推力在轨标定问题，对主要的摄动力进行了高精度建模。利用变步长龙格库塔积分RK78对卫星轨道进行了外推，通过与STK仿真结果对比验证了模型的精确性；使用无迹卡尔曼滤波器(UKF)同时估计小卫星轨道并标定卫星推力，证明了在轨估计算法的可行性。本文对轨道递推和推力标定算法的研究，可以应用于小卫星的实时在轨运算中，实时掌握小卫星自身的轨道状态和推力大小，同时算法的速度较快，有利于小卫星在线运算任务的实现。

2 微小推力下小卫星动力学模型

2.1 问题描述

小卫星在近地轨道上运动时，除了受到地球中心引力场的引力外，还受到地球非球形引力、太阳光压和第三体引力等摄动影响，小卫星采用微小推力进行轨道控制。对微小推力下小卫星的在轨运行情况进行预测，不仅需要精确的轨道运动模型，而且对轨道预测和估计的精度也有较高的要求。此外，由于微小推力量级小，且小卫星任务过程中存在不确定性，需要将微小推力进行在轨标定。

2.2 动力学模型

首先建立三个坐标系[24]，如图1所示。

图1 坐标系示意图Fig.1 Illustration of coordinate frames

地球惯性坐标系EXYZ(JDate)，坐标系原点在地心E，X轴指向当前时刻春分点，Z轴垂直于地球赤道平面向上，Y轴满足右手法则。

地球惯性坐标系J2 000，坐标系原点在地心E，X'轴指向J2 000时刻的平春分点，Z轴垂直于地球赤道平面向上，Y'轴满足右手法则。

卫星轨道坐标系Sxyz，坐标系原点在卫星质心S，x轴沿卫星矢径方向，从地心指向卫星，y轴垂直卫星矢径并指向速度方向，z轴垂直于轨道平面，遵循右手法则。

除此3个坐标系外，建立地球固联坐标系(ECEF)，其定义为：坐标系原点在地心，X''轴从地心指向赤道与国际参考子午线的交点，Z''轴沿着地球自转轴从地心指向北极点，根据右手准则Y''轴从地心指向赤道与90°东经子午线的交点。

另外，假设卫星本身采用三轴对地稳定，因此卫星的本体坐标系与轨道坐标系是相重合的。

卫星受到微小推力加速度aF的作用，且该微小作用力在本体坐标系中的方向不变，如图2所示。因此，aF在轨道坐标系下可以由θ和φ两个方位角表示，aF=[aFcosφcosθ,aFcosφsinθ,aFsinφ]T。

图2 微小推力在轨道坐标系下的示意Fig.2 Illustration of micro-thrust in frame Sxyz

卫星受到微小推力、地球扁率(J2、J22等)、太阳光压以及第三体(日月)引力等摄动，则卫星在地心惯性坐标系下的动力学方程可以描述为：

(1)

(2)

式中：ax,nm、ay,nm和az,nm分别表示地球形状摄动加速度在地球固联坐标系中三个方向上的分量：

(3)

式中：R⊕表示地球半径；地势系数Cnm,Snm由EGM96S引力模型给出；Vnm和Wnm分别为：

(4)

(5)

在地球惯性系下的太阳光压和太阳、月球等第三体引力引起的摄动力为：

(6)

(7)

式中：r⊕,rSun,rMoon,r分别表示地球惯性系下卫星到太阳的位置矢量；太阳的位置矢量；月球的位置矢量和卫星的位置矢量；G表示万有引力常数；AU表示地球到太阳的距离；P表示单位AU距离下的太阳光压力；ν表示地球的阴影函数；M表示太阳或者月球的质量；Cr表示辐射压力系数；m表示卫星的质量；A表示卫星的受照横截面积。

因此，其他摄动总加速度可以表示为：

(8)

式中：aFx、aFy、aFz、aJx、aJy、aJz、adx、ady、adz分别表示微小推力加速度，地球形状摄动加速度，其他摄动加速度在惯性坐标系下的三个坐标轴的分量。

3 卫星轨道递推算法

卫星常见的星上定轨方式是基于GNSS信号，但是由于GNSS输出的轨道信息一般在位置和速度上存在误差，直接应用GNSS信号进行轨道确定，会造成卫星轨道定位不精确。为了解决这一问题，提出一种变步长龙格库塔算法RK78进行轨道外推，具体的算法模型如下：

Runge-Kutta方法是一种间接采用泰勒级数展开而求解常微分方程初值问题的数值方法。该方法用某些点上函数值的不同组合来近似初值问题的解析解，从而避免函数高阶偏导数的计算，并获得高阶的方法[26]。将Runge-Kutta公式与初值问题解的泰勒展开式进行比较，使之有某些项相同，便可以确定组合系数，从而建立所需的方法。

由m个函数值线性组合构成的显式Runge-Kutta公式为：

(9)

式中：K1和Ki分别表示时间段开始时的斜率和每一时间段中点的斜率；ci,αi,βij均为待定系数；h为步长。可以看出，方程(9)含有m个函数值，所以称为m阶Runge-Kutta公式。本文使用的RK78算法表示采用七阶-八阶变步长Runge-Kutta算法，它用七阶方法提供候选解，八阶方法控制误差，是一种自适应步长的常微分方程数值解法，其整体截断误差为O(h8) 。其变步长的原理如下：

对于一步h1对应其误差估计Δ1，而h0对应其误差估计Δ0，两者关系为：

(10)

针对两个误差估计值，如果Δ1>Δ0，可以根据方程(10)确定所要减小步长的大小；如果Δ1<Δ0，可以根据方程(10)确定能够在安全范围内增加步长的大小。

针对微小推力轨道外推问题，根据上述Runge-Kutta公式，使用动力学模型(8)带入到七阶Runge-Kutta公式(9)中。在此基础上，设置合适的变步长，利用八阶Runge-Kutta公式设置截断误差，通过积分和迭代，对卫星轨道进行外推。

4 基于UKF的推力在轨标定算法

微小推力由于在轨保持时间较长，为了保证微小推力在轨工作的有效性，需要对推力进行在轨标定。本文采用UKF对微小推力进行高精度推力在轨标定。

4.1 UKF的基本原理

无迹卡尔曼滤波(Unscented Kalman Filter，UKF)，是无迹变换(Unscented Transform，UT)与标准卡尔曼滤波体系的结合，通过无迹变换使非线性系统方程适用于线性假设下的标准卡尔曼体系，该算法的核心思想是采用UT变换，利用一组Sigma采样点来描述随机变量的高斯分布，然后通过非线性函数的传递，再利用加权统计线性回归技术来近似非线性函数的后验均值和方差[27]。

UT变换的特点是对非线性函数的概率密度分布进行近似，而不是对非线性函数进行近似，即使系统模型复杂，也不增加算法实现的难度；而且所得到的非线性函数的统计量的准确性可以达到三阶；除此之外，它不需要计算雅可比矩阵，可以处理不可导非线性函数。

算法模型如下：

设非线性系统：

(11)

式中：Xk为nx维的系统状态向量；Zk为nz维的系统观测向量；vk为系统噪声；协方差矩阵为Pv；nk为观测噪声；协方差矩阵为Pn；式中：vk、nk都是高斯白噪声，且互不相关。算法具体计算步骤如下：

1)选定滤波初值。

(12)

2)计算k-1时刻的2n+1个σ(Sigma)样本点，构造(2n+1)维向量χi为σ向量：

(13)

λ=α2(n+k)-n

(14)

3)状态更新，计算k时刻的一步预测模型值。

(15)

4)计算k时刻的一步预测增广样本点。

(16)

5)观测更新，计算k时刻的观测变量的一步预测均值、方差和协方差。

(17)

6)计算增益矩阵。

(18)

7)计算滤波值。

(19)

4.2 在轨标定算法

利用UKF对微小推力的大小及方向进行在轨标定的过程如下：

(1)动力学模型改进

将笛卡尔动力学模型进行状态扩维，即包含带估计的参数。

(20)

(2)基本方程

利用上一部分的轨道递推算法，对一个测量周期内的航天器状态进行预测，写出预测方程；然后根据UKF基本原理中的公式分别写出状态估值方程、滤波增益方程、一步预测均方误差方程和估计均方误差方程等。

(3)更新

1)设置位置误差和速度误差，设计GNSS数值生成模拟器。

2)融合GNSS的测量值，在UKF滤波器内对状态预测进行更新。实时估计状态x、m、aF、φ、θ。

5 仿真校验

5.1 轨道递推仿真结果

针对微小推力轨道外推问题，使用变步长龙格库塔RK78算法进行卫星轨道预测。初始参数设置如下：

1)设置初始轨道为圆轨道，参数轨道半径为755 3 km，轨道倾角86.5°，偏心率0，升交点赤经40°，近地点幅角60°，真近点角30°，地球引力常数3.986×1014m3/s2。通过轨道要素变换，可以求解卫星的初始位置和速度，如表1所示。

表1 初始条件

2)设置卫星质量为194 kg，微小推力12 mN，平均分配到三轴推力[6.928 2,6.928 2,6.928 2]TmN，推力加速度为[3.571 24,3.571 24,3.571 24]T×10-5m/s2。此外，考虑到微小推力较小，卫星推进器的质量流量忽略不计。

3)设置太阳光压摄动参数：辐射压力系数为1.21，卫星表面积为3.88 m2，地球到太阳的距离AU=149 597 870 000 m，单位AU距离下的太阳光压力P=4.56×10-6N/m2。

4)设置太阳引力摄动太阳引力常数为：1.327 124 38×1020m3/s2，月球引力常数为：4.902 793 455×1012m3/s2。

利用C++软件编写整个推演过程的程序，和STK软件模拟卫星轨道推演24 h，比较不同情况下所编写的C++程序和STK的结果对比。根据对比结果，有如下结论：

1)由表2～表6可以看出，在完整考虑4种主要轨道摄动和微小推力作用下，本文所设计的C++程序和STK模拟结果具有较高的一致性，24 h轨道递推的位置误差均在22 m以下，速度误差小于0.012 m/s，几乎可以忽略不计。

2)由表2和表3对比可以看出，在仅考虑微小推力情况时(不考虑其他轨道摄动)，本文所设计的C++程序和STK模拟的结果具有高度一致性，两种方法得到的24 h轨道递推的位置差异小于0.3 m，速度差异小于0.001 m/s。而由表4和表2对比可以看出，当考虑4种主要轨道摄动时，两种方法得到的24 h轨道递推的位置差异增加至30 m量级，速度差异增加至0.01 m/s量级。

表2 无摄动+无推力结果对比

表3 无摄动+12 mN推力结果对比

表4 考虑4种主要轨道摄动+无推力结果对比

3)由表3与表6对比可以看出，地球形状摄动力是引起本文所设计的C++程序和STK模拟结果差异的最大因素，位置差异为22 m左右，主要原因可能是采用的地球引力模型之间存在轻微差异，其他三种主要摄动力并不会引起另种方法之间的较大位置差异。

表5 考虑4种主要轨道摄动+12 mN推力结果对比

表6 仅考虑地球扁率摄动+12 mN推力结果对比

总之，尽管本文所设计的C++程序对地球形状摄动力的处理与STK内置的地球形状摄动力有所差异，但本文所设计的C++程序仍可实现在微小推力作用下的高精度轨道预测。

为了验证其他影响因素对所设计的C++程序的影响，本文分别针对轨道高度、推力大小和推演时间对程序的结果进行分析。

图3和图4分别给出了不同轨道高度下(200 km～1 500 km)C++程序和STK模拟之间的位置误差和速度误差，其中除轨道高度外，其余初始参数均与本节第一部分一致，程序考虑12 mN推力且无摄动下的结果。通过结果可以看出，不同轨道高度下，程序的位置误差均在0.33 m以下，速度误差均在4×10-4m/s以下，这说明本文的程序在不同轨道高度下均可以保证高精度轨道递推。

图3 不同轨道高度下位置误差Fig.3 Position errors at different orbital altitudes

图4 不同轨道高度下速度误差Fig.4 Velocity errors at different orbital altitudes

图5和图6分别给出了不同推力大小下(0 mN-20 mN)C++程序和STK模拟之间的位置误差和速度误差，其中除推力大小外，其余初始参数均与本节第一部分一致，程序考虑存在推力且无摄动下的结果。通过结果可以看出，不同推力大小下，程序的位置误差均在0.6 m以下，速度误差均在6×10-4m/s以下，这说明本文的程序在不同推力大小下均可以保证高精度轨道递推。

图5 不同推力大小下位置误差Fig.5 Position errors with different thrusts

图6 不同推力大小下速度误差Fig.6 Velocity errors with different thrusts

图7和图8分别给出了不同推演时间下(1～14 d)C++程序和STK模拟之间的位置误差和速度误差，其中除推演时间外，其余初始参数均与本节第一部分一致，程序考虑存在4种摄动和12 mN推力下的结果。通过结果可以看出，随着推演时间的增加位置误差和速度误差都会变大，14 d内的位置误差在140 m以下，位置误差在0.35 m/s以下，这说明本文的程序随着推演时间的增大，误差会累积增大，但还能保证正常的精度要求。

图7 不同推演时间下位置误差Fig.7 Position errors with different propagation time

图8 不同推演时间下速度误差Fig.8 Velocity errors with different propagation time

5.2 在轨标定仿真结果

根据上文的算法，利用C/C++对微小推力大小和其本体坐标系下的方向在3 h内进行在轨标定的仿真验证。

针对本文中的微小推力在轨标定问题，使用无迹卡尔曼滤波器(UKF)进行卫星轨道估计和推力标定，通过预测和更新，完成对微小推力的标定。卫星轨道递推的初始参数与5.1节第一部分初始参数一致。

利用所编写的C++程序在3 h时间间隔内对卫星微小推力进行推力在轨标定，初始位置各分量误差为1 000 m，速度各分量误差为0.05 m/s，微小推力各方向初始误差为0.000 01 m/s2，测量精度为1 m，滤波器更新时间间隔为10 s。表7展示了一次轨道标定结果，可以看出位置误差、速度误差和推力误差均收敛，推力精度占比收敛至0.429 4%。

表7 推力在轨标定3 h结果

如图9～图12以图示的形式给出了UKF推力在轨标定的结果，图9展示了3 h内的位置估计误差，图10为速度误差，图11为推力误差，图12为推力精度占比，由上图可以看出，各误差和占比均收敛。

图9 推力在轨标定的位置误差Fig.9 Position errors of on-orbit thrust calibration

图10 推力在轨标定的速度误差Fig.10 Velocity errors of on-orbit thrust calibration

图11 推力在轨标定的推力误差Fig.11 Thrust errors of on-orbit thrust calibration

图12 推力在轨标定的推力精度占比Fig.12 Thrust accuracy ratio of on-orbit thrust calibration

此外，通过改变C++程序中的位置误差、速度误差和位置误差，该程序仍能保证各误差和占比均收敛。

6 结论

对于微小推力下的小卫星轨道递推和推力在轨标定问题，本文提出了基于变步长龙格库塔的轨道递推算法，利用高精度的轨道动力学模型，能够对小卫星轨道进行高精度的推演。在此基础上，设计基于UKF推力标定算法，对小卫星施加的微小推力进行在轨标定。仿真实例中对轨道递推算法与STK进行比较，位置误差均在22 m以下，速度误差均在0.02 m/s以下；同时推力在轨标定各误差均收敛，且精度得到保证，仿真数据说明了本文所设计算法的可行性。本文设计的算法基于C++程序，相较于其他程序算法具有更快的计算速率，有利于实际任务中实现星上的在线计算。在下一步的研究中，致力于消除地球形状摄动力对算法带来的误差。