A novel Gaussian sum quaternion constrained cubature Kalman filter algorithm for GNSS/SINS integrated attitude determination and positioning system

The quaternion cubature Kalman filter (QCKF) algorithm has emerged as a prominent nonlinear filter algorithm and has found extensive applications in the field of GNSS/SINS integrated attitude determination and positioning system (GNSS/SINS-IADPS) data processing for unmanned aerial vehicles (UAV). However, on one hand, the QCKF algorithm is predicated on the assumption that the random model of filter algorithm, which follows a white Gaussian noise distribution. The noise in actual GNSS/SINS-IADPS is not the white Gaussian noise but rather a ubiquitous non-Gaussian noise. On the other hand, the use of quaternions as state variables is bound by normalization constraints. When applied directly in nonlinear non-Gaussian system without considering normalization constraints, the QCKF algorithm may result in a mismatch phenomenon in the filtering random model, potentially resulting in a decline in estimation accuracy. To address this issue, we propose a novel Gaussian sum quaternion constrained cubature Kalman filter (GSQCCKF) algorithm. This algorithm refines the random model of the QCKF by approximating non-Gaussian noise with a Gaussian mixture model. Meanwhile, to account for quaternion normalization in attitude determination, a two-step projection method is employed to constrain the quaternion, which consequently enhances the filtering estimation accuracy. Simulation and experimental analyses demonstrate that the proposed GSQCCKF algorithm significantly improves accuracy and adaptability in GNSS/SINS-IADPS data processing under non-Gaussian noise conditions for Unmanned Aerial Vehicles (UAVs).

1 Introduction

Currently, Global Navigation Satellite Systems (GNSS) and Strapdown Inertial Navigation Systems (SINS) have experienced rapid development in both military and civilian fields. However, a single navigation technology, although it has its advantages in a specific environment, generally performs poorly in conventional occasions and may even fail to complete the task (Groves, 2008; Noureldin et al., 2013). Therefore, combining GNSS and SINS to form a GNSS/SINS integrated positioning and attitude determination system (GNSS/SINS-IADPS) can maximize their advantages and compensate for each other’s limitations (Teunissen and Montenbruck, 2017; Farrell and de Haag, 2020). At present, GNSS/SINS-IADPS has become one of the key technologies in aviation, aerospace and land navigation systems, and receives more attention in military strike, civil aviation, economic construction, and scientific research in various countries (Li and Chen, 2022; Boguspayev et al., 2023). The foundation of GNSS/SINS-IADPS relies on attitude representation methods, which include the attitude quaternion method, Rodrigues parameters and so on (Savage, 1998; Markley, 2003). The attitude quaternion method, in particular, has increasingly garnered attention in the realm of GNSS/SINS-IADPS data processing due to its advantages of avoiding singularities, requiring less computational effort, offering higher accuracy, and enabling complete attitude maneuvers. Consequently, filter algorithms predicated on the attitude quaternion have become a pivotal technology for the processing of GNSS/SINS-IADPS data processing (Ryzhkov, 2021; Michał et al., 2022). However, in practical applications, the mathematical model of GNSS/SINS-IADPS is frequently nonlinear. As such, investigating quaternion-based nonlinear filter algorithms is pivotal for enhancing the efficacy of GNSS/SINS-IADPS data processing (Ali and Mirza, 2010; Zhu et al., 2021).

For a long time in the past, the quaternion extended Kalman filter (QEKF) algorithm has been an important method of quaternion based nonlinear filter algorithm, but this approximately linearized suboptimal filter algorithm has the defect of high-order truncation errors (Gui and de Ruiter, 2018). To overcome the limitations of QEKF algorithm, the quaternion unscented Kalman filter (QUKF) algorithm approximates the probability density distribution of the nonlinear system through the sigma point set, avoiding the linearization errors and solving the Jacobi matrix. Although the filtering estimation performance is improved, frequent switching of quaternions and modification of Rodrigues parameters during iteration result in a large computational burden (Julier et al., 2000; Chang et al., 2013; Challa et al., 2016). The Cubature Kalman filter (CKF) algorithm employs a set of cubature points generated by the third-order sphere-phase diameter cubature rule to approximate the probability density distribution of nonlinear system, which has a more rigorous theoretical basis and better numerical stability than the unscented Kalman Filter algorithm (Arasaratnam et al., 2010; Zhang et al., 2019; Chang et al., 2021). By combining it with the attitude quaternion method, the obtained quaternion cubature Kalman filter (QCKF) algorithm shows the characteristics of simple implementation, good convergence, high precision, and suitability for high-dimensional systems (Geng et al., 2021; Swati, 2022; Wang et al., 2023). Furthermore, the quaternion augmented cubature Kalman filter (QACKF) algorithm enhances the estimation accuracy of QCKF algorithm to some extent, but its computational complexity is significantly higher than that of QCKF algorithm (Wang et al., 2017). For this reason, Huang et al. studied the quaternion state switching cubature Kalman filter (QSSCKF) algorithm, which ensured the filtering estimation accuracy while effectively reducing the computation amount (Huang et al., 2020). However, the use of quaternions as state variables is bound by normalization constraints in practical applications. If these constraints are overlooked during data processing, the resulting accuracy of the filter and the positive quality of covariance will be compromised (Huang et al., 2016; Qiu and Qian, 2018). In addition, these previous studies are based on the assumption that the random model is the white Gaussian noise, and the actual process noise and measurement noise in GNSS/SINS-IADPS deviate from ideal Gaussian distribution. That is to say, GNSS/SINS-IADPS is a nonlinear non-Gaussian system. In such case, if the random model is mismatched, it may affect the accuracy of filtering estimation and even lead to accuracy divergence in severe cases (Duong and Chiang, 2012; Wang et al., 2020; Taghizadeh and Safabakhsh, 2023). Therefore, if the QCKF algorithm can be extended to take into account both the effect of non-Gaussian noise and normalization constraints, it will show better estimation performance.

Researchers have proposed various adaptive filter algorithms to address the problem of non-Gaussian noise in state estimation, which are mainly divided into functional model based adaptive filter and stochastic model based adaptive filter (Mohamed and Schwarz, 1999; Chang, 2014; Elmezayen and El-Rabbany, 2021; Jiang et al., 2021; Yang et al., 2021; Wu et al., 2022). The Sage-Husa filter algorithm can obtain real-time statistical data on current epoch process noise and measurement noise. However, when the moving carrier generates a large disturbance, it is difficult for this kind of filter algorithm to distinguish between the model error and measurement noise, thus affecting the estimation results (Song et al., 2022; Chen et al., 2023). The fading filter algorithm makes the algorithm meet the optimality through a fading factor, but this method is limited to dealing with non-Gaussian process noise only (Sun et al., 2022; Wang et al., 2022). The robust adaptive filter algorithm can handle non-Gaussian noise in both process and measurement noise. When these two types of noise are present, the algorithm can achieve the purpose of stable state estimation results by adjusting the adaptive factor and equivalent weight matrix factor, but this method requires redundant measurement values (Dong et al., 2023). Nonetheless, in the context of the aforementioned research, when dealing with non-Gaussian noise problems, the stochastic model of noise is usually covered by a Gaussian distribution with greater variance.

The Gaussian mixture model (GMM) offers an alternative approach to address the issue of non-Gaussian noise through its multi-mode approximation technique. After being processed by its multi-mode approximation method, it has higher accuracy compared to the traditional extended variance Gaussian distribution approximation method (Alspach and Sorenson, 1972). In recent years, the filter random model has been optimized by GMM, which has gradually been recognized as a superior approach, attracting attention in target tracking, speech recognition, signal analysis, navigation integrity monitoring, and other aspects (Sun et al., 2020; Zickert and Yarman, 2021; Zhu et al., 2022; Yu et al., 2023). The literature sequentially delves into the nonlinear optimization problem associated with the Gaussian sum filter (GSF) algorithm grounded in GMM, such as the Gaussian sum extended Kalman filter (GSEKF) algorithm, Gaussian sum unscented Kalman filter (GSUKF) algorithm, Gaussian sum quadrature Kalman Filter (GSQKF) algorithm, and Gaussian sum cubature Kalman filter (GSCKF) algorithm, and these studies have to some extent optimized the integrated navigation information fusion algorithm (Wang and Cheng, 2015; Qian et al., 2021; Wang et al., 2021; Bai et al., 2022). While numerous algorithms have been put forward to optimize the random model of nonlinear filter algorithms under non-Gaussian noise conditions based on GMM, there are limited reports on research on quaternion-based algorithm in GNSS/SINS-IADPS data processing for UAVs.

To tackle the issue of random model mismatch under non-Gaussian noise conditions and quaternions normalization constraints affecting the QCKF algorithm, which leads to degradation in estimation accuracy during GNSS/SINS-IADPS data processing for UAVs, a novel Gaussian sum quaternion constrained cubature Kalman filter (GSQCCKF) algorithm is proposed in this paper. The algorithm combines the GMM principle with the two-step projection method and improves the QCKF algorithm. Firstly, multiple sub-filters are decomposed using GMM. Secondly, the quaternion is restricted by the two-step projection method to achieve the purpose of quaternion normalization in attitude determination. Finally, simulation and experiments in GNSS/SINS-IADPS data processing for UAVs are conducted to verify the improvement of estimation accuracy and adaptability for GSQCCKF. The outcomes demonstrate that the proposed GSQCCKF algorithm significantly mitigates the adverse effects of non-Gaussian noise on state estimation, substantially improving both accuracy and adaptability in the GNSS/SINS-IADPS data processing utilized on UAVs.

2 Preliminaries and problem formulation

2.1 Mathematical models for GNSS/SINS-IADPS

We have adopted an integrated navigation system composed of single antenna GNSS and SINS, which is tightly integrated in combination. This tightly integrated GNSS/SINS-IADPS has better navigation accuracy and anti-interference performance than the individual attitude determination and positioning technology of GNSS and SINS, so it is widely used in numerous scientific fields. In the GNSS/SINS-IADPS, the QEKF algorithm is generally used to fuse GNSS and SINS navigation information. To use the QEKF algorithm, it is customary to use the linearized GNSS/SINS-IADPS mathematical model. However, when the GNSS/SINS-IADPS works in high maneuverability, it will show obvious nonlinear characteristics. At this time, if the GNSS/SINS-IADPS mathematical model is linearized, the estimation accuracy will be reduced because of the linearization error. Therefore, in order to cope with linearization errors, need to establish the nonlinear mathematical model of GNSS/SINS-IADPS. The nonlinear mathematical models consist of two components: state-space equation and measurement equation.

The state estimate of GNSS/SINS-IADPS at epoch $k-1$ is defined in formula (1).

$\begin{array}{ll}{x}_{k-1|k-1}=\left[\delta q\mathrm{,}\delta v\mathrm{,}\delta p\mathrm{,}\epsilon \mathrm{,}\nabla \right]& (1)\end{array}$

here, $\delta q$ is attitude quaternion errors, $\delta v$ is velocity errors, $\delta p$ is position errors, $\epsilon$ is the gyroscope biases drift, $\nabla$ is the accelerometer biases drift. The state-space equation of GNSS/SINS-IADPS is described as follows:

$\begin{array}{ll}{x}_{k|k-1}=f\left(x_{k-1|k-1}\right)+g_k{w}_k& (2)\end{array}$

here, $x_{k|k-1}$ is the state prediction, $f(\cdot )$ is a nonlinear function, $g_k$ is the system noise driven matrix, $w_k$ is the process noise. Assuming that $w_k$ is characterized as white Gaussian noise, it can be represented mathematically as $w_k~\left(0\mathrm{,}{Q}_k\right)$ .

The measurement $z_k$ of the measurement equation is composed of the corrected pseudo-range and pseudo-range rate, which have been adjusted for the satellite clock bias, tropospheric delay, and ionospheric delay. Then, according to the state prediction $x_{k|k-1}$ determined in formula (2) and measurement $z_k$ , the measurement equation of GNSS/SINS-IADPS is established:

$\begin{array}{ll}{z}_k=h\left(x_{k|k-1}\right)+v_k& (3)\end{array}$

here, $h(\cdot )$ is a nonlinear measurement function; $v_k$ is the measurement noise, which is caused by receiver noise, multipath effects, and orbit prediction errors. Assuming that the measurement noise  $v_k$  is characterized as white Gaussian noise, it can be represented mathematically as $v_k~\left(0\mathrm{,}{R}_k\right).$

2.2 Quaternion cubature Kalman filter algorithm

The QCKF algorithm is a Gaussian filter algorithm that estimates the posterior distribution of the probability density function (PDF) of a nonlinear function by utilizing a set of cubature points, thereby circumventing the necessity for linearization of the nonlinear function (Geng et al., 2021; Swati, 2022; Wang et al., 2023). The concrete implementation procedures of the QCKF algorithm are outlined below:

Step 1: The selection of the cubature points ${\chi }_{c,k-1|k-1}$ is accomplished through the utilization of the third-order sphere-phase diameter cubature rule in formulas (4, 5).

$\begin{array}{ll}{\chi }_{c,k-1|k-1}=S_{k-1|k-1}{\xi }_c+x_{k-1|k-1}& (4)\end{array}$

${\xi }_c=\sqrt{\frac{m}{2}}{\left\{\left(\begin{array}{c}1\\ ⋮\\ 0\end{array}\right)\mathrm{,}\cdots \mathrm{,}\left(\begin{array}{c}0\\ ⋮\\ 1\end{array}\right)\mathrm{,}\left(\begin{array}{c}-1\\ ⋮\\ 0\end{array}\right)\mathrm{,}\cdots \mathrm{,}\left(\begin{array}{c}0\\ ⋮\\ -1\end{array}\right)\right\}}_i\text{ (5)}$

here, the Cholesky decomposition of the matrix $S_{k-1|k-1}$ can be obtained as $P_{k-1|k-1}=S_{k-1|k-1}{S}_{k-1|k-1}^{\mathrm{T}}$ , $P_{k-1|k-1}$ is the covariance pertains to the state estimation at epoch $k-1$ , and $c=1,2,\cdots ,m=2n$ , $n$ is the dimension of the state estimation, that is, the overall quantity of cubature points is double the dimensionality of the state estimation.

Step 2: Prediction.

The propagated cubature points $x_{c,k|k-1}^{\ast }$ are estimated through the nonlinear function $f(\cdot )$ , which can be expressed in formula (6).

$\begin{array}{ll}{x}_{c,k|k-1}^{\ast }=f\left({\chi }_{c,k-1|k-1}\right)& (6)\end{array}$

Then, the state prediction $x_{k|k-1}$ at epoch $k$ can be computed, as illustrated in formula (7).

$\begin{array}{ll}{x}_{k|k-1}=\frac{1}{m}{\displaystyle \sum }_{c=1}^m{x}_{c,k|k-1}^{\ast }& (7)\end{array}$

The state prediction covariance $P_{k|k-1}$ can be derived, as depicted in formula (8).

$\begin{array}{ll}{P}_{k|k-1}=\frac{1}{m}{\displaystyle \sum }_{c=1}^m{x}_{c,k|k-1}^{\ast }{x}_{c,k|k-1}^{\ast \mathrm{T}}-x_{k|k-1}{x}_{k|k-1}^{\mathrm{T}}+Q_k& (8)\end{array}$

Step 3: Update.

The cubature points  ${\chi }_{c,k|k-1}$  can be re-estimated in formula (9).

$\begin{array}{ll}{\chi }_{c,k|k-1}=S_{k|k-1}{\xi }_c+x_{k|k-1}& (9)\end{array}$

here, $P_{k|k-1}=S_{k|k-1}{S}_{k|k-1}^{\mathrm{T}}$ .

The propagated cubature points $z_{c,k|}^{\ast }{}_{k-1}$ can be estimated by applying the nonlinear function $h(\cdot )$ , which is expressed in formula (10).

$\begin{array}{ll}{z}_{c,k|}^{\ast }{}_{k-1}=h\left({\chi }_{c,k|k-1}\right)& (10)\end{array}$

Subsequently, the measurement prediction $z_{k|k-1}$ can be computed, as illustrated in formula (11).

$\begin{array}{ll}{z}_{k|k-1}=\frac{1}{m}{\displaystyle \sum }_{c=1}^m{z}_{c,k|}^{\ast }{}_{k-1}& (11)\end{array}$

The covariance $P_{zz,k|k-1}$ of the measurement prediction, cross-covariance $P_{xz,k|k-1}$ of the measurement prediction, and the filter gain $K_k$ can be derived, as depicted in formulas (12)–(14), respectively.

$\begin{array}{ll}{P}_{zz,k|k-1}=\frac{1}{m}{\displaystyle \sum }_{c=1}^m{z}_{c,k|}^{\ast }{}_{k-1}{z}_{c,k|}^{\ast \mathrm{T}}{}_{k-1}-z_{k|k-1}{z}_{k|k-1}^{\mathrm{T}}+R_k& (12)\end{array}$

$\begin{array}{ll}{P}_{xz,k|k-1}=\frac{1}{m}{\displaystyle \sum }_{c=1}^m{\chi }_{c,k|k-1}{z}_{c,k|}^{\ast \mathrm{T}}{}_{k-1}-x_{k|k-1}{z}_{k|k-1}^{\mathrm{T}}& (13)\end{array}$

$\begin{array}{ll}{K}_k=P_{xz,k|k-1}{P}_{zz,k|k-1}^{-1}& (14)\end{array}$

Subsequently, the state estimation $x_{k|k}$ and its covariance $P_{k|k}$ at epoch $k$ can be obtained, which are expressed in formulas (15) and (16).

$\begin{array}{ll}{x}_{k|k}=x_{k|k-1}+K_k\left(z_k-z_{k|k-1}\right)& (15)\end{array}$

$\begin{array}{ll}{P}_{k|k}=P_{k|k-1}-K_k{P}_{zz,k|k-1}{K}_k^{\mathrm{T}}& (16)\end{array}$

2.3 The limitation of QCKF algorithm

The wavelet transform method is employed to analyze the statistical properties of SINS errors, and Allan variance analysis is utilized to scrutinize the statistical characteristics of GNSS residuals (El-Sheimy et al., 2008; Wang et al., 2018; Zhang et al., 2020). From the analysis results, it was found that SINS errors and GNSS residuals do not conform to the distribution of zero-mean white Gaussian noise. Instead, the SINS errors and GNSS residuals remnants are a mixed distribution of Gaussian noise and non-Gaussian noise (El-Sheimy et al., 2008; Wang et al., 2018; Zhang et al., 2020). Since the QCKF algorithm is based on the white Gaussian noise hypothesis, this hypothesis may lead to suboptimal estimation results for GNSS/SINS-IADPS due to the random model mismatch inherent in the QCKF algorithm.

In addition, the quaternion normalization problem exists in the attitude quaternion of state estimation. Neglecting the quaternion constraint during filtering calculations may result in a decline in the estimation accuracy, potentially leading to covariance singularity.

Therefore, to enhance the QCKF algorithm estimation performance of the GNSS/SINS-IADPS under non-Gaussian noise environments, it is necessary to further refine the random model of QCKF algorithm and optimize its algorithm model employed in the GNSS/SINS-IADPS.

3 Gaussian sum quaternion constrained cubature Kalman filter algorithm

The QCKF algorithm based on the assumption of white Gaussian noise is difficult to obtain ideal performance in state estimation due to the influence of random model mismatch and the oversight of quaternion normalization. In this section, a novel GSQCCKF algorithm is proposed, which is based on the QCKF algorithm framework. The aim of this algorithm is to address the state estimation problem of QCKF algorithm used in non-Gaussian noise environment for GNSS/SINS-IADPS data processing. The steps involved in the GSQCCKF algorithm are described as follows.

3.1 Modeling of non-Gaussian probability density function by GMM

Non-Gaussian noise can be modeled as a multi-component system based on the degree of nonlinearity or the maximum eigenvalue of the covariance matrix (Maebashi et al., 2008; Liu et al., 2014; Qian et al., 2021). The PDF of $x_{k-1|k-1}$ in formula (1) can be approximately expressed as follows:

$\begin{array}{l}p\left(x_{k-1|k-1}\right)\approx {\displaystyle \sum }_{i=1}^I{\omega }_{k-1}\left(i\right)\mathrm{N}\left(\begin{array}{l}{x}_{k-1|k-1}\mathrm{;}{\overline{\mu }}_{k-1}\left(i\right)\mathrm{,}\\ P_{k-1|k-1}\left(i\right)\end{array}\right),\\ \mathrm{}{\displaystyle \sum }_{i=1}^I{\omega }_{k-1}\left(i\right)=1\end{array}\text{ (17)}$

here, ${\omega }_{k-1}\left(i\right)$ denotes the weight of the $ith$ Gaussian component, $\mathrm{N}\left(x_{k-1|k-1}\mathrm{;}{\overline{\mu }}_{k-1}\left(i\right)\mathrm{,}{P}_{k-1|k-1}\left(i\right)\right)$ represents the $ith$ Gaussian component with a mean of ${\overline{\mu }}_{k-1}\left(i\right)$ and a variance of $P_{k-1|k-1}\left(i\right)$ , $I$ signifies the total number of Gaussian components.

Correspondingly, the PDF of the process noise $p\left(w_k\right)$ in formula (2) and the PDF of the measurement noise $p\left(v_k\right)$ in formula (3) can be approximately expressed as follows:

$\begin{array}{ll}p\left(w_k\right)\approx {\displaystyle \sum }_{j=1}^J{\alpha }_k\left(j\right)N\left(w_k\mathrm{;}{\overline{w}}_k\left(j\right)\mathrm{,}{Q}_k\left(j\right)\right),{\displaystyle \sum }_{i=1}^J{\alpha }_k\left(j\right)=1& (18)\end{array}$

$\begin{array}{ll}p\left(v_k\right)\approx {\displaystyle \sum }_{l=1}^L{\beta }_k\left(l\right)\mathrm{N}\left(v_k\mathrm{;}{\overline{v}}_k\left(l\right)\mathrm{,}{R}_k\left(l\right)\right),{\displaystyle \sum }_{l=1}^L{\beta }_k\left(l\right)=1& (19)\end{array}$

here, ${\alpha }_k\left(j\right)$ represents the weight of the $jth$ Gaussian component, $N\left(w_k\mathrm{;}{\overline{w}}_k\left(j\right)\mathrm{,}{Q}_k\left(j\right)\right)$ represents the $jth$ Gaussian component with a mean of ${\overline{w}}_k\left(j\right)$ and a variance of $Q_k\left(j\right)$ , $J$ represents the total number of Gaussian components; ${\beta }_k\left(l\right)$ represents the weight of the $lth$ Gaussian component, $\mathrm{N}\left(v_k\mathrm{;}{\overline{v}}_k\left(l\right)\mathrm{,}{R}_k\left(l\right)\right)$ represents the $lth$ Gaussian component with a mean of ${\overline{v}}_k\left(l\right)$ and a variance of $R_k\left(l\right)$ , $L$ represents the total number of Gaussian components.

3.2 Gaussian sum quaternion cubature Kalman filter algorithm

Following the implementation of the non-Gaussian PDF through GMM, the prediction and update on each sub-filters are carried out.

Step 1: Prediction.

The cubature points ${\chi }_{c,k-1|k-1}$ are formulated utilizing the third-order sphere-phase diameter cubature rule.

$\begin{array}{ll}{\chi }_{c,k-1|k-1}\left(i\right)=S_{k-1|k-1}\left(r\right)\cdot {\xi }_c+x_{k-1|k-1}\left(i\right)& (20)\end{array}$

here, $S_{k-1|k-1}\left(r\right)$ is derived through Cholesky decomposition of $P_{k-1|k-1}\left(r\right)$ , which can be expressed as $P_{k-1|k-1}\left(r\right)=S_{k-1|k-1}\left(r\right)S_{k-1|k-1}^{\mathrm{T}}\left(r\right)$ , $r=\left(i-1\right)\cdot I+j$ .

By propagating the cubature points ${\chi }_{c,k-1|k-1}\left(i\right)$ through the nonlinear function $f(\cdot )$ , we can obtain ${\chi }_{c,k|k-1}^{\ast }\left(i\right)=f\left({\chi }_{c,k-1|k-1}\left(i\right)\right).$ The state prediction $x_{k|k-1}\left(r\right)$ is calculated as follows:

$\begin{array}{ll}{x}_{k|k-1}\left(r\right)={\displaystyle \sum }_{c=1}^m{\omega }_c{\chi }_{c,k|k-1}^{\ast }\left(i\right)+{\overline{w}}_k\left(j\right)& (21)\end{array}$

Subsequently, the covariance of state prediction $P_{k|k-1}\left(r\right)$ can be derived, as depicted in formula (22).

$\begin{array}{l}{P}_{k|k-1}\left(r\right)={\displaystyle \sum }_{c=1}^m{\omega }_c{\chi }_{c,k|k-1}^{\ast }\left(i\right){\chi }_{c,k|k-1}^{\ast }{}^{\mathrm{T}}\left(i\right)-\\ \mathrm{}[x_{k|k-1}\left(r\right)-{\overline{w}}_k\left(j\right)]\\ \mathrm{}{\left[x_{k|k-1}\left(r\right)-{\overline{w}}_k\left(j\right)\right]}^{\mathrm{T}}+Q_k\left(j\right)\end{array}\text{ (22)}$

here, ${\omega }_c$ denotes the weight of the cubature points.

Step 2: Update.

The cubature points are assessed:

$\begin{array}{ll}{\chi }_{c,k|k-1}\left(r\right)=S_{k|k-1}\left(r\right){\xi }_c+x_{k|k-1}\left(r\right)& (23)\end{array}$

Here,

$P_{k|k-1}\left(r\right)=S_{k|k-1}\left(r\right)S_{k|k-1}^{\mathrm{T}}\left(r\right).$

By propagating the cubature points ${\chi }_{c,k|k-1}\left(r\right)$ through the nonlinear function $h(\cdot )$ , the propagated cubature points can be obtained as: $z_{c,k|k-1}^{\ast }\left(r\right)=h\left({\chi }_{c,k|k-1}\left(r\right)\right).$  The measurement prediction $z_{k|k-1}\left(r\mathrm{,}l\right)$ is calculated as follows:

$\begin{array}{ll}{z}_{k|k-1}\left(r\mathrm{,}l\right)={\displaystyle \sum }_{c=1}^m{\omega }_c{z}_{c,k|k-1}^{\ast }\left(r\right)+{\overline{v}}_k\left(l\right)& (24)\end{array}$

The covariance matrix of measurement prediction $P_{zz,k|k-1}\left(r\mathrm{,}l\right)$ , the cross-covariance $P_{xz,k|k-1}\left(r\mathrm{,}l\right)$ , the filter gain $K_k\left(r\mathrm{,}l\right)$ , the state estimation $x_{k|k}\left(g\right)$ and its corresponding covariance $P_{k|k}\left(g\right)$ can be derived, as depicted in formulas (25)–(29), respectively.

$\begin{array}{l}{P}_{zz,k|k-1}\left(r\mathrm{,}l\right)={\displaystyle \sum }_{c=1}^m{\omega }_c{\chi }_{c,k|k-1}\left(r\right){\chi }_{c,k|k-1}^{\mathrm{T}}\left(r\right)-\\ \mathrm{}\left(z_{k|k-1}\left(r\mathrm{,}l\right)-\overline{v}\left(l\right)k\right)\\ \mathrm{}{\left(\begin{array}{l}{z}_{k|k-1}\left(r\mathrm{,}l\right)-\\ \overline{v}\left(l\right)k\end{array}\right)}^{\mathrm{T}}+R_k\left(l\right)\end{array}\text{ (25)}$

$\begin{array}{l}{P}_{xz,k|k-1}\left(r\mathrm{,}l\right)={\displaystyle \sum }_{c=1}^m{\omega }_c{\chi }_{c,k|k-1}^{\ast }\left(r\right)z_{c,k|k-1}^{\ast }{\left(r\right)}^{\mathrm{T}}-\\ \mathrm{}{x}_{k|k-1}\left(r\right){\left[z_{k|k-1}\left(r\mathrm{,}l\right)-\overline{v}\left(l\right)k\right]}^{\mathrm{T}}\end{array}\text{ (26)}$

$\begin{array}{ll}{K}_k\left(r\mathrm{,}l\right)=P_{xz,k|k-1}\left(r\mathrm{,}l\right)P_{zz,k|k-1}^{-1}\left(r\mathrm{,}l\right)& (27)\end{array}$

$\begin{array}{ll}{x}_{k|k}\left(g\right)=x_{k|k-1}\left(r\right)+K_k\left(r\mathrm{,}l\right)\left[z_k-z_{k|k-1}\left(r\mathrm{,}l\right)\right]& (28)\end{array}$

$\begin{array}{ll}{P}_{k|k}\left(g\right)=P_{k|k-1}\left(r\right)-K_k\left(r\mathrm{,}l\right)P_{zz,k|k-1}\left(r.l\right)K_k^{\mathrm{T}}\left(r\mathrm{,}l\right)& (29)\end{array}$

Here,

$g=\left(r-1\right)L+l.$

Step 3: Global point estimate.

The state estimation $x_{k|k}$ and its covariance matrix $P_{k|k}$ are computed as

$\begin{array}{ll}{x}_{k|k}={\displaystyle \sum }_{g=1}^{I\cdot J\cdot L}{\omega }_k\left(g\right)x_{k|k}\left(g\right)& (30)\end{array}$

$\begin{array}{l}{P}_{k|k}={\displaystyle \sum }_{g=1}^{I\cdot J\cdot L}{\omega }_k\left(g\right)(P_{k|k}\left(g\right)+\\ \mathrm{}\left(x_{k|k}\left(g\right)-x_{k|k}\right){\left(x_{k|k}\left(g\right)-x_{k|k}\right)}^{\mathrm{T}})\end{array}\text{ (31)}$

here, ${\omega }_k\left(g\right)$ denotes the weight of the $gth$ Gaussian component, which is calculated as formula (32).

$\begin{array}{ll}{\omega }_k\left(g\right)=\frac{{\alpha }_{k|k-1}\left(r\right){\beta }_k\left(l\right)p\left(z_{k|k}|x_{k|k}\mathrm{,}g\right)}{{\displaystyle \sum }_{r=1}^{I\cdot J}{\displaystyle \sum }_{l=1}^L{\alpha }_{k|k-1}\left(r\right){\beta }_k\left(l\right)p\left(z_{k|k}|x_{k|k}\mathrm{,}g\right)}& (32)\end{array}$

In this regard, the PDF of the $gth$ Gaussian component is given by, and its computation can be expressed as follows:

$\begin{array}{ll}p\left(z_k|x_{k|k}\mathrm{,}g\right)=\frac{1}{\sqrt{2\pi {\sigma }_g^2}}\exp \left(-\frac{1}{2}{\left(\frac{z_k-z\left(r\mathrm{,}l\right)k|k-1}{{\sigma }_g}\right)}^2\right)& (33)\end{array}$

From formulas (32) and (33), it is evident that the weight ${\omega }_k\left(g\right)$ of the Gaussian component adaptively modifies in response to the innovation $z_k-z\left(r\mathrm{,}l\right)k|k-1$ , thereby enhancing the robustness of the filter algorithm.

Step 4: Gaussian component reduction.

From formula (20) to (30), it can be found that after step 3 the number of Gaussian components reaches $I\cdot J\cdot L$ . If $I\cdot J\cdot L>I$ , there will be a mismatch between the number of Gaussian components at epoch k and epoch k − 1 when the recurrence operation is performed at epoch k. It can be seen that the total number of Gaussian components will increase with each iteration of filtering, which will eventually lead to an exponential increase in filtering recursion. Therefore, it is necessary to reduce the number of Gaussian components after each iteration, this ensures that the total number of Gaussian components in the state estimation remains $I$ . Gaussian component merging generally employs a quasi-Bayesian approximation, but the threshold selection depends on experience. Here, we adopt a different approach, by arranging the weight values of the Gaussian components in descending order, the Gaussian components are then sequentially labeled as $g=1,\cdots ,I\cdot J\cdot L$ , while retaining $I-1$ components, then the weight value ${\omega }_k\left(I\right)$ of the $I-1$ Gaussian component at epoch $k$ can be calculated:

$\begin{array}{ll}{\omega }_k\left(I\right)={\displaystyle \sum }_{g=1}^{J\cdot L}{\omega }_k\left(g\right)& (34)\end{array}$

The mean of state estimation $x_{k|k}\left(I\right)$ and its covariance $P_{k|k}\left(I\right)$ corresponding to the $Ith$ Gaussian component is shown in formulas (35) and (36), respectively.

$\begin{array}{ll}{x}_{k|k}\left(I\right)={{\displaystyle \sum }_{n=1}^{J\cdot L}}_k\tilde{\omega }\left(g\right)x_{k|k}\left(g\right)& (35)\end{array}$

$\begin{array}{l}{P}_{k|k}\left(I\right)={\displaystyle \sum }_{n=1}^{J\cdot L}{\tilde{\omega }}_k\left(g\right)[P\left(g\right)k|k+\left(x_{k|k}\left(g\right)-x_{k|k}\left(I\right)\right)\\ \mathrm{}{\left(x_{k|k}\left(g\right)-x_{k|k}\left(I\right)\right)}^T]\end{array}\text{ (36)}$

here, ${\tilde{\omega }}_k\left(g\right)={\omega }_k\left(g\right)/{\omega }_k\left(I\right)$ denotes the regularization weight. Through the aforementioned procedures, $I-1$ Gaussian components and the $Ith$ Gaussian component as described in formulas (35) and (36) are employed for the filtering recursion in the subsequent epoch. It is evident that following the merging of Gaussian components, the overall number of Gaussian components within the filter algorithm remains unchanged.

3.3 Two-step projection method

The two-step projection method (Tang et al., 2012; Huang, 2017) is employed herein to address quaternion constraint problem. In the initial step, the unconstrained state estimation distribution is projected onto the constrained surface so that the attitude estimation distribution complies with the quaternion constraint. However, this action may result in a reduction in the variance of attitude estimation. In the subsequent step, the constrained state estimation distribution is projected onto the constrained surface so that the mean of the attitude estimation satisfies the quaternion constraint. Simultaneously, the attitude estimation variance is compensated, to enhance the accuracy of attitude estimation while ensuring quaternion normalization. The specific processing steps are outlined as follows:

Step 1: the mean of constrained state estimation distribution ${\overrightarrow{x}}_{k|k}$ and its covariance ${\overrightarrow{P}}_{k|k}$ can be computed, which are detailed below:

$\begin{array}{ll}{\overrightarrow{x}}_{k|k}={\displaystyle \sum }_{c=1}^{2n}{\omega }_c{\tilde{\chi }}_{c,k|k}& (37)\end{array}$

$\begin{array}{ll}{\overrightarrow{P}}_{k|k}={\displaystyle \sum }_{c=1}^{2n}{\omega }_c\left({\tilde{\chi }}_{c,k|k}-{\overrightarrow{x}}_{k|k}\right){\left({\tilde{\chi }}_{c,k|k}-{\overrightarrow{x}}_{k|k}\right)}^{\mathrm{T}}& (38)\end{array}$

here, ${\tilde{\chi }}_{c,k|k}=\phi \left({\chi }_{c,k|k}^{\ast }\right),$ ${\chi }_{c,k|k}^{\ast }={\overline{S}}_{k|k}{\xi }_c+x_{k|k},$ $P_{k|k}={\overline{S}}_{k|k}{\overline{S}}_{k|k}^{\mathrm{T}}$ , $\phi (\cdot )$ is the constraint function.

Step 2: the ultimate state estimation ${\tilde{x}}_{k|k}$ and its covariance ${\tilde{P}}_{k|k}$ can be calculated:

$\begin{array}{ll}{\tilde{x}}_{k|k}=\phi \left({\overrightarrow{x}}_{k|k}\right)& (39)\end{array}$

$\begin{array}{ll}{\tilde{P}}_{k|k}={\overrightarrow{P}}_{k|k}+\left({\tilde{x}}_{k|k}-{\overrightarrow{x}}_{k|k}\right){\left({\tilde{x}}_{k|k}-{\overrightarrow{x}}_{k|k}\right)}^{\mathrm{T}}& (40)\end{array}$

3.4 GSQCCKF algorithm structure

The GSQCCKF algorithm is implemented in five stages, as illustrated in Figure 1.

Figure 1

Figure 1. The flow diagram of proposed GSQCCKF algorithm.

Stage 1: State estimation, process noise, and measurement noise are modeled by GMM, as shown in formulas (17)–(19).

Stage 2: The prediction of GSQCKF algorithm at epoch $k$ is calculated, as shown in formulas (20)–(22).

Stage 3: The Update of GSQCKF algorithm at epoch $k$ is calculated, as shown in formulas (23)–(29).

Stage 4: Global point estimate (formulas 30–33) and then perform Gaussian component reduction (formulas 34–36).

Stage 5: The two-step projection method is applied to deal with the quaternion normalization issue, and the output of state estimation at epoch $k$ , as shown in formulas (39) and (40).

4 Performance evaluation and discussion

In this section, the performance of the GSQCCKF algorithm was assessed, then, a comparative analysis was conducted on the performance of four different algorithms (GSQCCKF, QEKF, QCKF, and QCCKF).

4.1 Simulations and analysis

The performance of proposed GSQCCKF algorithm was evaluated by the Monte Carlo simulations, which were based on the design of a dynamic UAV equipped with a GNSS/SINS-IADPS. The flight trajectory of the UAV is depicted in Figure 2, which includes a variety of maneuvering states, such as climbing, level flight, turning and descending. The initial attitude of UAV was set as roll 0°, pitch 0°, and yaw 0°, while its initial speed was specified as 0 m/s in East, 120 m/s in North, and 0 m/s in Up. The initial position of UAV was established at 110.2° longitude, 34.0° latitude, and 2000 m altitude. Additionally, the initial attitude error of UAV was specified as roll 1 “, pitch 1 “, yaw 1.5 “, while its initial velocity error was set at 0.3 m/s in the East, 0.3 m /s in the North, 0.3 m/s in the Up direction, and the initial position error of UAV was defined as 10 m longitude, 10 m latitude, and 10 m altitude. The SINS parameter configurations of the GNSS/SINS-IADPS are presented in Table 1. The GNSS measurement was simulated based on the satellite constellation and epoch information on June 13, 2023. The pseudo range observation error of the GNSS receiver was $10\mathrm{m}$ , and the sampling frequency of the GNSS receiver was 1 Hz.

Figure 2

Figure 2. The UAV flight trajectory.

Table 1

Table 1. Parameters of SINS.

White Gaussian noise scenario and non-Gaussian noise scenario were simulated separately. For each scenario, 100 Monte Carlo simulations were carried out, and each simulation time was 1,500 s. The configuration of the computing platform used in the simulation is as follows. CPU: Inter Core i7-12700, 2.9GHZ; Internal memory: DDR4 16GB; Simulation software: Matlab R2020b. The precision of each filter estimation was assessed using the root mean square error (RMSE), which is mathematically expressed in formula (41).

$\begin{array}{ll}\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left({\overline{X}}_i-X_i\right)}^2}& (41)\end{array}$

here, $N$ signifies the total number of Monte Carlo simulations; $X_i$ represents the reference; ${\overline{X}}_i$ signifies the estimation.

4.1.1 White Gaussian noise scenario

In this scenario, both process noise and observation noise are simulated as white Gaussian noise, the covariance of which are defined in formula (42).

$Q=diag\left[\begin{array}{l}{\left({0.01}^{\circ }/\sqrt{\mathrm{h}}\right)}^2{I}_{3\times 3}\mathrm{,}{\left(1\times 0.0001g\cdot \sqrt{\mathrm{s}}\right)}^2\\ I_{3\times 3}\mathrm{,}{0}_{9\times 9}\mathrm{,}10\mathrm{m}/\sqrt{\mathrm{s}}\mathrm{,}5\mathrm{m}/{\mathrm{s}}^{3/2}\end{array}\right]\text{ (42)}$

$R={\left(25\mathrm{m}\right)}^2{I}_{4\times 4}$

The RMSE of attitude and position results calculated by four different algorithms (QEKF, QCKF, QCCKF, and GSQCCKF) under white Gaussian noise scenario are shown in Figure 3, where the epoch range ranges from 0 s to 1500s.

Figure 3

Figure 3. RMSEs of the attitude and position obtained by four different algorithms for the white Gaussian noise scenario. (A) Attitude. (B) Position.

Additionally, the average RMSEs (ARMSEs) for each algorithm are also computed and presented in Figure 4 for ease of comparison.

Figure 4

Figure 4. ARMSEs of the position and attitude obtained by four different algorithms for the white Gaussian noise scenario. (A) Attitude. (B) Position.

As is evident from Figures 3, 4, QEKF, QCKF, QCCKF, and GSQCCKF, these four algorithms all demonstrate convergence in the white Gaussian noise scenario. Our analysis reveals that the QEKF algorithm exhibits the highest estimation error. This can be attributed to the fact that when the initial attitude error is present in the QEKF algorithm, the high-order truncation error within the filter leads to a reduction in estimation accuracy and filtering stability. In contrast, algorithms of QCKF, QCCKF, and GSQCCKF are all derived from the QCKF algorithm framework. The QCKF algorithm computes the state estimation and its covariance using a set of cubature points, thus enabling algorithms of QCKF, QCCKF, and GSQCCKF to mitigate the impact of nonlinear function linearization errors inherent in the QEKF algorithm on estimation accuracy.

In comparison to the QCKF algorithm, both the QCCKF algorithm and the GSQCCKF algorithm exhibit a marked enhancement in estimation accuracy. This is because both of them consider the quaternion constraint in the state estimation, so the filter gain calculated by these two algorithms is also constrained, which further improves the estimation accuracy of attitude and position. It is noteworthy that the computational accuracy of the QCCKF algorithm and the GSQCCKF algorithm are similar to each other. This outcome is attributable to the fact that the QCCKF algorithm is based on the assumption of white Gaussian noise, while the GSQCCKF represents an enhancement of the QCCKF, grounded in the GMM approach. In white Gaussian noise scenario, both the QCCKF algorithm and the GSQCCKF algorithm are capable of converging, resulting in similar estimation accuracies.

4.1.2 Non-Gaussian noise scenario

In non-Gaussian noise scenario, non-Gaussian process noise is modeled and generated by the model $0.9\mathrm{N}\left(0\mathrm{,}{Q}_k\right)+\left(1-0.9\right)\mathrm{N}\left(0\mathrm{,}10Q_k\right),$ while the non-Gaussian measurement noise is also modeled and generated by the model $0.9\mathrm{N}\left(0\mathrm{,}{R}_k\right)+\left(1-0.9\right)\mathrm{N}\left(0\mathrm{,}10R_k\right).$ The RMSEs of attitude and position which result from four different algorithms (QEKF, QCKF, QCCKF, and GSQCCKF) are shown in Figure 5 and their corresponding ARMSEs are illustrated in Figure 6, where the epoch range ranges from 0 s to 1500s.

Figure 5

Figure 5. RMSEs of the attitude and position obtained by four different algorithms for the non-Gaussian noise scenario. (A) Attitude. (B) Position.

Figure 6

Figure 6. ARMSEs of the attitude and position obtained by four different algorithms for the non-Gaussian noise scenario. (A) Attitude. (B) Position.

As can be observed from Figures 5, 6, the estimation errors of QEKF algorithm, QCKF algorithm, and QCCKF algorithm increase significantly in the non-Gaussian noise scenario when compared to the white Gaussian noise scenario. Taking the yaw estimation results as example, the yaw ARMSEs calculated by algorithms of QEKF, QCKF, and QCCKF increased from 0.39′, 0.31′, 0.29′ to 0.61′, 0.47′, 0.44′, with error increase rates approximately 56.41, 52.61, and 51.72%, respectively. In contrast, the variation of the estimation error of GSQCCKF algorithm is the smallest. And the ARMSE calculated by the GSQCCKF algorithm changes from 0.28′ to 0.36′, resulting in an error increase rate of about 28.57%. It indicates that the GSQCCKF algorithm exhibits the highest estimation accuracy in the presence of non-Gaussian noise. The reason for this phenomenon lies in the fact that the GSQCCKF algorithm enhances the filtering estimation accuracy by refining the random model through GMM, and this optimization method can make the adaptability of the algorithm more effectively to mitigate the impact of non-Gaussian noise on the GNSS/SINS-IADPS data processing.

4.1.3 Computational performance

The computational time per epoch run, for four different algorithms (QEKF, QCKF, QCCKF, and GSQCCKF), as shown in Table 2. In addition, the computational efficiency of them was compared in Figure 7.

Table 2

Table 2. Computational time per epoch run of four different algorithms.

Figure 7

Figure 7. Relative computational efficiencies of four different algorithms. (A) White Gaussian noise scenario. (B) Non-Gaussian noise scenario.

As depicted in Table 2 and Figure 7, it reveals that the changes of computation time for all these four algorithms are relatively similar in white Gaussian noise scenario and non-Gaussian noise scenario. The QEKF algorithm exhibits the shortest computational time. However, the computational time of QCKF algorithm is at least 27.57% longer than that of QEKF algorithm, due to the complex and time-consuming process involved in cubature transformation calculations. Further, the computational time of the QCCKF algorithm is at least 15.19% larger than that of the QCKF algorithm owing to the calculation of the quaternion constraint. Due to the complexity of the computational process, the GSQCCKF algorithm takes the longest computational time, approximately twice that of QCCKF algorithm, as it performs distributed filtering and global point estimation at each iteration. Fortunately, the Gaussian components are trimmed and merged during the processing, preventing exponential growth in computational time (accounting for approximately 285.76% of QEKF). Furthermore, with the significant increase in computing power today, the millisecond-level operation time of GSQCCKF algorithm (16.43 ms) can still meet the real-time data processing requirements of GNSS/SINS-IADPS. In conclusion, despite its increased computational time, the GSQCCKF algorithm remains capable of handling high-dynamic navigation for GNSS/SINS-IADPS equipped on UAVs.

The simulation analysis and comparison conducted in sections 4.1.1, 4.1.2, and 4.1.3 reveal that the proposed GSQCCKF algorithm can effectively refine the random model of filter algorithm and mitigate the impact of non-Gaussian noise on the estimation performance in the GNSS/SINS-IADPS under non-Gaussian noise conditions. As a result, the GSQCCKF algorithm exhibits a higher level of computational accuracy and adaptability when compared to algorithms of QEKF, QCKF, and QCCKF. Despite the increased computational time, the GSQCCKF algorithm remains suitable for real-time solution of GNSS/SINS-IADPS under dynamic navigation states.

4.2 Experiments and analysis

The performance of GSQCCKF algorithm was evaluated by experiments using UAV that involved a diverse range of maneuvers. The experimental data were collected continuously for a duration of 70 min on September 15, 2023, in Zhengzhou, China (114.0°E, 34.3°N).

The UAV is equipped with a GNSS/SINS-IADPS, the parameters of which are detailed in Table 3. The GNSS reference station is located on the ground, with a maximum distance of 20 km from the UAV. The UAV is also equipped with a GNSS receiver, which processes the paired data between it and the GNSS reference station to obtain the differential GPS (DGPS) data with an accuracy of better than 0.1 m. This DGPS data serves as a reference value for assessing the performance of different algorithms.

Table 3

Table 3. System parameters of GNSS/SINS-IADPS.

The starting position of the UAV was 34.654° latitude, 109.193° longitude, and an altitude of 2,683 meters. The initial velocity for the eastern, northern, and up direction are 180 m/s, 60 m/s, and 40 m/s, respectively. Other initial parameters are the same as the simulation. Four different algorithms (QEKF, QCKF, QCCKF, and GSQCCKF) were, respectively, used for GNSS/SINS-IADPS data processing. The test accuracy is measured by 3D positioning error, calculated in formula (43).

$\begin{array}{ll}\Delta p=\sqrt{\Delta {\lambda }^2+\Delta L^2+\Delta h^2}& (43)\end{array}$

here, $\Delta \lambda$ is the positioning error in longitude; $\Delta L$ is the positioning error in latitude; $\Delta h$ is the positioning error in altitude.

To ascertain whether the process noise and measurement noise encountered during experiments are non-Gaussian noise, Allan variance analysis is carried out on the inertial element, the results of which are depicted in Figure 8. The findings indicate that the noise of the inertial element exhibited by the inertial element utilized in the experiments is not white Gaussian noise, but rather a complex noise term encompassing angle random walk, rate slope, quantization noise, and bias instability.

Figure 8

Figure 8. Accelerometer Allan variance results.

The statistical properties of the pseudo-range noise of satellites with different cutoff angle are shown in Table 4 and Figure 9. Notably, G14 has a higher cutoff angle while G22 possesses a lower cutoff angle. As observed from Table 4 and Figure 9, the kurtosis value of the pseudo noise of the G22 satellite is significantly less than 3, that is, it shows negative kurtosis. Consequently, it can be deduced that the pseudo-range noise of G22 exhibits pronounced non-Gaussian characteristics.

Table 4

Table 4. Statistical characteristics of the pseudo-range noise.

Figure 9

Figure 9. Non-Gaussian characteristics of the pseudo-range noise.

According to Figure 8 and Table 4 as well as Figure 9, the SINS errors and the pseudo-range noise of GNSS present in this experiment data are not white Gaussian noise, but rather non-Gaussian noise.

The 3D positioning error curves of four different algorithms (QEKF, QCKF, QCCKF, and GSQCCK) for the epoch range ranging from 100s to 1100s are illustrated in Figure 10. The 3D positioning ARMSEs of different algorithms based on 1,000 sets of data and 4,000 sets of data are shown in Figure 11. It is worth noting that the epoch ranges for the 1,000 sets of data encompass the epoch from 100 s to 1100s, while the epoch range for the 4,000 sets of data spans from 100 s to 4,100 s.

Figure 10

Figure 10. 3D positioning errors of UAV for experiments.

Figure 11

Figure 11. RMSEs of UAV 3D positioning for experiments with different datasets.

As depicted in Figures 10, 11, the period ranges from 100 s to 1100s reveals that the QEKF algorithm is susceptible to linearization errors, resulting in a substantial RMSE of approximately 17.01 m. In contrast, the QCKF algorithm employs a set of cubature points to compute the mean and its covariance, thereby mitigating the linearization error associated with the nonlinear function. Consequently, the QCKF algorithm exhibits a more accurate positioning performance compared to the QEKF algorithm, with the RMSE of approximately 16.38 m. Moreover, because the QCCKF algorithm considers the constraint condition of quaternion, the RMSE is slightly reduced compared with the QCKF algorithm, reaching 16.07 m. Notably, the maximum variation range of the 3D positioning error curves of GSQCCKF algorithm is narrower than that of observed in algorithms of QEKF, QCKF, and QCCKF. This can be attributed to the GSQCCKF algorithm’s use of GMM to accurately modeling the random model, resulting in the smallest RMSE of approximately 15.45 m. Therefore, the same conclusion as the simulation can be obtained, that is, compared with the other three algorithms (QEKF, QCKF, QCCKF), the GSQCCKF algorithm can achieve the best estimation accuracy and adaptability in GNSS/SINS-IADPS data processing.

It can also evident from Figure 11 that, an increase in the experiment data from 1,000 sets to 4,000 sets, the estimation accuracy of four different algorithms (QEKF, QCKF, QCCKF and GSQCCKF) is reduced (QEKF from 17.01to 17.76 m, QCKF from 16.38 to 16.82 m, QCCKF from 16.07 to 16.54 m, and GSQCCKF from 15.45 to 15.63 m). Notably, the estimation accuracy of GSQCCKF algorithm is always better than the other three filter algorithms. When the experiment data comprises 1,000 sets, the 3D positioning RMSEs of GSQCCKF algorithm are about 9.18, 5.545, and 3.70% higher than those of QEKF, QCKF and QCCKF, respectively. As the experiment data increases to 4,000 sets, the accuracy advantage of GSQCCKF algorithm becomes even more pronounced, which is about 11.01, 6.89, and 5.27% higher than those of QEKF, QCKF, and QCCKF, respectively. This shows that the GSQCCKF algorithm possesses robust processing capabilities for non-Gaussian noise and significantly enhances the GNSS/SINS-IADPS estimation accuracy. Furthermore, the GSQCCKF algorithm maintains high accuracy in long-sailing GNSS/SINS-IADPS applications.

5 Conclusion

This paper introduces a novel Gaussian sum quaternion constrained cubature Kalman filter algorithm to tackle the limitations of using QCKF algorithm in the non-Gaussian environments for GNSS/SINS-IADPS. The primary contributions of this research are summarized as follows:

1. The framework of GSQCCKF algorithm is set up. Firstly, the QCKF algorithm based on attitude quaternion for nonlinear/Gaussian systems is presented, followed by an analysis of its limitations. Secondly, the idea of GMM is introduced, which employs a set of Gaussian distributions to approximate the PDF of non-Gaussian variables, including process noise, measurement noise, and state estimation. Thirdly, the combination of QCKF algorithm and GMM produces the proposed GSQCCKF algorithm used for nonlinear non-Gaussian system estimation, which essentially consists of a set of parallel QCKFs.

2. In order to address the quaternion normalization problem in attitude estimation of the proposed GSQCCKF algorithm, a two-step projection method is proposed to resolve the quaternion constraint issue. This approach further enhances the accuracy and numerical stability of GSQCCKF algorithm for GNSS/SINS-IADPS data processing.

Results of simulation and experimentation demonstrate that the proposed GSQCCKF algorithm exhibits a remarkable ability to counteract the adverse effects of non-Gaussian noise on state estimation for GNSS/SINS-IADPS, and it demonstrates significantly enhanced estimation accuracy and adaptability in comparison with algorithms of QEKF, QCKF, and QCCKF.

The GSQCCKF algorithm proposed in this article also has limitations. Theoretically, the GSQCCKF algorithm are unable to adjust to time-varying non-Gaussian noise. In practical terms, due to the non-stationary nature of challenging operation environments of GNSS/SINS-IADPS data processing for UAVs, the non-Gaussian noise may vary over time. The inability to effectively mitigate the impact of such time-varying non-Gaussian noise in the GSQCCKF algorithm could potentially distort the performance of GNSS/SINS-IADPS. Future research efforts should focus on discussing modeling approaches for time-varying non-Gaussian noise in GNSS/SINS-IADPS data processing for UAVs.

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: the data that support the findings of this study are available from the corresponding author upon reasonable request. Requests to access these datasets should be directed to QD, daiqing@lypt.edu.cn.

Author contributions

QD: Funding acquisition, Writing – original draft, Writing – review & editing. G-RX: Methodology, Writing – original draft. G-HZ: Formal analysis, Validation, Writing – review & editing. Q-QY: Visualization, Writing – review & editing. S-YH: Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This research was funded by the National Natural Science Foundation of China (grant no. 42274045), the Henan Province Science and Technology Research Projects (grant no. 242102241067), the Key Research Funding Projects for Higher Education Institutions in Henan Province (grant no. 24A420003), and the Scientific Research Project of Wenzhou University of Technology (grant no. ky202208).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ali, J., and Mirza, M. R. U. B. (2010). Performance comparison among some nonlinear filters for a low cost SINS/GPS integrated solution. Nonlinear Dynam. 61, 491–502. doi: 10.1007/s11071-010-9665-y

Crossref Full Text | Google Scholar

Alspach, L. D., and Sorenson, W. H. (1972). Nonlinear Bayesian estimation using Gaussian sum approximations. IEEE Trans. Autom. Control 17, 439–448. doi: 10.1109/TAC.1972.1100034

Crossref Full Text | Google Scholar

Arasaratnam, I., Haykin, S., and Hurd, R. T. (2010). Cubature Kalman filtering for continuous-discrete systems: theory and simulations. IEEE Trans. Signal Process. 58, 4977–4993. doi: 10.1109/TSP.2010.2056923

Crossref Full Text | Google Scholar

Bai, J. G., Ge, Q. B., Li, H., Xiao, J. M., and Wang, Y. L. (2022). Aircraft trajectory filtering method based on Gaussian-sum and maxi-mum Correntropy Square-root cubature Kalman filter. Cogn. Comput. Syst. 4, 205–217. doi: 10.1049/ccs2.12049

Crossref Full Text | Google Scholar

Boguspayev, N., Akhmedov, D., Raskaliyev, A., Kim, A., and Sukhenko, A. (2023). A comprehensive review of GNSS/INS Integra-tion techniques for land and air vehicle applications. Appl. Sci. 13:4819. doi: 10.3390/app13084819

Crossref Full Text | Google Scholar

Challa, S. M., Moore, G. J., and Rogers, J. D. (2016). A simple attitude unscented Kalman filter: theory and evaluation in a Magne-tometer-only spacecraft scenario. IEEE Access. 4, 1845–1858. doi: 10.1109/ACCESS.2016.2559445

Crossref Full Text | Google Scholar

Chang, G. B. (2014). Robust Kalman filtering based on Mahalanobis distance as outlier judging criterion. J. Geod. 88, 391–401. doi: 10.1007/s00190-013-0690-8

Crossref Full Text | Google Scholar

Chang, L. B., Hu, B. Q., and Chang, G. B. (2013). Modified unscented quaternion estimator based on quaternion averaging. J. Guid. Control. Dyn. 37, 305–309. doi: 10.2514/1.61723

Crossref Full Text | Google Scholar

Chang, Y. Z., Wang, Y. Q., Shen, Y. Y., and Ji, C. G. (2021). A new fuzzy strong tracking cubature Kalman filter for INS/GNSS. GPS Solutions 25:120. doi: 10.1007/s10291-021-01148-5

Crossref Full Text | Google Scholar

Chen, Y., Yan, H., and Li, Y. C. (2023). Vehicle state estimation based on sage–Husa adaptive unscented Kalman filtering. World Electr. Veh. J. 14:167. doi: 10.3390/wevj14070167

Crossref Full Text | Google Scholar

Dong, J. Q., Lian, Z. Z., Xu, J. C., and Yue, Z. (2023). UWB localization based on improved robust adaptive cubature Kalman filter. Sensors 23:2669. doi: 10.3390/s23052669

PubMed Abstract | Crossref Full Text | Google Scholar

Duong, T. T., and Chiang, W. K. (2012). Non-linear, non-Gaussian estimation for INS/GPS integration. Sens. Lett. 10, 1081–1086. doi: 10.1166/sl.2012.2267

Crossref Full Text | Google Scholar

Elmezayen, A., and El-Rabbany, A. (2021). Real-time GNSS precise point positioning using improved robust adaptive Kalman filter. Surv. Rev. 53, 528–542. doi: 10.1080/00396265.2020.1846361

Crossref Full Text | Google Scholar

El-Sheimy, N., Hou, H. Y., and Niu, X. J. (2008). Analysis and modeling of inertial sensors using Allan variance. IEEE Trans. Instrum. Meas. 57, 140–149. doi: 10.1109/TIM.2007.908635

Crossref Full Text | Google Scholar

Farrell, J., and de Haag, M. U. (2020). Position, navigation, and timing technologies in the 21st century. Hoboken, NJ: Wiley.

Google Scholar

Geng, J. J., Xia, L. Y., and Wu, D. J. (2021). Attitude and heading estimation for indoor positioning based on the adaptive cubature Kalman filter. Micromachines 12:79. doi: 10.3390/mi12010079

PubMed Abstract | Crossref Full Text | Google Scholar

Groves, D. P. (2008). Principles of GNSS, inertial, and multisensor integrated navigation systems. London: Artech House.

Google Scholar

Gui, H. C., and de Ruiter, H. J. A. (2018). Quaternion invariant extended Kalman filtering for spacecraft attitude estimation. J. Guid. Control. Dyn. 41, 863–878. doi: 10.2514/1.G003177

Crossref Full Text | Google Scholar

Huang, W. (2017). Quaternion constrained cubature Kalman filter attitude estimation based on uncertain measurements. J. Harb. Inst. Technol. 49, 116–121. doi: 10.11918/j.issn.201509022

Crossref Full Text | Google Scholar

Huang, Y., Wu, L. H., and Yu, Q. (2020). Underwater square-root cubature attitude estimator by use of quaternion-vector switch-ing and geomagnetic field tensor. J. Syst. Eng. Electron. 31, 804–814. doi: 10.23919/JSEE.2020.000055

Crossref Full Text | Google Scholar

Huang, W., Xie, H. S., Shen, C., and Li, J. P. (2016). A robust strong tracking cubature Kalman filter for spacecraft attitude estimation with quaternion constraint. Acta Astronaut. 121, 153–163. doi: 10.1016/j.actaastro.2016.01.009

Crossref Full Text | Google Scholar

Jiang, C., Zhang, S. B., Li, H., and Li, Z. (2021). Performance evaluation of the filters with adaptive factor and fading factor for GNSS/INS integrated systems. GPS Solutions 25:130. doi: 10.1007/s10291-021-01165-4

Crossref Full Text | Google Scholar

Julier, S., Uhlmann, J., and Durrant-Whyte, F. H. (2000). A new method for the nonlinear transformation of means and covariances in filters and estimators. IEEE Trans. Autom. Control. 45, 477–482. doi: 10.1109/9.847726

Crossref Full Text | Google Scholar

Li, B. F., and Chen, G. E. (2022). Improving the combined GNSS/INS positioning by using tightly integrated RTK. GPS Solutions 26:144. doi: 10.1007/s10291-022-01331-2

Crossref Full Text | Google Scholar

Liu, Y., Dong, K., Wang, H. P., Liu, J., He, Y., and Pan, L. N. (2014). Adaptive Gaussian sum squared-root cubature Kalman filter with Split-merge scheme for state estimation. Chin. J. Aeronaut. 27, 1242–1250. doi: 10.1016/j.cja.2014.09.007

Crossref Full Text | Google Scholar

Maebashi, K., Suematsu, N., and Hayashi, A. (2008). Component reduction for Gaussian mixture models. IEICE Trans. Inf. Syst. E91-D, 2846–2853. doi: 10.1093/ietisy/e91-d.12.2846

Crossref Full Text | Google Scholar

Markley, F. L. (2003). Attitude error representations for Kalman filtering. J. Guid. Control. Dyn. 26, 311–317. doi: 10.2514/2.5048

Crossref Full Text | Google Scholar

Michał, G., Michał, W., Cezary, S., Mariusz, K., Albert, Z., and Krystian, B. (2022). Quaternion attitude control system of highly maneuverable aircraft. Electronics 11, 1–13. doi: 10.3390/electronics11223775

Crossref Full Text | Google Scholar

Mohamed, A. H., and Schwarz, K. P. (1999). Adaptive Kalman Filtering for INS/GPS. J. Geod. 73, 193–203. doi: 10.1007/s001900050236

Crossref Full Text | Google Scholar

Noureldin, A., Karamat, B. T., and Georgy, J. (2013). Fundamentals of inertial navigation, satellite-based positioning and their integration. Berlin Heidelberg: Springer.

Google Scholar

Qian, C., Song, C. Y., Li, S., Chen, Q. W., and Guo, J. (2021). Algorithm of Gaussian sum filter based on SGQF for nonlinear non-Gaussian models. Int. J. Control. Autom. Syst. 19, 2830–2841. doi: 10.1007/s12555-020-0490-x

Crossref Full Text | Google Scholar

Qiu, Z. B., and Qian, H. M. (2018). Modified multiplicative quaternion cubature Kalman filter for attitude estimation. Int. J. Adapt. Control Sign. Process. 32, 1182–1190. doi: 10.1002/acs.2895

Crossref Full Text | Google Scholar

Ryzhkov, L. M. (2021). Quaternion attitude determination by vector measurement. Int. Appl. Mech. 57, 613–617. doi: 10.1007/s10778-021-01111-4

Crossref Full Text | Google Scholar

Savage, G. P. (1998). Strapdown inertial navigation integration algorithm design part 1: attitude algorithms. J. Guid. Contr. Dynam. 21, 19–28. doi: 10.2514/2.4228

Crossref Full Text | Google Scholar

Song, H. R., Cheng, S. X., Xu, Z. X., and Zang, N. Research on PPP/INS algorithm based on sequential Sage-Husa adaptive filtering. China Satellite Navigation Conference (CSNC 2022) Proceedings, Beijing, China, 26–28 April 2022; Springer Nature Singapore: Singapore, (2022).

Google Scholar

Sun, L. L., Cao, Y. H., Wu, W. H., and Liu, Y. T. (2020). A multi-target tracking algorithm based on Gaussian mixture model. J. Syst. Eng. Electron. 31, 482–487. doi: 10.23919/JSEE.2020.000020

Crossref Full Text | Google Scholar

Sun, B., Zhang, Z. W., Liu, S. C., Yan, X. B., and Yang, C. X. (2022). Integrated navigation algorithm based on multiple fading factors Kalman filter. Sensors 22:5081. doi: 10.3390/s22145081

PubMed Abstract | Crossref Full Text | Google Scholar

Swati, (2022). Continuous discrete cubature quadrature Kalman filter. Asian J. Control 24, 483–493. doi: 10.1002/asjc.2505

Crossref Full Text | Google Scholar

Taghizadeh, S., and Safabakhsh, R. (2023). Low-cost integrated INS/GNSS using adaptive H∞ cubature Kalman filter. J. Navig. 76, 1–19. doi: 10.1017/S0373463322000583

Crossref Full Text | Google Scholar

Tang, X. J., Liu, Z. B., and Zhang, J. S. (2012). Square-root quaternion cubature Kalman filtering for spacecraft attitude estimation. Acta Astronaut. 76, 84–94. doi: 10.1016/j.actaastro.2012.02.009

Crossref Full Text | Google Scholar

Teunissen, J. G. P., and Montenbruck, O. (2017). Springer handbook of global navigation satellite systems. Berlin Heidelberg: Springer.

Google Scholar

Wang, J. W., Chen, X. Y., and Shi, C. F. (2023). A novel robust iterated CKF for GNSS/SINS integrated navigation applications. Eurasip J. Adv. Sign. Process. 1:83. doi: 10.1186/s13634-023-01044-9

Crossref Full Text | Google Scholar

Wang, L., and Cheng, X. H. (2015). Algorithm of Gaussian sum filter based on high-order UKF for dynamic state estimation. Int. J. Contr. Autom. Syst. 13, 652–661. doi: 10.1007/s12555-014-0114-4

Crossref Full Text | Google Scholar

Wang, D. J., Dong, Y., Li, Q. S., Li, Z. Y., and Wu, J. (2018). Using Allan variance to improve stochastic modeling for accurate GNSS/INS integrated navigation. GPS Solutions 22:53. doi: 10.1007/s10291-018-0718-x

Crossref Full Text | Google Scholar

Wang, L., Gao, W.X., Wang, L., and Hu, F.Z. Design and analysis of Gaussian sum high-order CKF for nonlinear/non-Gaussian dynamic state estimation. 2021 33rd Chinese Control and Decision Conference (CCDC). IEEE: Kun-ming (2021).

Google Scholar

Wang, Z. P., Li, X., Zhu, Y. B., Li, Q., and Fang, K. (2022). Integrity monitoring of global navigation satellite system/inertial Navi-gation system integrated navigation system based on dynamic fading filter optimisation. IET Radar Sonar Navig. 16, 515–530. doi: 10.1049/rsn2.12199

Crossref Full Text | Google Scholar

Wang, D. J., Lv, H. F., and Wu, J. (2017). Augmented cubature Kalman filter for nonlinear RTK/MIMU integrated navigation with non-additive noise. Measurement 97, 111–125. doi: 10.1016/j.measurement.2016.10.056

Crossref Full Text | Google Scholar

Wang, G. C., Xu, X. S., and Zhang, T. (2020). M-M estimation-based robust cubature Kalman filter for INS/GPS integrated navigation system. IEEE Trans. Instrum. Meas. 70, 1–11. doi: 10.1109/TIM.2020.3021224

Crossref Full Text | Google Scholar

Wu, Y. L., Chen, S. A., and Yin, T. T. (2022). GNSS/INS tightly coupled navigation with robust adaptive extended Kalman filter. Int. J. Automot. Technol. 23, 1639–1649. doi: 10.1007/s12239-022-0142-7

Crossref Full Text | Google Scholar

Yang, Z. H., Li, Z. K., Liu, Z., Wang, C. C., Sun, Y. W., and Shao, K. F. (2021). Improved robust and adaptive filter based on non-holonomic constraints for RTK/INS integrated navigation. Meas. Sci. Technol. 32:105110. doi: 10.1088/1361-6501/ac0370

Crossref Full Text | Google Scholar

Yu, B., Zhang, Y. Z., Xie, W. S., Zuo, W. J., Zhao, Y. M., and Wei, Y. L. (2023). A network traffic anomaly detection method based on Gaussian mixture model. Electronics 12:1397. doi: 10.3390/electronics12061397

Crossref Full Text | Google Scholar

Zhang, L., Li, S., Zhang, E., Chen, Q. W., and Guo, J. (2019). Improved square root adaptive cubature Kalman filter. IET Signal Process. 13, 641–649. doi: 10.1049/iet-spr.2018.5029

Crossref Full Text | Google Scholar

Zhang, Q., Niu, X. J., and Shi, C. (2020). Impact assessment of various IMU error sources on the relative accuracy of the GNSS/INS systems. IEEE Sensors J. 20, 5026–5038. doi: 10.1109/JSEN.2020.2966379

Crossref Full Text | Google Scholar

Zhu, T. G., Liu, Y., Li, W. K., and Li, K. L. (2021). The quaternion-based attitude error for the nonlinear error model of the INS. IEEE Sensors J. 21, 25782–25795. doi: 10.1109/jsen.2021.3118039

Crossref Full Text | Google Scholar

Zhu, W. Q., McBrearty, W. I., Mousavi, S. M., Ellsworth, L. W., and Beroza, C. G. (2022). Earthquake phase association using a Bayesian Gaussian mixture model. J. Geophys. Res. Solid Earth. 5:127. doi: 10.1029/2021JB023249

Crossref Full Text | Google Scholar

Zickert, G., and Yarman, C. E. (2021). Gaussian mixture model decomposition of multivariate signals. Signal Image Video Process. 16, 429–436. doi: 10.1007/s11760-021-01961-y

Crossref Full Text | Google Scholar