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).

here, δ q is attitude quaternion errors, δ v is velocity errors, δ p is position errors, ε is the gyroscope biases drift, ∇ is the accelerometer biases drift. The state-space equation of GNSS/SINS-IADPS is described as follows:

here, x k | k − 1 is the state prediction, f ( ⋅ ) 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 ~ ( 0 , Q k ) .

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:

here, h ( ⋅ ) 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 ~ ( 0 , R k ) .

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 χ c , k − 1 | k − 1 is accomplished through the utilization of the third-order sphere-phase diameter cubature rule in formulas (4, 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 T , P k − 1 | k − 1 is the covariance pertains to the state estimation at epoch k − 1 , and c = 1 , 2 , ⋯ , m = 2 n , 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 ∗ are estimated through the nonlinear function f ( ⋅ ) , which can be expressed in formula (6).

Then, the state prediction x k | k − 1 at epoch k can be computed, as illustrated in formula (7).

The state prediction covariance P k | k − 1 can be derived, as depicted in formula (8).

Step 3: Update.

The cubature points  χ c , k | k − 1  can be re-estimated in formula (9).

here, P k | k − 1 = S k | k − 1 S k | k − 1 T .

The propagated cubature points z c , k | ∗   k − 1 can be estimated by applying the nonlinear function h ( ⋅ ) , which is expressed in formula (10).

Subsequently, the measurement prediction z k | k − 1 can be computed, as illustrated in formula (11).

The covariance P z z , k | k − 1 of the measurement prediction, cross-covariance P x z , k | k − 1 of the measurement prediction, and the filter gain K k can be derived, as depicted in formulas (12)–(14), respectively.

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).

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.

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:

here, ω k − 1 ( i ) denotes the weight of the i t h Gaussian component, N ( x k − 1 | k − 1 ; μ ¯ k − 1 ( i ) , P k − 1 | k − 1 ( i ) ) represents the i t h Gaussian component with a mean of μ ¯ k − 1 ( i ) and a variance of P k − 1 | k − 1 ( i ) , I signifies the total number of Gaussian components.

Correspondingly, the PDF of the process noise p ( w k ) in formula (2) and the PDF of the measurement noise p ( v k ) in formula (3) can be approximately expressed as follows:

here, α k ( j ) represents the weight of the j t h Gaussian component, N ( w k ; w ¯ k ( j ) , Q k ( j ) ) represents the j t h Gaussian component with a mean of w ¯ k ( j ) and a variance of Q k ( j ) , J represents the total number of Gaussian components; β k ( l ) represents the weight of the l t h Gaussian component, N ( v k ; v ¯ k ( l ) , R k ( l ) ) represents the l t h Gaussian component with a mean of v ¯ k ( l ) and a variance of R k ( l ) , L represents the total number of Gaussian components.

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 χ c , k − 1 | k − 1 are formulated utilizing the third-order sphere-phase diameter cubature rule.

here, S k − 1 | k − 1 ( r ) is derived through Cholesky decomposition of P k − 1 | k − 1 ( r ) , which can be expressed as P k − 1 | k − 1 ( r ) = S k − 1 | k − 1 ( r ) S k − 1 | k − 1 T ( r ) , r = ( i − 1 ) ⋅ I + j .

By propagating the cubature points χ c , k − 1 | k − 1 ( i ) through the nonlinear function f ( ⋅ ) , we can obtain χ c , k | k − 1 ∗ ( i ) = f ( χ c , k − 1 | k − 1 ( i ) ) . The state prediction x k | k − 1 ( r ) is calculated as follows:

Subsequently, the covariance of state prediction P k | k − 1 ( r ) can be derived, as depicted in formula (22).

here, ω c denotes the weight of the cubature points.

Step 2: Update.

The cubature points are assessed:

Here,

By propagating the cubature points χ c , k | k − 1 ( r ) through the nonlinear function h ( ⋅ ) , the propagated cubature points can be obtained as: z c , k | k − 1 ∗ ( r ) = h ( χ c , k | k − 1 ( r ) ) .  The measurement prediction z k | k − 1 ( r , l ) is calculated as follows:

The covariance matrix of measurement prediction P z z , k | k − 1 ( r , l ) , the cross-covariance P x z , k | k − 1 ( r , l ) , the filter gain K k ( r , l ) , the state estimation x k | k ( g ) and its corresponding covariance P k | k ( g ) can be derived, as depicted in formulas (25)–(29), respectively.

Here,

Step 3: Global point estimate.

The state estimation x k | k and its covariance matrix P k | k are computed as

here, ω k ( g ) denotes the weight of the g t h Gaussian component, which is calculated as formula (32).

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

From formulas (32) and (33), it is evident that the weight ω k ( g ) of the Gaussian component adaptively modifies in response to the innovation z k − z ( r , l ) 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 ⋅ J ⋅ L . If I ⋅ J ⋅ 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 , ⋯ , I ⋅ J ⋅ L , while retaining I − 1 components, then the weight value ω k ( I ) of the I − 1 Gaussian component at epoch k can be calculated:

The mean of state estimation x k | k ( I ) and its covariance P k | k ( I ) corresponding to the I t h Gaussian component is shown in formulas (35) and (36), respectively.

here, ω ˜ k ( g ) = ω k ( g ) / ω k ( I ) denotes the regularization weight. Through the aforementioned procedures, I − 1 Gaussian components and the I t h 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.

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 x → k | k and its covariance P → k | k can be computed, which are detailed below:

here, χ ˜ c , k | k = φ ( χ c , k | k ∗ ) , χ c , k | k ∗ = S ¯ k | k ξ c + x k | k , P k | k = S ¯ k | k S ¯ k | k T , φ ( ⋅ ) is the constraint function.

Step 2: the ultimate state estimation x ˜ k | k and its covariance P ˜ k | k can be calculated:

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

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).

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 m , and the sampling frequency of the GNSS receiver was 1 Hz.

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).

here, N signifies the total number of Monte Carlo simulations; X i represents the reference; X ¯ i signifies the estimation.

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.

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).

here, Δ λ is the positioning error in longitude; Δ L is the positioning error in latitude; Δ 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.

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.

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.

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.