基于精确星历最小二乘拟合的轨道快速递推方法*

车 征

（陕西星邑空间技术有限公司 西安 710043）

1 引言

高精度的轨道快速递推模型可为卫星敏感器、控制器等部件的在轨自主工作提供位置输入信息。由于星上计算资源和运算效率限制，无法使用高精度数值积分方式进行轨道预报。传统的摄动分析方法采用一阶或二阶摄动理论［1～2］，虽然能够进行较快速的平均根数与瞬时根数的转换和轨道递推，但需要计算大量的公式推导工作，且在摄动模型项较多时运算仍然较为复杂。GPS 星历外推使用了一种参数拟合的方式，将摄动规律转为有限个参数描述的近似项［3～5］。参数个数从最初的16参数、发展出18 参数、22 参数等［6～9］。这种方法具有很高的拟合精度，但缺点是仅仅适用于短时间内的位置解算，当拟合时段变长则误差急剧增大，且没有包含速度信息，不适合用于长时间轨道递推［10～13］。还有一种方法是多项式插值方式［14～15］，在拟合区间内能够实现很高精度，但多项式阶数较高时计算量较大，且用于预报存在发散趋势。

本文从瞬根数变化规律中提取轨道根数各分量的频率信息，使用含少量待定系数的三角函数和线性函数组合关系式描述轨道根数的变化趋势，使用精确动力学预报结果通过最小二乘拟合的方式确定系数值。避免了复杂的解析公式运算，同时保证了与精确动力学模型预报结果的吻合度。仿真算例表明拟合结果可以保证24h 平均预报误差小于400m，且计算耗时小。

2 拟合模型

基本思想为根据瞬根周期规律进行傅里叶分解，找出每个根数变化的短周期频率，分别指定拟合函数，对每个轨道根数单独拟合。以第一类无奇点根数为变量，定义各个分量的变化规律如下：1）a的周期为T/2，还包括随时间线性变化项。2）ex和ey的周期为T/3和T。3）i的周期为T/2和1/2天。4）Ω的周期为T/2，还包括随时间线性变化项。5）λ的周期为T和1/2天，还包括随时间线性变化项。T为平均轨道周期。因此设置拟合函数如下式：

一共有30 个参数。对六根数分别拟合，设待拟合的参数符号为x（代表a，ex，ey，i，Ω，λ中的一个），待拟合参数的初值为x0到xn，则：

利用不同时间共m 组x 的精确动力学模型预报数据，有

其中L=xpre-x为按照式（1）和时间序列预报的xpre与真实动力学积分得到的x之间的偏差。

则根据L可以反算出x0到xn应修正量和修正后的待拟合变量为

其中（a，ex，ey，i，Ω，λ）对各自待拟合变量的偏导数为

重复式（3）～（6）直到计算的x0到xn应修正量接近0，即表示待拟合变量收敛。由于最小二乘方法需要一个较为准确的初值才能稳定收敛，下节给出初值的计算方法。

3 拟合初值计算方法

3.1 a的系数初值计算

对于给定的数值积分获得的一系列时刻的a（ti），i=1 到N，首先根据大气衰减较小的假设，取a3=0。然后根据式（1），可知

当tf远大于T时，可以认为

其中Δt为数据时间间隔，且(N-1)ΔtΔn是一个小量，则

T作为常值用于随后各个轨道根数系数的拟合过程中。根据式（1）的性质，可知：

则：

3.2 ex，ey 的系数初值计算

则

求解方程组得：

同理利用

可计算：

ey的初值计算与式（15）到式（18）完全相同。

3.3 i的系数初值计算

同样地，i0初值直接取i 的平均值，。根据和i(T)的特性，有：

3.4 Ω的系数初值计算

根据：

可得：

3.5 λ 的系数初值

首先λ3初值取为式（10）中的轨道角速度n。

根据下式：

得：

再根据下式：

得：

4 测试结果

初始瞬根为［7130km，0，98°，0°，0°，360°］。则拟合一天的参数为

a=7120.9688+9.0790 sin(2πt/T+1.5751)-3.3803e-08t

ex=-0.00046978+0.0007439 sin(6πt/T+1.5752)+0.0002941 sin(2πt/T-1.5749)

ey=2.2291e-5+0.0007435 sin(6πt/T+0.004545)+0.0009327 sin(2πt/T-3.1340)

i=1.7105+8.9871e-05 sin(4πt/T-1.5839)+2.7833e-05 sin(2πt/43200-0.6091)Ω=2.0126e-05+9.0666e-05 sin(4πt/T-3.1393)+1.9014e-07t

λ=0.0001174+0.0009410 sin(2πt/T-3.0928)+0.001049t+9.7748e-05 sin(2πt/43200-2.1751)+9.4583e-05 sin(4πt/T+0.008591)

对应的残差换算到三维位置偏差见图1，平均总误差小于400m。设去除3 倍周期以内的短周期项后的剩余长期项（常数项、长周期项和时间线性项）为平根数，与瞬根的对比如图2～图7。

图1 位置残差

图2 半长轴长期项提取

图3 ex长期项提取

图4 ey长期项提取

图5 i长期项提取图

图6 Ω长期项提取

图7 λ 长期项提取

可见拟合获得的平根能够很好地反映轨道根数长期变化。

5 结语

本文使用三角函数和时间线性函数来逼近低轨近圆轨道的轨道根数变化趋势，给出了最小二乘拟合函数系数的算法和流程，并设计了快速获得拟合系数初值的计算方法。仿真实例表明，拟合方法能够快速获得较高精度的轨道根数递推函数，24h数据拟合精度达到400m，且计算量很小。算法可用于星上快速轨道递推和瞬平根数转换。