2.1 Multi-GNSS real-time PPP model

Similar to the RTK algorithm, PPP utilises both code and carrier phase measurements in the filtering process, and also estimates the integer ambiguity. Differently, PPP eliminates measurement errors by corrections, precise models and estimation, while RTK eliminates them based on the double-difference between satellites and receivers. The observation equations for the un-differenced GNSS pseudorange and carrier phase can be expressed as follows:

(1) \begin{equation}P_{r, j}^{s}=\rho_{r}^{s}+c\left(d t_{r}-d T^{s}\right)+I_{r, j}^{s}+T_{r}^{s}+D_{r, P_{1}}-D_{P_{1}}^{s}+\varepsilon_{P},\end{equation}

(2) \begin{equation}\Phi_{r, i}^{s}=\rho_{r}^{s}+c\left(d t_{r}-d T^{s}\right)-I_{r, i}^{s}+T_{r}^{s}+\lambda_{i} B_{r, i}^{s}+d \Phi_{r, i}^{s}+\varepsilon_{\Phi, i},\end{equation}

where

(3) \begin{equation}B_{r, i}^{s}=\phi_{r, 0, i}-\phi_{0, i}^{s}+N_{r, i}^{s},\end{equation}

(4) \begin{equation}d \Phi_{r, i}^{s}=-d_{r, p c o, i}^{T} e_{r, e m u}^{s}+\left(E^{s} d_{p o, i}^{s}\right)^{T} e_{r}^{s}+d_{r, p c v, i}(E l)+d_{p c v, i}^{s}(\theta)-d_{r, d i s p}^{T} e_{r, e n u}^{s}+\lambda_{s} \phi_{p w},\end{equation}

and i and j represent two different frequencies, $ P^s_{r,i} $ is the pseudorange between satellite s and receiver r, and $ \Phi_{r, i}^{s}=\lambda_{t} \phi_{r, i}^{s} $ is the phase range, where $\phi_{r, i}^{s}$ is the carrier phase measurement. $\rho_{r}^{s}=\left|p_{r}-p^{s}\right|$ is the geometrical distance between the receiver antenna mass centre and satellite mass centre, c is the speed of light in vacuum, $d t_{r}$ and $d T^{s}$ represent, respectively, the receiver clock offset at the signal receiving time and the satellite clock offset at the signal transmitting time, $I_{r, i}^{s}$ is the ionospheric delay along the signal propagation path at frequency i, $T_{r}^{s}$ is the tropospheric delay of the signal path, $D_{r, P_i}$ and $ D^s_{P_i} $ are the differential code biases (DCB) for the receiver and satellite, $ \lambda_i $ is the wavelength at frequency i, $B_{r, i}^{s}$ indicates the non-integer phase ambiguity, which contains the integer phase ambiguity $N_{r, i}^{s}$ and the un-calibrated phase delay (UPD) $\phi_{r, 0, i}$ and $\phi_{0, i}^{s}$ , where $\phi_{r, 0, i}$ is the initial phase (cycle) of the receiver local oscillator and $\phi_{0, i}^{s}$ is the initial phase of the transmitted navigation signal at initial time. $d \Phi_{r, i}^{s}$ represents the carrier-phase correction terms, including the antenna phase centre offsets (PCO) $d_{r, pco, i}^{T}$ , antenna phase centre variations (PCV) $d_{r, pcv, i}$ , station displacement by earth tides $d_{r, {disp}}^{T}$ and phase windup effect $\phi_{p w}$ , and $e_{r, e n u}^{s}$ is the vector from satellite s to receiver r in the East-North-Up (ENU) frame. $ E^s $ is the coordinate transformation matrix from the satellite body-fixed frame to the Earth-Centred Earth-Fixed (ECEF) frame, $e_{r}^{s}$ is the vector from satellite to receiver in the ECEF frame, El is the satellite elevation angle, ad $ \theta $ is the satellite nadir angle. Finally, the symbols $\varepsilon_{p}$ and $\varepsilon_{\Phi, i}$ are the pseudorange and phase-range observation noise and un-modelled multipath error.

In the real-time PPP model, the satellite orbit error, satellite clock error, code bias and phase delay can be corrected by applying internet broadcast data, the PCO and PCV can be eliminated by using IGS antenna files, and the earth tide effect and phase windup can be modelled precisely. Besides, there are two dominant errors: ionospheric delay and tropospheric delay, which should be considered carefully in the data processing. The former is often eliminated by applying the ionospheric-free combination of pseudorange and carrier phase measurements. The combination can be expressed as

(5) \begin{equation}\left[\begin{array}{c}{P_{I F}} \\ {\Phi_{I F}}\end{array}\right]=\left[f_{i}^{2} /\left(f_{i}^{2}-f_{j}^{2}\right)-f_{j}^{2} /\left(f_{i}^{2}-f_{j}^{2}\right)\right]\left[\begin{array}{l@{\quad}l}{P_{i}} & {\Phi_{i}} \\ {P_{j}} & {\Phi_{j}}\end{array}\right]\!,\end{equation}

where $ P_{IF} $ and $ \Phi_{IF} $ are the ionospheric-free combination of pseudorange and phase-range, $ f_i $ and $ f_j $ are the frequency of two difference GNSS signals in the same constellation, and $ P_i $ , $ \Phi_i $ , $ P_j $ and $ \Phi_j $ are the pseudorange and phase range of the corresponding frequencies.

The tropospheric delay is often estimated. It is divided into the zenith hydro-static delay and zenith delay, and the relative mapping functions:

(6) \begin{equation}T_{r}^{s}=m_{H} Z_{H, r}+m_{W}\left(Z_{T, r}-Z_{H, r}\right)\!,\end{equation}

where $Z_{H, r}$ is the tropospheric zenith hydro-static delay, which can be calculated by applying the Saastamoinen model. $Z_{T, r}$ is the tropospheric Zenith Total Delay (ZTD), which should be estimated in the filtering process. $ m_H $ is the hydro-static mapping function, and $m_{W}$ is the wet mapping function. For the mapping functions, we use the Niell Mapping Function (NMF) [17] in this paper. The ZTD is usually modelled as a random walk process:

(7) \begin{equation}\delta Z_{T, r, k}=\delta Z_{T, r, k-1}+\omega_{Z},\end{equation}

where $\omega_{z} \sim N\left(0, q_{T}^{2} \Delta t\right)$ , $ q_T $ represents the power spectral density (PSD) of ZTD and $ \Delta{t} $ is the time duration between the k-th and $ k-1 $ -th time step.

By applying the corrections and ionospheric-free combination, the observation equations can be rewritten as follows:

(8) \begin{equation}\tilde{P}_{r, I F, M}^{s}=\rho_{r}^{s}+c d t_{r, M}+T_{r}^{s}+\varepsilon_{P, I F},\end{equation}

(9) \begin{equation}\tilde{\Phi}_{r, I F, M}^{s}=\rho_{r}^{s}+c d t_{r, M}+T_{r}^{s}+\lambda_{I F, M} B_{r, I F, M}^{s}+\varepsilon_{\Phi, I F},\end{equation}

where the subscript IF indicates the ionosphere-free combination, and M indexes the satellite systems, including GPS, GLONASS, BDS and Galileo. The equations show different receiver clock errors between the satellite systems, which is called inter-system bias (ISB), being due to the different signal structure and different hardware delay for each GNSS system. The ISB should be considered and estimated when multi-GNSS observations are used together in the PPP data process^{(Reference Gao18)}. It can be written as

(10) \begin{equation}\left[\begin{array}{c}{\mathrm{ISB}_{G, R}} \\[3pt] {\mathrm{ISB}_{G, C}} \\[3pt] {\mathrm{ISB}_{G, E}}\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{-1} & {1} & {0} & {0} \\[3pt] {-1} & {0} & {1} & {0} \\[3pt] {-1} & {0} & {0} & {1}\end{array}\right]\left[\begin{array}{c}{d t_{r, G}} \\[3pt] {d t_{r, R}} \\[3pt] {d t_{r, C}} \\[3pt] {d t_{r, E}}\end{array}\right]\!,\end{equation}

where ISB $ _{M,N} $ is the ISB of system N with respect to system M, and the subscripts G, R, C, E correspond to GPS, GLONASS, BDS and Galileo, respectively. Similar to the tropospheric ZTD, ISB is also modelled as a random walk process as shown in Equation (7).

For each observation, the noise is not the same, which is related to the atmospheric thickness along the signal path, multipath effects, the User Range Accuracy (URA) and the carrier-to-noise ratio. We use the following coefficients to represent the observation noise:

(11) \begin{equation}\sigma^{2}=F_{M}^{s} R_{r}\left(a_{\sigma}^{2}+b_{\sigma}^{2} / \sin E l_{r}^{s}\right)\!,\end{equation}

where $ F_M^s $ is the URA of system M, $ R_r $ is the inverse of the standardised carrier-to-noise error ratio, and $a_{\sigma}$ and $b_{\sigma}$ are error factors. Empirically, we set $a_{\sigma}$ as 0.003m for phase-range measurements and 0.3m for pseudorange measurements, and set $b_{\sigma}$ as the same as $a_{\sigma}$ .

Finally, the state parameters for four-constellation real-time PPP is

(12) \begin{equation}{\textbf{\textit{x}}}_{P P P}=\left[{\textbf{\textit{p}}}_{r}, {\textbf{\textit{v}}}_{r}, {\textbf{\textit{d}}} {\textbf{\textit{t}}}_{r}, {{\textbf{\textit{d}}}\dot{\textbf{\textit{t}}}}_{r}, Z_{T, r}, {\textbf{\textit{B}}}_{r, I F, M}^{s}\right]^{T},\end{equation}

where $p_{r}=\left[x_{r}, y_{r}, z_{r}\right]^{T}$ is position of the antenna mass centre in ECEF coordinates, and $d t_{r}=\left[d t_{r, G}, \operatorname{ISB}_{G, R}, \operatorname{ISB}_{G, C}, \operatorname{ISB}_{G, E}\right]^{T}$ is the clock vector containing the GPS receiver clock delay and ISBs of the other three systems. ${\textbf{\textit{v}}}_{r}$ represents the velocity, and ${{\textbf{\textit{d}}}\dot{\textbf{\textit{t}}}}_{r}$ is the frequency drift. Because $ {\textbf{\textit{v}}}_{r} $ and $ {{\textbf{\textit{d}}}\dot{\textbf{\textit{t}}}}_{r} $ can be derived respectively from $ {\textbf{\textit{p}}}_{r} $ and $ {\textbf{\textit{d}}} {\textbf{\textit{t}}}_{r} $ , the total number of states includes three positions, m clocks (for an m-constellation system), one tropospheric ZTD and n cycle ambiguities (assuming there are n satellites in view). Therefore, the minimum number of visible satellites to support PPP is $ 4+m $ .

2.2 IMU model

The IMU consists of an accelerometer and gyroscope, which can provide three-axis specific force and three-axis angular velocity measurements. These measurements can be employed in an Inertial Navigation System (INS) for dead reckoning and attitude determination. When using the ECEF frame as the navigation frame, the INS dynamic can be written as^{(Reference Liu, Zhan, Liu and Liu20)}

(13) \begin{equation}{\dot{\textbf{\textit{p}}}}^{e}={\textbf{\textit{v}}}^{e},\end{equation}

(14) \begin{equation}{\dot{\textbf{\textit{v}}}}^{e}={\textbf{\textit{C}}}_{m}^{e} \,{\textbf{\textit{f}}}^{m}-\left[2 {\boldsymbol{\omega}}_{i e}^{e}\right]^{\times} {\textbf{\textit{v}}}^{e}+G^{e},\end{equation}

(15) \begin{equation}\dot{\textbf{\textit{C}}}_{m}^{e}=-\left[\omega_{i e}^{e}\right] {\textbf{\textit{C}}}_{m}^{e}+{\textbf{\textit{C}}}_{m}^{e}\left[{\boldsymbol{\omega}}_{im}^{m}\right]^{\times},\end{equation}

where e is the ECEF frame, m is the IMU measurement reference frame, i is the inertial frame, $ {\boldsymbol{\omega}}_{xy} $ represents the angular velocity of the y -frame with respect to the x -frame, $ {\textbf{\textit{C}}}_{m}^{e} $ is the rotation matrix from the m-frame to the e-frame, $[*]^{\times}$ represents matrix cross product, $ {\boldsymbol{\omega}}_{im}^{m} $ and $\ {\textbf{\textit{f}}}^{m} $ is the angular velocity and specific force of the IMU output, respectively, and $ G^{e} $ is the gravity vector in the e -frame.

The IMU device has a series of error sources in both the specific force and angular velocity measurements, especially for a MEMS IMU. The error models can be written as follows:

(16) \begin{equation}\tilde{{\boldsymbol{\omega}}}_{i m}^{m}=\left({\textbf{\textit{I}}}_{3}+{\textbf{\textit{M}}}_{g}\right) {\boldsymbol{\omega}}_{i m}^{m}+{\textbf{\textit{b}}}_{g}+{\textbf{\textit{w}}}_{g},\end{equation}

(17) \begin{equation}\tilde{{\textbf{\textit{f}}}}^{m}=\left({\textbf{\textit{I}}}_{3}+{\textbf{\textit{M}}}_{a}\right) {\textbf{\textit{f}}}^{m}+{\textbf{\textit{b}}}_{a}+{\textbf{\textit{w}}}_{a},\end{equation}

where $ \tilde{{\boldsymbol{\omega}}}_{i m}^{m} $ and $\ \tilde{{\textbf{\textit{f}}}}^{m} $ are the true angular rate and acceleration, respectively. $ {\textbf{\textit{b}}}_i(i=g,a) $ is the bias error, and $ {\textbf{\textit{I}}}_3 $ is the three-dimensional identity matrix. $ {\textbf{\textit{M}}}_i(i=g,a) $ is the mixed matrix of the scale factor error ( $ {\textbf{\textit{S}}} $ ) and cross-coupling error ( ${\textbf{\textit{K}}}$ ), as shown in Equation (18). $ {\textbf{\textit{w}}}_i(i=g,a) $ represents other un-modelled errors, which can be regarded as white noise.

(18) \begin{equation}{\textbf{\textit{M}}}_{i}={\textbf{\textit{S}}}_{i}+{\textbf{\textit{K}}}_{i}=\left[\begin{array}{c@{\quad}c@{\quad}c}{s_{x}^{i}} & {k_{z}^{i}} & {-k_{y}^{i}} \\[5pt] {-k_{z}^{i}} & {s_{y}^{i}} & {k_{x}^{i}} \\[5pt] {k_{y}^{i}} & {-k_{x}^{i}} & {s_{z}^{i}}\end{array}\right]\end{equation}

The dynamic behaviour of the variation of the IMU cross-coupling error can be modelled as white noise and can be absorbed into $ {\textbf{\textit{w}}}_i(i=g,a) $ ^{(Reference Groves21)}. The scale factor and bias can be described by the first-order Gauss–Markov procedure^{(Reference Robert and Hwang22)}, which can be expressed as

(19) \begin{equation}\left[\begin{array}{c}{{\textbf{\textit{b}}}_{k}} \\[5pt] {{\textbf{\textit{S}}}_{k}}\end{array}\right]=\left[\begin{array}{c}{e^{-\Delta t / \tau_{b}} {\textbf{\textit{b}}}_{k-1}} \\[5pt] {e^{-\Delta t / \tau_{s}} {\textbf{\textit{S}}}_{k-1}}\end{array}\right]+\left[\begin{array}{c}{{\textbf{\textit{w}}}_{b}} \\[5pt] {{\textbf{\textit{w}}}_{s}}\end{array}\right]\!,\end{equation}

where ${\textbf{\textit{w}}}_{i} \sim N\left(0,2 \sigma^{2} \Delta t /\tau_{i}\right)(i=b, s)$ , $ \Delta{t} $ is the time difference between the $ k$ -th time step and the $ k-1 $ -th time step, $\tau_{i}(i=b, S)$ is the correlation time, $ {\textbf{\textit{w}}} $ is white noise and $ \sigma $ is the variance of $ {\textbf{\textit{w}}} $ determined by the instability of the IMU bias and scale factor error. Commonly, one can find the variance of the errors from the product specifications. However, since there are many other unexpected factors such as the temperature, device ageing etc. that may affect the stability of the biases and scale factor errors, we use larger measurement error variances than the specifications in our implementation to accommodate such uncertainties.

Finally, the estimated parameters for the IMU model are

(20) \begin{equation}{\textbf{\textit{x}}}_{I M U}=\left[{\textbf{\textit{p}}}_{r}, {\textbf{\textit{v}}}_{r}, {\textbf{\textit{b}}}_{a}, {\textbf{\textit{b}}}_{g}, {\textbf{\textit{s}}}_{a}, {\textbf{\textit{s}}}_{g}\right]^{T},\end{equation}

where ${\textbf{\textit{s}}}_{i}=[s_{x}^{i}, s_{y}^{i}, s_{z}^{i}]^{T}(i=a, g)$ is the diagonal element of the matrix ${\textbf{\textit{M}}}_{i}$ .

In IMU-only navigation, the user has to apply compensation or calibration algorithms to make the IMU outputs accurate. Meanwhile, in integrated GNSS/INS algorithms, the errors of the IMU can be estimated in real time with the assistance of GNSS. In this case, only one needs to obtain proper initial states for the IMU biases and scale factors, as well as their variances. In our implementation, the initial scale factors are set to zero, and the biases are estimated using the first 1s of IMU data. Hence, the UAV must remain stable for the first 1s after powering-on. The initial biases are estimated as

(21) \begin{equation}\bm{b}_{g,0} = \bar{{\boldsymbol{\omega}}}_{i m}^{m},\end{equation}

(22) \begin{equation}\bm{b}_{a,0} = \bar{{\textbf{\textit{f}}}}^{m} + \bm{C}_e^m G^e,\end{equation}

where $ \bar{{\boldsymbol{\omega}}}_{i m}^{m} $ and $\ \bar{{\textbf{\textit{f}}}}^{m} $ are the 1-s-averaged IMU outputs.

2.3 Multi-GNSS real-time PPP/IMU tight integration

The tight integration of PPP and IMU utilises raw observations of GNSS (pseudorange and phase-range) and IMU (specific force and angular velocity) to achieve higher navigation quality than each method individually. The Extended Kalman Filter (EKF) is used to deal with the nonlinear terms of the states and measurements and combine the two kinds of signal. The estimated parameters for EKF processing are

(23) \begin{equation}\delta {\textbf{\textit{x}}}=\left[\delta {\textbf{\textit{p}}}_{r}, \delta {\textbf{\textit{v}}}_{r}, \delta {\textbf{\textit{b}}}_{a}, \delta {\textbf{\textit{s}}}_{a}, \delta {\textbf{\textit{b}}}_{g}, \delta {\textbf{\textit{s}}}_{g}, \delta {\textbf{\textit{d}}} t_{r}, \delta {\textbf{\textit{d}}} {\textbf{\textit{t}}}_{r}, \delta Z_{T, r}, \delta {\textbf{\textit{B}}}_{r, I F, M}^{s}\right]^{T},\end{equation}

where $ \delta $ represents the error (update) between two consecutive measurements. The definitions of the parameters can be found in Equations (12) and (21).

To apply the EKF algorithm, one should first derive the linearised state-space representation, which is given as

(24) \begin{align}\delta {\textbf{\textit{x}}}_k &= {\textbf{\textit{F}}}_k \delta {\textbf{\textit{x}}}_{k-1} + {\textbf{\textit{G}}}_k {\textbf{\textit{w}}}_{k-1} \nonumber\\[5pt] \delta {\textbf{\textit{z}}}_k &= {\textbf{\textit{H}}}_k \delta {\textbf{\textit{x}}}_{k} + \boldsymbol{\varepsilon}_{k},\end{align}

where $ {\textbf{\textit{F}}}_k $ is the state transition matrix and $ {\textbf{\textit{H}}}_k $ is the measurement matrix.

The process model can be derived from the vehicle motion model, GNSS model, INS dynamics model and INS error model as follows:

(25) \begin{equation}\delta {\textbf{\textit{x}}}_{k}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{\textbf{\textit{F}}}_{m}} & {} & {} & {} & {} \\ {} & {{\textbf{\textit{F}}}_{g}} & {} & {} & {} \\ {} & {} & {{\textbf{\textit{F}}}_{t}} & {} & {} \\ {} & {} & {} & {1} & {} \\ {} & {} & {} & {} & {{\textbf{\textit{I}}}_{n}}\end{array}\right]\delta {\textbf{\textit{x}}}_{k-1}+\left[\begin{array}{cc}{G_{m}} & {} \\[4pt] {} & {I_{n+5}}\end{array}\right] {\textbf{\textit{w}}}_{k-1},\end{equation}

where $\delta \bm{x}$ is the state vector, which is given in Equation (21), and $ \bm{w} $ is the noise vector, which is

(26) \begin{equation}{\textbf{\textit{w}}}=\left[{\textbf{\textit{w}}}_{v}, {\textbf{\textit{w}}}_{b_{a}}, {\textbf{\textit{w}}}_{s_{a}}, {\textbf{\textit{w}}}_{b_{g}}, {\textbf{\textit{w}}}_{s_{g}}, {\textbf{\textit{w}}}_{d t}, {\textbf{\textit{w}}}_{d f}, {w}_{Z}, {\textbf{\textit{w}}}_{B}\right]^{T}.\end{equation}

The items represent the white noises from the vehicle velocity, accelerometer bias and scale factor, the gyroscope bias and scale factor, the receiver clock offset (from the different navigation systems) and frequency drift, the tropospheric delay and the receiver phase ambiguity, respectively. We define the covariance matrix of $ {\textbf{\textit{w}}} $ as $ {\textbf{\textit{W}}}_k $ , which is a diagonal matrix determined by the user dynamics. The matrices $ \bm{F} $ and $ \bm{G} $ are given as follows:

(27) \begin{equation}{\textbf{\textit{F}}}_{m}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\textbf{\textit{F}}}_{p p}} & {{\textbf{\textit{F}}}_{p v}} & {} & {} \\[3pt] {{\textbf{\textit{F}}}_{v p}} & {{\textbf{\textit{F}}}_{vv}} & {{\textbf{\textit{C}}}_{m}^{e}} & {{\textbf{\textit{C}}}_{m}^{e} {\,\textbf{\textit{f}}}^{m}} \\[3pt] {} & {} & {\beta_{b_{a}}} \\[3pt] {} & {} & {} & {\beta_{s_{a}}}\end{array}\right]\!, \bm{G}_{m}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{}{1} & {} & {} \\[3pt] {} & {1} & {} \\[3pt] {} & {} & {1}\end{array}\right]\!,\end{equation}

where ${\textbf{\textit{F}}}_{p p}, {\textbf{\textit{F}}}_{w}, {\textbf{\textit{F}}}_{v p}, {\textbf{\textit{F}}}_{p v}$ are system dynamics matrices representing the relationship between the position and velocity. $\beta_{i}=e^{-\Delta t / \tau_{i}}$ is the state transition term of the first-order Gauss–Markov process.

(28) \begin{equation}\bm{F}_{g}=\left[\begin{array}{c@{\quad}c}{\beta_{b_{g}}} & {} \\ {} & {\beta_{s_{g}}}\end{array}\right]\!,\end{equation}

where $ \beta $ has the same definition as in Equation (25). Finally, $ \bm{F}_t $ represents the relationship between the clock offset and frequency drift of the GPS satellite clock and ISB items, having the following form:

(29) \begin{equation}\bm{F}_{d t_{r, G}}=\bm{F}_{{IS\ B}_{M, N}}=\left[\begin{array}{l@{\quad}l}{1} & {\Delta t} \\[3pt] {0} & {1}\end{array}\right]\!.\end{equation}

The linearised measurement model of the PPP/IMU tight integration filter is formed by the IMU-predicted pseudorange and phase-range and the GNSS measurements. The measurement equation is

(30) \begin{equation} \left[\begin{array}{c} {\tilde{\textbf{\textit{P}}}}_{i, IF}^{s}-{\textbf{\textit{P}}}_{IMU} \\[6pt] {\tilde{\boldsymbol{\Phi}}}_{i, I F}^{s}-{\boldsymbol{\Phi}}_{IMU}\end{array}\right]= \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} {-{\textbf{\textit{e}}}_{r, i}^{s}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\textbf{\textit{I}}}_{4} & {\boldsymbol{0}} & {\textbf{\textit{m}}}_{w} & {\boldsymbol{0}} \\[6pt] {-{\textbf{\textit{e}}}_{r, i}^{s}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\textbf{\textit{I}}}_{4} & {\boldsymbol{0}} & {\textbf{\textit{m}}}_{w} & {\boldsymbol{\lambda}}_{I F}\end{array}\right] \delta {\textbf{\textit{x}}}_{k-1}+{\boldsymbol{\varepsilon}}_{i, k},\end{equation}

where $ \tilde{\textbf{\textit{P}}}_{i, IF}^{s} $ and $ {\tilde{\boldsymbol{\Phi}}}_{i, I F}^{s} $ are the corrected GNSS pseudorange and phase-range measurements of GPS, GLONASS, BDS and Galileo. ${\textbf{\textit{e}}}_{r, i}^{s}=[\left({\textbf{\textit{p}}}_{x}^{s}-{\textbf{\textit{p}}}_{r, x}\right) / \rho_{r}^{s} \quad({\textbf{\textit{p}}}_{y}^{s}-{\textbf{\textit{p}}}_{r, y}) / \rho_{r}^{s} ({\textbf{\textit{p}}}_{z}^{s}- {\textbf{\textit{p}}}_{r, z}) / \rho_{r}^{s}]^{T}$ is the vector from a specific receiver to the satellite, $ {\textbf{\textit{m}}}_w $ is a column vector containing four same tropospheric delay wet mapping functions $ m_w $ , $ \boldsymbol{\lambda}_IF $ is the frequency corresponding to each satellite system and $ \boldsymbol{\varepsilon}_{i, k} $ is the GNSS measurement noise. The covariance matrix of $ \boldsymbol{\varepsilon}_{i, k} $ is $ {\textbf{\textit{R}}}_{i, k} $ , which is a diagonal matrix with the elements described by Equation (11). $ {\textbf{\textit{P}}}_{IMU} $ and $ {\boldsymbol{\Phi}}_{IMU} $ are the IMU-predicted GNSS measurements, which can be written as

(31) \begin{equation}\left[\begin{array}{c}{{\textbf{\textit{P}}}_{IMU}} \\[7pt] {\boldsymbol{\Phi}_{ {IMU }}}\end{array}\right]=\left[\begin{array}{c}{\left|{\textbf{\textit{p}}}_{ {r,IMU}}-{\textbf{\textit{p}}}^{\text {s }}\right|+c d t_{\text {r}, {\textbf{\textit{M}}}}+T_{r}^{s}} \\[4pt] {\left|{\textbf{\textit{p}}}_{{r,IMU}}-{\textbf{\textit{p}}}^{s}\right|+c d t_{ {r,M}}+T_{r}^{s}+\lambda_{IF, {\textbf{\textit{M}}}} B_{r,IF,M}^s}\end{array}\right]\!,\end{equation}

where $ {\textbf{\textit{p}}}_{r,IMU} $ is the IMU-predicted position calculated by Equation (23).

The EKF process can be described in two steps, the first of which is the prediction step:

(32) \begin{equation}\begin{array}{l}\hat{{\textbf{\textit{x}}}}_{k}^{-}={\textbf{\textit{F}}}_{k} \hat{{\textbf{\textit{x}}}}_{k-1} \\[5pt] {\textbf{\textit{Q}}}_{k}^{-}={\textbf{\textit{F}}}_{k} {\textbf{\textit{Q}}}_{k-1} {\textbf{\textit{F}}}_{k}^{T}+{\textbf{\textit{G}}}_{k} {\textbf{\textit{W}}}_{k}{\textbf{\textit{G}}}_{k}^T,\end{array}\end{equation}

where $ \hat{\cdot} $ indicates estimated values and $ {\cdot}^- $ indicates predicted values. $ {\textbf{\textit{Q}}}_k $ is the covariance matrix of states at the k -th time step. Commonly, we set the initial $ {\textbf{\textit{Q}}} $ matrix as an infinite diagonal matrix.

The second step of EKF is the update step, which is given as follows:

(33) \begin{align}{\textbf{\textit{K}}}_{k}&={\textbf{\textit{Q}}}_{k}^{-} {\textbf{\textit{H}}}_{k}^{T}\left({\textbf{\textit{H}}}_{k} {\textbf{\textit{Q}}}_{k}^{-} {\textbf{\textit{H}}}_{k}^{T}+{\textbf{\textit{R}}}_{k}\right)^{-1} \nonumber\\[4pt] \hat{{\textbf{\textit{x}}}}_{k}&=\hat{{\textbf{\textit{x}}}}_{k}^{-}+{\textbf{\textit{K}}}_{k}\left({\textbf{\textit{z}}}_{k}-{\textbf{\textit{H}}}_{k} \hat{{\textbf{\textit{x}}}}_{k}^{-}\right) \\[6pt] {\textbf{\textit{Q}}}_{k}&=\left({\textbf{\textit{I}}}-{\textbf{\textit{K}}}_{k} {\textbf{\textit{H}}}_{k}\right) {\textbf{\textit{Q}}}_{k}^{-},\nonumber\end{align}

where $ {\textbf{\textit{K}}}_{k} $ is the Kalman gain.

The tight integration scheme is shown in Fig. 1. The raw data of the GNSS receiver output is first combined with precise corrections and then sent to the integration process. In the GNSS receiver and IMU combined system, IMU measurements always arrive more frequently than from the GNSS receiver. Therefore, the filter must process IMU-only measurements when GNSS data are not available, by just applying the prediction step using Equation (30). When GNSS data arrive, EKF measurement updates will be applied. After updating the integrated state, IMU errors, including the biases and scale factors of the accelerometer and gyroscope, are fed back to the IMU error models and employed to correct the measurements.

Figure 1. Description of the real-time multi-GNSS PPP/IMU tight integration algorithm.

3.1 Open-sky performance

To compare the PPP performance between single- and multi-constellation GNSS, we applied different data processing with the same dataset. Note that the on-board real-time processing was set as the four-constellation GNSS PPP/IMU combination mode, whereas all other modes were processed off-board. The available satellites during this process and the corresponding Position Dilution of Precision (PDOP) are shown in Fig. 4.

The three GNSS modes and three GNSS/IMU integration modes of data processing were conducted using the same settings. The performance comparison in the NEU frame is shown in Fig. 5, and the corresponding RMS value is presented in Table 1. Since PPP needs tens of minutes to achieve convergence whereas only the performance during the flight period is the main focus here, we just take the time range from 46 to 53min to calculate the bias and RMS. The biases of the three direction errors are subtracted to make the plot clearer. Since the PPP bias is related to the quality of the precise orbit and clock, which lies beyond the scope of this article, we do not discuss it in detail. The reader can refer to Ref. [Reference Kazmierski, Sośnica and Hadas23] for more detailed information. The relative statistics are shown in Fig. 6.

According to the results, real-time PPP can provide an accurate position output, with a RMS of less than 0.025m in horizontal and 0.1m in vertical direction, after convergence. By applying two-GNSS PPP, the RMS is reduced by 58.2%, 30.0% and 22.1% for the north, east and up component, respectively, with respect to single-GNSS PPP. The RMS is further reduced by 70.3%, 73.8% and 72.8% when using four-system observation. Besides, PPP/IMU tight integration can provide a smoother position output and achieve a lower RMS. Compared with GNSS-only processing, the combined system reduces the RMS by 2.7%, 8.4% and 3.1% in the three directions respectively, for BDS-only processing. For the two.system case, it decreases by 7.2%, 2.6% and 8.3%, respectively, and for the four-system case, it is reduced by 3.9%, 3.1% and 15.3%, respectively.

To evaluate the performance of the vehicle velocity determination, the velocity error relative to the PPK velocity output is shown in Fig. 7. The corresponding RMS value is presented in Table 2, and the relative statistics are given in Fig. 8. The results show a great improvement when applying multi-GNSS PPP and PPP/IMU tight integration. The two-GNSS PPP reduces the RMS by 40.0%, 43.5% and 43.4% in the north, east and up directions, respectively, with respect to the single-GNSS PPP. The RMS is further reduced by 65.2%, 72.3% and 51.9% when using four-system observation. Besides, by applying PPP/IMU tight integration, the RMS reduces by 40.0%, 35.2% and 57.5% in the three directions respectively, for BDS-only processing. For the two-system case, it decreases by 27.6%, 16.5% and 55.9%, respectively, and for the four-system case, it is reduced by 9.1%, 3.8% and 28.4%, respectively.

Consequently, the performance of real-time PPP can be improved significantly by applying multi-constellation GNSS systems. This is because of the improvement in satellite availability and geometric distribution. PPP/IMU tight integration can smooth the position and velocity output and achieve a lower disturbance, especially for kinematic vehicle motion because of the high short-term precision of the IMU integral.

Table 1 Position RMS of the six modes of data processing (in mm)

Figure 4. Number of satellites and corresponding PDOP, where 1 Sys, 2 Sys and 4 Sys refer to BDS, GPS+BDS and GPS+BDS+GLONASS+Galileo, respectively.

Figure 5. Positioning performance comparison of the six modes of data processing. The biases of the three direction errors are subtracted to make the plot clearer. (a) PPP only (b) PPP/IMU integration.

Table 2 Velocity RMS of the six modes of data processing (in mm/s)

Figure 6. Positioning performance comparison of PPP and PPP/IMU tight integration.

Figure 7. Velocity determination performance comparison of the six modes of data processing. (a) PPP only (b) PPP/IMU integration.

3.2 GNSS outage simulation

GNSS outage is a common scenario in GNSS positioning, especially for carrier-phase-based processing. Internally, this occurs when the GNSS receiver becomes unlocked for some specific satellites, which leads to cycle slips or phase jumps. Externally, it occurs when a vehicle passes through a bridge, crowded building, canyon etc., which can cause data interruption. Unfortunately, PPP is very vulnerable to GNSS outage. It may divergence after such outage and takes several minutes to re-converge, which could be disastrous for users.

Figure 8. Velocity determination performance comparison of PPP and PPP/IMU tight integration.

Figure 9. Positioning performance comparison of PPP and PPP/IMU positioning for outage simulation (10s). (a) PPP only (b) PPP/IMU integration.

Table 3 Position and velocity RMS of the GNSS outage simulation (in mm and mm/s)

Figure 10. Velocity determination performance comparison of PPP and PPP/IMU positioning for outage simulation (10s). (a) PPP only (b) PPP/IMU integration.

To test the sensitivity of the system to GNSS outage, we simulated an outage by removing all satellites for 10s. The data outage began at 49min and ended at 49.17min. Since this test focused on a comparison between PPP only and the PPP/IMU tight integration system, only the four-system mode and its integration with IMU was utilised. The positioning result is shown in Fig. 9, and the velocity determination result is shown in Fig. 10. The corresponding RMS value is given in Table 3. Note that just the time range from 40 to 65min is plotted.

According to these results, the PPP-only system diverged when suffering an outage, while the PPP/IMU system could maintain a stable position output. The RMS position value for the PPP-only system was 110.7, 259.2 and 887.2mm for the north, east and up direction, respectively, whereas the PPP/IMU tight integration system reduced this to 24.5, 21.4 and 23.3mm for each direction component, corresponding to a 77.9%, 91.8% and 97.4% performance improvement in the three directions, respectively. The RMS velocity value was reduced from 24.3, 41.8 and 156.3mm/s to 5.0, 5.1 and 13.0mm/s in the three directions, corresponding to a 79.6%, 87.9% and 91.7% performance improvement in the three directions, respectively. This is because the position and velocity states of PPP can be estimated and maintained through the IMU short-term accurate position, the tropospheric delay, clock error and frequency drift can be maintained by the system dynamic modelling, and only ambiguities are re-estimated after the outage.

The integration of the IMU measurements leads to position and velocity drift when the period of GNSS outage extends. To test the performance of the PPP/IMU system under a long period of GNSS outage, we further simulated 10–60s of GNSS outage in steps of 10s. The position and velocity drift with respect to the outage period are shown in Fig. 11.

Figure 11. Position and velocity drift under different GNSS outage durations.

According to Fig. 11, the position and velocity errors increase as the outage period gets longer. The solutions drift little when the outage period is short because the PPP/IMU integration provides a high resolution of IMU biases when the GNSS measurements are sufficient. The drift becomes greater because the variances of the biases increase without correction from the GNSS information. Hence, users should bear this solution drift effect in mind in city or canyon environments where GNSS is not available for long periods. From the control perspective, the speed controller and trajectory controller cannot work without accurate velocity and position feedback. So when large drift occurs, the autopilot system of UAVs may not be reliable. However, the UAV can still land safely with a stand-alone IMU because the basic flight control system can work with only attitude feedback.

3.3 Insufficient GNSS measurement simulation

Besides outages, insufficient GNSS observation also occurs in city and canyon environments, causing the PPP-only system to fail to provide navigation solutions. It is common for satellites to be partially shadowed by tall buildings or cliffs. Hence, it is interesting to discuss the superiority of the proposed algorithm in such situations.

To test the performance of the proposed algorithm in an insufficient observation situation, we applied several azimuth masks to simulation shadowing of satellites, including $ 240^\circ $ , $ 300^\circ $ and $ 360^\circ $ azimuth masks. The $ 360^\circ $ azimuth mask is equivalent to a GNSS outage. Because the position and velocity showed little drift in the 10-s outage simulation in Section 3.2, we simulate a 30-s azimuth mask in this section to reveal the different performance in the simulations more clearly. The masks began at 49min and ended at 49.50min. The number of satellites is shown in Fig. 12. According to Section 2.1, this situation is insufficient for PPP positioning during the azimuth mask period in all of the groups. The results are shown in Figs 13 and 14. The corresponding RMS value is given in Table 4.

Figure 12. Number of observed satellites under different azimuth masks (30s).

Figure 13. Positioning performance comparison of PPP/IMU positioning under different azimuth masks (30s).

Figure 14. Velocity performance comparison of PPP/IMU positioning under different azimuth masks (30s).

Table 4 Position and velocity RMS of the insufficient measurement simulation with different degrees of azimuth masks (in mm and mm/s)

According to Fig. 12, the number of observed satellites is 4, 2 and 0 for the $ 240^\circ $ , $ 300^\circ $ and $ 360^\circ $ mask, respectively. All of these cases are insufficient for PPP-only positioning. According to Figs 13 and 14 and Table 4, the position and velocity solutions drift dramatically when the $ 360^\circ $ mask is applied. The filter re-initialises after the 30-s mask because the position error is too large. The system performs better when the $ 300^\circ $ mask is applied, which provides two observed satellites. Although only two satellites cannot be used in a common PPP filter, they help to shrink the drift rate of the position and velocity solutions in the proposed PPP/IMU system, reducing the RMS of position by 41.9%, 5.6% and 22.1% in the three directions, respectively, and the RMS of velocity by 56.7%, 63.5% and 44.1%, respectively. For the $ 240^\circ $ mask case, when the number of satellites rises to four, the drift rate of the solutions is further reduced, being 93.2%, 61.2% and 93.9% for position and 97.1%, 90.2% and 86.1% for velocity, respectively.

3.4 Computational load

Computational time is an important characteristic for embedded real-time systems. To study the computation time of the proposed PPP/IMU integration system, we tested the one-, two- and four-system PPP/IMU integration on the raspberry pi 3B platform. The results are shown in Fig. 15.

Figure 15. Computation time of the proposed PPP/IMU integration algorithm.

According to Fig. 15, the computational time is stable over time steps. The average computation time is 0.55, 0.80 and 1.72ms for one-system, two-system and four-system PPP/IMU integration, respectively. This indicates that the computational load is strongly related to the number of observed satellites. Since the number of satellites influences the precision and convergence speed of the system, there is a trade-off between the computational time and the precision and convergence period. Users are free to decide the priority.