Single-Frequency GPS/BDS RTK and INS Ambiguity Resolution and Positioning Performance Enhanced with Positional Polynomial Fitting Constraint

Single-frequency GPS/BeiDou navigation satellite system (BDS) real-time kinematic (RTK) and inertial navigation system (INS) integration has wide range of application prospects due to the global deployment of GPS along with the rapid development of BDS. The instantaneous single-frequency ambiguity resolution will be significantly improved by the combined GPS/BDS and INS configuration. Owing to road conditions and an inertial measurement unit (IMU) on the carrier not being rigidly mounted, biased measurements in the IMU will occasionally emerge, leading to biased INS predictions. However, bias or inaccuracy from INS-predicted position can prevent the successful resolution of the whole set of ambiguities. This paper proposes the use of a positional polynomial fitting (PPF) constraint to compensate for the epochs with abnormal INS predictions. The aid from PPF is provided at two levels, i.e., at the ambiguity resolution (AR) level and at the solution level. In order to further increase the availability of ambiguity-fixed positioning solutions, a partial ambiguity resolution (PAR) strategy is introduced when full ambiguity resolution (FAR) fails. A field vehicular experiment was performed to show the validity of the proposed PPF-aided method by comparing different schemes regarding different INS-aided satellite system configurations, different AR strategies, and whether the PPF-aided method was adopted. The results show that the most attractive scheme is to combine the PAR with the PPF-aided dual-constellation and INS integration.


Introduction
Global navigation satellite systems (GNSS) technology has been widely researched and applied to many engineering applications, such as vehicle navigation, personal positioning, and mobile surveying. For centimeter-level demanding applications, the high-precision carrier phase observations-based real-time kinematic (RTK) positioning is attractive, which relies on a set of correctly fixed integer ambiguities [1]. In open-sky and short baseline conditions, a dual-frequency global positioning system (GPS) RTK technique can achieve rapid, even instantaneous, ambiguity resolution (AR) with a high success rate [2]. However, when it comes to single-frequency RTK for low-cost GPS receivers, particularly for single-epoch kinematic positioning, there is a low success rate for AR. Along with the rapid development of the Chinese BeiDou navigation satellite system (BDS), combining the observations of GPS and BDS improves the single-frequency single-epoch AR performance significantly [3].
The GNSS-only configuration may not provide continuous and reliable solutions in adverse circumstances due to the frequent signal blockages and multipath effects, even with the combined GPS and BDS constellations. Attracted by the reduced size and low-cost inertial navigation system (INS) using micro-electro-mechanical sensors (MEMS) technology, integrating the low-cost inertial sensors with a single-frequency receiver has broad application prospects in the civil market [4]. For GNSS RTK/INS integration is required reliable ambiguity resolution: incorrect ambiguities can cause position errors as well as sensor biases divergent and failed ambiguity resolution reduces availability.
In order to aid the AR with inertial sensors, the tightly coupled Kalman filter (TCKF) mechanization is commonly adopted, which is carried out in the raw measurements domain. To handle the ambiguity problem in GNSS/INS integration, two approaches can be used, i.e., using carrier-smoothing-code filtering (CSCF) and performing integer AR. These two techniques are also usually chosen by GNSS-only users [5]. For carrier-smoothing-code filtering, the epoch-differenced carrier phases and the code pseudoranges are simultaneously considered [6]. Although the ambiguities are canceled by epoch difference operations, the limitations of CSCF are that cycle slips should be absent or repaired in advance and the filter state estimates are essentially the delayed ones when considering the correlation between process and measurement noises [7]. For the integer AR technique, the AR process in the TCKF is usually resolved in two steps. The first step is to solve the real-valued ambiguities and its variance-covariance matrix through least squares or the Kalman filter. Then, the second step tries to fix the real-valued ambiguities with a chosen integer least-squares estimate. For the single-epoch AR mode, if the ambiguities can be correctly fixed during the test period, the centimeter-level and real-time positioning results will be maintained [8]. In this paper, the integer AR technique is taken care of. Biases in satellite observations or inertial measurement unit (IMU) measurements will inevitably hamper the AR performance or result in incorrectly fixed ambiguities. For the observation-related biases, which may come from residual atmospheric error, large multipath error, and abnormal observation noise, the corresponding methods to relieve or delete the bias effects in the AR process have been extensively researched and can be categorized as the following three kinds: (i) the innovation-based reweighting method, also called the robust estimation method, which uses an innovation-based inflating factor to downscale the weight of the observation matrix [9][10][11]. In this way, the effect of biased observations on float ambiguity will be reduced so as to improve the precision of the float ambiguity solution. (ii) There is the partial ambiguity resolution (PAR) strategy, which uses a predetermined subset ambiguities chosen strategy to fix an ambiguity subset with a high success rate. The subset ambiguities could be chosen based on the successively increased elevations [12], the bias-affected success rate [13], the maximum rounding fixing rate ordering method [4], and so on. The PAR strategy is beneficial to the availability of positioning results, even with the degraded position accuracy due to fewer observations. This situation will have less of an effect especially when using multiple GNSS constellations with many satellites visible [14,15]. (iii) There is the receiver autonomous integrity monitoring (RAIM) method, which can directly delete the biased observations by outlier identification, reliability testing, and isolation steps [16]. If all biased observations are identified and isolated, pure observations can be used for subsequent steps.
For biased measurements (outliers or gross errors) from inertial sensors, the INS predictions from integral operation in INS mechanization will be affected. Outliers in IMU measurements occasionally occur when the carrier of an IMU is traveling along a bumpy road or passing through a deceleration strip in an urban environment. The sudden bump will affect the outputs of the accelerometer triple and may affect the outputs of the gyroscope triple due to the fact that an IMU is often not firmly fixed on a carrier. Therefore, how to alleviate the negative effects of abnormal predictions on the AR and the final filter solutions are our main focus in this work.
With biased IMU measurements, abnormal predicted residuals (innovations) will be found during the AR process in the TCKF. In this case, the adaptive estimation method could be used to give less weight to the variance-covariance matrix of the predicted states and reduce the negative effect of the predicted position in estimating the float ambiguities [17,18]. The RAIM method could also be adopted to resist the abnormal dynamic predictions as shown in [19], which use an extended RAIM (eRAIM) procedure to detect faults in predictions. The eRAIM procedures are executed within the framework of the Gauss-Markov model, where the observation vector is combined by the actual observations and the predicted system states. Then, least-squares estimation is adopted and the abnormal predictions can be judged according to residual testing. The above observation vector is similar to the combined observation vector used in the AR process of the TCKF, where the virtual observations that are the components of the combined observation vector consist of the predicted position components. However, both the adaptive estimation and RAIM methods could be marginally helpful for successful AR. This is because either down-weighting the variance-covariance matrix of states or deleting the corresponding abnormal state components means the reduction of redundant information in estimating the reliable float ambiguity estimates. Even if the ambiguities are fixed by using the above two methods, biased predictions will also influence the final TCKF states due to the fact that the final states are obtained by subtracting the estimated state errors from the approximates (i.e., the biased predicted states). Since subdecimeter-level predicted position is needed for successful AR in GNSS/INS integration, one should maintain the prediction accuracy and minimize the impact of biased predictions at the same time. From this point of view, a multiple model adaptive estimation method could be considered, which originates from tracking applications [20] and has been extended to multi-sensor integration [21][22][23]. The multiple model approach shares the same observation model but assumes that the system obeys one of a finite number of dynamic models. The final solution of multiple model adaptive estimation is the fusion of all Kalman filter solutions with respect to each dynamic model and with the corresponding transition probabilities as weights. It is unquestionable that the computational expenses will increase linearly with the number of dynamic models. Moreover, the motion behavior of one epoch should be uniquely specified by a reasonable dynamic model rather than a fused model with empirical weights.
To ensure the optimization of the Kalman filter when the system experiences unexpected dynamic conditions, a dynamic model that is able to adapt to a wide range of dynamic conditions should be guaranteed [24]. In [25], a positional polynomial fitting (PPF) approach with model orders from 1 to 4 was used to construct the dynamic models. Furthermore, in [26], a model evaluation criterion was derived to choose a reliable dynamic model from a series of candidate models with different model orders in the given time window. As has been pointed out by [25,26], an object's motion could be characterized by a low-order polynomial constructed by the positions of multiple successive historical epochs in a given time window. Therefore, in [27], the PPF approach with a model order of 2 was used for single-frequency GNSS cycle slip estimation, and the estimation of cycle slips can be significantly improved compared with that from two traditional methods, i.e., measurement-based polynomial fitting and triple-differenced (TD) residual-based snooping. A main advantage of the constructed positional polynomial is that it does not depend on external measurements; it requires only that the positions within the time window are all available. However, the methods in [25][26][27] were only proposed to enhance the performance of the GNSS. For GNSS and INS integration, the INS predictions may be affected by gross errors existing in the IMU measurements. Thus, using PPF as an auxiliary model to aid the inertial dynamic model may improve the performance of GNSS/INS integration due to their complementary characteristics. Hence, when abnormal INS predictions are detected, one can use the position predicted by the constructed polynomial fitting model as an alternative. In this paper, we propose a PPF-aided GPS/BDS RTK and INS integration method to estimate the navigation solution in order to achieve continuous and reliable solutions for both ambiguities and positioning performance. We also propose the use of the model probabilities for both the INS and PPF models together with the ambiguity fixing state to detect an epoch with abnormal predictions. The PPF will be triggered once an epoch with abnormal predictions is detected. Furthermore, the full ambiguity resolution (FAR) and PAR strategies are respectively embedded into the proposed method while fixing the float ambiguities.
The typical characteristics of the proposed method when comparing with some existing methods in the domain of GNSS and INS integration are threefold: (i) with the aid of the predicted position from PPF, the proposed method has the ability to resist the abnormal predictions from INS. These abnormal predictions cannot be handled by the commonly used GNSS/INS integrated methods, e.g., [4,11], where no constraints on dynamic model error were introduced. (ii) The adaptive estimation and RAIM methods for GPS/INS integration as presented in [18,19], respectively, could alleviate the impact of abnormal INS predictions. However, these two kinds of methods cannot make full use of the a priori predicted information by down-weighting the variance-covariance matrix of states and deleting the corresponding abnormal state components. The proposed method uses PPF to compensate for the epochs with abnormal INS predictions, which will retain all the useful information from the INS predictions. (iii) Theoretically, the multiple model approaches [20][21][22][23] could achieve consistent results compared with the proposed method since one can incorporate all possible dynamic models, which will characterize the motion behavior of an object into the alternative dynamic models. However, the multiple model methods require nearly linear computation time as a function of the number of models. We should note that only code observations were considered in the existing methods mentioned in this paragraph. The effectiveness of these methods when considering the phase observations still needs further study.
The remainder of this paper is organized as follows: in Section 2, methods for a single-frequency GPS/BDS/INS tightly coupled Kalman filter, the construction of a positional polynomial motion constraint, and the proposed PPF-aided GPS/BDS RTK and INS integration are presented. The experimental results are presented and discussed in Sections 3 and 4, respectively. Finally, in Section 5, the conclusions are summarized.

Inertial Dynamic Model
For the GPS/INS integration problem, the inertial dynamic model of the INS error model could be described as the following psi-angle error equations [4]: where δr, δv, and δψ are the position, velocity, and orientation error vectors, respectively. The term f denotes the specific force; δg is the gravity uncertainty error vector; ω ie is the rotation rate of Earth with respect to the inertial frame; ω en is the rotation rate of the navigation frame with respect to Earth; and ∇ and ε are the accelerometer error and gyro error vectors, respectively. All these terms are parameterized in the north-east-down navigation (NED) frame. The entire error states, which consist of 27 states, can be written as: Here, for the convenience of symbolic representation, the vector x is used to indicate the entire error state vector. The vectors x Nav , x Acc , x Gyro , x Grav , and x Ant represent the specific error states corresponding to navigation, accelerometer, gyroscope, gravity, and lever arm, respectively. The subscripts N, E, and D represent the navigation frame, and the subscripts bx, by, and bz represent the body frame pointing in the forward, right, and down directions; ∇ b and ∇ f represent the accelerometer bias and scale factor; ε b and ε f are the gyro bias and scale factor; and δL b represents lever arm errors. In our processing strategy, the inertial sensor bias errors and gravity errors are modeled as a first-order Gauss-Markov (GM) processes, and the lever arm errors are modeled as random constants.

Observation Model
The observation model of the relative positioning describes the relationship between the double-differenced (DD) observations and the unknown parameters. In a general form, the DD code and carrier phase observation model can be written as: where ∆∇ denotes the symbol of the DD operator. The superscript " * " is the satellite system index, which, in this paper, can be replaced by the capitals "G" and "C" for GPS and BDS, respectively. The terms ρ and φ denote the code pseudorange and carrier phase observations, ρ 0 is the geometric distance between the receiver and satellite, and T and I stand for the tropospheric and ionospheric delay parameters. The terms λ and N represent the carrier phase wavelength and integer ambiguity, M and m are the code pseudorange and carrier phase multipath errors, and ε ρ and ε φ are the code and carrier phase noises. For the RTK positioning of short baseline, the between-satellite and between-receiver DD model is widely used to eliminate the receiver-and satellite-dependent errors. Therefore, residual DD tropospheric delays after corrections with the standard tropospheric model and the ionospheric delays in Equations (5) and (6) can be ignored. Nevertheless, for long distance kinematic positioning, the residual atmospheric uncertainties cannot be neglected, which should be estimated as unknown parameters.
With the dual-constellation configuration, the strength of the integrated observation model is improved. In our GPS/BDS model, system-specific double-differentiation is adopted. Moreover, compared to those with high elevation satellites, observations at lower elevation angles suffer from more larger atmospheric and multipath errors. To account for that problem, the elevation-dependent weighting scheme σ 2 = σ 2 0 + σ 2 0 / sin 2 (e) [4] is applied for the un-differenced observations with σ 0 , the standard deviation at zenith, and e, the elevation angle.

Single-Epoch GPS/BDS RTK Ambiguity Resolution with Inertial Aiding and Tightly Coupled Kalman Filter Solution
The whole process of the tightly coupled Kalman filter (TCKF) solution can be divided into two main steps. The first step uses the least squares method to estimate the float ambiguities and then tries to fix the float ambiguities to their integers. The second step executes the Kalman filter to obtain the final solution, i.e., the navigation parameters, INS sensor biases, gravity uncertainty, and lever arm.
For the first step, we first present the linearized DD GPS/BDS observation equation and the virtual observations corresponding to the INS-predicted position, which are as follows: where l k and H k are the observation vector and the coefficient matrix, which are determined by Equation (7). Here, the terms ∆∇φ and ∆∇ρ are the carrier phase and code observations corrected with the predicted atmospheric delays; ∆∇ρ 0 is the DD geometric ranges calculated according to the INS-predicted antenna position and satellite positions; H Nav is the receiver-satellite geometry matrix; x Amb indicates the DD carrier phase ambiguity error vector; ∆∇ε INS denotes the noise vector of the virtual observations; and the subscripts n and k are the number of DD ambiguities and epoch index, respectively. The combination of the DD observation equation with the virtual observations yields the following combined observation equation for float ambiguity resolution: Using the least-squares algorithm to estimate the parameter vector of Equation (9), the float ambiguities and the corresponding variance-covariance matrix can be obtained. Then, the Least-squares AMBiguity Decorrelation Adjustment (LAMBDA) method is used to obtain the integer-valued ambiguities [28]. When they have passed the ambiguity validation test, the ambiguities can be regarded as correct integers. Then, the ambiguity-fixed carrier phase observations are adopted to adjust the 27 states of Equation (4) using the Kalman filter.
In the second step, the coefficient matrix for the Kalman filter can be written as: where C e n is the rotation matrix from the navigation frame to the earth frame, C n b is the rotation matrix from the body frame to the navigation frame, l b is the lever arm expressed in the body frame, and the subscript m denotes the number of observations.
For the observation vector y k , the DD ambiguity-fixed carrier phase and code observations are used. If the ambiguities cannot be fixed, only the code observations are used to update the state vector. Then, the error state vectorx k is estimated by the following Kalman filter formulas: where K k is the so-called gain matrix, Q k is the DD variance-covariance matrix of the observations, and x − k is the predicted state vector. Since we use the error vector to indicate the state, x − k is indeed a 27 × 1 vector with all its elements equal to zero. The term Q x k denotes the a priori variance-covariance matrix associated with x − k . The minus symbol (−) as an exponent of a variable represents a "predicted" quantity; the symbol (ˆ) above a variable represents an "estimated" quantity.
Then, the 27 × 1 state vectorX k can be obtained by subtracting the estimated error statesx k from the INS-predicted states X − k , that is: Finally, the estimated sensor biases are fed back to compensate the raw INS measurements, and the gravity uncertainty is also fed back into the INS mechanization. Then, the navigation solution is updated based on the estimated navigation states of the previous epoch.

Construction of Positional Polynomial Motion Constraint
An INS can be used to predict the position, velocity, and attitude of the object on which an inertial measurement unit (IMU) is mounted by using the outputs of an accelerometer and gyroscope measurements. However, when the IMU outputs are contaminated by gross errors, the INS-predicted position will deviate from a reasonable one. This will undoubtedly bias the upcoming ambiguity resolution using Equation (7). Therefore, if another reliable predicted position can be provided, this may guarantee the acquisition of a continuous fixed solution so as to improve the reliability and effectiveness of the integrated system.
At epoch k, the position of the rover antenna is approximated as a polynomial with a q-order fitting coefficient [25,26]. Taking all the p positions within the time window, the PPF model reads as: where δt denotes the epoch difference between successive two epochs; a i (i = 0 · · · q) denotes a 3 × 1 coefficient vector. Here, p is a predefined time window. The matrix form of Equation (15) is: ; ς is the corresponding residual error vector. Accordingly, the coefficients a in Equation (16) can be estimated with the least squares method, i.e.:â = (M T Q −1 where Qˆr (k−p:k−1) is constructed by the position variance-covariance matrix from epoch k − p to k − 1. Hence, the position predicted by PPF and its associated variance-covariance matrix at epoch k can then be obtained by the following equations: Note from Equation (19) that, when the time window and model order of PPF are set in advance, the precision of the predicted position is mainly related to the variance-covariance matrix Qˆr (k−p:k−1) . If the positions in the time window are all estimated based on ambiguity-fixed carrier phase observations, centimeter-precision position estimates will be achieved. Then, a stable and even centimeter precision of r − k will maintained.

Proposed PPF-Aided GPS/BDS RTK and INS Method
With the construction of the PPF model, the position predicted by PPF can be employed to estimate the float ambiguities, which are then fixed to their integers. Therefore, the fixed solution could be obtained with the ambiguity-fixed carrier phase observations. Obviously, the above process is similar to the conventional GNSS-only estimation but uses the PPF-predicted position instead.
In the following, for simplicity, we call this the PPF-based solution (i.e., with the 3 × 1 position vector and its variance-covariance matrix), while the tightly coupled GPS/BDS RTK and INS integration as described in Section 2.1.3 is named as the INS-based solution (i.e., with the 27 × 1 state vector and its variance-covariance matrix).
Because of the existence of outliers in IMU measurements, using these measurements to predict the current states will degrade the accuracy of float ambiguity solution; this may subsequently affect the ambiguity fixing and therefore degrade the Kalman filter solution. At this point, the PPF-based solution is used to try to provide extra position information on the current epoch. The following should be noted while combining the PPF-based solution with the INS-based solution: (i) in what conditions the constructed PPF method will be triggered to aid the ambiguity fixing failed by an INS-aided system; (ii) when there are outliers in IMU measurements, the float ambiguities estimated based on INS-predicted positions may not be fixed or may be wrongly fixed to their integers, so no ambiguity-fixed carrier phase can be used for the second step of INS-based solutions; (iii) even with the fixed ambiguities provided by the PPF model, outliers in IMU measurements will bias the INS-based solution in the sense that the INS-predicted states X − k will, to some extent, deviate. Hence, in response to the former three problems, three strategies are proposed: (i) the model probabilities corresponding to the PPF model and inertial dynamic model are calculated, which are, in this paper, used as one of the indicators to trigger or not trigger the PPF method; (ii) the fixed ambiguities based on the PPF-predicted position can be directly used to replace the float/wrongly fixed ambiguities estimated based on the INS-predicted position; (iii) an additional Kalman filter taking the PPF-based solution as virtual observations is needed, with the aim of updating the INS-based solution. As a result, in the following subsubsections, the proposed PPF-aided GPS/BDS RTK and INS integration is presented in three parts, namely, the calculation of the model probability, the substitution of float/wrongly fixed ambiguities, and updating the INS-based solution.

Model Probability Calculation
Using Bayes' rule [29], the a posteriori probability of model j when the observation is up to k is given by: where the symbol P{·} denotes the probability; the symbol p{·} denotes the probability density function; the symbol M j indicates a specific model, which in this paper can be the PPF model (M 1 ) and inertial dynamic model (M 2 ); and the epoch indexes k and k − 1 in the upper right corner of L indicate the observation set up to k and k − 1, respectively. The term l k is the observation vector at epoch k. As seen from Equation (20), when there are several conditional random variables, Bayes' rule can be used to switch only some of them. In Equation (20), the first term on the right-hand side of the last identity is the likelihood function Λ j k for M j at epoch k, which can be calculated by: where S denotes the variance-covariance matrix of the innovation vector l k of Equation (7). As seen from Equation (7), we use the error state rather than the full state as the parameter vector; the observation vector l k is essentially equivalent to the innovation vector in this case. The term P{M j |L k−1 } is the a priori probability µ j k−1 of model M j . When k = 1, µ j 0 = P{M j |L 0 } is the a priori information of that model. Consequently, Equation (20) could also be written as: where the original superscripts j = 1, 2 in µ j k or µ j k−1 are symbolized by PPF and INS, which are used to represent the model probabilities of the PPF and INS models, respectively.
The calculation of the model probabilities of the PPF and INS models is detailed. After estimating the float ambiguities of the respective model, they try to fix them to their integers using the LAMBDA. Then, the two sets of fixed or unfixed (when not passing the ambiguity validation test) ambiguities are used to correct the carrier phase observations, forming the respective innovation vectors. With the innovation vector and its corresponding innovation variance matrix H k Q x k H T k + Q k , the respective model probability can be calculated by Equation (22).
The reason why Equation (22) can be regarded as an indicator for judging the validity of an ambiguity set is that the innovation in the likelihood function considers the fixed/unfixed/wrongly fixed ambiguity vector; thus, either the unfixed or wrongly fixed ambiguity will enlarge the element in the innovation vector, reduce the value of the corresponding likelihood function, and furthermore decrease the model probability. It needs to be emphasized that the validity of an ambiguity set, in this paper, denotes whether the ambiguity set determined based on the PPF model is more credible than the ambiguity set determined based on INS, and vice versa. The term "validity" has a different meaning in the ambiguity validation test of LAMBDA [28], which is used to evaluate the most credible integer ambiguity set chosen in a candidate group searched based on a float ambiguity set and its variance-covariance matrix. Furthermore, as can be seen from Equation (22), the model probability is determined by two factors, i.e., the likelihood function and the a priori probability of each model. In other words, the model probability at the current epoch contains not only the current ambiguity fixing information but also the ambiguity fixing information of the previous epoch.
We note that if the predictions by both of the two models can correctly or cannot fix their associated ambiguity sets, the values of the likelihood functions will be in the same order for the two models. Then, the model probability of the respective model is mainly dependent on the a priori model probability. In this paper, the initial model probabilities are set to 0.9 and 0.1 for the inertial dynamic model and PPF model, respectively. If one of the ambiguity sets can be correctly fixed by the corresponding model prediction, the model probability of this model will be significantly larger than the other one. This is why one can use the proposed model probability method to detect an abnormal prediction of the corresponding model.

Substitution of INS-Based Float/Wrongly Fixed Ambiguities
The substitution of the INS-based float/wrongly fixed ambiguities by the fixed ambiguities determined based on PPF prediction is described. We should emphasize that the discussions in the following two subsubsections are all based on the idea that the ambiguity set can be correctly fixed by the PPF-predicted position. Hence, when there are outliers in IMU measurements, the ambiguities will, to a great extent, not be fixed or be wrongly fixed to their actual integers based on the INS-predicted position, resulting in deteriorated navigation performance of the INS-based solution. However, since the PPF-predicted position does not depend on the external measurements, it requires only that the positions within the time window are all available [27]. If the ambiguity set can be correctly fixed based on PPF, one can simply use these fixed ambiguities to replace the unfixed/wrongly fixed ambiguities estimated based on INS. Then, the ambiguity-fixed carrier phase observations can still be provided to the INS-based solution. In this way, the continuous high-precision observations can still be maintained. From this point of view, the fixed ambiguities provided by the PPF-predicted position can be regarded as an alternative when INS-based ambiguity fixing fails.

Updating the INS-Based Solution
Updating the INS-based solution with the fixed PPF-based solution is detailed. As has been stated above, this step will be triggered only when the INS-based ambiguity fixing has failed or been wrong.
Although the high-precision carrier phase observations are aided by PPF, the estimated state parameters in Equation (14) wherer PPF,k andr INS,k represent the PPF-based solution and the position components of INS-based solution, respectively. The terms z k and B k denote the observation vector and coefficient matrix, which are determined by the above equation. The term η k is the corresponding error vector. The negative sign in B k is due to the error states being defined by the approximated states minus the actual ones.
The symbol x LGKF,k is the state error vector of the LCKF, which is used to distinguish the TCKF state error vector x k of Equation (11). The a priori state information for the additional LCKF and the observation variance-covariance matrix are chosen from the a posteriori variance-covariance matrix of the INS-based solution and PPF-based solution, respectively. Finally, the PPF-aided INS-based solution is calculated througĥ X k −x LCKF,k . An overview of the architecture of the PPF-aided GPS/BDS RTK and INS integration is shown in Figure 1.
As shown in Figure 1, the INS mechanization is processed and the PPF model is constructed based on the consecutive historical position information. The raw GPS/BDS observations are double differenced by the individual GPS and BDS system. Two ambiguity sets are resolved with the inertial aiding and the PPF aiding, respectively. A validation test is then followed to verify the ambiguity sets. If undergoing a normal INS prediction and the ambiguities pass the validation test, the ambiguity-fixed phase observations together with the code observations are processed by the TCKF to obtain the INS-based fixed solution. Meanwhile, a GNSS-only solution using the ambiguity-fixed high-precision phase and code observations is estimated by a least squares method, i.e., a PPF-based solution. In this time, no assistance from the PPF-related information (i.e., the fixed ambiguities and the PPF-based solution) is required. When there are outliers in the IMU measurements, the INS predictions may be biased and the biased predictions will result in a failed AR process (i.e., the AR validation test fails, or the AR validation test passes but the ambiguity set is wrongly fixed). Therefore, two criteria are used to judge whether to trigger the steps related to the PPF-aided INS-based solution. The first criterion guarantees that the PPF-based AR is valid. The second criterion states whether the model probability µ PPF k calculated based on PPF is larger than the INS counterpart µ INS k . The criterion of the model probability µ PPF k > µ INS k used here guarantees that the PPF-based AR is statistically significant compared to the INS-based AR. Hence, updating the INS-based solution with the LCKF will improve the INS-based solution. We should note that there is also a situation in which no gross errors exist in the IMU measurements and the ambiguity fixing states of the INS and PPF model are both failed. In this situation, the PPF-aided INS-based solution will not be triggered. This can be easily derived from the two criteria. The failure of the AR processes may be owing to the quality of the observations and/or the constellation of satellites. At this point, only the code observations are used. The proposed method is expected to provide continuous and highly reliable navigation performance.
solution, respectively. The terms k z and k B denote the observation vector and coefficient matrix, which are determined by the above equation. The term k η is the corresponding error vector. The negative sign in k B is due to the error states being defined by the approximated states minus the actual ones. The symbol LCKF,k x is the state error vector of the LCKF, which is used to distinguish the TCKF state error vector k x of Equation (11).
The a priori state information for the additional LCKF and the observation variance-covariance matrix are chosen from the a posteriori variance-covariance matrix of the INS-based solution and PPF-based solution, respectively. Finally, the PPF-aided INS-based solution is calculated through . An overview of the architecture of the PPF-aided GPS/BDS RTK and INS integration is shown in Figure 1.

Field Test Description
A vehicular test was carried out to evaluate the performance of the developed FFP-aided GPS/BDS RTK and INS tightly coupled method in an urban environment (see Figure 2). Seventy minutes of rover station data with 100 Hz INS data and 1 Hz dual-frequency GPS/BDS data were collected by a NovAtel SPAN-CPT system and a high-grade geodetic GNSS dual-frequency receiver on 19 September 2014. Another same receiver mounted on the roof of a tall building was set up as the reference station. The baseline length of the whole experiment was less than 5 km. The test covered the typical area characterized by foliage, a sub-dense urban canyon, a lake, and an open-sky condition. Hence, the satellite positioning performance was limited because of signal blocking and multipath problems. In this test, the single-frequency GPS/BDS (GPS L1 and BDS B1 phase and code observations) and INS data were processed to validate the proposed method using a single-epoch mode, while the dual-frequency ambiguity-fixed GPS/BDS RTK and INS solution was processed for reference. For evaluating the ambiguity, one can compare the position components of the proposed method with the counterparts of the reference trajectory through a predefined threshold to judge whether the corresponding ambiguities are correct or incorrect fixing. The threshold used in this paper was 0.05 m for three position components. Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 27 From Figure 2b, it can be seen that the number of visible satellites of BDS was marginally larger than that of the GPS counterparts during the whole test. From the perspective of position dilution of precision (PDOP) values, however, the PDOPs of BDS are larger than those of the GPS ones, which is due to the current deficiency of the BDS geometry configuration. By combining GPS with BDS, the number of visible satellites is doubled compared with in the GNSS-only system. The PDOPs are significantly improved for the combined GPS and BDS constellation. Figures 2c-d show the velocities in the navigation frame and the vehicle attitudes. It can be seen that the vehicle was moving in a moderate dynamic.
During the test, the standard deviation 0 σ at zenith was empirically set to 3 mm and 0.3 m for the phase and code observations, respectively, and the cutoff elevation was set to 15°. The specifications of the inertial sensor are listed in Table 1.  From Figure 2b, it can be seen that the number of visible satellites of BDS was marginally larger than that of the GPS counterparts during the whole test. From the perspective of position dilution of precision (PDOP) values, however, the PDOPs of BDS are larger than those of the GPS ones, which is due to the current deficiency of the BDS geometry configuration. By combining GPS with BDS, the number of visible satellites is doubled compared with in the GNSS-only system. The PDOPs are significantly improved for the combined GPS and BDS constellation. Figure 2c,d show the velocities in the navigation frame and the vehicle attitudes. It can be seen that the vehicle was moving in a moderate dynamic.

Positional Polynomial Fitting Performance by Simulated Vehicle Trajectories
During the test, the standard deviation σ 0 at zenith was empirically set to 3 mm and 0.3 m for the phase and code observations, respectively, and the cutoff elevation was set to 15 • . The specifications of the inertial sensor are listed in Table 1.

Positional Polynomial Fitting Performance by Simulated Vehicle Trajectories
In order to aid the ambiguity fixing so as to obtain continuous and highly reliable navigation solutions, the necessity of subdecimeter-precision-predicted positions for reliable ambiguity resolution should be confirmed by positional polynomial fitting. In this subsection, both the fitting errors and prediction errors of PPF are presented. For the fitting errors, it is calculated by the mean value of the differences between the fitted position coordinates Mâ of Equation (16) and the corresponding truth position coordinates within the time window. The prediction errors indicate the difference between the predicted position r − k of Equation (18) and the truth position of epoch k. How to select the time window p and model order q was studied in [25,26]. It was pointed out that the motion of a vehicle often obeys a low-order polynomial over a short period [27]. Therefore, in this subsection, we chose two sets of polynomial fitting model characterized by p = 4, q = 2 and p = 5, q = 3 to build the kinematic model of PPF, respectively. A straight line vehicle trajectory with variable accelerated motion similar to [27] and a combined straight and curved vehicle trajectory with accelerated and decelerated motion were simulated. All the position coordinates were expressed in the earth frame and then transformed to the navigation frame to show the vehicle trajectories in Figure 3.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 27 In order to aid the ambiguity fixing so as to obtain continuous and highly reliable navigation solutions, the necessity of subdecimeter-precision-predicted positions for reliable ambiguity resolution should be confirmed by positional polynomial fitting. In this subsection, both the fitting errors and prediction errors of PPF are presented. For the fitting errors, it is calculated by the mean value of the differences between the fitted position coordinates Ma of Equation (16)  How to select the time window p and model order q was studied in [25,26]. It was pointed out that the motion of a vehicle often obeys a low-order polynomial over a short period [27]. Therefore, in this subsection, we chose two sets of polynomial fitting model characterized by p = 4, q = 2 and p = 5, q = 3 to build the kinematic model of PPF, respectively. A straight line vehicle trajectory with variable accelerated motion similar to [27] and a combined straight and curved vehicle trajectory with accelerated and decelerated motion were simulated. All the position coordinates were expressed in the earth frame and then transformed to the navigation frame to show the vehicle trajectories in Figure 3. As seen from Figure 3, the vehicle is simulated to move with a moderate velocity in these two cases, to try to match the actual motion in an urban environment. The description of the straight line vehicle trajectory with variable accelerated motion is separated into the following nine scenarios: S1 (  As seen from Figure 3, the vehicle is simulated to move with a moderate velocity in these two cases, to try to match the actual motion in an urban environment. The description of the straight line vehicle trajectory with variable accelerated motion is separated into the following nine scenarios: S1

Full Ambiguity Resolution and Positioning Performance
In terms of the

Full Ambiguity Resolution and Positioning Performance
In terms of the

Full Ambiguity Resolution and Positioning Performance
In terms of the different satellite system combinations, the INS-based solutions, in this paper, can be classified into the following three schemes:  [25] in integrated navigation with phase observations. Other methods that could be used to handle the abnormal predictions from the dynamic system-such as the adaptive estimation method [18], RAIM method [19], and multiple model method [23]-are not considered here for numerical analyses due to their shortcomings as summarized in Section 1. Table 2 lists the number of epochs with fully and incorrectly fixed ambiguities and their fixing and incorrect fixing rates. In Table 2, the fixing rate (P FR ) and incorrect fixing rate (P IFR ) are defined by: P FR = number of epochs with fixed ambiguities total number of epochs × 100% P IFR = number of epochs with incorrectly fixed ambiguities total number of epochs with fixed ambiguities × 100% It should be noted that the number of epochs with fixed ambiguities mentioned here indicates the number of epochs for which the full set of ambiguities is fixed to their integers. This is used to distinguish the number of epochs with partially fixed ambiguities mentioned in the next subsection, which indicates the number of epochs for which only a subset of ambiguities is fixed to integers.
To further investigate the effectiveness of the schemes with the PPF-aided method, an error with a size of 0.05 g was added to the output of the accelerometer z-axis at GPS time 460,270 s. The position errors and ambiguity fixing states around this epoch are shown in Figure 6. Figure 7 shows the model probabilities calculated by the three PPF-based schemes.

Partial Ambiguity Resolution and Positioning Performance
To further improve the ambiguity fixing rate so as to achieve continuous high-precision positioning results, the partial ambiguity resolution (PAR) strategy is added to the ambiguity fixing processes of INS-based AR and PPF-based AR. The PAR used in this paper is conducted in an iterative way, which uses an already-fixed subset of ambiguities as precise ranges to resolve the

Partial Ambiguity Resolution and Positioning Performance
To further improve the ambiguity fixing rate so as to achieve continuous high-precision positioning results, the partial ambiguity resolution (PAR) strategy is added to the ambiguity fixing processes of INS-based AR and PPF-based AR. The PAR used in this paper is conducted in an iterative way, which uses an already-fixed subset of ambiguities as precise ranges to resolve the

Partial Ambiguity Resolution and Positioning Performance
To further improve the ambiguity fixing rate so as to achieve continuous high-precision positioning results, the partial ambiguity resolution (PAR) strategy is added to the ambiguity fixing processes of INS-based AR and PPF-based AR. The PAR used in this paper is conducted in an iterative way, which uses an already-fixed subset of ambiguities as precise ranges to resolve the remaining unfixed ambiguities until all the ambiguities are fixed or there are no more ambiguities that can be fixed.
The vital step for PAR is the selection of an ambiguity subset that can be reliably fixed. In this paper, we chose the commonly used strategy to determine the subset of ambiguities by successively increased elevations. The PAR is only activated if the full ambiguity resolution fails. If FAR fails, the cutoff elevation is then increased by 5 • . Therefore, the subset of ambiguities is chosen based on the new cutoff elevation threshold. The LAMBDA method is used to try to fix the subset ambiguities to their integers. If the fixed ambiguity subset does not pass the ratio test and the success rate is lower than the predefined value, this procedure is repeated until a selected ambiguity subset that can be successfully fixed or is empty is obtained. If the fixed ambiguity subset passes the ratio test and the success rate fulfills the predefined value, we think that the subset ambiguities are correctly fixed. Then, the ambiguity-fixed carrier phase observations can be used as precise ranges and aid the fixing of the remaining ambiguities. If none of the ambiguities can be fixed, only the code observations are used and the float solution of the state vector is achieved. The ratio test value and the success rate for PAR are set to 3 and 0.99, respectively. Table 3 lists the single-frequency single-epoch ambiguity fixing states of each scheme using PAR. To demonstrate a clear effect of PPF in aiding GNSS/INS solutions by using PAR, in this subsection, only the integrated solutions are considered. In Table 3, an additional column containing the number of epochs with partially fixed ambiguities and the partial fixing rate is presented, compared with Table 2. The partial fixing rate (P PFR ) is computed as: P PFR = number of epochs with partially fixed ambiguities total number of epochs × 100% (26) Table 3. Single-frequency single-epoch ambiguity fixing states with different schemes using partial ambiguity resolution (PAR). In order to manifest the advantage of the proposed PPF-aided method, errors with a size of 0.05 g were added to each output of the accelerometer axes at GPS time 460,270 s in the same way as in the second simulation in Section 3.2. The position errors and ambiguity fixing states around GPS time 460,270 s estimated by the schemes using PAR are shown in Figure 10. The corresponding model probabilities calculated by the schemes with the PPF-aided method using PAR are shown in Figure 11.  To give an overall comparison of the number of ambiguities fixed by PAR and FAR, the ratio of the number of ambiguities ( FR P Δ ) fixed by PAR with respect to those fixed by FAR is defined as:   To give an overall comparison of the number of ambiguities fixed by PAR and FAR, the ratio of the number of ambiguities ( FR P Δ ) fixed by PAR with respect to those fixed by FAR is defined as: To give an overall comparison of the number of ambiguities fixed by PAR and FAR, the ratio of the number of ambiguities (∆P FR ) fixed by PAR with respect to those fixed by FAR is defined as:

Number of Epochs with
∆P FR = number of ambiguities fixed by PAR-number of ambiguities fixed by FAR number of total ambiguities × 100% (27) The ∆P FR results corresponding to the integrated solutions are shown in Figure 12.
The FR P Δ results corresponding to the integrated solutions are shown in Figure 12. Finally, in Figures 13 to 15, the performance of each PPF-based scheme using PAR is evaluated by analyzing the positioning errors and the number of fixed DD ambiguities.   Finally, in Figures 13 to 15, the performance of each PPF-based scheme using PAR is evaluated by analyzing the positioning errors and the number of fixed DD ambiguities.

Discussion on the Performance of the Positional Polynomial Fitting Constraint
Generally, the fitting errors shown in the top two subplots of Figures 4 and 5 are at the centimeter level, which verifies that the PPF model can be well established by the respective model order and its associated time window. However, the prediction accuracy plays a key role in aiding the ambiguity resolution.
As can be seen from the bottom two subplots of Figure 4, the prediction errors are at the subdecimeter level, except for some epochs with large dynamic changes. Even though the prediction errors of the second-order PPF model are marginally smaller than those of the third-order PPF model when there is no dynamic change, the prediction errors of the second-order PPF model are larger than those of the third-order PPF model when undergoing dynamic change, especially for the epochs during S4 and S6, where variable accelerated motions are implemented. This can be seen from the bottom-left subplot of Figure 4; the prediction errors during S4 and S6 (labeled by black frames) fluctuate, while this phenomenon cannot be explicitly seen from the bottom-right subplot of Figure  4. This is because a third-order PPF model can better characterize the variable accelerated behavior than a lower order PPF model.
Similar results can also be seen from the bottom two subplots of Figure 5, where the combined straight and curved vehicle trajectory is considered. We should note that the position coordinates are

Discussion on the Performance of the Positional Polynomial Fitting Constraint
Generally, the fitting errors shown in the top two subplots of Figures 4 and 5 are at the centimeter level, which verifies that the PPF model can be well established by the respective model order and its associated time window. However, the prediction accuracy plays a key role in aiding the ambiguity resolution.
As can be seen from the bottom two subplots of Figure 4, the prediction errors are at the subdecimeter level, except for some epochs with large dynamic changes. Even though the prediction errors of the second-order PPF model are marginally smaller than those of the third-order PPF model when there is no dynamic change, the prediction errors of the second-order PPF model are larger than those of the third-order PPF model when undergoing dynamic change, especially for the epochs during S4 and S6, where variable accelerated motions are implemented. This can be seen from the bottom-left subplot of Figure 4; the prediction errors during S4 and S6 (labeled by black frames) fluctuate, while this phenomenon cannot be explicitly seen from the bottom-right subplot of Figure  4. This is because a third-order PPF model can better characterize the variable accelerated behavior than a lower order PPF model.
Similar results can also be seen from the bottom two subplots of Figure 5, where the combined straight and curved vehicle trajectory is considered. We should note that the position coordinates are

Discussion on the Performance of the Positional Polynomial Fitting Constraint
Generally, the fitting errors shown in the top two subplots of Figures 4 and 5 are at the centimeter level, which verifies that the PPF model can be well established by the respective model order and its associated time window. However, the prediction accuracy plays a key role in aiding the ambiguity resolution.
As can be seen from the bottom two subplots of Figure 4, the prediction errors are at the subdecimeter level, except for some epochs with large dynamic changes. Even though the prediction errors of the second-order PPF model are marginally smaller than those of the third-order PPF model when there is no dynamic change, the prediction errors of the second-order PPF model are larger than those of the third-order PPF model when undergoing dynamic change, especially for the epochs during S4 and S6, where variable accelerated motions are implemented. This can be seen from the bottom-left subplot of Figure 4; the prediction errors during S4 and S6 (labeled by black frames) fluctuate, while this phenomenon cannot be explicitly seen from the bottom-right subplot of Figure 4. This is because a third-order PPF model can better characterize the variable accelerated behavior than a lower order PPF model. Similar results can also be seen from the bottom two subplots of Figure 5, where the combined straight and curved vehicle trajectory is considered. We should note that the position coordinates are badly predicted by the second-order PPF model when carrying out a circle motion with the speed at nearly 10 m/s. Nevertheless, the third-order PPF model can still specify this complicated movement. Overall, considering the situation of our field test, a third-order PPF model with the time window set to five is employed.

Discussion on the Full Ambiguity Resolution and Positioning Performance
As can be seen from Table 2, the fixing rates for single-frequency single-epoch GPS and BDS are low, only 18.38% for GPS and 18.72% for BDS. These can also be seen from the number of epochs with fixed ambiguities for GPS and BDS. Moreover, the incorrect fixing rates of GPS and BDS reach 5.43% and 5.46%, respectively. All these results show that it is difficult for a single-frequency GNSS-only method to achieve reliable dynamic positioning results, considering the limited ambiguity resolutions of them. However, with the combined GPS and BDS system, the fixing rate is greatly improved, reaching 81.38%. This can be attributed to the enhanced satellite geometry and the redundancy of observations, which improve the precision of the float ambiguity estimates. The benefit of the GPS/BDS system is even more obvious in terms of the incorrect fixing rate, reaching 0% for the combined system.
If GPS or BDS is aided by INS, the fixing rate can be improved to more than 80%, exceeding the fixing rate of the GPS and BDS combined system in our experiment. The incorrect fixing rates of GPS/INS and BDS/INS are greatly decreased compared with the GNSS-only case. It is thus evident that the INS-predicted position information benefits the AR performance. This is because the additional INS-predicted position information enhances the variance-covariance matrix of the estimated float ambiguities, which defines the ambiguity search space and the quality of the float ambiguity estimates. However, the incorrect fixing rate of GPS/INS is higher than that of BDS/INS. We find from the GPS/INS solution that the epochs with incorrectly fixed ambiguities are concentrated in a successive period of epochs. This may be due to the quality of the observations during those epochs, and some quality control methods should be introduced to delete those deficient observations in future work. If PPF is constructed and is regarded as a dynamic model for GNSS, then the a priori information of the position could be provided so as to assist the ambiguity solution. In Table 2, slightly more epochs could be fixed by GPS and BDS-only schemes with PPF aid. Furthermore, the incorrect fixing rates of the two schemes are also improved compared with those of their respective counterparts without PPF aid. However, the combined GPS and BDS solution with PPF aid becomes worse than the GPS/BDS solution in terms of the ambiguity fixing rate. Based on these results, one can conclude that only taking the PPF as a dynamic model cannot stably improve the ambiguity fixing rate. This is because the PPF is vulnerable to the position accuracy within the time window and the motion of the object.
Hence, if the object's motion is mainly described by INS and the PPF is regarded as a subordinate model, which will only be triggered when INS model cannot work, the results are improved. We find that the fixing rates of PPF-aided GPS/INS, BDS/INS and GPS/BDS/INS are slightly higher than those of their counterparts without PPF aid. Furthermore, the incorrect fixing rates of PPF-aided GPS/INS and BDS/INS are lower than those of their counterparts without PPF aid. The incorrect fixing rate of PPF-aided GPS/BDS/INS still remains zero. The improvements of the schemes with PPF aid reflect that the INS-predicted position, in some epochs, does not contribute to successful ambiguity fixing, while, in these epochs, the PPF-aided method is triggered, providing correctly fixed ambiguities and updating the INS-based solutions. The effectiveness of the schemes with PPF aid may be more significant when a lower-cost IMU is equipped.
As seen from Figure 6, the position errors for the six schemes are all at the centimeter level before GPS time 460,270 s; in other words, with a reliable INS-predicted position, all six schemes could achieve reliable ambiguity-fixed solutions. However, when there are outliers in the IMU measurements, the navigation performances of the three INS-based solutions visibly deteriorated; at the same time, the ambiguity fixing states became unfixed. Furthermore, for the GPS/INS and BDS/INS schemes, it will take a long time to re-fix the unfixed ambiguities than with the GPS/BDS/INS scheme, as shown in the left three subplots of Figure 6.
When it is detected that there are outliers in the IMU measurements, the PPF-predicted position can still provide reliable position information so as to aid and update the unfixed ambiguities as well as the biased INS-based solutions. The right three subplots of Figure 6 show that continuous and reliable positioning results can be achieved by the three schemes with PPF aid. However, there are a few epochs whose ambiguity sets cannot be fully fixed. This may be due to the specific processing scheme and the quality of the observations of those epochs. However, the ambiguity sets can be immediately fixed right after the unfixed epochs. This is because the ambiguities, in this paper, are resolved in an epoch-by-epoch mode, and when there are no abnormal dynamic predictions, the ambiguity fixing state at previous epochs will not affect the AR performance at the current epoch. Furthermore, the positioning results for the epochs after a g error is introduced will be roughly equivalent to the results where no abnormal INS predictions occur at GPS time 460,270 s.
The model probabilities of the PPF and inertial dynamic models calculated by the three PPF-aided RTK/INS integrated schemes are shown in Figure 7. When the outputs of the IMU contain gross errors at GPS time 460,270 s, the model probability of the INS suddenly decreases to 0.008, while the model probability of PPF increases to 0.992. After one to two epochs, however, the model probability corresponding to INS returns to 0.999. This verifies that with the aid of PPF, the impact of outliers in the IMU measurements on subsequent epochs is greatly reduced. Consequently, the PPF-aided system quickly restores to an INS-aided integration system from the perspective of model probability. The proposed model probability method has the ability to sensitively sense the abnormal predictions caused by IMU measurements with gross errors.
When errors with a size of 0.05 g were added to each output of the accelerometer axes at GPS time 460,270 s, similar positioning performance could be seen as shown by Figure 8 compared with the subplots of Figure 6. The relatively larger outliers in the IMU measurements do not impact the overall positioning performance of the PPF-aided RTK/INS integrated schemes but have a slight impact on the AR performance. This can be seen from the right subplots of Figure 8, where more than one unfixed epoch is found in the GPS/INS-PPF and BDS/INS-PPF schemes, compared with the corresponding subplots of Figure 6. In Figure 8, the positioning results for the epochs after re-convergence are also similar to their counterparts in Figure 6. We note that once the results are convergent, the positioning results for converged epochs will be at the centimeter level if no other factors affect successful AR. From Figure 9, it can be seen that the variation of the curve still verifies the effectiveness of the calculated model probabilities. Theoretically, the larger the outliers in the IMU measurements, the more easily the model probability reflects the actual situation.

Discussion on the Partial Ambiguity Resolution and Positioning Performance
By comparing the second columns of Tables 2 and 3, it can be seen that the number of epochs with fully fixed ambiguities with each scheme using PAR is larger than those using FAR. The fixing rates of all the schemes are higher than 90% when the PAR technique is adopted. Moreover, with the PAR, there are epochs where the subset of ambiguities is fixed. The second row of each scheme lists the results by numerically adding the fully fixed results to the partially fixed results. We can find that more than 95% of the epochs can achieve fixed solutions by using PAR, especially for the GPS/INS, GPS/BDS/INS, GPS/INS-PPF, and GPS/BDS/INS-PPF schemes; the full and partial fixing rates of them exceed 99%. In terms of the incorrect fixing rate, the GPS/INS and GPS/INS-PPF schemes using PAR have a great improvement compared with their counterparts using FAR. However, the number of epochs with incorrectly fixed ambiguities obtained by BDS/INS-PPF using PAR is slightly more than that obtained by BDS/INS-PPF using FAR. However, considering that relatively many epochs of fixed solutions can be obtained by the BDS/INS-PPF scheme using PAR, the deficiency in the incorrectly fixed rate is moderate. The most reliable and continuous high-precision positioning results can be achieved by the GPS/BDS/INS and GPS/BDS/INS-PPF schemes using PAR. It is noteworthy that the performance of these two schemes, in this test, is equivalent. This could be attributed to three factors, namely, the strengthened satellite geometry of the combined GPS and BDS system, the utilization of the PAR strategy, and, especially, the fact that there were no significant outliers in the IMU measurements during the whole test. Therefore, the constructed PPF model could hardly be triggered.
When errors with a size of 0.05 g were added to each output of the accelerometer axes at GPS time 460,270 s, the three INS-based solutions were biased due to the outliers in IMU measurements, even with the PAR strategy, as shown by the left three subplots of Figure 10. The PAR can help to achieve continuous and reliable positioning results but is not essential for the adaptivity of an integrated system. Although there are some epochs between the last fixed epoch and the first re-fixed epoch where the ambiguity subsets can be fixed, it is not sufficient to converge the divergent positioning results. With the aid of the PPF model, the continuity and reliability of the positioning results are maintained. The centimeter-level positioning errors are shown in the right three subplots of Figure 10. We should note that the ambiguity fixing states have an obvious improvement compared with the right three subplots of Figure 8, where the FAR strategy is adopted. There are fewer epochs with the fixing states equal to zero when combining PAR with PPF-aided RTK/INS integrated schemes. Similar to the schemes using FAR, once the results are convergent after GPS time 460,270 s, the ambiguity-fixed positioning results will hold if no other factors affect successful AR.
As shown in Figure 11, the model probabilities calculated by the PPF-aided RTK/INS integrated schemes also verify the ability of the proposed model probability method to detect the outliers in the IMU measurements. It is interesting that there are two changes in the model probability curve in the first and third subplots of Figure 11, indicating that the PPF method is triggered twice during this time period. The first change takes place when there are outliers in the IMU measurements as has been explained in Section 4.2, while the second change is due to the use of PAR. The reasons for the second change in the model probability curves can be explained by the following aspects. Taking the GPS/BDS/INS-PPF scheme using FAR and the same scheme using PAR as examples, we can find that there is one epoch with its ambiguity set unfixed in the last subplot of Figure 8, while the ambiguity set of that epoch can be partially fixed as shown in the last subplot of Figure 10. This is because the likelihood function calculated based on the ambiguity-unfixed phase observations is numerically smaller than the counterpart calculated based on the partially fixed phase observations. Furthermore, as can be seen from the last subplot of Figure 11, the second changes for the two curves reflect that there are no ambiguities or wrong ambiguities that can be fixed by GPS/BDS/INS (PAR), while the PPF-aided GPS/BDS/INS (PAR) scheme can correctly fix a subset of ambiguities. Therefore, in that epoch, the model probability based on PPF is significantly larger than the one based on INS. This change in the model probability is found to be consistent with the actual experimental situation when checking the test date, which again verifies the sensitivity of the model probability method in detecting situations with abnormal predictions.
In Figure 12, the ∆P FR results of the six schemes are greater than or equal to zero in the majority of epochs. This means that the PAR strategy does have the ability to fix more ambiguities than FAR, particularly for the GPS/BDS/INS and PPF-aided GPS/BDS/INS schemes. For the dual-constellation and INS configurations, the number of epochs with ∆P FR = 100% is significantly larger than the that for the single-constellation and INS configurations, which shows that with an additional satellite system, the PAR strategy is rather beneficial for fixing the whole set of ambiguities. There are two to four epochs with ∆P FR < 0 when using the single-constellation and INS configurations, which indicates that in these epochs, the FAR can fix all the ambiguities while the PAR can only fix a subset of the ambiguities. However, when checking the positioning results from the epochs with ∆P FR < 0, the positioning errors using FAR reach the decimeter level, which reveals that wrongly fixed ambiguities exist in those epochs. The positioning errors using the PAR in those epochs are at the centimeter level with the exception of one epoch with ∆P FR = −14.29% at GPS time 460,679 s of the BDS/INS and PPF-aided BDS/INS schemes, whose positioning errors reach the decimeter level. All in all, the schemes using PAR tend to fix more reliable ambiguities than the FAR strategy, which is of great benefit when continuous and reliable high precision positioning solutions are needed. On the other hand, continuous and reliable high precision positioning solutions are also a prerequisite for the reliable prediction of the PPF model. In this aspect, the PAR strategy is more attractive when a PPF model is considered.
The positioning results of the three PPF-aided RTK/INS integrated schemes are shown in Figures 13-15. As can be seen from Figures 13 and 14, the positioning errors are generally at the centimeter level except for a few epochs where the errors are larger than 0.1 m. Epochs with larger positioning errors are usually found when fewer ambiguities can be fixed. The most attractive positioning results are achieved with the PPF-aided GPS/BDS/INS scheme using PAR as shown in Figure 15. The positioning errors are all within 0.06 m during the whole test period, and at least six DD ambiguities can be fixed at each epoch. We note that the better performance of the PPF-aided GPS/BDS/INS scheme by using PAR is attributed to three aspects: (i) the PAR strategy used in this paper can help to fix as many ambiguities as it could; (ii) an additional satellite system can strengthen the satellite geometry and improves the redundancy of the observations; and, furthermore, (iii) with the aid of the PPF model, the adaptivity of the integrated system is also guaranteed. Therefore, good AR performance and positioning results could be obtained from the PPF-aided dual-constellation and INS configuration by using PAR.

Conclusions
Reliable and continuous ambiguity resolution is a prerequisite for high-precision positioning using carrier phase observations. However, a challenge is efficient ambiguity and positioning resolution in a single-frequency GNSS/INS integrated system with abnormal INS predictions, especially for the use of a low-cost IMU. In this paper, we have proposed a PPF-aided GPS/BDS RTK and INS integration aiming to compensate for the epochs with abnormal INS predictions. The key to using PPF is to use its predicted current position information to fulfill a successful AR process so as to accomplish an independent GNSS-only solution with ambiguity-fixed phase observations (i.e., the PPF-based solution). The main purpose of PPF is to use its fixed ambiguities and the PPF-based solution to replace the unfixed ambiguities and update the biased float solution of GPS/BDS RTK and INS integration. Moreover, a partial ambiguity resolution strategy is adopted to replace the FAR in order to further improve the ambiguity fixing rate and the positioning solutions. The performance of PPF-based GPS/BDS RTK and INS integration has been evaluated in terms of ambiguity fixing states and positioning solutions through a field vehicular test. The following conclusions can be drawn: (i) With the aid of PPF, the fixing rates of the GPS/INS, BDS/INS, and GPS/BDS/INS schemes using FAR are all greater than 86%, higher than their counterparts without PPF aid. Furthermore, the incorrect fixing rates of PPF-aided GPS/INS and BDS/INS are lower than the ones without PPF aid. However, the incorrect fixing rate of PPF-aided GPS/BDS/INS remains 0%. All these results indicate that the proposed PPF method indeed benefits the improvements in AR performance. (ii) Further improvement of AR and positioning performance could be obtained by resolving a subset of ambiguities instead of the full set. A success rate greater than 90% can be obtained for all singleor dual-constellation and INS integrations with or without PPF aid. In terms of the incorrect fixing rate, it remains at 0% if dual-constellation systems are considered. (iii) If abnormal INS-predicted positions occur, the PPF-aided method can provide ambiguity-fixed phase observations and can help to update the float solution of the GNSS/INS integrated system. As a result, a continuous ambiguity-fixed resolution as well as high-precision and reliable positioning results can be maintained. (iv) In order to maintain the reliability of PPF prediction, the continuous and precise positions within the time window are needed. Furthermore, a more continuous and precise positioning solution can be provided with the help of PAR than with the FAR strategy. Therefore, the most attractive option is to combine the PAR with the PPF-aided dual-constellation and INS integration.
To further improve the applicability of the PPF method for a GNSS/INS integrated system, the construction of PPF with different model orders for different dynamic motion should be taken into account. Further research could be pursued by recommending some machine learning algorithms to detect the vehicle motion state, and then the corresponding model order of PPF could be selected on this basis. Moreover, the prediction accuracy of PPF is highly related to the sampling interval, high-rate GNSS data, therefore, will further improve the effect of the proposed method.