CN110146084B

CN110146084B - 面向卫星故障的多星编队系统分布式协同导航滤波方法

Info

Publication number: CN110146084B
Application number: CN201910437283.9A
Authority: CN
Inventors: 李文玲; 刘杨; 陈达
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2019-05-24
Filing date: 2019-05-24
Publication date: 2020-12-04
Anticipated expiration: 2039-05-24
Also published as: CN110146084A

Abstract

本发明公开了一种面向卫星故障的多星编队系统分布式协同导航滤波方法，采用无迹卡尔曼滤波‑扩展卡尔曼滤波算法作为基本的滤波算法，通过将上一时刻的计算结果作为反馈值，代入到当前时刻进行节点的计算结果验证，判断节点计算是否正常，并以此为依据，对一致性算法的权值进行自适应计算，最后使用协方差交集算法进行数据融合；该方法可以提升多星编队系统的容错性，能够有效应对多星编队系统中某颗环绕卫星出现量测故障或计算结果丢失、不准确的情况，保证最终滤波结果的精确度在上述情况出现时无明显偏差，对于环绕卫星测距误差较大，节点计算结果丢失、不准确的情况，都具备一定的鲁棒性。

Description

面向卫星故障的多星编队系统分布式协同导航滤波方法

技术领域

本发明涉及卫星导航技术领域，尤其涉及一种面向卫星故障的多星编队系统分布式协同导航滤波方法。

背景技术

随着卫星功能的日益庞大复杂，仅仅利用单颗卫星在轨工作的成本和风险逐渐增加。受到计算机网络技术发展的影响，多星编队系统的概念在20世纪90年代被提出。多星编队系统是按照一定要求将两颗或两颗以上的卫星分布在一种或多种轨道上，各个卫星间协同工作、相互联系，以完成各项飞行任务的卫星系统。

目前，多星编队系统分布式导航滤波算法，大都采用卡尔曼滤波算法及其相关算法(例如无迹卡尔曼滤波算法)作为基础的滤波算法。标准的卡尔曼滤波，每次迭代的末值会作为下一次迭代的初值状态，经过一步状态预测得到中值，并进行量测预测。最终根据状态预测和量测预测的差值，乘以适当的修正系数之后，反馈回中值，得到此刻的最终结果，并开始下一轮迭代。

在多星编队系统中，还可以利用一致性算法，一方面能够分散计算压力，另一方面可以得到更为准确的数据。一致性算法是一类施加于分布式网络各节点上具有相同形式且能使节点特定信息状态趋于一致的算法。基于一致性算法的卡尔曼滤波框架对上述中值加入了额外的处理过程，即在得到自身节点的中值后，利用邻居节点的中值，进行一致性算法求解，将求解的稳态值替代原节点中值进行后续计算。

在多星编队系统中，对于各个节点得到的目标卫星的状态信息，还可以利用协方差交集(CI)算法进行数据融合。当局部数据之间的相关程度准确可知时，可利用卡尔曼融合得到最优融合，而当局部数据之间的相关程度不可知时，协方差交集算法则能够不失一般性地代表融合结果。

综上，现有的多星编队系统分布式协同导航滤波方法，大都将关注点放在了计算精确度的提升上，而忽略了多星编队系统本身容错性较高的特点，未考虑到多星编队系统中部分节点工作在非正常状态时所带来的影响。

发明内容

有鉴于此，本发明提供了一种面向卫星故障的多星编队系统分布式协同导航滤波方法，用以使多星编队系统中某颗环绕卫星出现量测故障或计算结果丢失、不准确的情况时仍能保持多星编队系统得到较为准确的结果。

因此，本发明提供了一种面向卫星故障的多星编队系统分布式协同导航滤波方法，包括如下步骤：

S1：对各节点的滤波器进行初始化；其中，目标卫星的多颗环绕卫星构成一个节点；

S2：利用无迹变换进行一步状态预测，得到第一状态预测值；

S3：将得到的第一状态预测值进行分布式通信，通过自适应一致性算法达成一致状态，作为各滤波器后续处理的预测值；

S4：按照扩展卡尔曼滤波的方式进行量测更新，以卡尔曼滤波的标准形式进行后续的状态更新；

S5：利用上一时刻的滤波结果按照扩展卡尔曼滤波的方式进行一步状态预测，得到第二状态预测值，计算第二状态预测值与当前时刻各节点的状态估计值的误差比例，更新各节点的工作状态；

S6：利用协方差交集算法得到当前时刻的融合结果，作为当前时刻的滤波结果；

S7：预定时间间隔后，返回步骤S2，执行步骤S2～步骤S6，进行下一时刻的滤波计算，直至导航结束。

在一种可能的实现方式中，在本发明提供的上述多星编队系统分布式协同导航滤波方法中，步骤S3，将得到的第一状态预测值进行分布式通信，通过自适应一致性算法达成一致状态，作为各滤波器后续处理的预测值，具体包括：

S31：设定变量t_ij和变量t_ij的阈值范围；其中，t_ij表示环绕卫星i记录环绕卫星j正常运转与非正常运转的差值，t_ij的初始值为0；

S32：设定常数ρ₁、ρ₂；其中，ρ₁表示节点的初始权重，ρ₂表示节点权重趋向0或1的速度，0<ρ₁,ρ₂<1；

S33：判断t_ij是否小于零；若是，则执行步骤S34；若否，则执行步骤S35；

S34：按照下式计算自适应权值：

S35：按照下式计算自适应权值：

S36：对环绕卫星i和与环绕卫星i进行通信的环绕卫星j的集合进行权值的归一化：

S37：对各节点进行一致性计算；

步骤S5，利用上一时刻的滤波结果按照扩展卡尔曼滤波的方式进行一步状态预测，得到第二状态预测值，计算第二状态预测值与当前时刻各节点的状态估计值的误差比例，更新各节点的工作状态，具体包括：

S51：利用上一时刻的滤波结果按照扩展卡尔曼滤波的方式进行一步状态预测，得到第二状态预测值，计算第二状态预测值与当前时刻各节点的状态估计值的误差比例，计算方式如下：

其中，R_i表示第i个状态分量的误差比例，x_sub,i为节点输出的状态分量，x_pre,i为利用反馈值进行一步状态预测的状态分量；

S52：对R_i设定一个阈值，判断R_i是否超过阈值；若是，则执行步骤S53；若否，则执行步骤S54；

S53：确定卫星j对应的节点出现故障，变量t_ij减1；

S54：确定卫星j对应的节点正常，变量t_ij加1。

S55：判断变量t_ij是否超出变量t_ij的阈值范围；若是，则执行步骤S56；若否，则执行步骤S57；

S56：将变量t_ij限定在变量t_ij的阈值范围内，返回步骤S33；

S57：返回步骤S33。

本发明提供的上述多星编队系统分布式协同导航滤波方法，采用无迹卡尔曼滤波-扩展卡尔曼滤波算法作为基本的滤波算法，通过将上一时刻的计算结果作为反馈值，代入到当前时刻进行节点的计算结果验证，判断节点计算是否正常，并以此为依据，对一致性算法的权值进行自适应计算，最后使用协方差交集算法进行数据融合；该方法可以提升多星编队系统的容错性，能够有效应对多星编队系统中某颗环绕卫星出现量测故障或计算结果丢失、不准确的情况，保证最终滤波结果的精确度在上述情况出现时无明显偏差，对于环绕卫星测距误差较大，节点计算结果丢失、不准确的情况，都具备一定的鲁棒性。

附图说明

图1为本发明实施例提供的面向卫星故障的多星编队系统分布式协同导航滤波方法的流程示意图；

图2为本发明实施例提供的面向卫星故障的多星编队系统分布式协同导航滤波方法的流程图之一；

图3为本发明实施例提供的面向卫星故障的多星编队系统分布式协同导航滤波方法的流程图之二；

图4为本发明实例1中多星编队系统的通信拓扑结构图；

图5为本发明实例1中多星编队系统在仿真时间750s时的仿真图；

图6a为本发明实例1中本仿真的位置误差图；

图6b为本发明实例1中本仿真的速度误差图；

图7a为比较本仿真、扩展卡尔曼滤波算法与无迹卡尔曼滤波算法的位置误差图；

图7b为比较本仿真、扩展卡尔曼滤波算法与无迹卡尔曼滤波算法的速度误差图；

图8a为一致性扩展卡尔曼滤波算法的位置误差图；

图8b为一致性扩展卡尔曼滤波算法的速度误差图；

图9a为本发明实例2中本仿真在测距故障下的位置误差图；

图9b为本发明实例2中本仿真在测距故障下的速度误差图；

图10a为本发明实例2中一致性扩展卡尔曼滤波算法在测距故障下的位置误差图；

图10b为本发明实例2中一致性扩展卡尔曼滤波算法在测距故障下的速度误差图；

图11a为本发明实例2中无迹卡尔曼滤波结合协方差交集算法在测距故障下的位置误差图；

图11b为本发明实例2中无迹卡尔曼滤波结合协方差交集算法在测距故障下的速度误差图；

图12a为本发明实例3中本仿真在通信故障下的位置误差图；

图12b为本发明实例3中本仿真在通信故障下的速度误差图；

图13a为本发明实例3中一致性扩展卡尔曼滤波算法在通信故障下的位置误差图；

图13b为本发明实例3中一致性扩展卡尔曼滤波算法在通信故障下的速度误差图；

图14a为本发明实例3中无迹卡尔曼滤波结合协方差交集算法在通信故障下的位置误差图；

图14b为本发明实例3中无迹卡尔曼滤波结合协方差交集算法在通信故障下的速度误差图。

具体实施方式

下面将结合本申请实施方式中的附图，对本申请实施方式中的技术方案进行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限制本申请。

本发明实施例提供的一种面向卫星故障的多星编队系统分布式协同导航滤波方法，如图1和图2所示，图1和图2分别为多星编队系统分布式协同导航滤波方法的流程示意图和流程图，包括如下步骤：

本发明实施例提供的上述多星编队系统分布式协同导航滤波方法，采用无迹卡尔曼滤波(UKF)-扩展卡尔曼滤波(EKF)算法作为基本的滤波算法，通过将上一时刻的计算结果作为反馈值，代入到当前时刻进行节点的计算结果验证，判断节点计算是否正常，并以此为依据，对一致性算法的权值进行自适应计算，最后使用协方差交集算法进行数据融合；该方法可以提升多星编队系统的容错性，能够有效应对多星编队系统中某颗环绕卫星出现量测故障或计算结果丢失、不准确的情况，保证最终滤波结果的精确度在上述情况出现时无明显偏差，对于环绕卫星测距误差较大，节点计算结果丢失、不准确的情况，都具备一定的鲁棒性。

本发明实施例提供的上述多星编队系统分布式协同导航滤波方法，步骤S2采用无迹卡尔曼滤波的方案对各节点进行一步状态预测，可以在状态方程非线性强的情况下仍能得到精确的预测结果；步骤S3中的自适应一致性算法，一方面可以提升算法的计算精确度，另一方面是处理多星编队系统故障的主要环节，当有环绕卫星出现计算不准确、结果丢失的情况时，通过一致性算法的权值自适应，能够及时排除问题环绕卫星，增大正常环绕卫星计算结果所占权值；步骤S4利用扩展卡尔曼滤波进行非线性较弱的量测预测，在损失计算精度较小的情况下，可以保证计算速度，使多星编队系统在突发情况下具有更快速的反应能力；步骤S5利用上一时刻的计算结果作为当前时刻的反馈值，判断各节点是否出现问题，给下一时刻的权值自适应提供依据，同时，出现问题的节点会将反馈值的一步预测作为故障状态下的输出值，可以为多星编队系统替换或排除问题卫星提供反应时间，从而可以保证最终滤波结果不会有明显的波动；步骤S6对各节点输出结果进行融合，得到更加精确可靠的结果。

需要说明的是，在本发明实施例提供的上述多星编队系统分布式协同导航滤波方法中，步骤S7，预定时间间隔后，返回步骤S2，执行步骤S2～步骤S6，进行下一时刻的滤波计算，直至导航结束，该预定时间间隔可以为1s，即可以在1s后，返回步骤S2，执行步骤S2～步骤S6，进行下一时刻的滤波计算，也就是说，每秒都输出一个滤波结果。具体地，可以根据实际需要设定预定时间间隔，在此不做限定。

在具体实施时，在本发明实施例提供的上述多星编队系统分布式协同导航滤波方法中，步骤S3，将得到的第一状态预测值进行分布式通信，通过自适应一致性算法达成一致状态，作为各滤波器后续处理的预测值，如图3所示，具体可以包括如下步骤：

S34：按照下式计算自适应权值：

S35：按照下式计算自适应权值：

S37：对各节点进行一致性计算；

步骤S5，利用上一时刻的滤波结果按照扩展卡尔曼滤波的方式进行一步状态预测，得到第二状态预测值，计算第二状态预测值与当前时刻各节点的状态估计值的误差比例，更新各节点的工作状态，如图3所示，具体可以包括如下步骤：

S53：确定卫星j对应的节点出现故障，变量t_ij减1；

S54：确定卫星j对应的节点正常，变量t_ij加1；

S55：判断变量t_ij是否超出变量t_ij的阈值范围；若是，则执行步骤S56；若否，则返回步骤S33；

S56：将变量t_ij限定在变量t_ij的阈值范围内；返回步骤S33。

下面通过一个具体的实例对本发明实施例提供的上述多星编队系统分布式协同导航滤波方法的具体实施过程进行详细说明。

实例1：部分节点不可观的情况。

多星编队系统的分布式协同导航滤波方法，是在各环绕卫星状态(空间位置及速度)已知的情况下，协同求取目标卫星状态的方法。环绕卫星能够对目标卫星量测，得到部分或全部的状态量测值。各环绕卫星按照一定的拓扑结构形成通信链路，进行数据交互，组成多个计算节点。节点是能够进行算法计算的基本的单位，是一组具备星间链路的环绕卫星的集合。该集合中包含一颗环绕卫星，该环绕卫星与节点中的其他环绕卫星都有信息通道，以获取其他环绕卫星的量测数据，并在该卫星上进行计算。节点是虚拟出来的整体，同样的一组环绕卫星，如果负责计算的卫星不同，可看做构成了不同的节点。

建立一个具有四颗环绕卫星和一颗目标卫星的仿真实例，采用星间测距的方式建立量测方程，其通信拓扑结构如图4所示。每三颗环绕卫星可组成一个节点，例如，环绕卫星1、2、4可以构成一个节点进行计算，由于环绕卫星1可同时与环绕卫星2、4进行通信，因此，此节点的最终计算结果由环绕卫星1输出，称该节点的对应卫星为环绕卫星1。同理，环绕卫星1、2、3，环绕卫星2、3、4，环绕卫星1、3、4，皆可构成一个节点，其各自对应的环绕卫星分别为2、3、4。该多星编队系统共可组成4个节点进行滤波计算，最终将各个节点的计算结果进行融合，得到最后的滤波结果。

仿真中，考虑到地球J₂项非球形摄动是卫星所收到的最主要的摄动项，在此仿真实例中，使用J2000地心惯性坐标系，并考虑J₂项摄动影响。各卫星的初始位置如表1所示。

表1卫星初始位置

表1中，a为半长轴，e为偏心率，i_orb为卫星轨道倾角，ω为升交点赤经，Ω为近地点幅角，f是真近点角。

此时，卫星的状态方程如下：

其中，x,y,z分别为J2000地心惯性坐标系下目标卫星的位置；R_e为地球半径，取值为6378137m；μ为地球引力常数，取值为398600.44×10⁹m³/s²；J₂为带谐项系数，取其值为4.8416685×10^-4；f_x、f_y、f_z分别为真近点角在J2000地心惯性坐标系下的分量。

采用星间测距的方式对目标卫星进行观测时，观测模型(即量测方程)如下：

采用EKF的思想，对量测方程进行线性化：

至此，获得卡尔曼滤波所需的状态方程(式(1))和量测方程(式(3))。

方正采样时间为1s，对模型的计算时间为7501s(对应仿真中的0时刻)到10501s。在仿真时间750s(对应模型时间8250s)附近，多星编队系统中部分环绕卫星近似在一条直线上，存在节点状态不可观的情况，如图5所示。

设初始滤波值已知，为目标卫星0时刻状态，卡尔曼滤波的初始状态估计值：

状态估计值误差协方差矩阵：

P_x,0＝diag(10⁶ 10⁶ 10⁶ 10⁶ 10⁶ 10⁶)

过程噪声协方差矩阵：

Q₀＝diag(50 50 50 10 10 10)

测量噪声协方差矩阵：

R₀＝diag(10⁴ 10⁴ 10⁴)

量测误差服从高斯分布，均值为0，标准差在1m～5m之间随机生成，运行50次Monte-Carlo仿真。

利用上一时刻(第一次计算时，上一时刻为0时刻)的状态估计值，使用无迹(UT)变换计算当前时刻的一步状态预测，计算公式如下：

本仿真中，取α＝10^-3，β＝2，κ＝0。

求出状态预测值x_tmp与状态协方差预测值p_x,tmp。

对求出的协方差预测值加上过程噪声协方差矩阵：

p′_x,tmp＝p_x,tmp+Q_k (6)

之后再求取后续一致性计算所需的自适应加权系数，求取过程如下：

对于每个环绕卫星，以环绕卫星i为例，初始化卫星状态记录变量t_ij＝0。N_i为环绕卫星i和能够与环绕卫星i进行通信的环绕卫星的集合。对N_i中包含的各环绕卫星，其一致性算法的权值按如下公式计算：

其中，ρ₁和ρ₂为自定义值，0<ρ₁,ρ₂<1，ρ₁表示初始的权重因子，本仿真中ρ₁取值为0.5，ρ₂表示权重趋向0或1的速度，本仿真中ρ₂取值为0.5。

利用下式进行权值的归一化：

求得自适应权值之后，对状态估计值按照下式进行一致性计算：

协方差预测值一致性计算公式为：

p′_j,x,tmp表示卫星j的协方差预测值。

采用EKF的思想，利用式(3)求取量测预测，利用下式求取量测预测的协方差阵：

完成后，利用卡尔曼滤波的状态更新方程，即下式的向量形式，完成节点的状态求取：

对于节点的输出，在进行数据融合前，先对上一时刻的滤波结果进行EKF的一步状态预测，得到预测值x_pre。之后按照下式计算节点输出的各状态分量的误差比例：

对各状态分量的误差比例需要设定一个阈值，若误差比例大于此阈值则说明环绕卫星故障。仿真中取各状态分量的阈值为0.1，即，若R_i>0.1，则环绕卫星故障，替换环绕卫星输出结果为x_pre。

更新卫星状态记录变量t_ij，若环绕卫星故障，则t_ij减1，若环绕卫星正常，则t_ij加1。当t_ij超出指定阈值范围时，t_ij不再增加或减小。仿真中取-20≤t_ij≤20。即，若t_ij大于20，则t_ij不再增大，一直保持20；若t_ij小于-20，则t_ij不再减小，一直保持-20。

得到各节点输出后，使用协方差交集算法进行最终的数据融合，协方差交集算法如下式：

权值更新公式如下式：

融合后的状态值即为最终滤波结果。

进行新一时刻的计算时，将最终滤波结果作为新的初始状态，开始新一轮计算。

下面对本仿真的输出结果进行说明。取50s为间隔绘制仿真输出结果如图6a和图6b所示，图6a为本仿真中的位置误差，图6b为本仿真中的速度误差，如图6a和图6b所示，在仿真750s附近，环绕卫星的状态滤波依旧在收敛的过程中，位置误差和速度误差皆无明显波动。

下面比较EKF算法、UKF算法、一致性EKF算法和本仿真的计算精确度(速度误差、位置误差)。

对比使用EKF算法和UKF算法，EKF算法和UKF算法以使用环绕卫星1、2、4组成的节点进行计算为例，仿真结果如图7a和图7b所示。图7a为比较本仿真、EKF算法和UKF算法的位置误差，图7b为比较本仿真、EKF算法和UKF算法的速度误差，从图7a和图7b可以看出，本仿真在计算精确度上优于EKF算法和UKF算法。

再对EKF算法进行一致性改造，权值用最大度加权计算，其结果如图8a和图8b所示，图8a为一致性EKF算法的位置误差，图8b为一致性EKF算法的速度误差。如图8a和图8b所示，一致性EKF算法在750s附近仍然会因为部分节点不可观的情况，给位置误差带来一个峰值，且收敛速度慢，整体波动明显。

取2500s后，滤波结果稳定后的误差值进行统计，分别计算EKF算法、UKF算法、一致性EKF算法与本仿真的位置误差与速度误差的RMS(Root Mean Square)，结果为表2所示。

表2四种方案的RMS

方案	EKF	UKF	一致性EKF	本仿真
					位置RMS/m	5.0140	4.5772	5.0591	2.8923
速度RMS/(m/s)	2.6896	2.5823	2.2808	1.7598

综上可知，多星编队系统出现节点不可观的情况下，采用分布式计算能够有效消除这种情况所带来的影响，且本仿真在稳定后的位置误差和速度误差皆为最小值。

下面通过一个具体的实例比较本仿真、一致性EKF算法和UKF结合CI算法在测距故障的情况下的结果。

实例2：

仍然采用本仿真模型和各项数值，算法流程不变。选取模型9001s～11500s进行仿真，模型的9001s对应本仿真的1s。此时间段内，节点均可观。令环绕卫星1在2001s到2050s，这50s中距离量测出现故障，返回值为一较大数值，仿真中设为100000000m。比较本仿真与一致性EKF算法、UKF结合CI算法的结果，图9a和图9b为测距故障下本仿真的位置误差和速度误差，图10a和图10b为测距故障下一致性EKF算法的位置误差和速度误差，图11a和图11b为测距故障下UKF结合CI算法的位置误差和速度误差。从仿真结果中可以看出，本仿真对于环绕卫星的量测出现较大偏差时能够较好地对自身进行调整，状态估计值基本不受影响，而使用一致性EKF算法，状态估计值在量测出现误差之后全部跑偏，使用UKF结合CI算法，状态估计值在量测出现误差的时刻会出现较大偏差，一段时间后迅速缩减，但仍比原来的误差要大。

下面通过一个具体的实例比较本仿真、一致性EKF算法和UKF结合CI算法在节点结果丢失的情况下的结果。

实例3：

仍然采用本仿真模型和各项数值，算法流程不变。继续在时间段9001s～11500s内进行卫星通信故障仿真。将9000s作为仿真的起始时间点，环绕卫星1在2501s～2550s时间段内，结算结果丢失，结果通信值恒为0，比较本仿真、一致性EKF算法和UKF结合CI算法的结果。图12a和图12b为通信故障下本仿真的位置误差和速度误差，图13a和图13b为通信故障下一致性EKF算法的位置误差和速度误差，图14a和图14b为通信故障下UKF结合CI算法的位置误差和速度误差。从仿真结果中可以看出，在丢失了环绕卫星1的计算结果后，一致性EKF算法与UKF结合CI算法都会在丢失时刻产生巨大误差，在重新收悉到环绕卫星1的节点后，才回到正常工作的状态；而本仿真在丢失环绕卫星1的节点后，位置误差和速度误差都有一个微小的增大，此时丢失的卫星节点的权值迅速下降，待其输出正常后，权值逐渐增大，回到初始的权值配比中，整个过程保证了滤波结果的波动较小，依旧能够得到精确度同卫星无故障时相同的结果。

本发明实施例提供的上述多星编队系统分布式协同导航滤波方法，采用无迹卡尔曼滤波-扩展卡尔曼滤波算法作为基本的滤波算法，通过将上一时刻的计算结果作为反馈值，代入到当前时刻进行节点的计算结果验证，判断节点计算是否正常，并以此为依据，对一致性算法的权值进行自适应计算，最后使用协方差交集算法进行数据融合；该方法可以提升多星编队系统的容错性，能够有效应对多星编队系统中某颗环绕卫星出现量测故障或计算结果丢失、不准确的情况，保证最终滤波结果的精确度在上述情况出现时无明显偏差，对于环绕卫星测距误差较大，节点计算结果丢失、不准确的情况，都具备一定的鲁棒性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。