CN114417912B - 一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法 - Google Patents

Abstract

本发明公开了一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，包括：建立卫星姿态确定的非线性系统；根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，通过状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差；根据一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，通过测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、一步预测值与一步预测状态估计值的互协方差；基于中心误差熵准则建立卫星状态的线性化回归方程，获取当前时刻卫星的状态估计值和状态协方差。本发明能提高处理非高斯噪声时的姿态估计精度和鲁棒性。

Description

一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的 卫星姿态确定方法

技术领域

本发明涉及卫星姿态估计技术领域，具体涉及一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法。

背景技术

高精度高可靠的姿态确定是卫星开展空间在轨服务等任务的基础。目前卫星的姿态确定方法按姿态解算方法不同可以分为确定性方法和状态估计法两类，其中，状态估计法采用滤波的方法根据观测信息对卫星状态量进行估计，可以有效的克服参考矢量的不确定性影响。

在非线性姿态估计过程中，主要使用扩展卡尔曼滤波(EKF)算法进行姿态估计。然而，由于扩展卡尔曼滤波本身的局限性导致其在强非线性条件下滤波精度不高。为了克服使用扩展卡尔曼滤波算法所存在的问题，现有技术提出了一种Sigma点卡尔曼滤波(SPKF)算法。SPKF算法主要包括无迹卡尔曼滤波(UKF)算法和中心差分卡尔曼滤波(CDKF)算法，SPKF算法通过不同的采样准则进行Sigma点采样，然后基于这些Sigma点经过非线性传递后的结果来计算系统状态的均值和协方差，相比于EKF算法，在处理非线性问题上具有更高的精度，在高斯噪声的条件下有着良好的滤波效果。然而，在实际卫星姿态确定系统中，由于传感器故障、野值干扰等状况会导致噪声服从厚尾非高斯分布。此时，传统的CDKF算法的滤波精度会降低甚至出现滤波发散的现象，无法满足滤波的精度要求。

为处理非高斯噪声，目前主要使用非高斯滤波器，非高斯滤波器包括：粒子滤波器(PF)、Huber滤波器、最大相关熵卡尔曼滤波器(MCKF)、最小误差熵准则滤波器(MEEKF)和Student’s t滤波器。其中，粒子滤波器采用序贯重要性采样方法来近似计算后验密度，可以处理任意非线性模型及非高斯分布；Huber滤波器基于Huber代价函数对经典卡尔曼滤波算法进行了重构，可以处理非线性非高斯系统；最大相关熵卡尔曼滤波器和最小误差熵准则滤波器分别以最大相关熵准则和最小误差熵准则为最优准则推导得到，也可以处理复杂的非高斯噪声；Student’s t滤波器将非高斯噪声假设为student’s t分布，针对这种特定分布进行求解得到后验状态估计，可以处理非线性模型及非高斯分布。

然而，上述的非高斯滤波器中，粒子滤波器具有较高的计算复杂度，且其粒子退化和粒子贫化问题也是很难处理的问题；Huber滤波器在处理非高斯噪声时的精度较低；最大相关熵卡尔曼滤波器严重依赖核宽，在面对更复杂的非高斯噪声时会出现精度较低的问题；最小误差熵准则滤波器具有平移不变性，无法使误差概率密度函数趋于零，并且容易出现滤波发散，存在着一定的不稳定性；Student’s t滤波器只能对特定的非高斯噪声进行处理，环境适应性较差。

发明内容

为解决上述现有技术中存在的部分或全部技术问题，本发明提供一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法。

本发明的技术方案如下：

提供了一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，包括：

根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统；

根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差；

根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差；

基于中心误差熵准则建立卫星状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差。

在一些可能的实现方式中，根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统为：

其中，x_k表示k时刻卫星的n维状态向量，f(·)表示系统的状态方程，x_k-1表示k-1时刻卫星的n维状态向量，ω_k-1表示k-1时刻的n维过程噪声序列，z_k表示k时刻的m维测量向量，h(·)表示系统的测量方程，ν_k表示k时刻的m维测量噪声序列。

在一些可能的实现方式中，设定前一时刻为k-1时刻，当前时刻为k时刻，根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，包括：

基于确定性采样点对称采样策略，根据前一时刻卫星的状态估计值和状态协方差，利用以下公式三计算2n+1个确定性采样点；

利用以下公式四确定每个确定性采样点对应的加权系数；

其中，ξ_i,k-1,i＝0,1,...,2n表示由k-1时刻卫星的状态估计值和状态协方差矩阵计算得到的确定性采样点，表示k-1时刻卫星的状态估计值，P_k-1表示k-1时刻卫星的状态协方差矩阵，表示矩阵的第i列，n表示系统状态维数，h表示调节因子，表示第1个采样点均值的加权系数，W_i ^m表示第i+1个采样点均值的加权系数，W_i ^c1表示第i+1个采样点第1类加权系数，W_i ^c2表示第i+1个采样点第2类加权系数。

在一些可能的实现方式中，基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差，包括：

基于获取的确定性采样点，利用以下公式五进行状态传递，获取一步预测状态量的2n+1个采样点；

γ_i,k|k-1＝f(ξ_i,k-1),i＝0,1,2,...,2n公式五

基于一步预测状态量的2n+1个采样点，利用以下公式六进行状态预测估计，获取当前时刻卫星的一步预测状态估计值；

基于一步预测状态量的2n+1个采样点，利用以下公式七获取当前时刻卫星的一步预测状态协方差；

其中，γ_i,k|k-1表示一步预测状态量的第i+1个采样点，表示k时刻卫星的一步预测状态估计值，P_k|k-1表示k时刻卫星的一步预测状态协方差矩阵，Q表示过程噪声协方差矩阵。

在一些可能的实现方式中，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，包括：

基于确定性采样点对称采样策略，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用以下公式八计算2n+1个确定性采样点；

利用以下公式四确定每个确定性采样点对应的加权系数；

其中，ξ_i,k|k-1,i＝0,1,...,2n表示由k时刻卫星的一步预测状态估计值和一步预测状态协方差矩阵P_k|k-1计算得到的确定性采样点，表示矩阵的第i列。

在一些可能的实现方式中，基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差，包括：

基于获取的确定性采样点，利用以下公式九进行测量信息传递，获取2n+1个经测量方程传递的一步预测采样点；

χ_i,k|k-1＝h(ξ_i,k|k-1),i＝0,1,2,...,2n公式九

基于2n+1个经测量方程传递的一步预测采样点，利用以下公式十进行状态预测估计，获取卫星量测输出量的一步预测值；

基于2n+1个经测量方程传递的一步预测采样点，利用以下公式十一获取卫星量测输出量的一步预测值与一步预测状态估计值的互协方差；

其中，χ_i,k|k-1表示第i+1个经测量方程传递的一步预测采样点，表示卫星量测输出量的一步预测值，表示卫星量测输出量的一步预测值与航天器的一步预测状态估计值的互协方差阵，χ_1:n,k|k-1表示i从1到n的前n个经测量方程传递的一步预测采样点χ_i,k|k-1，χ_n+1:2n,k|k-1表示i从n+1到2n的后n个经测量方程传递的一步预测采样点χ_i,k|k-1。

在一些可能的实现方式中，设定一步预测误差为：

设定量测斜率矩阵为：

将测量向量z_k近似为：

建立卫星状态对应的线性化回归方程为：

其中，I表示单位矩阵，r_k表示高阶误差项。

在一些可能的实现方式中，中心误差熵准则滤波的代价函数为：

其中，λ表示权重系数，L＝m+n，表示核宽为σ₁的高斯核函数，e_i,k表示误差变量e_k的第i维状态，表示核宽为σ₂的高斯核函数，e_j,k表示误差变量e_k的第j维状态。

在一些可能的实现方式中，利用以下公式二十确定当前时刻卫星的状态；

其中，表示k时刻卫星的状态的最优估计值，

在一些可能的实现方式中，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差，包括以下步骤：

对代价函数计算梯度，将梯度计算公式表达成矩阵形式，并令梯度等于0；

基于矩阵形式的梯度计算公式，利用定点迭代算法和矩阵求逆引理，得到当前时刻卫星的状态估计值和状态协方差。

本发明技术方案的主要优点如下：

本发明的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法通过利用Stirling插值变换进行一步预测状态估计值和一步预测状态协方差的获取，然后构造线性化回归方程，利用中心误差熵准则进行卫星的后验状态的更新求解，能够有效应对卫星姿态确定的非线性系统中出现的野值非高斯噪声，提高处理野值非高斯噪声时的卫星姿态估计精度和鲁棒性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一实施例的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法的流程图；

图2为采用传统中心差分卡尔曼滤波算法、最大相关熵中心差分卡尔曼滤波算法、最小误差熵中心差分卡尔曼滤波算法和本发明一实施例的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法所得到的卫星滚转角的均方根误差对比示意图；

图3为采用传统中心差分卡尔曼滤波算法、最大相关熵中心差分卡尔曼滤波算法、最小误差熵中心差分卡尔曼滤波算法和本发明一实施例的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法所得到的卫星俯仰角的均方根误差对比示意图；

图4为采用传统中心差分卡尔曼滤波算法、最大相关熵中心差分卡尔曼滤波算法、最小误差熵中心差分卡尔曼滤波算法和本发明一实施例的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法所得到的卫星偏航角的均方根误差对比示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明具体实施例及相应的附图对本发明技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

以下结合附图，详细说明本发明实施例提供的技术方案。

参见图1，本发明一实施例提供了一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，该方法用于估计卫星姿态，包括以下步骤：

S1，根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统；

S2，根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差；

S3，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差；

S4，基于中心误差熵准则建立卫星状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差。

以下以前一时刻为k-1时刻，当前时刻为k时刻为例，对本发明一实施例提供的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法的步骤及原理进行具体说明。

步骤S1，根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统。

具体地，根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统为：

其中，x_k表示k时刻卫星的n维状态向量，f(·)表示系统的状态方程，x_k-1表示k-1时刻卫星的n维状态向量，ω_k-1表示k-1时刻的n维过程噪声序列，z_k表示k时刻的m维测量向量，h(·)表示系统的测量方程，ν_k表示k时刻的m维测量噪声序列，且ω_k和ν_k互不相关。

设定：卫星初始状态x₀与ω_k和ν_k相独立，ω_k与ν_k相互独立，ω_k和ν_k的统计特性如下：

其中，E(·)表示数学期望，Q_k表示过程噪声ω_k的协方差矩阵，ω_k表示k时刻的n维过程噪声序列，R_k表示测量噪声ν_k的协方差矩阵。

步骤S2，根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差。

具体地，以前一时刻为k-1时刻为例，根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，包括：

基于确定性采样点对称采样策略，根据前一时刻卫星的状态估计值和状态协方差，利用以下公式三计算2n+1个确定性采样点；

利用以下公式四确定每个确定性采样点对应的加权系数；

其中，ξ_i,k-1,i＝0,1,...,2n表示由k-1时刻卫星的状态估计值和状态协方差矩阵计算得到的确定性采样点，表示k-1时刻卫星的状态估计值，P_k-1表示k-1时刻卫星的状态协方差矩阵，表示矩阵的第i列，n表示系统状态维数，表示第1个采样点均值的加权系数，W_i ^m表示第i+1个采样点均值的加权系数，W_i ^c1表示第i+1个采样点第1类加权系数，W_i ^c2表示第i+1个采样点第2类加权系数，h表示调节因子。

进一步地，基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差，可以包括：

基于获取的确定性采样点，利用以下公式五进行状态传递，获取一步预测状态量的2n+1个采样点；

γ_i,k|k-1＝f(ξ_i,k-1),i＝0,1,2,...,2n公式五

基于一步预测状态量的2n+1个采样点，利用以下公式六进行状态预测估计，获取当前时刻卫星的一步预测状态估计值；

基于一步预测状态量的2n+1个采样点，利用以下公式七获取当前时刻卫星的一步预测状态协方差；

其中，γ_i,k|k-1表示一步预测状态量的第i+1个采样点，表示k时刻卫星的一步预测状态估计值，P_k|k-1表示k时刻卫星的一步预测状态协方差矩阵，Q表示过程噪声协方差矩阵。

步骤S3，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差。

具体地，以当前时刻为k时刻为例，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，包括：

基于确定性采样点对称采样策略，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用以下公式八计算2n+1个确定性采样点；

利用以下公式四确定每个确定性采样点对应的加权系数；

其中，ξ_i,k|k-1,i＝0,1,...,2n表示由k时刻卫星的一步预测状态估计值和一步预测状态协方差矩阵P_k|k-1计算得到的确定性采样点，表示矩阵的第i列。

进一步地，基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差，可以包括：

基于获取的确定性采样点，利用以下公式九进行测量信息传递，获取2n+1个经测量方程传递的一步预测采样点；

χ_i,k|k-1＝h(ξ_i,k|k-1),i＝0,1,2,...,2n公式九

基于2n+1个经测量方程传递的一步预测采样点，利用以下公式十进行状态预测估计，获取卫星量测输出量的一步预测值；

基于2n+1个经测量方程传递的一步预测采样点，利用以下公式十一获取卫星量测输出量的一步预测值与一步预测状态估计值的互协方差；

其中，χ_i,k|k-1表示第i+1个经测量方程传递的一步预测采样点，表示卫星量测输出量的一步预测值，表示卫星量测输出量的一步预测值与一步预测状态估计值的互协方差，χ_1:n,k|k-1表示i从1到n的前n个经测量方程传递的一步预测采样点χ_i,k|k-1，χ_n+1:2n,k|k-1表示i从n+1到2n的后n个经测量方程传递的一步预测采样点χ_i,k|k-1。

步骤S4，基于中心误差熵准则建立卫星状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差。

具体地，基于中心误差熵准则建立卫星状态对应的线性化回归方程，包括：

定义一步预测误差为：

定义量测斜率矩阵为：

将测量向量z_k近似为：

建立卫星状态对应的线性化回归方程为：

其中，I表示单位矩阵，r_k表示高阶误差项。

进一步地，设定：

其中，S_k、S_p,kk-1和S_r,k分别表示矩阵Ω_k、P_k|k-1和R_k的Cholesky分解。

对上述公式十五所示的线性化回归方程的等式两边同乘以将线性化回归方程变形为：

d_k＝W_kx_k+e_k公式十七

其中，

进一步地，设定：e_k＝[e_1,k,e_2,k,…,e_L,k]^T，d_k＝[d_1,k,d_2,k,…,d_L,k]^T，W_k＝[w_1,k,w_2,k,…,w_L,k]^T，e_i,k＝d_i,k-w_i,kx_k(i＝1,…,L)，L＝m+n，e_i,k表示e_k的第i维状态，d_i,k表示d_k的第i维状态，w_i,k表示W_k的第i行向量；

则CEECDKF算法的代价函数为：

其中，λ表示权重系数，表示核宽为σ₁的高斯核函数，表示核宽为σ₂的高斯核函数。

设定：则公式十八可以表述为：

在中心误差熵(centered error entropy，CEE)准则下，可以通过最大化上述代价函数得到当前时刻卫星的状态的最优估计值。

具体地，以当前时刻为k时刻为例，k时刻卫星的状态的最优估计值可以通过以下公式二十确定；

其中，表示k时刻卫星的状态的最优估计值，即k时刻卫星的状态估计值，表示J_L(x_k)取得最大值时对应的x_k值。

进一步地，本发明一实施例中，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差，可以包括以下步骤：

对代价函数计算梯度，将梯度计算公式表达成矩阵形式，并令梯度等于0；

基于矩阵形式的梯度计算公式，利用定点迭代算法和矩阵求逆引理，得到当前时刻卫星的状态估计值和状态协方差。

具体地，以当前时刻为k时刻为例，可以利用以下公式二十一对代价函数计算梯度，并令梯度等于0；

其中，

将上述公式二十一所示的代价函数梯度计算公式表达成如以下公式二十二所示的矩阵形式；

其中， Ω_k表示第i行第j列元素为(Ω_k)_ij的矩阵，

基于矩阵形式的代价函数梯度计算公式，采用定点迭代算法，得到k时刻卫星的状态估计值为：

其中，表示第t次迭代时的k时刻卫星的状态估计值，表示第t-1次迭代时的k时刻卫星的状态估计值，表示第t-1次迭代时的加权矩阵M_k。

进一步地，进行如以下公式二十四至公式三十四的设定：

其中，为n×n维矩阵，为m×n维矩阵，为n×m维矩阵，为m×m维矩阵，表示矩阵的第i行第j列组成的矩阵，表示矩阵的第i行第j列组成的矩阵， 表示矩阵的第i行第j列组成的矩阵，表示第i行第j列元素为的矩阵，

进一步地，基于上述公式十七、公式二十六、以及相关参数的设定，可以将第t次迭代的k时刻卫星的状态估计值表示为：

其中，

利用矩阵求逆引理，可以将公式三十五更新为：

其中，

进一步地，基于上述设定，当前k时刻卫星的状态协方差矩阵可以更新为：

其中，表示第t次迭代时的k时刻卫星的状态协方差矩阵，I表示单位矩阵。

基于上述限定的步骤和内容，本发明一实施例中，野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法的滤波步骤的迭代求解过程可以包括以下步骤：

S201，设定k＝1，选定核宽σ₁和σ₂，选定权重系数λ，设定迭代停止条件，选定正数ε的数值，给定卫星状态估计值和状态协方差的初始值和P₀；

S202，利用上述公式六计算得到卫星的一步预测状态估计值利用上述公式七计算得到卫星的一步预测状态协方差P_k|k-1；

S203，利用上述公式十计算得到卫星量测输出量的一步预测值利用上述公式十一计算得到卫星量测输出量的一步预测值与一步预测状态估计值的互协方差通过Choleskey分解得到参数S_k、S_p,k|k-1和S_r,k，利用上述公式十七得到参数d_k、W_k和e_k；

S204，设定：t＝1，其中，表示第t次迭代时的k时刻卫星的状态估计值；

S205，利用上述公式三十六计算第t次迭代时的k时刻卫星的状态估计值利用上述公式三十七计算第t次迭代时的k时刻卫星的状态协方差矩阵

S206，比较如下判别条件：

若上述判别条件成立，则设定并执行步骤S207，若上述判别条件不成立，则设定t＝t+1并返回步骤S205；

S207，设定k＝k+1，输出和P_k，并返回步骤S202。

进一步地，基于上述设定和计算更新过程，本发明一实施例中给出了如下表示1所述野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法的滤波步骤的迭代求解过程的伪代码。

表1

本发明一实施例提供的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法通过利用Stirling插值变换进行一步预测状态估计值和一步预测状态协方差的获取，然后构造线性化回归方程，利用中心误差熵准则进行卫星的后验状态的更新求解，能够有效应对卫星姿态确定的非线性系统中出现的野值非高斯噪声，提高处理野值非高斯噪声时的卫星姿态估计精度和鲁棒性。

以下结合具体示例，对本发明一实施例提供的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法的有益效果进行说明。

某卫星姿态确定系统的系统状态方程和测量方程如公式三十八所示：

其中，q_bo为四元数表示的卫星姿态，b为陀螺仪的常值漂移，也作为状态变量被估计，ω_g为陀螺仪测量值，ω_g为输入，η_g和η_b为零均值过程噪声，η_g和η_b的协方差为Q_g和Q_b，Ω_d为传递矩阵，可以表示为：

q_opt为卫星姿态敏感器的观测输出，q_N为测量噪声四元数，ω_oi＝[0 -ω₀ 0]表示惯性系下的轨道角速度，E_bo表示从轨道系到星体系的坐标转换矩阵，可以表示为：

用于卫星姿态确定的参数设定如表2所示；

表2

假设陀螺噪声为高斯分布，过程噪声协方差设为假设星敏感器量测噪声为混合高斯噪声，

附图2-4给出了传统中心差分卡尔曼滤波算法(CDKF)、最大相关熵中心差分卡尔曼滤波算法(MCCDKF)、最小误差熵中心差分卡尔曼滤波算法(MEECDKF)和本发明一实施例的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法(CEECDKF)得到的不同姿态估计结果对比示意图。附图中，RMSE of滚转角均方根误差，RMSE ofθ：俯仰角均方根误差，RMSE ofψ：偏航角均方根误差，time：时间。

其中，各个算法的核宽参数选择如表3所示；

表3

得到的不同算法下三轴姿态角的平均的均方根误差(ARMSE)如表4所示；

表4

可以看出，本发明一实施例提供的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法在非高斯噪声条件下具有最高的滤波精度，同时滤波收敛后具有最低的估计误差协方差，即最好的滤波稳定性，能够更好地应对复杂环境下的野值噪声干扰。

需要说明的是，在本文中，诸如“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。此外，本文中“前”、“后”、“左”、“右”、“上”、“下”均以附图中表示的放置状态为参照。

最后应说明的是：以上实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (8)

1.一种野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，包括：

根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统；

根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差；

根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差；

基于中心误差熵准则建立卫星状态对应的线性化回归方程，利用线性化回归方程确定中心误差熵准则滤波的代价函数，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差；

基于获取的确定性采样点，通过非线性系统中的状态方程进行状态传递，获取当前时刻卫星的一步预测状态估计值和一步预测状态协方差，包括：

基于获取的确定性采样点，利用以下公式五进行状态传递，获取一步预测状态量的2n+1个采样点；

γ_i,kk-1＝f(ξ_i,k-1),i＝0,1,2,...,2n公式五

基于一步预测状态量的2n+1个采样点，利用以下公式六进行状态预测估计，获取当前时刻卫星的一步预测状态估计值；

基于一步预测状态量的2n+1个采样点，利用以下公式七获取当前时刻卫星的一步预测状态协方差；

γ_i,kk-1表示一步预测状态量的第i+1个采样点，ξ_i,k-1,i＝0,1,...,2n表示由k-1时刻卫星的状态估计值和状态协方差矩阵计算得到的确定性采样点，f(·)表示系统的状态方程，表示k时刻卫星的一步预测状态估计值，W_i ^m表示第i+1个采样点均值的加权系数，P_k|k-1表示k时刻卫星的一步预测状态协方差矩阵，W_i ^c1表示第i+1个采样点第1类加权系数，W_i ^c2表示第i+1个采样点第2类加权系数，Q表示过程噪声协方差矩阵；

基于获取的确定性采样点，通过非线性系统中的测量方程进行测量信息传递，获取卫星量测输出量的一步预测值、以及一步预测值与一步预测状态估计值的互协方差，包括：

基于获取的确定性采样点，利用以下公式九进行测量信息传递，获取2n+1个经测量方程传递的一步预测采样点；

χ_i,k|k-1＝h(ξ_i,k|k-1),i＝0,1,2,...,2n 公式九

基于2n+1个经测量方程传递的一步预测采样点，利用以下公式十进行状态预测估计，获取卫星量测输出量的一步预测值；

基于2n+1个经测量方程传递的一步预测采样点，利用以下公式十一获取卫星量测输出量的一步预测值与一步预测状态估计值的互协方差；

χ_i,k|k-1表示第i+1个经测量方程传递的一步预测采样点，ξ_i,k|k-1,i＝0,1,...,2n表示由k时刻卫星的一步预测状态估计值和一步预测状态协方差矩阵P_k|k-1计算得到的确定性采样点，h(·)表示系统的测量方程，表示卫星量测输出量的一步预测值，表示卫星量测输出量的一步预测值与航天器的一步预测状态估计值的互协方差阵，χ_1:n,k|k-1表示i从1到n的前n个经测量方程传递的一步预测采样点χ_i,k|k-1，χ_n+1:2n,k|k-1表示i从n+1到2n的后n个经测量方程传递的一步预测采样点χ_i,k|k-1。

2.根据权利要求1所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，根据卫星的测量数据和卫星姿态动力学模型，建立卫星姿态确定的非线性系统为：

其中，x_k表示k时刻卫星的n维状态向量，x_k-1表示k-1时刻卫星的n维状态向量，ω_k-1表示k-1时刻的n维过程噪声序列，z_k表示k时刻的m维测量向量，ν_k表示k时刻的m维测量噪声序列。

3.根据权利要求2所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，设定前一时刻为k-1时刻，当前时刻为k时刻，根据前一时刻卫星的状态估计值和状态协方差，利用Stirling插值变换计算确定性采样点，包括：

基于确定性采样点对称采样策略，根据前一时刻卫星的状态估计值和状态协方差，利用以下公式三计算2n+1个确定性采样点；

利用以下公式四确定每个确定性采样点对应的加权系数；

其中，表示k-1时刻卫星的状态估计值，P_k-1表示k-1时刻卫星的状态协方差矩阵，表示矩阵的第i列，n表示系统状态维数，h表示调节因子，表示第1个采样点均值的加权系数。

4.根据权利要求3所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用Stirling插值变换计算确定性采样点，包括：

基于确定性采样点对称采样策略，根据当前时刻卫星的一步预测状态估计值和一步预测状态协方差，利用以下公式八计算2n+1个确定性采样点；

利用以下公式四确定每个确定性采样点对应的加权系数；

其中，表示矩阵的第i列。

5.根据权利要求4所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，设定一步预测误差为：

设定量测斜率矩阵为：

将测量向量z_k近似为：

建立卫星状态对应的线性化回归方程为：

其中，I表示单位矩阵，r_k表示高阶误差项。

6.根据权利要求5所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，中心误差熵准则滤波的代价函数为：

其中，λ表示权重系数，L＝m+n，表示核宽为σ₁的高斯核函数，e_i,k表示误差变量e_k的第i维状态，表示核宽为σ₂的高斯核函数，e_j,k表示误差变量e_k的第j维状态。

7.根据权利要求6所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，利用以下公式二十确定当前时刻卫星的状态；

其中，表示k时刻卫星的状态的最优估计值，

8.根据权利要求7所述的野值噪声干扰下基于中心误差熵中心差分卡尔曼滤波的卫星姿态确定方法，其特征在于，对代价函数进行最大化处理，得到当前时刻卫星的状态估计值和状态协方差，包括以下步骤：

对代价函数计算梯度，将梯度计算公式表达成矩阵形式，并令梯度等于0；

基于矩阵形式的梯度计算公式，利用定点迭代算法和矩阵求逆引理，得到当前时刻卫星的状态估计值和状态协方差。