*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 \text{\textquotedblleft}^s\_r-rounding \left(\Delta \text{\textquotedblright}^s\_r\right) \tag{15}$$

where Δ*N*ˆ *<sup>s</sup> <sup>r</sup>* 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. 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.

**Figure 4.** WL (**top row**) and NL (**bottom row**) residual distributions for GPS (**red left column**), Galileo (**green left column**), BDS-2 (**light blue right column**), and BDS-3 (**navy blue right column**).


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

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

**Figure 5.** Statistical diagram of mean 31-day kinematic results from 14 stations for the GPS-only (**top**), GPS+Galileo (**middle**), and GPS+Galileo+BDS (**bottom**) solutions (white diagonal line—float solution; solid line—fixed solution).


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

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.

**Figure 6.** Mean convergence time of the kinematic (**top row**) and static (**bottom row**) solutions for the GPS-only (**red left column**), GPS+Galileo (**green middle column**), and GPS+Galileo+BDS (**blue right column**) systems (white diagonal line—float solution; solid line—fixed solution).

#### *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:

$$\varphi\_i = a\_0 + a\_1(t\_i - t\_0) + a\_2(t\_i - t\_0)^2 + \dots + a\_n(t\_i - t\_0)^n, \ (i = 1, 2, 3, \dots, k; k > n+1) \tag{16}$$

where *n* represents the fitting order (second-order fitting was adopted in this study); *ai* denotes the *i*-order fitting coefficient; *t*<sup>0</sup> and *ti* represent the time of fitting and interpolated values, respectively; and *ϕ<sup>i</sup>* represents the interpolated values.

**Figure 7.** The kinematic positioning errors at the same 14 stations in Figure 5 for the GPS-only system from DOY 121 to 151 in 2021 (black—float solution; red—fixed solution).

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.

**Figure 8.** Kinematic positioning errors at ALIC for GPS-only on DOY 138 in 2021.

#### **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 realtime 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.
