CN106482728B - 基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法 - Google Patents

Abstract

本发明提供一种基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，利用运动学方程推导，构建通信保障航天器和通信故障航天器的相对运动方程模型。利用具有测量距离远、测量精度高，可全天候工作等优点的无线电设备为测量部件，测出两个航天器的相对测量信息；采用一种基于最大互相关熵准则无迹卡尔曼滤波，进行相对位置和相对速度信息的估计。本发明在具有冲激噪声的实际系统中，能有效提高相对位置和相对速度的估计精度，进而提高通信保障航天器天线指向精度，保障通信链路和质量，防止通信中断。

Description

基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相 对状态估计方法

技术领域

本发明属于通信领域，涉及一种通信保障航天器的相对状态估计方法。

背景技术

随着航天技术的快速发展，空间在轨服务技术已经成为在复杂的空间环境中亟待解决的重要问题。在轨服务技术能实现在轨维护、在轨装配和在轨修理等空间操控任务，引起了越来越多的关注，拥有广泛的市场应用前景。其中在轨服务航天器在空间网络中具有非常重要的功能，而通信保障航天器作为一种在轨服务航天器，其主要任务就是利用自身设备为通信故障航天器建立星地通信链路。在一般的情况下，通信保障航天器通过机动飞行到达通信功能受限的航天器通信覆盖区域，与其构成飞行编队，并为其提供与地面的通信中继服务。在这个过程中，无论是在轨服务还是航天器编队飞行任务的实现，都依赖于航天器之间相对位置、相对速度的确定精度。因此，相对位置、相对速度的确定是卫星编队飞行和通信有效保障的关键技术之一。

常用的编队卫星相对状态确定方法包括GPS测量法、激光测量法、红外测量法、可见光测量法和无线电测量法等。但GPS的S/A政策对我国有一定限制，因此其存在不安全性；激光测量抗干扰能力强，测量精度高，但其需要消耗的功率大，数据率低，波束狭窄，需要其他测量系统进行引导；红外测量仅仅能进行角度测量，不能单独实现卫星的相对状态测量；可见光测量仅仅适用于近距离不超过几十米的测量，且随着距离不断增加，测量精度不断变差；而星间无线电测量，具有作用距离远、测量覆盖高、测量精度高、实时性强，可全天候工作等优点，且具有信息通信功能，满足相对测量的要求。

选择采用无线电测量法进行通信保障航天器和故障航天器的相对位置和相对速度测量时，由于外部环境对系统存在扰动以及测量系统自身存在一定的噪声影响，同时，状态方程和测量方程均为非线性的，因此需要涉及非线性状态估计方法，而在非线性估计领域，扩展卡尔曼滤波(EKF)是广泛应用的方法。扩展卡尔曼滤波是对非线性方程进行一阶线性化，提供一阶逼近结果，依赖于局部非线性程度，因此可能会为后验均值和协方差带来较大的误差，甚至导致滤波发散。而无迹卡尔曼滤波(UKF)是一种更精确的估计方法，它继承了卡尔曼滤波框架，是一种切实可行便于应用的滤波方法。在无迹卡尔曼滤波中，非线性函数的概率分布是通过一个特定的最小样本点集(sigma点)来逼近的。这个最小样本点集是通过无迹变化(UT)得到的，这些样本点包含了均值和协方差信息，在实际的非线性方程传播时，得到的均值和协方差能达到三阶精度。与扩展卡尔曼滤波相比，无迹卡尔曼滤波不需要计算雅可比矩阵，只需要通过无迹变化获得均值和协方差，特别是对于比较复杂的系统，计算简便很多。

但无迹卡尔曼滤波只是对于高斯分布的噪声具有较好的性能，对于非高斯噪声，性能下降很快，鲁棒性差。近年来，信息论学习在许多领域受到了越来越多的关注，基于信息论的算法已经显示了很多的优越性，特别是互相关熵(correntropy)作为一种新的相似度度量，被用作目标函数来解决了许多问题，其对应的最大互相关熵准则具有坚实的理论基础，为非高斯噪声的处理提供了新的、鲁棒的解决方案。

为降低无线电测量设备误差的影响，特别是非高斯噪声的影响，因此，亟待提出一种新的无迹卡尔曼滤波方法对测量结果进行修正，提高相对位置和速度估计精度，从而为通信保障航天器的有效工作提供基础。

发明内容

本发明的目的在于提供一种基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法。

为达到上述目的，本发明采用了以下技术方案：

(1)利用安装在通信故障航天器上的无线电设备测出通信故障航天器和通信保障航天器的相对测量信息；

(2)通过运动学方程获取在地心惯性坐标系下通信故障航天器的绝对运动状态变量，并求得在希尔坐标系中通信保障航天器和通信故障航天器的相对运动方程；

(3)以互相关熵作为相似度度量并将最大互相关熵准则与无迹卡尔曼滤波结合，形成基于最大互相关熵准则的无迹卡尔曼滤波方法；

(4)将相对测量信息作为测量值，采用基于最大互相关熵准则的无迹卡尔曼滤波，计算出通信故障航天器和通信保障航天器之间的相对状态的估计值。

所述相对测量信息包括相对距离ρ、方位角θ和仰角φ：

其中，x、y和z分别为以通信故障航天器为中心的希尔坐标系的X_c、Y_c和Z_c轴的坐标，坐标原点固定于通信故障航天器所在轨道位置，X_c轴沿地心指向通信故障航天器质心的矢量方向，Z_c轴与通信故障航天器轨道角动量正方向一致，Y_c轴垂直于轨道平面并与X_c和Z_c轴构成右手正交坐标系。

所述步骤(2)具体包括以下步骤：

(2-1)在地心惯性坐标系下通信故障航天器的绝对运动状态变量的运动方程为：

其中，r_c和ω为地心惯性坐标系下通信故障航天器的绝对角加速度、绝对位置和绝对角速度，μ为地心引力常数；通过构造绝对状态向量为用函数等价表示上式，然后进行离散化，得到离散形式的运动方程；

(2-2)在以通信故障航天器为中心的希尔坐标系下通信保障航天器的相对运动方程为

其中，和分别为X_c、Y_c和Z_c轴上相应的加速度；通过构造相对状态向量为用函数等价表示上式，然后进行离散化，得到离散形式的运动方程。

所述基于最大互相关熵准则的无迹卡尔曼滤波包括以下步骤：

(3-1)根据无迹变化计算出先验的相对状态估计值和协方差矩阵P(k|k-1)，并采用Cholesky分解获得T(k)

其中，R(k)为测量噪声协方差；

(3-2)根据步骤(3-1)构造以下式子：

D(k)＝g(k,x(k))+e(k)

其中， x(k)为真实状态值，y(k)为测量值，h(·)为测量函数，r(k)为测量噪声；

(3-3)由步骤(3-2)得到基于互相关熵的代价函数

其中，d_i(k)表示向量D(k)的第i个元素，g_i(k,x(k))表示g(k,x(k))的第i个元素，L＝n+m其中n为状态的维数，m为测量值的维数；

(3-4)在最大互相关熵准则下，状态量的优化值由下式得到

其中，e_i(k)＝d_i(k)-g_i(k,x(k))；

(3-5)更新的协方差表示为

其中

取则得到

(3-6)将代替R(k)，根据无迹变化就可以计算出后验的相对状态估计值和协方差矩阵P(k|k)。

所述步骤(4)还包括以下步骤：在利用基于最大互相关熵准则的无迹卡尔曼滤波进行相对状态估计前，通过使用无迹卡尔曼滤波加快估计初期的收敛速度。

所述初期是指时间标签≤100。

本发明与现有技术相比，优点在于：

(1)本发明方法采用了无线电测量法，具有作用距离远、测量覆盖高、测量精度高、实时性强，可全天候工作等优点，为方法的正确实施奠定了基础。

(2)本发明方法中，将最大互相关熵准则和无迹卡尔曼滤波方法结合，增强了鲁棒性，提高了在实际系统存在冲激噪声的状况下的性能。

附图说明

图1为本发明的实时场景示意图；

图2为本发明的测量系统示意图；

图3为本发明方法的流程框图；

图4为本发明算法与其它算法的相对位置估计误差比较；MCUKF表示基于最大互相关熵准则的无迹卡尔曼滤波，σ表示基于最大互相关熵准则的无迹卡尔曼滤波中的核宽度；

图5为本发明算法与其它算法的相对速度估计误差比较。

具体实施方式

下面结合附图和实施例对本发明做进一步详细说明。

首先说明需要用到的两个坐标系。

地心惯性坐标系：坐标原点位于地球质心，X轴在赤道平面内指向春分点，Z轴垂直赤道面指向地球北极，Y轴由右手法则确定。其中，常用的地心惯性坐标系是J2000坐标系，它的坐标原点在地球质心，参考平面是J2000.0平赤道面，Z轴向北指向平赤道面北极，X轴指向J2000.0平春分点，Y轴与X轴和Y轴组成直角右手系；

希尔坐标系：即当地垂直水平坐标系，坐标原点O_c固定于航天器所在轨道位置，X_c轴沿地心指向航天器质心的矢量方向，Z_c轴与航天器轨道角动量正方向一致，Y_c轴垂直于轨道平面并与前两轴构成右手正交坐标系，参见图1以及图2。

本发明是一种基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，具体如下：

步骤1：在通信故障航天器上安装无线电设备，利用无线电设备测出通信故障航天器和保障航天器的相对测量信息，包括相对距离ρ，方位角θ和仰角φ。

无线电测量法具有作用距离远、测量覆盖高、测量精度高、实时性强，可全天候工作等优点，测量方程如下：

式中，x、y和z分别为以通信故障航天器为中心的希尔坐标系下的X_c、Y_c和Z_c轴的坐标。

对式(1)考虑噪声扰动项r_ρ，r_θ，r_φ，其协方差矩阵为R，即

其中，σ_ρ、σ_θ和σ_φ分别为测量系统中相对距离ρ、方位角θ和仰角φ的随机扰动的标准差。

步骤2：通过运动方程，可以获取在地心惯性坐标系下通信故障航天器的绝对运动状态变量，并求得在希尔坐标系中通信保障航天器和通信故障航天器的相对运动方程。

在地心惯性坐标系下通信故障航天器的绝对运动状态变量的运动方程为：

其中，r_c和ω为地心惯性坐标系下通信故障航天器的绝对角加速度、绝对位置和绝对角速度，μ为地心引力常数，和分别为绝对加速度和绝对速度。

构造绝对状态向量为用函数f_jd等价表示式(3)，即

表示绝对状态向量r的一阶导数，下标jd表示绝对的缩写；

对式(4)离散化得到相应的离散形式运动方程为

r_k+1＝F_jd(r_k) (5)

在以通信故障航天器为中心的希尔坐标系下通信保障航天器的相对运动方程为

其中，x、y、z分别为希尔坐标系下的X_c、Y_c、Z_c轴的坐标，和分别为相对运动模型中X_c、Y_c和Z_c轴上相应的加速度，和分别表示X_c和Y_c轴上相应的速度，ω、和r_c为地心惯性坐标系下通信故障航天器的绝对角速度、绝对角加速度和绝对位置，μ为地心引力常数；ω，r_c为利用步骤(2-1)得到的值；

对式(4)考虑噪声扰动项其协方差矩阵为Q，即

其中，和分别为相对运动模型中X_c、Y_c和Z_c轴上相应的加速度随机扰动的标准差。

构造相对状态向量为以函数f_xd等价表示式(6)，即有

表示相对状态向量X的一阶导数，下标xd表示相对的缩写；

对式(8)离散化得到相应的离散形式运动方程为

X_k+1＝F_xd(X_k) (9)

步骤3：基于最大互相关熵准则无迹卡尔曼滤波方法采用一种新的相似度度量——互相关熵(correntropy)及其对应的鲁棒准则，即最大互相关熵准则，来改进传统的无迹卡尔曼滤波在实际存在的冲激噪声环境中的性能。具体如下：

步骤3-1

该方法首先根据无迹变化(Unscented Transform)计算出先验的相对状态估计值和协方差矩阵P(k|k-1)，并采用Cholesky分解获得T(k)，其中

R(k)为测量噪声协方差。

步骤3-2

根据步骤3-1构造以下式子：

D(k)＝g(k,x(k))+e(k)

其中 x(k)为真实状态值，y(k)为测量值，h(·)为测量函数，r(k)为测量噪声。

步骤3-3

由步骤3-2，我们可以得到基于一种新的相似度度量——互相关熵(correntropy)的代价函数

d_i(k)表示向量D(k)的第i个元素，g_i(k,x(k))表示g(k,x(k))的第i个元素，L＝n+m，其中n为状态的维数，m为测量值的维数，其中σ表示核宽度。

步骤3-4

在最大互相关熵准则下，状态量的优化值可以由下式得到

其中，e_i(k)＝d_i(k)-g_i(k,x(k))。

步骤3-5

实际上该过程(即步骤3-4)是对残差e(k)的协方差进行加权，因此更新的协方差表示为

其中

取值则容易得到

步骤3-6

将代替原来的R(k)(步骤3-1中的R(k))，由步骤1中获取的相对测量信息作为测量值，根据无迹变化就可以计算出后验的相对状态估计值和协方差矩阵P(k|k)。

步骤4：采用基于最大互相关熵准则无迹卡尔曼滤波方法，对步骤1以及步骤2中所述的应用场合、运动模型和测量关系，构造相对状态估计器，从而得到两航天器间的相对位置和相对速度信息，参见图3。

步骤4-1

建立估计器的数据处理机制，并初始化估计器。具体流程为：将绝对状态向量和相对状态向量各自以时间标签进行标注，用下标符号‘_k’表示第k个时间标签的数据，用下标符号‘_k|k-1’表示第k-1个与第k个时间标签之间的时间片段。初始时间标签为下标‘₀’。

初始化绝对状态向量r₀、相对状态向量估计值以及协方差矩阵P₀。其方法为：通信故障航天器通过绝对运动方程求解出自身绝对状态向量r₀，根据轨道信息大致估计出通信保障航天器的相对状态向量根据初值相对于真值的误差程度，初始化P₀为较大的值。

由于互相关熵是一种局部优化的相似度测量，为克服初始估计值与实际值误差偏大时，收敛速度慢的缺点，前面一段时间(k＝0,1,…,100)将采用传统的无迹卡尔曼滤波。

步骤4-2

从时间标签0开始，到时间标签100，递归计算估计状态，按以下①～⑤进行。

①状态预测。在时间标签为k-1的时刻，根据当前的r_k-1、P_k-1，利用式(3)，式(6)进行状态预测以及无迹变化，得到P_k|k-1。

②测量预测。在时间标签为k-1的时刻，根据预测的状态和协方差矩阵P_k|k-1，利用式(1)以及无迹变化进一步预测测量值，得到P_y。

③准备更新。在时间标签为k-1的时刻，根据无迹变化，得到和之间的协方差P_xy，计算增益矩阵K_k。

④状态更新。在时间标签为k的时刻，得到测量值y_k，利用测量值y_k，按下式进行状态更新

⑤时间更新，递归计算。

将时间标签由k-1修改为k，重复步骤①～⑤。

步骤4-3

从时间标签101开始，根据步骤3的方法递归计算估计状态。

步骤4-4

根据步骤4-2和4-3的结果，计算出两个航天器间的相对位置和相对速度。具体方法为：在任意时间标签为k的时刻，由相对状态向量或中取出x、y、z作为在希尔坐标系中X_c、Y_c、Z_c轴上的相对位置，取出作为在希尔坐标系中X_c、Y_c、Z_c轴上的相对速度。

从图4可以看出，在冲激噪声的影响下，本方法(MCUKF)比传统的扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF)在相对距离估计误差(MSD_p)上具有更好的效果。

从图5可以看出，在冲激噪声的影响下，本方法(MCUKF)比传统的扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF)在相对速度估计误差(MSD_v)上具有更好的效果。

总之，本发明利用无线电测量设备测量原理，建立了通信保障航天器和通信故障航天器的相对测量方程，利用运动学方程，建立了通信保障航天器和通信故障航天器的相对状态方程(相对运动方程)，在此基础之上，结合无迹卡尔曼滤波和最大互相关熵准则的优势，提出一种基于最大互相关熵准则无迹卡尔曼滤波的方法来对相对状态进行估计，本发明在具有冲激噪声的实际系统中，能有效提高相对位置和相对速度的估计精度，进而提高通信保障航天器天线指向精度，保障通信链路和质量，防止通信中断。本发明不仅提高了传统算法的性能，还可以实现自主相对导航，减轻了地面测控的负担，还可以实现全天候、全天时工作。

Claims (6)

1.基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，其特征在于：包括以下步骤：

(1)利用安装在通信故障航天器上的无线电设备测出通信故障航天器和通信保障航天器的相对测量信息；

(2)通过运动学方程获取在地心惯性坐标系下通信故障航天器的绝对运动状态变量，并求得在希尔坐标系中通信保障航天器和通信故障航天器的相对运动方程；

(3)以互相关熵作为相似度度量并将最大互相关熵准则与无迹卡尔曼滤波结合，形成基于最大互相关熵准则的无迹卡尔曼滤波方法；

(4)将相对测量信息作为测量值，采用基于最大互相关熵准则的无迹卡尔曼滤波，计算出通信故障航天器和通信保障航天器之间的相对状态的估计值。

2.根据权利要求1所述基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，其特征在于：所述相对测量信息包括相对距离ρ、方位角θ和仰角φ：

其中，x、y和z分别为以通信故障航天器为中心的希尔坐标系的X_c、Y_c和Z_c轴的坐标，坐标原点固定于航天器所在轨道位置，X_c轴沿地心指向航天器质心的矢量方向，Z_c轴与航天器轨道角动量正方向一致，Y_c轴垂直于轨道平面并与X_c和Z_c轴构成右手正交坐标系。

3.根据权利要求2所述基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，其特征在于：所述步骤(2)具体包括以下步骤：

(2-1)在地心惯性坐标系下通信故障航天器的绝对运动状态变量的运动方程为：

其中，r_c和ω为地心惯性坐标系下通信故障航天器的绝对角加速度、绝对位置和绝对角速度，μ为地心引力常数；通过构造绝对状态向量为用函数等价表示上式，然后进行离散化，得到离散形式的运动方程，和分别为绝对加速度和绝对速度；

(2-2)在以通信故障航天器为中心的希尔坐标系下通信保障航天器的相对运动方程为

其中，和分别为X_c、Y_c和Z_c轴上相应的加速度；通过构造相对状态向量为 和分别表示X_c和Y_c轴上相应的速度，用函数等价表示上式，然后进行离散化，得到离散形式的运动方程。

4.根据权利要求1所述基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，其特征在于：所述基于最大互相关熵准则的无迹卡尔曼滤波包括以下步骤：

(3-1)根据无迹变化计算出先验的相对状态估计值和协方差矩阵P(k|k-1)，并采用Cholesky分解获得T(k)

其中，R(k)为测量噪声协方差；

(3-2)根据步骤(3-1)构造以下式子：

D(k)＝g(k,x(k))+e(k)

其中， x(k)为真实状态值，y(k)为测量值，h(·)为测量函数，r(k)为测量噪声；

(3-3)由步骤(3-2)得到基于互相关熵的代价函数

其中，d_i(k)表示向量D(k)的第i个元素，g_i(k,x(k))表示g(k,x(k))的第i个元素，L＝n+m其中n为状态的维数，m为测量值的维数；其中σ表示核宽度；

(3-4)在最大互相关熵准则下，状态量的优化值由下式得到

其中，e_i(k)＝d_i(k)-g_i(k,x(k))；

(3-5)更新的协方差表示为

其中

取则得到

(3-6)将代替R(k)，根据无迹变化就可以计算出后验的相对状态估计值和协方差矩阵P(k|k)。

5.根据权利要求1所述基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，其特征在于：所述步骤(4)还包括以下步骤：在利用基于最大互相关熵准则的无迹卡尔曼滤波进行相对状态估计前，通过使用无迹卡尔曼滤波加快估计初期的收敛速度。

6.根据权利要求5所述基于最大互相关熵准则无迹卡尔曼滤波的通信保障航天器相对状态估计方法，其特征在于：所述初期是指时间标签≤100。