Research on the Real-Time Ambiguity Resolution Algorithm of GPS/Galileo/BDS Based on CNES Real-Time Products

Abstract

Real-Time (RT) Precise Point Positioning (PPP) uses precise satellite orbits and clock corrections, and employs a separate receiver for positioning. With the growing demand, RT PPP is becoming an increasingly popular research topic. The ambiguity resolution (AR) can significantly improve the positioning accuracy and convergence time of PPP, so it is essential to study PPPAR in RT mode. In this paper, 37 MGEX stations from around the world are chosen, and the RT orbit, clock, and phase biases products broadcast by the Centre National d’Etudes Spatiales (CNES) are applied to PPPAR. Additionally, the residuals of the RT phase biases products, convergence time, and positioning accuracy are investigated. The results indicate that GPS products have the best quality of AR, with wide-lane (WL) and narrow-lane (NL) residuals of 98.9% and 95.3%, respectively, within ±0.25 cycles. Within ±0.25 cycles, the WL and NL residuals of the Galileo are 98.2% and 94.3%, respectively. Within ±0.25 cycles, the (Beidou Navigation Satellite System) BDS has a poor quality of AR, with WL and NL residuals of 97.3% and 73.1%, respectively. Due to the poor quality of the BDS AR, the convergence time of the BDS is not calculated in this paper. The convergence time of other systems is significantly reduced after AR processing, and the convergence time of the GPS/Galileo combination is the fastest, being 17.14 min in kinematic mode and only 11.85 min in static mode. The positioning accuracy of the GPS, Galileo, GPS/Galileo, and GPS/Galileo/BDS in the E and U directions is significantly improved after PPPAR.

1. Introduction

Precise Point Positioning (PPP) is a technology that does not require a reference station, has a flexible mode of operation, low cost, and high accuracy, and can provide users with centimeter- or even millimeter-level positioning services. As a result, PPP has a wide range of applications [1,2,3,4,5]. However, the release of post-PPP products suffers from a considerable lag, and it is challenging to meet the needs of RT users. Therefore, Muellerschoen et al. of JPL were the first to propose Real-Time (RT) PPP [6]. In recent years, an increasing number of users have switched from post-PPP to RT PPP [7]. The international GNSS Service (IGS) launched the RT project plan officially in 2007 and began to provide RT precision products officially in April 2013 [8].

The RT orbit and clock products are generated via the processing of data from 100 globally dispersed tracking stations [9,10,11]. Yao et al. [12] analyzed the RT and the final products of the GPS, and found that the difference in accuracy between the two was only 1.2 cm. Zhang et al. [13] examined the impact of RT products from various analysis centers on GPS data processing and found that CLK90 products are superior. Liu et al. [14] investigated a GPS RT PPP algorithm and model. Wang et al. [15] analyzed the RT PPP performance of the BDS. In order to prevent quality issues in the RT orbit and clock of satellites in RT mode, Ji et al. [16] monitored the orbit and clock in RT, and effectively detected problematic satellites.

Many scholars have conducted studies on multi-system RT PPP. However, de Bakker et al. found that the addition of GPS/GLONASS combined processing to the GLONASS did not improve the positioning performance in most cases [17]. Wang et al. [18] conducted kinematic experiments on RT PPP using BDS observation data. Abdi et al. [19] compared the positioning accuracy of BDS/GPS using RT products after incorporating BDS observation data, revealing an average decrease from 1.01/0.58 m to 0.70/0.50 m. Liu et al. [20] proposed a novel weighting method to constrain the ionospheric delay in RT GPS/Galileo-combined PPP. Based on the BDS PPP-B2b service, the BDS/GPS in the RT PPP kinematic mode can achieve a positioning accuracy of 5.9 cm, 3.6 cm, and 9.4 cm in the East (E)/North (N)/Up (U) directions, respectively [21]. Currently, CNES provides real-time products that support all systems. Kazmierski et al. [22] discovered that inter-system weighting can improve the coordinate repeatability and decrease the convergence time. In addition, Kazmierski et al. [23] compared the RT positioning performance of a four-system PPP model. However, when adding GPS or BDS, unexpected drops in the positioning accuracy were observed. Based on Ionospheric-Free (IF) observations [24], the combined RT PPP kinematic model employing GPS/BDS/Galileo/GLONASS can achieve E/N/U direction accuracies of 1.6 cm, 1.2 cm, and 3.4 cm.

Traditional RT PPP technology has faced challenges in meeting the current demand for high-precision RT positioning performance, thus necessitating research into AR technology. Gabor [25] first proposed a method to correct satellite-end wide-lane (WL) and narrow-lane (NL) Fractional Cycle Biases (FCB) by fixing the interstellar single-difference PPP ambiguity. Ge [26] uses the single-difference between-satellites model to fix the ambiguity. In contrast, Collins [27] proposed a clock decoupling model that uses distinct satellite clock corrections to restore the integer characteristics of zero-difference ambiguity through the pseudo-range and carrier phase, thereby successfully resolving the zero-difference ambiguity. Laurichesse [28] proposed estimating the integer satellite clock by recovering the integer characteristics of zero-difference phase ambiguity, and achieving PPP with fixing zero-difference ambiguity. Relevant studies indicate that these three methods are equivalent because they absorb the same error term but express it differently [29,30,31,32]. Based on Ge’s method, Geng et al. [33] explored single-difference narrow-lane using the Least-squares AMBiguity Decorrelation Adjustment (LAMBDA) technique. Zhang and Li conducted an analysis of the PPPAR in a zero-difference model. By incorporating the station FCB as a reference, they were able to restore single-difference FCB to zero-difference FCB, thereby achieving a solution for PPPAR in the zero-difference model. This successful application of GPS PPP ambiguity fixing was utilized for low-orbit satellite orbit determination [34,35,36]. Laurichesse et al. proposed a method to mitigate NL FCB by employing the satellite clock correction, which was implemented using the CNES analysis center in order to generate GRGS products [37]. Liu et al. [38] employed the integer-phase clock provided by the CNES in France to achieve PPPAR. Numerous kinematic PPP solution tests demonstrate that fixed PPP has a faster convergence rate and superior positioning accuracy and stability than float PPP. Laurichesse et al. [28] described the calculation of the CNES phase biases products. Liu et al. [39] evaluated the positioning performance of the CNES real-time phase biases product for GPS and Galileo in RT PPPAR. Gazzino et al. [40] conducted a thorough evaluation of the use and quality of phase biases products offered by the CNES. Du et al. [41] proposed a prediction method to improve the reliability of PPP-AR in the absence of real-time OSB products for a short time, and improve the positioning accuracy and the ambiguity fixed rate. Li et al. [7] presented a multi-system PPPAR solution comprising GPS, GLONASS, BDS, and Galileo. This joint system optimizes satellite space geometry to achieve greater redundancy and model strength than a single system. In comparison with a single GPS solution, the multi-frequency multi-system approach can enhance the AR rate [42,43,44]. Geng et al. [45] proposed a method for computing observable-specific signal bias (OSB) and conducted experiments using OSB products. Li et al. [46] investigated the effects of an IF model and an undifferenced and uncombined (UDUC) model on ambiguity resolution; their findings demonstrated that under static and kinematic modes, the positioning accuracy increased by 45% with IF model implementation, while the UDUC model implementation resulted in a 40% increase.

It is essential to conduct RT research on multi-system PPPAR in order to attain a higher precision positioning performance and a quicker convergence time in RT mode. However, there is a scarcity of studies focusing on the RT multi-system PPPAR technique. The emergence of multi-system RT phase biases products from CNES institutions in recent years has sparked a surge in research interest in multi-system RT PPPAR. With the completion and development of the BDS system in our country, more and more users now use the BDS system for positioning, but most of the current research on RT PPPAR is based on GPS and the Galileo system. Therefore, the BDS system is added in this paper, and all the systems are analyzed comprehensively. Due to the small number of BDS-2 satellites, the impact on positioning results is not obvious, so only BDS-3 satellites are used in this paper for RT PPP and RT PPPAR. This article can also provide an important data reference for other scholars’ research in this field. Furthermore, this paper adopts the IF combination model. The ambiguities of GPS, BDS and Galileo are fixed, respectively, in real-time mode. The residuals of WL and NL, the convergence time and the positioning accuracy of each system are analyzed in both the static and kinematic mode.

2. Materials and Methods

2.1. PPP Zero-Difference Observation Equation

The observation of Global navigation satellite systems (GNSS) is mainly pseudo-range (P) and phase (L), which can be written as follows (the unit of measurement is in meters):

$$P_{r,j}^s={\rho }_r^s+c(d{t}_r^{}-d{t}^s)+T_r^s+I_{r,j}^s+b_{r,j}^{}-b_j^s+e_{r,j}^s$$

(1)

$$L_{r,j}^s={\rho }_r^s+c(d{t}_r^{}-d{t}^s)+T_r^s-I_{r,j}^s+{\lambda }_j^{}(N_{r,j}^s+B_{r,j}^{}-B_j^s)+{\epsilon }_{r,j}^s$$

(2)

where $s$ is the satellite, $j$ is the frequency, $r$ represents the receiver, $\rho$ is the distance between the satellite and the receiver, $c$ represents the speed of light, $d{t}_r$ and $d{t}^s$ denote the receiver and satellite clock, $T_r^s$ is the inclined path tropospheric delay, $I_{r,j}^s$ represents the j-frequency inclined path first-order ionospheric delay, and $N_{r,j}^s$ is the phase ambiguity. $B_{r,j}$ and $B_j^s$ are the phase hardware delay of the j-frequency receiver and the satellite, respectively [47]. ${\lambda }_j$ is the j-frequency wavelength, $b_{r,j}$ is the pseudo-range hardware delay between the j-frequency receiver antenna and the signal correlator, and $b_j^s$ is the pseudo-range hardware delay between the j-frequency satellite end signal transmitter and the satellite antenna. $e$ and $\epsilon$ are the pseudo-range and carrier phase measurement errors between the satellite and the receiver.

2.2. Ionospheric-Free Combination with Ambiguity Fixing

In PPP, the IF model is the most frequently used function model. The observed combination equation is as follows [1,2]:

$$P_{r,IF}^s={\rho }_r^s+c(d{t}_r^{}-d{t}^s)+T_r^s+b_{r,IF}^{}-b_{IF}^s+e_{r,IF}^s$$

(3)

$$L_{r,IF}^s={\rho }_r^s+c(d{t}_r^{}-d{t}^s)+T_r^s+{\lambda }_{IF}^{}(N_{r,IF}^s+B_{r,IF}^{}-B_{IF}^s)+{\epsilon }_{r,IF}^s$$

(4)

Among them,

$$b_{r,IF}=(f_1^2{b}_1^s-f_2^2{b}_2^s)/(f_1^2-f_2^2)$$

(5)

$$N_{r,IF}^s=c(f_1^{}{N}_{r,1}^s-f_2^{}{N}_{r,2}^s)/(f_1^2-f_2^2)/{\lambda }_{IF}$$

(6)

$$B_{r,IF}^{}=c(f_1^{}{B}_{r,1}^{}-f_2^{}{B}_{r,2}^{})/(f_1^2-f_2^2)/{\lambda }_{IF}$$

(7)

$$B_{IF}^s=c(f_1^{}{B}_1^s-f_2^{}{B}_2^s)/(f_1^2-f_2^2)/{\lambda }_{IF}$$

(8)

where $f_1$ and $f_2$ denote the 1 and 2 frequency, respectively. When the RT PPP calculation is carried out, the satellite pseudo-range clock products released by IGS often use $cd{t}_{P_{IF}}^s=cd{t}^s+b_{IF}^s$, and the receiver pseudo-range clock is also estimated as $cd{t}_{r,P_{IF}}^s=cd{t}_r+b_{r,IF}^{}$. In addition, since the phase delay is related to the ambiguity parameter and we generally consider it to have high time stability, it can be absorbed when the ambiguity is resolved using a real-valued [29,48]. After clock correction, Equations (3) and (4) can be rewritten as follows:

$$P_{r,IF}^s={\rho }_r^s+cd{t}_{r,IF}-cd{t}_{P_{IF}}^s+T_r^s+e_{r,IF}^s$$

(9)

$$L_{r,IF}^s={\rho }_r^s+cd{t}_{r,IF}-cd{t}_{P_{IF}}^s+T_r^s+{\lambda }_{IF}{\overline{N}}_{r,IF}^s+{\epsilon }_{r.IF}^s$$

(10)

Among them,

$${\overline{N}}_{r,IF}^s=N_{r,IF}^s+d_{r,IF}-d_{IF}^s$$

(11)

$$d_{r,IF}^{}=B_{r,IF}-b_{r,IF}/{\lambda }_{IF}$$

(12)

$$d_{IF}^s=B_{IF}^s-b_{IF}^s/{\lambda }_{IF}$$

(13)

where ${\overline{N}}_{r,IF}^s$ is ionospheric-free ambiguity, it does not have integer characteristics, and the combination of the hardware delay and ambiguity parameters are linearly correlated, which is difficult to separate. In parameter estimation, the two parameters are typically combined into one parameter, and the real-valued solution is adopted. In fact, since the integer portion of the hardware delay has no effect on the integer characteristics of the ambiguity, the ambiguity parameter and the integer portion of the hardware delay can be considered equivalent, and only the fractional portion influence is considered. Therefore, as long as the fractional portion can be corrected in RT, the ambiguity in RT PPP can be fixed. Consequently, ionospheric-free ambiguity can be separated into the following forms [49]:

$${\lambda }_{IF}{\overline{N}}_{r,IF}^s=\left(\frac{{\lambda }_{NL}{f}_1^2}{f_1^2-f_2^2}{\overline{N}}_{r,NL}^s-\frac{{\lambda }_{WL}{f}_2^2}{f_1^2-f_2^2}{\overline{N}}_{r,WL}^s\right)$$

(14)

where ${\lambda }_{WL}$ and ${\lambda }_{NL}$ are the WL and NL wavelength, and ${\overline{N}}_{r,WL}^s$ and ${\overline{N}}_{r,NL}^s$ are the corresponding hardware delay. Equations (15) and (16) can be used to show that because the integer part of the hardware delay has no effect on the integer characteristic of the ambiguity unless otherwise specified, the integer ambiguity described below includes the integer part of the hardware delay.

$${\overline{N}}_{r,WL}^s=N_{r,WL}^s+b_{r,WL}-b_{WL}^s$$

(15)

$${\overline{N}}_{r,NL}^s=N_{r,NL}^s+b_{r,NL}-b_{NL}^s$$

(16)

where $N_{r,WL}^s$ and $N_{r,NL}^s$ are the WL and NL ambiguity, $b_{r,WL}$ and $b_{r,NL}$ are the corresponding hardware delays at the receiver, and $b_{WL}^s$ and $b_{NL}^s$ are the corresponding hardware delays at the satellite, respectively. Therefore, if we can obtain the RT wide and narrow-lane hardware delay or related products and fix them in real-time, followed by the ionospheric-free ambiguity in real-time, we can finally achieve RT PPP with ambiguity fixing.

Typically, the MW combination is used to resolve WL ambiguity [50,51], which is defined below:

$$L_{MW}^s=\frac{f_1{L}_1^s-f_2{L}_2^s}{f_1-f_2}-\frac{f_1{P}_1^s+f_2{P}_2^s}{f_1-f_2}={\lambda }_{WL}(N_{WL}^s+b_{r,WL}-b_{WL}^s)+{\epsilon }_{L_{MW}}^s$$

(17)

where $L_{MW}^s$ is the Melbourne–Wübbena (MW) observation quantity, and ${\epsilon }_{L_{MW}}^s$ represents the MW combined observation noise. The MW combination observation noise is too great, necessitating a multi-epoch smoothing process. After smoothing, the MW combination observation quantity is as follows:

$$\frac{{\overline{L}}_{MW}^s}{{\lambda }_{WL}}=\langle \frac{L_{MW}^s}{{\lambda }_{WL}}\rangle ={\tilde{N}}_{r,WL}^S=N_{r,WL}^S+b_{r,WL}-b_{WL}^s$$

(18)

where $\langle *\rangle$ is multi-epoch averaging. The actual process is as follows:

$${\tilde{N}}_{WL}^s(t)={\tilde{N}}_{WL}^s(t-1)+\frac{1}{n}\left[{\overline{N}}_{WL}^s(t)-{\tilde{N}}_{WL}^s(t-1)\right]$$

(19)

$${\sigma }_{{\tilde{N}}_{WL}^s(t)}^2={\sigma }_{{\tilde{N}}_{WL}^s(t-1)}^2+\frac{1}{n}\left[{\left({\overline{N}}_{WL}^s(t)-{\tilde{N}}_{WL}^s(t-1)\right)}^2-{\sigma }_{{\tilde{N}}_{WL}^s(t-1)}^2\right]$$

(20)

where ${\tilde{N}}_{WL}^s(t)$ is the mean value of WL ambiguity, n is the number of the continuous observation epoch, and ${\sigma }_{{\tilde{N}}_{WL}^s(t)}^2$ is the variance. After obtaining the average value, it can be resolved by correcting the product. Since it is the same for all satellites, the decimal part after adding the corrected products can be obtained by averaging the following:

$$b_{r,WL}=\frac{1}{m}{{\displaystyle \sum _{i=1}^m\left[(N_{r,WL}^s-b_{WL}^s)-{\lambda }_{WL}\langle \frac{N_{r,WL}^s-b_{WL}^s}{{\lambda }_{WL}}\rangle \right]}}_i$$

(21)

where m is the number of satellites. The integer property of WL ambiguity can be obtained by substituting the obtained mean value $b_{r,WL}$ back into Equation (18). To ensure that the obtained integer solution is correct and fixed, it is also necessary to perform a calculation using the following formula [52]:

$$\{\begin{array}{l}{\tilde{N}}_{r,WL}^s(t)-round({\tilde{N}}_{r,WL}^s)<{\delta }_1\\ P_0=1-{\displaystyle \sum _{i=1}^{\infty }\left[erfc\left(\frac{i-\left|{\tilde{N}}_{r,WL}^s-N_{r,WL}^s\right|}{\sqrt{2}{\sigma }_{{\tilde{N}}_{r,WL}^s}}\right)-erfc\left(\frac{i+\left|{\tilde{N}}_{r,WL}^s-N_{r,WL}^s\right|}{\sqrt{2}{\sigma }_{{\tilde{N}}_{r,WL}^s}}\right)\right]}<{\delta }_2\\ erfc(x)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_x^{\infty }{e}^{-t^2}}dt\end{array}$$

(22)

where $round(*)$ is the nearest integer function, $P_0$ is the fixation to the nearest integer probability, $N_{r,WL}^s$ represents the nearest full number of ${\tilde{N}}_{r,WL}^s$, and ${\delta }_1$ and ${\delta }_2$ are the fixed acceptance thresholds. Therefore, the closer the real-valued ambiguity solution is to the nearest integer, the smaller the variance in the real-valued ambiguity and the better the wide-lane fixing.

When the ionospheric-free ambiguity real-valued solution ${\overline{N}}_{r,IF}^s$ is known, and when the fixed WL ambiguity $N_{r,WL}^s$ is substituted into Equation (14), the real-valued solution of NL ambiguity ${\overline{N}}_{r,NL}^s$ is as follows:

$$\{\begin{array}{l}{\overline{N}}_{r,NL}^s=\frac{f_1+f_2}{f_1}{\overline{N}}_{r,IF}^s+\frac{f_2}{f_1-f_2}{N}_{r,WL}^s\\ Q_{{\overline{N}}_{r,NL}^s}=\frac{{\lambda }_1^2}{{\lambda }_n^2}{Q}_{{\overline{N}}_{r,IF}^s}\end{array}$$

(23)

where $Q_{{\overline{N}}_{r,NL}^s}$ and $Q_{{\overline{N}}_{r,IF}^s}$ are the variance–covariance matrices of NL ambiguity and ionospheric-free ambiguity, and ${\lambda }_n$ is the wavelength of the n-frequency. Because the NL wavelength is shorter, the optimal solution is obtained using the LAMBDA algorithm, and the integer characteristic of IF combination ambiguity is restored by fixing WL and NL ambiguity sequentially.

3. Experimental Results and Analysis

3.1. Experimental Data and Processing Scheme

In this study, 37 Multi-GNSS Experiment (MGEX) stations were randomly selected globally for RT PPP processing and RT PPPAR processing using the CNES real-time phase biases product. The phase biases products are broadcasted in real-time with an updating rate of 5 s on the mountpoint CLK93 from the CNES caster (ntrip.gnsslab.cn: 2101); these real-time products are also saved as files and available from the following website: (http://www.ppp-wizard.net/products/REAL_TIME/). The observation period for the experimental data ranged from the 150th to the 169th day of 2023. The precise coordinates of the stations are available in file ‘‘IGS0OPSSNX_2023DDD0000_01D_01D_CRD.SNX’’ provided by IGS; we compared the station coordinates calculated using different methods with the coordinates supplied via IGS weekly solutions. The filter was considered to have converged when kept within 10 cm for 10 consecutive epochs in both the static and kinematic modes. In terms of data quality control, we detected and labeled cycle slips using a combined MW+(Geometry-Free) GF approach [53,54]. In the experimental section, first, the ionosphere-free ambiguity was calculated. Then, it was decomposed into WL and NL ambiguity, and the WL ambiguity was corrected, smoothed and fixed. The floating-point single difference subtracted the WL deviation to obtain the floating-point NL. Finally, the NL ambiguity was fixed, the fixed floating-point ionospheric ambiguity was obtained, and then the equation constraint was solved. The well-known ratio test was used for ambiguity validation. In this study, the critical criterion for the ratio test was selected as 3. Figure 1 and

Table 1. Basic information of MGEX stations.

Station	Latitude	Longitude	Station	Latitude	Longitude
ALIC	−23.67	133.886	KERG	−49.351	70.256
ASCG	−7.916	−14.333	LAMA	53.892	20.67
BOGT	4.64	−74.081	LCK3	53.892	20.67
CEDU	−31.867	133.81	MARS	43.279	5.354
CORD	−31.528	−64.47	MAW1	−67.605	62.871
CZTG	−46.432	51.855	MIKL	46.973	31.973
DAV1	−68.577	77.973	MIZU	39.135	141.133
EBRE	40.821	0.492	NKLG	0.354	9.672
ENAO	39.091	−28.026	NIST	39.995	−105.263
GAMB	−23.13	−134.965	NLIB	41.772	−91.575
GAMG	35.59	127.92	POVE	−8.709	−63.896
GLSV	50.364	30.497	PTBB	52.296	10.46
GRAZ	47.067	15.493	ULAB	47.865	107.052
HOLB	50.64	−128.135	URUM	43.808	87.601
HUEG	47.834	7.596	USUD	36.133	138.362
INEG	21.856	−102.284	WIND	−22.575	17.089
IISC	13.021	77.57	XMIS	−10.45	105.689
IITK	26.521	80.232	YEL2	62.481	−114.481
JDPR	26.207	73.024

provide the specific global distribution of the 37 stations,

Table 2. PPP processing strategies.

Item	Setting
Observations	Pseudo-range and phase observations
Solution model	Static/Kinematic
Frequency	GPS: L1/L2; Galileo: E1/E5a; BDS: B1I/B3I
Orbits and clocks	CNES RT orbit and clock products
Elevation cut off	7°
Sampling offset	30 s
Phase windup	Phase polarization effects applied
Parameter estimation method	Kalman filtering
Inter-system bias	Estimated as a random walk
Receiver coordinates	Parameters estimation
Receiver clock error	Parameters estimation
Earth tides	IERS2010
Weighting scheme	Elevation dependent weight
Ionosphere	Ionosphere-free
A priori troposphere delay	Saastamoinen
Tropospheric mapping function	VMF1
Satellite phase center	igs20.atx
Phase ambiguity	WL + NL
Cycle slip method	GF + MW

details the specific PPP processing strategies.

According to Figure 2 and Figure 3, the number of satellites in Galileo is eight, with the PDOP close to 2.4. In contrast, GPS has an average satellite count of 10 with a PDOP below 2, whereas BDS has the lowest position dilution of precision (PDOP) consistently; this is below 1.2 in the Asia–Pacific region, where it possesses a significantly higher number of available satellites. This is up to 12 due to differences in its satellite constellation design, operating altitude, and orbit inclination parameters compared to the BDS configuration. Consequently, BDS surpasses other systems in the Asia–Pacific region in terms of its satellite quantity and spatial distribution [55].

3.2. Ambiguity Resolution Residuals Analysis

According to Figure 4, within ±0.25 cycles, GPS exhibits WL and NL residuals of 98.9% and 95.3%, respectively, while within ±0.15 cycles, these values are 87.1% and 83.5%, respectively. Moreover, the WL residuals for GPS within ±0.25 cycles are 98.2% and the NL residuals within ±0.25 cycles are 83.9%, while the NL residuals for GPS within ±0.15 cycles are observed to be approximately 81%. In general, the GPS residual distribution outperforms that of Galileo. This is due to the extensive RT tracking of GPS signals by receivers and the high quality of GPS orbit determination. However, due to the inadequate treatment of RT deviation products in BDS, both WL and NL residuals exhibit poor distribution characteristics. As shown in the figure, the WL residuals for BDS at the ±0.25 and ±0.15 cycles are 97.3% and 79.2%, respectively, while the NL residuals only reached 73.1% and 53.7%.

As shown in Figure 4, the RT PPP method exhibits the highest ambiguity fixing rate for GPS, while Galileo also performs well in terms of ambiguity fixing rate. In this regard, however, the BDS yields poor results. In this paper, therefore, only the RMS of the combined BDS/GPS/Galileo is considered, and the convergence time is ignored.

3.3. Convergence Performance Analysis

As depicted in Figure 5, the convergence time in static mode is the slowest for Galileo. Comparatively, the GPS has a faster convergence time than Galileo, whereas the GPS/Galileo has the quickest convergence time. Upon ambiguity resolution, significant improvements can be expected in the convergence time for all systems. In RT PPP, the convergence time for Galileo is 20.89 min. After fixing, the convergence time can reach 16.38 min. Prior to fixing, the GPS has a convergence time of 17.46 min, which decreases to 14.87 min after fixing. When the combined GPS/Galileo position is not fixed, the convergence time is 14.62 min, with 62% of stations converging between 10 and 20 min and less than 20% achieving fixation within 10 min. After ambiguity resolution, the convergence time improves to a remarkable 11.45 min, with over 40% of stations converging within 10 min. These results demonstrate that using a dual system combination observation mode significantly reduces the convergence time of positioning, while ambiguity fixing also plays a crucial role in enhancing station convergence.

As depicted in Figure 6, the Galileo convergence time is the slowest in kinematic mode. In RT PPP, it takes approximately 27.43 min to converge, which reduces to 23.14 min after ambiguity resolution. Following Galileo, the GPS has a longer convergence time of approximately 24.19 min without ambiguity resolution, which decreases to 21.68 min after the ambiguity is resolved. On the other hand, the GPS/Galileo demonstrates the quickest convergence time, at approximately 20.47 min without ambiguity resolution and 17.14 min after ambiguity resolution. After ambiguity resolution, it becomes evident that the proportion of stations converging within 10 increases significantly. This indicates that ambiguity resolution increases the convergence time in kinematic mode intuitively.

3.4. Positioning Accuracy Analysis

Figure 7 and Figure 8 depict the time series of the HOLB and GRAZ stations in kinematic mode on the 150th of 2023, displaying the RT PPP and RT PPPAR modes. The time-series of position accuracy in the static mode is not given due to its small fluctuations. As the convergence time has already been calculated, here only the positioning errors in three directions are counted.

As shown in Figure 7 and Figure 8, for GRAZ stations, it can be clearly observed that the error fluctuation of the Galileo system in the E and U directions is significantly higher than that of GPS, but that the positioning accuracy of Galileo is significantly improved after AR processing. For GPS/Galileo, it can be observed that the accuracy of RT PPPAR is significantly better than RT PPP in the E direction, which shows that AR is very helpful with regard to improving the positioning accuracy. For HOLB stations, the RT PPP AR position error using GPS and the Galileo system decreases in both the E and N directions, and the time-series fluctuations become smoother as well. Despite the insignificant or even opposite performance in the U direction, the holistic position accuracy is still improved. Moreover, the GPS/Galileo is greater than that of a single system; however, significant fluctuations still persist. However, when BDS is incorporated into the combined positioning of all three systems, both the E and N direction fluctuations are significantly reduced, indicating enhanced stability and precision in the horizontal components.

As shown in Figure 9a, 60% of the Galileo stations had improved localization accuracy in the E, N, U directions after PPPAR, and the U direction of the CZTG station had the best improvement effect, reaching 4.9 cm; in addition, the E, N, U directions of only three stations were worse after PPPAR, which may indicate that the same strategy is not practical for all the MGEX stations. It is clear in Figure 9b that, for the GPS system, 54% of the stations had an improved positioning accuracy in all three directions after PPPAR, while the N direction was significantly better than the E direction. A 2 cm improvement was observed in the E direction for the GAMB stations. In Figure 9c,d, it can be seen that regardless of whether GPS/Galileo or GPS/Galileo/BDS was used, the positioning accuracies in the E, N, U directions were significantly better than the GPS and Galileo systems, which is importantly related to the fact that multi-systems combine a large number of satellites and a good satellite spatial geometric configuration. However, the positioning accuracy of GPS/Galileo/BDS in RT PPPAR was worse than that of GPS/Galileo in the both the E direction and U direction, which is very much related to the poor solving of the BDS-phase deviation products.

From Figure 10a, it can be seen that for the Galileo system, all the stations improved in one or more directions in E, N and U. It is also clear that 78% of the stations improved in the E direction and 59% and 70% of the stations improved in the N and U directions. From Figure 10b, it can be seen that for the GPS system, after PPPAR, the improvement in the E direction was the best, followed by the U direction, and lastly the N direction, but that the positioning accuracy in the N direction was the best. The E direction of the HOLB station was improved by 0.7 cm, but the accuracy of the USUD station was reduced by 1.2 cm after PPPAR. From Figure 10c,d, it is evident that, no matter the RT PPP or RT PPPAR, the combination of multi-systems can still obviously improve the positional accuracy, although it is possible that the accuracy of individual stations is not as good as those without AR treatment after adding BDS. Although the accuracy of individual stations after PPPAR may not be as good as those without AR processing, the positioning accuracy after adding the BDS system is still better than that of GPS/Galileo.

Table 3. RMS values under different processing modes (unit: cm).

System	Kinematic						Static
	RT PPPAR			RT PPP			RT PPPAR			RT PPP
	E	N	U	E	N	U	E	N	U	E	N	U
Galileo	3.23	2.43	5.67	3.66	2.69	6.26	1.20	0.58	1.69	1.35	0.63	1.80
GPS	1.75	1.56	4.09	2.33	1.63	4.32	0.86	0.58	1.37	1.12	0.62	1.45
GPS/Galileo	1.05	0.95	3.08	1.77	1.22	3.51	0.75	0.50	1.17	1.00	0.54	1.32
GPS/Galileo/BDS	1.33	1.32	2.94	1.36	1.17	2.98	0.79	0.35	0.88	0.91	0.33	0.89

shows the average RMS value of each system under different processing modes. It is evident from the table that, regardless of whether the static or kinematic mode is used, after AR processing, the positioning accuracy is significantly improved. For a single system, GPS is superior to Galileo, due to its enhanced number of satellites and constellation configuration. At the same time, the GPS/Galileo combination improves the positioning results. The poor solution of the GPS/Galileo/BDS to resolve BDS phase products results in suboptimal positioning accuracy for the combined BDS system. However, the final average positioning accuracy is 1 cm or better in static mode and better than 3 cm in kinematic mode, which fully satisfies the requirements for precision in the daily RT mode.

4. Conclusions

To conduct a comprehensive analysis of the RT PPPAR algorithm, we employed the CNES’ precise clock, orbit, and phase biases products. Specifically, we selected the data for the 150th to 169th day of 2023 from 37 MGEX stations around the world. Subsequently, RT ambiguity fixing was performed separately for GPS, Galileo, GPS/Galileo, and GPS/Galileo/BDS, and the following conclusions were derived:

(1)

The RT phase biases products were analyzed. GPS has the best ambiguity fixing, followed by BDS and Galileo. Within a range of ±0.25 cycles, the GPS WL and NL residuals were 98.9% and 95.3%, respectively, while they were 87.1% and 83.5%, respectively, within a range of ±0.15 cycles. For the Galileo system, the WL residuals within ±0.25 cycles and ±0.15 cycles were 98.2% and 83.9%, while the residuals within ±0.15 cycles for NL were 81.4%; in addition, the corresponding value for within ±0.25 cycles is 94.3%. The distribution of the WL and NL residuals of BDS is poor. Within ±0.25 and ±0.15 cycles, the residuals of WL and NL were 97.3% and 79.2%, respectively, whereas the residuals of NL were only 73.1% and 53.7%.

(2)

PPPAR was effective in reducing convergence time. After applying AR processing, the GPS convergence time in static mode decreased from 17.46 min to 14.87 min, while it decreased from 20.89 min to 16.38 min for Galileo. At the same time, the GPS/Galileo decreased from 14.62 min to 11.85 min. In kinematic mode, the GPS decreased from 24.19 min to 21.68 min, while the average convergence of Galileo time decreased from 27.43 min to 23.14 min. After applying AR processing, the convergence time of the combined GPS/Galileo decreased from 20.47 min to 17.14 min, indicating a significant increase in the convergence speed.

(3)

After PPPAR processing, the improvement in the positioning accuracy in dynamic mode was more obvious than that in static mode. Meanwhile, the positioning effect of the Galileo system was the worst, reaching 6.2 cm in the U direction of dynamic RT PPP. The positioning effect of GPS was obviously better than Galileo, and the positioning accuracy of the multi-system combination was also better. Although the RT PPPAR results of individual stations are not as good as RT PPP after adding the BDS system, which is caused by the unsatisfactory solving effect of the BDS system, the positioning accuracy of the multi-system was still very good. Therefore, we should choose a combination of multiple systems for positioning in our daily real-time applications.

This paper analyzed the positioning performance of RT PPPAR with GPS/Galileo/BDS in detail. The results show that RT PPPAR can improve the position accuracy and shorten the convergence time of conventional RT PPP results. The research in this paper makes up for the lack of research on BDS RT PPPAR, and also shows that it can better meet higher demands for precision and timeliness. Further research on the RT PPPAR of the BDS system will be carried out in the future.