An Efficient Nonlinear Filter for Spacecraft Attitude Estimation

Increasing the computational efficiency of attitude estimation is a critical problem related to modern spacecraft, especially for those with limited computing resources. In this paper, a computationally efficient nonlinear attitude estimation strategy based on the vector observations is proposed. The Rodrigues parameter is chosen as the local error attitude parameter, to maintain the normalization constraint for the quaternion in the global estimator. The proposed attitude estimator is performed in four stages. First, the local attitude estimation error system is described by a polytopic linear model. Then the local error attitude estimator is designed with constant coefficients based on the robust H 2 filtering algorithm. Subsequently, the attitude predictions and the local error attitude estimations are calculated by a gyro based model and the local error attitude estimator. Finally, the attitude estimations are updated by the predicted attitude with the local error attitude estimations. Since the local error attitude estimator is with constant coefficients, it does not need to calculate thematrix inversion for the filter gainmatrix or update the Jacobianmatrixes online to obtain the local error attitude estimations. As a result, the computational complexity of the proposed attitude estimator reduces significantly. Simulation results demonstrate the efficiency of the proposed attitude estimation strategy.


Introduction
Attitude determination is a very important part for a spacecraft to achieve its designed mission.There are various methods for the spacecraft attitude determination.They can be divided into two classes: deterministic method and optimal estimation method [1].The deterministic method, the TRIAD algorithm, uses a minimal set of date and then solves three possibly nonlinear equations to obtain the attitude [1].It is simple and elegant; however, it is suboptimal and of limited use because it makes use of only two unit-vector measurements and ignores one piece of information from one of the unit vectors [2].The optimal estimation method is based on the solutions to the Wahba's problems, which obtains the optimal attitude estimation by minimizing an appropriate loss function.
There are many nonlinear estimation algorithms for the spacecraft attitude estimation since it is essentially a nonlinear problem.The most widely used algorithm for real time attitude estimation is the EKF.The EKF is recursive and easy to implement, but the accuracy can be surprisingly bad in the cases that the dynamic and the measurement models have highly nonlinearities or the system is with large process noise [3,4].The poor performance has driven several nonlinear filters for attitude estimation, among which the sigma point filters have attracted much attention, such as the Unscented Kalman Filter (UKF) [5], the Cubature Kalman Filter (CKF) [6], the Gauss-Hermite Quadrature Filter (GHQF) [7], and the Particle Filter (PF) [8].They deal with the nonlinear functions directly by choosing some points to approximate the probability density function of the nonlinear functions according to certain rules.It is generally believed that the sigma point filters are more accurate than the EKF; nonetheless, the computational cost of the sigma point filter seems high for engineering implementation [9].Even the implementation of the EKF is also computationally complex [10], because the Jacobian matrixes are required to update online which can be a very cumbersome and errorprone process, and it needs to calculate the matrix inversion for the gain matrix, resulting in heavy burden of the onboard computer especially for the systems with high dimension.
Several alternatives on the computational cost of the filter have been developed for attitude estimation.Wei used the optimal-REQUEST to estimate the attitude and the UKF to 2 International Journal of Aerospace Engineering estimate the gyro drifts [11].The computational cost of the attitude estimator was reduced by setting the state dimension to three rather than six.Fan and Kiani improved the realtime performance of the attitude estimator by minimizing the number of the required sigma points [9,12].Tang et al., Miao et al., and Choukroun et al. presented a reduced quaternion measurement model without losing information to reduce the computational complexity of the attitude estimator [6,13,14].The improved methods mentioned above have reduced the computational cost of the attitude estimator efficiently; nevertheless, the improved sigma point methods remain to have high computational cost for spacecraft which is small and with limited computational abilities.The calculation of the matrix inversion remains in the attitude estimator on the bias of reducing the dimension of the measurement model.
In this paper, an efficient nonlinear filtering method is developed for spacecraft attitude estimation.By introducing the polytopic linear differential inclusion (PLDI) theory given by Boyd et al. [15], the local attitude estimation error system is represented by an uncertain polytopic linear model.This leads to the local error attitude estimator designed with constant filter coefficients, without calculating matrix inversion or updating the Jacobian matrixes online.Thus the computational cost is sharply reduced.
The rest of the paper is organized as follows: Section 2 briefly introduces the attitude kinematics and the sensor models.Section 3 presents the implementation of the efficient nonlinear attitude estimator in detail.Section 4 demonstrates the performance of the attitude estimator and compares the results of this method with the Multiplicative EKF (MEKF) and other filters.Section 5 gives the conclusion remarks.

Attitude Kinematics and Sensor Models
In this section, the attitude kinematics and the sensor models are briefly introduced.
The spacecraft attitude can be described by various parameters, such as the direction cosine matrix, the principal axis and angle, the Euler angels, quaternion, Rodrigues parameters (RPs), and the modified Rodrigues parameters (MRPs).The most widely used attitude parameter is the quaternion because of its nonsingular character for any arbitrary rotation angle and its bilinear kinematic equation.
The quaternion is defined as where ⃗  is principal rotation axis and  is the corresponding rotation angle.
The attitude kinematic equation in the quaternion form is given by where and  is the angular velocity of the spacecraft.
The gyro is commonly used to measure the angular velocity of the spacecraft; a general gyro model is given by where   is the measurement noise and  is the drift rate bias driven by a white noise where   is assumed to be white noise.
In most practical applications, a typical attitude estimation system for the spacecraft comprises several gyros and vector sensors, such as the sun sensor, star sensor, and magnetometer.Therefore, the vector observation is chosen as the attitude measurement for the most general case.The measurement model for a signal vector observation is described as where () is the attitude matrix,    is the observation vector, and    is the known reference vector.V  is the measurement noise assumed to be white noise, and the superscript  denotes the index of the observation vector.

The Computationally Efficient Attitude Estimator
The spacecraft attitude state is given by the attitude quaternion and the gyro drift rate bias Then the dynamic model for the attitude estimation system is Denote the error attitude state as where q and b are the predicted attitude quaternion and gyro drift rate bias and q−1 is the inversion of q.The dynamic model for the attitude estimation error system is [16] International Journal of Aerospace Engineering where ω is the predicted angular velocity; the angular velocity error is expressed as follows: Substitute the above equation into (10); the dynamic model for the attitude estimation error system can be rewritten as The quaternion must obey a normalization constraint.In order to ensure that the quaternion maintains the normalization constraint, the most common attitude quaternion estimation method uses an unconstrained three-component vector to represent the local attitude error.In this paper, the RPs are chosen as the local attitude error parameter.The error RPs is defined in the terms of the error quaternion by and the inverse transformation is given by If the attitude error is quite small or tends to be zero, (13) can be approximated as Therefore, the dynamic model ( 12) can be easily approximated to the first order form by using the error RPs where The corresponding measurement equation of the attitude estimation error system is where . . .

The Multiplicative EKF for Attitude Estimation.
The most well-known nonlinear filter for spacecraft attitude estimation is the EKF.There are several different implementations of the attitude EKF, depending on both the attitude parameter used in the state vector and the form in which observations are input [3].The best known and most widely used attitude EKF is MEKF.The attitude MEKF is derived from the following equations.
Attitude predictions [5]: Predicted measurement error: Covariance matrix for the attitude predicted errors: Gain matrix: Covariance matrix for the attitude estimation errors: (24) Local error attitude estimations: Error quaternion estimation: Attitude estimations: where where   and  V are the covariance matrixes for the process noise and the vector observation noise, respectively, and Δ is the discretization step size.
International Journal of Aerospace Engineering

The Computationally Efficient Attitude Estimator.
The evaluation for the gain matrix in attitude MEKF requires to calculate the inverse of a 3 × 3 matrix, resulting in heavy computational burden when  is large, especially for small spacecraft with limited computational source.A computationally efficient attitude MEKF based on the reduced vector observation model (RMEKF) is developed to solve this problem [13,17].
According to Cayley transformation, the attitude matrix  can be mapped to a minimum-element attitude parameterization, expressed by the skew symmetric Rodrigues matrix  [18]: and the inverse transformation is expressed as follows: where  is the skew symmetric matrix generated from the RPs :  = [×].Equation ( 6) can be rewritten as Substitute the second term of (29) into the above equation; one can obtain the following equation: where Then, the observation model can be expressed as where where . ., .In order to reduce the computational cost of attitude MEKF, the dimension of the observation model equation ( 33) can be reduced to 3 by multiplying both sides of the equation by   1 [13,17].The weighted factor is firstly designed on the bias of the information conservation principle, to ensure the information for each vector observation without losing after the dimension of the observation model reduced to 3. The weighted factor is designed as where . Then, (33) can be rewritten in the following form: where . . .
Multiply both sides of the above equation by   2 , one can get where The local attitude estimation error system is composed of ( 16) and (38).It is obvious that the dimension of the vector observation model described by the above equation is 3.It only requires evaluating a 3 × 3 matrix inversion for the gain matrix   in the MEKF, rather than the 3 × 3 matrix inversion.As a result, the computational burden is reduced.
According to the MEKF, it is easy to get the following implementation of the RMEKF [13,17].The equations for attitude predictions, error attitude estimations, and estimations are the same in the RMEKF and MEKF.For the above reason, only the equations for local error attitude estimations in the RMEKF are shown here.

The Improved Computationally Efficient Attitude Estimator.
A 3-dimension vector observation model is used in attitude RMEKF instead of the original 3-dimension model.
Only the calculation of a 3 × 3 matrix inversion is required for the gain matrix, resulting in much less computational cost than that of attitude MEKF.However, the Jacobian matrixes need to update online in attitude RMEKF too, and the process of the matrix inversion remains in the calculation for the gain matrix.In order to further reduce the computational burden of attitude EKF, a new attitude estimation strategy is developed in this section.
The local attitude error estimation system composed of ( 16) and ( 18) is rewritten as follows: where It is assumed that the process noise   and the measurement noise V are uncorrelated white noise.The parameter  is bounded by the 6-dimension space and the value set of  belongs to a compact set, since the angular velocity varies in a finite interval in most practical applications and each element of attitude quaternion takes value on the interval [0, 1].That is to say,  ∈ Ω ⊂ R 6 and Ω is a compact set. Denote Then, the above equation can be approximated as a convex combination of the  constant linear system matrixes   ,  = 1, . . .,  [19], namely, where   () is the bias of the convex combination, which satisfies ∑  =1   () = 1, 0 ≤   () ≤ 1.Thus, the local attitude estimation error system equation (45) can be described with the following form: where the system matrixes , , , and  are denoted as (48).The equivalent discrete-time form of (49) can be approximated as where  , is the equivalent discrete-time noise The state transition matrix can be approximated by Since   () satisfies ∑  =1   () = 1, 0 ≤   () ≤ 1, the above equation can be easily reexpressed as For convenience of notation, let then (53) can be rewritten as It is now a straightforward matter to show that the discrete-time form of the local attitude error estimation system can be described by an uncertain discrete linear polytopic model in the form of (50).With this model, the local error attitude estimation problem is converted to a robust linear one, that is to find a stable local error attitude estimator in the form such that the local attitude estimation error variance, [(Δ  − Δx  )  (Δ  − Δx  )], is minimized, where (  ,   ,   ,   ) are constant matrixes to be determined.According to the robust  2 filtering algorithm given by Duan et al. [19], the matrixes (  ,   ,   ,   ) can be obtained by solving an optimization given in the following Lemma.
Lemma 1 (see [19]).Consider the system (50); a filter of the form (56) that achieves a suboptimal guaranteed filtering error covariance bound can be derived from the following optimization: where  1 and  2 are fixed parameters, and The suboptimal filter is given by As a consequence, the improved computationally efficient attitude estimator can be presented as shown in Table 1.
The filter coefficients (  ,   ,   ,   ) for the local error attitude estimator can be solved before the recursive process of the attitude estimation.Since they are constant, the matrix inversion is not required for the gain matrix and the updating of the Jacobian matrixes online neither.As a result, the computational complexity of the proposed efficient attitude estimator (PEF) reduces extensively.As can be seen in the above attitude estimation strategy, the key to implement the PEF is to determine the vertexes of the polytopic linear model (49).There are two kinds of methods: parameter boundary method and TP model transformation method [20].The polytopic linear model for an affine parameter system obtained by the former is quite accurate, but, for other systems, the polytopic linear model obtained by the latter is with less conservativeness.The TP model transformation algorithm is described in Algorithm 1.

Simulation Results
In this section, the comparisons of the computational cost and the accuracy between the PEF and other attitude estimators are given.

Computational Complexity.
The computer complexity of the attitudes PEF, RMEKF, and MEKF is shown in Table 2.Only the computational cost of the local error attitude estimation is given since the equations of the attitude predictions and updatings are the same in the three attitude filters.
According to the statistics, it is obvious to find that the computational cost of the MEKF is much more than that of the other two filters.The additional cost of the PEF is less than that of the RMEKF if  < 28, while the multiplicational cost of the RMEKF is much more than that of the PEF if  ≤ 70.It happens that in a realistic situation the observation vector Step 1. Determine the bounded domain Ω where the parameter  varying on; Step 2. Set the transformation space Ω into  uniformly distributed grid: Step 3. Sample the functions () over the hyper rectangular grid and store the sample matrixes in the tensor F; Step 4. Execute the higher order singular value decomposition (HOSVD) on tensor F and extract the minimal basis, the result of this step is F ≈ V ⊗    ; Step 5. Normalize the basis matrix   , one can obtain () ≈ V ⊗    ; Step 6. Extract the vertexes of the polytopic linear model (49) from V.
Algorithm 1: TP model transformation algorithm [20].number is less than 28; therefore, one can get the conclusion that the PEF has the least computational cost in the three attitude filters and the MEKF has the most.
Denote the computational cost error between the PEF and other attitude filters as where  PEF and   are the computational cost of the multiplication or addition for the PEF and other filters (MEKF or RMEKF), respectively.The computational cost error for the MEKF and PEF is shown in Figure 1.The increasing rates of the multiplicational cost error and the additional cost error are nearly the same for the two filters.They both increase rapidly as the observation vector number increases, indicating that the computational complexity of the PEF is much less than the MEKF, especially when the observation vector number is large.
The computational cost errors for the RMEKF and PEF are shown in Figure 2. The advantage of less computation complexity is not evident as the observation vector number increases in the PEF, while the decrescent rate of the multiplicational cost error for the two filters reduces more slowly than that of the additional cost error.Since the multiplication is much more complex than the addition, it can be concluded that the computational burden of the PEF reduces much more than the other two filters for the application of spacecraft attitude estimation in the practical engineering.The initial gyro drift rate is 5 deg/h and the standard covariance of the driven noise is 0.003 deg/ √ h.The standard covariances of the gyro measurement noise and the vector observation noise are 0.5 deg/h and 5 arcsec, respectively.The initial covariance of attitude estimation error is 10 −8  6 and the discretization step size is 200 ms.The initial attitude estimations are set as their true avoiding the poor performance or divergence of the three estimators by lacking a good a priori estimate of attitude.The diagrammatic representation of the attitude estimation system (attitude PEF or the other two estimators) designed in the simulation is shown in Figure 3; the attitude estimations are given by the attitude estimator (attitude PEF or the other attitude estimators) based on the gyro and the vector observations.The observation vectors can be given by unit sun, star, and Earth's magnetic field vectors.Two unit star vectors are chosen in the simulation, because several star vector observations can be obtained from the star tracker at a time.The Monte Carlo simulation is computed over an ensemble of 100 independent runs.In order to make the simulation results more intuitive, the attitude estimation errors are shown by the error attitude 3-1-2 Euler angles instead of quaternion, because the quaternion does not have intuitive physical meanings.
The average estimation errors of the attitude and the gyro drift rate bias with their respective 3 bounds in the PEF are shown in Figure 4.As can be seen in the figure, the estimation errors of the attitude and gyro drift rate bias in the PEF all converge to within their respective 3, indicating that the PEF performs well for the attitude estimation.
The attitude principal rotation angle is used to scale the attitude estimation error, which is expressed in the form of quaternion as follows: The average initial attitude estimation errors of the three filters are shown in Figure 5.The maximal estimation error of the attitude angle is the smallest in the PEF during the transient process of the filters.The initial attitude estimation errors in the RMEKF and MEKF are nearly identical.
The steady attitude estimation errors of attitudes PEF and MEKF are given in Figure 6.The steady attitude estimation error in attitude RMEKF is not shown here, because it is nearly the same as that in attitude MEKF.From the figure, one can get that the steady attitude error angle in attitude PEF varies around 4.13 arcsec, while it varies around 1.41 arcsec in attitude MEKF.The magnitude order of the accuracy achieved by attitude PEF is 1 arcsec and it is identical to that of attitude MEKF.The maximum values of the absolute steady estimation errors for the triaxial angular velocities and gyro drift rate bias are shown in Table 3.The magnitude order of the steady estimation errors for the triaxial angular velocities and gyro drift rate bias are 10 −4 deg/s and 10 −6 deg/s in attitude PEF, respectively.They are nearly identical to those of attitude MEKF.Therefore, one can conclude that the proposed attitude PEF performs well for the spacecraft attitude estimation.Since its computational cost is much lower than that of the MEKF and RMEKF, it should be a good choice for the spacecraft attitude estimation application with limited computing resources and low accuracy demand.

Conclusion
This paper develops a computationally efficient attitude estimation method for the spacecraft based on the gyro measurements and the vector observations.The global attitude is described by the quaternion, while the local attitude error is represented by the Rodrigues parameter to ensure that the quaternion satisfies the normalization constraint.Rather than reducing the computational complexity of attitude estimator by reducing the number of the sigma points or the measurement model dimension, the local error attitude estimator is designed with constant coefficients in the proposed attitude estimator.It does not need to calculate the matrix inversion for gain matrix or update the Jacobian matrixes online during  the recursive process of the attitude estimator.As a result, the computational cost reduces extensively and the convergence speed for the attitude is much faster than other attitude filters.Simulation results show that the proposed attitude estimator improves the computational efficiency, though the accuracy for the attitude is a little larger than that of MEKF.It should be preferred over the MEKF in the practical spacecraft attitude estimation application with limited computing resources and low accuracy demand.

𝑞:
A t t i t u d e q u a t e r n i o n ⃗ : The The space of 6 dimensional vectors with real components.

Table 2 :
The statistics of the computational cost for three filters.
a The matrix inversions in MEKF and RMEKF are calculated by QR decomposition method.

Table 3 :
The maximum absolute steady estimation errors.
,   : Known reference vector and the observation vector , , ,  : Error quaternion, error angular velocity, error gyro drift rate, and error Rodrigues parameters , : Jacobian matrixes International Journal of Aerospace Engineering : A t t i t u d em a t r i x   ,   : Gyro drift driven noise and gyro measurement noise   , V: Process and measurement noise  , : The equivalent discrete-time noise Δ: Discretization step size q, ω, b: Estimations of quaternion, angular velocity, and gyro drift rate   : G a i nm a t r i xi nt h eM E K F  /−1 ,   : Covariance matrixes for predicted and estimated attitude error   ,  V : Covariance matrixes for the process and measurement noise    ,   : Filtercoefficientsofthelocalerrorattitude estimator   ,   : Filtercoefficientsofthelocalerrorattitude estimator   ,   ,   : Matrix variables in the optimization   ,  11 ,  21 : Matrix variables in the optimization  2 ,  11 ,  21 : Matrix variables in the optimization  11 ,  12 ,  22 : Matrix variables in the optimization R 6 :