Modeling and Assessment of GPS/BDS Combined Precise Point Positioning

Precise Point Positioning (PPP) technique enables stand-alone receivers to obtain cm-level positioning accuracy. Observations from multi-GNSS systems can augment users with improved positioning accuracy, reliability and availability. In this paper, we present and evaluate the GPS/BDS combined PPP models, including the traditional model and a simplified model, where the inter-system bias (ISB) is treated in different way. To evaluate the performance of combined GPS/BDS PPP, kinematic and static PPP positions are compared to the IGS daily estimates, where 1 month GPS/BDS data of 11 IGS Multi-GNSS Experiment (MGEX) stations are used. The results indicate apparent improvement of GPS/BDS combined PPP solutions in both static and kinematic cases, where much smaller standard deviations are presented in the magnitude distribution of coordinates RMS statistics. Comparisons between the traditional and simplified combined PPP models show no difference in coordinate estimations, and the inter system biases between the GPS/BDS system are assimilated into receiver clock, ambiguities and pseudo-range residuals accordingly.


Introduction
To prepare for the next phase of generating products for all GNSS available, the IGS [1] initiated the IGS Multi-GNSS Experiment (MGEX) campaign. The MGEX campaign focuses on tracking the newly available GNSS signals including the Chinese BeiDou Navigation Satellite System (BDS). BDS had been providing regional service since December 2012 with five Geostationary Earth Orbit (GEO) satellites, five Inclined Geosynchronous Satellite Orbit (IGSO) satellites and four Medium Earth Orbit (MEO) satellites. Many studies are being carried out and much progress has been achieved in BDS-only and GPS/BDS-combined precise data analysis. Precise BDS orbits/clocks are currently available at many institutes and several IGS analysis centers [2][3][4][5].
Many publications demonstrate the daily static Precise Point Positioning (PPP, [6] using BDS observations, and precise orbits/clocks could reach a precision of a few cm [5,7,8]. With the additional tracking of BDS constellations and therefore, a significant increase of observed satellites, GPS/BDS combined PPP could improve solution availability and accuracy by improving tracking geometry with a reduction of the position dilution of precision (PDOP) [9], especially in environments like urban canyons and ravines where there are limits in sky view. Chen et al. [10] demonstrated that near real time GPS/BDS combined PPP solution show a higher degree of precision and better robustness, which

PPP Model for Single System
Taking a GPS system as an example, pseudo-range and phase observation functions of the ionosphere-free combination between a receiver and a satellite G can be written as: where P G , L G are the pseudo-range and carrier phase observation; ρ G is geometrical distance; c is light speed and c¨dt G is receiver clock correction; b G and B G are receivers pseudo-range and carrier phase hardware delay bias; N G is ambiguity, and ZPD G is the slant tropospheric delay, ς G and ε G are noise terms.
The pseudo-range hardware delay biases b G is assimilated into clock correction c¨dt G . The carrier phase hardware delay biases B G is satellite dependent and stable over time, and it will be grouped into ambiguity [27,28]. Thus, the ionosphere-free combination PPP observation can be rewritten as: where dt G and N G are reformed station clock and ambiguity with: For BDS phase and pseudo-range observations the same Equations (2) and (3) also apply.

Traditional GPS/BDS Combined PPP Model
The ionosphere-free pseudo-range and phase observations for the GPS/BDS combined PPP can be written as: In Equation (4), the superscript and subscript G, C refers to GPS and BDS satellites. The slant tropospheric delay term ZPD could be modelled as the hydrostatic and wet parts with following: ZPD " m f h¨Z HD`m f w¨Z WD In Equation (5), ZHD, ZWD are the hydrostatic and wet zenith path delays, m f h , m f w are the mapping functions. The hydrostatic part of ZPD is normally corrected using empirical models and the term ZTD w is estimated.
The traditional GPS/BDS combined PPP model requires the estimation of an additional inter-system bias parameter. Defining GPS as the reference system, the inter-system bias could be written as, Considering Equations (5) and (6), the traditional GPS/BDS combined PPP equations are written as,

Simplified GPS/BDS Combined PPP Model
Based on the scaled sensitivity matrix (SSM, [19,29]) method, Chen et al. [19] proves that there is no correlation between ISB and coordinate parameters in multi-GNSS combined PPP, and ISB parameter can be removed in multi-GNSS PPP. Based on this conclusion, we remove the ISB term in Equation (7) and rewrite observation equations as: where dt 1 is defined as new receiver clock, N are new ambiguity terms. Compared to Equation (7), the ISB term is absorbed in the new terms and pseudo-range residuals with: where λ is called the element of the SSM, which quantitatively defines the ratio of ISB parameter assimilated into the receiver clock. As described in Chen et al. [19], the benefits of the new model are: (1) under some extreme circumstances, the traditional Multi-GNSS PPP model may final when there are only four satellites observed, while the new approach will still give coordinate solution; (2) the correlation between the clock and ambiguity parameters is reduced, thus making solution more stable using the new model; (3) Multi-GNSS PPP realization is simplified and unified. In this new model, observations of different GNSS systems are treated in a unified way as they were of the same satellite system. In this simplified model, the ISB is assimilated into a receiver clock and ambiguity parameter. Due to the much smaller weight assigned on pseudo-range observations, all pseudo-range observations will show big residuals. The GPS pseudo-range residuals are close to the amount of ISB value that assimilates into station clock parameter. The BDS pseudo-range residuals are close to the amount of ISB value that assimilates into BDS ambiguity.

Data Processing
To assess GPS/BDS PPP performance under different processing models, and to validate the simplified GPS/BDS combined PPP model, daily static and kinematic PPP data analysis is performed. Data of 11 IGS MGEX stations in January 2014 is used. Figure 1 shows the geographic distribution of these tracking stations.

Position Differences between Static PPP and IGS Daily Solutions
We compared GPS/BDS daily position estimates with the IGS daily solutions covering the same period through a seven-parameter Helmert transformation [28,36]. The RMS of each station is showed in Figure 2 and the mean value is shown in Table 1. GPS/BDS combined PPP achieves the best performance compared with GPS-only PPP and BDS-only PPP. The improvement of combined PPP over GPS-only PPP is not very obvious. BDS-only PPP has the worst accuracy, which may be due to the less accurate BDS orbit and clock products and a regional satellite constellation with less satellites compared with GPS. We notice that stations JFNG (Wuhan, China), CUT0 (Perth, Australia), GMSD (Guts Masda, Japan), and REUN (Le Tampon, La Reunion, France) have the best accuracy in BDS PPP, which is due to the fact that more BDS satellites are tracked for these sites under the current BDS constellation. Also, the GPS/BDS combined PPP using a simplified and traditional model is nearly For results evaluation and comparison, data processing is performed in the following four scenarios: GPS-only PPP, BDS-only PPP, and GPS/BDS combined PPP with traditional and simplified models, where station-wise ISB is estimated as a daily constant in a traditional model.
Precise satellite orbits and satellite clocks from SHA [17,30] are used for data analysis. We apply the IGS absolute antenna phase center model for GPS observations [31], the phase-wind up modeling [32] and the station displacement are modeled according to the IERS conventions 2010 [33]. A cut-off angle of 7˝is set for usable measurements and data sampling is set to 30 s. An elevation-dependent weighting strategy is applied to measurements at low elevations and the noise ratio between pseudo-range and phase observation is 500:1. Moreover, we estimate wet zenith path delay every hour by applying the most recently developed GPT2 empirical slant delay model [34]. An improved version of the LTW_BS software is used [35].

Position Differences between Static PPP and IGS Daily Solutions
We compared GPS/BDS daily position estimates with the IGS daily solutions covering the same period through a seven-parameter Helmert transformation [28,36]. The RMS of each station is showed in Figure 2 and the mean value is shown in Table 1. GPS/BDS combined PPP achieves the best performance compared with GPS-only PPP and BDS-only PPP. The improvement of combined PPP over GPS-only PPP is not very obvious. BDS-only PPP has the worst accuracy, which may be due to the less accurate BDS orbit and clock products and a regional satellite constellation with less satellites compared with GPS. We notice that stations JFNG (Wuhan, China), CUT0 (Perth, Australia), GMSD (Guts Masda, Japan), and REUN (Le Tampon, La Reunion, France) have the best accuracy in BDS PPP, which is due to the fact that more BDS satellites are tracked for these sites under the current BDS constellation. Also, the GPS/BDS combined PPP using a simplified and traditional model is nearly the same and their values differ at the level of µm, which will be discussed in detail in Sections 3.3 and 4.

Position Differences between Static PPP and IGS Daily Solutions
We compared GPS/BDS daily position estimates with the IGS daily solutions covering the same period through a seven-parameter Helmert transformation [28,36]. The RMS of each station is showed in Figure 2 and the mean value is shown in Table 1. GPS/BDS combined PPP achieves the best performance compared with GPS-only PPP and BDS-only PPP. The improvement of combined PPP over GPS-only PPP is not very obvious. BDS-only PPP has the worst accuracy, which may be due to the less accurate BDS orbit and clock products and a regional satellite constellation with less satellites compared with GPS. We notice that stations JFNG (Wuhan, China), CUT0 (Perth, Australia), GMSD (Guts Masda, Japan), and REUN (Le Tampon, La Reunion, France) have the best accuracy in BDS PPP, which is due to the fact that more BDS satellites are tracked for these sites under the current BDS constellation. Also, the GPS/BDS combined PPP using a simplified and traditional model is nearly the same and their values differ at the level of µm, which will be discussed in detail in Sections 3.3 and 4.   We make 3D RMS statistics of the 11 stations over the whole month and plot the histogram of the three scenarios in Figure 3. The GPS/BDS combined PPP achieves the best accuracy, with a mean 3D We make 3D RMS statistics of the 11 stations over the whole month and plot the histogram of the three scenarios in Figure 3. The GPS/BDS combined PPP achieves the best accuracy, with a mean 3D RMS of 1.5 cm, while those of GPS-only and BDS-only is 1.6 cm and 3.7 cm, respectively. The percentage of 3D RMS within the range of [0, 1.5] cm is 26%, 40%, and 59% for BDS-only, GPS-only and combined solutions, which demonstrates a more accurate solution is obtained by adding observations from different systems.

Position Differences between Kinematic PPP and IGS Daily Solutions
To assess kinematic position estimates differences among PPP scenarios, epoch-wise kinematic coordinate estimates are compared with the IGS daily estimates. In each daily kinematic solution, the epoch-wise coordinates after the first 1 h (120 epochs) are used for statistics analysis.

Position Differences between Kinematic PPP and IGS Daily Solutions
To assess kinematic position estimates differences among PPP scenarios, epoch-wise kinematic coordinate estimates are compared with the IGS daily estimates. In each daily kinematic solution, the epoch-wise coordinates after the first 1 h (120 epochs) are used for statistics analysis.

Position Differences between Traditional and New Model
Assessing the differences between the position estimates can directly illustrate to what extent the traditional and simplified models agree in their positioning results. In the following we compare the GPS/BDS combined static PPP results between the traditional and simplified models.
For each of the 11 IGS MGEX stations, we computed GPS/BDS combined static PPP position differences between the two models over 1 month. The biggest position difference is less than 1 mm and most points have position differences of less than 0.1 mm. Figure 6 shows the magnitude distribution of these coordinate differences in the North, East, and Up component.

Position Differences between Traditional and New Model
Assessing the differences between the position estimates can directly illustrate to what extent the traditional and simplified models agree in their positioning results. In the following we compare the GPS/BDS combined static PPP results between the traditional and simplified models.
For each of the 11 IGS MGEX stations, we computed GPS/BDS combined static PPP position differences between the two models over 1 month. The biggest position difference is less than 1 mm and most points have position differences of less than 0.1 mm. Figure 6 shows the magnitude distribution of these coordinate differences in the North, East, and Up component.

Position Differences between Traditional and New Model
Assessing the differences between the position estimates can directly illustrate to what extent the traditional and simplified models agree in their positioning results. In the following we compare the GPS/BDS combined static PPP results between the traditional and simplified models.
For each of the 11 IGS MGEX stations, we computed GPS/BDS combined static PPP position differences between the two models over 1 month. The biggest position difference is less than 1 mm and most points have position differences of less than 0.1 mm. Figure 6 shows the magnitude distribution of these coordinate differences in the North, East, and Up component. The mean biases are´0.6,´2.5 and 0.4 µm and RMS are 1.1, 3.9, and 1.9 µm for each coordinate component. In addition, about 91.5% in the North, 96.9% in the East, and 98.7% in the Up component of all deviations are within twice the standard deviations. These results verify the same position estimates of these two models.

Parameter Differences between Combined PPP Models
In the following, we compare the station clocks, ambiguities and pseudo-range residuals between the traditional and simplified GPS/BDS combined PPP models. Data used is the kinematic solution of station JFNG on DOY 028, 2014.

Kinematic Coordinate Difference
As discussed above, in the simplified GPS/BDS PPP model where ISB is not considered, the coordinate solution is indentical with that of traditional model. Here the coordinate difference between traditional and new model of JFNG's kinematic solution is presented in Figure 7. After 8 epochs the cooridnate differences is less than 1 cm, and after 30 epochs the difference reduces to less than 1 mm. After 500 epochs the difference is less than 5 µm and remains stable until the last epoch.

Station Clock and Ambiguity
In the simplified GPS/BDS PPP model, station clock, and ambiguity estimates absorb the ISB parameters together, and the percentage which goes into each parameter depends on the normal

Parameter Differences between Combined PPP Models
In the following, we compare the station clocks, ambiguities and pseudo-range residuals between the traditional and simplified GPS/BDS combined PPP models. Data used is the kinematic solution of station JFNG on DOY 028, 2014.

Kinematic Coordinate Difference
As discussed above, in the simplified GPS/BDS PPP model where ISB is not considered, the coordinate solution is indentical with that of traditional model. Here the coordinate difference between traditional and new model of JFNG's kinematic solution is presented in Figure 7. After 8 epochs the cooridnate differences is less than 1 cm, and after 30 epochs the difference reduces to less than 1 mm. After 500 epochs the difference is less than 5 µm and remains stable until the last epoch.

Parameter Differences between Combined PPP Models
In the following, we compare the station clocks, ambiguities and pseudo-range residuals between the traditional and simplified GPS/BDS combined PPP models. Data used is the kinematic solution of station JFNG on DOY 028, 2014.

Kinematic Coordinate Difference
As discussed above, in the simplified GPS/BDS PPP model where ISB is not considered, the coordinate solution is indentical with that of traditional model. Here the coordinate difference between traditional and new model of JFNG's kinematic solution is presented in Figure 7. After 8 epochs the cooridnate differences is less than 1 cm, and after 30 epochs the difference reduces to less than 1 mm. After 500 epochs the difference is less than 5 µm and remains stable until the last epoch.

Station Clock and Ambiguity
In the simplified GPS/BDS PPP model, station clock, and ambiguity estimates absorb the ISB parameters together, and the percentage which goes into each parameter depends on the normal

Station Clock and Ambiguity
In the simplified GPS/BDS PPP model, station clock, and ambiguity estimates absorb the ISB parameters together, and the percentage which goes into each parameter depends on the normal equation. In Equation (9), the ISB assimilation shows the opposite sign for station clock and GPS Sensors 2016, 16, 1151 9 of 13 ambiguity parameters, thus the sum quantity of station clock and GPS ambiguity for the traditional and simplified model is theoretically the same at each epoch. The upper subplot of Figure 8 shows the magnitude distribution of this sum quantity differences between the two models using all GPS satellite/station pairs at all epochs. The mean difference is 2.0 µm and all the differences are below 1 mm with around 0.35% bigger than 0.1 mm.
According to Equation (9), the epoch-wise sum of station clock and BDS ambiguity in the simplified model is theoretically the same as the sum quantity of station clock, ISB, and BDS ambiguity in the traditional model. The bottom subplot of Figure 8 shows the magnitude distribution of this sum quantity differences between the two models using all BDS satellite/station pairs at all epochs. The mean difference is 4.6 µm and all the differences are below 1 mm at around 0.72% bigger than 0.1 mm. Both plots show that the differences are so small, which further proves the same coordinate estimates in these two models.
Sensors 2016, 16, 1151 9 of 13 equation. In Equation (9), the ISB assimilation shows the opposite sign for station clock and GPS ambiguity parameters, thus the sum quantity of station clock and GPS ambiguity for the traditional and simplified model is theoretically the same at each epoch. The upper subplot of Figure 8 shows the magnitude distribution of this sum quantity differences between the two models using all GPS satellite/station pairs at all epochs. The mean difference is 2.0 µm and all the differences are below 1 mm with around 0.35% bigger than 0.1 mm. According to Equation (9), the epoch-wise sum of station clock and BDS ambiguity in the simplified model is theoretically the same as the sum quantity of station clock, ISB, and BDS ambiguity in the traditional model. The bottom subplot of Figure 8 shows the magnitude distribution of this sum quantity differences between the two models using all BDS satellite/station pairs at all epochs. The mean difference is 4.6 µm and all the differences are below 1 mm at around 0.72% bigger than 0.1 mm. Both plots show that the differences are so small, which further proves the same coordinate estimates in these two models.

Station Clock and Pseudo-Range Residual
In Equation (9), the ISB assimilation shows opposite signals for station clock and GPS pseudorange observation residuals, thus the epoch-wise sum quantity of station clock and GPS pseudorange residual is theoretically the same for the traditional and simplified model. The upper subplot of Figure 9 shows the station clock differences between the traditional and simplified PPP models, where the differences are between the range of −31 and −28 m. The bottom subplot of Figure 9 shows the magnitude distribution of the epoch-wise sum quantity differences between the two models for all BDS satellite/station pairs at all epochs. The mean difference is −0.7 µm and all the differences are below 0.1 mm.

Station Clock and Pseudo-Range Residual
In Equation (9), the ISB assimilation shows opposite signals for station clock and GPS pseudo-range observation residuals, thus the epoch-wise sum quantity of station clock and GPS pseudo-range residual is theoretically the same for the traditional and simplified model. The upper subplot of Figure 9 shows the station clock differences between the traditional and simplified PPP models, where the differences are between the range of´31 and´28 m. The bottom subplot of Figure 9 shows the magnitude distribution of the epoch-wise sum quantity differences between the two models for all BDS satellite/station pairs at all epochs. The mean difference is´0.7 µm and all the differences are below 0.1 mm.

Ambiguity and Pseudo-Range Residual
In Equation (9), the ISB assimilation is the same for ambiguity and pseudo-range observation residual in both the GPS and BDS systems. Thus, the ambiguity differences and pseudo-range residual differences between the two models should theoretically be the same at the same epoch for the traditional and simplified model. In the following, we define the epoch-wise "double difference" term by making difference between the ambiguity difference and pseudo-range residual difference between the two models, which by definition should theoretically be zero. Figure 10. (a,c) GPS and BDS Pseudo-range residuals between traditional and simplified models; (b,d) magnitude distribution of the double differences for GPS and BDS observations. The left corner of right subplots (subplot b) and d)) shows the bias and the standard deviation (σ), whereas the right corner of right subplots (subplot b) and d)) shows the percentages of deviations that are within 2σ and 3σ.

Ambiguity and Pseudo-Range Residual
In Equation (9), the ISB assimilation is the same for ambiguity and pseudo-range observation residual in both the GPS and BDS systems. Thus, the ambiguity differences and pseudo-range residual differences between the two models should theoretically be the same at the same epoch for the traditional and simplified model. In the following, we define the epoch-wise "double difference" term by making difference between the ambiguity difference and pseudo-range residual difference between the two models, which by definition should theoretically be zero.

Ambiguity and Pseudo-Range Residual
In Equation (9), the ISB assimilation is the same for ambiguity and pseudo-range observation residual in both the GPS and BDS systems. Thus, the ambiguity differences and pseudo-range residual differences between the two models should theoretically be the same at the same epoch for the traditional and simplified model. In the following, we define the epoch-wise "double difference" term by making difference between the ambiguity difference and pseudo-range residual difference between the two models, which by definition should theoretically be zero. Figure 10. (a,c) GPS and BDS Pseudo-range residuals between traditional and simplified models; (b,d) magnitude distribution of the double differences for GPS and BDS observations. The left corner of right subplots (subplot b) and d)) shows the bias and the standard deviation (σ), whereas the right corner of right subplots (subplot b) and d)) shows the percentages of deviations that are within 2σ and 3σ. Figure 10. (a,c) GPS and BDS Pseudo-range residuals between traditional and simplified models; (b,d) magnitude distribution of the double differences for GPS and BDS observations. The left corner of right subplots (subplot b) and d)) shows the bias and the standard deviation (σ), whereas the right corner of right subplots (subplot b) and d)) shows the percentages of deviations that are within 2σ and 3σ. The left subplots of Figure 10 show the GPS and BDS pseudo-range residual differences between the two models, where we see that they present a similar shape but with different signs as presented in Equation (9). The right subplot of Figure 10 shows the magnitude distribution of the double difference quantity using all satellite/station pairs at all epochs. The mean double differences are´2.5 and 5.0 µm for GPS and BDS. All the double differences are below 0.1 mm for GPS, while all the BDS double differences are below 1 mm with around 0.5% bigger than 0.1 mm.

Conclusions and Suggestions
In this study, we discuss the different parameterization of GPS/BDS PPP, where combined PPP between the traditional and a simplified model are compared. GPS/BDS PPP were performed with different models and scenarios using 1 month of GPS/BDS data from 11 IGS MGEX stations. Comparing daily static coordinate estimates of BDS-only, GPS-only and GPS/BDS combined solutions with the IGS daily estimates, the mean 3D RMSes were 3.7, 1.6 and 1.5 cm, respectively. In the kinematic case, the 3D RMSes were 17.0, 8.2, 6.5 cm. Magnitude distribution of the 3D RMS shows higher percentage of more accurate solutions in GPS/BDS combined solutions, which is due to the inclusion of data from different systems.
Comparing GPS/BDS combined PPP coordinate estimates between the traditional and simplified models, mean biases of the differences are´0.6,´2.5, and 0.4 µm and RMSes are 1.1, 3.9, and 1.9 µm in the North, East, and Up components, respectively. The analysis of parameter and pseudo-range residual differences between the two models show the same quantity with deviation at the µm level for: