An Improved Adaptive Unscented Kalman Filter With Application in the Deeply Integrated BDS/INS Navigation System

The integration of the BeiDou Navigation Satellite System(BDS) and the Inertial Navigation System(INS) can provide a more reliable and accurate navigation service than either system alone. The deeply coupled architecture for BDS and INS integration has more superior performance than the loosely coupled or the tightly coupled. Owing to the complicated dynamic scenario and the nonlinear system’s noise uncertainty, the adaptive Kalman filter(AKF) algorithm is often adopted in the deep integration(DI) system. The adaptive Sage window methods including innovation-based adaptive estimation(IAE) and residual-based adaptive estimation(RAE) are widely applied in AKF algorithms, but they have several limitations. We propose an improved adaptive unscented Kalman filter(AUKF) based on forgetting-factor-weight smoothing and multi-factor adaptation to overcome these limitations. Compared with the Sage window methods, the improved AUKF algorithm is immune to the quantity change of the satellites concerning integration and more sensitive to present dynamic. Furthermore, it can reduce the computation and storage burden in implementation. A simulation test based on a software platform and the deeply integrated BDS/INS navigation system is carried out to evaluate the performance of the improved AUKF. Simulation results show that the improved AUKF algorithm outperforms the extended Kalman filter(EKF) and has a similar performance with the RAE-AUKF.


I. INTRODUCTION
The BeiDou Navigation Satellite System(BDS) plays an important role in the Global Navigation Satellite System(GNSS), and it has made great progress in recent years. The BDS has been able to provide positioning and navigation service to global users since December 28, 2018. The signal of satellites is not subject to error accumulation by nature, but it is susceptible to interference, jamming, multipath and blocking obstacles. In contrast, the Inertial Navigation System(INS) is self-contained and inherently immune to exterior interference and environmental disturbances [1]. However, caused by biases and random noise, error accumulation of the INS easily leads the navigation solution to diverge with time. Therefore, the INS needs to be calibrated by the INS alignment and/or integration algorithms [2].
The associate editor coordinating the review of this manuscript and approving it for publication was Rajeeb Dey .
Considering both the advantages and drawbacks of GNSS and INS, integration of GNSS and INS can supply more accurate and reliable navigation information than either system alone. There are mainly three integrated architectures named loosely coupled, tightly coupled, and deeply coupled, respectively. The loosely coupled integration uses position and velocity estimations from GNSS and INS to form measurements. This architecture is sub-optimal in preventing outages (i.e., with less than four satellites in sight) [3], [4]. The tightly coupled architecture adopts pseudo range and doppler differences between GNSS and INS as the integration filter's input information. The performance of the tightly coupled architecture is improved [4] compared with the loosely coupled one, especially in degraded satellite constellations. In the deep integration(DI) architecture, the dynamic assistant information from INS is fed back to tracking loops. Substantial improvements in the performance of noise suppression can be achieved in the deeply coupled architecture, compared to VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a traditional loop filter [5]. Generally, the deeply coupled integration can be classified into two architectures named centralized filtering architecture and federated filtering architecture [6], [7]. The centralized filtering architecture combines the in-phase and quadrature correlator components as measurements [8]- [10], and the federated one utilizes outputs from discriminators or pre-filters as measurements [11]- [13].
The system architecture proposed in this paper is a federated one based on [13]. The deep integration of BDS and INS introduces nonlinearity into the system model so that the traditional linear Kalman filter(KF) is not suitable for the integration filtering. The Extended Kalman Filter(EKF) is one of the most popular methods for sensor integration. The EKF was first applied in sensor fusion by Willner et al. [14], and it has become one of the most popular filtering methods because of its efficiency and accuracy. In the EKF algorithm, the first-order Taylor series expansion is adopted. In this expansion, nonlinear models' Jacobian matrices are calculated, and the high order terms are abandoned. Caused by the first-order linearization process, the EKF easily suffers from performance degradation, especially in strong nonlinear systems. On account of high maneuverability in high-speed applications, it is difficult to neglect the high-order nonlinear items. The propagation process of the system error should be modeled by the nonlinear system model. As a result, the filtering method for strong nonlinear systems is required for the integration of GNSS and INS [15].
The Unscented Kalman Filter which was first proposed by Julier et al. [16] is regarded as a better alternative for the integration of BDS and INS. The UKF introduces a minimal set of sigma points which are chosen deterministically to approximate the probability distribution of the nonlinear model by ''unscented transformation''(UT) instead of relying on analytical linearization of the system. Unlike the linearization process of the EKF, the UT does not need to calculate the partial derivatives. Additionally, the posterior mean and covariance of the Gaussian random variable can be approximated with second-order accuracy [17]. The UKF has higher estimation accuracy and convergence than EKF [18], but the UKF filtering solution will be unstable if stochastic errors are involved in the pre-defined kinematic and measurement model [19].
The uncertainty of the process noise and measurement noise reduces the system positioning performance and dynamic adaptability. If the statistic properties of the system are set unchanged, it may result in performance degradation or even divergence. One of the effective approaches is applying the adaptive Kalman filters(AKFs). Sage and Husa proposed adaptive Bayes estimation algorithms with unknown prior statistics for filtering which are called the innovation-based adaptive estimation(IAE) and residual-based adaptive estimation(RAE) [20]. The Sage window methods reevaluate the noise covariance matrixes of the predicted state and the measurement vector, but it requires the number of satellites tracked to stay constant.
Yang et al. developed a flexible AKF method based on the adaptive factor calculated by utilizing innovation vectors or difference between predicted and estimated states [21], [22]. In [23], an adaptive federated filter method is applied to the PPP/INS integrated system to improve filter efficiency and adaptivity. The experiment results show that the performance of the filter is improved. It should be noted that all the adaptive Kalman filtering methods above are related to the linear state-space model.
Many adaptation methods related to the UKF have been put forward by researchers in recent years. Gao et al. [24] proposed an adaptive UKF(AUKF) by combining the windowing and random weighting concepts. The proposed adaptive UKF outperforms the conventional one under the condition without precise knowledge of system noise statistics. In [25], an adaptation strategy is provided to develop a modified Masreliez-Martin UKF for discrete-time nonlinear stochastic systems. The proposed filter is easy to be implemented and derivative-free. Gao et al. [26] presented a novel adaptive UKF by employing the maximum likelihood principle(MLP) and moving horizon estimation(MHE). The proposed adaptive UKF can achieve the online estimation of system noise and improve robustness. Cho and Choi [27] reported a sigma-point Kalman filter based on receding-horizon strategy. The proposed filter can work even in the case of exiting the unmodeled random walk of the inertial sensors. Various situations indicate that the proposed filter can be implemented in a low-cost DR/GPS integration system for seamless land navigation.
In this paper, we propose an improved AUKF algorithm to overcome the limitations of the Sage window methods. The improved AUKF algorithm has similar performance with the Sage window methods, but it is immune to the quantity change of the visible satellites and more sensitive to dynamic change compared with the Sage window method. Additionally, it can reduce the computation and storage burden. This algorithm is also applied in a deeply integrated BDS/INS system in which the performance of the proposed AUKF algorithm is evaluated. The following sections are organized as follows. In section II, we first introduce the deeply integrated BDS/INS model including dynamic model and measurement model. Then in section III, we interpret the basic procedures of the standard UKF, and an improved AUKF algorithm is proposed. In section IV, we carry out a simulation test based on a software receiver to compare the performance of the proposed AUKF with the other two types of Kalman filters. Simulation results and discussions are presented. Finally, a conclusion comes as the last section.  carrier is mixed with the IF signal to obtain in-phase and quadrature-phase signals. The correlator then correlates the processed signals with the one-half chip early, prompt,and one-half chip late replica codes in each channel. Six I-Q integration values for each channel can be obtained in every integration interval, and they are further processed by prefilters. The pre-filters can estimate tracking errors and their covariances to prepare for the integration filtering. Next, the integration filter based on the adaptive UKF estimates an error-state vector that can periodically correct the INS outputs, biases of the Inertial Measurement Unit(IMU), and receiver clock status. Specific force and angular rate measurements from the IMU are delivered to the Inertial Navigation Processor, and high frequency position and velocity estimations are calculated. Finally, utilizing ephemerides, the receiver's position and velocity estimations are projected along the line of sight of each satellite. Numerically Controlled Oscillator(NCO) commands are generated with the assistance of dynamic information from the INS. The loop of the whole system is closed. Because all the channels share dynamic information, which is called the vector tracking architecture, the tracking and positioning performance of the deeply integrated system is improved.

B. INTEGRATION FILTER MATHEMATICAL MODEL
The principle of the deeply coupled BDS/INS integration is to employ pre-filters' outputs as measurements to estimate the error of the INS and the receiver clock status. Details about the pre-filter model can be found in [6], [13]. In this section, the mathematical model consisting of a dynamic model and a measurement model is established for the deep integration system. The dynamic model is derived from the inertial navigation equations and the IMU error equations. The measurement model is constructed by utilizing the pre-filters' output information.

1) SYSTEM DYNAMIC MODEL
The state vector of the integration filter is a 17-state one which includes INS errors and receiver clock errors. The components of the state vector X k are defined as follows, and they are explained in Table 1.
where the superscript ''T '' means matrix transposition, and the subscript ''k'' means at the epoch of k. δϕ = [δϕ E δϕ N δϕ U ] T represents the vector of pitch, roll, and yaw misalignment angles. The term δv n = [δv E δv N δv U ] T is the velocity error vector expressed in n-frame. The position error vector is expressed as δp = [δL δλ δh] T where L, λ, and h represent latitude, longitude, and height respectively. In this paper, the body frame(b-frame) is selected as ''Right-Front-Up'', and the local navigation frame(n-frame) is defined as ''East-North-Up(ENU)''. The n-frame means the actual navigation frame, and n is the navigation frame estimated by the receiver. i-frame represents the geocentric inertial frame, e-frame is the Earth-Centered Earth-Fixed(ECEF) frame.
The classic INS error propagation model is based on the assumption that the misalignment angles are small. Nevertheless, this assumption may be broken by environmental disturbances and sensor errors [28].
Stochastic nonlinearities from high-maneuverability or sudden environment changes need to be considered. A nonlinear dynamic model is adopted to avoid small angle assumption and improve the robustness of the filter [29]. The nonlinear system dynamic model of INS is written as where C n n is the rotation matrix from n-frame to n -frame. Similarly, C n n and C n b are also rotation matrixes. ω n in is the rotation velocity of the n-frame with respect to i-frame solved in n-frame. δω n in is the error vector of the angular velocity ω n in , and δω b ib is the error vector of ω b ib which is measured by the gyro. δω b ib is composed of gyro constant drift error ε b and Gaussian white noise. ω n ie is the rotation velocity vector of e-frame with respect to i-frame solved in n-frame, and δω n ie is the error vector of it. ω n en stands for the rotation velocity vector of n-frame relative to e-frame, and δω n en is the corresponding error. f b ib is the specific force sensed by the accelerometer resolved about the b-frame, and δf b ib is the error vector of f b ib . δf b ib consists of the accelerometer bias error ∇ b and Gaussian white noise. R M is the median radius, and R N is the normal radius. C −1 ω is the transformation matrix from the relative angular velocity to the Euler angle error. When the small angle assumption is valid, we can get C ω ≈ I and I − C n n ≈ (δϕ×), and the linear system dynamic model can be established.
The gyro bias error and accelerometer bias error are set as random constants, which is expressed as The receiver clock model is shown as follows where w b and w d are Gaussian white noise for clock bias and drift. δd clk is modeled as first-order Markov process, and τ clk is the time constant. Combining (2) to (4), the discrete-time form after discretization of the nonlinear dynamic model of the deeply integrated BDS/INS system can be summarized as follows where f (·) represents the nonlinear function and W k−1 is the process noise vector.

2) MEASUREMENT MODEL
Pre-filters can provide tracking error estimations and their covariances. The discrete form of the measurement model is given by (6) where Z k is the measurement vector, H k is the observation matrix, and V k is the measurement noise vector.
In the deeply coupled integration, the NCO commands are generated using INS status information. Therefore, the tracking errors of tracking channels contain the information of the INS's residual errors. Code phase error and carrier frequency error estimations from pre-filters are taken as measurements of the integration filter. The measurement vector Z k at epoch k is defined as where N is the tracking channel number, and i = 1, 2, . . . , N . i means the ith channel of all N tracking channels. The pseudo range error vector δρ and the doppler error vector δρ are the converting results of (8) and (9). In (8) and (9), c is the speed of light, δτ i is the code phase error estimation from the ith channel, and δω i is the carrier frequency error estimation. f code and f carr represent code frequency and carrier frequency, respectively. The noise covariance matrix R k is directly related to the pre-filters' state covariance matrixes [13]. The measurement matrix H k is given in (10).
and H ρ1 , H ρ2 , Hρ1 and Hρ2 are defined as follows where O M×N is a zero matrix of size M × N , and similarly is the unit vector of line-of-sight direction from user to the ith satellite. C 1 and C 2 are transformation matrixes between different frames.

III. IMPROVED ADAPTIVE UNSCENTED KALMAN FILTER A. STANDARD UKF
Considering the following nonlinear system model with additive noises (13) where X k and Z k are the state vector and measurement vector at epoch k. f (·) is the nonlinear function of the process model, and H k is the measurement transposition matrix. W k−1 and V k are the process and measurement noises. They are assumed as additive zero-mean Gaussian white noise with covariances Q k and R k , which is shown as follows where δ kf is the Kronecker-δ function.
The UKF uses unscented transform technique to approximate the probability distribution of the nonlinear system state. A set of weighted sigma points whose sample mean and covariance match those of the priori distribution is deterministically chosen. Then the nonlinear system function is applied to the sigma points to yield the transformed samples. By using the weighted mean and covariance of the transformed samples, the predicted mean and covariance are calculated [30]. The calculation steps of the UKF are summarized as follows [31], [32] step 1: Initialization of the state vector and covariance matrix.
step 2: Calculation of the sigma points. In this step, the symmetric sampling strategy based on the Cholesky decomposition is adopted [33]. Assuming the state estimateX k−1 and the error covariance matrixP k−1 at epoch k − 1 are given, the sigma points can be selected as where (P k−1 ) 1 2 i denotes the ith column of the square root of matrixP k−1 . n is the element number of the state vector. γ is a scale factor which is defined as where λ = α 2 (n + κ) − n is a scale factor, and κ is a secondary parameter which is usually equal to or greater than zero [31]. The parameter α determines the spread of the sigma points aroundX k−1 , and it is usually a positive value smaller than 1(e.g., 0.0001 ≤ α ≤ 1).
where W n i and W c i represent the weight of the mean and covariance, respectively. β incorporates prior knowledge of the distribution of X, and β = 2 when X is normally distributed [34].
The sigma points are instantiated through the dynamic model f (·) to get the set of the transformed samples The predicted mean X − k and covariance P − k are updated as   where the superscript ''−'' means the prior. step 4: Measurement Update. Since the measurement model is a linear one, the measurement update process can be performed as the standard Kalman Filter. The measurement update process at epoch k is shown as followŝ where the superscript ''+'' represents the posteriori. K k is the Kalman filter gain at epoch k.

B. IMPROVED ADAPTIVE UKF
The Sage window methods [20], [35] including innovation-based adaptive estimation(IAE) and residual-based adaptive estimation(RAE) estimate the covariance matrixes by a moving window. Because the noise matrixes are reevaluated and adjusted according to the real data, the Sage window method improves the performance of the Kalman filter and it is widely applied. However, the Sage window method has several limitations itemized as follows [36], [37] • Since the Sage window method evaluates the covariance matrix by a moving window, the response of the adaptive evaluation is delayed. The historical innovation or residual sequences contribute equally to the evaluation, which makes the adaptive mechanism insensitive to the dynamic change.
• The smoothing method implies that the number of the measurements should have the same dimension. VOLUME 8, 2020 Once the number of measurements changes(e.g. some tracking channels lose lock due to building blocks), the method will fail to work. This essential problem may lead the IAE and RAE unusable in the dynamic condition.
• To estimate the innovation or residual sequence covariance matrix, we have to use the innovation or residual sequences and covariance matrixes of m epoches where m is the window width. It causes a significant storage and computation burden especially when the dimension of the matrix is high. An improved AUKF algorithm is proposed in this section to overcome these limitations. This adaptive filtering algorithm contains two steps including Q-Adaptation and R-Adaptation. The adaptation process of the system noise matrix Q k and the measurement noise matrix R k is presented next.

The innovation vector v k at epoch k is expressed as
The estimation of the system noise matrix only using the innovation vector at the current epoch is written as follows. To guaranteeQ cur,k to be positive definite, an approximation is adopted here [38]. (27) where the subscript ''cur'' means ''current'', and the innovation covariance matrixĈ cur,v − k at epoch k is estimated using the current innovation vector v − k as followŝ In order to eliminate the noise correlation and reduce computation burden,Ĉ cur,v − k is supposed to be a diagonal matrix as followsĈ where ''•'' represents the Hadamard product operation, and diag{a} means a diagonal matrix whose diagonal elements consist of the vector a. However, the estimationQ cur,k is noisy, becauseQ cur,k is estimated only using the innovation sequence at current epoch without smoothing. Hence,Q cur,k is not suitable for application in AUKF. A smoothing window method whose window width is m based on forgetting factor weights is implemented here. The smoothing estimation of the system noise matrix at epoch k + 1 is given where µ j is the jth element of the weight vector. µ j can be prepared before the filter begins to work, and it is defined as 1, 2, . . . , m) (31) Generally, the constant d is satisfied with 0.90 < d < 0.99, and d = 0.95 is adopted here. The forgetting factor weights are plotted in Figure 2. As can be seen from the figure, the weight decreases with the increase of index j, which means that the termQ cur,k−j+1 closer to epoch k contributes more to the estimationQ k+1 . Utilizing this forgetting-factor-weight-based smoothing method, the sensitivity to dynamic change is improved compared with the average smoothing method. Additionally, instead of expensive storage and computation consumption in the smoothing process of matrixĈ cur,v − k , the object of the smoothing process is transformed into the diagonal elements of Q. Hence, the computation and storage burden is reduced in theQ k estimation process compared with the IAE and RAE method. Furthermore, the window smoothing process of the matrixĈ v − k in IAE and RAE requires that the number of observations remain the same in the smoothing window. The proposed AUKF algorithm avoids this limitation by smoothing the diagonal elements of Q whose dimension is constant.

2) R-ADAPTATION
The residual vector r k at epoch k is defined as whereẐ + k is the posteriori estimation of the measurement vector Z k .
The estimation of the residual covariance matrixĈ v + k is derived from [35] as follows, and an adaptive scale factor matrix S k is adopted here for the adaptation of the measurement noise covariance matrix R whose dimension is (2N × 2N ) where R − ρ,k and R − ρ,k are related to pseudo range error and pseudo range rate error, respectively. R − ρ,k and R − ρ,k are obtained from the pre-filters' outputs before the integration filter update [13].
Since R ρ and R P ρ have different units, it is not appropriate to adapt the matrix R as a whole. For precise filtering, the components of the scale factor S k for R ρ and R P ρ should be estimated separately. S k is defined as (35) where s 1,k and s 2,k are the adaptive coefficients at epoch k, and the posteriori estimation of R at epoch k is given To estimate the adaptive coefficients s 1,k and s 2,k of S k , calculate the trace of both sides of (33) and we can get where tr(·) a:b represents the sum of the ath diagonal element to the bth diagonal element of the related matrix. Finally, the estimation of adaptive coefficients for R-adaptation can be calculated from (37) as follows Because the R-adaptation process is executed after UKF measurement update,ŝ 1,k+1 andŝ 2,k+1 estimated at epoch k are prepared for epoch (k + 1). When epoch k = 1, bothŝ 1,k and s 2,k are equal to 1 for initialization. On the purpose of overcoming the limitation 2(i.e. the measurement number should remain the same) mentioned above, the following smoothing method for the trace estimations of residual covariance matrix tr(Ĉ v + where N j is the channel number at epoch j, the operation (·) a:b stands for the ath element to the bth element of the related vector. µ j is the weight factor defined in (31). Using the smoothing method in (39), the estimation process is not subject to the change of the measurement number which is related to channel number N k . Furthermore, similarly with the Q-adaptation process, the later residual sequence contributes more than the earlier one does to the S k estimation because of the property of µ j . The storage and computation burden is also reduced compared with the Sage window method, because only the traces of the residual covariance matrixesĈ v + k in the smoothing window are computed and stored.
The adaptation process of the proposed algorithm can be summarized as Table 2.

IV. PERFORMANCE EVALUATION AND DISCUSSIONS A. SIMULATION SETUP
In this section, a simulation is carried out to evaluate the performance of the AUKF proposed in comparison with EKF and RAE-based AUKF. The EKF is the most widely applied in the nonlinear navigation system, and the RAE-based adaptive method is one of the Sage window methods. These three algorithms are employed in a deeply integrated BDS/INS system. In order to get better insights into the ability of the proposed AUKF, the simulation test is implemented on a Matlab-based software receiver on a PC, in which way the navigation solution can be compared with the designed trajectory, and the navigation error can be calculated precisely. Since the dynamic information is fed back to the tracking loops in the deep integration system, which is the main distinction between the tightly coupled and the deeply coupled architecture, the simulation is designed at the Intermediate Frequency(IF) signal level rather than depending on pseudo range and doppler measurements.
Before showing the simulation results, a summary of assumptions is made as follows.
• IF signals are generated with standard atmospheric models. Satellite clock errors, multipath errors, and ephemeris errors are not considered.
• Only scale factor errors, constant bias, and random walk noise are considered for INS simulation. Askew installation errors, lever arm errors, and correlated bias errors are not modeled in the simulation.
• The INS has been calibrated properly before the integration filter begins to work. The simulation framework is shown in Figure 3. This figure mainly illustrates the generation procedure of IMU measurements and IF signals. Each block stands for a module with a explicit function. Different colors represent different types of functions. The relationships among the modules are also described in this figure.
The user dynamic is designed as shown in Figure 4. Figure 4(a), 4(b), and 4(c) indicate the trajectory, velocity, VOLUME 8, 2020  and attitude reference, respectively. Table 3 illustrates the receiver initial status. The motion states of the receiver contain acceleration, climbing, turning, circling, and uniform moving.
In this simulation, raw IMU measurements without noise are firstly generated with the trajectory, velocity, and attitude references mentioned above. Then artificial measurement noise is added to the raw measurements. The noise model consists of scale factor errors, constant bias, and random walk noise. The parameters of both accelerometer and gyroscope of IMU are in accordance with the MEMS grade IMU [39] and listed in Table 4.
Satellite position and velocity are calculated based on actual BDS satellite ephemerides achieved from International GNSS Service(IGS) products. The broadcast ephemerides are recorded on 29th, January 2020 from BCEmerge. The ephemerides are also utilized for positioning and integration of the DI system. We select seven BDS visible satellites for this simulation. The aeroview of these satellites at the beginning epoch are expressed in Figure 5. The elevation angles and azimuth angles are marked in turn in this figure. The PRN numbers of each satellite are labeled on the right side of the corresponding markers.
With the user's and the satellites' dynamic information, the distance and doppler can be calculated. Then atmospheric delay computed by standard atmospheric model are added to obtain carrier phase and code phase in the simulation. We adopts the BDS B3I signal format in the IF signal generation process, and the related parameters are listed in Table 5. White gaussian noise is injected in the digital IF data based on a prepared C/N 0 value for all satellites. Finally, the digital IF data is quantized and recorded in a text file to process later.
The DI BDS/INS system based on a software receiver processes the IMU measurements and IF signals using  EKF, RAE-AUKF, and the proposed AUKF, respectively. Simulation results and discussions are presented in the following section.

B. RESULTS AND DISCUSSIONS
The prepared IF signals and IMU measurements are transmitted into the deeply integrated BDS/INS system as shown in Figure 3. We apply the Kalman filter methods including   EKF, RAE-AUKF, and Improved AUKF for signal processing and integrating. The navigation solutions consisting of attitude, velocity, and position estimations for each method are then recorded.
The navigation solution errors computed by comparing the estimations with the references of trajectory, velocity, and attitude are shown in Figure 4. The attitude errors, velocity errors, and position errors are illustrated in Figure 6, Figure 7, Figure 8, respectively. The IMU processing frequency is 200Hz, and the system solution output sample rate is 5Hz. To verify that the improved AUKF is immune to the number change of the satellites concerning integration, we control the changing trend of the visible satellite number as shown in Figure 9 when the improved AUKF is implemented.
The Figure 6 presents the attitude estimation errors of the three Kalman filter algorithms. This figure reveals that the attitude errors fluctuate with the changing trend of user dynamics. As can be seen from this figure, the noisy parts(t = 20∼28s and 36∼44s) correspond to the two circling motion states expressed in Figure 4. Similarly, changing of the user motion state (e.g. accelerating, climbing, and turning) results in increase of attitude error more or less. Compared with EKF, the adaptation strategy of RAE-AUKF and the improved AUKF helps to alleviate the attitude estimation error influenced by the changes of user's movement direction. Additionally, the Yaw estimation errors of the improved  AUKF are closer to zero compared with RAE-AUKF, which indicates that the improved AUKF is more sensitive to the change of heading direction. Moreover, it should be noted that the mean values of the attitude estimation errors are not equal to zero as expected. The main reason is that the grade of the IMU in the simulation is selected as MEMS whose large  scale factor errors and constant bias have a strong influence on accuracy of the IMU measurements. Figure 7 illustrates the velocity estimation errors in the ENU navigation frame, and Figure 8 shows the position estimation errors in the ECEF frame. From these two figures, we can draw the conclusion that the velocity and position estimation errors of RAE-AUKF and improved AUKF are on the same level. The EKF has larger estimation errors than the RAE-AUKF and the improved AUKF.
The overall estimation of navigation solution error at epoch k is defined as where the subscripts ''1, 2, 3'' represent the three orthogonal components' indexes of the error vector x k . x k represents the attitude error, velocity error, or position error here. The RMSE of the navigation solution errors is expressed as where M is the total number of epochs.
The overall estimation errors of the navigation solutions are compared next. RMSEs and maximum absolute(Max Abs) errors are computed and recorded in Table 6. From the error statistics of the navigation solutions, we can conclude that the improved-AUKF outperforms the EKF and has similar performance with the RAE-AUKF in the accuracy of attitude, velocity, and position estimations. There is no adaptation process in the EKF, that is why it has the poorest performance among the three algorithms. Although the RAE-AUKF and the improved AUKF implement the adaptation process in different manners, the two AUKF methods both extract noise covariance information from the residual sequences of the UKF, and that is the reason why they have similar performance.
The total-time-cost statistics of adaptation processes are also collected and recorded for RAE-AUKF and the improved AUKF, respectively. With the same processor and simulation conditions, the time cost is 0.2935s for RAE-AUKF and 0.1851s for the improved AUKF. The computation burden is reduced when the improved AUKF is implemented.

V. CONCLUSION
The conventional Sage window method for Kalman filter adaptation has limitations in application. We present an improved AUKF algorithm to overcome these limitations with application in the deeply integrated BDS/INS system. A forgetting-factor-weight smoothing method and a multi-factor adaptation method are adopted in the proposed AUKF algorithm. Without loss of smoothness, the forgettingfactor-weight smoothing method is more sensitive to the present motion compared with the average one. The smoothing failure problem caused by number changing of tracking channels is avoided. Additionally, using the improved AUKF algorithm, the storage and computation burden during the adaptation process is reduced. A simulation test based on a software receiver is carried out to evaluate the performance of the proposed algorithm. The results show that the proposed AUKF algorithm outperforms EKF and has similar performance with RAE-AUKF in terms of attitude, velocity and position accuracy. In future work, we will concentrate on its real-time implementation on a hardware platform.