**1. Introduction**

Since Precise Point Positioning (PPP) was first introduced [1,2], it has been a popular tool in numerous applications, for instance, GNSS (Global Navigation Satellite System) seismology [3,4], GNSS meteorology [5,6], deformation monitoring [7,8], etc. PPP utilizes only single-station GNSS data and orbit/clock products with high precision supplied via the International GNSS Service (IGS) centers to provide users with centimeter-level positioning solutions, and the errors contained in observations are remedied by appointed models, treated as unknown parameters, or parameterized after model correction. The research presented in [9–11] has proved the advantages of PPP, especially after applying the ambiguity fixing model [12,13] and the slant ionospheric delay and receiver Differential

**Citation:** Lv, J.; Gao, Z.; Xu, Q.; Lan, R.; Yang, C.; Peng, J. Assessment of Real-Time GPS/BDS-2/BDS-3 Single-Frequency PPP and INS Tight Integration Using Different RTS Products. *Remote Sens.* **2022**, *14*, 4367. https://doi.org/10.3390/rs14174367

Academic Editors: Chuang Shi, Shengfeng Gu, Yidong Lou and Xiaopeng Gong

Received: 16 July 2022 Accepted: 31 August 2022 Published: 2 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Code Bias (DCB) constraint models [14,15]. However, most of the works conducted thus far mainly consider conditions of dual-frequency observations. In contrast, Single-Frequency PPP (SF-PPP) is not as widely applied as dual-/triple-frequency PPP [16]. This is chiefly owing to the fact that no effective model can be utilized in SF-PPP to decrease or eliminate the effect of the ionospheric delay. Hence, this may directly lead to low positioning accuracy and long convergence time, especially in kinematic conditions.

Recently, GNSS satellite availability has significantly improved, along with the evolution of multi-constellation GNSS, and obtaining reliable and continuous positioning solutions using SF-PPP has become possible. The performance of three SF-PPP models, namely the GRoup And PHase Ionosphere Correction (GRAPHIC) model, GRAPHIC with code observation model, and an ionosphere-constrained model were analyzed in [17]. The contribution of the QZSS (Quasi-Zenith Satellite System) to GPS/BDS/GLONASS/Galileo SF-PPP was also presented. The outcomes proved that QZSS significantly increased the positioning accuracy of BDS and GLONASS when using simulated dynamic data. Such improvement was invisible while using the GPS-only, Galileo-only, and GPS/BDS/GLONASS/Galileo data. Furthermore, in the work of [17], the SF-PPP with ionosphere constraints presented a superior convergence performance, while the positioning accuracies supplied by the three SF-PPP models were close to each other. Moreover, with the completion of the BDS-3 constellation [18] on 23 June 2020, the modeling and performance provided by BDS-2/BDS-3 SF-PPP have become a new topic. Based on data from global distribution stations, Shi et al. [18] presented a BDS-2/BDS-3 SF-PPP model and studied the improvement of the combination based on BDS-2 and BDS-3 in precise positioning. The results indicated that the positioning accuracy was relevant to the receiver type adopted in different GNSS networks. The positioning performance based on the B1I signal is superior to that based on the B1C signal from the stations of the international GNSS Monitoring and Assessment System (iGMAS), while it is opposite at the stations of the Multi-GNSS EXperiment (MGEX). Contrasted with BDS-2 SF-PPP and BDS-3 SF-PPP, the convergence performance of BDS-2 + BDS-3 SF-PPP is significantly increased. As a major error in SF-PPP, ionospheric delay must be considered carefully. For ionospheric studies, Su et al. [19] proposed a novel single-frequency ionospheric-free-half PPP technique. To validate the reliability of the method, the performances of two conventional PPP models, and the proposed method, were evaluated. Results showed that slant ionospheric observables can be extracted with submeter-level accuracy towards the proposed novel method.

However, the continuous satellite availability will decrease in most dynamic applications when the receiver passes through overbridges, tunnels, and near towering constructions. Under these circumstances, the availability of GNSS satellites will be obstructed, and it is difficult to provide sufficient observations for PPP calculation. Therefore, high positioning accuracy with SF-PPP is difficult to achieve. Fortunately, according to [20], GPS positioning performance can be improved by fusing GPS with an Inertial Navigation System (INS), especially in challenging environments. In the last few years, profiting from the rapid development of PPP, PPP/INS integration has been regarded as one of the most efficient approaches to provide positions and attitudes. However, there are just a small number of works on SF-PPP/INS tight integration. Gao et al. [21] presented a multi-GNSS tight integration technique for single-frequency measurements and IMU outputs, by which the performance supplied by the ionospheric-delay- and receiver-DCB-constrained model was ameliorated significantly. The results indicated that the proposed model can obtain more accurate, continuous, and credible solutions in open-sky and GNSS-denied environments compared with traditional SF-PPP. The work in [22] presented a GNSS/INS tight integration with augmentations on ionospheric delay and tropospheric delay. Results illustrated that tight integration has a great effect in urban environments. Regional atmospheric augmentation increases the positioning accuracy of SF-PPP significantly on account of the sensibility of single-frequency signals to the ionospheric delay.

With the development of society, the requirement for real-time positioning is becoming urgent. To satisfy this demand, the Real-Time Services (RTS) of the IGS afford

precise orbit and clock products to global users through the internet. The work in [23] verified that the RTS products are, with high accuracy, comparing with the final products of the ESA (European Space Agency, Paris, France). Subsequently, the work in [24] appraised the positioning performance of BDS-2-only real-time PPP based on RTS products (CLK93) and MGEX stations. Results showed that the positioning accuracy acquired can be centimeter-/decimeter-level in static/kinematic environments. Kazmierski et al. [25] further investigated the quality of real-time products of CNES (Centre National d'Etudes Spatiales) for GPS, GLONASS, Galileo, and BDS-2. Results illustrated that the orbits and clocks of the GPS are the most accurate currently. Based on these backgrounds, this paper introduces a real-time GPS/BDS-2/BDS-3 SF-PPP/INS tight integration model. It is assessed via a set of vehicle-borne data and real-time orbit and clock products afforded by the IGS centers of CAS (Chinese Academy of Sciences), GFZ (Deutsche GeoForschungsZentrum), and WHU (Wuhan University). The paper is organized as follows: The methodology of the GPS/BDS-2/BDS-3 SF-PPP/INS tight integration model is described in Section 2. Then, the evaluations and the conclusions are presented in Sections 3 and 4.

#### **2. Methodology**

In this section, the methods of undifferenced and uncombined GPS/BDS-2/BDS-3 SF-PPP, recovery of real-time orbit and clock products, GPS/BDS-2/BDS-3 SF-PPP/INS tight integration, and parameter adjustments are described in detail.

## *2.1. Real-Time SF-PPP Model*

The linearized observation functions for the undifferenced and uncombined pseudorange and carrier phase are expressed as [21]

$$P\_{r,j}^{s} = \mathbf{u}\_r^s \cdot \mathbf{x} + \mathbf{c}t\_r + \mathbf{M}\mathbf{w}\_r^s \cdot \mathbf{Z}\mathbf{V}\mathbf{D}\_r + \gamma\_j \cdot I\_{r,1}^s + \mathbf{c}d\_{r,j} + e\_{r,j}^s \tag{1}$$

$$L\_{r\_{\eta}}^{s} = \mathbf{u}\_{r}^{s} \cdot \mathbf{x} + ct\_{r} + \mathbf{M} \mathbf{w}\_{r}^{s} \cdot \mathbf{Z} \mathbf{V} \mathbf{D}\_{r} - \gamma\_{\bar{\jmath}} \cdot I\_{r,1}^{s} - \lambda\_{\bar{\jmath}} N\_{r,\bar{\jmath}}^{s} + \boldsymbol{\varepsilon}\_{r,\bar{\jmath}}^{s} \tag{2}$$

where *P<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* and *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* refer to observed minus computed values of the original pseudorange and carrier phase observations; *s* represents satellite, *r* is the receiver, and *j* is signal frequency; **x** is the receiver position increments vector; **u***<sup>s</sup> <sup>r</sup>* is the direction cosine between receiver and satellite; *tr* is the receiver clock offset; ZWD*<sup>r</sup>* and Mw*<sup>s</sup> <sup>r</sup>* are the residual of tropospheric delay and the relevant mapping function; *I<sup>s</sup> <sup>r</sup>*,1 is the ionospheric delay at frequency 1 when the signal is passing through the propagation path; *γ<sup>j</sup>* is the multiplier factor dependent on signal frequency; *dr*,*<sup>j</sup>* represents receiver hardware time delay; *c* is the velocity of light; *λ<sup>j</sup>* is the wavelength and *N<sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* is integer ambiguity; *<sup>e</sup><sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* and *<sup>ε</sup><sup>s</sup> <sup>r</sup>*,*<sup>j</sup>* represent multipath biases, which are unmodeled, and the observation noise from the pseudorange and carrier phase, respectively.

The ionospheric delay is a key section that influences the performance of SF-PPP. Therefore, several models, for instance, the GIM (Global Ionosphere Map) data correction model [26], the single-frequency GRAPHIC model [27], and the single-frequency ZIDE (Zenith Ionospheric Delay Estimation) PPP model [28], have been proposed to limit such influence. For example, Øvstedal [26] recommended using the GIM to amend the ionospheric delay of the pseudorange and carrier phase. However, the drawback of the GIM-correction-based SF-PPP is that the residual of each satellite's ionospheric delay and the undisposed receiver DCB will degrade the performance of SF-PPP. To overcome this drawback, Montenbruck [29] tried to use the GRAPHIC model [27] instead of GIM data to eliminate the ionospheric delay in the carrier phase. In the GRAPHIC model, the characteristic that the ionospheric delays on pseudorange and carrier phase are equal in magnitude and opposite in sign for the same satellite is utilized. Hence, a linear ionospheric-free combination between pseudorange and carrier phase is adopted to form a new ionospheric-free carrier phase observation. Even so, the SF-PPP positioning accuracy still cannot satisfy the accuracy requirement because of the large noise of the new carrier phase [30]. Therefore, based on the GRAPHIC model, Beran et al. [28,31] proposed the ZIDE model to estimate

the zenith ionospheric delay of each satellite as a parameter, which has been verified to be effective to increase the performance of SF-PPP [30]. According to the works of [32–34], the position accuracy and convergence time of ZIDE SF-PPP can be ameliorated by adopting precise global and regional ionospheric models. However, the receiver DCB on pseudorange is ignored in the ZIDE model. Hence, the SF-PPP model based on Ionospheric delay and Receiver DCB Constraint (IRC) [21] is applied in this contribution. Here, to separate the receiver DCB from the ionospheric delay [14], the real-time ionospheric product provided by Wuhan University is utilized to generate a virtual external observation for each ionospheric delay. The basic observation equation can be written as

$$\begin{array}{l}P\_{r,f\_1}^{\mathcal{G}} = \mathbf{u}\_r^{\mathcal{G}} \cdot \mathbf{x} + \delta t\_r + \mathbf{M} \mathbf{w}\_r^{\mathcal{G}} \cdot \mathbf{Z} \mathbf{W} \mathbf{D}\_r + I\_{r,f\_1}^{\mathcal{G}} + d\_{r,f\_1}^{\mathcal{G}}\\P\_{r,f\_1}^{\mathcal{B}} = \mathbf{u}\_r^{\mathcal{B}} \cdot \mathbf{x} + \delta t\_r + \mathbf{M} \mathbf{w}\_r^{\mathcal{B}} \cdot \mathbf{Z} \mathbf{W} \mathbf{D}\_r + I\_{r,f\_1}^{\mathcal{B}} + d\_{r,f\_1}^{\mathcal{B}}\end{array} \tag{3}$$

$$\begin{array}{l} \mathbf{L}\_{r\_{\!\!\!\!=1}^{G}}^{\mathbf{G}} = \mathbf{u}\_{r}^{\mathbf{G}} \cdot \mathbf{x} + \delta t\_{\!\!r} + \mathbf{M} \mathbf{w}\_{r}^{\mathbf{G}} \cdot \mathbf{Z} \mathbf{W} \mathbf{D}\_{r} - I\_{r\_{\!\!\!=1}^{G}}^{\mathbf{G}} - \lambda\_{f\_{1}}^{\mathbf{G}} \mathbf{N}\_{r\_{\!\!\!=1}^{G}}^{\mathbf{G}} \\ \mathbf{L}\_{r\_{\!\!\!=1}^{B}}^{\mathbf{B}} = \mathbf{u}\_{r}^{\mathbf{B}} \cdot \mathbf{x} + \delta t\_{\!\!r} + \mathbf{M} \mathbf{w}\_{r}^{\mathbf{B}} \cdot \mathbf{Z} \mathbf{W} \mathbf{D}\_{r} - I\_{r\_{\!\!=1}^{G}}^{\mathbf{B}} - \lambda\_{f\_{1}}^{\mathbf{B}} \mathbf{N}\_{r\_{\!\!=1}^{B}}^{\mathbf{B}} \end{array} \tag{4}$$

$$I\_{r,f\_1}^s = 40.28 \cdot \text{VTEC} / \left(f\_1^2 \cos(Z\_\theta) \right) + \varepsilon\_{I\_{r,f\_1}^s}^s \varepsilon\_{I\_{r,f\_1}^s} \sim N\left(0, \sigma\_{\varepsilon\_{f\_1^s}}^2 \right) \tag{5}$$

$$
\begin{bmatrix} d\_{DCB} \\ 0 \end{bmatrix} = \begin{bmatrix} 1 & -1 \\ \frac{f\_1^2}{f\_1^2 - f\_2^2} & -\frac{f\_2^2}{f\_1^2 - f\_2^2} \end{bmatrix} \begin{bmatrix} d\_1 \\ d\_2 \end{bmatrix} \tag{6}
$$

where VTEC and *Z<sup>θ</sup>* stand for vertical total electron content acquired from the WHU's real-time ionospheric product and the zenith angle at the Ionospheric Punctuation Point (IPP); *εI<sup>s</sup> r*, *f* 1 is the accuracy of the real-time ionospheric product with prior variance of *σ*<sup>2</sup> *εIs r*, *f* 1 ; *dDCB* is the DCB between pseudoranges *P*<sup>1</sup> and *P*<sup>2</sup> in the same GNSS system; *f*<sup>1</sup> and *f*<sup>2</sup> are signal frequencies; *d*<sup>1</sup> and *d*<sup>2</sup> are hardware time delays in frequencies 1 and 2; the meanings of other parameters are the same as those in Equations (1) and (2).

In real-time data processing, products from IGS are transmitted into the RTCM-SSR (Radio Technical Commission for Maritime Services State Space Representation) form. The real-time corrections are the differences between precise products and broadcast products. The real-time products cannot be applied in positioning directly. Therefore, the products must be recovered to satisfy the requirement of PPP. SSR corrections for real-time orbit at the present epoch *t* can be calculated by corrections at the reference epoch *t*<sup>0</sup> using [35]

$$
\begin{bmatrix}
\Delta\_r \\
\Delta\_\mathbf{d} \\
\Delta\_\mathbf{c}
\end{bmatrix}\_t = \begin{bmatrix}
\Delta\_r \\
\Delta\_\mathbf{d} \\
\Delta\_\mathbf{c}
\end{bmatrix}\_{t\_0} + \begin{bmatrix}
\dot{\Delta}\_r \\
\dot{\Delta}\_\mathbf{d} \\
\dot{\Delta}\_\mathbf{c}
\end{bmatrix} (t - t\_0) \tag{7}
$$

where <sup>Δ</sup>*r*, <sup>Δ</sup>*a*, and <sup>Δ</sup>*<sup>c</sup>* are satellite position corrections in radial, along, and cross directions; . Δ*r*, . <sup>Δ</sup>*a*, and . Δ*<sup>c</sup>* are the corresponding velocities.

The rotation matrix *R* between the orbital coordinate system and Earth-Centered Earth-Fixed (ECEF) frame can be expressed as

$$\mathbf{R} = \begin{bmatrix} \frac{\dot{r}}{|\dot{r}|} \times \frac{r \times \dot{r}}{|r \times \dot{r}|} & \frac{\dot{r}}{|\dot{r}|} & \frac{r \times \dot{r}}{|r \times \dot{r}|} \end{bmatrix} \tag{8}$$

where *r* and . *r* stand for position and velocity vector for satellites calculated from the broadcast ephemeris, respectively.

The orbit corrections can be transferred from the orbital coordinate system to the ECEF frame via

$$
\begin{bmatrix}
\Delta\_{\mathbf{x}} \\
\Delta\_{\mathbf{y}} \\
\Delta\_{z}
\end{bmatrix}\_{t} = \mathbf{R} \cdot \begin{bmatrix}
\Delta\_{r} \\
\Delta\_{d} \\
\Delta\_{c}
\end{bmatrix}\_{t} \tag{9}
$$

where Δ*x*, Δ*y*, and Δ*<sup>z</sup>* are satellite orbit corrections in the ECEF frame.

Then, the precise products can be recovered by combining orbit corrections with satellite positions calculated by broadcast ephemeris

$$
\begin{bmatrix} X\_{prc} \\ Y\_{prc} \\ Z\_{prc} \end{bmatrix}\_t = \begin{bmatrix} X\_{brd} \\ Y\_{brd} \\ Z\_{brd} \end{bmatrix}\_t - \begin{bmatrix} \Delta\_x \\ \Delta\_y \\ \Delta\_z \end{bmatrix}\_t \tag{10}
$$

where *Xpre*, *Ypre*, and *Zpre* are satellite coordinates in ECEF frame at the time *t* after SSR corrections; *Xbrd*, *Ybrd*, and *Zbrd* are satellite coordinates calculated from broadcast ephemeris.

For real-time clock offset, the SSR corrections at the present epoch *t* can be calculated by corrections at the reference epoch *t*<sup>0</sup> by [35]

$$
\Delta t\_{\mathbb{C}} = \mathbb{C}\_0 + \mathbb{C}\_1(t - t\_0) + \mathbb{C}\_2(t - t\_0)^2 \tag{11}
$$

where Δ*tC* is the real-time clock correction in distance form calculated from SSR clock offsets. Then, the real-time precise clock can be obtained by

$$t\_{pre}^{sat} = t\_{brd}^{sat} - \frac{\Delta t\_C}{c} \tag{12}$$

where *t* sat pre is the satellite offset at epoch *t* after SSR corrections.

The signal–noise ratio, GNSS satellite elevation, as well as GNSS measurement environment are chiefly factors to influence the observation quality. Considering these terms, the conventional satellite elevation-dependent weight function is utilized here to figure out the prior variance [36]

$$\sigma^2 = \begin{cases} \sigma\_0^2 & E \ge 30^\circ \\ \sigma\_0^2 / (2 \cdot \sin(E))^2 & \text{else} \end{cases} \tag{13}$$

where *E* and *σ*<sup>2</sup> <sup>0</sup> are the elevation of satellites and the corresponding prior variance.
