Multi-GNSS Combined Precise Point Positioning Using Additional Observations with Opposite Weight for Real-Time Quality Control

Abstract

The emergence of multiple global navigation satellite systems (multi-GNSS), including global positioning system (GPS), global navigation satellite system (GLONASS), Beidou navigation satellite system (BDS), and Galileo, brings not only great opportunities for real-time precise point positioning (PPP), but also challenges in quality control because of inevitable data anomalies. This research aims at achieving the real-time quality control of the multi-GNSS combined PPP using additional observations with opposite weight. A robust multiple-system combined PPP estimation is developed to simultaneously process observations from all the four GNSS systems as well as single, dual, or triple systems. The experiment indicates that the proposed quality control can effectively eliminate the influence of outliers on the single GPS and the multiple-system combined PPP. The analysis on the positioning accuracy and the convergence time of the proposed robust PPP is conducted based on one week’s data from 32 globally distributed stations. The positioning root mean square (RMS) error of the quad-system combined PPP is 1.2 cm, 1.0 cm, and 3.0 cm in the east, north, and upward components, respectively, with the improvements of 62.5%, 63.0%, and 55.2% compared to those of single GPS. The average convergence time of the quad-system combined PPP in the horizontal and vertical components is 12.8 min and 12.2 min, respectively, while it is 26.5 min and 23.7 min when only using single-GPS PPP. The positioning performance of the GPS, GLONASS, and BDS (GRC) combination and the GPS, GLONASS, and Galileo (GRE) combination is comparable to the GPS, GLONASS, BDS and Galileo (GRCE) combination and it is better than that of the GPS, BDS, and Galileo (GCE) combination. Compared to GPS, the improvements of the positioning accuracy of the GPS and GLONASS (GR) combination, the GPS and Galileo (GE) combination, the GPS and BDS (GC) combination in the east component are 53.1%, 43.8%, and 40.6%, respectively, while they are 55.6%, 48.1%, and 40.7% in the north component, and 47.8%, 40.3%, and 34.3% in the upward component.

1. Introduction

Real-time precise point position (PPP) is an innovative global positioning technique using a single receiver on the basis of real-time satellite orbit and clock products [1,2,3,4]. With the assistance of the uncalibrated phase delay (UPD) products and some atmospheric products such as tropospheric and ionospheric delay corrections, the PPP real-time kinematic (RTK) method achieves position accuracy at the centimeter level with an initialization time of 1 min or better [5,6,7]. Compared to the traditional reference station-based RTK positioning approach, the real-time PPP technique overcomes the limitation of the distances between reference stations and could be performed in any place where satellite correction data are available [8,9,10,11,12]. In addition, the real-time PPP is experiencing rapid development towards multiple systems and multiple frequencies for further realizing high-precision and high-reliability positioning.

Real-time global positioning system (GPS) PPP has been widely used in the real-time determination of integrated water vapor, precision agriculture applications, and natural hazards monitoring [13,14,15,16,17]. However, it takes much longer convergence time with only GPS in some cases, depending on satellite geometry and prevailing atmospheric conditions [18,19]. With the newly emerging Galileo and Beidou navigation satellite system (BDS) as well as the ongoing modernization of GPS and global navigation satellite system (GLONASS), the prosperity of multiple global navigation satellite systems (multi-GNSS) brings great opportunities for PPP to accelerate initialization time and improve positioning performance [20]. Having more available satellites increases the number of redundant observations and improves satellite geometry remarkably [21]. Unfortunately, the quality control problem also arises inevitably as the number of observations increases. As a result, the parameter estimation of the multi-GNSS combined PPP will be influenced and it is difficult to take full advantage of the combined systems.

Statistical testing procedures and robust estimation methods have been widely investigated in the adjustment of GNSS data in order to control the influences of outliers. Different strategies are proposed for quality control in different scenarios, such as the data snooping of a single outlier [22], the statistical testing procedure for multiple outliers [23,24], and the quasi-accurate detection method and robust estimation scheme for correlated observations [25,26]. The robust estimation method not only constructs the equivalent weights in the observation model, but also adjusts the updated parameters in the dynamic model to resist the influences of outlying kinematic model errors [27]. For statistical testing methods, the assumption of multiple alternative hypotheses is established using the test statistic to detect and identify outliers according to the critical region and the probability density function [28]. These methods are usually based on rigorous and perfect theory. However, the computation is complex and even practically impossible, although possible in theory.

For these quality control methods, the outliers are first detected and identified according to the strict analysis of measurement residuals. They are then removed or downweighted until the normal distributions of a posteriori residuals are satisfied. The readjustment is carried out by adding model parameters or updating the weights of observations. To simplify the quality control procedure of the multi-GNSS combined PPP, an approach is proposed that uses the experienced thresholds of the maximum residual and the maximum unit weight standard deviation (STD) as the criteria to judge whether large abnormal measurements exist. The influences of outliers can be removed from the normal equation by the traditional method, which gives the outliers zero weights and then makes readjustments. In order to avoid the recomputation and readjustment of the normal equation, an equivalent method is used by assuming the additional observations have opposite weight and offsetting their contribution to the normal equation against the influence of the corresponding outlier [29].

In this paper, Section 2 firstly discusses the observation model and the estimation equation of multi-GNSS combined PPP, which is needed by the robust model of the additional observations with opposite weight and real-time PPP implementation. Section 3 describes the data processing strategy of the multi-GNSS combined PPP. The experimental validation is then presented and the proposed method evaluated in Section 4. Finally, we come to some conclusions in Section 5.

2. Methodology

For a better understanding of the combined PPP robust estimation, the observation model and the equation of multi-GNSS PPP are first discussed and then the sequential least squares solution is presented for completeness, although it is common knowledge. On this basis, the quality control is studied based on the robust model using additional observations with opposite weight. Finally, its implementation in the PPP estimation is introduced in detail.

2.1. Observation Model of Multi-GNSS PPP

For a satellite s observed by a receiver r, the basic equations for raw GNSS pseudorange and carrier phase observations are expressed as follows:

$$\begin{array}{l}{p}_{r,j}^s={\mu }_r^s\Delta x_r+t_r-t^s+d_{r,j}^s-d_j^s+{\beta }_{r,j}^s{I}_{z,r}+m_r^s{Z}_r+{\epsilon }_{r,j}^s\\ l_{r,j}^s={\mu }_r^s\Delta x_r+t_r-t^s-{\beta }_{r,j}^s{I}_{z,r}+m_r^s{Z}_r+{\lambda }_j^s\left(B_{r,j}^s+b_{r,j}^s-b_j^s\right)+{\xi }_{r,j}^s\end{array}$$

(1)

where $p_{r,j}^s$ and $l_{r,j}^s$ are the observed minus computed (OMC) pseudorange and phase observations on frequency band j in meters; ${\mu }_r^s$ is the unit row vector from the receiver to the satellite after linearization; $\Delta x_r$ is the column vector of the receiver’s position increments relative to the a priori position; $t_r$ and $t^s$ are the clock offsets for the receiver and satellite, respectively; $d_{r,j}^s$ and $d_j^s$ are the frequency-dependent uncalibrated code delays (UCDs) for the receiver and satellite; $I_{z,r}$ is the zenith total ionospheric content and ${\beta }_{r,j}^s={\gamma }_r^s\cdot 40.3/f_j^2$ is a frequency-dependent scalar factor, in which ${\gamma }_r^s$ is the ionospheric mapping function and $f_j$ is the frequency [30]; $Z_r$ is the zenith tropospheric delay and $m_r^s$ is the corresponding mapping function; $B_{r,j}^s$ is the integer phase ambiguity in unit of cycles and ${\lambda }_j^s$ is the corresponding wavelength; $b_{r,j}^s$ and $b_j^s$ are the frequency-dependent UPDs for the receiver and satellite; and ${\epsilon }_{r,j}^s$ and ${\xi }_{r,j}^s$ are the measurement noise of pseudorange and carrier phase observations.

The ionospheric-free pseudorange and carrier phase observables are usually utilized to remove the first-order ionospheric delay errors [31]. The satellite-dependent UCDs can be also eliminated from the ionospheric-free observation model as it is involved in the precise satellite clock products, which are generated with the ionospheric-free undifferenced observations [32]. Because of different signal frequencies and signal structures, the intersystem biases (ISBs) and interfrequency biases (IFBs) should be considered for multi-GNSS combined PPP [33]. However, the introduction of the GLONASS IFB parameter will increase the complexity of the solution. Considering that it does not influence the validation of this proposed robust approach, the GLONASS IFB parameter is approximately treated as one ISB parameter. If GPS system is chosen as the reference system, the ionospheric-free model of the quad-system combination may be expressed as

$$\begin{array}{l}{p}_{r,IF}^{s,G}={\mu }_r^{s,G}\Delta x_r+{\overline{t}}_r+m_r^{s,G}{Z}_r+{\epsilon }_{r,IF}^{s,G}\\ l_{r,IF}^{s,G}={\mu }_r^{s,G}\Delta x_r+{\overline{t}}_r+{\overline{B}}_{r,IF}^{s,G}+m_r^{s,G}{Z}_r+{\xi }_{r,IF}^{s,G}\end{array}$$

(2)

$$\begin{array}{l}{p}_{r,IF}^{s,SYS}={\mu }_r^{s,SYS}\Delta x_r+{\overline{t}}_r+IS{B}_r^{SYS-G}+m_r^{s,SYS}{Z}_r+{\epsilon }_{r,IF}^{s,SYS}\\ l_{r,IF}^{s,SYS}={\mu }_r^{s,SYS}\Delta x_r+{\overline{t}}_r+IS{B}_r^{SYS-G}+{\overline{B}}_{r,IF}^{s,SYS}+m_r^{s,SYS}{Z}_r+{\xi }_{r,IF}^{s,SYS}\end{array}$$

(3)

where “SYS” represents the GLONASS (R), BDS (C), and Galileo (C) systems, respectively, and the symbol “G” refers to the GPS system; the receiver clock can be expressed as

$${\overline{t}}_r=t_r^G+d_{r,IF}^G$$

(4)

and the biases are defined as

$$IS{B}_r^{SYS-G}=d_{r,IF}^{SYS}-d_{r,IF}^G$$

(5)

where $d_{r,IF}^G$ and $d_{r,IF}^{SYS}$ denote the ionospheric-free receiver UCDs of the GPS reference system and other systems. Considering the satellite clock products contain the satellite-dependent UCDs, the ambiguity can be expressed as

$$\begin{array}{l}{\overline{B}}_{r,IF}^{s,G}=B_{r,IF}^{s,G}+b_{r,IF}^G-b_{IF}^{s,G}+d_{IF}^{s,G}-d_{r,IF}^G\\ {\overline{B}}_{r,IF}^{s,SYS}=B_{r,IF}^{s,SYS}+b_{r,IF}^{SYS}-b_{IF}^{s,SYS}+d_{IF}^{s,SYS}-d_{r,IF}^{SYS}\end{array}$$

(6)

where $B_{r,IF}^{s,G}$ and $B_{r,IF}^{s,SYS}$ denote the ionospheric-free float ambiguity of the GPS reference system and other systems; $b_{r,IF}^G$ and $b_{r,IF}^{SYS}$ denote the ionospheric-free receiver UPDs; $b_{IF}^{s,G}$ and $b_{IF}^{s,SYS}$ denote the ionospheric-free satellite UPDs; and $d_{IF}^{s,G}$ and $d_{IF}^{s,SYS}$ denote the ionospheric-free satellite UCDs.

2.2. Multi-GNSS Combined PPP Estimation

Assuming that the receiver tracks n satellites, where the number of GPS, GLONASS, BDS, and Galileo satellites is $n_G$, $n_R$, $n_C$, and $n_E$, respectively, and the matrix $P$ is the weight matrix of the observations, the observation equation for the quad-system combined PPP can be written as

$$L=AX+e,P$$

(7)

where $L$ is the OMC vector, $A$ is the design matrix, $X$ is the estimated parameter vector, and $e$ is the measurement noise vector; for the elements of matrices, see the reference [30].

Using the sequential least squares adjustment, the solution of multi-GNSS PPP estimation is as follows:

$$\begin{array}{l}\widehat{X}={\left(A^{T}PA+\overline{N}\right)}^{-1}\left(A^{T}PL+\overline{W}\right)=N^{-1}W\\ V=A\widehat{X}-L\\ {\widehat{\sigma }}_{STD}=\sqrt{\frac{V^{T}PV}{n}}\end{array}$$

(8)

where $\overline{N}$ and $\overline{W}$ are the existing normal equation and will be replaced by $N$ and $W$, which are the normal equation after contributing the observations at this epoch. $V$ is the residual vector of dimension n, $\widehat{X}$ is the parameter vector, and ${\widehat{\sigma }}_{STD}$ is the unit weight STD.

2.3. Robust Model Based on Additional Observations with Opposite Weight

Let the residual of the $i^{th}$ observation be the largest, and if it is defined as one outlier, its contribution can be removed by combining Equation (7) and an additional observation equation, which is same as the $i^{th}$ observation equation but with the opposite weight. The additional observation can be expressed as

$$l_i=a_{i}X+e_i,-p_i$$

(9)

where $a_i$ is the $i^{th}$ row vector of the design matrix $A$, and $l_i$ is the $i^{th}$ element of the OMC vector $L$. $p_i$ is the $i^{th}$ diagonal element of the weight matrix $P$.

By combining Equation (7) and Equation (9), the robust model of real-time PPP estimation becomes the following formulas [34]:

$$\begin{array}{l}{N}_{new}^{-1}=N^{-1}-\left(N^{-1}{a}_i^T\right)\beta \left(a_i{N}^{-1}\right)\\ X_{new}=X+\left(N^{-1}{a}_i^T\right)\beta \left(l_i-a_{i}X\right)\\ \beta ={\left(1/p_i+a_i{N}^{-1}{a}_i^T\right)}^{-1}\end{array}$$

(10)

Substituting the new estimates $X_{new}$ into Equation (8), the updated residuals can be expressed as $V_{new}$. The unit weight STD can be obtained and replaced by

$${\widehat{\sigma }}_{STD}=\sqrt{\frac{V_{new}^{T}P{V}_{new}}{n-1}}$$

(11)

where the number of observations reduces by one because of the removal of the outlier observation.

2.4. Real-Time Robust PPP Implementation

The real-time detection of cycle slips usually takes the TurboEdit algorithm, which uses the differences between epochs of both the geometry-free and Melbourne–Wübbena (MW) observations to identify cycle slips [35]. Because the differences are influenced by the observation interval and the elevation angle, the empirical thresholds are generally determined by some tests in advance [36]. Since it is difficult to detect all possible cycle slips and some outliers occur inevitably, especially for observations at low elevation angle, the PPP estimation would be abnormal without the robust control. As a result, the unit weight STD and the maximum residual are larger than the theoretical observation noise. Therefore, the unit weight STD can be used as the overall quality control and the maximum residual is imposed simultaneously to identify all possible outliers. The criteria of the quality control can be expressed as

$$\begin{array}{l}{\widehat{\sigma }}_{STD}<{\sigma }_0\\ \max \left|V\right|<c{\widehat{\sigma }}_{STD}\end{array}$$

(12)

where ${\sigma }_0$ is the experience value of the unit weight STD; and $c$ is usually chosen as 3.0-6.0 [37]. In our experiments, ${\sigma }_0$ is set as 0.020 m and $c$ uses 3.0 based on the observation noise and the model error.

After the PPP estimation without quality control is finished at one epoch, the normal equation and a posteriori residuals are available. The observation corresponding to the maximum residual will be identified according to Equation (9) if the unit weight STD or the maximum residual does not satisfy the condition in Equation (12). According to Equation (10), the influence of this abnormal observation will be offset against the normal equation by the use of one additional observation with the opposite weight. The influence of the unmodeled error can be easily eliminated without the complicated model assumption. After one outlier is removed, the unit weight STD and the residual are recomputed with the updated normal equation. The above procedure is carried out through adding the identified observation with opposite weight one by one until the condition is satisfied. The PPP robust estimation is realized at this epoch and the procedure can be repeated for next epoch. This quality control method is implemented into our real-time PPP software package, which was designed for real-time applications with the C++ programming language. The flow chart of the quality control using the additional observation with opposite weight and the real-time robust PPP implementation is presented in Figure 1.

3. Data Processing Strategy

The GPS L1/L2, GLONASS G1/G2, Galileo E1/E5a, and BDS B1/B2 observations were used for the ionosphere-free combination. The estimation parameters of the multi-GNSS combined PPP included the station coordinate, receiver clock offset, tropospheric zenith wet delay, ISB, and ionospheric-free combination ambiguity. The differential code biases (DCBs) of GPS and GLONASS were corrected with products from the center for orbit determination in Europe (CODE) and was not corrected for BDS and Galileo, since there are not yet official values available at present, to our knowledge.

Table 1. The processing strategy and error model in the multiple global navigation satellite systems (multi-GNSS) combined PPP.

Items	Models
Observations	Undifferenced phase and code ionosphere-free combination observations
Elevation mask	7°
Observation weight	Elevation-dependent weight, 1 for elevation angle (E) > 30°, otherwise 2sin(E)
Adjustment	Sequential least square
Phase-windup effect	Corrected [38]
Satellite antenna phase center and variation	Global positioning system (GPS), global navigation satellite system (GLONASS), Galileo: IGS14.atx; Beidou navigation satellite system (BDS): European space agency (ESA) Model [39,40]
Receiver antenna phase center and variation	GPS: IGS14.atx; GLONASS, Galileo, BDS: same as GPS
Station displacement	Solid Earth tide, pole tide, ocean tide loading: IERS Convention 2010 [41]
Relativistic effects	Corrected (IERS 2003)
Satellite orbit	German research centre for geosciences (GFZ) multi-GNSS (GBM) final product
Satellite clock	GFZ GBM final product
Receiver clock	Estimated, white noise
Station coordinate	Estimated, white noise
Tropospheric delay	Saastamoinen model and global mapping function (GMF) [42], and estimated as piecewise constants with a step size of one hour
the intersystem biases (ISBs)	Estimated as constant parameters
Phase ambiguities	Estimated as a constant for each continuous arc

shows the processing strategy and error model in the multi-GNSS combined PPP.

To validate the proposed robust approach, one week’s data from 32 stations from 5–11 March (day of year (DOY) 064–070) 2017 were chosen from the multi-GNSS experiment (MGEX) network. These stations are globally distributed and equipped with quad-system (GPS, BDS, GLONASS, and Galileo) receivers. The distribution is shown in Figure 2. The simulated real-time kinematic PPP was carried out with a sampling interval of 30 s for a single GPS system and multi-GNSS combination using the German research centre for geosciences (GFZ) multi-GNSS (GBM) final product, which provides GPS, BDS, GLONASS, and Galileo orbit- and clock-based daily solutions.

4. Experimental Validations

In order to avoid the influence of day boundary error caused by the orbit and clock product, only 20 hours’ data from 02:00:00 to 22:00:00 were processed for each day. The positioning details at the station BRST on 5 March 2017 are first presented as an example, in which the influence of cycle slips on PPP estimation is depicted. The positioning errors and satellite numbers of the combined solutions are also presented. Afterwards, the positioning performance is evaluated for all stations based on the proposed robust PPP approach.

4.1. Influence of Cycle Slips on PPP Estimation

The positioning performance of the PPP estimation will be decreased by the unmodeled error due to the failure in detecting cycle slips. To analyze the influence of cycle slips on PPP estimation and validate the proposed quality control based on additional observations with opposite weight, the cycle detection algorithm uses the relatively loose threshold condition where the cycle slip will be defined if the epoch-difference geometry-free phase combination is more than three times larger than the difference between wavelengths on two frequencies, or the epoch-difference MW combination is larger than four cycles. By applying these experienced thresholds, the precision of PPP estimation is investigated in terms of the unit weight STD and the maximum weighted residual at each epoch. The weighted residuals are used to process the phase and pseudorange residuals uniformly. Figure 3 shows the unit weight STD with and without quality control, and Figure 4 shows the corresponding maximum weighted residuals.

It can be seen from Figure 3 and Figure 4 that there are several large anomalies on the unit weight STD curve and the maximum weighted residual curve of the single GPS PPP solution without quality control. For example, the maximum weighted residual occurs on the G07 phase data between the epochs 1142 and 1151 (11:31:00 and 11:35:00, coordinated universal time (UTC)) and reaches about 0.26 m, which is 0.84 m in terms of the original residual. The corresponding unit weight STD is about 0.18 m at this epoch. Similarly, the phase weighted residuals of the G05 satellite are also abnormal, with the variation from 0.11 m to 0.07 m at epochs 2180 to 2297 (20:09:30 to 21:08:00, UTC time) and the unit weight STD is from 0.06 m to 0.04 m. The corresponding original residual is from 0.32 m to 0.28 m. At the last few epochs after 21:46:00, the weighted residual of the G29 satellite is also abnormal. There are also some other small anomalies for the solution without quality control. The anomalies might be caused by the ignorance of cycle slips, where one is about 0.24 m and 0.19 m for L1 and L2, respectively. The magnitude is enlarged by the ionospheric-free combination and other observations’ residuals might be influenced through their common parameters. In order to analyze these anomalies, the times of the difference between wavelengths on two frequencies for the epoch-difference geometry-free observations and the elevation angle are presented in Figure 5.

Figure 5 shows that the epoch difference of the geometry-free phase combination is abnormal at epochs around 11:31:00, 20:09:30, and 21:50:30 for G07, G05, and G29 satellites, respectively. It indicates that cycle slips are easy to occur at the period when the elevation angle is low and they are difficult to be detected by the data preprocessing procedure. The larger residual and worse PPP estimation performance will be caused by these unmodeled cycle slips. The geometry-free phase variation of the G07 satellite is largest and reflects the largest cycle slip, which results in the greatest impact on the positioning error curve. Figure 6 shows the positioning error of the single GPS PPP solution with quality control and without quality control.

Without quality control, the positioning error fluctuates more widely and the convergence time is longer than that of the solution with quality control. The additional observations with opposite weight are used to remove the influence of the undetected cycle slips if the unit weight STD or the maximum residual does not satisfy the threshold condition. The improvements of 4.32 cm, 0.25 cm, and 3.09 cm on the positioning RMS are achieved in the east, north, and upward components, respectively. Although the positioning performance can be improved from the quality control, there are still small fluctuations on the positioning error curve since the single GPS observations with low elevation angle are used.

4.2. Comparisons of Real-Time Multi-GNSS PPP with and without Quality Control

To test the performance of the quality control based on additional observations with opposite weight in the multi-GNSS combined PPP, the dual-system, triple-system, and quad-system combined PPP experiments were first carried out at the station BRST. The solution without quality control is also presented as a comparison. On the basis of the GPS system, seven different combined schemes were used, defined as follows: GPS/BDS (GC), GPS/Galileo (GE), GPS/GLONASS (GR), GPS/BDS/Galileo (GCE), GPS/GLONASS/BDS (GRC), GPS/GLONASS/Galileo (GRE), and GPS/GLONASS/BDS/Galileo (GRCE).

Figure 7 shows the positioning errors of the seven schemes without and with quality control, respectively. There are more spikes on the positioning error curve without quality control than the solution with quality control due to the unremoved outliers. Although the combination significantly increases the number of observations, the quad-system combination rarely improves the positioning accuracy comparing to the dual-system or triple-system combination because of the influence of outliers on the PPP solutions without quality control. The positioning errors of the robust GC combined PPP, which are denoted by the black curve, achieve obvious improvements. For example, the biggest jump on the positioning curve without quality control is well fixed at epochs around the thirteenth hour. The results of the other combined schemes are also improved remarkably and better positioning accuracy can be achieved. It should be noted that a small positioning error also exists in the GC combined solutions. With the number of observations increasing, it becomes small in the GE and GCE combinations and completely disappears in the GR, GRC, GRE, and GRCE combinations.

Table 2. The RMS statistic of the positioning error in east, north, and upward components after convergence.

System	No Control			Quality Control
	E (cm)	N (cm)	U (cm)	E (cm)	N (cm)	U (cm)
GRCE	4.66	4.87	10.17	0.79	0.86	4.24
GRC	4.86	5.23	11.35	0.99	0.88	4.44
GRE	4.33	1.45	7.55	0.90	0.85	4.43
GCE	5.09	7.22	11.50	1.14	1.77	4.96
GR	4.37	1.51	7.37	0.98	0.82	4.35
GE	4.71	1.61	7.94	1.32	1.79	4.72
GC	6.68	8.78	12.80	1.81	1.82	5.64

provides the RMS statistic of the positioning error in east, north, and upward components after convergence, and

Table 3. The convergence time statistic in horizontal (H) and vertical (V) components.

System	No Control		Quality Control
	H (min)	V (min)	H (min)	V (min)
GRCE	5.5	45.5	8.0	6.5
GRC	31.0	49.0	9.0	9.5
GRE	22.5	4.5	8.5	4.5
GCE	48.5	49.0	13.0	31.5
GR	31.5	49.0	9.0	6.5
GE	45.5	49.5	14.5	31.5
GC	64.0	80.5	29.5	40.5

shows the convergence time, which is defined as the period from the first epoch to the converged epoch. The positioning is considered to have converged when the positioning errors reach 0.1 m and remain constant for the next 20 epochs. Obviously, the positioning accuracy and the convergence time without quality control are poorer than the solution with quality control for almost all schemes. Although the positioning performance cannot be improved for one single component, in a few cases where the positioning accuracy or the convergence time is already good, the quality control can still improve the results as a whole. Among the schemes with quality control, the positioning accuracy of the GC combined PPP is worst, since the number of the available satellites is the lowest and the BDS geostationary earth orbit (GEO) satellite orbit and clock products have the poorest precision because of the particular orbit design for covering the Asia–Pacific region. Unfortunately, it is also vulnerable to outliers if it lacks quality control. The proposed quality control is the most effective for this case and the best improvement can be achieved. As a result, the positioning RMS errors are 1.81 cm, 1.82 cm, and 5.64 cm in the east, north, and upward components, respectively, and the convergence time in the horizontal and vertical components is 29.5 min and 40.5 min, respectively. For the triple-system combination scheme, the GCE combined PPP is similar: the corresponding RMS is 1.14 cm, 1.77 cm, and 4.96 cm, and the convergence time is 13.0 min and 31.5 min.

Figure 8 describes the number of available satellites for each scheme and the corresponding position dilution of precision (PDOP). It can be seen that the number of available satellites is generally more than 20 for the quad-system combination and almost 30 at some epochs. The quad-system PDOP is the best and the average is 0.65. For the GRC and GRE combination, the PDOP value is about 0.70, which is slightly larger than the quad-system constellation since the latter has more visible satellites. Therefore, the quality control can remove all possible outliers one by one, benefiting from the best geometry. Table 3 shows that with the quality control, the GR, GRC, and GRE combinations achieve an accuracy of about 1 cm in the east and north components, while it is about 4 cm in the upward component, which is only about 2 mm worse than the GRCE combined PPP. Their convergence time is less than 10 min in horizontal and vertical components. The number of visible satellites for the GCE combination is almost equal to that of the GR combination. The average satellite numbers are 15 and 16, respectively, while the corresponding PDOP values are 0.81 and 0.76. Therefore, the GR combination can achieve comparable positioning performance with the triple-system combinations. However, the number of visible satellites is only about 12 and the PDOP value is about 0.90 for both GC and GE combination, which leads to relatively worse positioning accuracy and much longer convergence time.

4.3. Assessment of Positioning Performance with the Proposed Quality Control Method

In order to assess the positioning performance of the multi-GNSS combined PPP based on the proposed quality control method, the observation data collected at 32 stations from March 5 to 11 (DOY 064–070) 2017 were used to calculate the positioning results of the seven combination schemes. The performance was analyzed based on daily solutions, for which the average positioning accuracy and the convergence time statistic are depicted in Figure 9, Figure 10 and Figure 11.

The bottom panel in Figure 9 denotes the result of the GPS system and the top panel is the GRCE combined result. It shows that the positioning RMS of the GRCE combination is better than the single GPS result. The positioning RMS of the GPS system and the GRCE combination in the east component is 3.2 cm and 1.2 cm, respectively, and it is 2.7 cm and 1.0 cm in the north component, while the RMS is 6.7 cm and 3.0 cm in the upward component. The corresponding positioning accuracy is improved by 62.5%, 63.0%, and 55.2% in the east, north, and upward components, respectively. The single GPS PPP takes a much longer time to converge. For example, only 35.7% and 42.9% of stations converge in 20 min in the horizontal and vertical components, respectively. With the integration of the other three systems, the percentage is improved to 75.5% and 77.7% correspondingly and the average time is 12.8 min and 12.2 min. Due to the significantly increased number of satellites, the GRCE combined PPP can still provide the reliable positioning service, although there are some observations with low elevation angle.

Figure 10 shows the results of the GRC, GRE, and GCE PPP. The positioning RMS values of the GRC and GRE PPP are comparable to the GRCE combined PPP and they are better than the accuracy of the GCE combined PPP. The further integration with another system for the triple-system combination, however, does not make large improvement in the positioning accuracy and the convergence time. Compared to the GRC combined PPP, the GRE combined PPP achieves slightly better positioning accuracy. For the horizontal and vertical components, 71.0% and 77.2% of stations for the GRE combination converge in 20 min, respectively, while 69.7% and 76.3% are obtained for the GRC combination. The average convergence time is about 13 min for both the GRE and GRC PPP. The GCE combined PPP takes a much longer time to converge than the GRE and GRC PPP and the average time is about 17 min.

The positioning RMS and the distribution of the convergence time of the GR, GC, and GE PPP are shown in Figure 11. It can also be seen that the combination of two systems not only improves the positioning accuracy, but also shortens the convergence time compared to the single GPS PPP. The improvement of the positioning RMS in the east component for the GR, GE, and GC combination is 53.1%, 43.8%, and 40.6%, respectively, while for the north component, it is 55.6%, 48.1%, and 40.7%. It is 47.8%, 40.3%, and 34.3% in the vertical component. The GE and GC PPP accuracy is slightly poorer than the GR PPP because the number of satellites is lower. The best dual-system combination is the GR combination, where the positioning performance is close to that of the triple-system combination. The convergence time of the GE and GC combination is about 22 min, while it is about 15 min for the GR combination. For all cases, the north component shows higher positioning accuracy than the east component, which is related to the satellite constellation configuration.

5. Discussions and Conclusions

In this paper’s work, the multi-GNSS combined PPP model was first built for processing arbitrary combinations as well as the single system positioning. For the quality control of the multi-GNSS combined PPP, a robust method was proposed based on additional observations with opposite weight. The corresponding model and the implementation of quality control were discussed in detail and then the performance was analyzed. Experimental results show that the positioning performance with quality control is better than the solution without quality control benefitting from the control of outliers.

The positioning accuracy and the convergence time of the multi-GNSS combined PPP with the proposed quality control were evaluated using one week’s data from 32 globally distributed stations. The GPS positioning accuracy in the east, north, and upward components is 3.2 cm, 2.7 cm, and 6.7 cm, respectively, and the convergence time is 26.5 min and 23.7 min in the horizontal and vertical components, respectively. The corresponding positioning RMS of the quad-system combined PPP is 1.2 cm, 1.0 cm, and 3.0 cm with an improvement of 62.5%, 63.0%, and 55.2% compared to GPS, and the convergence time is improved by 12.8 min and 12.2 min, respectively.

For the triple-system combined PPP, the positioning performance of the GRC and GRE PPP is close to that of the quad-system combined PPP and better than that of the GCE combined PPP. For the dual-system combined PPP, the best one is the GR combination, where the positioning accuracy and the convergence time are close to that of the triple-system combination. The convergence time of the GC and GE combinations is about 22 min. The positioning accuracy of the former is 1.9 cm, 1.6 cm, and 4.4 cm in the east, north, and upward components, respectively, and it is 1.8 cm, 1.4 cm, and 4.0 cm for the latter.

Considering the weak geometry caused by few visible satellites and the complexity of multi-GNSS observation data, it is sometimes difficult to solve all abnormal situations due to the limitation of the single quality control model. The integration of multiple quality control methods might help to improve the stability of the solutions. Therefore, the next step of our work is to study the combination of the statistical testing method, the adaptively equivalent weights method, and this proposed method.