CN113792411B - 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法 - Google Patents

一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法 Download PDF

Info

Publication number
CN113792411B
CN113792411B CN202110929245.2A CN202110929245A CN113792411B CN 113792411 B CN113792411 B CN 113792411B CN 202110929245 A CN202110929245 A CN 202110929245A CN 113792411 B CN113792411 B CN 113792411B
Authority
CN
China
Prior art keywords
spacecraft
state
representing
attitude
covariance
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110929245.2A
Other languages
English (en)
Other versions
CN113792411A (zh
Inventor
曹璐
杨宝健
冉德超
肖冰
蒋臣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Defense Technology Innovation Institute PLA Academy of Military Science
Original Assignee
National Defense Technology Innovation Institute PLA Academy of Military Science
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Defense Technology Innovation Institute PLA Academy of Military Science filed Critical National Defense Technology Innovation Institute PLA Academy of Military Science
Priority to CN202110929245.2A priority Critical patent/CN113792411B/zh
Publication of CN113792411A publication Critical patent/CN113792411A/zh
Application granted granted Critical
Publication of CN113792411B publication Critical patent/CN113792411B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Feedback Control In General (AREA)

Abstract

本发明公开了一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，包括：根据航天器测量数据和航天器姿态动力学模型，建立航天器姿态确定的非线性系统；根据前一时刻航天器状态和状态协方差，利用预设采样方式生成多个Sigma点，建立时间更新传递公式，获取当前时刻航天器的一步预测状态估计值和一步预测状态协方差；根据一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，获取航天器量测输出量的一步预测值、自协方差和互协方差；基于中心误差熵准则建立航天器状态的线性化回归方程，确定中心误差熵准则滤波的代价函数，得到当前时刻航天器的状态和状态协方差。本发明能提高处理非高斯噪声时的姿态估计精度和鲁棒性。

Description

一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定 方法
技术领域
本发明涉及航天器姿态估计技术领域，具体涉及一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法。
背景技术
高精度高可靠的姿态确定是航天器开展空间在轨服务等任务的基础。目前航天器的姿态确定方法按姿态解算方法不同可以分为确定性方法和状态估计法两类，其中，状态估计法采用滤波的方法根据观测信息对航天器状态量进行估计，可以有效的克服参考矢量的不确定性影响。
在非线性姿态估计过程中，主要使用扩展卡尔曼滤波(EKF)算法进行姿态估计。然而，由于扩展卡尔曼滤波本身的局限性导致其在强非线性条件下滤波精度不高。为了克服使用扩展卡尔曼滤波算法所存在的问题，目前提出了一种无迹卡尔曼滤波(UKF)算法，UKF算法基于无损变换(Unscented transform)，相比于EKF算法，在处理非线性问题上具有更高的精度，在高斯噪声的条件下有着良好的滤波效果。然而，在实际工程实践中往往会出现非高斯噪声，在面对工程中经常出现的野值、脉冲等诱导的非高斯噪声时，由于最小均方误差(MMSE)不具有鲁棒性，传统的UKF算法的滤波精度会降低甚至出现滤波发散的现象，无法满足滤波的精度要求。
为处理非高斯噪声，目前主要使用非高斯滤波器，非高斯滤波器包括：粒子滤波器(PF)、高斯和滤波器(GSP)、Huber无迹滤波器、Student’s t滤波器、最大相关熵无迹滤波器(MCUKF)和最小误差熵无迹滤波器(MEEUKF)。其中，粒子滤波器采用序贯重要性采样方法来近似计算后验密度，可以处理任意非高斯噪声；高斯和滤波器将非高斯噪声建模为高斯和分布以处理非高斯噪声；Huber无迹滤波器由Unscented变换与Huber代价函数结合而成，可以处理非线性非高斯系统；Student’s t滤波器将非高斯噪声假设为student’s t分布以处理非高斯噪声；最大相关熵无迹滤波器和最小误差熵无迹滤波器分别以最大相关熵准则和最小误差熵准则为最优准则，相比传统的最小均方误差准则有着更好的非高斯噪声的处理效果。
然而，上述的非高斯滤波器中，粒子滤波器具有较高的计算复杂度，且其粒子退化和粒子贫化问题也是很难处理的问题；高斯和滤波器假设噪声的概率密度是已知的，且有着较大的计算量；Huber无迹滤波器的精度有限；Student’s t滤波器只能对特定的非高斯噪声进行处理，环境适应性较差；最大相关熵无迹滤波器虽然可以将概率密度函数的峰值固定到零，但是单独使用时算法的精度有限；最小误差熵无迹滤波器具有平移不变性，不能保证估计误差收敛到零。
发明内容
为解决上述现有技术中存在的部分或全部技术问题，本发明提供一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法。
本发明的技术方案如下：
提供了一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，所述方法用于估计航天器姿态，包括：
根据航天器的测量数据和航天器姿态动力学模型，建立航天器姿态确定的非线性系统；
根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，建立时间更新传递公式，获取当前时刻航天器的一步预测状态估计值和一步预测状态协方差；
根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差；
基于中心误差熵准则建立航天器状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差。
在一些可能的实现方式中，建立航天器姿态确定的非线性系统为：
Figure BDA0003210645090000021
其中，xk表示k时刻航天器的n维状态向量，f(·)表示系统的状态方程，xk-1表示k-1时刻航天器的n维状态向量，ωk-1表示k-1时刻的n维系统噪声序列，zk表示k时刻的m维测量向量，h(·)表示系统的测量方程，νk表示k时刻的m维测量噪声序列。
在一些可能的实现方式中，设定前一时刻为k-1时刻，当前时刻为k时刻，根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，包括：
利用以下公式六生成2n+1个Sigma点；
Figure BDA0003210645090000031
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure BDA0003210645090000032
Figure BDA0003210645090000033
其中，ξ0,k-1、ξi,k-1和ξi+n,k-1表示由k-1时刻航天器的状态和状态协方差矩阵生成的Sigma点，
Figure BDA0003210645090000034
表示k-1时刻航天器的状态，Pk-1表示k-1时刻航天器的状态协方差矩阵，
Figure BDA0003210645090000035
表示(n+λ)Pk-1平方根矩阵的第i列，n表示系统状态维数，λ＝α2(n+κ)-n，α表示正的比例缩放因子，κ表示比例系数，Wi m表示一阶加权系数，Wi c表示二阶加权系数，β表示非负的权系数。
在一些可能的实现方式中，基于非线性系统和多个Sigma点，建立时间更新传递公式为：
Figure BDA0003210645090000036
其中，
Figure BDA0003210645090000037
表示k时刻航天器的一步预测状态估计值，Pk|k-1表示k时刻航天器的一步预测状态协方差矩阵，Qk-1表示系统噪声ωk-1的协方差矩阵。
在一些可能的实现方式中，根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，包括：
利用以下公式八生成2n+1个Sigma点；
Figure BDA0003210645090000041
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure BDA0003210645090000042
Figure BDA0003210645090000043
其中，ξ0,k|k-1、ξi,k|k-1和ξi+n,k|k-1表示由k时刻航天器的一步预测状态估计值和一步预测状态协方差矩阵生成的Sigma点，
Figure BDA0003210645090000044
表示(n+λ)Pk|k-1平方根矩阵的第i列。
在一些可能的实现方式中，利用以下公式九至公式十二获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差；
χi,k|k-1＝h(ξi,k|k-1)i＝0,1,…,2n公式九
Figure BDA0003210645090000045
Figure BDA0003210645090000046
Figure BDA0003210645090000047
其中，
Figure BDA0003210645090000048
表示航天器量测输出量的一步预测值，
Figure BDA0003210645090000049
表示航天器量测输出量的一步预测值的自协方差阵，
Figure BDA00032106450900000410
表示航天器量测输出量的一步预测值与航天器的一步预测状态估计值的互协方差阵，Rk表示测量噪声νk的协方差矩阵。
在一些可能的实现方式中，设定一步预测误差为：
Figure BDA00032106450900000411
设定伪观测矩阵为：
Figure BDA00032106450900000412
将测量向量zk近似为：
Figure BDA00032106450900000413
建立航天器状态对应的线性化回归方程为：
Figure BDA0003210645090000051
Figure BDA0003210645090000052
其中，
Figure BDA0003210645090000053
I表示单位矩阵，rk表示高阶误差项。
在一些可能的实现方式中，中心误差熵准则滤波的代价函数为：
Figure BDA0003210645090000054
其中，η表示权重系数，L＝m+n，
Figure BDA0003210645090000055
表示核宽为σ1的高斯核函数，ei,k表示误差变量ek的第i维状态，
Figure BDA0003210645090000056
表示核宽为σ2的高斯核函数，ej,k表示误差变量ek的第j维状态。
在一些可能的实现方式中，利用以下公式二十一确定当前时刻航天器的状态；
Figure BDA0003210645090000057
其中，
Figure BDA0003210645090000058
表示k时刻航天器的状态的最优估计值，
Figure BDA0003210645090000059
在一些可能的实现方式中，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差，包括以下步骤：
对代价函数计算梯度，将梯度计算公式表达成矩阵形式，并令梯度等于0；
基于矩阵形式的梯度计算公式，采用定点迭代算法，得到当前时刻航天器的状态和状态协方差。
本发明技术方案的主要优点如下：
本发明的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法通过利用无损变换进行一步预测状态估计值和一步预测状态协方差的获取，然后构造线性化回归方程，利用中心误差熵准则进行航天器的后验状态的求解，能够有效应对航天器姿态确定的非线性系统中出现的非高斯噪声，提高处理非高斯噪声时的航天器姿态估计精度和鲁棒性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明一实施例的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法的流程图；
图2为采用传统无迹卡尔曼滤波算法、最大相关熵无迹卡尔曼滤波算法、最小误差熵无迹卡尔曼滤波算法和本发明一实施例的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法所得到的航天器滚转角的均方根误差对比示意图；
图3为采用传统无迹卡尔曼滤波算法、最大相关熵无迹卡尔曼滤波算法、最小误差熵无迹卡尔曼滤波算法和本发明一实施例的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法所得到的航天器俯仰角的均方根误差对比示意图；
图4为采用传统无迹卡尔曼滤波算法、最大相关熵无迹卡尔曼滤波算法、最小误差熵无迹卡尔曼滤波算法和本发明一实施例的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法所得到的航天器偏航角的均方根误差对比示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明具体实施例及相应的附图对本发明技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
以下结合附图，详细说明本发明实施例提供的技术方案。
参见图1，本发明一实施例提供了一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，该方法用于估计航天器姿态，包括以下步骤：
S1，根据航天器的测量数据和航天器姿态动力学模型，建立航天器姿态确定的非线性系统；
S2，根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，建立时间更新传递公式，获取当前时刻航天器的一步预测状态估计值和一步预测状态协方差；
S3，根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差；
S4，基于中心误差熵准则建立航天器状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差。
以下以前一时刻为k-1时刻，当前时刻为k时刻为例，对本发明一实施例提供的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法的步骤及原理进行具体说明。
S1，根据航天器的测量数据和航天器姿态动力学模型，建立航天器姿态确定的非线性系统。
具体地，根据航天器的测量数据和航天器姿态动力学模型，建立航天器姿态确定的非线性系统为：
Figure BDA0003210645090000071
其中，xk表示k时刻航天器的n维状态向量，f(·)表示系统的状态方程，xk-1表示k-1时刻航天器的n维状态向量，ωk-1表示k-1时刻的n维系统噪声序列，zk表示k时刻的m维测量向量，h(·)表示系统的测量方程，νk表示k时刻的m维测量噪声序列，且ωk和νk互不相关。
设定：航天器初始状态x0与ωk和νk相独立，ωk与νk相互独立，ωk和νk的统计特性如下：
Figure BDA0003210645090000072
其中，E(·)表示数学期望，Qk表示系统噪声ωk的协方差矩阵，ωk表示k时刻的n维系统噪声序列，Rk表示测量噪声νk的协方差矩阵。
S2，根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，建立时间更新传递公式，获取当前时刻航天器的一步预测状态估计值和一步预测状态协方差。
本发明一实施例中，预设采样方式包括：
利用以下公式三生成2n+1个Sigma点；
Figure BDA0003210645090000081
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure BDA0003210645090000082
Figure BDA0003210645090000083
其中，ξ0、ξi和ξi+n表示生成的Sigma点，x表示输入的状态均值，Px表示输入的状态均值对应的状态协方差矩阵，
Figure BDA0003210645090000084
表示(n+λ)Px平方根矩阵的第i列，n表示系统状态维数，λ＝α2(n+κ)-n，α表示正的比例缩放因子，κ表示比例系数，Wi m表示一阶加权系数，Wi c表示二阶加权系数，β表示非负的权系数。
当根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点时，则将公式三中的状态均值x和状态协方差矩阵Px替换为前一时刻航天器的状态和状态协方差矩阵。
具体地，以前一时刻为k-1时刻为例，根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，包括：
利用以下公式六生成2n+1个Sigma点；
Figure BDA0003210645090000085
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure BDA0003210645090000086
Figure BDA0003210645090000091
其中，ξ0,k-1、ξi,k-1和ξi+n,k-1表示由k-1时刻航天器的状态和状态协方差矩阵生成的Sigma点，
Figure BDA0003210645090000092
表示k-1时刻航天器的状态，Pk-1表示k-1时刻航天器的状态协方差矩阵，
Figure BDA0003210645090000093
表示(n+λ)Pk-1平方根矩阵的第i列。
进一步地，基于非线性系统和多个Sigma点，建立时间更新传递公式为：
Figure BDA0003210645090000094
其中，
Figure BDA0003210645090000095
表示k时刻航天器的一步预测状态估计值，Pk|k-1表示k时刻航天器的一步预测状态协方差矩阵，Qk-1表示系统噪声ωk-1的协方差矩阵。
利用上述的时间更新传递公式进行一步预测，能够获取k时刻航天器的一步预测状态估计值和一步预测状态协方差，即当前时刻航天器的一步预测状态估计值和一步预测状态协方差。
S3，根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差。
具体地，以当前时刻为k时刻为例，根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，包括：
利用以下公式八生成2n+1个Sigma点；
Figure BDA0003210645090000096
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure BDA0003210645090000101
Figure BDA0003210645090000102
其中，ξ0,k|k-1、ξi,k|k-1和ξi+n,k|k-1表示由k时刻航天器的一步预测状态估计值和一步预测状态协方差矩阵生成的Sigma点，
Figure BDA0003210645090000103
表示(n+λ)Pk|k-1平方根矩阵的第i列。
进一步地，基于非线性系统和多个Sigma点，利用以下公式九至公式十二获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差；
χi,k|k-1＝h(ξi,k|k-1)i＝0,1,…,2n 公式九
Figure BDA0003210645090000104
Figure BDA0003210645090000105
Figure BDA0003210645090000106
其中，
Figure BDA0003210645090000107
表示航天器量测输出量的一步预测值，
Figure BDA0003210645090000108
表示航天器量测输出量的一步预测值的自协方差阵，
Figure BDA0003210645090000109
表示航天器量测输出量的一步预测值与航天器的一步预测状态估计值的互协方差阵。
本发明一实施例中，通过利用非线性测量方程h(·)将采样得到的Sigma点ξi,k|k-1(i＝0,1,…,2n)传播为χi,k|k-1，利用χi,k|k-1获取航天器量测输出量的一步预测值及其自协方差、以及其与航天器的一步预测状态估计值的互协方差。
S4，基于中心误差熵准则建立航天器状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差。
具体地，基于中心误差熵准则建立航天器状态对应的线性化回归方程，包括：
定义一步预测误差为：
Figure BDA0003210645090000111
定义伪观测矩阵为：
Figure BDA0003210645090000112
将测量向量zk近似为：
Figure BDA0003210645090000113
建立航天器状态对应的线性化回归方程为：
Figure BDA0003210645090000114
其中，rk表示高阶误差项，
Figure BDA0003210645090000115
I表示单位矩阵。
进一步地，设定：
Figure BDA0003210645090000116
其中，Sk、Sp,k|k-1和Sr,k分别表示矩阵
Figure BDA0003210645090000117
Pk|k-1和Rk的Cholesky分解。
对上述公式十六所示的线性化回归方程的等式两边同乘
Figure BDA0003210645090000118
以将线性化回归方程变形为：
dk＝Wkxk+ek 公式十八
其中，
Figure BDA0003210645090000119
进一步地，设定：ek＝[e1,k,e2,k,…,eL,k]T，dk＝[d1,k,d2,k,…,dL,k]T，Wk＝[w1,k,w2,k,…,wL,k]T，ei,k＝di,k-wi,kxk(i＝1,…,L)，L＝m+n，ei,k表示ek的第i个元素，di,k表示dk的第i个元素，wi,k表示Wk的第i行向量；
则中心误差熵准则滤波(CEEKF)的代价函数为：
Figure BDA0003210645090000121
其中，η表示权重系数，
Figure BDA0003210645090000122
表示核宽为σ1的高斯核函数，
Figure BDA0003210645090000123
表示核宽为σ2的高斯核函数。
设定：
Figure BDA0003210645090000124
则公式十九可以表述为：
Figure BDA0003210645090000125
在中心误差熵(centered error entropy，CEE)准则下，可以通过最大化上述代价函数得到当前时刻航天器的状态的最优估计值，该最优估计值即为估计得到的当前时刻航天器的状态。
具体地，以当前时刻为k时刻为例，k时刻航天器的状态可以通过以下公式二十一确定；
Figure BDA0003210645090000126
其中，
Figure BDA0003210645090000127
表示k时刻航天器的状态的最优估计值，即k时刻航天器的状态，
Figure BDA0003210645090000128
表示JL(xk)取得最大值时对应的xk值。
进一步地，本发明一实施例中，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差，可以包括以下步骤：
对代价函数计算梯度，将梯度计算公式表达成矩阵形式，并令梯度等于0；
基于矩阵形式的梯度计算公式，采用定点迭代算法，得到当前时刻航天器的状态和状态协方差。
具体地，以当前时刻为k时刻为例，可以利用以下公式二十二对代价函数计算梯度，并令梯度等于0；
Figure BDA0003210645090000131
其中，
Figure BDA0003210645090000132
Figure BDA0003210645090000133
将上述公式二十二所示的代价函数梯度计算公式表达成如以下公式二十三所示的矩阵形式；
Figure BDA0003210645090000134
其中，
Figure BDA0003210645090000135
Figure BDA0003210645090000136
(Ωk)ij表示Ωk的第i行第j列元素，
Figure BDA0003210645090000137
基于矩阵形式的代价函数梯度计算公式，采用定点迭代算法，得到k时刻航天器的状态为：
Figure BDA0003210645090000141
其中，
Figure BDA0003210645090000142
表示第t+1次迭代的k时刻航天器的状态。
进一步地，进行如以下公式二十五至公式三十四的设定：
Figure BDA0003210645090000143
Figure BDA0003210645090000144
Figure BDA0003210645090000145
Figure BDA0003210645090000146
Figure BDA0003210645090000147
Figure BDA0003210645090000148
Figure BDA0003210645090000149
Figure BDA00032106450900001410
Figure BDA00032106450900001411
Figure BDA00032106450900001412
其中，Λx,k表示n×n维矩阵，Λxy,k为m×n维矩阵，Λyx,k为n×m维矩阵，Λy,k为m×m维矩阵，Λi,j；k表示Λk矩阵的第i行第j列组成的矩阵，Ξi,j；k表示Ξk矩阵的第i行第j列组成的矩阵，Ωi,j；k表示Ωk矩阵的第i行第j列组成的矩阵。
根据公式十八和设定的公式二十五至公式三十四，可以将k时刻航天器的状态表示为：
Figure BDA0003210645090000151
其中，
Figure BDA0003210645090000152
Figure BDA0003210645090000153
Figure BDA0003210645090000154
利用公式三十五获取的
Figure BDA0003210645090000155
为k时刻航天器状态xk的最优估计值，可作为估计得到的当前k时刻航天器的状态。
相应的，基于上述设定，当前k时刻航天器的状态协方差矩阵Pk可以更新为：
Figure BDA0003210645090000156
其中，I表示单位矩阵。
本发明一实施例提供的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法通过利用无损变换进行一步预测状态估计值和一步预测状态协方差的获取，然后构造线性化回归方程，利用中心误差熵准则进行航天器的后验状态的求解，能够有效应对航天器姿态确定的非线性系统中出现的非高斯噪声，提高处理非高斯噪声时的航天器姿态估计精度和鲁棒性。
以下结合具体示例，对本发明一实施例提供的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法的有益效果进行说明。
某航天器姿态确定系统的系统状态方程和测量方程如公式三十七所示：
Figure BDA0003210645090000157
其中，qbo为四元数表示的航天器姿态，b为陀螺仪的常值漂移，也作为状态变量被估计，ωg为陀螺仪测量值，ωg为输入，ηg和ηb为零均值系统噪声，ηg和ηb的协方差为Qg和Qb，Ωd为传递矩阵，可以表示为：
Figure BDA0003210645090000161
Figure BDA0003210645090000162
qopt为航天器姿态敏感器的观测输出，qN为测量噪声四元数，ωoi＝[0 -ω0 0]表示惯性系下的轨道角速度，Ebo表示从轨道系到星体系的坐标转换矩阵，可以表示为：
Figure BDA0003210645090000163
用于航天器姿态确定的参数设定如表1所示；
表1
Figure BDA0003210645090000164
更具体地，陀螺常值漂移：b＝[30 30 30]T(°)/h，常值漂移白噪声均方差σb＝0.5(°)/h，陀螺测量噪声σg＝0.5(°)/h，滤波初值选取为：qbo(0)＝[0,0,0,1]T，b(0)＝[30 3030]T(°)/h，ωbo＝10-4×[cos(10ω0t) cos(8ω0t) cos(5.7ω0t)]，轨道角速度ω0＝0.0012rad/s，初始协方差P0＝diag(I3×3 0.04I3×3)，星敏感器量测噪声为混合高斯噪声，
Figure BDA0003210645090000171
附图2至附图4给出了利用传统无迹卡尔曼滤波算法(UKF)、最大相关熵无迹卡尔曼滤波算法(MCUKF)、最小误差熵无迹卡尔曼滤波算法(MEEUKF)和本发明一实施例提供的基于中心误差熵准则无迹卡尔曼滤波(CEEUKF)的航天器姿态确定方法得到的不同姿态估计结果对比示意图。附图中，RMSE of
Figure BDA0003210645090000172
滚转角均方根误差，RMSE ofθ：俯仰角均方根误差，RMSE ofψ：偏航角均方根误差，time：时间。
其中，各个算法的核宽参数选择如表2所示；
表2
Figure BDA0003210645090000173
得到的不同算法下三轴姿态角的平均的均方根误差(ARMSE)如表3所示；
表3
Algorithms UKF MEEUKF MCUKF CEEUKF
Row(deg) 18.829743 28.723249 10.161122 4.049165
Pitch(deg) 19.290676 28.752157 6.039377 3.564183
Yaw(deg) 10.095340 20.991913 4.576413 3.650843
可以看出，本发明一实施例提供的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法在非高斯噪声条件下具有最高的滤波精度，同时滤波收敛后具有最低的估计误差协方差，即最好的滤波稳定性，能够更好地应对非高斯噪声。
需要说明的是，在本文中，诸如“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。此外，本文中“前”、“后”、“左”、“右”、“上”、“下”均以附图中表示的放置状态为参照。
最后应说明的是：以上实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (10)

1.一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，所述方法用于估计航天器姿态，包括：
根据航天器的测量数据和航天器姿态动力学模型，建立航天器姿态确定的非线性系统；
根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，建立时间更新传递公式，获取当前时刻航天器的一步预测状态估计值和一步预测状态协方差；
根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，基于非线性系统和多个Sigma点，获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差；
基于中心误差熵准则建立航天器状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差。
2.根据权利要求1所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，建立航天器姿态确定的非线性系统为：
Figure FDA0003210645080000011
其中，xk表示k时刻航天器的n维状态向量，f(·)表示系统的状态方程，xk-1表示k-1时刻航天器的n维状态向量，ωk-1表示k-1时刻的n维系统噪声序列，zk表示k时刻的m维测量向量，h(·)表示系统的测量方程，νk表示k时刻的m维测量噪声序列。
3.根据权利要求2所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，设定前一时刻为k-1时刻，当前时刻为k时刻，根据前一时刻航天器的状态和状态协方差，利用预设采样方式生成多个Sigma点，包括：
利用以下公式六生成2n+1个Sigma点；
Figure FDA0003210645080000012
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure FDA0003210645080000021
Figure FDA0003210645080000022
其中，ξ0,k-1、ξi,k-1和ξi+n,k-1表示由k-1时刻航天器的状态和状态协方差矩阵生成的Sigma点，
Figure FDA0003210645080000023
表示k-1时刻航天器的状态，Pk-1表示k-1时刻航天器的状态协方差矩阵，
Figure FDA0003210645080000024
表示(n+λ)Pk-1平方根矩阵的第i列，n表示系统状态维数，λ＝α2(n+κ)-n，α表示正的比例缩放因子，κ表示比例系数，Wi m表示一阶加权系数，Wi c表示二阶加权系数，β表示非负的权系数。
4.根据权利要求3所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，基于非线性系统和多个Sigma点，建立时间更新传递公式为：
Figure FDA0003210645080000025
其中，
Figure FDA0003210645080000026
表示k时刻航天器的一步预测状态估计值，Pk|k-1表示k时刻航天器的一步预测状态协方差矩阵，Qk-1表示系统噪声ωk-1的协方差矩阵。
5.根据权利要求4所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，根据当前时刻航天器的一步预测状态估计值和一步预测状态协方差，利用预设采样方式生成多个Sigma点，包括：
利用以下公式八生成2n+1个Sigma点；
Figure FDA0003210645080000027
利用以下公式四和公式五确定每个Sigma点对应的加权系数；
Figure FDA0003210645080000031
Figure FDA0003210645080000032
其中，ξ0,k|k-1、ξi,k|k-1和ξi+n,k|k-1表示由k时刻航天器的一步预测状态估计值和一步预测状态协方差矩阵生成的Sigma点，
Figure FDA0003210645080000033
表示(n+λ)Pk|k-1平方根矩阵的第i列。
6.根据权利要求5所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，利用以下公式九至公式十二获取航天器量测输出量的一步预测值、一步预测值的自协方差、以及一步预测值与一步预测状态估计值的互协方差；
χi,k|k-1＝h(ξi,k|k-1)i＝0,1,…,2n公式九
Figure FDA0003210645080000034
Figure FDA0003210645080000035
Figure FDA0003210645080000036
其中，
Figure FDA0003210645080000037
表示航天器量测输出量的一步预测值，
Figure FDA0003210645080000038
表示航天器量测输出量的一步预测值的自协方差阵，
Figure FDA0003210645080000039
表示航天器量测输出量的一步预测值与航天器的一步预测状态估计值的互协方差阵，Rk表示测量噪声νk的协方差矩阵。
7.根据权利要求6所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，设定一步预测误差为：
Figure FDA00032106450800000310
设定伪观测矩阵为：
Figure FDA00032106450800000311
将测量向量zk近似为：
Figure FDA00032106450800000312
建立航天器状态对应的线性化回归方程为：
Figure FDA0003210645080000041
其中，
Figure FDA0003210645080000042
I表示单位矩阵，rk表示高阶误差项。
8.根据权利要求7所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，中心误差熵准则滤波的代价函数为：
Figure FDA0003210645080000043
其中，η表示权重系数，L＝m+n，
Figure FDA0003210645080000044
表示核宽为σ1的高斯核函数，ei,k表示误差变量ek的第i维状态，
Figure FDA0003210645080000045
表示核宽为σ2的高斯核函数，ej,k表示误差变量ek的第j维状态。
9.根据权利要求8所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，利用以下公式二十一确定当前时刻航天器的状态；
Figure FDA0003210645080000046
其中，
Figure FDA0003210645080000047
表示k时刻航天器的状态的最优估计值，
Figure FDA0003210645080000048
10.根据权利要求9所述的基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法，其特征在于，对代价函数进行最大化处理，得到当前时刻航天器的状态和状态协方差，包括以下步骤：
对代价函数计算梯度，将梯度计算公式表达成矩阵形式，并令梯度等于0；
基于矩阵形式的梯度计算公式，采用定点迭代算法，得到当前时刻航天器的状态和状态协方差。
CN202110929245.2A 2021-08-13 2021-08-13 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法 Active CN113792411B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110929245.2A CN113792411B (zh) 2021-08-13 2021-08-13 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110929245.2A CN113792411B (zh) 2021-08-13 2021-08-13 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法

Publications (2)

Publication Number Publication Date
CN113792411A CN113792411A (zh) 2021-12-14
CN113792411B true CN113792411B (zh) 2022-10-11

Family

ID=79181719

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110929245.2A Active CN113792411B (zh) 2021-08-13 2021-08-13 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法

Country Status (1)

Country Link
CN (1) CN113792411B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114417912B (zh) * 2021-12-20 2024-04-12 中国人民解放军军事科学院国防科技创新研究院 一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法
CN114880874B (zh) * 2022-06-07 2024-03-12 东南大学 一种水面无人船参数自适应鲁棒估计方法与系统
CN115166635B (zh) * 2022-06-24 2023-03-28 江南大学 基于风险敏感fir滤波的机器人定位方法
CN114858166B (zh) * 2022-07-11 2022-10-11 北京神导科技股份有限公司 基于最大相关熵卡尔曼滤波器的imu姿态解算方法
CN117705108A (zh) * 2023-12-12 2024-03-15 中铁第四勘察设计院集团有限公司 一种基于最大相关熵的滤波定位方法、系统及滤波器

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108520107A (zh) * 2018-03-19 2018-09-11 山西大学 基于最大似然准则鲁棒卡尔曼滤波的系统状态估计方法
CN109631913A (zh) * 2019-01-30 2019-04-16 西安电子科技大学 基于非线性预测强跟踪无迹卡尔曼滤波的x射线脉冲星导航定位方法及系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108520107A (zh) * 2018-03-19 2018-09-11 山西大学 基于最大似然准则鲁棒卡尔曼滤波的系统状态估计方法
CN109631913A (zh) * 2019-01-30 2019-04-16 西安电子科技大学 基于非线性预测强跟踪无迹卡尔曼滤波的x射线脉冲星导航定位方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
改进的平方根UKF及其在交会对接中的应用;李鹏等;《电机与控制学报》;20101115(第11期);全文 *

Also Published As

Publication number Publication date
CN113792411A (zh) 2021-12-14

Similar Documents

Publication Publication Date Title
CN113792411B (zh) 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法
Liu et al. Maximum correntropy square-root cubature Kalman filter with application to SINS/GPS integrated systems
Huang et al. A novel robust student's t-based Kalman filter
CN108225337B (zh) 基于sr-ukf滤波的星敏感器和陀螺组合定姿方法
CN106599368B (zh) 基于改进粒子提议分布和自适应粒子重采样的FastSLAM方法
Tang et al. Square-root quaternion cubature Kalman filtering for spacecraft attitude estimation
CN108844536B (zh) 一种基于量测噪声协方差矩阵估计的地磁导航方法
CN110109470A (zh) 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制系统
CN111798491A (zh) 一种基于Elman神经网络的机动目标跟踪方法
Dang et al. Cubature Kalman filter under minimum error entropy with fiducial points for INS/GPS integration
Tang et al. Square-root adaptive cubature Kalman filter with application to spacecraft attitude estimation
CN108562290B (zh) 导航数据的滤波方法、装置、计算机设备及存储介质
CN111983927A (zh) 一种新型的最大协熵椭球集员滤波方法
Qiu et al. Adaptive robust nonlinear filtering for spacecraft attitude estimation based on additive quaternion
CN110706265A (zh) 一种改进srckf强跟踪滤波的机动目标跟踪方法
CN113449384A (zh) 一种基于中心误差熵准则扩展卡尔曼滤波的姿态确定方法
Penenko et al. Sequential data assimilation algorithms for air quality monitoring models based on a weak-constraint variational principle
CN113587926A (zh) 一种航天器空间自主交会对接相对导航方法
Wang et al. Variational Bayesian cubature RTS smoothing for transfer alignment of DPOS
CN110006462B (zh) 基于奇异值分解的星敏感器在轨标定方法
CN113792412B (zh) 一种基于中心误差熵准则容积卡尔曼滤波的航天器姿态确定方法
Ma et al. Jet transport particle filter for attitude estimation of tumbling space objects
CN115453534B (zh) 一种顾及解缠误差的序贯InSAR时序形变解算方法
CN110045363B (zh) 基于相对熵的多雷达航迹关联方法
CN108225274A (zh) 基于不完全分解预处理的共轭梯度法光束法平差方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant