Get Our e-AlertsSubmit Manuscript
Space: Science & Technology / 2022 / Article

Research Article | Open Access

Volume 2022 |Article ID 9854601 |

Baojian Yang, Hao Huang, Lu Cao, "Centered Error Entropy-Based Sigma-Point Kalman Filter for Spacecraft State Estimation with Non-Gaussian Noise", Space: Science & Technology, vol. 2022, Article ID 9854601, 13 pages, 2022.

Centered Error Entropy-Based Sigma-Point Kalman Filter for Spacecraft State Estimation with Non-Gaussian Noise

Received28 Mar 2022
Accepted17 Jul 2022
Published29 Jul 2022


The classical sigma-point Kalman filter (SPKF) is widely used in a spacecraft state estimation area with the Gaussian white noise hypothesis. The actual sensor noise is often disturbed by outliers in the harsh space environment, and the SPKF algorithm will reduce the filtering accuracy or even diverge. In this study, to enhance the robustness under non-Gaussian noise condition, the outlier-robust SPKF algorithm based on a centered error entropy (CEE) criterion is derived. Unscented Kalman filter (UKF) is typical of SPKF; combining the deterministic sampling criterion with the centered error entropy criterion, a robust centered error entropy UKF (CEEUKF) algorithm is proposed. The CEEUKF uses the unscented transformation (UT) method to perform time update step and then uses the robust regression model and CEE criterion to reconstruct the measurement update step. The effectiveness of the proposed CEEUKF is verified by a spacecraft attitude determination system.

1. Introduction

The star sensor is a high-precision attitude sensor. The spacecraft’s inertial attitude can be directly determined by using the measurement data of the star’s direction from the star sensor fixed to the satellite body [13]. The gyroscope has complementary characteristics with the star sensor, its dynamic response is fast, the random measurement error of the high-precision gyroscope is small, and its drift can be calibrated by the star sensor. Therefore, the combined use of star sensors and gyroscopes is the guarantee for high-precision attitude determination [4]. Nanosatellites often use low-precision sensors, which are easily disturbed by harsh environments. So, it is of great significance to study high-precision robust state estimation algorithms [57].

A spacecraft attitude kinematics model, attitude measurement model, and filter algorithm are three important parts in spacecraft attitude determination, and a high-precision filtering algorithm is the key to attitude determination. The Kalman filter (KF) is the optimal filter with the linear Gaussian framework. Once it came out, it has been widely used in various fields of industry [8, 9]. However, actual systems are often nonlinear systems, and there is no optimal filtering algorithm for nonlinear systems. Only approximate methods can be used for the nonlinear Gaussian systems. For nonlinear filtering problems, approximating the nonlinear model to a linear model is a widely used method at present. The extended Kalman filter (EKF) adopted this idea, which linearized the nonlinear model and used the traditional Kalman filter to estimate the state [912]. The nonlinear filtering algorithm based on a deterministic sampling criterion can approximately solve the a posteriori integral under Gaussian hypothesis, so as to determine the first two moments of the a posteriori state. Compared with the linearization of nonlinear function, a deterministic sampling method has high precision and moderate amount of calculation. It is a powerful tool to solve the problem of nonlinear filtering. Unscented transformation, cubature transformation, and Stirling interpolation formula are three widely used deterministic sampling methods. The resulting nonlinear Gaussian filtering methods are called unscented Kalman filter (UKF), cubature Kalman filter (CKF), and central differential Kalman filter (CDKF). Since these methods involve the sampling of deterministic points, they are called SPKF methods in this paper. In the Gaussian system, these sampling methods can approach the first three moments of the state and have similar filtering accuracy [1318].

Although the SPKF algorithm performs well in ideal Gaussian white noise, the actual operating conditions of the spacecraft in orbit are complicated. Space environmental interference, solar panel jitter, and flicker noise will make the noise no longer meet the Gaussian distribution and present a heavy-tailed non-Gaussian situation [1921]. In the nonlinear non-Gaussian system, the classical SPKF filtering method is no longer applicable, and there will be obvious accuracy degradation or even filtering divergence. The essential reason is that the classical nonlinear filtering algorithms are all based on the minimum mean square error (MMSE) criterion to calculate the mean and covariance of a posteriori state, while the MMSE criterion is sensitive to outliers and is not robust to non-Gaussian noise. Robust algorithms can effectively deal with non-Gaussian noise situations [22, 23].

A particle filter (PF) has no restrictions on the characteristics and noise conditions of the system and can deal with non-Gaussian noise. Its filtering accuracy is positively correlated with the number of particles. In the high-dimensional model, the number of particles exponentially increasing with the dimension will lead to the rapid increase in the computational complexity of the PF algorithm, which makes it difficult to carry out practical engineering application [24]. A Gaussian sum filter (GSF) models the system process noise and measurement noise in the form of finite weighted sum of Gaussian distribution and then runs multiple Gaussian Kalman filters to calculate the weighted sum, so as to approximate the posterior probability density of the system and obtain the posterior state estimation [25]. When applying the GSF method, the probability density of noise is required to be known, and the increase in Gaussian summation term will lead to the exponential increase in calculation [26]. Student’s -based Kalman filter models non-Gaussian noise as special Student’s distribution [26]; it is not suitable for general non-Gaussian noise. The robust estimation method based on generalized maximum likelihood estimation (M-estimation) is also an effective way to deal with non-Gaussian noise [27]. The M-estimation method establishes the regression model by combining the system equation and the measurement equation and then obtains the final state estimation value based on the influence function. This method assigns different weights to different sample points, so it can reduce the influence of outliers on the estimation value. The robust Huber filtering (HF) algorithm is a classical M-estimation method based on Huber’s cost function, and it is a combination of minimum and norm [28, 29]. The accuracy of the HF algorithm is limited because the Huber cost function cannot be reduced back. Filters based on the entropy criterion in information theory learning are also M-estimation methods. The maximum correntropy (MC) criterion-based Kalman filter (MCKF) suitable for a linear non-Gaussian system is developed first [30], and then, the nonlinear non-Gaussian filters such as maximum correntropy unscented Kalman filter (MCUKF) and maximum correntropy square-root cubature Kalman filter (MCSCKF) are derived [3133]. The minimum error entropy (MEE) criterion may be more robust than MC. Several filters based on the MEE criterion, such as MEEKF and MEEUKF, can effectively deal with heavy-tailed non-Gaussian noise and multimodal distribution noises [34].

However, the existing MEE-based filters are biased estimators because the MEE criterion is shift-invariant [35]. Besides, the MEEKF may have numerical instability [36, 37]. Recently, the MEEKF algorithm has been corrected and the robust MEEKF (R-MEEKF) algorithm is formed [36], but this cannot change the characteristics of the MEE. The CEE criterion has been verified to be more robust than MEE and MC criteria [38, 39]. The CEEKF for the linear non-Gaussian system has been developed [40]. In this paper, we are committed to extending this algorithm to nonlinear and non-Gaussian fields. By referring to the research idea of nonlinear filtering, this paper combines the deterministic sampling criterion with the CEE criterion and develops several sigma-point filtering methods based on the CEE criterion. Among our study, the sigma-point sampling methods are used to perform time update step. Then, the posterior state is redesigned by CEE. Because the higher-order information of the error is captured by the CEE criterion, CEESPKFs should be more robust to deal with non-Gaussian noise. Because several sampling methods have similar accuracy, the typical UKF algorithm is used for algorithm development and numerical simulation. The application to a tracking problem and spacecraft attitude determination system verifies this theory. The simulation results show that the proposed filter has higher accuracy than the classical UKF and existing robust MCUKF and MEEUKF algorithms.

The outline of this paper is as follows. The classical UKF algorithm and CEE criterion are introduced in Section 2. The CEEUKF algorithm is derived in Section 3. The simulation example is shown in Section 4. The conclusion is introduced in Section 5.

2.1. Centered Error Entropy Criterion
2.1.1. Error Entropy Criterion

The quadratic Renyi’s entropy of error is expressed as where the exact probability density function (PDF) can be approximated by samples [41]: where is the Gaussian kernel function and is the kernel size. So, Renyi’s entropy can be approximated as

The information potential (IP) can be denoted by

So, minimizing error entropy means maximizing the information potential.

2.1.2. Correntropy Criterion

For variables and , the correntropy has the following definition [30]: where is the Gaussian kernel function and is the kernel size; is the Gaussian function. The samples are used to approximate the correntropy.

2.1.3. Centered Error Entropy Criterion

We take the weighted combination of MC and MEE as the expression of CEE [40, 42], and the CEE criterion is shown as where is the weight constant. The CEE becomes the MEE when which degenerates to the MC when . and are the kernel sizes with the same definition as (4) and (6).

2.2. Sigma-Point Kalman Filter

The nonlinear discrete time system is where and are the state and measurement vectors, respectively. and are the progress function and measurement function, respectively. is the system process noise, and is the sensor measurement noise. Noises have zero mean, and their statistical characteristics are where and are process noise and measurement noise covariances, respectively.

According to Bayesian filtering theory, the filtering process includes two steps: time update and measurement update. The time update step contains computing the prediction PDF, and the measurement update step contains computing the posterior PDF [43]. When the prior PDF is known, the state one-step prediction PDF can be expressed as where represents the measurement value at the previous times. When the observation at time is known, the state posterior PDF can be recursively transferred as where is the likelihood function. The posterior PDFs in equations (11) and (12) are only ideal solution, which cannot be solved accurately in practical application.

The Gaussian assumption is often used to study the implementation of Bayesian filtering theory. The Gaussian distribution can be perfectly described by the mean and variance of the system state estimation. Therefore, equations (11) and (12) can be expressed as follows: where represents the Gaussian distribution with mean and variance . The Gaussian approximate filter contains time update step and measurement update step [43]:

The time update step contains computing the prediction state mean and error covariance:

The measurement update step contains computing the posterior state mean and error covariance: where

The above five Gaussian weighted integrals can be written in a unified form: where is an arbitrary integral function and is the integral region, the analytical solution of the integral cannot be obtained, and the approximate solution is usually obtained by the quadrature criterion. where and are the integration points and corresponding weights of the Gaussian density, respectively, and is the number of integration points.

Different sigma-point Kalman filters are obtained by using different approximate integration methods. UKF, CKF, and CDKF are all sigma-point Kalman filters. Next, the typical UT method is used in this paper, and the UKF is reviewed.

2.3. Classical UKF

Classical UKF uses the UT method to get sampling points and approximate the state mean and error covariance of a PDF [13]. The UKF method considers that it is easier to approximate PDF than a nonlinear function. Time update and measurement steps are also contained in it.

Time update step: firstly, the sigma points are generated by and . where is a positive scaling factor which can be chosen as . When the nonlinearity of the system is severe, is usually taken as a small positive value, such as . is also a scaling parameter, and it is usually set as . represents the i-th column of the matrix .

Then, the one-step prediction sigma-points are got by , and the one-step prediction values and are given by where where the parameter is a nonnegative weight function used to describe the prior information of the state. For Gaussian distribution, it is usually chosen as . and are weighting factors of mean and covariance respectively.

Measurement update step: similarly, use the one-step prediction values and to get the sigma points , and then, propagate the sigma points to through the measurement function. Then, the output prediction and cross-covariance matrix are obtained:

When getting the measurement value , the posterior state mean and covariance can be got by


In this section, the centered error entropy-based UKF (CEEUKF) is derived by the CEE criterion. The CEEUKF also has time and measurement update steps. For the nonlinear system (8) and (9), the time update of the CEEUKF algorithm is the same as the classical UKF algorithm, and it is calculated by equations (21) and (22). The new measurement update step is designed based on the CEE criterion. It includes two steps: the establishment of an augmented model and the posterior state estimation.

3.1. Establishment of the Augmented Model

The approximate measurement function based on statistical linearization [33] can be expressed as where can be calculated by

The one-step prediction error is

As stated by (29) and (31), the regression model is established as where , with where , , and ; “chol” represents the Cholesky decomposition operation. Multiply both sides of equation (33) by , and the regression equation can be got by where

The error is used in solving optimal states, and , .

3.2. Posterior State Estimation by the CEE Criterion

Similar to paper [40], the optimal state estimation can be got by maximizing the CEE criterion: where and .

Calculate the gradient of the CEE criterion and set it to zero:

By setting and , the above equation is rewritten as where

Since and are coupled, the fixed-point iteration [30] is used to get the exact solution of equation (40): where represents the -th iteration and is written as with

Using equations (26), (27), and (38), equation (37) can also be written as where

The recursive solution of state in equation (44) can be got by the matrix inversion lemma: where

The posterior covariance is

The CEEUKF algorithm is summarized in Algorithm 1.

Inputs:, , , ,
1. ; choose small , kernel sizes , and weight ; set initial value and ;
Time update
2. Use Eq. (21) and (22) to get and ;
Measurement update
3. Calculate and by using Eq. (25) and Eq. (26); calculate by using Eq. (30); Get , and from the Cholesky decomposition of ; Get and by using Eq. (35) and Eq. (36);
Fixed point iteration
4. Initialization: ,
 Get , , , and from Eq. (43)
6. if , set , and go to step 7. If not, let
, and return to step 5;
7. and return to step 2.
Outputs: ,

Next, two theorems are given to show the properties of the CEEUKF algorithm with special parameters.

Theorem 1. When , CEEUKF degenerates to MCUKF. When , CEEUKF degenerates to MEEUKF. This can be easily got according to the characteristics of CEE.

Theorem 2. When and , CEEUKF becomes the classical UKF.

Proof of Theorem 3. On the basis of Theorem 1, when , CEEUKF becomes MCUKF. So, if , then and the state is ; this value equals to the estimated value of classical UKF.

Remark 4. Three parameters , , and affect the filtering effect. By adjusting them, the proposed CEEUKF can have a better effect on processing non-Gaussian noise. According to our experience, when the parameter is set at [0.8, 1), the CEEUKF algorithm will have better filtering results; is usually used; then, and can be tried by users in different applications.

Remark 5. The convergence analysis of the proposed CEEUKF algorithm is similar to paper [40], and it will not be shown here.

4. Simulation

4.1. Attitude Estimation System

Gyro and star sensor are used for the spacecraft attitude estimation, and the orbital frame is chosen as the reference frame. The quaternion is defined by , with . The quaternion satisfies . In this section, the gyro model, attitude determination system model, and measurement model are introduced [14, 44]. In addition, the steps of attitude determination by using UKF are introduced briefly.

4.1.1. Gyro Model

The mathematical model of gyro is usually modeled as [45] where is the measurement output of gyro, is the angular velocity of the satellite relative to the inertial space, is the constant drift of the gyro, is the measurement noise of gyro, and is the gyro constant drift random noise. and are uncorrelated white Gaussian noise. where is the Dirac delta function.

4.1.2. System Model

In this paper, the combination of the star sensor and gyro is used as the model of attitude determination, and the attitude quaternion and gyro constant drift are selected as state vectors [12]. The system model for attitude determination consists of an attitude kinematics equation and gyro model: where represents the quaternion vector of the orbital coordinate frame rotated to the satellite body frame. is the transfer matrix and is defined as where is the cross-product matrix, and it is defined as

is the angular velocity vector of the satellite relative to the orbital frame, and it is calculated by where is the orbital angular velocity and is the conversion matrix from the orbit frame to the satellite body frame [46] and is expressed as

To facilitate the computer simulation, the discrete model of the system equation is given as where is the attitude quaternion at time , is the gyro constant drift at time , the is the angle velocity at time , and is the sampling interval of the gyro.

4.1.3. Measurement Model

The attitude quaternion is selected as the measurement output vector of the star sensor, and the measurement equation of the star sensor is [45] where is the measurement output of the star sensor, is the measurement noise quaternion, and represents the quaternion multiplication.

4.1.4. Unscented Filtering for Attitude Estimation

In this section, the core steps of attitude determination by using the UKF algorithm are given; the different robust UKF algorithms have similar steps with UKF. When the quaternion is directly used to determine the attitude, there will be a problem of normalization of the weighted operation [47]. For this reason, modified Rodrigues parameters (MRPs) are introduced, and the attitude determination process is completed by switching between MRPs and quaternions [14, 45]. The generalized Rodrigues parameters (GRPs) are defined as follows: where are are scalars, is the GRPs, is the vector part of the quaternion, and is the scalar part of the quaternion.

The quaternion can be represented by GRPs as follows: When and , the GRPs are the MRPs.

In UKF attitude determination, the MRPs are chosen as states and they are converted to quaternion for state transfer. When calculating the mean and error covariance, the quaternion is converted to MRPs, which avoids the problem of quaternion weighted calculation and covariance singularity.

The specific UKF attitude determination steps are given. The form of state variables with MRPs as parameters is , where is the attitude with MRP form. (a)Initialization(b)Time update

Use equations (19) and (20) to get sigma points ; in the time update step, quaternion updating is involved. The sigma points are divided into MRP part and gyro part:

Convert MRP sigma points to quaternion sigma points: where

The one-step prediction value of quaternion can be got by the first line of equation (63): where and is the estimation of gyro drift at time .

Next, the quaternions are converted to MRPs:

The time update of gyro drift is

So, the predicted state mean and variance are calculated by equations (21) and (22). (c)Measurement updateThe sigma points in the measurement step are generated by the predicted mean and covariance. These sigma points are still divided into MRP part and non-MRP part .

Only the attitude quaternion is contained in the measurement equation, and the MRP sigma points are converted to the quaternion type: where