An Improved Invariant Kalman Filter for Lie Groups Attitude Dynamics with Heavy-Tailed Process Noise

: Attitude estimation is a basic task for most spacecraft missions in aerospace engineering and many Kalman type attitude estimators have been applied to the guidance and navigation of spacecraft systems. By building the attitude dynamics on matrix Lie groups, the invariant Kalman ﬁlter (IKF) was developed according to the invariance properties of symmetry groups. However, the Gaussian noise assumption of Kalman theory may be violated when a spacecraft maneuvers severely and the process noise might be heavy-tailed, which is prone to degrade IKF’s performance for attitude estimation. To address the attitude estimation problem with heavy-tailed process noise, this paper proposes a hierarchical Gaussian state-space model for invariant Kalman ﬁltering: The probability density function of state prediction is deﬁned based on student’s t distribution, while the conjugate prior distributions of the scale matrix and degrees of freedom (dofs) parameter are respectively formulated as the inverse Wishart and Gamma distribution. For the constructed hierarchical Gaussian attitude estimation state-space model, the Lie groups rotation matrix of spacecraft attitude is inferred together with the scale matrix and dof parameter using the variational Bayesian iteration. Numerical simulation results illustrate that the proposed approach can signiﬁcantly improve the ﬁltering robustness of invariant Kalman ﬁlter for Lie groups spacecraft attitude estimation problems with heavy-tailed process uncertainty.


Introduction
The navigation and control operations in aerospace engineering usually require adequate accurate knowledge of spacecraft orientation in space, and attitude determination from vector observations is usually employed in astronautically applications [1][2][3][4]. The widely used Kalman type filters can infer the state in Euclidean vector space by fusing attitude dynamic and observation sensor according to their probabilities [5][6][7][8]. Recently, the building of attitude dynamics on matrix Lie groups has been investigated actively for the navigation and control of spacecraft targets because it can significantly improve the estimation and control performance to make full use of the geometrical properties of Lie groups models [9,10].
For attitude estimation, there have already been significant research interests in the estimation and control of spacecraft targets such as the invariant observers and filters for attitude dynamics on the special orthogonal group SO (3), or the navigation systems defined on the special Euclidean group SE(3) [11,12]. Accounting for the geometrical invariance of symmetry attitude dynamics, the resulted attitude description with rotation matrix on SO(3) can avoid the trouble of singularities and unwinding encountered when using unit Machines 2021, 9,182 2 of 19 quaternion [13] and, to some extent, also contribute to better estimation performance for attitude determination algorithms [14]. Exploiting the mapping relation between matrix Lie groups and Lie algebra, invariant Kalman filter (IKF) can convert the attitude dynamics on SO(3) to Euclidean vector space and mimic classical Kalman filtering steps [15]. IKF is based on the invariance and log-linearity of symmetry dynamics, which then leads to better convergence on a much bigger set of state trajectories [16]. In these available researches, the process noise of attitude dynamical model was assumed to be a concentrated Gaussian distribution [17,18] and the statistics parameter are also assumed to be accurately given in advance as required in classical Kalman theory [19,20]. However, due to the presence of server maneuvers in space missions, there might be some unpredictable disturbances or outliers that would contaminate the Gaussian distribution and violate the required assumptions for Kalman filtering on matrix Lie groups.
For Euclidean space filtering problems, some adaptive methods are commonly used to deal with the trouble of inaccurate/unknown noise parameters [21][22][23]. Some adaptive techniques have been designed to estimate online unknown/inaccurate noise parameters, mainly including innovation-based adaptive methods and multiple model methods [24][25][26]; nevertheless, the convergence of innovation-based estimation methods cannot be guaranteed while the parallel-running multiple filter banks might increase the computational complexities. The variational Bayesian(VB) iteration-based adaptive methods have also been employed to jointly determine the system state and inaccurate noise statistics [27][28][29]. However, these adaptive methods all assume the noise distribution to be Gaussian and cannot be applied to the filtering cases with heavy-tailed noise. For heavy-tailed system noises, the Student's t distribution assumes the distribution of the latent state and system noises as a joint Student's t-distribution, and it has been proved to be more suitable to represent the non-Gaussian probability density function (PDF) [30,31].
To deal with the heavy-tailed process noise caused by severe maneuvering of rigid bodies, a robust Student's t-based hierarchical Gaussian system model for the SINS/GPS integration is presented by modeling the prior probability density function as a Student's t distribution and the conjugate prior distributions [21]. To improve the filtering performance of Student's t-based methods for cases with both heavy-tailed process and observation noises, a robust Student's t-based Kalman filter is proposed employing the variational Bayesian-based iterative approach in [32]; in the work of [33], by applying the Student's t approximate distributions to posterior filtering and smoothing probability density functions, the robust nonlinear Student's t-based filter and smoother are also presented as the generalization and extension of the linear Student's t filtering methods. Note that the above researches provide a sound theoretical foundation for attitude estimation problems with heavy-tailed noises, but they are originally designed in Euclidean space systems for filtering of vector variable states. To the best of our knowledge, for attitude estimation problem in SO (3), until now there is still no special study on the robust Student's t-based adaptive methods for invariant Kalman filtering on matrix Lie groups.
Note that although the available invariant Kalman filter also obeys the basic theory of classical Kalman filter, the above researches about robust Student's t-based methods cannot be directly applied to matrix Lie groups systems. The uncertainty definition of concentrated Gaussian distribution on matrix Lie groups requires the assumption of noisy errors being small enough [34,35], which is the condition for first order approximation to the Baker-Campbell-Hausdorff (BCH) formula [36,37]. For attitude estimation on SO(3) in aerospace engineering, the presence of heavy-tailed noises and disturbances imposes a critical constraint on the application of robust Student's t-based methods to improve the performance of invariant Kalman filter.
According to the above discussions, it is rather meaningful to investigate the robust Student's t-based invariant Kalman filtering method especially for the matrix Lie groups attitude estimation problems in the presence of heavy-tailed outliers in system process noise. In this work, for attitude estimation problem with heavy-tailed process noise, this paper proposes a hierarchical Gaussian state-space model for invariant Kalman filtering, and the Lie groups rotation matrix of spacecraft attitude is inferred together with the scale matrix and dof parameter using the variational Bayesian iteration. Numerical simulation results illustrate that the proposed approach can significantly improve the filtering robustness of invariant Kalman filter for Lie groups spacecraft attitude estimation problems with heavy-tailed process uncertainty. Attitude determinition from two vector observations is an important task in most aerospace engineering. The attitude dynamic on special orthogonal group SO(3) ⊂ R 3×3 with the associated observation vectors can be formulated as the following system model [1,14]

Primaries and
where the special orthogonal group SO(3) := RR T = I 3×3 , det(R) = 1 is the group of 3 × 3 rotation matrices and R k ∈ SO(3) is the mapping rotation from the body frame of the spacecraft to the earth-fixed frame at time instant k; Ω k ∈ SO(3) denotes the control rotation; Y k ∈ R 6 is a composition of the vector observations y k , y k ∈ R 3 with their respective Gaussian white and mutually independent observation noises with N c denoting the Gaussian distribution. Note that the vector observation based attitude determination is a common application scenario because the Sun sensor, magnetometer, horizon sensor, and the gravitational accelerometer are usually equipped by most spacecrafts in astronautical engineering [1]. A graphical presentation of target application scenario is given in Figure 1. Note that the additive forms of vector noises in Euclidean space are not supported in Lie Groups space. The vector process noise w k ∈ R 3 is injected through Lie exponential operator Exp G (·) and, for any w = [w 1 w 2 w 3 ] T ∈ R 3 , the projection between the vector noise w = [w 1 w 2 w 3 ] T ∈ R 3 and the associated groups element R ∈ SO(3) are where exp(·) is the matrix exponential operator and w × is the skew-symmetric matrix. The concept of concentrated Gaussian distribution (CGD) was employed to define the probability distribution for random variables on matrix Lie groups [19,20,37]. For attitude estimation, in concentrated Gaussian distribution the groups variable R ∈ SO(3) with mean R ∈ SO(3) and covariance Σ ∈ R 3×3 is defined as where G(·, ·) denotes the concentrated Gaussian distribution (CGD). Note that, in the definition of CGD, the eigenvalues of Σ are assumed to be small enough, i.e., almost all the mass of the distribution should be concentrated in a small neighborhood around the mean of the rotation variable [20,37]. A graphical illustration of the relation between the concentrated Gaussian distribution on matrix Lie groups SO(3) and the classical Gaussian distribution relation is presented in Figure 2.

Model Projection Based on the Invariance Property of Attitude Estimation System
For attitude estimation on SO(3), the usual Kalman filtering methods in Euclidean vector space could not be applied directly. However, by properly defining the system variable, the Lie groups dynamic model (1) can be projected into the usual Euclidean space based on the invariance property. LetR k−1 be an estimate for R k−1 at instant k − 1 and R k =R k−1 Ω k−1 be the prediction for R k . Define the invariant error ξ k = Exp −1 G R k R −1 k , then the evolution model of invariant error ξ k obeys where the BCH(·, ·) of the last step is a first order form of the Baker-Campbell-Hausdorff (BCH) formula [34][35][36]. Note that if the norm of the error term ξ k and w k are small enough so that the second order term in (7) could be negligible, then the matrix Lie groups attitude dynamic (1) for rotation matrix R k could be approximately converted into the following linear state space model in Euclidean space for the invariant error ξ k As for observation model (2), the observation sequence Y k can also be reformulated into a function of the invariant error ξ k , i.e., the new innovation sequence z k ∈ R 6 [14] where are the projected observation matrix and noise, respectively; note that, since the observation noises v k , v k are isotropic and mutually independent noises and withR kR Obviously, based on the definition of the invariant error ξ k , the Lie groups dynamics (1) and observation model (2) could be converted as the Euclidean space model (8) and (9) for the invariant error ξ k (as in Figure 3). Note that both the true rotation state R k−1 and the control rotation Ω k−1 do not appear in the converted model (8) and (9); therefore, the evolving motion of ξ k is actually independent of the true trajectory of rotation R k−1 and control rotation Ω k−1 , which is called the invariance property of symmetry model [15,16].

The Invariant Kalman Filter for Attitude Estimation
As the consequence of invariance property, for the attitude estimation system, the probability distribution of rotation state R k based on the observation Y k is equivalent to that of the invariant error ξ k in R 3 , constituting the basis of invariant Kalman filter. Then, if the process noise obeys classical Gaussian distribution, i.e., w k ∼ N c (w k ; 0 3×1 , Σ w ) ∈ R 3 , the filtering steps for attitude estimation system on Lie group SO(3) can be implemented by mimicking that of the classical Kalman filter for invariant error ξ k in Euclidean space.
For the attitude dynamic (1) and associated observation (2), invariant Kalman filter (IKF) aims to recursively predict and correct the estimate for the rotation state. Given the initial estimateR 0|0 = Exp G ξ 0|0 R 0 for true R 0 with errorξ 0|0 ∼ N c ξ 0|0 ; 0 3×1 , Σ 0|0 , as the prediction for prior error estimateξ k|k−1 based on (8), the prior error rotation estimatê R k|k−1 based on the dynamic model (1) can be implemented as followŝ whereR k|k−1 andR k−1|k−1 are the prior and posterior error estimate for true rotation state R k , respectively;ξ k|k−1 andξ k−1|k−1 are associated with the prior and posterior error estimate for true error ξ k . The covariance propagation of rotation group R k is same as that of the classical Kalman filter for invariant error ξ k based on the evolution model (8), where Σ k−1|k−1 stands for the covariance of the posterior estimateξ k|k−1 andR k|k−1 . Therefore, based on the converted innovation (9), the correction steps of IKF for the rotation statê R k|k−1 can be implemented by updating the estimateR k|k−1 with innovation sequence z k where K k is the Kalman gain matrix, and the posterior estimateR k|k is then updated with Lie exponential operation of the correction term K k z k as in (12) [16,17]. A graphical presentation of the filtering diagram of invariant Kalman filter is presented in Figure 4.

The Attitude Estimation Problem with the Trouble of Heavy-Tailed Process Noise
Obviously, the invariant Kalman filter (IKF) is a general extension of the Euclidean space Kalman filter to matrix Lie groups systems based on the invariance property of symmetry dynamics [16]. However, for the attitude estimation in aerospace engineering, there is still some difficulty for practical implementation of the IKF-based attitude estimation because the presence of heavy-tailed process noise in attitude dynamic might introduce great challenge and trouble to the theory of IKF: Firstly, in the filtering step (10)~(14) of IKF, the process noise of attitude was assumed to be a concentrated Gaussian distribution and the statistics parameter is also required to be accurately known (as in Figure 4). However, due to the presence of severely maneuvering operations in space tasks, there might be some unpredictable disturbances or substantial outliers that would contaminate the Gaussian distribution and violate the assumptions for Kalman filtering on matrix Lie groups. Therefore, the theoretical foundation for optimal Kalman filtering is sure to be destroyed, which will degrade the final estimation.
Secondly, for the invariant projection of the Lie groups dynamics (1) into the linear vector form (8) in Section 2.2, the first order approximation to BCH formula is employed based on the assumption that both ξ k and w k were small enough and the second and higher order residual could be dropped from (7) to (8). However, for aerospace engineering applications, the unpredictable disturbances or substantial outliers caused by severely maneuvering operations might violate this condition and lead to inaccurate approximation to the system model, which is prone to a potential source of error and influences the estimation accuracy of IKF.
For attitude determination in astronautical engineering, it is rather meaningful to further study Lie groups estimation problems with heavy-tailed process noise and investigate new robust adaptive version of invariant Kalman filtering methods to improve the filtering precisions for cases with heavy-tailed process noises.

Probability View of Attitude Estimation with Heavy-Tailed Process Noise
For attitude estimation, if the process noise w k in dynamic model (1) is heavy-tailed and contaminated by larger outliers, the Gaussian distribution and concentrated Gaussian distribution cannot accurately describe the probability distribution of heavy-tailed process noise. In this work, to account for the larger outliers due to severely maneuvering, the probability of the heavy-tailed noise is described as the student's t distribution where p(w k ) denotes the probability density of the process noise w k ; S t (·; µ, Σ, τ) denotes the student's t distribution with mean vector µ, scale matrix Σ, and degree of parameter τ; γ w is an auxiliary random variable in (15); and Γ(γ; α, β) represents the Gamma probability density function of the scalar γ with α and β as the shape and rate parameter, i.e., 3. Robust Student's t Based Invariant Kalman Filter for Attitude Estimation on SO(3)

Probability View of Attitude Estimation with Heavy-Tailed Process Noise
For attitude estimation, if the process noise k w in dynamic model (1) is heavy-tailed and contaminated by larger outliers, the Gaussian distribution and concentrated Gaussian distribution cannot accurately describe the probability distribution of heavy-tailed process noise. In this work, to account for the larger outliers due to severely maneuvering, the probability of the heavy-tailed noise is described as the student's t distribution where   k p w denotes the probability density of the process noise k w ; notes the student's t distribution with mean vector  , scale matrix  , and degree of parameter  ; w  is an auxiliary random variable in (15); and   ; , Gamma probability density function of the scalar  with  and  as the shape and rate parameter, i.e., where e is the natural constant and     represents the Gamma function. Obviously, the probability density function of the student's t distribution can be regarded as the infinite mixture of the classical Gaussian probability density function [21,30].
where e is the natural constant and

Probability View of Attitude Estimation with Heavy-Tailed Process Noise
For attitude estimation, if the process noise k w in dynamic model (1) is heavy-tailed and contaminated by larger outliers, the Gaussian distribution and concentrated Gaussian distribution cannot accurately describe the probability distribution of heavy-tailed process noise. In this work, to account for the larger outliers due to severely maneuvering, the probability of the heavy-tailed noise is described as the student's t distribution where   k p w denotes the probability density of the process noise k w ; notes the student's t distribution with mean vector  , scale matrix  , and degree of parameter  ; w  is an auxiliary random variable in (15) where e is the natural constant and     represents the Gamma function. Obviously, the probability density function of the student's t distribution can be regarded as the infinite mixture of the classical Gaussian probability density function [21,30].
(·) represents the Gamma function. Obviously, the probability density function of the student's t distribution can be regarded as the infinite mixture of the classical Gaussian probability density function [21,30].
If given the Gaussian approximation of the posterior probability density function p( ξ k−1 |z 1:k−1 ) and concentrated Gaussian approximation of posterior p( R k−1 |Y 1:k−1 ) for instant k−1, the prior probability density function of p( ξ k |z 1:k−1 ) and p( R k |Y 1:k−1 ) for time instant k can be obtained according to the Chapman-Kolmogorov equation [38,39], where p( ξ k−1 |z 1:k−1 ) denotes the posterior probability density function of the ξ k−1 based on the innovation sequence z k from time instant 1 to k − 1 and p( R k−1 |Y 1:k−1 ) denotes the posterior probability density function of rotation R k based on observation Y k from instant 1 to k − 1. Note that resorting to the equivalence of the R k and ξ k in probability density function, is the student's t distribution p(w k ) for system process noise and so we have p( R k |Y 1:k−1 ) = p( ξ k |z 1:k−1 ). Further, it is evident that the prior probability p( ξ k |z 1:k−1 ) for instant k is not Gaussian but a student's t distribution with mean vector parameterξ k|k−1 , scale matrix Σ k|k−1 and parameter τ k , i.e., whereξ k|k−1 , Σ k|k−1 , and τ k are the mean vector, scale matrix, and degree of freedom parameter of the prior student's t distribution p( R k |Y 1:k−1 ); and γ k is the auxiliary random variable. Employing (15)~ (19), the prior student's t distribution p( R k |Y 1:k−1 ) can be reformulated into a hierarchical composition of the Gaussian distributions Attitude estimation aims to propagate the posterior probability density distribution based on the observation Y k . Since the observation noise in (2) and its equivalent form V k in (9) are zero mean Gaussian white sequence with covariance Σ V , the likelihood probability density function of the innovation sequence z k based on the error ξ k is equal to that of observation Y k given the rotation R k , i.e., p( Y k |R k ) = p( z k |ξ k ) = N c z k ; H ξ ξ k , Σ V where p( Y k |R k ) is the probability density function of Y k given rotation R k , and p( z k |ξ k ) is the probability density function of z k based on error ξ k .
Note that the accurate propagation of the probability density function requires accurate knowledge of scale matrix Σ k|k−1 and degree of freedom τ k . However, for real engineering applications in the presence of unpredictable disturbances or substantial outliers, these two distribution parameters are usually not correct or unavailable in advance. Therefore, from the view of probability density, the main troubles of attitude estimation with heavy-tailed process noise are threefold: (1) With the student's t distribution process noise, the probability density function p( ξ k |z 1:k−1 ) in classical invariant Kalman filter depends on the auxiliary random variable γ k and becomes the hierarchical form p( ξ k |γ k , z 1:k−1 ) in (20) and (21). (2) For the unpredictable disturbances or outliers induced by severely maneuvering operations, the accurate scale matrix Σ k|k−1 and degree of freedom γ k are usually unavailable but essential to propagate the posterior estimates.
As a consequence, the heavy-tailed attitude estimation problem cannot be solved with the invariant Kalman filter in Section 2 since the distribution parameter Σ k|k−1 and τ k should be estimated together with the invariant error ξ k and the rotation state R k .

Prior Probability Definition for the Parameters of Student's t Distribution
To solve Lie groups filtering problems with Student's t distribution process noise, this work aims to estimate the rotation group state R k together with Σ k|k−1 and τ k by variational Bayesian iterations as in Figure 5. The prior distributions of parameter Σ k|k−1 and τ k based on all historical observation sequences should be conjugate prior distributions so that their posterior distribution would maintain the same form as their prior distribution. According to (20) and (21), the estimated scale matrix Σ k|k−1 is proportional to the covariance matrix of Gaussian distribution. Therefore, the commonly used inverse Wishart distribution and the Gamma distribution of Bayesian statistics are employed as the conjugate prior distributions for parameters Σ k|k−1 and γ k , respectively.
The distribution of scale matrix Σ k|k−1 is defined as inverse Winshart distribution [29,39,40] which is often used as conjugate prior for the covariance matrix of the Gaussian variable. For attitude estimation problem, the inverse Winshart function of Σ k|k−1 can be written as where IW denotes the inverse Winshart distribution; λ k denotes the degrees of freedom; Ψ k is the 3 × 3 inverse scale matrix; |·| is the determinant operator; and rameter  ; w  is an auxiliary random variable in (15); and    Gamma probability density function of the scalar  with  and rate parameter, i.e., where e is the natural constant and     represents the Gamma the probability density function of the student's t distribution can be nite mixture of the classical Gaussian probability density function [2 3 (·) is the 3-variate gamma function [29]. For the scale matrix Σ k|k−1 ∼ IW Σ k|k−1 ; λ k , Ψ k , the expectation of the The distribution parameter τ k can be defined as the Gamma distribution [21], where a k , b k are respectively the shape and rate parameter of the Gamma distribution τ k , respectively.
In summary, Equations (20)~(25) actually constitute a new student's t-based hierarchical Gaussian state-space model for the heavy-tailed attitude estimation problem, whose diagram is presented in Figure 6. The initialization of the required distribution parameter is based on the nominal parameter Σ k|k−1 , which will be discussed later on. Using the above distribution definition, the rotation state R k and ξ k , the auxiliary random variable γ k , the scale matrix Σ k|k−1 , and the parameter τ k will be jointly estimated based on the newly constructed hierarchical Gaussian state-space model by variational Bayesian iterations, which then leads to the new robust student's t-based invariant Kalman filter for the heavy-tailed attitude estimation problem.

Variational Beayesian Approximations of Posterior Probability Density Function
To jointly estimate the rotation state R k , auxiliary random variable γ k , scale matrix Σ k|k−1 , and parameter τ k requires calculating the posterior probability density function p ξ k , γ k , Σ k|k−1 , τ k |z 1:k . Since the accurate analytical expression of p ξ k , γ k , Σ k|k−1 , τ k |z 1:k is not available, this work tries to obtain its factored approximate based on variational Bayesian iterations [38]. i.e., where q(ξ k ), q(γ k ), q Σ k|k−1 , and q(τ k ) represent, respectively, the approximate posterior probability density function of ξ k , γ k , Σ k|k−1 , τ k ; q(ξ k ), q(γ k ), q Σ k|k−1 and q(τ k ), which can be determined by minimizing the Kullback-Leibler divergence (KLD) [38,41] between the factored approximate q(ξ k )q(γ k )q Σ k|k−1 q(τ k ) and the true p ξ k , γ k , Σ k|k−1 , τ k |z 1:k [29,38], i.e., p(x) dx denotes the KLD between q(x) and p(x). According to the variational Bayesian iteration, these approximate posterior q(ξ k ), q(γ k ), q Σ k|k−1 , and q(τ k ) to (27) satisfy the following equations [21,38] where log(·) denotes the natural logarithm function; E −ξ k [·] represents the mathematical expectation operation for all variables except ξ k , while , and E −τ k [·] are with similar definitions; c ξ k is the constant with respect to ξ k , while c γ k , c Σ k|k−1 , and c τ k are similarly defined. Obviously, Equations (28)~(31) are coupled and so no analytical solutions to the above equations exist. Therefore, only the approximate solutions are available and in this work the fixed-point iterations are employed to determine these parameters [38].
According to the conditional independence properties, the joint probability density function p ξ k , γ k , Σ k|k−1 , τ k , z 1:k can be factored as follows Then, using (20)~(26), we have Using (16) and (24), log p ξ k , γ k , Σ k|k−1 , τ k , z 1:k can be formulated as tion for optimal Kalman filtering is sure to be destroyed, which will degrad estimation. Secondly, for the invariant projection of the Lie groups dynamics (1) int vector form (8) in Section 2.2, the first order approximation to BCH formula is based on the assumption that both k  and k w were small enough and the s higher order residual could be dropped from (7) to (8). However, for aerospac ing applications, the unpredictable disturbances or substantial outliers caused maneuvering operations might violate this condition and lead to inaccurate a tion to the system model, which is prone to a potential source of error and inf estimation accuracy of IKF.
For attitude determination in astronautical engineering, it is rather me further study Lie groups estimation problems with heavy-tailed process noise tigate new robust adaptive version of invariant Kalman filtering methods to im filtering precisions for cases with heavy-tailed process noises.

Probability View of Attitude Estimation with Heavy-Tailed Process Noise
For attitude estimation, if the process noise k w in dynamic model (1) is h and contaminated by larger outliers, the Gaussian distribution and concentrat distribution cannot accurately describe the probability distribution of heavy-tai noise. In this work, to account for the larger outliers due to severely maneu probability of the heavy-tailed noise is described as the student's t distribution where   k p w denotes the probability density of the process noise k w ; notes the student's t distribution with mean vector  , scale matrix  , and de rameter  ; w  is an auxiliary random variable in (15); and   ; , Gamma probability density function of the scalar  with  and  as the rate parameter, i.e., where e is the natural constant and     represents the Gamma function.
the probability density function of the student's t distribution can be regarded nite mixture of the classical Gaussian probability density function [21,30].

Fixed-Point Iteration of the System State and Distribution Parameters
Based on the conditional independence property of the invariant error ξ k and parameters γ k , Σ k|k−1 , and τ k , the Equations (28)~(31) can be rewritten with above Equation (34) and fixed-point methods can be used to obtain the approximate estimates to the invariant error ξ k as well as the parameters γ k ,Σ k|k−1 , and τ k [21].

Fixed-Point Iteration of the Invariant Error ξ k
With (34) and the independence property of γ k and Σ k|k−1 , (28) can be written as where q (i+1) (·) is the iterative approximation to q(·) at i + 1th iteration; E (i) −ξ k [·] represents the mathematical expectation at ith iteration for variables except the invariant error ξ k Machines 2021, 9, 182 11 of 19 while E (i) γ k [·] and E (i) [·] denote the mathematical expectation at ith iteration for γ k and Σ k|k−1 , respectively; Σ (i) k|k−1 is the modified prior error covariance matrix calculated as Then, q (i+1) (ξ k ) can also be updated as a Gaussian distribution, i.e., whereξ (i+1) k|k and Σ (i+1) k|k can be determined by following the update stepŝ where K (i+1) k denotes the calibrated Kalman gain matrix at the i + 1th iteration.

Fixed-Point Iteration of the Auxiliary Random Variable γ k
Using (34) and conditional independency, the log q (i+1) (γ k ) in (29) can be rewritten as where E (i) −γ k [·] represents the mathematical expectation at ith iteration for variables except the γ k while E (i) τ k [·] denotes the mathematical expectation at ith iteration for the parameter τ k . Note that the approximateξ (i+1) k|k and Σ (i+1) k|k updated in (35)~ (40) can be used to propagate Equation (41) and the term E (i) Then (41) could be deduced into which is actually also a Gamma distribution as (16), i.e., where η (i+1) and θ (i+1) are the shape and rate parameter determined by With (34), (37), (44) and the independence of ξ k and γ k , (30) can be rewritten as Therefore, using (24) and (47), q (i+1) Σ k|k−1 is also an inverse Wishart distribution as

Fixed-Point Iteration of the Prior Estimate for Parameter τ k
With (34), (37), (44), and (49), (33) can be rewritten as Firstly, in the filtering step (10)~(14) of IKF, the process noise of attitude was assumed to be a concentrated Gaussian distribution and the statistics parameter is also required to be accurately known (as in Figure 4). However, due to the presence of severely maneuvering operations in space tasks, there might be some unpredictable disturbances or substantial outliers that would contaminate the Gaussian distribution and violate the assumptions for Kalman filtering on matrix Lie groups. Therefore, the theoretical foundation for optimal Kalman filtering is sure to be destroyed, which will degrade the final estimation.
Secondly, for the invariant projection of the Lie groups dynamics (1) into the linear vector form (8) in Section 2.2, the first order approximation to BCH formula is employed based on the assumption that both k  and k w were small enough and the second and higher order residual could be dropped from (7) to (8). However, for aerospace engineering applications, the unpredictable disturbances or substantial outliers caused by severely maneuvering operations might violate this condition and lead to inaccurate approximation to the system model, which is prone to a potential source of error and influences the estimation accuracy of IKF.
For attitude determination in astronautical engineering, it is rather meaningful to further study Lie groups estimation problems with heavy-tailed process noise and investigate new robust adaptive version of invariant Kalman filtering methods to improve the filtering precisions for cases with heavy-tailed process noises.

The Variational Beayesian Iteration Based Robust Student's t Invariant Kalman Filter
To propagate the approximate posterior distributions, using (44), (49) and (54), the in Lie algebra of 5000 random runs is used to calculate the root mean square error RMSE k during the 5000 runs' filtering time and average root mean square error ARMSE, i.e., where ξ k,l denotes the true estimate at the k-th time instant of the l-th simulation run andξ k,l is the corresponding estimate for the true ξ k,l ; · denotes the Euclidean vector norm.
In this work, three cases of heavy tailed process noises are simulated with α = 5, 20 and β = 100, 1000 to investigate the filtering performance of different methods. The ARMSE result data of different filtering methods are given in Figure 7 while the RMSE k result data during the first 500 runs' filtering processes are presented in Figures 8 and 9.   The ARMSE result of IKF in Figure 7 and RMSE k data of IKF in Figures 8-11 are always the largest in all methods for all cases, which certifies that the estimation performance of IKF is severely corrupted by heavy tailed outliers.  In STF, the noise distribution is assumed as a given Student's t distribution, and the ARMSE result in Figure 7 and RMSE k in Figures 8-11 is smaller than that of IKF. However, since it still requires accurate prior knowledge of the noise distribution parameter, the unknown or inaccurate parameters still mislead the estimate results.
The QeIKF aims to estimate, online, the unknown noise covariance parameter and also output better results than IKF and STF; however, it assumes the distribution as an Gaussian one so its ARMSE result in Figure 7 and RMSE k in Figure 8 for the case α = 5 become worse than the case of α = 20 in Figure 9; while, as shown in Figures 10 and 11, its filtering precision would be significantly degraded by large outliers with β = 1000.
HKF shows some robustness to heavy tailed outliers IKF, STK, and QeKF for the cases of β = 1000 as in Figures 10 and 11. Note that, the estimate given by HKF is rather conservative and its robustness is at the cost of inferior precision.
In all simulation cases of α and β, the ARMSE and RMSE k result of proposed RSIKF is always the least one in Figures 7-11, demonstrating its superiority filtering performance to the other available methods. For RSIKF, different setting of N from 1~15 are conducted with all cases of α with the ARMSE result given in Figure 12; obviously, N = 10 already guarantees the iterative convergence and a larger N would significantly increase the computational cost.

Conclusions
For aerospace, satellite, and robotics engineering, the matrix Lie groups attitude estimation problem with heavy-tailed process noise is investigated. An improved robust Student's t invariant Kalman filtering is proposed based on a hierarchical Gaussian statespace model. The probability density function of state prediction is defined based on student's t distribution, while the conjugate prior distributions of the scale matrix and degrees of freedom (dofs) parameter are respectively formulated as the inverse Wishart and Gamma distribution. For attitude estimation state-space model, the Lie groups rotation matrix of spacecraft attitude is inferred together with the scale matrix and dof parameter using the variational Bayesian iteration. The new approach can improve the filtering robustness of the invariant Kalman filter for Lie groups spacecraft attitude estimation problems with heavy-tailed process uncertainty. Some further applications of the proposed methods can also be extended to other problems such as state smoother, parameter identification, state observers, and so on.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available as the data also forms part of an ongoing study.