低轨航天器多点多参数轨道上行数据计算方法与流程

本发明属于航天测量与控制领域，涉及地面测控中心对低轨道航天器加工上注轨道数据的计算方法。

背景技术：

航天器在空间飞行过程中，为完成对地观测任务，需要有自身相对于地球的位置信息；同时，航天器的太阳帆板为了最大程度地利用太阳能，必须使太阳帆板的法向指向太阳，因此需要有航天器相对于太阳的位置信息；此外，航天器为了规划未来的使命任务，也需要具备自主的轨道预报能力。

但低轨航天器精确地完成自主轨道预报比较困难。首先，航天器在轨运行过程中受复杂摄动力的作用，建立精确的模型比较困难。运行于低地球轨道的航天器受大气阻力的影响，轨道高度逐渐衰减，轨道半长轴逐渐降低，轨道偏心率减小，轨道周期变短。航天器的大气阻尼力主要受大气密度的影响，大气密度模型与太阳辐射流量、地磁指数有关，航天器在轨运行过程中很难实时获得空间环境信息。航天器还受地球非球形引力摄动的影响，升交点赤经随时间变化，在惯性空间中表现为椭圆轨道面不断旋转，其平近点角在轨道面内的相位角也不断变化，在惯性空间表现为椭圆轨道近地点的进动。建立准确的引力模型需要考虑考虑地球引力场势函数中的高阶项，计算量大，运算复杂。其次，航天器的计算资源有限，星上计算机要负责数传、星务管理、姿控、载荷数据处理等任务，受星上计算机性能的限制，分配给轨道计算的资源不足以完成高精度的建模与轨道预报任务。

因此，通常低轨航天器需要地面上注加工好的轨道信息，星上计算机采用简单的动力学模型完成轨道预报。这种方法注入轨道数据的精度较差，轨道误差很容易随时间快速发散，因此要求轨道注入的频率高，不能满足卫星长期轨道预报的需求。

技术实现要素：

为了克服现有技术的不足，本发明提供一种多点多参数轨道上行数据计算方法，将地面测控中心完成的一组轨道预报结果以函数、函数初值和函数参数的形式上注给航天器，提高了航天器自主轨道预报的精度。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤一，地面测控站接收航天器的遥外测数据，通过数据网发送给地面测控中心，地面测控中心完成数据处理后进行精密轨道确定，生成一分钟一点的弹道文件，以格式[yyyymmddhhmisecxyzvxvyvz]存储，其中，yyyy、mm、dd、mi、sec分别为协调世界时的年、月、日、时、分、秒，x、y、z、vx、vy、vz分别为航天器在该时刻j2000坐标系下三个方向轴的位置与速度；

步骤二，将步骤一生成的弹道文件转换为j2000坐标系瞬时轨道根数，每一点弹道记录对应一组瞬时轨道根数，格式为[yyyymmddhhmisecaeiωωm]，其中，a、e、i、ω、ω、m分别为航天器轨道在j2000坐标系的半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角、平近点角；

步骤三，把步骤二生成的j2000坐标系瞬时轨道根数转换为拟平均轨道根数，包含二阶长期变化项的平根数构造如下：

生成的拟平根数存储为[yyyymmddhhmisec]；其中，j2、j3、j4分别为航天器环绕行星的引力场球谐函数中带谐项第二、三、四项的系数，为二体引力下的轨道根数，σ1(t-t0)、σ2(t-t0)分别为一阶、二阶长期变化项，为轨道平半通径，为归一化的轨道平角速度，μ为地球引力常数，t0为初始时刻，a0、e0、i0、ω0、ω0、m0、n0分别为参考椭圆轨道的半长轴、偏心率、倾角、升交点赤经、近地点幅角、平近点角、轨道角速度，a2、a3、a4分别为引力场球谐函数系数的简化符号；

步骤四，对步骤三生成的拟平根数进行处理，得到纬度幅角采用纬度幅角替换平近点角，将轨道倾角、升交点赤经、近地点幅角、纬度幅角转换为弧度的格式；

步骤五，判断纬度幅角的变化趋势，若纬度幅角变化趋势为增加，则纬度幅角的值穿越2π时在后续所有纬度幅角值上加上2π，若纬度幅角趋势变化趋势为减小，则纬度幅角的值穿越0时在后续所有纬度幅角值减去2π；每穿越一次均处理穿越点后面的所有记录点；

步骤六，采用差值拟合法计算q参数及其变率q参数是指纬度幅角的变化率除包括以下步骤：

计算l＝m+ω；

对n进行多项式拟合，得到n＝pn(1)t+pn(2)；

对n求积分，得到

将l与n的积分做差，得到

对q进行拟合，得到q＝pq(1)t²+pq(2)t+pq(3)，则即

步骤七，分别对轨道平半长轴和轨道角速度n进行线性拟合，得到平半长轴变率和平均轨道角速度轨道角速度变率分别对升交点赤经和近地点幅角进行二次多项式拟合，得到升交点赤经的一阶变率二阶变率和近地点幅角的一阶变率二阶变率对纬度幅角进行归一化差值拟合，得到q参数的一阶变率二阶变率

本发明的有益效果是：首先对地面精密轨道平均化，然后用函数拟合的方法得到平均轨道的拟合公式，以上注轨道参数的方式使航天器得到较高精度的轨道。示例分析表明，对低轨航天器每12小时上注一次轨道参数，轨道计算误差可保持在20米以内。因此，本发明既利用了地面轨道预报精度较高的优点，又减少了上注的数据量，同时满足航天器快速计算和精度较高的需求，具有较高的实用价值。

附图说明

图1是多点多参数轨道上行数据计算方法总体框图；

图2是纬度幅角的周跳修正示意图；

图3是拟合的各项指标曲线及拟合残差示意图，其中，(a)为半长轴随时间变化曲线，(b)为半长轴线性拟合残差，(c)为近地点幅角随时间变化率曲线，(d)为近地点幅角拟合残差，(e)为近地点幅角随时间变化率曲线，(f)为近地点幅角拟合残差，(g)为平近点角拟合曲线，(h)为平近点角拟合残差，(i)为纬度幅角随时间变化曲线，(j)为纬度幅角的拟合残差，(k)为升交点赤经随时间变化曲线，(l)为升交点赤经拟合残差，(m)为平均角速度随时间变化曲线，(n)为平均角速度的拟合残差；

图4是48小时弧段拟合残差示意图，其中，(a)为48h弧段纬度幅角的拟合残差，(b)为48h弧段近地点幅角与平近点角的拟合残差；

图5是24小时弧段纬度幅角拟合残差示意图；

图6是12小时弧段纬度幅角拟合残差示意图；

图7是6小时弧段纬度幅角拟合残差示意图；

图8是6小时弧段半长轴拟合残差示意图；

图9是6小时弧段升交点赤经拟合残差示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明提供的低轨航天器多点多参数轨道上行数据计算方法，总体流程如图1所示，包括以下步骤：

步骤一：地面测控站接收航天器的遥外测数据，通过数据网发送给地面测控中心，地面测控中心完成数据处理后进行精密轨道确定，生成一分钟一点的弹道文件，以j2000坐标系位置、速度的格式存储，格式如下：

[yyyymmddhhmisecxyzvxvyvz](1)

yyyy、mm、dd、mi、sec分别为协调世界时的年、月、日、时、分、秒，x、y、z、vx、vy、vz分别为航天器该时刻j2000坐标系的位置与速度。

步骤二：把步骤一生成的弹道文件转换为j2000坐标系瞬时轨道根数，每一点弹道记录对应一组瞬时轨道根数，格式如下：

[yyyymmddhhmisecaeiωωm](2)

a、e、i、ω、ω、m分别为航天器轨道在j2000坐标系的半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角、平近点角。

步骤三：把步骤二生成的j2000坐标系瞬时轨道根数转换为拟平均轨道根数，瞬平转换的算法主要考虑地球非球形引力摄动与大气阻尼摄动，通过求解航天器轨道运动的摄动方程，把解展开为小级数的形式，舍去其中的短周期项，保留长周期项、长期项。平根数可以表示为零阶项、一阶长期变化项、二阶长期变化项，乃至更高阶长期变化项之和，包含二阶长期变化项的平根数的构造如下：

生成的拟平根数存储为

其中各参数的含义如表1所示。

表1各参数含义

步骤四：对步骤三生成的拟平根数进行处理，得到纬度幅角其计算公式如下

并采用纬度幅角替换平近点角，并将轨道倾角、升交点赤经、近地点幅角、纬度幅角转换为弧度的格式。

步骤五：纬度幅角的归一化。由于纬度幅角在0～2π之间变化，为计算纬度幅角的变率，需要将纬度幅角归一化，主要的方法是判断纬度幅角的变化趋势，若纬度幅角变化趋势为增加，则纬度幅角的值穿越2π时在后续所有纬度幅角值上加上2π，若纬度幅角趋势变化趋势为减小，则纬度幅角的值穿越0时在后续所有纬度幅角值减去2π。每穿越一次均处理穿越点后面的所有记录点，处理前后的纬度幅角如图2所示。

步骤六：采用差值拟合法计算q参数及其变率q参数是指纬度幅角的变化率去除主要的常数项其计算方法为

子步骤6.1：计算l＝m+ω；

子步骤6.2：对n进行多项式拟合，得到n＝pn(1)t+pn(2)；

子步骤6.3：对n求积分，得到

子步骤6.4：将l与n的积分做差，得到

子步骤6.5：对q进行拟合，得到q＝pq(1)t²+pq(2)t+pq(3)，则q的导数为由子步骤6.4可知，

步骤七：采用下表的方法得到航天器轨道上行数据的各参数变率。

表2各参数变率的处理方法

本发明以采用某超低轨道航天器的上行轨道数据加工为例，进行计算验证。

接收航天器下传的星载gps数据，进行精密定轨。精密定轨的动力学模型中重力场模型采用32阶jgm-3地球重力场模型，大气模型采用msis-90，地影模型采用锥形地影，并考虑日月等第三体摄动，采取如下的初始轨道：

表3初始轨道

在精密定轨软件中生成时长为两天的轨道，并转换为平根数，画出平根数随时间的变化曲线，如图2所示，横坐标单位时间为30分钟，角度单位为rad，半长轴单位为m。对要求提供变化率的平根数进行拟合，其中，半长轴、平均角速度采用一次多项式函数拟合，升交点赤经、近地点幅角采用二次多项式函数拟合，纬度幅角采用归一化差值拟合算法。程序采用matlab-m语言，运算环境为matlab7.0。

通过以上数据的拟合与分析，得到如下的拟合公式：

选择弹道第一点为平根数注入需求点，相应的参数变率输出结果为

表4加工的轨道上行数据中的参数变率

由图2可知，在两天的时间长度内，近地点幅角拟合残差最大值为1.1×10^-3rad，平近点角的拟合残差最大值也近似为1.1×10^-3rad，两者相差正负号，曲线关于y＝0近似对称，因此两者的拟合残差绝大部分相互抵消，纬度幅角的拟合残差最大值为5.323×10^-5rad。因此采用二阶多项式的方法对纬度幅角进行推算，其位置误差最大为355m，主要是沿迹向。升交点赤经拟合的最大残差为2.5×10^-6rad，引起的位置误差最大为16.6米，主要是法向。半长轴拟合的最大残差为17米，主要是径向。

总体而言，采用二次多项式对卫星进行拟合所引起的位置误差为：径向20米以内，法向20米以内，沿迹向400m以内。最大误差发生的时刻大致在轨道参数注入的48小时，若每24小时注入一次轨道根数，其拟合残差曲线如图5所示，由于参与拟合的弧段变短，纬度幅角的拟合残差最大值为10^-5rad，可将星上拟合的沿迹向误差降低至48弧段误差的五分之一，约66.6m，若每12小时注入一次轨道根数，其拟合残差曲线如图5所示，可将纬度幅角拟合误差降低至3×10^-6rad以内，即星上拟合的沿迹向误差将保持在20m以内。若每6小时注入一次轨道根数，其拟合残差曲线如图6所示，可将纬度幅角拟合误差降低至2×10^-6rad，星上拟合的沿迹向误差将保持在13.3m以内，半长轴拟合误差如图7所示，径向误差将保持在4米以内，升交点赤经拟合误差如图9所示，法向误差将保持在6米以内。

可以看出采用本方法加工的轨道上行数据48小时内预报偏差小于400米，具有预报精度高，计算快捷，注入参数简单等优点，适合用于航天器中等精度的自主轨道预报。