一种分布式信息融合的航天器实时轨道确定算法

黄静琪，何雨帆，孙山鹏

(1. 西安卫星测控中心宇航动力学国家重点实验室，西安 710043；2. 火箭军工程大学，西安 710025)

0 引 言

航天器实时定轨在空间态势感知、航天器实时跟踪监视、碰撞预警、航天器机动效果快速评估等方面有着重要作用。实时定轨的测量数据一般由多个相同或者不同类型的传感器组成的传感器网络提供，因此如何融合多传感器的测量数据快速准确地进行轨道确定十分关键。

目前，我国航天测控网主要采用集中式滤波算法[1-3]进行实时定轨，各测量数据均传输至轨道计算中心进行融合计算。随着航天技术应用的不断发展，对实时定轨系统可靠性的要求也越来越高。集中式算法过分依靠中心节点使的鲁棒性差。同时，各传感器测量数据均要传输至中心节点，导致通信量大、中心节点计算负担高、扩展性较差，分布式算法可有效解决以上问题[4-6]。

分布式状态估计方法主要分为两类：基于一致性算法和基于融合算法。对于一致性算法，Olfati-Saber在文献[7]中将平均一致性算子加入节点的局部卡尔曼滤波中，提出了卡尔曼一致性滤波(Kalman consistency filter,KCF)算法，该算法要求每个传感器节点和邻居节点对目标有联合可观性。随后，Olfati-Saber在文献[4]中首次提出了基于状态一致性的最优分布式卡尔曼滤波，但该方法要求每个节点已知邻居的邻居节点信息，且计算量大，无法应用到实际中。为了降低该方法的通信和计算复杂度，文献[4]中给出了一种次优的分布式卡尔曼滤波器。以文献[4]为基础，针对不同问题发展出许多分布式卡尔曼滤波算法。上述一致性算法通常要求各节点在局部更新前与邻居节点进行多次通信。除了一致性算法外，文献[8]提出了基于融合方法的卡尔曼滤波(Distributed Kalman filtering,DKF)，文献[9]将协方差交叉算法[10](Consensus-plus-innovations,CI)与DKF算法结合，提出了一种基于CI算法的DKF算法(CI-DKF)，但该方法普适性不强，只能用来解决文中提到的两种特定场景下的分布式滤波问题。

上述分布式算法针对不同的问题提出，常常有特定的假设，对节点及其邻居节点的联合可观性、通信拓扑的强连通性、动力学模型和通信拓扑的时不变性等都有严格的要求，有的还要求局部节点需要知道与通信拓扑结构相关的全局信息，不利于实际应用。文献[11]在上述研究的基础上，提出了一种分布式混合信息滤波算法(Distributed hybrid information filtering,DHIF)，该算法以扩展信息滤波为基础，与协方差交叉算法(CI)结合，是一种普适性较强的分布式滤波算法，通信代价低，并且可以使每个节点的估计值逼近全局最优。

本文以DHIF算法为基础，结合我国航天器测控系统现状，设计了一种分布式实时轨道确定算法，并对测量网络中不同传感器类型和不同网络拓扑结构进行了仿真计算，结果表明该算法使得各测站定轨精度均优于单站滤波精度，并逼近集中式滤波定轨精度。最后，本文分析了该分布式实时定轨算法收敛速度与通信拓扑的关系。

1 信息滤波算法

1.1 扩展信息滤波

在实施扩展卡尔曼滤波(Extended Kalman filter,EKF)的过程中，若迭代更新的不是协方差矩阵P，而是协方差矩阵P的逆矩阵，则该滤波过程称作扩展信息滤波(Extended information filter,EIF)。

设有非线性系统：

x(k)=f(x(k-1))+w(k-1)

(1)

z(k)=h(x(k))+v(k)

(2)

式中:x(k)为k时刻的状态向量;z(k)为k时刻传感器的测量值。w(k)和v(k)是零均值不相关白噪声，且协方差矩阵分别为Q(k)和R(k)。

1) 状态预测

(3)

Pk|k-1=F(k-1)Pk-1|k-1(F(k-1))T+Q(k-1)

(4)

Ik|k-1=(Pk|k-1)-1

(5)

2) 状态更新

Ik|k=Ik|k-1+(H(k))T(R(k))-1H(k)

(6)

K(k)=(Ik|k)-1(H(k))T(R(k))-1

(7)

(8)

(9)

式(9)乘Ik|k后，将(6)代入式(9)右侧，整理后可得：

(10)

1.2 多传感器的集中式信息滤波

对于多传感器的集中式信息融合滤波，每个传感器节点有各自的测量方程，k时刻第i个节点的量测方程可以表示为：

(11)

集中式信息滤波算法如下：

1)信息状态预测

状态、协方差矩阵和信息矩阵的预测同式(3)、(4)、(5)。信息状态为：

(12)

2)信息状态更新

根据文献[2]的集中式融合方法，系统信息状态及信息矩阵的更新为：

(13)

(14)

其中,i=1，…，N，N为网络中传感器的数量

(15)

(16)

2 基于传感器网络的分布式滤波

在航天器实时轨道确定这类问题中，往往有多个测站，各测站可能使用不同类型的传感器对目标进行跟踪测量，可以归为基于传感器网络的融合滤波定轨问题。目前航天器实时轨道确定系统中使用的集中式融合滤波方法依赖于计算中心，一方面需要各传感器向中心传递大量数据，增加了通信负载和计算量，另一方面中心节点故障将直接导致整个系统瘫痪。使用分布式滤波方法定轨，各测站既是测量节点也是计算节点，均可以对目标状态进行估计。

现有的大多数分布式算法中，针对不同背景问题，有较为严格的应用条件。例如，单个节点需要知道某些全局信息(整个网络中的节点数量或通信拓扑的最大入度等)、通信拓扑不能发生改变、每个节点和邻居节点是联合可观的等。DHIF算法则克服了上述缺点：1)每次局部状态更新前，邻居节点之间只需进行一次信息交换，通信量低；2)节点间交换的信息只有局部状态估计及量测信息，单个节点无需知道整个网络的全局信息，方便随时增减传感器节点，实现了传感器在网络中的即插即用；3)DHIF算法应用条件宽松，不需要每个节点及其邻居都是联合可观的，只要满足传感器网络是联合可观的，即可保证各节点的局部估计误差有统一下界，并且误差的期望值渐进为零[10]；4)算法支持时变的通信拓扑。

DHIF算法如下：

1) 局部状态预测

(17)

(18)

2) 局部状态更新

在k时刻，节点i将自身以及接收到的邻居节点的信息融合进行局部状态更新。为了保证各节点只通过接收邻居节点信息就可以获取到网络中其他节点的信息，网络的通信拓扑要求为连通图。

融合的信息包括两部分，第一部分为自身及邻居节点信息矩阵和信息状态的局部先验估计，第二部分为自身及邻居节点的y(k)和s(k)值。

首先，第一部分的融合采用了CI[10]算法，如下：

(19)

(20)

(21)

(22)

3 分布式滤波实时轨道确定

3.1 动力学模型

(23)

式中：μ为地球引力常数，Re为地球半径。式(23)可以表示为如下形式

(24)

将其离散化，其他摄动看作系统噪声w，T为离散时间间隔，可得

x(k)=x(k-1)+A(x(k-1))T+w(k-1)

(25)

式(25)可写为状态方程(1)的形式，线性化可得状态转移矩阵为

(26)

其中，I为3×3的单位矩阵，且

(27)

3.2 量测模型

(28)

(29)

在测站坐标系中有：

(30)

(31)

(32)

(33)

式(30)～(33)可以表示为量测方程(2)的形式，线性化后可得量测矩阵为：

(34)

3.3 分布式实时轨道确定算法

在分布式实时轨道确定中，所有测站构成一个传感器网络，每个观测站看作一个节点，算法采用DHIF，各测站分别对航天器进行实时定轨。计算流程如下：

图1 分布式实时轨道确定算法流程Fig.1 Distributed real-time orbit determination algorithm

1)步骤1:测站状态初始化

2)步骤2:状态预测

3)步骤3:获取量测值

4)步骤4:测站i计算状态信息和量测信息

(35)

(36)

(37)

(38)

5)步骤5:测站信息交换

7)步骤7:量测更新

(39)

(40)

4 仿真结果与分析

本节对分布式实时定轨算法进行仿真计算。假设有4个测站构成的传感器网络对一航天器进行分布式实时轨道确定。航天器初始轨道见表1，定轨弧段为2007-7-1 12∶32∶19至2007-7-1 12∶43∶59，使用STK高精度轨道外推模型(HPOP)将初轨外推700 s数据作为精密星历，观测站坐标见表2，各地面站跟踪航天器弧段见表3，可以看出航天器不同时对所有测站可见。

表1 目标航天器初始轨道Table 1 Target spacecraft initial orbit

表2 测站坐标Table 2 Coordinate of the observation stations

表3 各测站跟踪弧段Table 3 The observation arc of each station s

计算采用200次蒙特卡洛仿真，分别使用图2两种通信拓扑结构计算。

图2 两种不同的通信拓扑结构Fig.2 Two different communication topologies

4.1 测站测量方式相同

设4个测站均采用统一S频段系统(USB)设备测量，观测数据间隔1 s。

设初始状态为

初始状态协方差矩阵为

系统噪声方差矩阵为

Q=diag(0.12,0.12,0.12,0.012,0.012,0.012)

量测噪声方差矩阵为

Ri=diag(102,0.052,0.022,0.022)

试验结果见图3～图4。其中，图3 (a)、(b)分别为在图2两种通信拓扑下，各测站80 s内分布式滤波定轨位置和速度的均方根误差(Root mean squared error,RMSE)[1]曲线。

图3 各测站80 s内分布式滤波定轨的位置、速度RMSE曲线Fig.3 80 s position and velocity RMSE of DHIF for each station

结果表明，在初轨误差较大的情况下(位置误差13.69 km，速度误差52.1 m/s)，该分布式滤波定轨算法可以使各测站的计算结果快速收敛，且定轨结果趋于一致。根据图3，在前65 s只有测站S1、S2有量测数据(见表3)的情况下，各测站50 s定轨位置误差达到200 m，速度误差下降到10 m/s。在65 s时，测站S3开始跟踪目标，此时传感器网络中三个测站具有量测数据，这些量测信息通过相邻测站之间的通信在网络中传播。从RMSE曲线斜率的变化趋势可以看出，各测站定轨收敛速度加快。80 s时，定轨位置误差达到10 m以内，速度误差在0.1 m/s以内。

从式(39)、(40)可以看出，对某一测站，其相邻测站的状态信息和量测信息直接影响该测站的滤波结果，不同的通信拓扑会导致测站的邻居节点发生变化，继而影响该测站计算收敛的速度。对比图3(a)、(b)，在各测站观测弧段相同的情况下，拓扑a的收敛速度要快于拓扑b。

图4为各测站使用不同滤波方法在整个定轨弧段的RMSE曲线。图4中，DHIF-a，DHIF-b为拓扑a，b下的分布式滤波，IF为测站单站信息滤波，CIF为4个测站集中式信息滤波。表4和表5列出了集中式滤波和分布式滤波RMSE的均值。

结果表明，使用该分布式算法，各测站的定轨精度均优于单站信息滤波定轨精度，且在不同的通信拓扑下均逼近于集中式信息滤波定轨精度。该算法中，测站每次状态更新前与相邻测站进行信息传递，单个测站融合了其他测站的信息，因此定轨结果优于单站信息滤波。同时，虽然不像集中式算法计算中心可以直接获取所有测站的测量值，但通信拓扑的连通性使得各个测站间接获得了其他测量信息，得到了与集中式定轨精度相当的结果。

4.2 测站测量方式不同

设4个测站采用不同测量设备，见表6，不同测量设备的观测量见表7。

R1=R4=diag(202,0.022,0.022)，

R2=diag(0.0022,0.0022)，

R3=diag(102,0.052,0.022,0.022)。

图5为两种拓扑下，各测站80 s内分布式滤波定轨位置和速度的RMSE曲线。

表5 各测站集中式与分布式滤波定轨速度平均RMSETable 5 Velocity average RMSE of DHIF and CIF for each station

表6 各测站测量设备Table 6 Measure equipment of each observation station

80 s时，定轨位置误差在10 m以内，速度误差在0.1 m/s以内。图6为各测站不同拓扑结构下定轨比较。

表7 不同设备的观测量Table 7 The observation type of different equipment

图5 各测站80 s内分布式滤波定轨的位置、速度RMSE曲线Fig.5 80 s position and velocity RMSE of DHIF for each station

图6(a)、(b)中，测站S1和S2在拓扑b下收敛速度快于拓扑a。这是因为拓扑b的结构使得测站S1和S2的1和2类信息能更快的传输至整个网络中。S1和S2在起始阶段，收到更多的是1和2类信息。而在拓扑a中，起始阶段S1和S2很快收到了来自S3和S4的第3类信息，特别对于S2，开始几秒只接收到了第3类信息，因此收敛速度最慢。

图6(c)、(d)中，结果则相反，测站S3，S4在拓扑a下收敛速度快于拓扑b。因为虽然通信拓扑a的信息交换速度比拓扑b慢，但拓扑a中的S3和S4在前几秒收到的更多是第1和2类信息，尤其S4，开始只收到了第1和2类信息，因此很快收敛。而拓扑b中的S3和S4开始时收到的1和2类信息少，第3类信息多，因此收敛速度相对较慢。

图5也验证了上述结论，图5 (a)中，根据拓扑a中信息的传递方向，收敛速度最快的为S4，最慢的为S2。图5(b)中，因为拓扑b使得各类信息在整个网络中可以快速传递，四个测站的收敛速度基本接近。

图6 测站在不同拓扑下80 s内分布式实时定轨位置、速度RMSE曲线对比Fig.6 Comparison of 80 s position and velocity RMSE for each station in different topologies

5 结 论

针对我国航天器测控系统，本文结合DHIF设计了一种分布式实时轨道确定算法。该算法使得单站分布式定轨精度优于单站滤波定轨精度并逼近集中式滤波方法，可以扩展至天地基多种传感器联合观测的分布式实时定轨中。该算法要求各测站的量测数据在同一时刻，在实际工程中，由于不同传感器的量测模型和数据修正模型不同，存在不同测站修正后的量测数据不在同一时刻的问题，因此如何进行多传感器异步数据的分布式融合，是需进一步研究的课题。