基于EKF和多信息源融合的空间交会对接过程实时定轨方法

马鹏斌 王 丹

（1中国西安卫星测控中心 2宇航动力学国家重点实验室）

1 引言

我国天宫一号目标飞行器和神舟八号飞船将要进行空间交会对接。其测量手段有地（海）基测量站USB测量（测距、测速、测角），中继星4程测距数据，目标飞行器和神舟八号飞船间的激光测距和雷达测距等外测手段。另外在飞船上还搭载有高精度的加速度计。可通过这些数据融合来快速确定天宫一号目标飞行器和神舟八号飞船的轨道。

2 系统方程

基于J2000惯性坐标系，对于航天器的运动方程，采用位置矢量r和速度矢量来描述。运动方程为=f（r，，t），可将多个航天器的轨道同时作为待估量，如目标飞行器（r1，）和飞船（r2，）两个航天器可同时求解。另外还可增加中继星作为待估量。

在轨道确定中需要求解的量往往不限于航天器状态r、，如求解大气阻尼系数CD，测量系统误差等。设其它需要求解的模型参数为c，因c为常数，所以有=0。

令

航天器运动状态方程可扩展成如下形式：

代表了一个由n个非线性一阶常微分方程组成的系统[1]。

3 动力学模型

航天器运动方程基于J2000惯性坐标系，位置矢量r和速度矢量描述。运动方程即为初值问题[1]

右函数f可分为f0和fε两部分，f0为二体问题下航天器的加速度，fε为其它各种摄动力加速度之和，可表示如下：

针对序贯处理的需求，对运动方程采用了单步法RKF7（8）方法求解。天宫一号和神舟八号可同时积分求解。

一般情况下，各种摄动力的量级如下表所示[2]：

表1 摄动力量级表

兼顾计算效率，只考虑米级以上的摄动影响，包括地球引力场、日月引力、太阳光压、大气阻尼摄动。由于神舟八号飞船搭载有高精度的加速度计，可以只计算地球引力场和日月引力，交会对接过程中的轨道机动推力和其他非引力摄动力可直接使用加速度计实测数据。

系统动力学模型的状态转移矩阵为：

在计算传统动力学模型的状态转移矩阵时，考虑地心引力和J2项摄动即可完全满足精度要求。考虑地心引力和J2项摄动，则有：

4 测量模型

建立量测方程，将系统在某时刻的测量量和系统状态量联系起来。状态量为J2000惯性坐标系下的位置矢量、速度矢量 r，。测量数据有：陆（海）基测量站USB测量（测距、测速、测角），目标飞行器和神舟八号飞船间的激光测距和雷达测距，中继星4程距离和等。

4.1 陆（海）基测量站USB测量

陆（海）基测量站USB测量在测站地平坐标系中的观测量为：斜距R、斜距变化率、方位角A、俯仰角E。设测站地平坐标系中的直角坐标为（x，y，z，，，），则：

观测量对状态量的偏导数为：

上述各式中出现的B和L分别为观测站的大地纬度和经度，R测站为J2000惯性坐标系下测站地心向径矢量。

4.2 目标飞行器和飞船间距离测量

令航天器 1 的位置矢量为r1=（x1，y1，z1），航天器 2的位置矢量为r2=（x2，y2，z2），则航天器 1 对航天器2的测距为：

ρ对航天器1，2状态矢量的偏导分别为：

4.3 中继星四程距离和测量

已知中继星精密星历，可将四程距离和处理为中继星对所求解航天器的测距数据，认为中继星是一个地面站，其测量模型同USB测距模型相同。

5 序贯估计方法

5.1 扩展卡尔曼滤波基本原理[3]

设状态方程为（t）=f[x（t），t]+w（t），w为模型噪声；量测方程为y（t）=h[x（t），t]+R（t），R为观测噪声，线性化后可得xk=Φ（tk，tj）xj+wj，yk=Hkxk+Rk。

已知估值及其协方差矩阵Pk，对tk+1时刻的一个新观测数据,卡尔曼滤波基本方法如下：

计算tk+1时刻的预测轨道；

计算tk+1时刻的预测协方差矩阵

进行状态更新，得到tk+1时刻的估计值

计算tk+1时刻的误差协方差矩阵Pk+1=

对运动方程重新初始化，把运动方程中原参考轨道换成最新得出的tk时刻的最优估值。把运动方程积分至tk+1时刻以处理tk+1时刻的观测数据。

在卡尔曼滤波实际应用中很可能出现的一个现象是处理了一定量的数据后滤波值出现发散现象。每次在协方差矩阵上加一个小的噪音项Q（即过程噪声）可以抑制矩阵变得越来越小[4]。则有：

5.2 野值判断

野值剔除是序贯定轨中很重要的环节，一个大的野值可能导致序贯定轨的发散。

设得到t时刻的测量量为O，其测量方差为R。预报至t时刻的轨道状态矢量为X，先验方差矩阵为P，C为测量量的预报值。

先验方差矩阵P的对角线元素可以近似认为是轨道精度，则

可近似认为：

测距预报值精度SR为Rp

测角预报值精度SAE为arctan（Rp/CR）

测速预报值精度SD为（P11P44+P22P55+P33P66）/CR

O-C为观测值减计算值。若满足：

O-C＞n1·R 且 O-C＞n2·S

则测量值为野值。其中，R为测量方差，S为测量值预报精度。n1，n2为人为设定门限，若使用5σ剔除，则n1，n2为5。如果某一时刻的测量量只有某一项为野值，则剔除该量，其他量仍可使用。

6 实时软件流程

实时软件流程如下：

第一步：计算起步初值：

第二步：得到新的tk+1时刻观测数据后，进行轨道外推，计算tk+1时刻的预测轨道；

第三步：计算tk到tk+1时刻的状态转移矩阵Φ（tk+1，tk）；

第四步：计算tk+1时刻的预测协方差矩阵=Φ（tk+1，tk）PkΦT（tk+1，tk）；

第五步：计算观测量计算值；

第六步：计算观测偏导数；

第七步：进行数据质量检验，判断O-C是否满足规定门限要求，或O-C与预测O-C相差满足规定值，如不满足则剔除此值，返回第二步；

第九步：进行状态更新，得到tk+1时刻的估计值

第十步：计算tk+1时刻的误差协方差矩阵Pk+1=

第十一步：比较Pk-1、Pk与Pk+1，判断是否发散，如连续发散则返回第一步；

第十二步：输出航天器轨道状态量，也可预报下一时刻引导值并输出；

第十三步：返回第二步，如数据结束则转到下一步；

第十四步：进行结束处理。

7 结论

本文针对我国空间交会对接任务的快速轨道计算，融合使用USB测量数据、中继星测量数据、目标飞行器和飞船间测量数据以及船载高精度加速度计数据等，利用扩展卡尔曼滤波方法来快速确定目标飞行器和飞船的轨道。可用于我国空间交会对接任务的实时轨道监视。 ◇

［1］李济生.人造卫星精密轨道确定.北京：解放军出版社,1995

［2］刘林.人造地球卫星轨道力学.北京:高等教育出版社,1992

［3］贾沛漳，朱征桃.最优估计理论.科学出版社.1984.

［4］A.C.Long et al.Goddard Trajectory Determination System.GSFC.1989