Event-Driven Maximum Correntropy Filter Based on Cauchy Kernel for Spatial Orientation Using Gyros/Star Sensor Integration

Abstract

Gyros/star sensor integration provides a potential method to obtain high-accuracy spatial orientation for turntable structures. However, it is subjected to the problem of accuracy loss when the measurement noises become non-Gaussian due to the complex spatial environment. This paper presents an event-driven maximum correntropy filter based on Cauchy kernel to handle the above problem. In this method, a direct installation mode of gyros/star sensor integration is established and the associated mathematical model is derived to improve the turntable’s control stability. Based on this, a Cauchy kernel-based maximum correntropy filter is developed to curb the influence of non-Gaussian measurement noise for enhancing the gyros/star sensor integration’s robustness. Subsequently, an event-driven mechanism is constructed based on the filter’s innovation information for further reducing the unnecessary computational cost to optimize the real-time performance. The effectiveness of the proposed method has been validated by simulations of the gyros/star sensor integration for spatial orientation. This shows that the proposed filtering methodology not only has strong robustness to deal with the influence of non-Gaussian measurement noise but can also achieve superior real-time spatial applications with a small computational cost, leading to enhanced performance for the turntable’s spatial orientation using gyros/star sensor integration.

1. Introduction

The turntable structure can realize the orientation motion of the payload (camera, radar and so on) independently of the spacecraft platform. It plays a vital role for the space tracking and detection system in the aerospace field [1,2]. A two-dimensional turntable is one of the most widely used structures in the above systems. It consists of two rotating platforms, in which their shafting is perpendicular to each other [1,3]. The highly precise orientation measurement of the two-dimensional turntable is an important prerequisite for achieving high stability and enhancing the stability of the space tracking and detection system.

Spatial orientation accuracy refers to the error between the ideal spatial orientation and the actual spatial orientation associated with the inner frame of the turntable after the rotation and transformation of each axis. Due to the inevitable control error, the installation error and other factors in the practical application of the system, there is a deviation between the actual orientation and the ideal orientation of the turntable, which significantly affects the control stability of the turntable [4,5]. Therefore, highly precise orientation measurements of the turntable are particularly important.

Currently, in spatial orientation applications, the star sensor is a commonly used measurement device that can provide high-precision attitude and orientation information [6,7]. However, the operating frequency of the star sensor is quite low and it is easily limited by the time and environmental conditions [7,8]. In contrast, gyros, as an autonomous inertial measurement method, can continuously output high-frequency angular/attitude information with good real-time performance [9,10]. However, gyros have the weakness of error drifts, leading to poor long-term attitude/orientation accuracy [11]. It can be seen that the above two sensors have a good complementary characteristic in terms of attributes, so the gyros/star sensor integration scheme can use the high-precision attitude or orientation information obtained from the star sensor to compensate for the drift of the gyros, achieving high-precision turntable orientation information [12,13].

When the gyros/star sensor integration is implemented for the two-dimensional turntable’s spatial orientation, there are usually two schemes: one is the strap-down installation mode and the other is the direct installation mode [13,14]. In the former, both the gyros and star sensor are installed on the fixed base of the turntable, while in the latter, the gyros are installed in the inner frame of the two-dimensional turntable and firmly connected with the load camera/radar, and the star sensor is installed on the fixed base of the turntable. It should be noted that in many applications, due to the limited space, the star sensor is hardly installed in the inner frame of the turntable because of its large size. In the strap-down installation mode, the turntable’s control stability is easily affected by the turntable error since the gyros/star sensor integration cannot directly perceive the turntable angle. Thus, the direct installation mode has gained increasing attention by scholars in recent years. The research of this paper is also based on the direct installation mode of gyros/star sensor integration.

Nowadays, the data fusion of gyros/star sensor integration is mostly achieved based on Kalman filtering [15,16]. Under the assumption of linear and Gaussian noise models, Kalman filtering can obtain the optimal estimation of a system state [15,16]. However, due to the influence of the spatial light environment, the measurement of star sensors will be contaminated. This makes the measurement noise of gyros/star sensor integration no longer follow a Gaussian distribution, resulting in non-Gaussian noise. In this situation, the performance of Kalman filtering will be significantly deteriorated [17,18].

Research efforts have been dedicated to dealing with non-Gaussian noise involved in the dynamic system. Farahmand et al. proposed a doubly robust smoothing method to provide double robustness that enhances the accuracy and reliability of estimations, particularly in the presence of bias measurement noise covariance and model misspecification [18]. However, the smoothing method is not suitable for a real-time application. Further, Soken et al. established a robust Kalman filter by using a scalar factor to adjust the measurement noise covariance for curbing the influences of inaccurate measurement noise covariances [19]. However, its improvement is very limited since the scalar factor is chosen empirically. Particle filter (PF) and Gaussian sum filter (GSF) are also two typical forms of non-Gaussian filtering [20,21]. Kim et al. developed PF to further extend the non-Gaussian capabilities of filtering by employing sampling and resampling techniques to approximate posterior distributions [20]. Sorenson et al. constructed a GSF by using random sampling to obtain the global approximation of the system state’s posterior probability density, which is not limited to Gaussian noise [21]. However, global approximation requires a large amount of computation, making the above filters difficult to apply for the spatial orientation [22].

As an emerging filtering theory, information theoretic learning has garnered attention for its optimization criteria, which rely on information theoretic measures such as maximum correntropy criterion (MCC) derived directly from data, instead of conventional second-order statistical metrics [23,24]. It exhibits an excellent robustness to handle non-Gaussian noise since it uses the signals’ high-order moments rather than the second-order moments in the minimum mean square error principle. Chen et al. [23] first established the maximum correntropy Kalman filter by using the maximum correntropy criterion. It has superior performance to resist the influences of outliers and non-Gaussian noise for a linear system. Compared to PF and GSF, the filter based on the MCC shows a quite low computational burden. Further, Mohammadi et al. [24] presented a maximum–minimum correntropy criterion for robust and stable gene selection. Gao et al. [25] developed the theory of MCC to handle outliers and non-Gaussian noise to improve the performance of spectral redshift estimation for spectral redshift navigation. In the traditional maximum correntropy (MC) filter, the Gaussian kernel is always used to define the distance between distinct vectors. Nevertheless, it may not always be an optimal selection for the kernel function. An issue arises when the measurement is affected by multi-dimensional non-Gaussian noise, such as heavy-tailed noise in real-world applications. In such cases, the aforementioned MC filters may fail due to the occurrence of singular matrices, causing certain characteristics of the noise to be overlooked and resulting in reduced estimation accuracy [26,27].

To address the drawbacks of Gaussian kernel-based MC filters, a novel Cauchy kernel-based MCC filter is proposed [27,28]. It uses the Cauchy kernel instead of the original Gaussian kernel to describe the distance between different vectors. The Cauchy kernel originated from the Cauchy distribution. Its function curve has a long tail, indicating that the kernel function is very broad and can be used for data with high dimensions. Because of the above reason, the Cauchy kernel-based MC filter shows a superior performance when handling multi-dimensional non-Gaussian noise. It also behaves in a consistently stable manner when a different Cauchy kernel bandwidth is selected [28,29]. Unfortunately, the existing MC filters or Cauchy kernel-based MC filters mostly utilize the fixed-point iteration method to update the estimates of the state, which will burden the computational efficiency [27,30]. For the application of spatial orientation, the real-time performance is a very important property index. Thus, if the MCC is used for the whole filtering process of gyros/star sensor integration, it will increase the unnecessary computational burden, leading to a deteriorated real-time performance. This is because in the time point without measurement fault, i.e., non-Gaussian measurement noise does not exist, the Kalman filter can obtain the optimal estimation result with a lower computational cost than MC filters, while MC filters need to use the fixed-point iteration to update the estimates of the state, requiring a large computational cost.

To address the above problem for a spatial orientation application, this paper presents a novel method of a Cauchy kernel-based maximum correntropy filter that is event-driven to curb the influence of non-Gaussian measurement noise on the solution of gyros/star sensor integration. It first establishes a direct installation mode of gyros/star sensor integration to achieve a high turntable control stability and then derives its associated mathematical model. Further, an event-driven maximum correntropy filter based on the Cauchy kernel is proposed to handle the non-Gaussian noise involved in the measurement of gyros/star sensor integration caused by the complex spatial environment. It also constructs an event-driven mechanism according to the innovation information to reduce the redundant computational cost for the optimization of the filter’s real-time performance. The proposed methodology can achieve robustness to curb the influence of non-Gaussian measurement noise with a superior computational cost compared to the existing MC filters. Simulations and a comparison analysis based on the gyros/star sensor integration for a two-dimensional turntable’s spatial orientation have been carried out to comprehensively evaluate the effectiveness of the proposed filtering method.

2. Model of Gyros/Star Sensor Integration for Spatial Orientation

The direct installation mode is used for the integration of the gyros and star sensor, which is described in Figure 1. It can be seen from Figure 1 that the gyros are installed in the inner frame of the two-dimensional turntable and firmly connected with the load camera/radar, and the star sensor is installed on the fixed base of the turntable.

Based on the direct installation mode, the framework of the gyros/star sensor integration used for the spatial orientation is illustrated in Figure 2, which has two loops: a filtering loop and error calculation loop. In the filtering loop, the gyros’ measurement information, the turntable control quantities, and the attitude of the star sensor serve as the input of filter for the error state estimation of gyros/star sensor integration. Then, the estimated gyros’ attitude error is fed back to correct the gyros and the gyros’ attitude is updated for the next time step. In the error calculation loop, the orientation error between the tracking vector (the normalized vector between the spacecraft and target) and the attitude reference of the spacecraft is calculated by using the turntable’s control quantities and the spacecraft’s attitude at each time step to reflect the tracking accuracy.

It should be noted that as shown in Figure 1 and Figure 2, the star sensor realizes the three-dimensional attitude information of turntable base. Since the rotation angles of the turntable’s inner frame and outside frame are known through the pitch and azimuth code discs, the three-dimensional attitude information of the turntable base can be transformed into the turntable’s inner frame (E-frame) and inertial frame (I-frame). Based on this attitude information, the measurements and measurement model for the gyros/star sensor integration are established. Further, combining the process model with the measurement model, gyros’ platform misalignment angle can be estimated through the filtering process. This estimation will be fed back to the gyros to correct the attitude error. Thus, the star sensor can restrain and correct the gyros’ drift and further enhance the turntable’s spatial orientation accuracy.

Based on the above, the model of gyros/star sensor integration can be further derived and established, which includes the process and measurement models.

2.1. Process Model

Denote the inertial frame as the I-frame, body frame as the B-frame, turntable’s base frame as the G-frame, turntable’s outside frame as the A-frame, turntable’s inner frame as the E-frame and gyros’ composite frame as the P-frame.

The I-frame is chosen as the navigation frame; thus, the angular relation between the I-frame and P-frame can be expressed by the attitude matrix [31,32]

$${\dot{\mathit{C}}}_P^I={\mathit{C}}_P^I\left({\mathit{\omega }}_{IP}^P\times \right)$$

(1)

where ${\mathit{C}}_P^I$ is the attitude matrix from the P-frame and I-frame, ${\dot{\mathit{C}}}_P^I$ is the differential form of ${\mathit{C}}_P^I$, ${\mathit{\omega }}_{IP}^P$ represents the gyros’ output, and $\left({\mathit{\omega }}_{IP}^P\times \right)$ stands for the skew symmetric matrix of ${\mathit{\omega }}_{IP}^P$, which can be calculated as [33]

$$\left({\mathit{\omega }}_{IP}^P\times \right)=\left[\begin{array}{ccc}0& -{\mathit{\omega }}_{IPz}^P& {\mathit{\omega }}_{IPy}^P\\ {\mathit{\omega }}_{IPz}^P& 0& -{\mathit{\omega }}_{IPx}^P\\ -{\mathit{\omega }}_{IPy}^P& {\mathit{\omega }}_{IPx}^P& 0\end{array}\right]$$

(2)

where ${\mathit{\omega }}_{IPi}^P$ ($i=x,y,z$) represents the components of ${\mathit{\omega }}_{IP}^P$ in the directions of x, y and z.

Considering the measurement error, the above attitude matrix should be written as ${\mathit{C}}_P^{I^{\prime }}$, and the $I^{\prime }$-frame and I-frame have a platform misalignment angle. Thus, we have the following relationship [32,34]

$${\mathit{C}}_P^{I^{\prime }}={\mathit{C}}_I^{I^{\prime }}{\mathit{C}}_P^I=\left(\mathbf{I}-{\mathit{\varphi }}^I\times \right){\mathit{C}}_P^I$$

(3)

where ${\mathit{\varphi }}^I$ stands for the platform misalignment angle.

Further, according to the definition of platform misalignment angle ${\mathit{\varphi }}^I$ in [35], it has a differential relation with the gyros’ zero-bias error, which can be described as follows. The specifics can be seen in [35].

$${\dot{\mathit{\varphi }}}^I=-{\mathit{C}}_P^I{\mathit{\epsilon }}^P$$

(4)

where ${\mathit{\epsilon }}^P$ is the gyros’ zero-bias error.

It can be seen that the attitude update of the gyros diverges with time due to the existence of gyros’ zero-bias error; thus, we can use the star sensor’s output to estimate and correct the gyros’ attitude errors.

Denote the system state of gyros/star sensor integration as

$$\mathit{X}={\left[\begin{array}{cc}{({\mathit{\varphi }}^I)}^{\mathrm{T}}& {({\mathit{\epsilon }}^P)}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$$

(5)

Based on (5), the process model of gyros/star sensor integration can be described as

$$\dot{\mathit{X}}=\left[\begin{array}{cc}{\mathbf{0}}_{3\times 3}& -{\mathit{C}}_P^I\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\end{array}\right]\mathit{X}+\left[\begin{array}{c}{\mathit{w}}_P\\ {\mathbf{0}}_{3\times 1}\end{array}\right]$$

(6)

Discretizing (6), we can obtain

$${\mathit{X}}_k=\left[\begin{array}{c}{\mathit{\varphi }}_k^I\\ {\mathit{\epsilon }}_k^P\end{array}\right]=\left[\begin{array}{cc}{\mathit{I}}_{3\times 3}& -{\mathit{C}}_P^I\\ {\mathbf{0}}_{3\times 3}& {\mathit{I}}_{3\times 3}\end{array}\right]\left[\begin{array}{c}{\mathit{\varphi }}_{k-1}^I\\ {\mathit{\epsilon }}_{k-1}^P\end{array}\right]+\left[\begin{array}{c}{\mathit{I}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}\end{array}\right]{\mathit{W}}_{k-1}$$

(7)

2.2. Measurement Model

In the direct installation mode, there is a turntable angle between the gyros and the star sensor. Thus, according to Figure 1, the error matrix can be constructed based on the output of the star sensor

$${\mathit{M}}_{Zb}={\tilde{\mathit{C}}}_I^{P_i}{\tilde{\mathit{C}}}_B^I{\mathit{C}}_G^B{\mathit{C}}_A^G{\mathit{C}}_E^A={\tilde{\mathit{C}}}_E^{P_i}$$

(8)

where $\mathit{C}$ is the rotation matrix from one frame to another frame.

In the real application, since the installation angle can be calibrated to compensate for its error, the E-frame can be seen coincident with the P-frame. Further, by using the Rodrigues formula [34], (8) can be further written as

$$\begin{array}{l}{\mathit{M}}_{Zb}={\tilde{\mathit{C}}}_E^{P_i}={\tilde{\mathit{C}}}_P^{P_i}=({\tilde{\mathit{C}}}_{P_i}^P{)}^{\mathrm{T}}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}={\left(\mathit{I}+{\displaystyle \frac{\sin {\varphi }^P}{{\varphi }^P}}({\mathit{\varphi }}^P\times )+{\displaystyle \frac{1-\cos {\varphi }^P}{{\left({\varphi }^P\right)}^2}}{\left({\mathit{\varphi }}^P\right)}^2\right)}^{\mathrm{T}}\approx {\left(\mathit{I}+({\mathit{\varphi }}^P\times )\right)}^{\mathrm{T}}=\left(\mathit{I}-({\mathit{\varphi }}^P\times )\right)\end{array}$$

(9)

where ${\mathit{\varphi }}^P={[\begin{array}{ccc}{\mathit{\varphi }}_x^P& {\mathit{\varphi }}_y^P& {\mathit{\varphi }}_z^P\end{array}]}^{\mathrm{T}}$ is the platform misalignment angle between the P-frame and P_i-frame computed by the output of star sensor, which can be seen as a small amount; ${\varphi }^P=\left|{\mathit{\varphi }}^P\right|$ represents its modulus.

According to the definition of the skew symmetric matrix [33], we have

$$\left({\mathit{\varphi }}^P\times \right)=\left[\begin{array}{ccc}0& -{\mathit{\varphi }}_z^P& {\mathit{\varphi }}_y^P\\ {\mathit{\varphi }}_z^P& 0& -{\mathit{\varphi }}_x^P\\ -{\mathit{\varphi }}_y^P& {\mathit{\varphi }}_x^P& 0\end{array}\right]$$

(10)

Thus, expanding (9) yields

$${\mathit{M}}_{Zb}=\mathit{I}-\left[\begin{array}{ccc}0& -{\mathit{\varphi }}_z^P& {\mathit{\varphi }}_y^P\\ {\mathit{\varphi }}_z^P& 0& -{\mathit{\varphi }}_x^P\\ -{\mathit{\varphi }}_y^P& {\mathit{\varphi }}_x^P& 0\end{array}\right]=\left[\begin{array}{ccc}1& {\mathit{\varphi }}_z^P& -{\mathit{\varphi }}_y^P\\ -{\mathit{\varphi }}_z^P& 1& {\mathit{\varphi }}_x^P\\ {\mathit{\varphi }}_y^P& -{\mathit{\varphi }}_x^P& 1\end{array}\right]$$

(11)

Then, the measurement of gyros/star sensor integration can be constructed as

$${\mathit{Z}}_k=-{\displaystyle \frac{1}{2}}\left[\begin{array}{c}{M}_{32}-M_{23}\\ M_{13}-M_{31}\\ M_{21}-M_{12}\end{array}\right]=\left[\begin{array}{c}{\mathit{\varphi }}_x^P\\ {\mathit{\varphi }}_y^P\\ {\mathit{\varphi }}_z^P\end{array}\right]={\mathit{\varphi }}^P$$

(12)

where $M_{ij}$ represents the ith row and jth column element of the matrix ${\mathit{M}}_{Zb}$.

Using the relationship ${\mathit{\varphi }}^P={\mathit{C}}_I^P{\mathit{\varphi }}^I$, the measurement model can be established by

$${\mathit{Z}}_k={\displaystyle \frac{1}{2}}\left[\begin{array}{c}{M}_{32}-M_{23}\\ M_{13}-M_{31}\\ M_{21}-M_{12}\end{array}\right]=\left[\begin{array}{cc}{\mathit{C}}_I^P& {\mathbf{0}}_{3\times 3}\end{array}\right]{\mathit{X}}_k+{\mathit{V}}_k$$

(13)

where ${\mathit{V}}_k$ is the measurement noise.

Thus, the model of gyros/star sensor integration for spatial orientation, including the process and measurement models, is achieved.

3. Event-Driven Maximum Correntropy Filter Based on Cauchy Kernel

As discussed previously, the star sensor’s measurement is subject to the spatial environment’s interference due to low light intensity. This will make the measurement noise of gyros/star sensor integration in (12) no longer adhere to a Gaussian distribution. In this situation, the traditional Kalman filtering has no ability to deal with the non-Gaussian measurement noise, leading to divergent filtering results. This may cause ineffective correction for the gyros by using the star sensor. Thus, this section establishes an event-driven maximum correntropy filter based on the Cauchy kernel to address the above problem.

3.1. Maximum Correntropy Filter Based on Cauchy Kernel

3.1.1. Maximum Correntropy Criterion

In information theory, correlation entropy is often defined as a statistic that describes the similarity between two random variables. Cross-correlation entropy (correntropy) is an extension of the correlation between random variables. It is not only used to measure the second-order information of the joint probability density function but can also measure higher-order moments.

Firstly, the correntropy between the random variables X and Y is defined as [23]

$$C(X,Y)=E\left[G_{\sigma }(X,Y)\right]={\displaystyle \int G_{\sigma }(X,Y)d{F}_{XY}(x,y)}$$

(14)

where $E\left[\cdot \right]$ represents the calculation of expectation, $G_{\sigma }(\cdot )$ is a shift-invariant Mercer kernel, $F_{XY}(x,y)$ stands for the joint probability distribution of variables $X$ and $Y$, and ${\displaystyle \int \cdot }$ represents the integral operation.

The Gaussian kernel is the most widely used kernel function in Mercer theory, and it is described as [29]

$$G_{\sigma }(X,Y)=G_{\sigma }(X-Y)=\exp \left(-{\displaystyle \frac{{‖X-Y‖}^2}{2{\sigma }^2}}\right)$$

(15)

where $\sigma >0$ represents the bandwidth of the Gaussian kernel.

However, the Gaussian kernel may not always be an optimal selection due to the appearance of singular matrices, making the filter break down. In this paper, we use the Cauchy kernel instead of the Gaussian kernel as the kernel function for the correntropy to better capture the heavy-tailed features in the noise. The Cauchy kernel has merits such as being less sensitive to kernel bandwidth and having better stability compared to the Gaussian kernel. It is defined as [27]

$$G_{\sigma }(X,Y)=G_{\sigma }(X-Y)={\displaystyle \frac{1}{1+{‖X-Y‖}^2/\sigma }}$$

(16)

where $\sigma >0$ is the bandwidth of the Cauchy kernel, which is used to adjust the degree of influence of the kernel function on non-Gaussian noise.

In real applications, the true joint distribution function of random variables is usually inaccessible; thus, the sample mean estimation of correntropy is often calculated to approximate the true joint distribution function [28]

$$\widehat{C}(X,Y)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{G}_{\sigma }(x_i,y_i)}$$

(17)

where $N$ is the sample number, and $x_i$ and $y_i$ are the sample data of the random variables $X$ and $Y$.

Correntropy can capture the higher-order moment information of the observed signal, which can be described by the Taylor series expansion of the kernel function [28]

$$\widehat{C}(X,Y)={\displaystyle \sum _{n=0}^{\infty }\frac{{(-1)}^n}{{\sigma }^n}}E\left[{(X-Y)}^{2n}\right]$$

(18)

It can be seen from (15) and (16) that when $X$ = $Y$, the correntropy of two random variables reaches its maximum, indicating the strongest similarity between them, which is the core principle of the maximum correntropy criterion. In addition, it can be seen from (18) that the correntropy is a weighted sum of the even-order moments of $X$ − $Y$. Unlike the minimum mean square error criterion in Kalman filtering, only the second-order moment is considered. Thus, the maximum correntropy criterion can be used to handle the influence of non-Gaussian noises.

Further, from (15) and (16), we can see that the Gaussian kernel gets exponentially closer to 0 as $X$ − $Y$ increases, which leads to an increase in the possibility of matrix singularity. As a comparison, the Cauchy kernel approaches zero much more slowly. Thus, it can effectively reduce the probability of matrix singularity. This is the reason why the Cauchy kernel instead of the Gaussian kernel is used as the kernel function for the correntropy in this paper.

3.1.2. Process of the Maximum Correntropy Filter Algorithm Based on Cauchy Kernel

According to the above-described maximum entropy criterion based on the Cauchy kernel in Section 3.1.1, we can further combine it with Kalman filtering to establish the maximum correntropy filter based on the Cauchy kernel.

As shown in (7) and (13), the system process and measurement models of gyros/star sensor integration can be summarized as the following form

$$\left\{\begin{array}{l}{\mathit{X}}_k={\mathit{F}}_{k/k-1}{\mathit{X}}_{k-1}+{\mathit{W}}_k\hfill \\ {\mathit{Z}}_k={\mathit{H}}_k{\mathit{X}}_k+{\mathit{V}}_k\hfill \end{array}\right.$$

(19)

where

$${\mathit{F}}_{k/k-1}=\left[\begin{array}{cc}{\mathit{I}}_{3\times 3}& -{\mathit{C}}_P^I\\ {\mathbf{0}}_{3\times 3}& {\mathit{I}}_{3\times 3}\end{array}\right]$$

(20)

The calculation procedure of traditional Kalman filter can be described as

Step 1: Initialization parameter.

$$\left\{\begin{array}{l}{\widehat{\mathit{X}}}_0=E[{\mathit{X}}_0]\\ {\mathit{P}}_0=E[({\mathit{X}}_0-{\widehat{\mathit{X}}}_0)({\mathit{X}}_0-{\widehat{\mathit{X}}}_0{)}^{\mathrm{T}}]\end{array}\right.$$

(21)

Step 2: Time update. Using the previous state estimate and its error covariance matrix in time k − 1, the state prediction and its error covariance matrix are calculated as [17]

$${\widehat{\mathit{X}}}_{k/k-1}={\mathit{F}}_{k/k-1}{\widehat{\mathit{X}}}_{k-1}$$

(22)

$${\mathit{P}}_{k/k-1}={\mathit{F}}_{k/k-1}{\mathit{P}}_{k-1}{\mathit{F}}_{k/k-1}^{\mathrm{T}}+{\mathit{Q}}_{k-1}$$

(23)

where ${\mathit{F}}_{k/k-1}$ is the state transition matrix described by (20); ${\widehat{\mathit{X}}}_{k/k-1}$ and ${\mathit{P}}_{k/k-1}$ are the state prediction and its error covariance at time k; ${\widehat{\mathit{X}}}_{k-1}$ and ${\mathit{P}}_{k-1}$ are the state estimation and its error covariance at the previous time k − 1; and ${\mathit{Q}}_{k-1}$ is the process noise covariance.

Step 3: Measurement update. The processes for measurement updates are as follows [17]

$${\mathit{K}}^k={\mathit{P}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}{({\mathit{H}}_k{\mathit{P}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}+{\mathit{R}}_k)}^{-1}$$

(24)

$${\widehat{\mathit{X}}}_k={\widehat{\mathit{X}}}_{k/k-1}+{\mathit{K}}_k\left[{\mathit{Z}}_k-{\mathit{H}}_k{\widehat{\mathit{X}}}_{k/k-1}\right]$$

(25)

$${\mathit{P}}_k=(\mathit{I}-{\mathit{K}}_k{\mathit{H}}_k){\mathit{P}}_{k/k-1}$$

(26)

where ${\mathit{R}}_k$ is the measurement noise covariance; ${\mathit{K}}_k$ is the filtering gain; and ${\mathit{H}}_k$ is the measurement matrix.

However, due to the weak light environment in space, the star sensor’s measurement information is subject to be disturbed, leading to the occurrence of non-Gaussian measurement noise. This further reduces the estimation accuracy of the Kalman filtering and may even cause divergence, resulting in the loss of its ability to correct the gyros.

According to the maximum correntropy criterion based on the Cauchy kernel, we construct a new filtering method to modify the measurement update process of Kalman filtering to handle the influence of non-Gaussian measurement noise. Firstly, based on the system model shown in (19), we have

$$\left[\begin{array}{c}{\widehat{\mathit{X}}}_{k/k-1}\\ {\mathit{Z}}_k\end{array}\right]=\left[\begin{array}{c}\mathit{I}\\ {\mathit{H}}_k\end{array}\right]{\mathit{X}}_k+{\mathit{A}}_k$$

(27)

where

$${\mathit{A}}_k=\left[\begin{array}{c}{\widehat{\mathit{X}}}_{k/k-1}-{\mathit{X}}_k\\ {\mathit{V}}_k\end{array}\right],\mathrm{and}\mathit{E}\left[{\mathit{A}}_k{\mathit{A}}_k^T\right]=\left[\begin{array}{cc}{\mathit{P}}_{k/k-1}& 0\\ 0& {\mathit{R}}_k\end{array}\right]={\mathit{B}}_k{\mathit{B}}_k^{\mathrm{T}}$$

(28)

where ${\mathit{B}}_k$ is the Cholesky decomposition of $\mathit{E}\left[{\mathit{A}}_k{\mathit{A}}_k^T\right]$.

Multiplying both sides of (27) by ${\mathit{B}}_k^{-1}$ yields

$${\mathit{Y}}_k={\mathit{M}}_k{\mathit{X}}_k+{\mathit{E}}_k$$

(29)

where ${\mathit{Y}}_k={\mathit{B}}_k^{-1}\left[\begin{array}{c}{\widehat{\mathit{X}}}_{k/k-1}\\ {\mathit{Z}}_k\end{array}\right]$, ${\mathit{M}}_k={\mathit{B}}_k^{-1}\left[\begin{array}{c}\mathit{I}\\ {\mathit{H}}_k\end{array}\right]$, and ${\mathit{E}}_k={\mathit{B}}_k^{-1}{\mathit{A}}_k$.

Then, a correntropy cost function can be established as

$$\mathit{J}={\displaystyle \frac{1}{L}}{\displaystyle \sum _{i=1}^L{G}_{\sigma }\left(E_k^i\right)}={\displaystyle \frac{1}{L}}{\displaystyle \sum _{i=1}^L{G}_{\sigma }\left(y_k^i-{\mathit{m}}_k^i{\mathit{X}}_k\right)}$$

(30)

where L = n + m represents the dimension of ${\mathit{Y}}_k$, and m and n are the dimensions of ${\mathit{Z}}_k$ and ${\mathit{X}}_k$, respectively; $E_k^i$ and $y_k^i$ represent the ith elements of ${\mathit{E}}_k$ and ${\mathit{Y}}_k$; and ${\mathit{m}}_k^i$ is the ith row of ${\mathit{M}}_k$.

The optimal state estimation can be obtained using the following equation.

$${\displaystyle \frac{\mathsf{\partial }\mathit{J}}{\mathsf{\partial }{\mathit{X}}_k}}={\displaystyle \frac{2}{\sigma }}{\displaystyle \frac{1}{L}}{\displaystyle \sum _{i=1}^L{G}_{\sigma }{(E_k^i)}^2{\mathit{m}}_k^{i,T}\left(y_k^i-{\mathit{m}}_k^i{\mathit{X}}_k\right)}=0$$

(31)

After simplification, it is evident that the following result can be obtained [36]

$$\begin{array}{cc}{\mathit{X}}_k& ={\left({\displaystyle \sum _{i=1}^L{G}_{\sigma }{(E_k^i)}^2{\mathit{m}}_k^{i,T}{\mathit{m}}_k^i}\right)}^{-1}\left({\displaystyle \sum _{i=1}^L{G}_{\sigma }{(E_k^i)}^2{\mathit{m}}_k^{i,T}{y}_k^i}\right)\hfill \\ & ={\left({\mathit{M}}_k^T{\mathit{D}}_k{\mathit{M}}_k\right)}^{-1}\left({\mathit{M}}_k^T{\mathit{D}}_k{\mathit{Y}}_k\right)\hfill \end{array}$$

(32)

where ${\mathit{D}}_k=\left[\begin{array}{cc}{\mathit{D}}_k^x& \mathbf{0}\\ \mathbf{0}& {\mathit{D}}_k^z\end{array}\right]$ with

$$\begin{array}{l}{\mathit{D}}_k^x=diag\left(\begin{array}{ccc}{G}_{\sigma }{(E_k^1)}^2& ,\dots ,& G_{\sigma }{(E_k^n)}^2\end{array}\right)\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}=diag\left(\begin{array}{ccc}{\left({\displaystyle \frac{1}{1+{\left({\widehat{X}}_{k/k-1}^1-{\mathit{m}}_k^1{\widehat{\mathit{X}}}_k\right)}^2/\sigma }}\right)}^2& ,\dots ,& {\left({\displaystyle \frac{1}{1+{\left({\widehat{X}}_{k/k-1}^n-{\mathit{m}}_k^n{\widehat{\mathit{X}}}_k\right)}^2/\sigma }}\right)}^2\end{array}\right)\end{array}$$

(33)

$$\begin{array}{l}{\mathit{D}}_k^z=diag\left(\begin{array}{ccc}{G}_{\sigma }{(E_k^{n+1})}^2& ,\dots ,& G_{\sigma }{(E_k^{n+m})}^2\end{array}\right)\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}=diag\left(\begin{array}{ccc}{\left({\displaystyle \frac{1}{1+{\left(Z_k^1-{\mathit{m}}_k^{n+1}{\widehat{\mathit{X}}}_k\right)}^2/\sigma }}\right)}^2& ,\dots ,& {\left({\displaystyle \frac{1}{1+{\left(Z_k^m-{\mathit{m}}_k^{n+m}{\widehat{\mathit{X}}}_k\right)}^2/\sigma }}\right)}^2\end{array}\right)\end{array}$$

(34)

where ${\widehat{X}}_{k/k-1}^i(i=1,2,\cdots ,n)$ and $Z_k^i(i=1,2,\cdots ,m)$ are the ith elements of ${\widehat{\mathit{X}}}_{k/k-1}$ and ${\mathit{Z}}_k$.

By expanding (32), we can obtain

$$\begin{array}{l}{\mathit{X}}_k ={\left({\mathit{M}}_k^{\mathrm{T}}{\mathit{D}}_k{\mathit{M}}_k\right)}^{-1}\left({\mathit{M}}_k^{\mathrm{T}}{\mathit{D}}_k{\mathit{Y}}_k\right)\\ \hspace{1em} ={\left({\left({\mathit{B}}_k^{P,\mathrm{T}}\right)}^{-1}{\mathit{D}}_k^x{\left({\mathit{B}}_k^P\right)}^{-1}+{\mathit{H}}_k^{\mathrm{T}}{\left({\mathit{B}}_k^{R,\mathrm{T}}\right)}^{-1}{\mathit{D}}_k^z{\left({\mathit{B}}_k^R\right)}^{-1}{\mathit{H}}_k\right)}^{-1}\\ \hspace{1em} \left({\left({\mathit{B}}_k^{P,\mathrm{T}}\right)}^{-1}{\mathit{D}}_k^x{\left({\mathit{B}}_k^P\right)}^{-1}{\widehat{\mathit{X}}}_{k/k-1}+{\mathit{H}}_k^{\mathrm{T}}{\left({\mathit{B}}_k^{R,\mathrm{T}}\right)}^{-1}{\mathit{D}}_k^z{\left({\mathit{B}}_k^R\right)}^{-1}{\mathit{Z}}_k\right)\end{array}$$

(35)

where ${\mathit{B}}_k^P$ and ${\mathit{B}}_k^R$ are the Cholesky decomposition of ${\mathit{P}}_{k/k-1}$ and ${\mathit{R}}_k$.

Applying the matrix inversion lemma (${(A+B{D}^{-1}C)}^{-1}=A^{-1}-A^{-1}B{(D+C{A}^{-1}B)}^{-1}C{A}^{-1}$), (35) can be become

$${\widehat{\mathit{X}}}_k={\widehat{\mathit{X}}}_{k/k-1}+{\overline{\mathit{K}}}_k\left[{\mathit{Z}}_k-{\mathit{H}}_k{\widehat{\mathit{X}}}_{k/k-1}\right]$$

(36)

where

$${\overline{\mathit{K}}}_k={\overline{\mathit{P}}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}{({\mathit{H}}_k{\overline{\mathit{P}}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}+{\overline{\mathit{R}}}_k)}^{-1}$$

(37)

where ${\overline{\mathit{P}}}_{k/k-1}={\mathit{B}}_k^P{\left({\mathit{D}}_k^x\right)}^{-1}{\mathit{B}}_k^{P,\mathrm{T}}$, ${\overline{\mathit{R}}}_k={\mathit{B}}_k^R{\left({\mathit{D}}_k^z\right)}^{-1}{\mathit{B}}_k^{R,\mathrm{T}}$. The specific derivation from (35)–(37) can be found in [36].

To estimate the system state through (36), it is necessary to compute the filter gain ${\tilde{\mathit{K}}}_k$ in (35) firstly. The calculation of the gain inherently relies on the posterior state estimate ${\widehat{\mathit{X}}}_k$, which is involved in (33) and (34). This clearly indicates an iterative process for state estimation. Therefore, we outline the measurement update process of the maximum correntropy filter based on the Cauchy kernel as follows:

(i) Let the iteration index t = 1 and ${\widehat{\mathit{X}}}_k^1={\widehat{\mathit{X}}}_{k/k-1}$, ${\widehat{\mathit{X}}}_k^{t+1}$ can be achieved by the following process

$${\widehat{\mathit{X}}}_k^{t+1}={\widehat{\mathit{X}}}_k^t+{\tilde{\mathit{K}}}_k\left[{\mathit{Z}}_k-{\mathit{H}}_k{\widehat{\mathit{X}}}_k^t\right]$$

(38)

$${\tilde{\mathit{K}}}_k={\tilde{\mathit{P}}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}{({\mathit{H}}_k{\tilde{\mathit{P}}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}+{\tilde{\mathit{R}}}_k)}^{-1}$$

(39)

$${\tilde{\mathit{P}}}_{k/k-1}={\mathit{B}}_k^P{\left({\tilde{\mathit{D}}}_k^x\right)}^{-1}{\mathit{B}}_k^{P,T}$$

(40)

$${\tilde{\mathit{R}}}_k={\mathit{B}}_k^R{\left({\tilde{\mathit{D}}}_k^z\right)}^{-1}{\mathit{B}}_k^{R,T}$$

(41)

where ${\tilde{\mathit{D}}}_k^x=diag\left(\begin{array}{c}{G}_{\sigma }({\tilde{E}}_k^1{)}^2,\dots ,G_{\sigma }({\tilde{E}}_k^n{)}^2\end{array}\right)$ and ${\tilde{\mathit{D}}}_k^z=diag\left(\begin{array}{c}{G}_{\sigma }({\tilde{E}}_k^{n+1}{)}^2,\dots ,G_{\sigma }({\tilde{E}}_k^{n+m}{)}^2\end{array}\right)$; ${\tilde{E}}_k^i$ is the element of ${\tilde{\mathit{E}}}_k$ and ${\tilde{\mathit{E}}}_k={\mathit{Y}}_k-{\mathit{M}}_k{\widehat{\mathit{X}}}_k^t$.

(ii) Compare the posterior estimation from the current iteration ${\widehat{\mathit{X}}}_k^{t+1}$ with those from the previous iteration ${\widehat{\mathit{X}}}_k^t$. If it satisfies the following condition, the iterative process is terminated; otherwise, set the final state estimation ${\widehat{\mathit{X}}}_k={\widehat{\mathit{X}}}_k^{t+1}$ and return to Step (i).

$${\displaystyle \frac{‖{\widehat{\mathit{X}}}_k^{t+1}-{\widehat{\mathit{X}}}_k^t‖}{‖{\widehat{\mathit{X}}}_k^t‖}}\le \epsilon$$

(42)

where $\epsilon$ is a relatively small number.

(iii) Compute the posterior error covariance matrix as

$${\mathit{P}}_k=(\mathit{I}-{\tilde{\mathit{K}}}_k{\mathit{H}}_k){\mathit{P}}_{k/k-1}$$

(43)

At this time, the measurement update process at time k of the maximum correntropy filter based on the Cauchy kernel can be terminated.

It can be seen from (33) to (37) that as the bandwidth $\sigma \to \infty$, ${\mathit{D}}_k^x$ and ${\mathit{D}}_k^z$ will tend to unity matrices. In this case, the maximum correntropy filter based on the Cauchy kernel will be degraded to the Kalman filter and its robustness will gradually vanish. If $\sigma \to 0$, ${\mathit{D}}_k^x$ and ${\mathit{D}}_k^z$ will become 0, the calculation of filter gain ${\tilde{\mathit{K}}}_k$ will be out of operation and the filter will be disabled. Thus, an appropriate kernel bandwidth is vital to affect the performance of the maximum correntropy filter based on the Cauchy kernel. Currently, it is usually determined by prior information or a multiple test computation prior to the filtering solution.

3.2. Event-Driven Conditions

From the above filtering algorithm process, it can be seen that multiple iterations are involved. This means it needs a large computational cost for the filtering calculation at each time. Thus, in the times without non-Gaussian measurement noise, the maximum correntropy filter based on the Cauchy kernel will involve too many redundant calculations. Due to this reason, this section establishes an adaptive event-driven mechanism to effectively recognize the non-Gaussian measurement noise, which can significantly improve the computational efficiency of the maximum correntropy filter based on the Cauchy kernel.

The constructed event-driven conditions are related to the filter’s innovation vector information, which is calculated as [37]

$${\tilde{\mathit{Z}}}_k={\mathit{Z}}_k-{\mathit{H}}_k{\widehat{\mathit{X}}}_{k/k-1}$$

(44)

In the normal case, the innovation ${\tilde{\mathit{Z}}}_k$ should adhere to a Gaussian distribution and its covariance is [37]

$${\mathit{P}}_{ZZ,k/k-1}={\mathit{H}}_k{\mathit{P}}_{k/k-1}{\mathit{H}}_k^{\mathrm{T}}+{\mathit{R}}_k$$

(45)

Consider that ${\mathit{P}}_{ZZ,k/k-1}$ is a positive matrix, we can construct a unitary matrix ${\mathit{L}}_k$ to make the following relationship hold

$${\mathit{L}}_k^{\mathrm{T}}{\mathit{P}}_{ZZ,k/k-1}{\mathit{L}}_k={\mathrm{\Sigma }}_k$$

(46)

where ${\mathrm{\Sigma }}_k=diag\left\{\begin{array}{c}{\lambda }_{k1},\hspace{0.33em}{\lambda }_{k2},\dots ,{\lambda }_{k{n}_i}\end{array}\right\}\in n_i\times n_i$, $\begin{array}{c}{\lambda }_{k1},\hspace{0.33em}{\lambda }_{k2},\dots ,{\lambda }_{k{n}_i}\end{array}$ are the eigenvalues of ${\mathit{P}}_{ZZ,k/k-1}$ and $n_i$ represents the number of eigenvalues.

Further, define ${\mathit{F}}_k={\mathit{L}}_k{\mathrm{\Sigma }}_k^{-1/2}$, the innovation vector can be normalized as

$${\overline{\mathit{Z}}}_k={\mathit{F}}_k^{\mathrm{T}}{\tilde{\mathit{Z}}}_k$$

(47)

Equation (47) is called Mahalanobis transformation. The elements of ${\overline{\mathit{Z}}}_k$ should be uncorrelated and obey the standard Gaussian distribution [22,37]. Thus, the event-driven conditions for the detection of the non-Gaussian measurement noise can be constructed as

$${\tau }_k=\left\{\begin{array}{l}1\hspace{1em}\hspace{1em}\hspace{1em}{\kappa }_{\beta }\le {‖{\overline{\mathit{Z}}}_k‖}_{\infty }\le {\kappa }_{\alpha }\\ 0\hspace{1em}\hspace{1em}\hspace{1em}Otherwise\end{array}\right.$$

(48)

where ${‖{\overline{\mathit{Z}}}_k‖}_{\infty }$ represents the infinity norm of ${\overline{\mathit{Z}}}_k$, i.e., ${‖{\overline{\mathit{Z}}}_k‖}_{\infty }=\max \left\{{\overline{\mathit{Z}}}_{k,1},\hspace{0.33em}{\overline{\mathit{Z}}}_{k,2},\dots ,{\overline{\mathit{Z}}}_{k,q}\right\}$; ${\kappa }_{\beta }$ and ${\kappa }_{\alpha }$ are the lower threshold and the higher threshold, respectively.

Based on (48), there are three cases that will appear:

Case 1: ${‖{\overline{\mathit{Z}}}_k‖}_{\infty }<{\kappa }_{\beta }$ and ${\tau }_k=0$. This case considers that the integration of the gyros/star sensor can be predictable and the measurement does not need to be applied to update the estimation. The purpose of this case is to reduce the communication burden between the gyros and star sensor so as to decrease the whole system’s computational cost.

Case 2: ${\kappa }_{\beta }\le {‖{\overline{\mathit{Z}}}_k‖}_{\infty }\le {\kappa }_{\alpha }$ and ${\tau }_k=1$. This case considers that the measurement is normal and the non-Gaussian noise is not involved;

Case 3: ${‖{\overline{\mathit{Z}}}_k‖}_{\infty }>{\kappa }_{\alpha }$ and ${\tau }_k=0$. This case considers that the measurement is deviant and the non-Gaussian noise is involved.

It should be noted that the selection of thresholds ${\kappa }_{\beta }$ and ${\kappa }_{\alpha }$ is very important for the detection of the non-Gaussian measurement noise. ${\kappa }_{\beta }$ is chosen to balance the communication burden and estimation performance, which can be set according to the required communication rate. ${\kappa }_{\alpha }$ is always set according to the confidence level. Through the simulation test, a 97.5% confidence level can obtain a relatively satisfactory detection performance.

3.3. Algorithm Description

Based on the maximum correntropy Kalman filter based on the Cauchy kernel described in Section 3.1 and event-driven conditions constructed in Section 3.2, we can establish the algorithm of the event-driven maximum correntropy Kalman filter based on the Cauchy kernel, whose computational process is illustrated in Algorithm 1. It can balance the computational efficiency and robustness of the filter, reducing unnecessary computational load while ensuring the algorithm’s anti-interference performance.

Algorithm 1. Event-driven- maximum correntropy filter based on Cauchy kernel (ED-MCFCK)
1	Initialization: ${\widehat{\mathit{x}}}_0$, ${\mathit{P}}_0$
2	For k = 1,2,…
	{
3	Compute ${\widehat{\mathit{X}}}_{k/k-1}$ and ${\mathit{P}}_{k/k-1}$ as (22) and (23) of Kalman filtering procedures;
4	Calculate filter’s innovation vector as well as its covariance as (44) and (45), and further construct the normalized innovation vector as (47).
5	Compute and conduct the event driven condition by (48).
6	If ${‖{\overline{\mathit{Z}}}_k‖}_{\infty }<{\kappa }_{\beta }$ and ${\tau }_k=0$, set the state prediction ${\widehat{\mathit{X}}}_{k/k-1}$ and ${\mathit{P}}_{k/k-1}$ as the final state estimation.
7	If ${\kappa }_{\beta }\le {‖{\overline{\mathit{Z}}}_k‖}_{\infty }\le {\kappa }_{\alpha }$ and ${\tau }_k=1$, execute (24)–(26) to obtain system state esti-mation ${\widehat{\mathit{X}}}_k$ and ${\mathit{P}}_k$.
8	If ${‖{\overline{\mathit{Z}}}_k‖}_{\infty }>{\kappa }_{\alpha }$ and ${\tau }_k=0$, (indicating the presence of non-Gaussian noise in the measurement information), the MCFCK is driven.
	{
9	Let the iteration index t = 1 and ${\widehat{\mathit{X}}}_k^1={\widehat{\mathit{X}}}_{k/k-1}$;
10	Compute ${\tilde{\mathit{P}}}_{k/k-1}$ and ${\tilde{\mathit{R}}}_k$ according to (40) and (41), and further calcu late filter gain ${\tilde{\mathit{K}}}_k$ by (39).
11	Compute ${\widehat{\mathit{X}}}_k^t$ by (38) and test the Condition (42);
12	If Condition $\ge \epsilon$
13	Let t = t + 1 (iterations) and ${\widehat{\mathit{X}}}_k={\widehat{\mathit{X}}}_k^{t+1}$, and conduct the next iteration.
14	Else, the iterative process is terminated and ${\widehat{\mathit{X}}}_k={\widehat{\mathit{X}}}_k^t$.
	}
15	Compute the posterior error covariance matrix by (43).
	}

4. Simulation Evaluations and Discussions

Simulation evaluations have been conducted to comprehensively validate the efficacy of the proposed event-driven maximum correntropy filter based on the Cauchy kernel (ED-MCFCK) for a two-dimensional turntable’s spatial orientation using the gyros/star sensor integration. A comparison analysis of the proposed ED-MCFCK with the Kalman filter and MCF based on the Gaussian kernel (abbreviated as MCF) as well as MCFCK is also conducted in this section for the handling of non-Gaussian noise involved in the measurements.

As shown in Figure 3, a trajectory of spacecraft is designed according to its actual flight process, which is a classical parabolic curve motion. Further, the attitude of the spacecraft at every time step, including yaw, pitch, and roll is depicted in Figure 4. The initial attitude of the spacecraft was (142.16°, −89.78°, 63.16°) and the total simulation time was 3600 s. Accordingly, a two-dimensional turntable was installed on the spacecraft, and the tracking vector between the spacecraft and the target star is also shown in Figure 3.

For the filtering loop shown in Figure 2, the gyros’ sampling rates are both 20 Hz, the star sensor’s sampling rate is 1 Hz, and the other simulation parameters are listed in

Table 1. Simulation parameters of gyros/star sensor integration.

Information Sources	Simulation Parameters	Values
Gyros	Gyros’ constant drift	0.1°/h
	Gyros’ random walk coefficient	0.01°/$\sqrt{h}$
	Scale factor error	100 PPM
	Misalignment error	5″
Star sensor	Goniometric error	5″ (3σ)

. The initial state and its error covariance, as well as the initial process noise covariance, are set as

$$\left\{\begin{array}{l}{\mathit{X}}_0={[{\mathbf{0}}_{6\times 1}]}^{\mathrm{T}}\hfill \\ {\mathit{P}}_0=diag[\begin{array}{ccc}{(1rad)}^2& {(1rad)}^2& \begin{array}{cccc}{(1rad)}^2& {(1″)}^2& {(1″)}^2& {(1″)}^2\end{array}\end{array}]\hfill \\ {\mathit{Q}}_0=diag[\begin{array}{ccc}{({0.01}^{\circ }/\sqrt{h})}^2& {({0.01}^{\circ }/\sqrt{h})}^2& {({0.01}^{\circ }/\sqrt{h})}^2\end{array}]\hfill \end{array}\right.$$

(49)

A.

Accuracy Evaluation and Analysis

To verify the accuracy and robustness of the proposed ED-MCFCK in the presence of non-Gaussian measurement noise, two typical non-Gaussian measurement noise scenarios were considered, i.e., non-Gaussian noise with outliers in measurements and $\alpha$ stable noise in measurements.

By using the trial and error method [38], the bandwidth parameters $\sigma$ in the above-mentioned maximum correntropy filters were determined to be 13. For the designed event-driven conditions, the confidence level is set as 97.5% (α = 0.025) with 3 degrees of freedom (m = 3).

(1)

Non-Gaussian noise with outliers in measurements

To simulate this situation, a random noise with Gaussian distribution was directly injected into the measurement during the time interval (1500 s, 2500 s) with a probability of $\rho =0.5$. Thus, the measurement becomes the following

$${\mathit{Z}}_k=\left\{\begin{array}{l}{\mathit{Z}}_k,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} others\hfill \\ (1-\rho ){\mathit{Z}}_k+\rho [{\mathit{Z}}_k+N(0,{(4\times {10}^{-4}rad)}^2)],(1500s,2500s)\hfill \end{array}\right.$$

(50)

It can be seen that the measurement noise in (50) is mixed Gaussian, which is a non-Gaussian and strong impulsive noise (outliers).

For the case of non-Gaussian noise with outliers in measurements, Figure 5 and Figure 6 illustrate the orientation errors in azimuth and pitch achieved by the above-mentioned Kalman filter, MCF, MCFCK and ED-MCFCK, respectively. Further, the corresponding 3$\sigma$ error boundary for the ED-MCFCK, which is calculated based on the 10 runs of Monte Carlo simulations as in [39], is also provided in Figure 5 and Figure 6. As shown in Figure 5 and Figure 6, during the time interval (1500 s, 2500 s) with non-Gaussian noise, the Kalman filter has no ability and fails to handle the abnormal non-Gaussian noise. This further deteriorates the obtained state estimate during the filtering measurement update process, resulting in a significant degradation in orientation accuracy. The orientation errors in azimuth and pitch exhibit severe fluctuations. In contrast, MCF, MCFCK and ED-MCFCK can leverage the robustness of the Gaussian/Cauchy kernel to effectively handle the non-Gaussian noise information, leading to stable orientation errors in azimuth and pitch.

The global orientation accuracy and the corresponding Root Mean Square Error (RMSE) are further used to evaluate the orientation accuracy, where the global orientation error is defined as

$$‖\delta \theta ‖=\sqrt{\delta {\theta }_A^2+\delta {\theta }_E^2}$$

(51)

where $\delta {\theta }_A$ and $\delta {\theta }_E$ are the azimuth error and pitch error of the turntable, respectively.

The RMSE at time k is defined as

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{\left|\delta {\theta }_k\right|}^2}}$$

(52)

where N is the total simulation time, and $\delta {\theta }_k$ represents the orientation error at time k.

Table 2. Mean RMSEs of the orientation errors by Kalman filter, MCF, MCFCK and ED-MCFCK for the case of non-Gaussian noise with outliers.

Filtering Methods	Global Orientation Error (″)	Azimuth Error (″)	Pitch Error (″)
Kalman Filter	24.60	16.22	18.49
MCF	10.87	5.23	9.53
MCFCK	10.76	4.81	9.62
ED-MCFCK	10.6	4.79	9.46

shows the mean RMSEs of global orientation error as well as the azimuth and pitch errors for each filter. Due to the ability of robustness on non-Gaussian noise, the mean RMSEs of global orientation error as well as the azimuth and pitch errors achieved by the MCF are 55.8%, 67.76% and 48.46% lower than those of the Kalman filter. Meanwhile, the mean RMSEs of global orientation error as well as the azimuth and pitch errors achieved by the MCFCK are 56.26%, 70.35% and 47.97% lower than those of the Kalman filter. Further, the mean RMSEs of global orientation error as well as the azimuth and pitch errors achieved by the ED-MCFCK are 56.91%, 70.47% and 48.84% lower than those of the Kalman filter. These indicate that the MCF, MCFCK and ED-MCFCK exhibit a similar estimation performance for the case of non-Gaussian noise with outliers in measurements and they have a better sensitivity to the non-Gaussian measurement noises compared to the Kalman filter.

(2)

$\alpha$ stable noise in measurements

The $\alpha$ stable noise is a significant type of non-Gaussian random distribution, and its two most important features are the stability of the probability distribution and the heavy-tailed probability density function. The characteristic function of the $\alpha$ stable noise is given by [40]

$$\phi (t)=\exp \{j\delta t-\alpha {\left|t\right|}^{\beta }[1+j\zeta \mathrm{sgn}(t)\omega (t,\beta )]\}$$

(53)

with

$$\omega (t,\beta )=\left\{\begin{array}{c}\tan ({\displaystyle \frac{\beta \pi }{2}}),\beta \ne 1\\ {\displaystyle \frac{2}{\pi }}\log \left|t\right|,\beta =1\end{array}\right.$$

(54)

$$\mathrm{sgn}(t)=\left\{\begin{array}{l}1,\hspace{1em} \mathrm{t}\text{>}0\hfill \\ 0,\hspace{1em} \mathrm{t}=0\hfill \\ -1, \mathrm{t}\text{<}0\hfill \end{array}\right.$$

(55)

where $\beta \in (0,2]$ is the characteristic index, $\zeta \in [-1,\hspace{0.33em}1]$ is the symmetry parameter, $\delta$ is the location factor and $\alpha >0$ is the dispersion factor.

To test the performance of the four filters under $\alpha$ stable noise, the above noise form was injected into the measurement during the time intervals (1000 s, 1500 s) and (2500 s, 3000 s), where $(\beta ,\hspace{0.33em}\zeta ,\hspace{0.33em}\alpha ,\hspace{0.33em}\delta )=(1.8,\hspace{0.33em}0,\hspace{0.33em}1.6,\hspace{0.33em}0)$ is set in this simulation.

Figure 7 and Figure 8 provide the orientation errors in azimuth and pitch achieved by the Kalman filter, MCF, MCFCK and ED-MCFCK as well as the 3$\sigma$ error boundary of ED-MCFCK for the case of $\alpha$ stable noise in measurement. It can be seen that as stated previously, MCF is unable to handle heavy-tailed noise anomalies due to the limitation of the Gaussian kernel function; thus, its estimation error is divergent. In contrast, MCFCK and ED-MCFCK can effectively handle the influence of $\alpha$ stable noise with a heavy-tailed characteristic through the robustness of the Cauchy kernel. They show similar performance with the case of non-Gaussian noise with outliers in measurements and have a stronger robustness than the Kalman filter.

Table 3. Mean RMSEs of the orientation errors by Kalman filter, MCF, MCFCK and ED-MCFCK for the case of $\sigma$ stable noise.

Filtering Methods	Global Orientation Error (″)	Azimuth Error (″)	Pitch Error (″)
Kalman Filter	26.27	18.08	19.06
MCF	185.93	120.61	141.5
MCFCK	9.67	3.86	8.87
ED-MCFCK	8.82	3.26	8.20

shows the mean RMSEs of global orientation error as well as the azimuth and pitch errors by the above four filters. Due to the robustness of the Cauchy kernel of MCFCK and ED-MCFCK, the mean RMSEs of the global orientation error as well as the azimuth and pitch errors achieved by the MCFCK are 63.19%, 78.65% and 53.46% lower than those of the Kalman filter, and the mean RMSEs of ED-MCFCK are 66.43%, 81.97% and 56.98% lower than those of the Kalman filter. However, due to the limitation of the Gaussian kernel in handling the heavy-tailed noise anomalies, the mean RMSEs of global orientation error as well as the azimuth and pitch errors by the MCF increased by 2008.05%, 3599.69% and 1625.61% compared the ED-MCFCK. Therefore, it is evident that the proposed ED-MCFCK has a superior robust performance to curb heavy-tailed non-Gaussian noise such as $\alpha$ stable noise, leading to improved orientation accuracy in the complex space environment.

B.

Computational Performance Evaluation

The computational efficiency of the proposed ED-MCFCK is also discussed by comparison with the Kalman filter, MCF and MCFCK via a Monte Carlo simulation. It was conducted 10 times on an Intel(R) Core(TM) i5-9300H 2.40 GHz PC with 8 GB memory to study the computational performance of the proposed ED-MCFCK based on the above two types of non-Gaussian measurement noise scenarios. The computational times for each Monte Carlo run of the Kalman filter, MCF, MCFCK and ED-MCFCK are listed in

Table 4. The computational times for each Monte Carlo run of the Kalman Filter, MCF, MCFCK and ED-MCFCK based on the simulation cases.

Filtering Methods	Non-Gaussian Noise with Outliers	$\mathit{\alpha }$ Stable Noise
	Computational Time	Computational Time
Kalman Filter	2.78 s	2.80 s
MCF	5.22 s	6.98 s
MCFCK	5.14 s	5.62 s
ED-MCFCK	2.84 s	2.89 s

. Further, the relative efficiencies of all the filters compared to the ED-MCFCK are also computed and shown in Figure 9.

It can be seen from Table 4 and Figure 9 that the Kalman filter has optimal real-time performance and its computational time for each Monte Carlo run is smallest among all the four filters. This is the reason why it is widely used in engineering practices. Further, the computational times of MCF and MCFCK are significantly longer than those of the Kalman filter. The filtering computational time of MCF is 87.8% (for the first case) and 149.3% (for the second case) longer than that of the Kalman filter. The filtering computational time of MCFCK is 84.9% (for the first case) and 100.7% (for the second case) longer than that of the Kalman filter. They have a similar computational performance since multiple iterations are involved in filtering every time and result in unnecessary redundant calculations. Both MCF and MCFCK’s real-time performance are poor; thus, they are not suitable for the practice applications. As a comparison, the proposed ED-MCFCK establishes a flexible event-driven mechanism to detect and process the case of non-Gaussian measurement noise. This effectively avoids the unnecessary redundant calculations and significantly improves the computational efficiencies of MCF and MCFCK. The ED-MCFCK has an approximated computational efficiency compared to the Kalman filter (97.9% for the first case and 96.9% for the second case). Thus, through the designed event-driven mechanism, ED-MCFCK can simultaneously obtain a good robustness on non-Gaussian measurement noise and a strong real-time performance in practical engineering.

The results of the above simulation evaluations and analysis demonstrate that the proposed ED-MCFCK not only has a superior ability to deal with the influence of non-Gaussian measurement noise but can also achieve superior real-time spatial applications by avoiding unnecessary computational costs in time points without non-Gaussian measurement noise, leading to improved spatial orientation performance for the gyros/star sensor integration used in the two-dimensional turntable structure.

5. Conclusions

This paper proposes an event-driven maximum correntropy filter based on the Cauchy kernel to improve the performance of gyros/star sensor integration used for spatial orientation under the environment of non-Gaussian measurement noise. The main contributions of this paper can be concluded as follows: (i) it establishes and derives a direct installation mode of gyros/star sensor integration as well as its mathematical model to ensure the turntable control stability; (ii) it further proposes an event-driven maximum correntropy filter based on the Cauchy kernel to handle the influence of non-Gaussian measurement noise on the above gyros/star sensor integration, in which an event-driven mechanism is also constructed based on the innovation information to reduce the unnecessary computational costs involved in the filtering solution. The results of the simulations and comparison analysis indicate that the proposed methodology can effectively curb the non-Gaussian measurement noise through the maximum correntropy criterion based on the Cauchy kernel, and it also can avoid the unnecessary computational cost of MC filters such as MCF and MCFCK, leading to improved robustness for the gyros/star sensor integration of spatial orientation in comparison with the conventional Kalman filter and MC filters.

Future research will focus on the combination of the proposed ED-MCFCK with advanced intelligent methods such as neural networks or deep learning to improve its robustness by adaptively adjusting the kernel bandwidth in real-time to accommodate the complex spatial environment. Further, it will also concentrate on the practical/extension application of the ED-MCFCK in the fields of vehicle navigation or target tracking.