**1. Introduction**

One method utilized by the global navigation satellite system (GNSS) to accomplish high-precision positioning is known as precise point positioning (PPP) [1], which has many advantages, including the need for only one receiver configuration as well as its flexible operation and wide-area application. To achieve centimeter-level positioning, PPP requires a convergence time of around 30 min, which makes it challenging to meet user demands in real time and severely restricts the marketing of PPP technology applications. Therefore, PPP ambiguity resolution (AR) technology has been suggested as a solution to the PPP

**Citation:** Du, S.; Shu, B.; Xie, W.; Huang, G.; Ge, Y.; Li, P. Evaluation of Real-time Precise Point Positioning with Ambiguity Resolution Based on Multi-GNSS OSB Products from CNES. *Remote Sens.* **2022**, *14*, 4970. https://doi.org/10.3390/ rs14194970

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

Received: 2 September 2022 Accepted: 4 October 2022 Published: 6 October 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/).

convergence problem. The single-difference (SD) between satellites approach was first suggested by Ge et al. in 2008. The main idea was to estimate uncalibrated phase delay (UPD) products using the fractional portions of the float wide-lane (WL) and narrow-lane (NL) ambiguities. According to their experimental findings, more than 80% of the SD ambiguities from 450 stations over 14 days were fixed under the data test, and the positioning accuracy was increased by 30% in comparison to the float solution [2]. Meanwhile, the integer recovery clocks method and the decoupled clock model were proposed by Laurichese et al. [3,4] and Collins et al. [5,6], respectively, from the perspective of satellite clocks. Using the products provided by these methods, PPP-AR can also be implemented at the user end. Through theoretical derivation and a large amount of data analysis, Geng et al. and Shi et al. verified the equivalence of the three ambiguity resolution processes and the positioning performance [7,8]. The above and related studies enabled the implementation of PPP-AR, effectively shortening the convergence time and enhancing the positioning accuracy [9,10]. As a result, different forms of post-ambiguity resolution products have been released to users by several analysis centers. Chen et al. conducted an analysis of the performance of PPP-AR based on various ambiguity resolution products and demonstrated that these products could significantly increase positioning accuracy and reduce convergence times; however, the positioning performance varied, which suggests that the usefulness of ambiguity resolution products plays a role in PPP-AR performance [11].

The abovementioned research results indicate that PPP-AR technology can, to a certain extent, address the shortcomings of traditional PPP; however, these analyses were based on the post-processing mode. Various academics have conducted extensive studies on real-time PPP, including studies on the accuracy of real-time satellite orbit and clock offset products, the positioning performance of real-time PPP, real-time PPP for time transfer or tropospheric delay retrieval, and more [12–19]. Their findings have shown that realtime satellite orbit and clock offset products can be accurate enough to satisfy user needs, and real-time PPP positioning performance can achieve centimeter-level accuracy after convergence [13,16–19]. However, for the kinematic mode, the four-system positioning still requires longer than 15 min to converge to 0.1, 0.1, and 0.2 m for the east (E), north (N), and up (U) components, respectively [19]. Moreover, El-Mowafy presented a method that can guarantee real-time PPP with a 3D accuracy of less than a decimeter while real-time satellite orbit and clock offset products are unavailable [20]. With the demand for real-time AR products from users, CNES currently broadcasts real-time observable-specific signal bias (OSB) products to the world with open access, which allows for the user implementation of PPP-AR. However, there are a number of issues with real-time PPP-AR technology, including the inferior precision of real-time products compared to post-products, missing or outage real-time OSB corrections, etc. Thus, it is important to evaluate the quality of real-time OSB products and their performance in PPP-AR positioning in order to promote the application of real-time PPP.

In this contribution, the data quality of real-time multi-GNSS OSB products from CNES and the performance of the real-time PPP-AR were evaluated. In Section 2 of this paper, a dual-frequency mathematical model for PPP-AR based on OSB products is introduced. Then, in Section 3, the quality of the OSB products are analyzed using three metrics, the positioning accuracy and convergence time of the real-time PPP-AR are evaluated using various system combinations and positioning modes, and a method that can effectively avoid the interruption of short-term OSB products is validated. Finally, the conclusions are summarized.

#### **2. Methodology**

The dual-frequency ionosphere-free PPP model is first introduced in this section. Then, a method for recovering the integer features of ambiguities using OSB products is described. Finally, an AR method for an ionosphere-free combination is given.

#### *2.1. Dual-frequency Ionosphere-Free PPP Model*

The effects of first-order ionospheric delay can be removed using the dual-frequency ionosphere-free (IF) observation. The following are expressions for the code and carrier phase IF observations (*P<sup>s</sup> <sup>r</sup>*,*IF* and *<sup>L</sup><sup>s</sup> <sup>r</sup>*,*IF*, respectively) at a specific epoch:

$$\begin{cases} P\_{r,IF}^s = \rho\_r^s + t\_r - t^s + T\_r + (a\_{12}b\_{r,1} + \beta\_{12}b\_{r,2}) - \left(a\_{12}b\_1^s + \beta\_{12}b\_2^s\right) + e\_{r,IF}^s\\ L\_{r,IF}^s = \rho\_r^s + t\_r - t^s + T\_r + \left(a\_{12}B\_{r,1} + \beta\_{12}B\_{r,2}\right) - \left(a\_{12}B\_1^s + \beta\_{12}B\_2^s\right) + \lambda\_{IF}N\_{r,IF}^s + \varepsilon\_{r,IF}^s \end{cases} \tag{1}$$

$$\begin{cases} \lambda\_{IF} \mathcal{N}\_{r,IF}^s = \alpha\_{12} \lambda\_1 \mathcal{N}\_{r,1}^s + \beta\_{12} \lambda\_2 \mathcal{N}\_{r,2}^s\\ \mathcal{a}\_{12} = \frac{f\_1^2}{f\_1^2 - f\_2^2}, \mathcal{B}\_{12} = -\frac{f\_2^2}{f\_1^2 - f\_2^2} \end{cases} \tag{2}$$

where *s* and *r* denote the satellite and receiver, respectively; *ρ<sup>s</sup> <sup>r</sup>* is the geometric distance between the satellite and the receiver (m); *tr* and *t <sup>s</sup>* are the receiver and satellite clock (m), respectively; *Tr* is the slant tropospheric delay (m); *br*,*<sup>i</sup>* and *b<sup>s</sup> <sup>i</sup>* are the code hardware delay of the receiver and the satellite (m), respectively; *Br*,*<sup>i</sup>* and *B<sup>s</sup> <sup>i</sup>* are the phase hardware delay of the receiver and the satellite (m), respectively; *fi* is the frequency; *λIF* is the IF wavelength; *Ns <sup>r</sup>*,*IF* is the IF ambiguity (cycle); *<sup>e</sup><sup>s</sup> <sup>r</sup>*,12 and *<sup>ε</sup><sup>s</sup> <sup>r</sup>*,12 represent the sum of the measurement noises and the multipath effects for the code and carrier phase IF observations (m), respectively.

The satellite clock (*t s* ) in the PPP model is corrected using the clock offset products. Currently, the clock offset products are estimated via the IF combination, which causes the code IF hardware delay to be absorbed by the generated clock offset [21–23]. Therefore, the relationship between the clock offset products (*t s IF*) and the satellite clock (*t s* ) can be expressed as follows:

$$t\_{IF}^{s} = t^{s} + \left(\alpha\_{12}b\_{1}^{s} + \beta\_{12}b\_{2}^{s}\right) \tag{3}$$

By combining Equations (1) and (3), the linearized IF combined function model can be obtained as follows:

$$\begin{cases} \begin{array}{c} p\_{r,IF}^{s} = \mu\_{r}^{s} \cdot \mathbf{x} + t\_{r,IF} - t\_{IF}^{s} + T\_{r} + e\_{r,IF}^{s} \\\ I\_{r,IF}^{s} = \mu\_{r}^{s} \cdot \mathbf{x} + t\_{r,IF} - t\_{IF}^{s} + T\_{r} + \lambda\_{IF} \hat{\mathcal{N}}\_{r,IF}^{s} + \varepsilon\_{r,IF}^{s} \end{array} \tag{4}$$

where *p<sup>s</sup> <sup>r</sup>*,*IF* and *l s <sup>r</sup>*,*IF* denote the observed minus computed (OMC) IF observations of the code and carrier phase, respectively; *μ<sup>s</sup> <sup>r</sup>* is the unit vector from the receiver to the satellite, and *x* is the coordinate of the estimated parameter. The receiver clock will absorb the code hardware delay of the receiver, which can be represented by the following expression:

$$t\_{r,IF} = t\_r + (a\_{12}b\_{r,1} + \beta\_{12}b\_{r,2})\tag{5}$$

Since the code and carrier phase observations in the PPP model share the same receiver clock parameter, where the parameter is based on the code observation, the code hardware delay will be added to the carrier phase observation [24,25]. Thus, the estimated ambiguity absorbs both the code and phase hardware delays, which can be written as follows:

$$
\hat{N}\_{r,IF}^{s} = N\_{r,IF}^{s} + \left[ a\_{12} (B\_{r,1} - B\_1^s + b\_1^s - b\_{r,1}) + \beta\_{12} (B\_{r,2} - B\_2^s + b\_2^s - b\_{r,2}) \right] / \lambda\_{IF} \tag{6}
$$

where *N*ˆ *<sup>s</sup> <sup>r</sup>*,*IF* and *<sup>N</sup><sup>s</sup> <sup>r</sup>*,*IF* represent the float and integer ambiguities, respectively. The integer feature of the IF ambiguity (*N<sup>s</sup> <sup>r</sup>*,*IF*) is destroyed, as can be seen from the equation above. Therefore, the integer feature for the IF ambiguity should first be recovered before PPP-AR.

#### *2.2. Method for Recovering the Integer Feature of Ambiguity*

The same function model between the server and the user is required for traditional AR products, which significantly limits the application of the AR products. With the development of AR technology, Laurichesse et al. proposed the use of undifferenced and uncombined OSB products, which comprise both code bias and phase bias [26]. Based on these OSB products, a variety of function models can be used by users to conduct single-frequency or multi-frequency PPP-AR, which helps to meet the current demand for multi-frequency, multi-system, and multi-model PPP-AR. By directly utilizing the OSB products in the code and phase observations, the satellite's code and phase hardware delays can be corrected. The following is an expression for the dual-frequency code bias (*b s* <sup>1</sup> and *b s* 2) and phase bias (*B<sup>s</sup>* <sup>1</sup> and *<sup>B</sup><sup>s</sup>* 2) products [27]:

$$
\begin{pmatrix} \overline{b\_1^s} \\ \overline{b\_2^s} \end{pmatrix} = \begin{pmatrix} \beta\_{12} \begin{pmatrix} b\_1^s - b\_2^s \\ -a\_{12} \left( b\_1^s - b\_2^s \right) \end{pmatrix}, \begin{pmatrix} \overline{B\_1^s} \\ \overline{B\_2^s} \end{pmatrix} = \begin{pmatrix} a\_{12}b\_1^s + \beta\_{12}b\_2^s - B\_1^s \\ a\_{12}b\_1^s + \beta\_{12}b\_2^s - B\_2^s \end{pmatrix} \tag{7}
$$

Since there is one code and phase bias for every code and phase observation, respectively, it is only necessary to add the OSB products directly to the raw observations when they are being employed. Additionally, the IF code and phase biases are formed, which can be expressed as follows:

$$\begin{cases} \overline{\boldsymbol{b}}\_{IF}^{s} = \boldsymbol{a}\_{12} \overline{\boldsymbol{b}}\_{1}^{s} + \beta\_{12} \overline{\boldsymbol{b}}\_{2}^{s} = 0\\ \overline{\boldsymbol{B}}\_{IF}^{s} = \boldsymbol{a}\_{12} \overline{\boldsymbol{B}}\_{1}^{s} + \beta\_{12} \overline{\boldsymbol{B}}\_{2}^{s} = \boldsymbol{a}\_{12} \boldsymbol{b}\_{1}^{s} + \beta\_{12} \boldsymbol{b}\_{2}^{s} - \boldsymbol{a}\_{12} \boldsymbol{B}\_{1}^{s} - \beta\_{12} \boldsymbol{B}\_{2}^{s} \end{cases} \tag{8}$$

When *b s IF* and *<sup>B</sup><sup>s</sup> IF* are combined with Equation (1), the new IF ambiguity can be stated as follows:

$$
\hat{N}\_{r,IF}^{s} = N\_{r,IF}^{s} + \left[a\_{12}(B\_{r,1} - b\_{r,1}) + \beta\_{12}(B\_{r,2} - b\_{r,2})\right] / \lambda\_{IF} \tag{9}
$$

Because of the negative impact of the code and phase hardware delays on the receiver, the IF ambiguity still does not have the integer feature. Usually, the SD between satellites can be used to eliminate the negative impact of the receiver, and the IF-SD ambiguity will recover the integer feature.

## *2.3. PPP-AR Process*

The integer wide-lane (WL) and float narrow-lane (NL) ambiguities can be used to decompose the IF ambiguity in the IF combination, which translates as follows:

$$
\hat{N}\_{r,IF}^{s} = \left(\frac{f\_2}{f\_1 + f\_2} \lambda\_{WL} N\_{r,WL}^{s} + \lambda\_{NL} \hat{N}\_{r,NL}^{s}\right) / \lambda\_{IF} \tag{10}
$$

where *N<sup>s</sup> <sup>r</sup>*,*WL* is the integer WL ambiguity; *<sup>N</sup>*<sup>ˆ</sup> *<sup>s</sup> <sup>r</sup>*,*NL* is the float NL ambiguity; and *λWL* and *λNL* are the wavelengths of the WL and NL ambiguities, respectively, which can be formulated as follows: 

$$\begin{cases} \lambda\_{WL} = \frac{\mathbb{C}}{f\_1 - f\_2} \\ \lambda\_{NL} = \frac{\mathbb{C}}{f\_1 + f\_2} \end{cases} \tag{11}$$

where *c* denotes the speed of light.

The float NL ambiguity can be derived from Equation (10) when the float WL ambiguity is correctly fixed; the float WL ambiguity can be calculated through the Melbourne– Wübbena (MW) combination and fixed using the rounding method [2]. Then, the float WL and NL ambiguities can be summarized as follows:

$$\begin{cases} \begin{array}{c} \text{\(\hat{N}\_{r,WL}^{s} = \left(\frac{L\_{r,1}^{s}}{\lambda\_{1}} - \frac{L\_{r,2}^{s}}{\lambda\_{2}} - \frac{\lambda\_{2}P\_{r,1}^{s} + \lambda\_{1}P\_{r,2}^{s}}{\lambda\_{\text{WL}}\left(\lambda\_{2} + \lambda\_{1}\right)}\right) \\\ \text{\(\hat{N}\_{r,NL}^{s} = \frac{\lambda\_{1}P\_{r}^{s} \| \text{F}}{\lambda\_{NL}} - \frac{\lambda\_{1}N\_{r,1}^{s} \| \text{F} \|}{\lambda\_{2} - \lambda\_{1}} \end{array}} \tag{12}$$

The WL ambiguity is affected by measurement noise and observation errors to a lesser extent due to its long wavelength, and it can be calculated with high accuracy using multiepoch smoothing [28]. As a result, the rounding method can be used to fix it directly [2]. On the contrary, the float NL ambiguities from different satellites in the PPP model are highly correlated; hence, the LAMBDA method should be used for fixing. After inserting the fixed WL and NL ambiguities into (10), the fixed IF ambiguities are obtained. Then, a virtual observation is employed to constrain the filtering state, and the fixed solution is obtained.

It is important to note that when using the MW combination to obtain the float WL ambiguity, the antenna phase center correction of the receiver and satellite must be taken into account [29,30]. The formula is as follows:

$$\mathcal{N}\_{r,WL}^{s} = \left(\frac{L\_{r,1}^{s} + z\_{r,1}^{s}}{\lambda\_{1}} - \frac{L\_{r,2}^{s} + z\_{r,2}^{s}}{\lambda\_{2}} - \frac{\lambda\_{2} \cdot (P\_{r,1}^{s} + z\_{r,1}^{s}) + \lambda\_{1} \cdot (P\_{r,2}^{s} + z\_{r,2}^{s})}{(\lambda\_{2} + \lambda\_{1}) \cdot \lambda\_{WL}}\right) \tag{13}$$

$$\begin{cases} z\_{r,1}^s = z\_{r,1} \cdot \sin \theta\_r^s - z\_1^s \\ z\_{r,2}^s = z\_{r,2} \cdot \sin \theta\_r^s - z\_2^s \end{cases} \tag{14}$$

where *zr*,1 and *zr*,2 are the vertical phase center offsets of the receiver antenna for the two frequencies, and similarly, *z<sup>s</sup>* <sup>1</sup> and *<sup>z</sup><sup>s</sup>* <sup>2</sup> are those of the satellite antenna. *<sup>θ</sup><sup>s</sup> <sup>r</sup>* denotes the elevation angle of the satellite *s* with respect to the receiver *r*.

#### **3. Real-Time PPP-AR Performance**

In this section, we will first introduce the experimental data and processing strategies. Then, the quality of the real-time OSB products is analyzed. In addition, the performance of various PPP-AR combinations and modes based on the OSB products is assessed. Finally, a prediction method is proposed to effectively avoid the influence of short-term missing OSB products.

#### *3.1. Data and Strategy*

The post-store OSB products from day of year (DOY) 121 to 151 in 2021 were downloaded from http://www.ppp-wizard.net/products/REALTIME/ (accessed on 2 August 2022) to make it easier to analyze the real-time OSB products and PPP-AR.

The details of the dual-frequency and multi-GNSS real-time OSB products from CNES as of DOY 151 in 2021 are displayed in Table 1. It should be noted that the code bias was stable during the day [31]; therefore, the code bias was not analyzed.


**Table 1.** OSB product information provided by CNES.

**Note:** The sampling interval of post-store OSB products was 30 s. The "Frequency Number" only represents the index of the frequency for each system, which is convenient for the description below.

Figure 1 depicts the distribution of the 90 Multi-GNSS Experiment (MGEX) stations used to perform the PPP float solution and PPP-AR around the globe. These stations can receive GPS, Galileo, and BDS dual-frequency signals. The experiment time was from DOY 121 to 151 in 2021, which is a total of 31 days. To fully evaluate the real-time PPP-AR performance with different satellite systems, several combinations of static and kinematic positioning experiments were carried out. The combinations were as follows: 1) GPS-only; 2) GPS+Galileo; 3) GPS+Galileo+BDS. In the experimental analysis, incomplete observational data or data that did not pass quality checks (including detection, identification, and adaption (DIA) [32,33] and the Inter Quartile Range (IQR) method [34]) were excluded, where the excluded data were about 4% of the total data.

**Figure 1.** The 90 globally distributed MEGX stations used for the PPP-AR experiment.

All the PPP and PPP-AR experiments were performed with in-house software based on the secondary development of GAMP [35]. The phase windup was applied by the phase polarization effects. For the GPS L1/L2, Galileo E1/E5a, and BDS B1I/B3I frequencies, the igs14.atx file was utilized to correct the satellite phase center offset (PCO) and phase center variation (PCV); the receiver PCO and PCV of GPS were used for Galileo and BDS because those of Galileo and BDS were unavailable. It should be noted that the satellite orbit, clock offset, and OSB products from CNES were used in PPP-AR. The elevation and observable arc were used to define which ambiguity subsets should be employed, and the ratio test and bootstrapping success rate were used to validate whether the fixed subsets could be trusted. The detailed strategy for positioning is summarized in Table 2.


#### *3.2. Quality Analysis of Phase Bias*

The two key metrics for analyzing the quality of real-time OSB products are data availability (DA) and stability, where DA can be defined as the ratio of available epochs to all epochs in a day. Figure 2 displays the DA results for the GPS, Galileo, and BDS satellites based on the mean DA of each satellite at two frequencies during a period of 31 days. As can be observed, all GPS and Galileo satellites had a DA of more than 90%, whereas the majority of BDS satellites had DAs of less than 60%. The major reason for the inferior DAs of BDS satellites can be inferred as the subpar quality of the real-time satellite orbit and the offset products for BDS [37,38], which hinder the generation of the phase bias products in real time.

**Figure 2.** Mean DA of GPS (**left**), Galileo (**middle**), and BDS (**right**) satellites.

The mean DA of each satellite for each system is displayed in Table 3, which numerically reflects the DA of each system. Since frequency 1 and frequency 2 are two frequencies that are utilized to estimate the clock offset products of each system [39], the DAs of frequency 1 and frequency 2 are almostsame for each system. The overall DAs of the GPS and Galileo satellites were greater than 97%, while that of the BDS satellites was smaller than 50%. In conclusion, although real-time PPP-AR can be implemented using the OSB products, a consistent and reliable real-time PPP-AR service cannot be guaranteed due to its unstable DA, particularly in the case of BDS satellites. Therefore, a critical issue for real-time PPP-AR is how to avoid or weaken the influence of the absent phase bias products.

**Table 3.** Overall DA of GPS, Galileo, and BDS satellites.


The maximum fluctuation value (MAX) and standard deviation (STD), where MAX is the absolute value of the difference between the maximum and minimum phase biases within one day, can be applied to evaluate the stability of real-time OSB products. Given that the results for the two frequencies were nearly identical, the 31-day data were counted to obtain the mean MAX and STD of frequency 1 for each satellite; the results are shown in Figure 3. With the exception of G01 and G21, all other GPS satellites had MAX and STD values that were better than 0.08 and 0.025 cycles, respectively. The MAX and STD of Galileo were only marginally inferior to those of GPS, with values at around 0.12 and 0.05 cycles, respectively. Apart from C07, C13, and C24, the MAX and STD of the remaining BSD MEO and IGSO satellites were approximately 0.15 and 0.05 cycles, respectively, whereas the MAX and STD of the GEO satellites were noticeably greater than those of the other satellites, exceeding 0.5 and 0.2 cycles, respectively. This was due to the fact that the quality of the satellite orbit and clock offset products for GEO satellites is poorer than that for MEO and IGSO satellites, leading to more satellite orbit and clock offset errors being absorbed by the phase bias. It is worth mentioning that the MAX and STD trends for each satellite were consistent, indicating that fewer outliers were included in the phase bias, which is beneficial to achieving a more robust PPP-AR performance.

**Figure 3.** Mean MAX and STD of GPS (**top**), Galileo (**middle**), and BDS (**bottom**) satellites for frequency 1.

Similarly, the mean MAX and STD values of each system were calculated in order to depict the overall stability of each system numerically; the results are displayed in Table 4. The MAX and STD of the GPS, Galileo, and BDS systems were 0.045 and 0.012; 0.081 and 0.028; and 0.292 and 0.085 cycles, respectively, among which the BDS system was obviously inferior to the GPS and Galileo systems.


**Table 4.** Overall MAX and STD of GPS, Galileo, and BDS satellites.
