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

Abstract

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.

Keywords

Ambiguity resolution robust adaptive filtering tightly coupled BeiDou navigation satellite system/inertial navigation system relative positioning

Introduction

Traditional global navigation satellite system (GNSS) real-time kinematic (RTK) technology using differential carrier phase observations can achieve centimeter-level 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 vehicle (UAV) aerial refueling, carrier-based aircraft 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–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–Wü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–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–10 However, it remains a challenging problem for low-cost 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–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–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;¹⁸ 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–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–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.³⁰ 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.³³

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 single-frequency 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. The section “The EKF-based BDS KKRP model” presents the EKF-based BDS KKRP model. Then, a novel IRAFAR model for BDS single frequency is elaborated in the section “IRAFAR model for BDS single frequency.” Next, the implementation of the BDS/INS TC integrated system is described in the section “BDS/INS TC integrated model.” In the section “Experiment and discussion,” algorithm feasibility is tested and results are discussed in detail. Finally, the section “Conclusion” presents our conclusions and future work.

The EKF-based BDS KKRP model

In general, the code pseudorange, carrier phase and Doppler observations are used in the single-frequency single-epoch precise KKRP. The observation equations of double-difference (DD) code pseudorange, carrier phase, and Doppler can be written as follows

\begin{matrix} Δ \nabla P_{F, ur}^{sk} = Δ \nabla ρ_{ur}^{sk} + Δ \nabla I_{F, ur}^{sk} + Δ \nabla T_{ur}^{sk} \\ + Δ \nabla M_{P_{F}, ur}^{sk} + ε_{Δ \nabla P_{F}, ur}^{sk} \end{matrix}

(1)

\begin{matrix} Δ \nabla Φ_{F, ur}^{sk} = Δ \nabla ρ_{ur}^{sk} + λ_{F} Δ \nabla N_{F, ur}^{sk} \\ - Δ \nabla I_{F, ur}^{sk} + Δ \nabla T_{ur}^{sk} + Δ \nabla M_{Φ_{F}, ur}^{sk} + ε_{Δ \nabla Φ_{F}, ur}^{sk} \end{matrix}

(2)

\begin{matrix} λ_{F} Δ \nabla D_{F, ur}^{sk} = Δ \nabla {\overset{\cdot}{ρ}}_{ur}^{sk} + Δ \nabla {\overset{\cdot}{I}}_{F, ur}^{sk} + Δ \nabla {\overset{\cdot}{T}}_{ur}^{sk} \\ + Δ \nabla {\overset{\cdot}{M}}_{D_{F}, ur}^{sk} + ε_{Δ \nabla D_{F}, ur}^{sk} \end{matrix}

(3)

where $Δ \nabla$ 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$ , $Φ$ , and $D$ are DD code pseudorange, carrier phase (in distance), and Doppler observations, respectively. $ρ$ and $\overset{\cdot}{ρ}$ are the geometric distance and its changing rate between station and satellite, respectively; $λ$ and $N$ represent the carrier wavelength and initial integer ambiguity, repectively; $I$ , $T$ , $M$ and $\overset{\cdot}{I}$ , $\overset{\cdot}{T}$ , $\overset{\cdot}{M}$ represent ionospheric delay, tropospheric delay, and multipath error and their changing rate, respectively. $ε_{P}$ , $ε_{Φ}$ , and $ε_{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,10

State equation

An EKF method is used for parameter estimation of KKRP. The system model is constructed by BDS with constant velocity (CV) kinematic model and ambiguity parameters are modeled as a random walk process. The state parameter vector is defined by

X = {[\underset{X_{b}}{\underset{︸}{δ x, δ y, δ z,}} \underset{X_{v}}{\underset{︸}{δ v_{x}, δ v_{y}, δ v_{z},}} \underset{Δ \nabla N}{\underset{︸}{Δ \nabla N_{1}, Δ \nabla N_{2}, \dots, Δ \nabla N_{m}}}]}^{T}

(4)

where $δ x$ , $δ y$ , $δ z$ , $δ v_{x}$ , $δ v_{y}$ , $δ v_{z}$ represent the increments of baseline and relative velocity components, respectively and $Δ \nabla N_{i} (i = 1, 2, \dots, m)$ represent the DD ambiguities with $m + 1$ observed satellites.

Discrete state equation is written as below

X_{k} = Φ_{k, k - 1} X_{k - 1} + W_{k - 1}

(5)

where $Φ_{k, k - 1}$ is the one-step state transition matrix, which is defined by

Φ_{k, k - 1} = [\begin{matrix} F_{6 \times 6} & 0_{6 \times m} \\ 0_{m \times 6} & I_{m \times m} \end{matrix}], F_{6 \times 6} = [\begin{matrix} I_{3 \times 3} & T I_{3 \times 3} \\ 0_{3 \times 3} & I_{3 \times 3} \end{matrix}]

(6)

$W_{k - 1}$ is system noise vector, whose covariance matrix is $Q_{k - 1}$ . They are constructed as follows

W_{k - 1} = {[0, 0, 0, 0, 0, 0, ω_{Δ \nabla N_{1}}, \dots, ω_{Δ \nabla N_{m}}]}^{T}

(7)

Q_{k - 1} = diag (0, 0, 0, 0, 0, 0, T σ_{Δ \nabla N_{1}}^{2}, \dots, T σ_{Δ \nabla N_{m}}^{2})

(8)

where $σ_{Δ \nabla N_{i}} = σ_{Δ \nabla N} (i = 1, \dots, m)$ are the standard variance of DD ambiguities process noise and $T$ is the time interval.

Observation equation

The observation equation is obtained by linearizing equations (1)–(3), then it can be expressed in matrix form as follows

\begin{matrix} [\begin{matrix} Δ \nabla P - Δ \nabla ρ \\ Δ \nabla Φ - Δ \nabla ρ \\ λ Δ \nabla D - Δ \nabla \overset{\cdot}{ρ} \end{matrix}] = [\begin{matrix} DE & 0_{m \times 3} & 0_{m \times m} \\ DE & 0_{m \times 3} & λ I_{m \times m} \\ 0_{m \times 3} & DE & 0_{m \times m} \end{matrix}] \\ [\begin{matrix} X_{b} \\ X_{v} \\ Δ \nabla N \end{matrix}] + [\begin{matrix} ε_{Δ \nabla P} \\ ε_{Δ \nabla ϕ} \\ ε_{Δ \nabla D} \end{matrix}] \end{matrix}

(9)

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

\begin{matrix} D = [\begin{matrix} 1 & - 1 & 0 & \dots & 0 \\ 1 & 0 & - 1 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & 0 & 0 & \dots & - 1 \end{matrix}], \\ E = {(1_{u}^{1}, 1_{u}^{2}, \dots, 1_{u}^{m + 1})}^{T} \end{matrix}

(10)

Then, the discrete observation equation is written as

Z_{k} = H_{k} X_{k - 1} + V_{k}

(11)

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³⁴

σ^{2} = σ_{0}^{2} + σ_{0}^{2} / \sin^{2} (el)

(12)

where $el$ is satellite elevation angle and $σ_{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 is

\begin{matrix} {\hat{X}}_{k, k - 1} = Φ_{k, k - 1} {\hat{X}}_{k - 1, k - 1} \\ P_{k, k - 1} = Φ_{k, k - 1} P_{k - 1, k - 1} Φ_{k, k - 1}^{T} + Q_{k - 1} \end{matrix}

(13)

where ${\hat{X}}_{k - 1, k - 1}$ is the posteriori state vector at epoch $k - 1$ with covariance matrix of $P_{k - 1, k - 1}$ , ${\hat{X}}_{k, k - 1}$ is the priori state vector at epoch $k$ with covariance matrix $P_{k, k - 1}$ .

The update step is

\begin{matrix} K_{k} = P_{k, k - 1} H_{k}^{T} {(H_{k} P_{k, k - 1} H_{k}^{T} + R_{k})}^{- 1} \\ {\hat{X}}_{k, k} = {\hat{X}}_{k, k - 1} + K_{k} (Z_{k} - H_{k} {\hat{X}}_{k, k - 1}) \\ P_{k, k} = (I - K_{k} H_{k}) P_{k, k - 1} \end{matrix}

(14)

where $K_{k}$ is the filtering gain matrix, ${\hat{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 INS-derived position. Meanwhile, the virtual baseline increment from INS-derived position is used to enhance the model strength, which is written as

\begin{matrix} [\begin{matrix} Δ \nabla P - Δ \nabla P_{INS} \\ Δ \nabla Φ - Δ \nabla P_{INS} \\ 0_{3 \times 1} \end{matrix}] = [\begin{matrix} D E & 0_{m \times 3} 0_{m \times m} \\ D E & 0_{m \times 3} λ I_{m \times m} \\ I_{3 \times 3} & 0_{3 \times 3} 0_{3 \times m} \end{matrix}] \\ [\begin{matrix} X_{b} \\ X_{v} \\ Δ \nabla N \end{matrix}] + [\begin{matrix} ε_{Δ \nabla P} \\ ε_{Δ \nabla ϕ} \\ ε_{INS} \end{matrix}] \end{matrix}

(15)

where $Δ \nabla P_{INS}$ is the predicted DD satellite-to-ground distance by INS and $ε_{INS}$ is the observation noise of INS system obtained by the TC filtering prior equation.

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

\begin{matrix} P_{k, k} = {(H_{k}^{T} R_{k}^{- 1} H_{k})}^{- 1} = \\ {(\begin{matrix} H_{bk}^{T} W_{P} H_{bk} + H_{bk}^{T} W_{Φ} H_{bk} + W_{INS} & λ H_{bk}^{T} W_{Φ} \\ λ W_{Φ} H_{bk} & λ^{2} W_{Φ} \end{matrix})}^{- 1} \end{matrix}

(16)

where $H_{bk} = DE$ ; and $W_{P}$ , $W_{Φ}$ , $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

P_{k, k, \hat{N}}^{*} = λ^{- 2} [W_{Φ}^{- 1} + H_{bk} {(H_{bk}^{T} W_{P} H_{bk} + W_{INS})}^{- 1} H_{bk}^{T}]

(17)

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

{\bar{R}}_{k} = \frac{R_{k}}{r_{k}}

(18)

where $r_{k}$ is the robust factor matrix, whose elements are constructed by an improved IGG-III model³⁵ considering the correlation of DD observations as

r_{ij, k} = \sqrt{r_{ii, k} r_{jj, k}}

(19)

with

r_{i i, k} = {\begin{matrix} 1 & | {\tilde{v}}_{i, k} | \leq c_{0} \\ \frac{c_{0}}{| {\tilde{v}}_{i, k} |} {(\frac{c - | {\tilde{v}}_{i, k} |}{c_{1} - c_{0}})}^{2} & c_{0} < | {\tilde{v}}_{i, k} | \leq c_{1} \\ 10^{- 12} & | {\tilde{v}}_{i, k} | > c_{1} \end{matrix}

(20)

where ${\tilde{v}}_{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

{\bar{P}}_{k, k - 1} = \frac{P_{k, k - 1}}{γ_{k}}

(21)

where $γ_{k}$ is the adaptive factor matrix, whose elements are constructed by Huber’s $ψ$ function³⁶ as

γ_{k} = {\begin{matrix} 1 | Δ {\tilde{X}}_{k} | \leq c_{2} \\ \frac{c_{2}}{| Δ {\tilde{X}}_{k} |} | Δ {\tilde{X}}_{k} | > c_{2} \end{matrix}

(22)

where $c_{2}$ is the regulation factor with a value of 0.85∼3.0 usually; $Δ {\tilde{X}}_{k}$ is state discrepancy statistic expressed as

Δ {\tilde{X}}_{k} = \frac{‖ {\hat{X}}_{k, k} - {\hat{X}}_{k, k - 1} ‖}{\sqrt{tr (P_{k, k - 1})}}

(23)

Afterwards, ${\bar{R}}_{k}$ , ${\bar{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³⁷ without considering satellite system errors

{\begin{matrix} δ {\overset{\cdot}{r}}^{n} = - δ ω_{en}^{n} \times r^{n} - ω_{en}^{n} \times δ r^{n} + δ v^{n} \\ δ {\overset{\cdot}{v}}^{n} = f^{n} \times Φ - (2 δ ω_{ie}^{n} + δ ω_{en}^{n}) \times v^{n} - (2 ω_{ie}^{n} + ω_{en}^{n}) \times δ v^{n} + C_{b}^{n} \nabla^{b} + δ g^{n} \\ \overset{\cdot}{Φ} = Φ \times (ω_{ie}^{n} + ω_{en}^{n}) + (δ ω_{ie}^{n} + δ ω_{en}^{n}) - C_{b}^{n} ε^{b} \end{matrix}

(24)

where, the subscripts i, e, n, and b represent inertial, Earth, navigation (east–north–up) and body frame, respectively; $δ r$ , $δ v$ , and $Φ$ represent position, velocity, and attitude angle error vector, respectively; $f^{n}$ is specific force vector in n-frame; $ω_{ie}^{n}$ is the Earth rotation rate with respect to the inertial frame; $ω_{en}^{n}$ is the rotation rate of the navigation frame with respect to Earth; $δ ω_{ie}^{n}$ and $δ ω_{en}^{n}$ are the corresponding angular rate errors; $C_{b}^{n}$ is the transformation matrix from body frame to navigation frame; $δ g^{n}$ is the gravity error vector; $\nabla^{b}$ represents accelerometer error vector, which is modeled as first-order Gauss–Markov process; $ε^{b}$ represents gyro drift error vector, which is modeled as random constant and first-order Gauss–Markov process.³⁸

Therefore, the system state equation for EKF is expressed as follows

\overset{\cdot}{X} (t) = F (t) X (t) + G (t) W (t)

(25)

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

\begin{array}{l} X = [δ r_{E} δ r_{N} δ r_{U} δ v_{E} δ v_{N} δ v_{U} φ_{E} φ_{N} φ_{U} \\ ε_{b x} ε_{b y} ε_{b z} ε_{r x} ε_{r y} ε_{r z} \nabla_{a x} \nabla_{a y} \nabla_{a z}]^{T} \end{array}

(26)

where, $δ r_{E}, δ r_{N}, δ r_{U}$ , $δ v_{E}, δ v_{N}, δ v_{U}$ , $φ_{E}, φ_{N}, φ_{U}$ are the east–north–up coordinate components of position, velocity, and attitude angle error; $ε_{bx}, ε_{by}, ε_{bz}, ε_{rx}, ε_{ry}, ε_{rz}$ are the three axis components of gyro drift error; $\nabla_{ax}, \nabla_{ay}, \nabla_{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

\begin{matrix} Δ \nabla Φ - Δ \nabla P_{INS} - λ Δ \nabla N = {DEC}_{n}^{e} δ r + ε_{Δ \nabla_{Φ}} \\ λ \cdot Δ \nabla D - Δ \nabla v_{INS} = {DEC}_{n}^{e} δ v + ε_{Δ \nabla D} \end{matrix}

(27)

where $v_{INS}$ is the INS-predicted satellite-rover relative velocity, $C_{n}^{e}$ 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

Z (t) = H (t) X (t) + V (t)

(28)

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.

Figure 1.

The algorithm roadmap of the novel IRAFAR strategy of BDS/INS tightly coupled for KKRP. BDS: BeiDou navigation satellite system; DD: double difference; IRAFAR: INS-aided robust adaptive filtering ambiguity resolution; KKRP: kinematic-to-kinematic relative positioning; INS: inertial navigation system.

Experiment and discussion

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

Table 1.

The performance specifications of the IMU sensor.

IMU	Bias	Scale factor	Random walk
Gyroscopes	10°/h	1500 ppm	0.2°/ $\sqrt{h}$
Accelerometers	40 µg	4000 ppm	0.1 m/s/ $\sqrt{h}$

IMU: inertial measurement unit.

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–41 are defined as follows to evaluate the efficiency of AR, respectively

Ratio = \frac{Ω_{\sec}}{Ω_{min}}

(29)

where $Ω_{\sec}$ and $Ω_{min}$ are the sub-optimal and optimal quadratic form of ambiguity residual determined in LAMBDA

ADOP = {(\sqrt{det (Q_{{\hat{X}}_{a}})})}^{\frac{1}{m}}

(30)

where $Q_{{\hat{X}}_{a}}$ is the covariance matrix of ambiguity floats and $m$ is the number of ambiguity

P_{SE} = \frac{number of correctly fixed epochs}{total number of epochs}

(31)

The following four schemes are performed for AR and KKRP, respectively.

Scheme 1. Extended Kalman filtering method for AR and carrier-based BDS KKRP (EKFAR-BDRP).

Scheme 2. INS-aided extended Kalman filtering method for AR and carrier-based BDS/INS TC KKRP (IEKFAR-TCRP).

Scheme 3. INS-aided robust Kalman filtering method for AR and carrier-based BDS/INS TC KKRP (IRKFAR-TCRP).

Scheme 4. INS-aided robust adaptive filtering method for AR and carrier-based BDS/INS TC KKRP (IRAFAR-TCRP).

In Schemes 2, 3, and 4, the test thresholds are taken as empirical value $c_{0} = 1.0$ , $c_{1} = 3.5$ , $c_{2} = 1.5$ . The skyplot of visible satellites are shown in Figure 2. The number of visible satellites and the position dilution of precision (PDOP) value are shown in Figure 3.

Figure 2.

The skyplot of visible satellites.

Figure 3.

The number of visible satellites and PDOP value. PDOP: position dilution of precision.

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.

Figure 4.

The ratio value for original data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system–aided extended Kalman filtering ambiguity resolution.

Figure 5.

The ADOP value for original data. ADOP: ambiguity dilution of precision; EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system–aided extended Kalman filtering ambiguity resolution.

Figure 6.

The baseline error for original data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system–aided extended Kalman filtering ambiguity resolution; TCRP: tightly coupled relative position.

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 predictions with EKFAR and IEKFAR. However, the IRKFAR algorithm can eliminate the gross error completely at the epochs 200∼210 s, 600∼610 s and 1000∼1010 s, but it cannot decrease the abnormalities in model predictions at the epochs 350∼360 s, 500∼510 s, 750∼760 s, and 900∼910 s. Ulteriorly, the IRAFAR algorithm can fully eliminate both the measurement gross error and model outliers at all the epochs. The ratio value keeps a high level in the whole period. The ratio means of the four algorithms are 52.8, 87.5, 91.7, and 99.4, respectively. The success rates of the four algorithms as shown in 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.

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.

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.

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.

Table 2.

The errors statistics of different algorithms.

Schemes	Baseline errors (m)
	Maximum	Minimum	Mean	RMS
EKFAR-BDRP	0.88	−1.564	0.019	0.156
IEKFAR-TCRP	0.867	−1.46	0.018	0.141
IRKFAR-TCRP	0.658	−0.105	0.003	0.042
IRAFAR-TCRP	0.011	−0.076	0.0003	0.006

RMS: root mean square; 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.

Practical experiments with field vehicle data

To test the actual performance of IRAFAR-TCRP algorithm further, the practical experiments with field vehicle data are conducted. The field vehicle test was run on 5 September 2021 by two cars with the same route, and a NovAtel OEM-615 card navigation timing receiver and a Qianxun Cors RTK module were installed on the top of each car with a sampling frequency of 1 Hz. Moreover, an MTI-G-700 INS navigation sensor was installed in rover station, which is an IMU-enhanced GPS/INS integrated system. The raw IMU data sampling rate was 200 Hz. The data collection time was about 43 min, and the baseline length was less than 200 m. The position result output by Qianxun Cors RTK was used as a reference value. The satellite cut-off angle is set to 15°. The experiment scene of field vehicle is shown in Figure 10, and the trajectories of the three results are shown in Figure 11.

Figure 10.

The experiment scene of field vehicle.

Figure 11.

The trajectories of the three results. BDS: BeiDou navigation satellite system; TC: tightly coupled.

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. The ratio value and ADOP value for original data are shown in Figures 15 and 16. The ratio means of EKFAR and IEKFAR are 21.1 and 27.5, respecitively. The success rates of the two algorithm are 85.8% and 90.2%. The baseline errors are shown in Figure 17. The RMS errors of the two algorithms are 4.3 and 2.4 cm, which is improved by 44.2% with INS-aided.

Figure 12.

The skyplot of visible satellites.

Figure 13.

The number of visible satellites and PDOP value. PDOP: position dilution of precision.

Figure 14.

The baseline length and relative velocity.

Figure 15.

The ratio value for original data. EKFAR: extended Kalman filtering ambiguity resolution; IEKFAR: inertial navigation system–aided extended Kalman filtering ambiguity resolution.

Figure 16.

Figure 17.

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. The RMS errors of the BDS/INS TC system are within 0.2 m in the horizontal direction and 0.6 m in the up direction.

Figure 18.

Figure 19.

Figure 20.

The baseline error of 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.

Figure 21.

Position error of tightly coupled system.

Conclusion

In this article, a novel AR model of BDS/INS TC for KKRP is investigated. Especially, the INS-aided robust adaptive filtering AR algorithm is elaborated. A hierarchical filtering strategy is developed to utilize the “best” ambiguity parameters in BDS/INS TC integrated system. Both static experiment and field vehicle test were carried out to validate the AR and relative positioning performance under different data configurations. According to the theoretical analysis and experimental results, the following favorable conclusions can be drawn:

Compared with EKFAR algorithm, IEKFAR algorithm can decrease the ambiguity search space and improve the accuracy of ambiguity floats effectively so as to raise the success rate of AR. The ratio values obviously increase and ADOP values decrease.

Compared with EKFAR and IEKFAR algorithm, IRKFAR algorithm can suppress the influence of measurement gross error rather than model prediction outliers. However, IRAFAR algorithm can reduce the influence of observation outliers and model prediction abnormalities simultaneously.

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, multi-frequency 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.

Footnotes

Handling Editor: Lyudmila Mihaylova

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

ORCID iDs

Ling Wu

Yongrong Sun

Xiyu Fu

References

Chen

Shu

Ding

, et al. Kinematic single-frequency relative positioning for LEO formation flying mission. GPS Solut 2015; 19: 525–535.

Shao

. Performance analysis of BDS/GPS single epoch short baseline motion-to-motion positioning. Sci Surv Map 2021; 46(1): 1–6.

Duan

Zhao

Pang

, et al. A novel GPS precise relative positioning algorithm for moving base station. J Sichuan Univ Eng Sci Ed 2015; 47(3): 130–136.

Zhao

Pang

, et al. New method for Beidou dual-frequency kinematic-to-kinematic relative positioning. J Chin Inert Technol 2016; 24(6): 758–762.

Gao

Pan

, et al. Improving ambiguity resolution for medium baselines using combined GPS and BDS dual/triple-frequency observations. Sensors 2015; 15(11): 27525–27542.

Huang

Zhai

, et al. A new triple-frequency cycle slip detecting algorithm validated with BDS data. GPS Solut 2016; 20(4): 761–769.

Lou

Gong

, et al. Assessment of code bias variations of BDS triple-frequency signals and their impacts on ambiguity resolution for long baselines. GPS Solut 2017; 21(1): 177–186.

Dai

, et al. Accuracy and reliability of multi-GNSS real-time precise positioning: GPS, GLONASS, BeiDou, and Galileo. J Geod 2015; 89(6): 607–635.

Odolinski

Teunissen

PJG

Odijk

. Combined BDS, Galileo, QZSS and GPS single-frequency RTK. GPS Solut 2015; 19(1): 151–163.

10.

Yang

, et al. Performance assessment of single- and dual-frequency BeiDou/GPS single-epoch kinematic positioning. GPS Solut 2014; 18: 393–403.

11.

Yang

Liu

. An improved relative GNSS tracking method utilizing single frequency receivers. Sensors 2020; 20(15): 4073.

12.

Tsakiri

Sioulis

Piniotis

. The use of low-cost, single-frequency GNSS receivers in mapping surveys. Surv Rev 2018; 50(358): 46–56.

13.

, et al. Performance analysis of positioning solution using low-cost single-frequency u-blox receiver based on baseline length constraint. Sensors 2019; 19(19): 4352.

14.

Shi

Wang

, et al. Effect analysis of GNSS/INS processing strategy for sufficient utilization of urban environment observations. Sensors 2021; 21(2): 620.

15.

Wagner

. GNSS/INS integration: still an attractive candidate for automatic landing systems? GPS Solut 2005; 9: 179–193.

16.

Petovello

Cannon

Lachapelle

. Benefits of using a tactical-grade IMU for high-accuracy positioning. Navigation 2004; 51: 1–12.

17.

Grejner-Brzezinska

Toth

. GPS error modeling and OTF ambiguity resolution for high-accuracy GPS/INS integrated system. J Geod 1998; 72(11): 626–638.

18.

Falco

Pini

Marucco

. Loose and tight GNSS/INS integrations: comparison of performance assessed in real urban scenarios. Sensors 2017; 17: 255.

19.

Zhang

Niu

, et al. Tightly-coupled integration of multi-GNSS single-frequency RTK and MEMS-IMU for enhanced positioning performance. Sensors 2017; 17(11): 2462.

20.

Hinüber

ELV

Reimer

Schneider

, et al. INS/GNSS integration for aerobatic flight applications and aircraft motion surveying. Sensors 2017; 17(5): 941.

21.

Lan

Elsheikh

Abdelfatah

, et al. Integrated RTK/INS navigation for precision agriculture. In: Proceedings of the 32nd international technical meeting of the satellite division of the institute of navigation, Miami, FL, 16–20 September 2019. Manassas, VA: Institute of Navigation.

22.

Cox

Jr Brading

JDW

. Integration of LAMBDA ambiguity resolution with Kalman filter for relative navigation of spacecraft. Navigation 2000; 47: 205–210.

23.

Chang

Yang

Zhou

. MLAMBDA: a modified LAMBDA method for integer least-squares estimation. J Geod 2005; 79(9): 552–565.

24.

Chen

Zheng

Ramanandan

, et al. INS aided GPS integer ambiguity resolution. In: Proceedings of the IEEE international conference on control applications, Denver, CO, 28–30 September 2011, pp.567–572. New York: IEEE.

25.

Han

Wang

, et al. Reliable partial ambiguity resolution for single-frequency GPS/BDS and INS integration. GPS Solut 2017; 21(1): 251–264.

26.

Zhang

Ban

Niu

, et al. Improving the design of MEMS INS-aided PLLs for GNSS carrier phase measurement under high dynamics. Micromachines 2017; 8: 135.

27.

Cui

, et al. A tightly coupled RTK/INS algorithm with ambiguity resolution in the position domain for ground vehicles in harsh urban environments. Sensors 2018; 18: 2160.

28.

Dong

Wang

Zhang

, et al. Tightly coupled GNSS/INS integration with robust sequential Kalman filter for accurate vehicular navigation. Sensors 2020; 20: 561.

29.

Medina

Vilà-Valls

, et al. Robust filtering techniques for RTK positioning in harsh propagation environments. Sensors 2021; 21(4): 1250.

30.

Han

Wang

, et al. Single-frequency GPS/BDS RTK and INS ambiguity resolution and positioning performance enhanced with positional polynomial fitting constraint. Rem Sens 2020; 12(15): 2374.

31.

Yang

. Adaptively robust filtering for kinematic geodetic positioning. J Geod 2001; 75(2–3): 109–116.

32.

Zhang

Zhao

. A two-step robust adaptive filtering algorithm for GNSS kinematic precise point positioning. Chin J Aeronaut 2021; 34: 210–219.

33.

Zhang

Guo

Zhang

. Data fusion based on adaptive interacting multiple model for GPS/INS integrated navigation system. Appl Sci 2018; 8(9): 1682.

34.

Zhang

Zhao

Zhou

. A novel weighting approach for variance component estimation in GPS/BDS PPP. IEEE Sens J 2019; 19(10): 3763–3771.

35.

Yang

Song

. Robust estimator for correlated observations based on bifactor equivalent weights. J Geod 2002; 76(6–7): 353–358.

36.

Yang

Cui

. Adaptively robust filter with multi adaptive factors. Surv Rev 2008; 40(309): 260–270.

37.

Savage

. Strapdown inertial navigation integration algorithm design part 2: velocity and position algorithms. J Guid Control Dyn 1998; 21(2): 208–221.

38.

Wang

Xiong

Liu

, et al. A new tightly-coupled INS/CNS integrated navigation algorithm with weighted multi-stars observations. Proc IMechE, Part G: J Aerospace Engineering 2015; 230(4): 698–712.

39.

Wang

Feng

Guo

. Reliability control of single-epoch RTK ambiguity resolution. GPS Solut 2017; 21(2): 591–604.

40.

Teunissen

PJG

. Success probability of integer GPS ambiguity rounding and bootstrapping. J Geod 1998; 72(10): 606–612.

41.

Liu

Zhang

, et al. An extended ADOP for performance evaluation of single-frequency single-epoch positioning by BDS/GPS in Asia-Pacific region. Sensors 2017; 17(10): 2254.