Evaluation of Real-time Precise Point Positioning with Ambiguity Resolution Based on Multi-GNSS OSB Products from CNES

Abstract

Ambiguity resolution (AR) is a core technology that helps to speed up convergence time and increase positioning accuracy for precise point positioning (PPP), and the performance of PPP-AR is based on the quality of ambiguity resolution products. Real-time PPP-AR becomes a reality as users can now obtain publicly accessible real-time observable-specific signal bias (OSB) products from the Centre National d’Etudes Spatiales (CNES). Therefore, an analysis of the quality of OSB products and an evaluation of the performance of PPP-AR are required to promote the application of real-time positioning. For a total of 31 days between day of year (DOY) 121 and 151 in 2021, observation data were collected from 90 stations, and the OSB products were used to assess the experiments. As for the quality of the OSB products, the data availability (DA) of the GPS and Galileo satellites was greater than 97%, whereas that of BDS was less than 60%; the maximum fluctuation value (MAX) and standard deviation (STD) of the GPS, Galileo, and BDS satellites were 0.045 and 0.012; 0.081 and 0.028; and 0.292 and 0.085 cycles, respectively. In terms of ambiguity residuals using the OSB products, the wide-lane (WL) residual percentages within ±0.25 cycles for the GPS, Galileo, BDS-2, and BDS-3 systems were more than 92%, and the narrow-lane (NL) residual percentages within ±0.25 cycles for the four systems were 92%, 89%, 79%, and 60%, respectively. With regard to the performance of PPP-AR, the GPS+Galileo solution showed the best performance in the kinematic positioning mode, in which the mean root mean square (RMS) of positioning accuracy was 1.06, 1.27, and 2.85 cm for the east (E), north (N), and up (U) components, respectively, and the mean convergence time reached 9.6 min. In the static positioning mode, the mean convergence times of the GPS-only and GPS+Galileo solutions were 11.4 min and 8.0 min, respectively, and both of their mean RMS for positioning accuracy reached 0.79, 0.95, and 1.48 cm for the E, N, and U components, respectively. However, the addition of BDS did not further enhance the performance of multi-GNSS PPP-AR in either the kinematic or static positioning mode due to the poor quality of the real-time BDS products. More importantly, a prediction method was proposed to avoid re-convergence and to enhance the reliability of PPP-AR in the event of short-time missing real-time OSB products and to improve the positioning accuracy and the ambiguity fixed rate.

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 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,13,14,15,16,17,18,19]. Their findings have shown that real-time 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,17,18,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_{r,IF}^s$ and $L_{r,IF}^s$, respectively) at a specific epoch:

$$\{\begin{array}{l}{P}_{r,IF}^s={\rho }_r^s+t_r-t^s+T_r+\left({\alpha }_{12}{b}_{r,1}+{\beta }_{12}{b}_{r,2}\right)-\left({\alpha }_{12}{b}_1^s+{\beta }_{12}{b}_2^s\right)+e_{r,IF}^s\hfill \\ L_{r,IF}^s={\rho }_r^s+t_r-t^s+T_r+\left({\alpha }_{12}{B}_{r,1}+{\beta }_{12}{B}_{r,2}\right)-\left({\alpha }_{12}{B}_1^s+{\beta }_{12}{B}_2^s\right)+{\lambda }_{IF}{N}_{r,IF}^s+{\epsilon }_{r,IF}^s\hfill \end{array}$$

(1)

$$\{\begin{array}{l}{\lambda }_{IF}{N}_{r,IF}^s={\alpha }_{12}{\lambda }_1{N}_{r,1}^s+{\beta }_{12}{\lambda }_2{N}_{r,2}^s\hfill \\ {\alpha }_{12}=\frac{f_1^2}{f_1^2-f_2^2},{\beta }_{12}=-\frac{f_2^2}{f_1^2-f_2^2}\hfill \end{array}$$

(2)

where $s$ and $r$ denote the satellite and receiver, respectively; ${\rho }_r^s$ is the geometric distance between the satellite and the receiver (m); $t_r$ and $t^s$ are the receiver and satellite clock (m), respectively; $T_r$ is the slant tropospheric delay (m); $b_{r,i}$ and $b_i^s$ are the code hardware delay of the receiver and the satellite (m), respectively; $B_{r,i}$ and $B_i^s$ are the phase hardware delay of the receiver and the satellite (m), respectively; $f_i$ is the frequency; ${\lambda }_{IF}$ is the IF wavelength; $N_{r,IF}^s$ is the IF ambiguity (cycle); $e_{r,12}^s$ and ${\epsilon }_{r,12}^s$ 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,22,23]. Therefore, the relationship between the clock offset products ($t_{IF}^s$) 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)$$

(3)

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

$$\{\begin{array}{l}{p}_{r,IF}^s={\mu }_r^s·x+t_{r,IF}-t_{IF}^s+T_r+e_{r,IF}^s\hfill \\ l_{r,IF}^s={\mu }_r^s·x+t_{r,IF}-t_{IF}^s+T_r+{\lambda }_{IF}{\widehat{N}}_{r,IF}^s+{\epsilon }_{r,IF}^s\hfill \end{array}$$

(4)

where $p_{r,IF}^s$ and $l_{r,IF}^s$ denote the observed minus computed (OMC) IF observations of the code and carrier phase, respectively; ${\mu }_r^s$ 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+\left({\alpha }_{12}{b}_{r,1}+{\beta }_{12}{b}_{r,2}\right)$$

(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:

$${\widehat{N}}_{r,IF}^s=N_{r,IF}^s+\left[{\alpha }_{12}\left(B_{r,1}-B_1^s+b_1^s-b_{r,1}\right)+{\beta }_{12}\left(B_{r,2}-B_2^s+b_2^s-b_{r,2}\right)\right]/{\lambda }_{IF}$$

(6)

where ${\widehat{N}}_{r,IF}^s$ and $N_{r,IF}^s$ represent the float and integer ambiguities, respectively. The integer feature of the IF ambiguity ($N_{r,IF}^s$) 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 (${\overline{b}}_1^s$ and ${\overline{b}}_2^s$) and phase bias (${\overline{B}}_1^s$ and ${\overline{B}}_2^s$) products [27]:

$$\left(\begin{array}{c}{\overline{b}}_1^s\\ {\overline{b}}_2^s\end{array}\right)=\left(\begin{array}{c}{\beta }_{12}\left(b_1^s-b_2^s\right)\\ -{\alpha }_{12}\left(b_1^s-b_2^s\right)\end{array}\right),\left(\begin{array}{c}{\overline{B}}_1^s\\ {\overline{B}}_2^s\end{array}\right)=\left(\begin{array}{c}{\alpha }_{12}{b}_1^s+{\beta }_{12}{b}_2^s-B_1^s\\ {\alpha }_{12}{b}_1^s+{\beta }_{12}{b}_2^s-B_2^s\end{array}\right)$$

(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{array}{l}{\overline{b}}_{IF}^s={\alpha }_{12}{\overline{b}}_1^s+{\beta }_{12}{\overline{b}}_2^s=0\hfill \\ {\overline{B}}_{IF}^s={\alpha }_{12}{\overline{B}}_1^s+{\beta }_{12}{\overline{B}}_2^s={\alpha }_{12}{b}_1^s+{\beta }_{12}{b}_2^s-{\alpha }_{12}{B}_1^s-{\beta }_{12}{B}_2^s\hfill \end{array}$$

(8)

When ${\overline{b}}_{IF}^s$ and ${\overline{B}}_{IF}^s$ are combined with Equation (1), the new IF ambiguity can be stated as follows:

$${\widehat{N}}_{r,IF}^s=N_{r,IF}^s+\left[{\alpha }_{12}\left(B_{r,1}-b_{r,1}\right)+{\beta }_{12}\left(B_{r,2}-b_{r,2}\right)\right]/{\lambda }_{IF}$$

(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:

$${\widehat{N}}_{r,IF}^s=\left(\frac{f_2}{f_1+f_2}{\lambda }_{WL}{N}_{r,WL}^s+{\lambda }_{NL}{\widehat{N}}_{r,NL}^s\right)/{\lambda }_{IF}$$

(10)

where $N_{r,WL}^s$ is the integer WL ambiguity; ${\widehat{N}}_{r,NL}^s$ is the float NL ambiguity; and ${\lambda }_{WL}$ and ${\lambda }_{NL}$ are the wavelengths of the WL and NL ambiguities, respectively, which can be formulated as follows:

$$\{\begin{array}{c}{\lambda }_{WL}=\frac{C}{f_1-f_2}\\ {\lambda }_{NL}=\frac{C}{f_1+f_2}\end{array}$$

(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{array}{l}{\widehat{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 }_{WL}\left({\lambda }_2+{\lambda }_1\right)}\right)\hfill \\ {\widehat{N}}_{r,NL}^s=\frac{{\lambda }_{IF}{\widehat{N}}_{r,\mathrm{IF}}^s}{{\lambda }_{NL}}-\frac{{\lambda }_1{N}_{r,\mathrm{WL}}^s}{{\lambda }_2-{\lambda }_1}\hfill \end{array}$$

(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 multi-epoch 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:

$${\widehat{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·(P_{r,1}^s+z_{r,1}^s)+{\lambda }_1·(P_{r,2}^s+z_{r,2}^s)}{\left({\lambda }_2+{\lambda }_1\right)·{\lambda }_{WL}}\right)$$

(13)

$$\{\begin{array}{l}{z}_{r,1}^s=z_{r,1}·\sin {\theta }_r^s-z_1^s\hfill \\ z_{r,2}^s=z_{r,2}·\sin {\theta }_r^s-z_2^s\hfill \end{array}$$

(14)

where $z_{r,1}$ and $z_{r,2}$ are the vertical phase center offsets of the receiver antenna for the two frequencies, and similarly, $z_1^s$ and $z_2^s$ are those of the satellite antenna. ${\theta }_r^s$ 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. OSB product information provided by CNES.

System	Frequency Number	Frequency Type	Code Bias	Phase Bias
GPS	1	L1	C1C C1P C1W	L1C
	2	L2	C2C C2S C2L C2X C2W	L2W
GLONASS	1	G1	C1C C1P
	2	G2	C2C C2P
Galileo	1	E1	C1C	L1C
	2	E5a	C5Q	L5Q
BDS	1	B1I	C2I	L2I
	2	B3I	C6I	L6I

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.

. It should be noted that the code bias was stable during the day [31]; therefore, the code bias was not analyzed.

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.

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. Parameters for PPP configurations.

Parameter	Configurations
Estimator	Extended Kalman filter
Observations model	Ionosphere-free combinations
Frequency	GPS: L1/L2 Galileo: E1/E5a BDS: B1I/B3I
Observation weighting	0.3 m and 0.003 m for GPS/Galileo/BDS code and phase, respectively
Sampling rate and cutoff angle	30 s and 7°
Orbits and clocks	Real-time products from CNES
Code and phase biases	Corrected with the OSB products from CNES
Station coordinate	Kinematic mode: white noise Static mode: constant
Receiver clock	Estimated as white noise for each system
Tropospheric delay	A priori troposphere delay: Saastamoinen model [36] Zenith wet tropospheric delay: estimated as a random walk
AR validation	Bootstrapping success rate threshold: 95% Ratio test threshold: 2.0

.

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.

The mean DA of each satellite for each system is displayed in

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

System	Frequency 1	Frequency 2
GPS	97.6%	97.6%
Galileo	98.8%	98.8%
BDS	49.5%	49.5%

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

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.

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. Overall MAX and STD of GPS, Galileo, and BDS satellites.

System (Units: Cycles)	Frequency 1		Frequency 2
	MAX	STD	MAX	STD
GPS	0.045	0.012	0.045	0.012
Galileo	0.081	0.028	0.083	0.028
BDS	0.292	0.085	0.293	0.085

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

3.3. Distribution of Ambiguity Residuals Using OSB Products

One of the criterion indices for evaluating the quality of AR products is the distribution of ambiguity residuals [40,41]. The ambiguity should be close to an integer after applying the OSB product to conduct the SD between satellites; the rounding method can be utilized to determine the closest integers. The residuals denote the difference between the float values and the closest integers for the SD ambiguities, and these are specifically expressed as follows:

$$Residual=\Delta {\widehat{N}}_r^s-rounding\left(\Delta {\widehat{N}}_r^s\right)$$

(15)

where $\Delta {\widehat{N}}_r^s$ denotes the SD float ambiguity. The converged float ambiguities were used in these statistics to reduce the negative effects of other errors. Furthermore, the BDS system was separated into BDS-2 and BDS-3 for statistical analysis in this subsection.

The WL and NL residual distributions for GPS, Galileo, BDS-2, and BDS-3 were plotted using the same 31-day data; they are illustrated in Figure 4, while the specific values are shown in

Table 5. WL and NL residual statistic results for the four systems.

System	WL		NL
	With ± 0.15	With ± 0.25	With ± 0.15	With ± 0.25
GPS	82.8%	95.8%	80.4%	92.1%
Galileo	84.3%	92.3%	74.3%	88.7%
BDS-2	85.7%	95.9%	61.5%	79.0%
BDS-3	89.3%	97.2%	39.8%	60.1%

. It can be seen from the figure that besides the NL residual distribution of the BDS-3, the WL and NL residual distributions for the four systems followed a normal distribution. Regarding the WL residuals for the four systems, the residual percentages within ±0.15 cycles were 83%, 84%, 86%, and 90%, respectively, and those within ±0.25 cycles were greater than 92%. As for the NL residuals, the residual percentages within ±0.25 for the GPS and Galileo systems were 92% and 89%, which were considerably better than those of 79% and 60% for BDS-2 and BDS-3, respectively. Inferior satellite orbit and clock offset products, model residuals, and other factors can be blamed for the poor NL residuals of BDS-2 and BDS-3. It is evident that the ambiguity residuals of GPS showed the best performance, while those of Galileo were slightly worse than those of GPS; those of BDS-2 and BDS-3 were the worst, which is consistent with the quality analysis of the OSB products. Considering the DA, MAX, STD, and ambiguity residuals, it can be concluded that the positioning accuracy and convergence time will be impacted by BDS when performing multi-GNSS PPP-AR based on the OSB products.

3.4. Performance Analysis of PPP-AR

The convergence of real-time positioning is defined as having positioning errors in the east (E), north (N), and up (U) components smaller than 10 cm for 10 consecutive epochs. After convergence, the root mean square (RMS) of the ENU components represents the positioning accuracy, where the reference coordinates for each station are obtained from the IGS weekly solutions.

The average kinematic positioning accuracy results for the three combinations at 14 typical stations during the 31-day period are shown in Figure 5, where the float and fixed solutions are represented in the figure by the white diagonal and solid lines, respectively. The positioning accuracy of PPP-AR, especially for GPS-only solutions, is greatly increased when compared to that of float solutions, as can be seen in the figure. The mean positioning accuracy results of the kinematic and static PPP solutions for the three combinations were calculated using data from the 90 stations over 31 days to further reflect the improvement of PPP-AR in positioning accuracy; the results are given in

Table 6. Statistical results for the kinematic and static position accuracy of the three combinations.

Mode	System	Float (cm)			Fixed (cm)			Fixed Rate
		E	N	U	E	N	U
Kinematic	G	2.44	2.03	4.13	1.15 (53%)	1.38 (32%)	3.14 (24%)	97.8%
	GE	1.90	1.64	3.46	1.06 (44%)	1.27 (23%)	2.85 (17%)	98.5%
	GEC	1.89	1.64	3.43	1.06 (44%)	1.27 (22%)	2.85 (17%)	98.3%
Static	G	1.42	1.14	1.73	0.79 (44%)	0.95 (16%)	1.48 (15%)	98.8%
	GE	1.18	1.09	1.61	0.78 (34%)	0.94 (14%)	1.42 (12%)	98.7%
	GEC	1.18	1.09	1.59	0.77 (34%)	0.94 (14%)	1.41 (11%)	98.3%

.

When comparing the float solution in the GPS-only, the kinematic positioning accuracy of PPP-AR improved by 53%, 32%, and 24% (from 2.44, 2.03, and 4.13 cm to 1.15, 1.38, and 3.14 cm) for the E, N, and U components, respectively; on the other hand, it improved by 44%, 16%, and 15% (from 1.42, 1.14, and 1.73 cm to 0.79, 0.95, and 1.48 cm) for the E, N, and U components, respectively, in the static mode. With respect to the GPS+Galileo solution, the kinematic and static fixed solutions for the E, N, and U components improved from 1.90 to 1.06, 1.64 to 1.27, and 3.46 to 2.85 cm, and from 1.18 to 0.78, 1.09 to 0.94, and 1.61 to 1.42 cm, respectively. As demonstrated above, PPP-AR enhanced the positioning accuracy, most notably for the E component, while the GPS+Galileo solution improved both the float and fixed solutions to a certain extent in comparison to the GPS-only solution, particularly in the kinematic mode. However, the improvement in positioning accuracy was limited after the BDS satellites were involved in the PPP solutions; this can be attributed to the inferior quality of the DA of the BDS OSB products, which resulted in the participation of only a few BDS satellites in positioning. Meanwhile, the real-time satellite orbit and clock offset products for the BDS satellites were inaccurate, and the positioning accuracy could not be effectively improved under equal weight processing for all systems in this study. To summarize, in kinematic mode, the positioning accuracy of the float solutions after convergence reached about 2, 2, and 4 cm for the E, N, and U components, respectively, while the fixed solutions further improved the positioning accuracy to 1, 1, and 3 cm for the E, N, and U components, respectively. In static mode, the fixed solution improved the float solution from 1~2 cm to 7~9 mm for the E and N components, while the U component was greater than 1.5 cm. The mean ambiguity fixed rates for the kinematic and static modes of the three combinations were more than 97%.

Using the same data from 90 stations over 31 days, the kinematic and static convergence time results for the three combinations were counted using a frequency histogram, where the white diagonal and solid lines represent the convergence times for the float and fixed solutions, respectively. Each panel in Figure 6 represents the convergence time of the float solutions after promotion to fixed solutions; they are divided into the mean convergence time, the 10 min convergence ratio, the 20 min convergence ratio, and the 30 min convergence ratio. Since 10 epochs were used for MW smoothing in this study, the statistical analysis of convergence time was started after 5 min.

The mean convergence time of the kinematic and static float solutions for the GPS-only, GPS+Galileo, and GPS+Galileo+BDS solutions were 28.8, 19.7, and 20.4 min and 18.9, 14.7, and 15.0 min, respectively, while the results of the fixed solutions were 16.8, 9.6, and 9.89 min and 11.4, 8.0, and 8.1 min, respectively. The figure shows that the portion of each convergence period for the GPS+Galileo solution was higher than that of the GPS-only solution, especially at 10 min. This is because the satellite space geometry configuration is improved when more satellites are involved in PPP-AR. However, the convergence time of the GPS+Galileo and GPS+Galileo+BDS solutions was almost the same because the BDS real-time products were of poor accuracy and were frequently missing. To summarize, the convergence time of the fixed solution was much faster than that of the float solution for each convergence period, which demonstrates the effectiveness of PPP-AR in increasing the speed of convergence; in addition, the multi-GNSS PPP-AR further hastened the convergence time.

3.5. Dealing with the Missing Phase Bias

It should be pointed out that during the statistical analysis of positioning accuracy and convergence time, we found a problem where all stations re-converged in some periods—as shown in Figure 7, where the float and fixed solutions are represented in black and red, respectively—which greatly reduced the reliability of real-time PPP-AR services.

After the data analysis, it was found that this was mainly due to the missing real-time OSB products. Therefore, a polynomial fitting method was proposed to compensate for the temporary absence of the OSB products, which utilized the previous data in order to guarantee that the positioning result was reliable. The polynomial fitting method can be applied to fit an $n$-order polynomial based on $k$ known values, allowing for the value at the next epoch to be predicted, which can be expressed as follows:

$${\phi }_i=a_0+a_1\left(t_i-t_0\right)+a_2{\left(t_i-t_0\right)}^2+\dots +a_n{\left(t_i-t_0\right)}^n,\left(i=1,2,3,\dots ,k;k>n+1\right)$$

(16)

where $n$ represents the fitting order (second-order fitting was adopted in this study); $a_i$ denotes the $i$-order fitting coefficient; $t_0$ and $t_i$ represent the time of fitting and interpolated values, respectively; and ${\phi }_i$ represents the interpolated values.

By substituting the $k$ values into Equation (16) to perform polynomial fitting and using the least square method to obtain the coefficient of the polynomial, the polynomial can be employed to predict the missing value. It should be noted that the principle of the OSB and UPD is the same, and they can be converted into each other [42,43]. The missing phase bias products were only predicted within 15 min in this study because it has been demonstrated that the NL UPD products are stable within 15 min [2].

For instance, on DOY 138 in 2021, the phase bias products of all satellites were missing from 4:00:30 to 4:01:30 and from 11:43:00 to 11:44:30; the positioning errors at ALIC are plotted in Figure 8, where the upper figure represents the original results, and the lower figure represents the results using the proposed method. In comparison to the original positioning results, the compensated positioning results showed restored reliability while simultaneously avoiding re-convergence. Of course, the positioning accuracy and ambiguity fixed rate were also improved. It should be mentioned that the performance analysis in Section 3.4 was based on the method in Section 3.5.

4. Discussion

In this study, it can be seen that the BDS system still plays an auxiliary role in multi-GNSS PPP, although the visible BDS satellite number is very considerable. The main factors affecting the PPP performance are the accuracy and stability of precise orbit, clock, and OSB products. In fact, the quality of the OSB products depends on the upstream orbit and clock products. The average influence of GFZ real-time orbit error on NL UPD was from 0.06 to 0.78 cycles for four different GNSS systems [44]. In terms of CNES real-time orbit products, the 3D orbit RMS error is typically 5, 10, 18, 18, and 36 cm for GPS, GLONASS, Galileo, BDS MEO, and IGSO satellites, respectively [45]; the BDS orbit error was obviously higher than that of GPS. As a result, the effective solution for improving the BDS PPP ambiguity fixed rate is to reduce the effects of orbital errors. Estimating OSB products that can compensate for the BDS real-time orbit error may be a feasible solution in our future work.

Moreover, the stochastic model setting is also a vital point for multi-GNSS positioning as the errors of precise products vary in different GNSS systems or in different satellites ofthe same system. The errors of the CNES real-time orbit, clock, and OSB products were unknown, and equal weight was thus set for all satellite systems in this study. Therefore, due to the low quality of the BDS real-time products currently provided by CNES, the positioning accuracy became even worse after jointing the BDS into GPS PPP. Consequently, the integrity monitoring of satellite precise products is very important for a reliable real-time PPP service [46,47]. In addition, the post-store real-time precise products from CNES were used in this study, which prevented the possible communication delay in real-time applications. The extrapolated orbit and clock offset should be used when a communication delay occurs, which would allow the performance of PPP to drop rapidly. However, the extrapolation error will be different for different GNSS systems and different satellites [48]. This should also be considered in the stochastic model setting for multi-GNSS real-time PPP applications.

5. Conclusions

As a key product for real-time PPP-AR, real-time OSB products are currently broadcasted by CNES. Based on the dual-frequency ionospheric-free observation model, the integer feature of ambiguity was recovered by using the OSB products and PPP-AR could be performed. The quality of the OSB products was analyzed, and the performance of the dual-frequency PPP-AR was evaluated. The experimental results showed that:

(1) The mean DA of the GPS and Galileo satellites was above 97% according to the 31-day statistical analysis of the OSB products, whereas that of BDS was less than 60%. The mean MAX and STD of the GPS, Galileo, and BDS satellites were 0.045 and 0.012; 0.081 and 0.028; and 0.292 and 0.085 cycles, respectively. By analyzing the residual distribution of the OSB products, the WL residual percentages within ±0.25 cycles for the GPS, Galileo, BDS-2, and BDS-3 satellites were found to be more than 92%, while the NL residual percentages within ±0.25 cycles for the four systems were 92%, 89%, 79%, and 60%, respectively. Considering the DA, MAX, STD, and ambiguity residuals, we concluded that the reliability of PPP-AR, including the BDS system, may not be guaranteed using the OSB products from CNES, especially for BDS-3. In fact, the experiment described in this paper also confirmed that the advancement of positioning accuracy was constrained after BDS satellites had been incorporated into multi-GNSS PPP-AR.

(2) Due to the accuracy restrictions of real-time products, real-time PPP-AR cannot greatly increase positioning accuracy; however, it can dramatically speed up convergence time. In the kinematic mode, the GPS+Galileo fixed solution reduced the mean convergence time from 19.7 to 9.6 min and increased the convergence ratio within 10 min from 22.2% to 70.2% when compared to the float solution. After convergence, the positioning accuracy of the GPS+Galileo float solution reached about 2, 2, and 4 cm for the E, N, and U components, respectively, while that of the GPS+Galileo fixed solution was further enhanced to 1, 1, and 3 cm for the E, N, and U components, respectively. In the static mode, the overall fixed solution improved the float solution from 1~2 cm to 7~9 mm for the E and N components, and the U component was better than 1.5 cm. The mean convergence times of the GPS-only and GPS+Galileo solutions were 11.4 min and 8.0 min, respectively. The mean ambiguity fixed rates for the kinematic and static modes for the three combinations were more than 97%. In conclusion, the GPS+Galileo solution showed the best performance in the kinematic mode when using the OSB products from CNES, while the GPS-only solution can achieve almost the same performance as the GPS+Galileo solution in the static mode.

(3) A prediction method was proposed to compensate for the short-term missing OSB products, and our experiment showed that the positioning results were more reliable and avoided re-convergence.