*2.1. PPP Positioning Based on PPP-B2b Service*

As is known, precise satellite orbit and clock products play a key role in PPP processing. For PPP-B2b service, the satellite orbit corrections and clock offset corrections are broadcast by BDS-3 GEO satellites. Up to now, the PPP-B2b service has provided corrections only for BDS-3/GPS satellites. According to the PPP-B2b Interface Control Document (ICD) from China Satellite Navigation Office (CSNO) [26], the PPP-B2b orbit corrections are given in the radial (*δOradial*), along-track - *δOalong* , and cross-track (*δOcross*) directions. Hence, the orbit correction vector *<sup>δ</sup>OB*2*<sup>b</sup>* <sup>=</sup> *<sup>δ</sup>Oradial <sup>δ</sup>Oalong <sup>δ</sup>Ocross<sup>T</sup>* should first be transformed to the Earth-Center Earth-Fixed (ECEF) frame. This is because the satellite positions derived from CNAV1/LNAV broadcast ephemeris are based on the ECEF frame. This transformation can be described as follows:

$$
\delta \mathbf{X}^{sat} = \begin{bmatrix} \mathbf{e}\_{radial} & \mathbf{e}\_{along} & \mathbf{e}\_{cross} \end{bmatrix} \cdot \delta \mathbf{O}\_{B2b} \tag{1}
$$

with

$$\begin{array}{c} \mathbf{e}\_{radial} = \frac{r}{|r|}\\ \mathbf{e}\_{cross} = \frac{r \times \bar{r}}{|r \times \bar{r}|}\\ \mathbf{e}\_{along} = \mathbf{e}\_{cross} \times \mathbf{e}\_{radial} \end{array} \tag{2}$$

where *<sup>δ</sup>Xsat* <sup>=</sup> *δOx δOy δOz <sup>T</sup>* represents the PPP-B2b orbit correction vector in the ECEF frame. *r* and . *r* are the satellite position and velocity vectors derived from broadcast ephemeris. Then, the precise satellite position vector *<sup>X</sup>prec* can be computed by

⎧ ⎪⎪⎨ ⎪⎪⎩

$$\mathbf{X}\_{\text{precc}} = \mathbf{X}\_{\text{brdc}} - \delta \mathbf{X}^{\text{sat}} \tag{3}$$

where *<sup>X</sup>brdc* is the satellite position vector derived from broadcast ephemeris. As described in ICD document, the PPP-B2b clock correction parameter is defined as offset to the broadcast ephemeris clock in meters. The precise satellite clock offset is given by

$$dt\_{proc}^{sat} = dt\_{brdc}^{sat} - \frac{\mathbf{C}\_0}{c} \tag{4}$$

where *C*<sup>0</sup> represents the PPP-B2b clock correction parameter; *dtsat brdc* is the satellite clock offset derived from broadcast ephemeris; *dtsat prec* denotes the precise PPP-B2b clock offset; *c* is the velocity of light in a vacuum.

The GNSS raw code and carrier-phase measurements can read as [27,28]

$$\begin{cases} \begin{aligned} P\_i &= \rho + c \cdot (dt\_r - dt^s) + T + \frac{f\_1^2}{f\_i^2} \cdot I\_1 + B\_{r,i} - B\_i^s + \varepsilon\_{P\_i} \\ L\_i &= \rho + c \cdot (dt\_r - dt^s) + T - \frac{f\_1^2}{f\_i^2} \cdot I\_1 + \lambda\_i \left( N\_i + b\_{r,i} - b\_i^s \right) + \varepsilon\_{L\_i} \end{aligned} \end{cases} \tag{5}$$

where *i* represents the frequency number; *Pi* and *Li* are raw code and phase measurements; *ρ* denotes the geometric distance from satellite to receiver; *dtr* and *dt<sup>s</sup>* are the clock offsets at the receiver and satellite, respectively; *fi* is the frequency value, *T* is the tropospheric delay, and *I*<sup>1</sup> denotes the ionospheric delay for *L*1; *Br*,*<sup>i</sup>* and *B<sup>s</sup> <sup>i</sup>* represent the code bias at the receiver/satellite end; *br*,*<sup>i</sup>* and *b<sup>s</sup> <sup>i</sup>* denote the phase bias at the receiver/satellite end; *εPi* and *εLi* are the unmodelled errors of code/phase measurements. The relativistic effect, Sagnac effect, Shapiro time delay [29], site displacements [30], and phase windup [31] should be corrected according to the corresponding models.

The ionosphere-free code and phase combinations are usually adopted by PPP to eliminate the first-order ionospheric delay. The ionospheric-free code and phase combinations read as [32]:

$$\begin{cases} \begin{aligned} P\_{IF} &= \mathfrak{a} \cdot P\_{\tilde{t}} + (1 - \mathfrak{a}) \cdot P\_{\tilde{j}} = \mathfrak{\rho} + \mathfrak{c} \cdot d\tilde{t}\_{\tilde{r}} - \mathfrak{c} \cdot dt^{s} - B\_{IF}^{s} + T + \mathfrak{c}\_{P\_{IF}} \\ L\_{IF} &= \mathfrak{a} \cdot L\_{\tilde{i}} + (1 - \mathfrak{a}) \cdot L\_{\tilde{j}} = \mathfrak{\rho} + \mathfrak{c} \cdot d\tilde{t}\_{\tilde{r}} - \mathfrak{c} \cdot dt^{s} + T + \lambda\_{IF} \cdot \overline{\mathcal{N}}\_{IF} + \mathfrak{c}\_{L\_{IF}} \end{aligned} \end{cases} \tag{6}$$

where *α* = *f* <sup>2</sup> *i* / - *f* 2 *<sup>i</sup>* − *<sup>f</sup>* <sup>2</sup> *j* and *<sup>λ</sup>IF*·*NIF* = *<sup>λ</sup>IFNIF* − *<sup>b</sup><sup>s</sup> IF* + *br*,*IF* − *Br*,*IF*, *<sup>b</sup><sup>s</sup> IF* and *br*,*IF* are phase bias of ionospheric-free phase combination at the satellite/receiver end, *B<sup>s</sup> IF* and *Br*,*IF* represent code bias of ionospheric-free code combination at the satellite/receiver end, *λIF* and *NIF* are the ionospheric-free wavelength and the ionospheric-free ambiguity, respectively; *dtr* is the recombined receiver clock offset, which absorbs the receiver code bias of ionospheric-free code combination; *εPIF* and *εLIF* denote the unmodelled errors of

ionospheric-free code and phase combinations. It is noted that the ionospheric-free code bias at the satellite end should be corrected by applying PPP-B2b differential code bias corrections [15].

When the recovered PPP-B2b precise satellite orbits and clock offsets have been applied, the satellite orbit and clock errors are considered eliminated. The tropospheric delay can be divided into the dry and wet parts. The Saastamoinen model is usually used to correct the dry part, and the wet part must be estimated as unknown. Then the ionospheric-free code and phase combinations can be linearized as

$$\begin{cases} \begin{aligned} p\_{IF} &= -\mathbf{e} \cdot \mathbf{x} + \mathbf{c} \cdot d\overline{t}\_{I} + M\_{W} \cdot zwd + \varepsilon\_{P\_{IF}}\\ l\_{IF} &= -\mathbf{e} \cdot \mathbf{x} + \mathbf{c} \cdot d\overline{t}\_{I} + M\_{W} \cdot zwd + \lambda\_{IF} \cdot \overline{N}\_{IF} + \varepsilon\_{L\_{IF}} \end{aligned} \end{cases} \tag{7}$$

where *e*, *x* represent the unit vector from receiver to satellite, and the vector of position increment, respectively; *zwd* denotes the zenith wet delay and *MW* is the corresponding mapping function. In this equation, the remaining unknown parameters include only the receiver position increment vector *x*, the receiver clock offset *dtr*, the zenith wet delay *zwd* and the ionosphere-free ambiguity *λIFNIF*. By adopting a Kalman filter, the unknown parameters can be exactly estimated. It should be noted that in this paper the ionosphericfree Doppler combination is used to derive the velocity [23].

#### *2.2. INS Model*

The mechanization of INS in the navigation frame (*n*-frame) can be expressed as an integral process for the following equation [33,34]:

$$
\begin{bmatrix}
\dot{\boldsymbol{p}}\_{INS}^{\boldsymbol{n}} \\
\dot{\boldsymbol{\omega}}\_{INS}^{\boldsymbol{n}} \\
\dot{\boldsymbol{\mathcal{C}}}\_{b}^{\boldsymbol{n}}
\end{bmatrix} = \begin{bmatrix}
\boldsymbol{D}^{-1} \boldsymbol{v}\_{INS}^{\boldsymbol{n}} \\
\boldsymbol{\mathcal{C}}\_{b}^{\boldsymbol{n}} \boldsymbol{f}^{b} - \left(2\boldsymbol{\omega}\_{\dot{ic}}^{\boldsymbol{n}} + \boldsymbol{\omega}\_{c\boldsymbol{n}}^{\boldsymbol{n}}\right) \times \boldsymbol{v}\_{INS}^{\boldsymbol{n}} + \mathbf{g}^{\boldsymbol{n}} \\
\boldsymbol{\mathcal{C}}\_{b}^{\boldsymbol{n}} \left(\boldsymbol{\omega}\_{\dot{ib}}^{\boldsymbol{b}} \times\right) - \left(\boldsymbol{\omega}\_{\dot{ic}}^{\boldsymbol{n}} + \boldsymbol{\omega}\_{c\boldsymbol{n}}^{\boldsymbol{n}}\right) \times \boldsymbol{\mathcal{C}}\_{b}^{\boldsymbol{n}}
\end{bmatrix} \tag{8}
$$

where *D*−<sup>1</sup> is the diagonal matrix to transform the rectangular coordinates into geodetic coordinates; *<sup>p</sup><sup>n</sup> INS* and *<sup>v</sup><sup>n</sup> INS* denote the INS position vector and INS velocity vector in *<sup>n</sup>*-frame; *C<sup>n</sup> <sup>b</sup>* represents the transformation matrix from body-frame (*b*-frame) to *<sup>n</sup>*-frame; *<sup>f</sup> b* is the inertial measurement unit (IMU) specific force measurement projected in *b*-frame; *ωb ib* denotes the IMU angular rate measurement expressed in *<sup>b</sup>*-frame; *<sup>ω</sup><sup>n</sup> ie* is the vector of earth rotation rate expressed in *<sup>n</sup>*-frame; *g<sup>n</sup>* is the gravity vector presented in *<sup>n</sup>*-frame; *ω<sup>n</sup> en* represents the rotation rate of ECEF frame relative to that of *n*-frame expressed in *n*-frame; The symbol "×" represents the cross-product operator. Obviously, the INS update of position, velocity, and attitude can be implement based on Equation (8).

The perturbation of position, velocity and attitude must be considered in the INS data processing. In this paper, the error model of attitude, velocity, and position is defined as the psi-angle error model, which reads as

$$
\begin{bmatrix}
\delta \dot{p}\_{INS}^{\text{u}} \\
\delta \dot{\boldsymbol{w}}\_{INS}^{\text{u}} \\
\dot{\Phi}
\end{bmatrix} = \begin{bmatrix}
\end{bmatrix} \tag{9}
$$

where *<sup>δ</sup>p<sup>n</sup> INS* and *<sup>δ</sup>v<sup>n</sup> INS* are the error vector of INS position and velocity, respectively; *φ* indicates the vector of misalignment angles; *<sup>δ</sup>g<sup>n</sup>* is the gravity error; *<sup>δ</sup>f <sup>b</sup>* is error vector of specific force measurement; *δω<sup>b</sup> ib* denotes error vector of the IMU angular rate. Generally, the error vectors of specific force and angular rate measurement include scale factor errors, bias errors, and white noise [35]. Here, only bias errors are considered and the bias errors of IMU sensors are described as the first Gauss-Markov procedure [34]:

$$
\begin{bmatrix}
\delta \dot{\mathcal{B}}\_a \\
\delta \dot{\mathcal{B}}\_\mathcal{S}
\end{bmatrix} = -\frac{1}{\pi} \begin{bmatrix}
\delta \mathcal{B}\_a \\
\delta \mathcal{B}\_\mathcal{S}
\end{bmatrix} + \begin{bmatrix}
w\_a \\ w\_\mathcal{S}
\end{bmatrix} \tag{10}
$$

where *<sup>δ</sup>B<sup>g</sup>* and *<sup>δ</sup>B<sup>a</sup>* represent the error vector of gyro biases and accelerometer biases, respectively; *<sup>τ</sup>* is the correction time; and *<sup>w</sup>a*, *<sup>w</sup><sup>g</sup>* represents the driving noise for accelerometer and gyro, respectively.

#### *2.3. PPP-B2b/INS Loosely Coupled Integration Model*

A typical state vector of 15 states is used in PPP-B2b/INS loosely coupled integration, which can be expressed as

$$\mathbf{X}\_k = \begin{bmatrix} \delta p\_{INS}^n & \delta \boldsymbol{\sigma}\_{INS}^n & \boldsymbol{\Phi} & \delta \mathbf{B}\_a & \delta \mathbf{B}\_\circ \end{bmatrix}^T \tag{11}$$

The system model of PPP-B2b/INS loosely coupled integration in the discrete form can be simplified as

$$\mathbf{X}\_{k} = (\mathbf{I} + \mathbf{F} \cdot \mathbf{\hat{\boldsymbol{\omega}}} t) \mathbf{X}\_{k-1} + \mathbf{\hat{\Gamma}}\_{k-1} \mathbf{w}\_{k-1}, \mathbf{w}\_{k-1} \sim N(\mathbf{0}, \mathbf{Q}\_{k-1}) \tag{12}$$

where *<sup>X</sup>k*−<sup>1</sup> and *<sup>X</sup><sup>k</sup>* are the vector of state parameter at epoch *<sup>k</sup>* <sup>−</sup> 1 and *<sup>k</sup>*, respectively; *<sup>I</sup>* is the identify matrix; <sup>Δ</sup>*<sup>t</sup>* denotes the time interval between two adjacent epochs; *<sup>w</sup>k*−<sup>1</sup> represents the vector of the process noise at epoch *k* − 1, and **Γ***k*−<sup>1</sup> denotes the matrix of the noise distribution; *<sup>Q</sup>k*−<sup>1</sup> represents the covariance matrix of process noise; *<sup>F</sup>* is the dynamic matrix, which can be derived from the INS error model as follows:

$$F = \begin{bmatrix} F\_{rr} & F\_{rv} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ F\_{vv} & F\_{vv} & \left(\mathbf{f}^b \times \right) & \mathbf{C}\_b^n & \mathbf{0} \\ F\_{\phi r} & F\_{\phi v} & -\left(\omega\_{ic}^n + \omega\_{cn}^n \times \right) & \mathbf{0} & -\mathbf{C}\_b^n \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & -\frac{1}{\tau}I & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & -\frac{1}{\tau}I \end{bmatrix} \tag{13}$$

where **<sup>0</sup>** denotes the zero matrix; *<sup>F</sup>rr*, *<sup>F</sup>rv*, *<sup>F</sup>vr*, *<sup>F</sup>vv*, *<sup>F</sup>φr*, and *<sup>F</sup>φ<sup>v</sup>* are sub-matrices of dynamic matrix, which can be derived from Equation (9). Then the Kalman time update can be expressed as [36]

$$\begin{cases} X\_k^- = \Phi\_{k,k-1} X\_{k-1}^+ \\ P\_k^- = \Phi\_{k,k-1} P\_{k-1}^+ \Phi\_{k,k-1}^T + Q\_{k-1} \end{cases} \tag{14}$$

where the superscript "+" and "−" denote the estimated and predicted information; **<sup>Φ</sup>***k*,*k*−<sup>1</sup> <sup>=</sup> (*<sup>I</sup>* <sup>+</sup> *<sup>F</sup>*·Δ*t*) represents the state transition matrix; *<sup>P</sup><sup>k</sup>* represents the covariance matrix of the states.

When the PPP-B2b positioning solution and INS predicted solution are available at the same epoch, the Kalman filter measurement update can be operated. However, the positions and velocities from INS are referred to the IMU center, while those from PPP-B2b positioning are typically based on the GNSS antenna phase center. Therefore, the corresponding lever-arm offsets must be compensated. Here, the accurate lever-arm offsets are measured and directly applied to the position and velocity from INS. Then, the observation model for PPP-B2b/INS loosely coupled integration can be expressed as

$$\mathbf{Z}\_k = H\_k \mathbf{X}\_k + \begin{bmatrix} \mathfrak{e}\_p \\ \mathfrak{e}\_v \end{bmatrix}\_{\prime} \begin{bmatrix} \mathfrak{e}\_p \\ \mathfrak{e}\_v \end{bmatrix} \sim N(0, \mathbf{R}\_k) \tag{15}$$

with

$$\mathbf{Z}\_{k} = \begin{bmatrix} \mathbf{p}\_{GNSS}^{\mathrm{n}}\\ \mathbf{v}\_{GNSS}^{\mathrm{n}} \end{bmatrix} - \begin{bmatrix} \mathbf{p}\_{INS}^{\mathrm{n}}\\ \mathbf{v}\_{INS}^{\mathrm{n}} \end{bmatrix} \tag{16}$$

$$H\_k = \begin{bmatrix} I\_{3 \times 3} & \mathbf{0}\_{3 \times 3} & \left(\mathbf{C}\_b^n I^b \times \right) & \mathbf{0}\_{3 \times 3} & \mathbf{0}\_{3 \times 3} \\ \mathbf{0}\_{3 \times 3} & I\_{3 \times 3} & H\_\phi & \mathbf{0}\_{3 \times 3} & -\mathbf{C}\_b^n \left(I^b \times \right) \end{bmatrix} \tag{17}$$

where *<sup>H</sup><sup>φ</sup>* <sup>=</sup> *ω<sup>n</sup> ie* + *<sup>ω</sup><sup>n</sup> en* <sup>×</sup> *C<sup>n</sup> b* - *l b*× <sup>+</sup> *C<sup>n</sup> b* - *l <sup>b</sup>* <sup>×</sup> *<sup>ω</sup><sup>b</sup> ib* × ; *pn GNSS*, *v<sup>n</sup> GNSS* denote the vector of position and velocity obtained from PPP-B2b positioning; *ε<sup>p</sup> ε<sup>v</sup> <sup>T</sup>* denotes the vector of measurement noise of position and velocity and *<sup>R</sup><sup>k</sup>* represents its related covariance matrix; *l <sup>b</sup>* represents the vector of lever-arm offsets; *<sup>H</sup><sup>k</sup>* is the design matrix; *<sup>Z</sup><sup>k</sup>* is defined as the innovation vector, which consists of the difference of positions and velocities between PPP-B2b positioning solution and INS predicted solution. Finally, the Kalman measurement update can be given as

$$\begin{cases} \mathbf{K}\_k = \mathbf{P}\_k^- \mathbf{H}\_k^T \left( \mathbf{H}\_k \mathbf{P}\_k^- \mathbf{H}\_k^T + \mathbf{R}\_k \right)^{-1} \\ \mathbf{X}\_k^+ = \mathbf{X}\_k^- + \mathbf{K}\_k (\mathbf{Z}\_k - \mathbf{H}\_k \mathbf{X}\_k^-) \\ \mathbf{P}\_k^+ = (\mathbf{I} - \mathbf{K}\_k \mathbf{H}\_k) \mathbf{P}\_k^- \end{cases} \tag{18}$$

where *<sup>K</sup><sup>k</sup>* is the Kalman gain.

### **3. Description of Experiments**

In order to assess the performance of BDS-3 PPP-B2b/INS loosely coupled integration, two land vehicle experiments in open-sky scenery and urban canyon scenery (hereinafter referred to as Experiment A and Experiment B, respectively) were carried out in Qingdao, China. In these two experiments, a tactical-grade IMU (ISA100C) [37] and a MEMS-grade IMU (ADIS-16505) [38] were equipped to evaluate the impact of IMU grade on BDS-3 PPP-B2b/INS loosely coupled integration. Two GNSS antennas were mounted on the roof of the land vehicle, and connected to a NovAtel PwrPak7 receiver [39] and a FRII-PLUS PPP-B2b receiver (http://www.femtomes.com, accessed on 20 October 2021). The FRII-PLUS receiver was only used to collected PPP-B2b messages, while all GNSS observations used in the evaluation were collected by the PwrPak7 receiver. The sampling rate of GNSS observations was 1 Hz. Figure 1 shows the installation of experimental equipment. Relevant specific information for the IMU sensors is listed in Table 1.

**Figure 1.** Experimental equipment mounted on the vehicle.


**Table 1.** Parameters of the IMU sensors used in the experiments.

During these two experiments, a GNSS base station was established on the roof of Engineering-C building in the China University of Petroleum (East China). It should be noted that the distance between the base and rover stations was less than 5 km. Using GNSS observations collected by PwrPak7, and IMU data from ISA100C, together with GNSS measurements at the base station, the smoothed RTK/INS tightly coupled solution was obtained by commercial Inertial Explorer software (IE 8.90) from the NovAtel company, and this solution was used as reference.

In PPP-B2b positioning, the B1C and B2a signals of BDS-3, L1 and L2 signals of GPS were selected to form the ionospheric-free combination, and the cut-off elevation angle was set to 7 degrees. The remaining processing strategies are shown in detail in Table 2.

**Table 2.** Strategies of PPP-B2b positioning.


#### **4. Result and Discussion**

In this section, the performance of BDS-3 PPP-B2b/INS loosely coupled integration was investigated in both open-sky and complex urban environments.

#### *4.1. Experiment A*

Experiment A was conducted from 16:30:40 to 16:59:59 on 3 December 2021 in GPS time. The trajectory of the land vehicle is shown in Figure 2.

During Experiment A, the velocities of the land vehicle were within ±20 m/s in the east and north directions and within ±1 m/s in the vertical direction. The number of visible BDS-3/GPS satellites and the position dilution of precision (PDOP) values are shown in Figure 3. The number of visible satellites ranged from 9 to 15 with an average value of 12.5, and the PDOP varied from 1.28 to 2.70 with an average value of 1.59.

Figure 4 shows the positioning errors of PPP-B2b only, PPP-B2b/MEMS-IMU loosely coupled integration, and PPP-B2b/tactical-IMU loosely coupled integration, respectively. As with traditional PPP processing, a convergence period is indispensable in PPP-B2b positioning and PPP-B2b/INS loosely coupled integration. In this study, the convergence time was defined for horizontal/vertical positioning accuracy better than 30 cm/60 cm, with such a positioning accuracy for at least 60 continuous epochs. The statistics of the positioning results are summarized in Table 3. It should be noted that positioning results before full convergence were excluded when calculating positioning accuracy. The biases

in the east, north and up directions were 6.7/23.7/26.5 cm for PPP-B2b only, and the corresponding values for PPP-B2b/tactical-IMU and PPP-B2b/MEMS-IMU loosely coupled integration schemes were 6.6/23.3/26.0 cm and 6.2/23.5/25.8 cm, respectively. In terms of the root mean square (RMS) value of positioning errors, the improvement of PPP-B2b/INS loosely coupled integration was not significant compared to PPP-B2b only. This is because Experiment A was carried out within an open-sky environment and thus PPP-B2b on its own was already able to obtain high-precision positioning accuracy.

**Figure 2.** The trajectory of land vehicle in Experiment A.

**Figure 3.** Number of visible BDS-3/GPS satellites and PDOPs in Experiment A.

**Figure 4.** Positioning errors of PPP-B2b only (**top** panel), PPP-B2b/MEMS-IMU loosely coupled integration (**middle** panel), and PPP-B2b/tactical−IMU loosely coupled integration (**bottom** panel), for Experiment A.

**Table 3.** Statistics of the positioning errors (unit: cm).


It is recognized that PPP-B2b only is unable to obtain a high-precision solution during GNSS signal outages, while PPP-B2b/INS loosely coupled integration continues to work by INS mechanization. Four periods of GNSS signal outage were simulated to present the superiority of PPP-B2b/INS loosely coupled integration. The periods of GNSS outage were simulated from the 600th/900th/1200th/1500th epoch, and each period lasted 30 s. The corresponding positioning errors of PPP-B2b/MEMS-IMU and PPP-B2b/tactical-IMU loosely coupled integration are plotted in Figures 5 and 6, respectively. Clearly, the positioning accuracy become progressively worse with increased time of GNSS outages in both PPP-B2b/MEMS-IMU and PPP-B2b/tactical-IMU loosely coupled integration schemes. When the GNSS signal outage time increased to 30 s, the positioning errors in the east, north, and up directions decreased to 300.0 cm, 498.0 cm, and 41.0 cm, respectively, for PPP-B2b/MEMS-IMU loosely coupled integration. The results for the PPP-B2b/tactical-IMU loosely coupled integration scheme dropped to 18.6 cm, 21.8 cm, and 6.1 cm in the east, north, and up components, respectively. Table 4 lists the statistics results, showing that the positioning performance of PPP-B2b/tactical-IMU loosely coupled integration is significantly better than that of PPP-B2b/MEMS-IMU loosely coupled integration.

**Figure 5.** Positioning errors of the PPP-B2b/MEMS−IMU loosely coupled integration during the simulated GNSS signal outages.

**Figure 6.** Positioning errors of the PPP-B2b/tactical−IMU loosely coupled integration during the simulated GNSS signal outages.


**Table 4.** Statistics of positioning errors for PPP-B2b/INS loosely coupled integration during GNSS outages (unit: cm).
