A novel ambiguity resolution model of BeiDou navigation satellite system/inertial navigation system tightly coupled for kinematic-to-kinematic precise relative positioning

The carrier-based kinematic-to-kinematic relative positioning can obtain the precise baseline between two moving stations, which greatly expands the application field of dynamic relative positioning. However, the relative positioning performance is degraded greatly with low fixation rate of ambiguity with low-cost receivers. Especially, in the complex dynamic environment, ambiguity resolution effect is influenced by the satellite signal blocked, multipath outlier, and abnormal state prediction. Aiming at the problems, a novel inertial navigation system–aided robust adaptive filtering ambiguity resolution model is proposed. In addition, a hierarchical filtering strategy is developed to eliminate ambiguity parameters in BeiDou navigation satellite system/inertial navigation system tightly coupled integrated system. Finally, the precise relative position can be calculated with the “best” ambiguity solution. Both experiments with static data and field vehicle test were carried out to evaluate the algorithm efficiency in different data configurations. The results indicate that IRAFAR-TCRP method can effectively suppress the influence of observation outliers and model prediction abnormalities, which improves the success rate of ambiguity resolution, raises the accuracy as well as the continuity of relative positioning. The success rate of ambiguity resolution with single-frequency BeiDou navigation satellite system can reach 90% in the gross error and abnormal disturbance environments and centimeter-level accuracy can be achieved.


Introduction
Traditional global navigation satellite system (GNSS) real-time kinematic (RTK) technology using differential carrier phase observations can achieve centimeterlevel accuracy with correct and reliable ambiguity solution. However, it is limited by the fixed position of the reference station and limited coverage. When the rover station is far away from the reference station, the traditional RTK technology is difficult to meet the requirement, such as aircraft formation flying, unmanned air landing, and so on. Aiming at this problem, the relative positioning method based on moving reference station (also called kinematic-to-kinematic relative positioning, KKRP) is used. [1][2][3] In the carrier-based GNSS precision relative positioning, the ambiguity resolution (AR) is a key issue. Dual/triple-frequency-based RTK technology utilized for AR strategy has attracted much attention in recent years. Melbourne-Wu¨bbena (MW) or wide-lane (WL) combination method and three-carrier ambiguity resolution (TCAR) method are the most popular frequency solution algorithms, which can achieve rapid and reliable AR in various baseline range. [4][5][6][7] With the rapid development of multiple GNSS constellations, the multi-GNSS-based RTK method is also a research hotspot, which brings benefits for GNSS positioning. [8][9][10] However, it remains a challenging problem for lowcost receivers with single frequency for single-epoch solving in dynamic positioning. Especially, in GNSS-challenging surroundings, the ambiguity fixing effect is reduced significantly due to the GNSS signal occlusion, interruption, and multipath outliers, which cause poor geometry, decreased number of visible satellites, and low ambiguity float accuracy. [11][12][13] To solve the above problems, GNSS is often integrated with inertial navigation system (INS), which can provide short-term high-precision position information and strong autonomy. Therefore, the GNSS/INS integrated system can maintain high-precision navigation under complex dynamic environments. [14][15][16][17] There are two representative GNSS/INS integration methods with loosely coupled (LC) and tightly coupled (TC). Among them, the LC method is simple implementation but unreliability with less than four visible satellites; 18 while the TC method offers theoretical soundness and high accuracy, which is used in various fields such as airborne navigation, intelligent transportation, and precision agriculture. [19][20][21] In the carrier-based GNSS/INS TC integrated system, ambiguity parameter estimation is usually implemented in two ways. One way is to augment the ambiguity parameter into the state equation of integrated system and estimate the ambiguity floats in a single integration filter, which imposes heavy computational burden due to the large number of system states to be estimated. In addition, it may lead to longer solution time and poor positioning accuracy in the harsh environments, even causing filter divergence with cycle slip. The other way is to separate the ambiguity parameter from the integrated system, in which the least squares (LS) or extended Kalman filter (EKF) is always utilized to estimate the ambiguity floats first. Then, the LS ambiguity decorrelation adjustment (LAMBDA) algorithm is used to search and verify the integer ambiguity with corresponding covariance matrix. 22,23 Subsequently, the integrated system is performed with pseudorange and carrier phase of known ambiguity observations. In this method, an instantaneous high-precision position information derived by INS can be used to aid GNSS AR in turn. It has been theoretically and experimentally verified that the adding of INS information can aid AR in two aspects. [24][25][26][27] On one hand, the more precise satellite-to-ground distance predicted by INS is used to replace the geometrical distance by single-point positioning, which can improve the accuracy of ambiguity float. On the other hand, a virtual observation provided by INS is used to enhance the positioning model strength, which can reduce the ambiguity search space. Since the model reliability depends on the accuracy of INS-derived position and AR efficiency relies on observations quality as well, various parameter estimation methods are developed to improve the system performance. For the biases in GNSS observations, such as atmospheric residual error, pseudorange multipath gross errors, and large observation noise, the robust estimation method is always utilized to relieve the effects, in which the outlier observations are reweighted by the robust factor designed by predicted (innovation)-based residual test. 28,29 To alleviate the negative effects by occasional outlier in INS measurements, positional polynomial fitting constraint method is used to compensate for the epochs with abnormal INS predictions. 30 For the biases in predicted states in EKF, such as inaccuracy of motion model and abnormal disturbance in states, the adaptive filtering estimation method is generally used to reduce the effects, in which the abnormal predicted states are decreased weight by the adaptive factor designed by states discrepancy-based residual test. 31,32 To minimize the impact of predicted outlier but meanwhile maintain the prediction accuracy, multiple model adaptive estimation method is exploited, in which a fusion solution by multi-sub-filters with respect to different dynamic models and with weights determined by corresponding transition probabilities. 33 Therefore, to improve the navigation performance of KKRP with common BDS receiver in complex environments. A hierarchical filtering architecture of BDS/INS TC for KKRP algorithm is proposed in this research. Especially, aiming at the essential issue of carrier-based high-precision positioning, an INS-aided BDS singlefrequency AR strategy is developed. Moreover, the robust adaptive filtering AR model is exploited to enhance the AR success rate in BDS-challenging environment. Both simulation experiment with static data and field vehicle running test are conducted to verify the algorithm practicability.
The article is organized as follows.
where Dr denotes the symbol of DD operator; subscripts F, u and r represent the observed frequency (BDS B1) and the two moving stations, respectively; superscripts s and k are the observed satellites consisting DD. P, F, and D are DD code pseudorange, carrier phase (in distance), and Doppler observations, respectively. r and _ r are the geometric distance and its changing rate between station and satellite, respectively; l and N represent the carrier wavelength and initial integer ambiguity, repectively; I, T , M and _ I, _ T, _ M represent ionospheric delay, tropospheric delay, and multipath error and their changing rate, respectively. e P , e F , and e D are the measurement noise of pseudorange, carrier phase, and Doppler, respectively. For short baseline, the ionospheric and tropospheric delay in the DD are too small to be neglected. The multipath error can be modeled or processed within measurement noise. Since the position error of reference station has little effect on the baseline calculation, which is less than 1 cm with the reference station position error less than 10 m in the baseline of 20 km, the reference station can be selected with one of the moving station (called the moving reference station) with single-point pseudorange positioning. 2 where dx,dy,dz,dv x ,dv y ,dv z represent the increments of baseline and relative velocity components, respectively and DrN i (i = 1, 2, . . . , m) represent the DD ambiguities with m + 1 observed satellites. Discrete state equation is written as below where F k, kÀ1 is the one-step state transition matrix, which is defined by W kÀ1 is system noise vector, whose covariance matrix is Q kÀ1 . They are constructed as follows where D is the single difference operator, E is the direction matrix with unit sight vector between satellite and rover, and I is the identity matrix. It is expressed as follows 2   6  6  6  6  4   3   7  7  7  7 5 , Then, the discrete observation equation is written as where Z k is observation vector, H k is the observation matrix, and V k is the observation noise vector. Among them, the covariance matrix of V k is R k .
Here, the random function model based on the elevation angle is used for R k . The variance of un-differenced observation is calculated as 34 where el is satellite elevation angle and s 0 is the standard deviation factor, which is valued as 0.003 m, 0.3 m, and 0.01 m/s for carrier phase, code pseudorange, and Doppler measurement, respectively.

EKF fusion algorithm
An EKF fusion algorithm is adopted in this research for KKRP. It is a recursive solving procedure which contains two steps of prediction and update, which are written as follows. The prediction step iŝ whereX kÀ1, kÀ1 is the posteriori state vector at epoch k À 1 with covariance matrix of P kÀ1, kÀ1 ,X k, kÀ1 is the priori state vector at epoch k with covariance matrix P k, kÀ1 . The update step is where K k is the filtering gain matrix,X k, k is the posteriori state vector at epoch k with covariance matrix P k, k .

IRAFAR model for BDS single frequency
Carrier phase-based precise relative positioning depends on the fast and reliable resolution of integer ambiguity. However, the AR on-fly with single frequency is always a challenging question, especially in the dynamic complex environments such as signal blocked and severe multipath appeared as well as abnormal interference. INS is of short-term high-precision navigation performance, which can provide a prior information to reduce the ambiguity search space to improve the efficiency of AR. Therefore, an INS-aided robust adaptive filtering AR model is performed in this research to obtain high accuracy float solution of ambiguity, which is fixed with LAMBDA afterwards.

INS-aided AR model
In the observation equation (9), the geometric distance calculated by pseudorange single-point positioning is low precision, which is generally not enough for reliable AR, especially for single-frequency carrier. Therefore, the observation equations can be reconstructed by the predicted DD satellite-to-ground distance with INSderived position. Meanwhile, the virtual baseline increment from INS-derived position is used to enhance the model strength, which is written as Assuming that the priori covariance matrix P k, kÀ1 = 0 in EKF, the covariance matrix of weighted average solution of equation (15) with parameters of baseline and ambiguity, which are decoupled with X v can be given as where H bk = DE; and W P , W F , W INS are the weight matrices for the pseudorange, carrier phase, and INS virtual observation, respectively. Then, the variance matrix of ambiguity floats can be obtained as Compared with the variance matrix of ambiguity floats without INS-aided, the ambiguity search space is reduced with INS-aided, which improves the accuracy of the float ambiguities.

Robust adaptive filtering AR model
Although the AR reliability can be improved greatly by INS-aided in a good observation environment, it is also reduced significantly in the harsh surroundings including multipath gross error and abnormal interference. Hence, a robust adaptive AR model is performed to solve the problems. The robust factor and adaptive factor are designed, respectively, to adjust the weight between observation and prediction information. Namely, when there are gross errors in observation vector, robust factor is constructed to decrease the effects of outliers on the state estimation. On the contrary, when there are anomalous predictions in state vector, adaptive factor is constructed to reduce the effects of abnormal predicted states on the state estimation. Consequently, it is essential to design the robust factor and adaptive factor reasonably.
Design of robust factor. The observation noise covariance matrix is reconstructed by robust factor matrix to reduce the weight of observations with gross errors in the presence of outliers in observations as where r k is the robust factor matrix, whose elements are constructed by an improved IGG-III model 35 considering the correlation of DD observations as whereṽ i, k is the standardized residual statistic; c 0 ,c 1 are the threshold parameters, which are generally valued as c 0 = 1:0;2:5 and c 1 = 3:5;8:0, respectively.
Design of adaptive factor. The system covariance matrix is reconstructed by adaptive factor to decrease the weight of abnormal state predictions in the advent of improper disturbance as where g k is the adaptive factor matrix, whose elements are constructed by Huber's c function 36 as where c 2 is the regulation factor with a value of 0.85;3.0 usually; D e X k is state discrepancy statistic expressed as Afterwards, R k , P k, kÀ1 are used to update R k , P k, kÀ1 in filtering equation (14). When robust adaptive filtering is realized, the more accurate ambiguity floats are obtained. Then LAMBDA algorithm is adopted to fix ambiguities and ratio test is used to check the success rate to obtain the ''best'' ambiguity solution.

BDS/INS TC integrated model
For BDS/INS TC navigation system, the system model is derived from INS instrument, and the measurement model is composed of the difference between BDS DD observations and INS-derived counterparts. The integrated system is implemented by EKF with feedback correction to INS sensor biases.

State model of integrated system
System state space model depends on INS error model and system error description of inertial sensor. The integrated system model can be described as the following psi-angle error equations 37 without considering satellite system errors d _ r n = À dv n en 3 r n À v n en 3 dr n + dv n d _ v where, the subscripts i, e, n, and b represent inertial, Earth, navigation (east-north-up) and body frame, respectively; dr, dv, and F represent position, velocity, and attitude angle error vector, respectively; f n is specific force vector in n-frame; v n ie is the Earth rotation rate with respect to the inertial frame; v n en is the rotation rate of the navigation frame with respect to Earth; dv n ie and dv n en are the corresponding angular rate errors; C n b is the transformation matrix from body frame to navigation frame; dg n is the gravity error vector; $ b represents accelerometer error vector, which is modeled as first-order Gauss-Markov process; e b represents gyro drift error vector, which is modeled as random constant and first-order Gauss-Markov process. 38 Therefore, the system state equation for EKF is expressed as follows where F is state coefficient matrix; G is noise coefficient matrix; W is noise vector; and X is state vector with augmented gyro and accelerometer errors, which is written as e by e bz e rx e ry e rz r ax r ay r az T where, dr E , dr N , dr U ,dv E , dv N , dv U ,u E , u N , u U are the east-north-up coordinate components of position, velocity, and attitude angle error; e bx , e by , e bz , e rx , e ry , e rz are the three axis components of gyro drift error; r ax , r ay , r az are three axis components of the accelerometer error.

Measurement model of integrated system
Measurement model describes the relationship between the observations and the unknown parameters. In conventional TC measurement model, it is constituted by the differences between BDS DD code pseudorange, carrier phase, and Doppler observations and corresponding observations calculated by INS. As known, there are unknown parameters of ambiguity except for the navigation and INS sensor error parameters, which are augmented to state vector to be estimated. However, the number of ambiguity parameter is increased as the number of observed satellites, which makes the calculation too heavy, and even causes filtering divergence with cycle slip occurrence. In this research, a separated AR strategy is utilized in the BDS/INS TC to enhance the reliability of the integrated system, which is analyzed in section 2 and section 3. When the ambiguity parameters are solved reliably, the new measurement equations can be established with DD carrier phase and Doppler observations only, which are written as follows where v INS is the INS-predicted satellite-rover relative velocity, C e n is the transformation matrix from navigation frame to Earth frame, and other parameters are consistent with previous definitions. It is worth mentioning that the lever-arm correction from the inertial measurement unit (IMU) to the antenna phase center should be included. It can be compensated in observations by measurements with laser rangefinder or total station.
Write the observation equations into matrix form as follows where Z is the measurement vector, H is the measurement coefficient matrix, and V is the measurement noise matrix. An overview of the algorithm roadmap of the novel IRAFAR strategy of BDS/INS TC integrated system for KKRP is shown in Figure 1.

Simulation experiment with static data
To verify the IRAFAR algorithm proposed in the article, a set of static measured data is first used to perform static simulating dynamic for KKRP. The static BDS data were collected on the roof of the College of Automation Engineering Building of Nanjing University of Aeronautics and Astronautics on 6 September 2021. In the experiment, the BDS data with 1 Hz were collected by two NovAtel OEM-615 card navigation timing receivers. The IMU raw data were collected by MTI-G-700 sensor, and the data sampling rate was 200 Hz. The performance specifications of the IMU sensor are shown in Table 1. The data collection time was about 18 min, and the baseline length was 3.98 m. In the experiment, the B1 single-frequency data of BDS was used for algorithm analysis, and the position result output by Qianxun Cors RTK was used as a reference value. The measurement accuracy of pseudorange, carrier phase, and Doppler observations are 0.3 m, 0.003 m, and 0.01 m/s, respectively, and INS accuracy is 0.1 m. The satellite cut-off angle is set to 15°.
The data were processed as follows: (1) gross error of 15 m is added to the pseudorange observations randomly at epochs from 200 to 210 s, 600 to 610 s, and 1000 to 1010 s to simulate the outliers in observations; (2) velocity error of 5 m/s is added to the velocity predicted parameter at epochs from 500 to 510 s and 900 to 910 s to simulate the kinematic model abnormalities; (3) observation gross errors and velocity errors are added at epoch from 350 to 360 s and 750 to 760 s to simulate the two types of outliers exist simultaneously.
Three indicators of ratio value, ambiguity dilution of precision (ADOP) value, and empirical success rate [39][40][41] are defined as follows to evaluate the efficiency of AR, respectively where O sec and O min are the sub-optimal and optimal quadratic form of ambiguity residual determined in LAMBDA where QX a is the covariance matrix of ambiguity floats and m is the number of ambiguity P SE = number of correctly fixed epochs total number of epochs ð31Þ The following four schemes are performed for AR and KKRP, respectively.  Figure 3.
The ratio value and ADOP value of EKFAR and IEKFAR for original data is shown in Figures 4 and 5. The ratio test threshold is 3. Compared with EKFAR, the ratio value of IEKFAR algorithm is larger and ADOP value of is smaller. That is because the accuracy of ambiguity floats is improved with INS-aided, and  the search space is decreased. The ratio means of EKFAR and IEKFAR are 65.5 and 106.5, respectively. The success rates of the two algorithms are 98.7% and 99.6%, respectively. The baseline error is shown in Figure 6. The root mean square (RMS) errors of the two algorithms are 6 and 4.8 mm, which is improved by 20% with INS-aided.
Moreover, the processed data are utilized to test the performance of the EKFAR, IEKFAR, IRKFAR, and IRAFAR algorithms. As shown in Figure 7, the ratio value is decreased significantly at the epochs that are outliers in observations or abnormalities in model Ratio mean Figure 7. The ratio value for processed data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system-aided extended Kalman filtering ambiguity resolution; IRKFAR: inertial navigation system-aided robust extended Kalman filtering ambiguity resolution; IRAFAR: inertial navigation system-aided robust adaptive extended Kalman filtering ambiguity resolution.   Figure 8 are 94%, 96.5%, 98.1%, and 99.6%, respectively. The baseline errors are shown in Figure 9, and the errors statistics are shown in Table 2. The baseline errors are increased greatly at the epochs in the presence of outliers in observations and abnormalities in state parameters with EKFAR-BDRP and IEKFAR-TCRP algorithms. The baseline error is decreased with IRKFAR-TCRP algorithm, which can suppress the gross errors in the observations and improve the success rate in the presence of measurements outliers. However, the accuracy of IRAFAR-TCRP algorithm is improved significantly. The maximum error of the IRAFAR-TCRP algorithm is within 8 cm, and the RMS error is 6 mm, which is improved by 96% compared to EKFAR-BDRP. Success rate (%) Figure 8. The success rate for processed data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system-aided extended Kalman filtering ambiguity resolution; IRKFAR: inertial navigation system-aided robust extended Kalman filtering ambiguity resolution; IRAFAR: inertial navigation system-aided robust adaptive extended Kalman filtering ambiguity resolution. Baseline error (m) Figure 9. The baseline error for processed data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system-aided extended Kalman filtering ambiguity resolution; IRKFAR: inertial navigation system-aided robust extended Kalman filtering ambiguity resolution; IRAFAR: inertial navigation system-aided robust adaptive extended Kalman filtering ambiguity resolution; BDRP: BeiDou relative positioning; TCRP: tightly coupled relative position. of field vehicle is shown in Figure 10, and the trajectories of the three results are shown in Figure 11. The skyplot of visible satellites is shown in Figure  12. It is worth mentioning that there are no obvious outliers in observation data and no significant kinematic disturbances during the running. As shown in Figure 13, the visible satellite number is more than 10 and the PDOP value is less than 3. The baseline length and relative velocity are shown in Figure 14    Relative velocity (m/s) Figure 14. The baseline length and relative velocity.

Practical experiments with field vehicle data
To verify the performance of IRAFAR-TCRP algorithm further, the data are processed as follows: (1) gross error of 15 m is added to the pseudorange observations randomly at epochs from 700 to 710 s and 1400 to 1410 s to simulate the outliers in observations; (2) velocity error of 5 m/s is added to the velocity predicted parameter at epochs from 500 to 510 s and 1000 to 1010 s to simulate the kinematic model abnormalities; (3) satellite observation outage is utilized at epoch from 800 to 805 s, 1600 to 1610 s, and 2000 to 2020 s to simulate satellite signal blocked.
As shown in Figure 18, the ratio value is decreased largely at the epochs in the presence of outliers with EKFAR and IEKFAR algorithms. The IRKFAR algorithm can eliminate the gross error in the observations at 700;710 s and 1400;1410 s rather than the model abnormalities at the epochs of 500;510s and 1000;1010 s. The IRAFAR algorithm can eliminate both outliers in observations and model-predicted states. The ratio means of the four algorithms are 18.5, 24.4, 26.9, and 27.7, respectively. As shown in Figure  19, the success rates of the four algorithms are 82.4%, 86.6%, 88.4%, and 89.4%, respectively. The baseline errors are shown in Figure 20. It can be seen that the baseline errors are increased significantly at the epochs in the presence of the outliers. However, the baseline error of the IRAFAR-TCRP algorithm can keep a high accuracy in the presence of both measurement outliers and model abnormalities in the continuous period of satellite signal that the BDS/INS TC can remain high accuracy INS-derived position. Therefore, centimeter precision level can be achieved. However, the baseline error is increased rapidly with the outage duration of BDS signal. That is because the position information is provided by the INS mechanization alone during the satellite signal outage period, and the position error drift is mainly dependent on the quality of IMU sensors. In this field test, meter position precision can be obtained within 10 s outage duration, and the AR recovery is within 2 s when the BDS signal is recovery. The position error of the BDS/INS TC system is shown in Figure 21  Ratio mean Figure 18. The ratio value for processed data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system-aided extended Kalman filtering ambiguity resolution; IRKFAR: inertial navigation system-aided robust extended Kalman filtering ambiguity resolution; IRAFAR: inertial navigation system-aided robust adaptive extended Kalman filtering ambiguity resolution. Success rate (%) Figure 19. The success rate for processed data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system-aided extended Kalman filtering ambiguity resolution; IRKFAR: inertial navigation system-aided robust extended Kalman filtering ambiguity resolution; IRAFAR: inertial navigation system-aided robust adaptive extended Kalman filtering ambiguity resolution. 3. Compared with conventional EKFAR-BDRP algorithm, IRAFAR-TCRP algorithm can improve the accuracy of KKRP greatly. In particular, in the BDS-challenging environments, IRAFAR-TCRP not only increase the precision of relative positioning but also enhance the continuity of positioning when the GNSS signal is blocked.
With the rapid development of GNSS system, multifrequency multi-constellation can obtain more observations and virtual combination frequencies, which will provide a new idea for AR. Therefore, our future work will focus on the AR and precise relative positioning with multi-frequency multi-constellation GNSS/INS TC system under an actual urban GNSS-challenging environment.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The work described in this article was partially supported by the National Natural Science Foundation of China (grant nos. 61973160 and 61873125) and Funding of Jiangsu Innovation Program for Graduate Education (grant no. KYLX15_0276).