Laplace ℓ 1 robust Student’s T-filter for attitude estimation of satellites

Abstract

A Laplace ℓ₁ robust student’s T-filter is presented for satellites to estimate their attitude state despite severe measurement noise and modeling error. Although the student’s T-filter (STF) has the capability of handling measurement noise, it cannot address the unknown modeling error. It is further sensitive to the degree of freedom (DOF). Hence, the measurement covariance is updated by using the maximum correntropy criterion to accommodate the covariance of unknown modeling error in robust filtering design. Moreover, the Laplace distribution is derived to reduce the influence of the DOF parameter by forming a surrogate function of an optimization problem. Then, the majorization minimization approach is formulated to solve such an optimization problem and present the proposed filtering algorithm in the STF framework. Numerical simulation of applying the proposed attitude estimation scheme is performed by comparing it with the third/fifth-order cubature Kalman filters.

Introduction

Nonlinear filters play a vital role in signal processing, ¹ target tracking,^2–4 and navigation.^5–7 In the Bayesian estimation framework, filters are designed by solving the posterior probability density function (PDF). However, there does not exist an optimal solution to the posterior PDF for nonlinear systems, since the closed PDF does not exist. ⁸ The posterior PDF is thus assumed to satisfy the Gaussian distribution. Accordingly, many approximation methods are applied to design the suboptimal nonlinear filters such as the classical Kalman filter (KF) and its extensions.^9–14 For those filters, both the system states and the measurement predictions are assumed to have Gaussian characteristics. Considering space missions such as the orbit tracking or the attitude determination, the noise and the unknown perturbations may be induced. The measurement outliers may also exist. ¹⁵ The noise of the system and its measurement may have a heavier tail than the distribution with the Gaussian characteristics. ¹⁶ Moreover, unknown perturbations bring satellites with modeling error. This error will further increase the nonlinearity and the uncertainty of the noise. This complicates the state estimation further and deteriorates the estimation performance of KF.

The past decade has witnessed significant development of the filter theory in dealing with the heavy-tailed noise. Many investigations are available and they can be categorized into two types. The first is the regression-based robust filters. ¹⁷ The Huber theory^18,19 is applied to design a robust KF. The maximum correntropy criterion^20,21 is applied to maximize the prediction error and the residual correntropy with the non-Gaussian noise handled. However, their performance is sensitive to the estimation gains. This drawback drives the researchers to design alternative solutions to the heavy-tailed noise problem.

Some achievements focused on the PDF distributions describing the heavy-tailed noise. For example, the particle filter (PF) adopts a stochastic approximation method to address the non-Gaussian noise. For this filter, random sampling points are adopted to approximate the PDF via the sequence of the Monte Carlo sampling technique.^22,23 However, the PF has a computational burden because the considered particles’ number will increase exponentially along with the state’s dimension. The mixture Gaussian KF was derived by assuming that the sum of the Gaussians can approximate any density with arbitrary accuracy.^24,25 However, it is not easy to determine the number of filters and the corresponding computation burden is large.

An efficient approach to solve the above issue is the application of the Student’s T distribution to approximate the posterior PDF. This scheme is named as the Student’s T-filter (STF), in which the predicted and the posterior PDF are supposed to satisfy the Student’s T distribution. ²⁶ Therefore, the STF has significant advantages when dealing with the heavy-tailed noise. Many investigations focusing on this method have been presented in recent years.^27–29 To deal with the unknown modeling error, the maximum correntropy criterion is considered here to update the STF’s measurement covariance. This can compensate for the covariance of the unknown modeling error to suppress it for the robust estimation design. Besides, the estimation performance is also dependent on the degree-of-freedom (DOF) of the Student’s T distribution. The selection of DOF for STF is not easy. However, the Laplace distribution can solve this drawback. It can be implemented by optimizing a hybrid ℓ₂/ℓ₁ norm issue. Then, the surrogate function in a quadratic format is derived, ³⁰ which is employed to get the approximation of cost function. Then, the majorization minimization (MM) technique is used to obtain the closed-form expression of state estimation by an iterative algorithm in the STF framework. This approach is called as the Laplace ℓ₁ Robust STF (L1-RSTF).

The rest of the work is organized as follows. In section “Problem formulation,” the theory of the Student’s T-filter is introduced. In section “Development of the Laplace ℓ₁ robust Student’s T-filter,” the proposed L1-RSTF is presented based on the STF by using the maximum correntropy criterion and the Laplace distribution. Simulation results of the attitude estimation for small satellites are shown in section “Numerical simulation.” Some concluding remarks are given in section “Conclusion.”

Problem formulation

In this paper, the mathematical notations are standard defined as follows. The set of $m\times n$ matrices is represented by ${\Re }^{m\times n}$ . The set of all non-negative integers is denoted as $N_0$ . The Euclidean norm of vectors or matrices is denoted by $||·||$ . $(·{)}^{\mathrm{T}}$ denotes the transpose of matrix or vector. For any vector $\mathit{a}=[\begin{array}{cccc}{a}_1\hfill & a_2\hfill & \cdots \hfill & a_n\hfill \end{array}{]}^{\mathrm{T}}\in {\Re }^n$ , $||\mathit{a}|{|}_1=\sum _{i=1}^n|a_i|$ and $\mathit{diag}(\mathit{a})=\left[\begin{array}{cccc}{a}_1\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& a_2\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & a_n\hfill \end{array}\right]$ are defined; moreover, $(\mathit{a}{)}_i$ denotes the ith element of $\mathit{a}$ , $i=1,2,\cdots ,n$ . $\mathit{I}$ is the identity matrix with appropriate dimension. $\mathit{0}$ is a zero vector or matrix having an appropriate dimension. $\mathit{E}(·)$ is the mathematical expectation of a random signal. $\det (·)$ denotes the determinant operation. $\mathit{\Gamma }(·)$ is the Gamma function. For any vector $\mathit{b}=[\begin{array}{ccc}{b}_1\hfill & b_2\hfill & b_3\hfill \end{array}{]}^{\mathrm{T}}\in {\Re }^3$ , ${\mathit{b}}^{\times }=\left[\begin{array}{ccc}0\hfill & -b_3\hfill & b_2\hfill \\ \multicolumn{1}{c}{b_3}& 0\hfill & -b_1\hfill \\ \multicolumn{1}{c}{-b_2}& b_1\hfill & 0\hfill \end{array}\right]$ is defined. “w.p.” represents “with probability.” $\mathbf{\chi }({\mathit{\mu }}_0,{\mathbf{\Sigma }}_{\mathbf{\chi }})$ is the normal distribution with mean ${\mathit{\mu }}_0$ and covariance ${\mathbf{\Sigma }}_{\mathbf{\chi }}$ .

Mathematical model of a class of nonlinear systems

The following discrete-time nonlinear system is considered in this study:

$$\{\begin{array}{c}{\mathit{x}}_{k+1}=\mathit{f}({\mathit{x}}_k)+{\mathit{w}}_k\hfill \\ \multicolumn{1}{c}{{\mathit{y}}_k=\mathit{h}({\mathit{x}}_k)+{\mathit{v}}_k}\end{array}$$

(1)

where $k\in N_0$ is the discrete-time index. ${\mathit{x}}_k\in {\Re }^n$ is the system state at the discrete-time $t_k$ , ${\mathit{y}}_k\in {\Re }^m$ is the system’s measurement at $t_k$ , ${\mathit{w}}_k\in {\Re }^n$ is the process noise at $t_k$ , and ${\mathit{v}}_k\in {\Re }^m$ is measurement noise at $t_k$ . $\mathit{f}({\mathit{x}}_k)$ and $\mathit{h}({\mathit{x}}_k)$ are known, and they are the process and the measurement functions, respectively.

Assuming that those noises ${\mathit{w}}_k$ and ${\mathit{v}}_k$ satisfy the heavy-tailed distribution, then the PDFs of those noises are supposed to be the following Student’s-t distribution:

$$p({\mathit{w}}_k)=\mathbf{St}({\mathit{w}}_k;0,{\mathit{Q}}_k,{\upsilon }_1)$$

(2)

$$p({\mathit{v}}_k)=\mathbf{St}({\mathit{v}}_k;0,{\mathit{R}}_k,{\upsilon }_2)$$

(3)

where $\mathbf{St}(·;\mathit{\mu },\mathbf{\Sigma },\upsilon )$ is the Student’s T PDF with mean $\mathit{\mu }$ , scale matrix $\mathbf{\Sigma }$ , and DOF $\upsilon$ . Moreover, the specific expression of the Student’s T PDF with any variable $\mathit{z}$ could be given by

$$\begin{array}{c}p(z)=\mathrm{St}(z;\mathit{\mu },\Sigma ,\upsilon )=\frac{r(\frac{\upsilon +2}{2})}{r(\frac{\upsilon }{2})}\frac{1}{\sqrt{\left|\upsilon \pi \Sigma \right|}}\hfill \\ {[1+\frac{1}{\upsilon }{(z-\mathit{\mu })}^T{\Sigma }^{-1}(z-\mathit{\mu })]}^{-(\frac{\upsilon +n_z}{2})}\hfill \end{array}$$

(4)

For $\upsilon >2$ , the mean and the variance of $\mathit{z}$ can be calculated by $\mathit{E}(\mathit{z})=\mathit{\mu }$ and $\mathrm{var}(\mathit{z})=\mathit{E}((\mathit{z}-\mathit{\mu })(\mathit{z}-\mathit{\mu }{)}^{\mathrm{T}})=\frac{\upsilon }{\upsilon -2}\mathbf{\Sigma }$ , respectively.

It is known from the above definition that ${\mathit{Q}}_k$ and ${\upsilon }_1$ are the scale matrix and the process noise’s DOF, respectively. ${\mathit{R}}_k$ and ${\upsilon }_2$ are the scale matrix and the measurement noise’s DOF, respectively. Then, the initial value of ${\mathit{x}}_0$ can be assumed as the distribution of the Student’s T characteristics with mean ${\hat{\mathit{x}}}_{0|0}$ , scale matrix ${\mathit{P}}_{0|0}$ , and DOF ${\upsilon }_3$ as

$$p({\mathit{x}}_0)=\mathbf{St}({\mathit{x}}_0;{\hat{\mathit{x}}}_{0|0},{\mathit{P}}_{0|0},{\upsilon }_3)$$

(5)

where ${\mathit{x}}_0$ , ${\mathit{w}}_k$ , and ${\mathit{v}}_k$ are mutually uncorrelated.

Student’s T filtering theory

When applying the Bayesian theory framework to develop the Student’s T-filter, the posterior PDF $p({\mathit{x}}_k|{\mathit{Y}}_k)$ is calculated via the Bayesian relation

$$\begin{array}{c}\multicolumn{1}{c}{p({\mathit{x}}_k|{\mathit{Y}}_k)=\frac{p({\mathit{x}}_k|{\mathit{Y}}_{k-1})p({\mathit{y}}_k|{\mathit{x}}_k)}{{\int }_{{\Re }^n}p({\mathit{x}}_k|{\mathit{Y}}_{k-1})p({\mathit{y}}_k|{\mathit{x}}_k)d{\mathit{x}}_k}}\\ =\frac{p({\mathit{x}}_k,{\mathit{y}}_k|{\mathit{Y}}_{k-1})}{{\int }_{{\Re }^n}p({\mathit{x}}_k,{\mathit{y}}_k|{\mathit{Y}}_{k-1})d{\mathit{x}}_k}\hfill \end{array}$$

(6)

which represents the measurement update for the estimation.

Based on (6), it is ready to get the predicted PDF $p({\mathit{x}}_k|{\mathit{Y}}_{k-1})$ by the Chapman-Kolmogorov equation for the time update of the estimation, which is specified as

$$p({\mathit{x}}_k|{\mathit{Y}}_{k-1})={\int }_{{\Re }^n}p({\mathit{x}}_k|{\mathit{x}}_{k-1})({\mathit{x}}_{k-1}|{\mathit{Y}}_{k-1})d{\mathit{x}}_{k-1}$$

(7)

The PDFs $p({\mathit{x}}_k|{\mathit{x}}_{k-1})$ and $p({\mathit{y}}_k|{\mathit{x}}_k)$ are got from the process function and the measurement function in (1). The PDFs of the correlated noises are updated as

$$p({\mathit{x}}_k|{\mathit{x}}_{k-1})=\mathbf{St}({\mathit{x}}_k;\mathit{f}({\mathit{x}}_{k-1}),{\mathit{Q}}_k,{\upsilon }_1)$$

(8)

$$p({\mathit{y}}_k|{\mathit{x}}_k)=\mathbf{St}({\mathit{y}}_k;\mathit{h}({\mathit{x}}_k),{\mathit{R}}_k,{\upsilon }_2)$$

(9)

To obtain a solution to the posterior PDF (7), the jointly predicted PDF $p({\mathit{x}}_k,{\mathit{y}}_k|{\mathit{Y}}_{k-1})$ based on the state ${\mathit{x}}_k$ and the measurement ${\mathit{y}}_k$ is assumed as the Student’s-t distributed ²⁶ :

$$\begin{array}{c}p({\mathit{x}}_k,{\mathit{x}}_k|{\mathit{Y}}_{k-1})=\hfill \\ \mathbf{St}(\left[\begin{array}{c}{\mathit{x}}_k\hfill \\ \multicolumn{1}{c}{{\mathit{y}}_k}\end{array}\right];\left[\begin{array}{c}{\hat{\mathit{x}}}_{k|k-1}\hfill \\ \multicolumn{1}{c}{{\hat{\mathit{y}}}_{k|k-1}}\end{array}\right],\left[\begin{array}{cc}{\mathit{P}}_{k|k-1}\hfill & {\mathit{P}}_{k|k-1}^{\mathrm{xy}}\hfill \\ \multicolumn{1}{c}{{\mathit{P}}_{k|k-1}^{\mathrm{yx}}}& {\mathit{P}}_{k|k-1}^{\mathrm{yy}}\hfill \end{array}\right],{\upsilon }_3)\hfill \end{array}$$

(10)

where ${\hat{\mathit{x}}}_{k|k-1}$ is the mean vector, and ${\mathit{P}}_{k|k-1}$ is the scale matrix of the predicted PDF $p({\mathit{x}}_k|{\mathit{Y}}_{k-1})$ . ${\hat{\mathit{y}}}_{k|k-1}$ and ${\mathit{P}}_{k|k-1}^{\mathrm{yy}}$ are the mean and the scale matrix of the PDF $p({\mathit{y}}_k|{\mathit{Y}}_{k-1})$ , respectively. ${\mathit{P}}_{k|k-1}^{\mathrm{xy}}$ is the cross-scale matrices of the system state, and ${\mathit{P}}_{k|k-1}^{\mathrm{yx}}$ is the cross-scale matrices of the measurement.

To this end, the posterior PDF $p({\mathit{x}}_k|{\mathit{Y}}_k)$ in (7) can be calculated by the Student’s T distribution. In accordance, the generalized algorithm of STF can be summarized as follows.

• (1) Time update:

$$\begin{array}{c}{\hat{\mathit{x}}}_{k|k-1}=\mathit{E}({\mathit{x}}_k|{\mathit{Y}}_{k-1})=\hfill \\ {\int }_{{\Re }^n}\mathit{f}({\mathit{x}}_{k-1})\mathbf{St}({\mathit{x}}_{k-1};{\hat{\mathit{x}}}_{k-1|k-1},{\mathit{P}}_{k-1|k-1},{\upsilon }_3)d{\mathit{x}}_{k-1}\hfill \end{array}$$

(11)

$$\begin{array}{c}{\mathit{P}}_{k|k-1}=\frac{{\upsilon }_3-2}{{\upsilon }_3}\mathit{E}(({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}{)}^{\mathrm{T}}|{\mathit{Y}}_{k-1})\hfill \\ =\frac{{\upsilon }_3-2}{{\upsilon }_3}{\int }_{{\Re }^n}\mathit{f}({\mathit{x}}_{k-1}){(\mathit{f}({\mathit{x}}_{k-1}))}^{\mathrm{T}}\mathbf{St}({\mathit{x}}_{k-1};{\hat{\mathit{x}}}_{k-1|k-1},\hfill \\ {\mathit{P}}_{k-1|k-1},{\upsilon }_3)d{\mathit{x}}_{k-1}-\frac{{\upsilon }_3-2}{{\upsilon }_3}{\hat{\mathit{x}}}_{k|k-1}{\hat{\mathit{x}}}_{k|k-1}^{\mathrm{T}}+\frac{{\upsilon }_1({\upsilon }_3-2)}{({\upsilon }_1-2){\upsilon }_3}{\mathit{Q}}_{k-1}\hfill \end{array}$$

(12)

• (2) Measurement update:

$$\begin{array}{c}{\hat{\mathit{y}}}_{k|k-1}=\mathit{E}({\mathit{y}}_k|{\mathit{Y}}_{k-1})=\hfill \\ {\int }_{{\Re }^n}\mathit{h}({\mathit{x}}_k)\mathbf{St}({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1},{\upsilon }_3)d{\mathit{x}}_k\hfill \end{array}$$

(13)

$$\begin{array}{c}{\mathit{P}}_{k|k-1}^{\mathrm{yy}}=\frac{{\upsilon }_3-2}{{\upsilon }_3}\mathit{E}(({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1})({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1}{)}^{\mathrm{T}}|{\mathit{Y}}_{k-1})\hfill \\ =\frac{{\upsilon }_3-2}{{\upsilon }_3}{\int }_{{\Re }^n}\mathit{h}({\mathit{x}}_k){(\mathit{h}({\mathit{x}}_k))}^{\mathrm{T}}\mathbf{St}({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1},{\upsilon }_3)d{\mathit{x}}_k\hfill \\ -\frac{{\upsilon }_3-2}{{\upsilon }_3}{\hat{\mathit{y}}}_{k|k-1}{\hat{\mathit{y}}}_{k|k-1}^{\mathrm{T}}+\frac{{\upsilon }_1({\upsilon }_3-2)}{({\upsilon }_1-2){\upsilon }_3}{\mathit{R}}_{k-1}\hfill \end{array}$$

(14)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_{k|k-1}^{\mathit{x}y}=\frac{{\upsilon }_3-2}{{\upsilon }_3}\mathit{E}(({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1}{)}^{\mathrm{T}}|{\mathit{Y}}_{k-1})}\\ =\frac{{\upsilon }_3-2}{{\upsilon }_3}{\int }_{{\Re }^n}\mathit{f}({\mathit{x}}_{k-1}){(\mathit{h}({\mathit{x}}_k))}^{\mathrm{T}}\mathbf{St}({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1},{\upsilon }_3)d{\mathit{x}}_k\hfill \\ -\frac{{\upsilon }_3-2}{{\upsilon }_3}{\hat{\mathit{x}}}_{k|k-1}{\hat{\mathit{y}}}_{k|k-1}^{\mathrm{T}}\hfill \end{array}$$

(15)

where ${\mathit{P}}_{k|k-1}^{\mathrm{xy}}=({\mathit{P}}_{k|k-1}^{\mathrm{yx}}{)}^{\mathrm{T}}$ .

• (3) State update:

$${\mathrm{\Delta }}_k=\sqrt{{({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1})}^{\mathrm{T}}{\mathit{P}}_{k|k-1}^{\mathrm{yy}}({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1})}$$

(16)

$${\mathit{K}}_k={\mathit{P}}_{k|k-1}^{\mathrm{xy}}({\mathit{P}}_{k|k-1}^{\mathrm{yy}}{)}^{-1}$$

(17)

$${\hat{\mathit{x}}}_{k|k}={\hat{\mathit{x}}}_{k|k-1}+{\mathit{K}}_k({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1})$$

(18)

$${\mathit{P}}_{k|k}=\frac{({\upsilon }_3-2)({\upsilon }_3+{\mathrm{\Delta }}_k^2)}{{\upsilon }_3({\upsilon }_3+m-2)}({\mathit{P}}_{k|k-1}-{\mathit{K}}_k{\mathit{P}}_{k|k-1}^{\mathrm{yy}}{\mathit{K}}_k^{\mathrm{T}})$$

(19)

The generalized STF could be obtained by calculating (11)–(15) with the Student’s T weighted integral, and then by updating the estimated state and the scale matrix (16)–(19). However, the DOF parameters of the STF are difficult to select and they will have an important effect on the estimation accuracy. This issue will be addressed in the next section.

Development of the Laplace ℓ₁ robust Student’s T-filter

Maximum correntropy criterion design

In this section, the maximum correntropy criterion is adopted to update the measurement covariance of the unknown modeling error for the robust estimation scheme design. Therefore, the maximum correntropy criterion is introduced first. For any random variables ${\mathit{G}}_1$ and ${\mathit{G}}_2$ , their correntropy is denoted as ³¹

$$\mathit{V}({\mathit{G}}_1,{\mathit{G}}_2)=\mathit{E}(\mathit{\kappa }({\mathit{G}}_1,{\mathit{G}}_2))=\int \int \mathit{\kappa }({\mathit{g}}_1,{\mathit{g}}_2)d{\mathit{f}}_{{\mathit{G}}_1,{\mathit{G}}_2}({\mathit{g}}_1,{\mathit{g}}_2)$$

(20)

where $\mathit{\kappa }(·)$ denotes the kernel function and it is positive-definite. The function ${\mathit{f}}_{{\mathit{G}}_1,{\mathit{G}}_2}({\mathit{g}}_1,{\mathit{g}}_2)$ denotes the density of the joint probability.

In practical engineering, the joint probability cannot be obtained directly. It is usually represented by a series of data samples $\{(g_{1i},g_{2i}),i=1,2,\cdots N\}$ with finite numbers. Then, one can write the correntropy as

$$\mathit{V}({\mathit{G}}_1,{\mathit{G}}_2)\approx \frac{1}{N}\sum _{i=1}^N\mathit{\kappa }({\mathit{g}}_1,{\mathit{g}}_2)$$

(21)

According to Wang et al., ³² the following Gaussian kernel function can be employed for (19).

$${\mathit{\kappa }}_{\mathit{G}}({\mathit{e}}_i)={\mathit{\kappa }}_{\mathit{G}}({\mathit{g}}_{1i}-{\mathit{g}}_{2i})=\frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{{\mathit{e}}_i^2}{2{\sigma }^2})$$

(22)

where $\sigma$ is the correntropy’s kernel width. It is found that the kernel ${\mathit{\kappa }}_{\mathit{G}}({\mathit{e}}_i)$ will have its maximum value when ${\mathit{e}}_i=0$ , which can be used as the cost function of the maximum correntropy criterion.

$${\mathit{J}}_{\mathrm{MCC}}=min\sum _{i=1}^N({\mathit{\kappa }}_{\mathit{G}}(0)-{\mathit{\kappa }}_{\mathit{G}}({\mathit{e}}_i))$$

(23)

Considering that the investigation focuses on updating the measurement covariance to suppress the unknown modeling error, the new cost function is proposed by combining (23) and the STF framework with the definition of ${\mathit{e}}_k=(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\mathit{R}}_k{)}^{-\frac{1}{2}}({\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k)$ as

$$\begin{array}{c}{\hat{\mathit{x}}}_{k|k}=\arg min\hfill \\ (\frac{1}{2}||{\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}|{|}_{{(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1})}^{-1}}^2+\frac{1}{2}\sum _{i=1}^m({\mathit{\kappa }}_{\mathit{G}}(0)-{\mathit{\kappa }}_{\mathit{G}}({\mathit{e}}_{k,i})))\hfill \end{array}$$

(24)

with $||{\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}|{|}_{{(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1})}^{-1}}^2=({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}{)}^{\mathrm{T}}(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1}{)}^{-1}({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})$ . ${\hat{\mathit{x}}}_{k|k}$ is the system state estimations after the measurement update. ${\hat{\mathit{x}}}_{k|k-1}$ is system state estimations before measurement update. ${\mathit{P}}_{k|k-1}$ is the likelihood of the scalar matrix of the PDF $p({\mathit{x}}_k|{\mathit{Y}}_{k-1})$ . ${\mathit{H}}_k$ is the Jacobian matrix of the nonlinear function $\mathit{h}(·)$ for the following theoretical derivation.

Differentiating (24) with respect to ${\mathit{x}}_k$ yields

$$(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1}{)}^{-1}({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})-\sum _{i=1}^m\frac{\partial {\mathit{\kappa }}_{\mathit{G}}({\mathit{e}}_{k,i})}{\partial {\mathit{e}}_{k,i}}\frac{\partial {\mathit{e}}_{k,i}}{\partial {\mathit{x}}_k}=\mathit{0}$$

(25)

Inserting (22) into (23), one has

$$(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1}{)}^{-1}({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})-\sum _{i=1}^m{\psi }_{k,i}{\mathit{e}}_{k,i}\frac{\partial {\mathit{e}}_{k,i}}{\partial {\mathit{x}}_k}=\mathit{0}$$

(26)

where ${\psi }_{k,i}=\frac{1}{\sqrt{2\pi }{\sigma }^3}\exp (-\frac{{\mathit{e}}_{k,i}^2}{2{\sigma }^2})$ , $i=1,2,\cdots ,m$ .

Define the matrix ${\mathit{\psi }}_k$ as ${\mathit{\psi }}_k=\mathit{diag}([\begin{array}{cccc}{\psi }_{k,1}\hfill & {\psi }_{k,2}\hfill & \cdots \hfill & {\psi }_{k,m}\hfill \end{array}{]}^{\mathrm{T}})$ , (25) can be rewritten as

$$\begin{array}{c}(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1}{)}^{-1}({\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1})-{\mathit{H}}_k^{\mathrm{T}}((\frac{{\upsilon }_2}{{\upsilon }_2-2}{\mathit{R}}_k{)}^{-\frac{1}{2}}{)}^{\mathrm{T}}\hfill \\ {\mathit{\psi }}_k(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\mathit{R}}_k{)}^{-\frac{1}{2}}({\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k)=\mathit{0}\hfill \end{array}$$

(27)

The maximum posterior estimation of ${\mathit{x}}_k$ is then obtained by solving the following cost function, which is derived and transformed from (27).

$$\begin{array}{c}{\hat{\mathit{x}}}_{k|k}=\arg min\hfill \\ (\frac{1}{2}||{\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}|{|}_{{(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1})}^{-1}}^2+\frac{1}{2}||{\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k|{|}_{{(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\overline{\mathit{R}}}_k)}^{-1}}^2)\hfill \end{array}$$

(28)

where ${\overline{\mathit{R}}}_k={\mathit{R}}_k^{\frac{1}{2}}{\mathit{\psi }}_k^{-1}({\mathit{R}}_k^{\frac{1}{2}}{)}^{\mathrm{T}}$ . The modified cost function can be treated as the integration of the STF and the maximum correntropy criterion with the updated covariance. This can update the covariance according to the influence of unknown modeling error for robustness design.

The Laplace ℓ₁ distribution design

Based on the above result, the Laplace ℓ₁ distribution is introduced to formulate the proposed L1-RSTF. Firstly, the normalized residual vector is denoted as ${\mathit{\xi }}_k=(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\overline{\mathit{R}}}_k{)}^{-\frac{1}{2}}({\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k)=[\begin{array}{cccc}{\xi }_{k,1}\hfill & {\xi }_{k,2}\hfill & \cdots \hfill & {\xi }_{k,m}\hfill \end{array}{]}^{\mathrm{T}}$ . Then, the new cost function (28) can be simplified by the Laplace ℓ₁ distribution as

$$\begin{array}{c}{\hat{\mathit{x}}}_{k|k}=\arg min\hfill \\ (\frac{1}{2}||{\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}|{|}_{{(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1})}^{-1}}^2+\sqrt{2}||{\mathit{\xi }}_k|{|}_1)\hfill \end{array}$$

(29)

Because the differentiation of the ℓ₁ norm function in (29) does not exist at ${\mathit{\xi }}_k=0$ , it fails to apply the traditional gradient techniques to solve (29) directly.

To solve this above difficulty, an iterative surrogate function is to be presented for re-describing the cost function (29). The majorization minimization (MM) method will be designed further to solve this issue by calculating the minimum value of the surrogate function. Consider a positive function $g(·)$ and a positive scalar $a$ , it follows that

$$\sqrt{g(\mathit{x})}\le \frac{a}{2}+\frac{g(\mathit{x})}{2a}$$

(30)

Moreover, $\sqrt{g(\mathit{x})}=\frac{a}{2}+\frac{g(\mathit{x})}{2a}$ holds only when $g(\mathit{x})=a^2$ . Then, one has

$$\begin{array}{c}||{\mathit{\xi }}_k|{|}_1=\sum _{i=1}^m|{\xi }_{k,i}|\le \sum _{i=1}^m(\frac{{\Xi }_i}{2}+\frac{{\xi }_{k,i}^2}{2{\Xi }_i}),\hfill \\ \forall {\Xi }_i>0,i=1,\cdots ,m\hfill \end{array}$$

(31)

Substituting (31) into (29), it leaves the surrogate function as

$$\{\begin{array}{c}{\hat{\mathit{x}}}_{k/k}=\arg min\mathit{J}({\mathit{x}}_k,\mathit{\Xi })\hfill \\ \multicolumn{1}{c}{\mathit{J}({\mathit{x}}_k,\mathit{\Xi })=\frac{1}{2}||{\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}|{|}_{{(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1})}^{-1}}^2+\frac{1}{2}||{\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k|{|}_{{(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\hat{\mathit{R}}}_k)}^{-1}}^2+\frac{\sqrt{2}}{2}\sum _{i=1}^m{\Xi }_i}\end{array}$$

(32)

where ${\hat{\mathit{R}}}_k=\frac{\sqrt{2}}{2}{\overline{\mathit{R}}}_k^{\frac{1}{2}}\mathit{\Xi }({\overline{\mathit{R}}}_k^{\frac{1}{2}}{)}^{\mathrm{T}}$ and $\mathit{\Xi }=\mathit{diag}([\begin{array}{cccc}{\Xi }_1\hfill & {\Xi }_2\hfill & \cdots \hfill & {\Xi }_m\hfill \end{array}{]}^{\mathrm{T}})$ .

Now, the original optimization problem has been changed into solving the surrogate function (32) via the MM technique. It is described as

$${\hat{\mathit{x}}}_{k|k}^{(c+1)}=\arg {\min }_{{\mathit{x}}_k}\mathit{J}({\mathit{x}}_k,{\mathit{\Xi }}^{(c)})$$

(33)

$${\mathit{\Xi }}^{(c+1)}=\arg {\min }_{\mathit{\Xi }}\mathit{J}({\hat{\mathit{x}}}_{k|k}^{(c+1)},{\mathit{\Xi }}^{(c)})$$

(34)

Aiming of applying the equivalent condition (31), the solution of (34) could be approximated by

$${\Xi }_i^{(c+1)}=\left|{({(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\hat{\mathit{R}}}_k)}^{-\frac{1}{2}}({\mathit{y}}_k-{\mathit{H}}_k{\hat{\mathit{x}}}_{k|k}^{(c)}))}_i\right|+{\epsilon }_i$$

(35)

where ${\epsilon }_i$ is a very small non-zero constant (e.g. 10⁻⁶∼10⁻⁸) to avoid ${\Xi }_i^{(c)}=0$ with $c$ representing the iterations’ number. Generally, a starting weighting matrix ${\mathit{\Xi }}^{(0)}$ could be selected first to obtain ${\hat{\mathit{x}}}_{k|k}^{(1)}$ through (33). Then, the estimation ${\hat{\mathit{x}}}_{k|k}^{(1)}$ could be used to update and obtain a new weighting matrix ${\mathit{\Xi }}^{(1)}$ . In this work, the identity matrix $\mathit{I}$ is treated as the initial weighting matrix for simplification.

Note that the surrogate function (32) has some similarities as the cost function of STF for a given ${\Xi }_i^{(c)}$ . Hence, the function (32) is equivalent to the following form

$$\begin{array}{c}{\hat{\mathit{x}}}_{k|k}^{(c)}=\arg min\hfill \\ (\frac{1}{2}\Vert {\mathit{x}}_k-{\hat{\mathit{x}}}_{k|k-1}{\Vert }_{{(\frac{{\upsilon }_3}{{\upsilon }_3-2}{\mathit{P}}_{k|k-1})}^{-1}}^2+\frac{1}{2}||{\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k|{|}_{{(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\hat{\mathit{R}}}_k)}^{-1}}^2)\hfill \end{array}$$

(36)

It obtains the upper bound of this cost function in the STF framework. Once getting the estimation state ${\hat{\mathit{x}}}_{k|k}^{(c)}$ , one can update and iterate a new weighting matrix ${\mathit{\Xi }}^{(c+1)}$ for obtaining ${\hat{\mathit{x}}}_{k|k}^{(c+1)}$ .

The Laplace ℓ₁ robust Student’s T-filter

Based on section “The Laplace ℓ₁ distribution design,” we are ready to present the L1-RSTF in the following algorithm. The proposed L1-RSTF has a similar structure to the STF in the measurement updating. The major difference is in the state updating. Therefore, the new process of state updating for proposed L1-RSTF is described as follows.

• (1) State updating for L1-RSTF:

The initial value of the iteration and the scalar matrix for the measurement covariance is given by

$${\mathit{\Xi }}^{(0)}=\mathit{I},c=0$$

(37)

$$\begin{array}{c}{\overline{\mathit{P}}}_{k|k-1}^{\mathrm{yy}}=\frac{{\upsilon }_3-2}{{\upsilon }_3}\hfill \\ {\int }_{{\Re }^m}\mathit{h}({\mathit{x}}_k){(\mathit{h}({\mathit{x}}_k))}^{\mathrm{T}}\mathbf{St}({\mathit{x}}_k;{\hat{\mathit{x}}}_{k|k-1},{\mathit{P}}_{k|k-1},{\upsilon }_3)d{\mathit{x}}_k\hfill \\ -\frac{{\upsilon }_3-2}{{\upsilon }_3}{\hat{\mathit{y}}}_{k|k-1}{\hat{\mathit{y}}}_{k|k-1}^{\mathrm{T}}\hfill \end{array}$$

(38)

• (2) Repeating:

$${\mathit{e}}_k=(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\mathit{R}}_k{)}^{-\frac{1}{2}}({\mathit{y}}_k-{\mathit{H}}_k{\mathit{x}}_k)$$

(39)

$${\mathit{\psi }}_k=\mathit{diag}([\begin{array}{cccc}{\psi }_{k,1}({\mathit{e}}_{k,1})\hfill & {\psi }_{k,2}({\mathit{e}}_{k,2})\hfill & \cdots \hfill & {\psi }_{k,m}({\mathit{e}}_{k,m})\hfill \end{array}{]}^{\mathrm{T}})$$

(40)

$${\overline{\mathit{R}}}_k={\mathit{R}}_k^{\frac{1}{2}}{\mathit{\psi }}_k^{-1}({\mathit{R}}_k^{\frac{1}{2}}{)}^{\mathrm{T}}$$

(41)

$${\hat{\mathit{R}}}_k=\frac{\sqrt{2}}{2}{\overline{\mathit{R}}}_k^{\frac{1}{2}}{\mathit{\Xi }}^{(c)}({\overline{\mathit{R}}}_k^{\frac{1}{2}}{)}^{\mathrm{T}}$$

(42)

$${\mathit{P}}_{k|k-1}^{\mathrm{yy}}={\overline{\mathit{P}}}_{k|k-1}^{\mathrm{yy}}+\frac{{\upsilon }_2({\upsilon }_3-2)}{({\upsilon }_2-2){\upsilon }_3}{\hat{\mathit{R}}}_k$$

(43)

$${\mathit{K}}_k^{(c)}={\mathit{P}}_{k|k-1}^{\mathrm{xy}}({\mathit{P}}_{k|k-1}^{\mathrm{yy}}{)}^{-1}$$

(44)

$${\hat{\mathit{x}}}_{k+1|k+1}^{(c)}={\hat{\mathit{x}}}_{k+1|k}+{\mathit{K}}_{k+1}^{(c)}({\mathit{y}}_{k+1}-{\hat{\mathit{y}}}_{k+1|k})$$

(45)

$${\mathrm{\Delta }}_k=\sqrt{{({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1})}^{\mathrm{T}}{\mathit{P}}_{k|k-1}^{\mathrm{yy}}({\mathit{y}}_k-{\hat{\mathit{y}}}_{k|k-1})}$$

(46)

$${\Xi }_i^{(c+1)}=\left|{({(\frac{{\upsilon }_2}{{\upsilon }_2-2}{\hat{\mathit{R}}}_k)}^{-\frac{1}{2}}({\mathit{y}}_k-{\mathit{H}}_k{\hat{\mathit{x}}}_{k|k}^{(c)}))}_i\right|+{\epsilon }_i$$

(47)

$${\mathit{\Xi }}^{(c+1)}=\mathit{diag}([\begin{array}{cccc}{\Xi }_1^{(c+1)}\hfill & {\Xi }_2^{(c+1)}\hfill & \cdots \hfill & {\Xi }_m^{(c+1)}\hfill \end{array}{]}^{\mathrm{T}})$$

(48)

$$c=c+1$$

(49)

Until the desired convergence performance index is satisfied.

Up to now, the final estimation for the states and their scalar matrix of the covariance can be determined as

$${\hat{\mathit{x}}}_{k+1|k+1}={\hat{\mathit{x}}}_{k+1|k+1}^{(c)}$$

(50)

$${\mathit{P}}_{k+1|k+1}=\frac{({\upsilon }_3-2)({\upsilon }_3+{\mathrm{\Delta }}_k^2)}{{\upsilon }_3({\upsilon }_3+m-2)}({\mathit{P}}_{k+1|k}-{\mathit{K}}_{k+1}^{(c)}{\mathit{P}}_{k+1|k}^{\mathrm{yy}}{\mathit{K}}_{k+1}^{(c)})$$

(51)

It can be known that the proposed L1-RSTF is similar to the standard STF via the Laplace ℓ₁ theory. However, it is difficult to choose the convergence criterion. Hence, one can determine the number of iterations according to the precision requirements.

Numerical simulation

In this section, the application of the L1-RSTF developed in section “Development of the Laplace ℓ₁ robust Student’s T-filter” to achieve attitude estimation for small satellites will be done.

Mathematical model of a rigid satellite

The unit quaternion $\mathit{q}=[\begin{array}{cc}{q}_1\hfill & {\mathit{q}}_{24}\hfill \end{array}{]}^{\mathrm{T}}=[\begin{array}{cccc}{q}_1\hfill & q_2\hfill & q_3\hfill & q_4\hfill \end{array}{]}^{\mathrm{T}}$ , ${\mathit{q}}_{24}=[\begin{array}{ccc}{q}_2\hfill & q_3\hfill & q_4\hfill \end{array}{]}^{\mathrm{T}}$ , is adopted to represent the attitude with respect to the inertial frame. Then, the kinematics of a rigid satellite is given by

$$\stackrel{·}{\mathit{q}}=\frac{1}{2}\mathit{\Phi }(\mathit{q})\mathit{\omega }$$

(52)

where $\mathit{\omega }=[\begin{array}{ccc}{\omega }_1\hfill & {\omega }_2\hfill & {\omega }_3\hfill \end{array}{]}^{\mathrm{T}}$ is the angular velocity in the body frame with respect to the inertial frame. $\mathit{\Phi }(\mathit{q})$ is given by

$$\mathit{\Phi }(\mathit{q})=\left[\begin{array}{c}-{\mathit{q}}_{24}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{q_1\mathit{I}+{({\mathit{q}}_{24})}^{\times }}\end{array}\right]$$

(53)

For small satellites, their attitude dynamic model is established as

$$\stackrel{·}{\mathit{\omega }}={\mathit{J}}^{-1}(\mathit{\tau }-{\mathit{\omega }}^{\times }\mathit{J}\omega +\mathit{d})$$

(54)

where $\mathit{\tau }\in {\Re }^3$ is the total torque and $\mathit{J}\in {\Re }^{3\times 3}$ is the inertia of the satellite. Based on (50), the attitude dynamics without any gyro measurement could be written as:

$$\stackrel{·}{\mathit{x}}=\left[\begin{array}{c}0.5\mathit{\Phi }(\mathit{q})\mathit{\omega }\hfill \\ \multicolumn{1}{c}{{\mathit{J}}^{-1}(\mathit{\tau }-{\mathit{\omega }}^{\times }\mathit{J}\omega )+{\mathit{J}}^{-1}\mathit{d}}\end{array}\right]+\mathit{w}$$

(55)

where $\mathit{w}$ denotes the heavy-tailed noise and $\mathit{x}=[\begin{array}{cc}{\mathit{q}}^{\mathrm{T}}\hfill & {\mathit{\omega }}^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ .

Moreover, considering cost limitation, a three-axis magnetometer and a sun sensor are fixed to measure the states of the attitude control system. For these sensors, the measurement model could be given by

$$\mathit{y}=\mathit{h}(\mathit{X})=\left[\begin{array}{c}{\mathit{C}}_i^b(\mathit{q}){\mathit{S}}_1\hfill \\ \multicolumn{1}{c}{{\mathit{C}}_i^b(\mathit{q}){\mathit{B}}_2}\end{array}\right]+\left[\begin{array}{c}\mathrm{\Delta }{\nu }_1\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{\nu }_2}\end{array}\right]$$

(56)

where ${\mathit{C}}_i^b\in {\Re }^{3\times 3}$ represents the transition matrix from the inertial frame to the body-fixed frame, and it describes the attitude information. ${\mathit{S}}_1\in {\Re }^3$ and ${\mathit{B}}_2\in {\Re }^3$ are the sun vector and the magnetic field, respectively. $\mathrm{\Delta }{\nu }_1\in {\Re }^3$ and $\mathrm{\Delta }{\nu }_2\in {\Re }^3$ are the heavy-tailed measurement noise of the sensors, respectively.

To this end, the satellite’s attitude control system described by (51) and (52) can be written into the form of (1) after discretization. Consequently, the proposed L1-RSTF can be applied to estimate the attitude state with $\hat{\mathit{x}}=[\begin{array}{cc}{\hat{\mathit{q}}}^{\mathrm{T}}\hfill & {\hat{\mathit{\omega }}}^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}$ denoting the estimation of the state $\mathit{x}$ .

Simulation results

To validate the performance of the proposed L1-RSTF, simulations are conducted by comparing it with the CKF. In the simulation, the initial states of the satellites are $\mathit{q}(0)={\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}$ and $\mathit{\omega }(0)=[\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \end{array}{]}^{\mathrm{T}}$ deg/s. The states’ initial estimation is set to be $\hat{\mathit{q}}(0)=\frac{\sqrt{3}}{3}{\left[\begin{array}{cccc}\sqrt{3}\cos (0.5\Theta )\hfill & \sin (0.5\Theta )\hfill & \sin (0.5\Theta )\hfill & \sin (0.5\Theta )\hfill \end{array}\right]}^{\mathrm{T}}$ , $\Theta =5°$ , and $\hat{\mathit{\omega }}(0)=[\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \end{array}{]}^{\mathrm{T}}$ deg/s. The inertia of the satellite has a value of $\mathit{J}=\mathit{diag}([\begin{array}{ccc}49.96\hfill & 55.40\hfill & 63.0\hfill \end{array}{]}^{\mathrm{T}})$ kg m². The processing and the measurement noises are chosen as

$$\{\begin{array}{c}{\mathit{w}}_k=\{\begin{array}{c}\mathbf{\chi }(0,\mathit{Q}),\hfill \\ \multicolumn{1}{c}{\mathbf{\chi }(0,100\mathit{Q}),}\end{array}\begin{array}{cc}\mathrm{w}.p.\hfill & p1\hfill \\ \multicolumn{1}{c}{\mathrm{w}.p.}& 1-p1\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}{\mathit{v}}_k=\{\begin{array}{c}\mathbf{\chi }(0,\mathit{R}),\hfill \\ \multicolumn{1}{c}{\mathbf{\chi }(0,100\mathit{R}),}\end{array}\begin{array}{cc}\mathrm{w}.p.\hfill & p2\hfill \\ \multicolumn{1}{c}{\mathrm{w}.p.}& 1-p2\hfill \end{array}\hfill \end{array}}\\ \multicolumn{1}{c}{\mathit{Q}=\mathit{diag}([\begin{array}{cccccc}10\hfill & 10\hfill & 10\hfill & 1\hfill & 1\hfill & 1\hfill \end{array}{]}^{\mathrm{T}})\times {10}^{-11}}\\ \multicolumn{1}{c}{{\mathit{R}}_B=(5\times {10}^{-4})\mathit{I}}\\ \multicolumn{1}{c}{\begin{array}{c}{\mathit{R}}_S=\frac{1}{3}{\sin }^2(\frac{1}{900}\pi )\mathit{I}\hfill \end{array}}\end{array}$$

(57)

where $\mathit{Q}$ represent the covariance of the process noise with $p1=0.9$ . $\mathit{R}$ is the measurement noise’s covariance with $p2=0.9$ . The measurement accuracy of the magnetometer and the sun sensor is 50 nT and $0.2$ degrees, respectively. Their covariances are represented by ${\mathit{R}}_B$ and ${\mathit{R}}_S$ , respectively. Based on Cao et al., ³³ the 13th-order Gaussian sphere harmonic function of IGRF is considered to simulate the real ${\mathit{y}}_B$ . The 7th-order of IGRF is applied to simulate ${\mathit{B}}_2$ in the satellite system. The measurement error between the 13th-order and 7th-order of the IGRF is about 208 nT. This error is also considered. In the simulation, the following two scenarios are considered.

Remark: Let $\tilde{\mathit{x}}=[\begin{array}{cc}{\mathit{\omega }}_{\mathrm{error}}^{\mathrm{T}}\hfill & {\tilde{\mathit{q}}}^{\mathrm{T}}\hfill \end{array}{]}^{\mathrm{T}}=\mathit{x}-\hat{\mathit{x}}$ denote the state estimation error, where $\tilde{\mathit{q}}=\mathit{q}-\hat{\mathit{q}}$ and ${\mathit{\omega }}_{\mathrm{error}}=[\begin{array}{ccc}{W}_{\mathrm{xerror}}\hfill & W_{\mathrm{yerror}}\hfill & W_{\mathrm{zerror}}\hfill \end{array}{]}^{\mathrm{T}}=\mathit{\omega }-\hat{\mathit{\omega }}$ . In satellite engineering, the unit quaternion $\mathit{q}$ can be mathematically transformed into the Euler attitude angles to explicitly describe the satellite’s attitude, that is,

$$\mathit{\Theta }=[\begin{array}{ccc}\varphi \hfill & \theta \hfill & \psi \hfill \end{array}{]}^{\mathrm{T}}={\mathit{f}}_t(\mathit{q})$$

(58)

where ${\mathit{f}}_t(·)$ is the transformation function. ³⁴ $\varphi$ , $\theta$ , and $\psi$ are the roll, the pitch, the yaw attitude Euler angles with respect to $\mathit{q}$ , respectively. Hence, the attitude estimation error $\tilde{\mathit{q}}$ can be transformed into the attitude Euler angles estimation error denoted as ${\mathit{\Theta }}_{\mathrm{error}}=[\begin{array}{ccc}{\varphi }_{\mathrm{error}}\hfill & {\theta }_{\mathrm{error}}\hfill & {\psi }_{\mathrm{error}}\hfill \end{array}{]}^{\mathrm{T}}$ .

(1) The first scenario: This scenario is to validate the developed filter’s effectiveness in handling the heavy-tailed noise. Figures 1 and 2 show the resulted estimation errors. Because the CKF is widely used in practice, the 3rd-CKF and the 5th-CKF are also employed to compare with the proposed L1-RSTF. It is found in Figures 1 and 2 that the convergence rate of L1-RSTF is much faster and more stable than the CKFs without the obvious fluctuations. Although the advantage of L1-RSTF is not very obvious in the angular velocity estimation, this method still has a certain precision superiority. Moreover, there are a few estimated mutations for the L1-RSTF, which is mainly due to the mutation of heavy-tailed distribution with a small probability. However, due to the theoretical advantages of proposed L1-RSTF, the estimation error divergence can be rapidly suppressed and return to the normal estimate trajectory.

Figure 1.

The attitude estimation error in the first scenario.

Figure 2.

The angular velocity estimation error in the first scenario.

Table 1 lists the RMSE results of the estimation error between the 3rd-CKF, the 5th-CKF, and the L1-RSTF. It is listed that the developed L1-RSTF has higher estimation precision than the 3rd-CKF and the 5th-CKF, which is especially true in attitude estimation. In the aspect of angular velocity, the theoretical advantage of L1-RSTF is not obvious, but it still has the highest accuracy. Regarding the computational burden, the 3rd-CKF and the 5th-CKF need 12 and 73 sampling points to calculate and propagate the state estimation, respectively. However, the proposed L1-RSTF only needs 12 sampling points and the number of iterations is five. Therefore, the L1-RSTF also has obvious advantages in computational complexity.

Table 1.

The RSME comparison in the first scenario.

RSME	$\varphi$ (degrees)	$\theta$ (degrees)	$\psi$ (degrees)	${\omega }_x$ (rad/s)	${\omega }_y$ (rad/s)	${\omega }_z$ (rad/s)
3rd-CKF	0.0101	0.0154	0.0111	0.000349	0.000309	0.000302
5th-CKF	0.0091	0.0139	0.0100	0.000314	0.000279	0.000272
L1-RSTF	0.0039	0.0047	0.0052	0.000280	0.000280	0.000255

(2) The second scenario: The attitude dynamics of small satellites are characterized by modeling error. Hence, the modeling error in the attitude dynamics is considered. It is represented by $\mathit{d}=\mathrm{\Delta }{\mathit{N}}_c+\mathrm{\Delta }{\mathit{N}}_e$ , where $\mathrm{\Delta }{\mathit{N}}_c$ and $\mathrm{\Delta }{\mathit{N}}_e$ are assumed as

$${\mathit{N}}_c=[\begin{array}{ccc}0.0001\hfill & 0.0002\hfill & 0.0004\hfill \end{array}{]}^{\mathrm{T}}\mathrm{Nm}$$

(59)

$$\mathrm{\Delta }{\mathit{N}}_e=[\begin{array}{ccc}0.0003\cos (10{\omega }_{°}t)\hfill & 0.0002\cos (10{\omega }_{°}t)\hfill & -0.0001\cos (10{\omega }_{°}t)\hfill \end{array}{]}^{\mathrm{T}}\mathrm{Nm}$$

(60)

Therefore, the modeling errors and the heavy-tailed noise (57) are both taken into consideration in this case and introduced into the simulation system, which can be used to test the robust ability of presented L1-RSTF to the modeling error and the noise.

Figures 3 and 4 show that the estimation accuracy of the second scenario is inferior to the accuracy in the first scenario. This coincides with the preceding theoretical result. Especially for the CKFs, the estimated convergence curve deviates seriously from the zero-value due to the existence of model error in Figure 3. However, the L1-RSTF still fluctuates slightly around zero and maintains a high precision compared with the CKFs. The angular velocity estimate accuracy of L1-RSTF also has better convergence properties than the CKFs, as we can see in Figure 4.

Figure 3.

The attitude estimation error in the second scenario.

Figure 4.

The angular velocity estimation error in the second scenario.

Moreover, the comparison results listed in Table 2 show that the estimation accuracy obtained from the 3rd-CKF and the 5th-CKF is severely decreased due to the modeling error. The reason is that the nonlinearity and the uncertainty of the satellite system are increased by modeling error and noise. Only the proposed L1-RSTF maintains acceptable high precision. This validates the superiority of L1-RSTF.

Table 2.

The RSME comparison in the second scenario.

RSME	$\varphi$ (degrees)	$\theta$ (degrees)	$\psi$ (degrees)	${\omega }_x$ (rad/s)	${\omega }_y$ (rad/s)	${\omega }_z$ (rad/s)
3rd-CKF	0.0527	0.1081	0.0748	0.00120	0.000924	0.001100
5th-CKF	0.0475	0.0973	0.0674	0.00100	0.000832	0.000991
L1-RSTF	0.0074	0.0101	0.0103	0.000613	0.000467	0.000608

The above simulation results verify that the presented scheme is robust to the heavy-tailed noise and the modeling error with low computation required. For the CKFs, only the angular velocity estimation can get acceptable performance. It is concluded that the proposed scheme is less dependent on the model of the uncertainty and the noise. This lets the designed filter be superior to the traditional Kalman filter because the latter relies heavily on the known noise level. Therefore, it can be anticipated that this method has important application value and will be more widely used in the future.

Conclusion

This work presented a Laplace ℓ₁ robust Studen’s T-filter to solve the attitude estimation problem of small satellites subject to heavy-tailed noise and modeling errors. The maximum correntropy criterion and the Laplace distribution were used to model the measurement process. The majorization minimization technique was applied to solve the maximum posterior estimation via an iterative algorithm. Numerical simulation with the application of the proposed L1-RSTF was conducted by comparing it with the 3rd-CKF and the 5th-CKF. The effectiveness of the proposed approach in both the accuracy of system estimation and the system robustness with less calculation burden was validated. As one of future work, using the estimated states provided by the proposed filtering approach, the attitude controller design should be carried out to ensure high-accuracy attitude control. This could be done by invoking the robust methods reported in.^35–38