Precise positioning with current multi-constellation Global Navigation Satellite Systems: GPS, GLONASS, Galileo and BeiDou

The world of satellite navigation is undergoing dramatic changes with the rapid development of multi-constellation Global Navigation Satellite Systems (GNSSs). At the moment more than 70 satellites are already in view, and about 120 satellites will be available once all four systems (BeiDou + Galileo + GLONASS + GPS) are fully deployed in the next few years. This will bring great opportunities and challenges for both scientific and engineering applications. In this paper we develop a four-system positioning model to make full use of all available observations from different GNSSs. The significant improvement of satellite visibility, spatial geometry, dilution of precision, convergence, accuracy, continuity and reliability that a combining utilization of multi-GNSS brings to precise positioning are carefully analyzed and evaluated, especially in constrained environments.

collect and analyze data of GPS, GLONASS, BeiDou and Galileo and also serves as a framework for increasing the overall awareness of multi-GNSS within the scientific and engineering communities 4 . As a backbone of the MGEX project, over the past one to two years a new network of multi-GNSS monitoring stations has been deployed around the globe in parallel to the legacy IGS network for GPS and GLONASS. The MGEX network has grown to more than 90 stations now and it provides an excellent opportunity to demonstrate the potential benefits from multi-constellation GNSS (http://igs.org/ mgex/).
The fusion of multiple GNSSs will significantly increase the number of observed satellites, optimize the spatial geometry and improve continuity and reliability of positioning. However, in the past the data fusion of multi-GNSS only focused on GPS and GLONASS [8][9][10] . Thanks to the completion of the BeiDou regional system, GPS 1 BeiDou combined precise orbit determination 7,11-13 (POD) and precise point positioning [14][15] (PPP) have been investigated recently. Some initial results about GPS 1 Galileo combined POD [16][17] and real time kinematic 18 (RTK) positioning were also presented.
In this contribution, we develop a BeiDou 1 Galileo 1 GLONASS 1 GPS four-system model to fully exploit all the observations from four navigation satellite systems and demonstrate its significant contribution to precise positioning with current multi-constellations. The data processing model and strategy for multi-GNSS precise positioning are described in detail in the 'Methods' section. We perform a rigorous multi-GNSS analysis with careful consideration of inter-system and inter-frequency biases. All the observations from different GNSS are processed together in one common parameter estimation procedure with appropriate ionospheric constraints. In the 'Results' section, the satellite visibility, positional dilution of precision (PDOP), signal-to-noise ratio (SNR), multipath, and phase and code residuals of different constellations are carefully analyzed and compared to assess their individual performance. Furthermore, a comprehensive analysis, including satellite visibility, spatial geometry, dilution of precision, convergence, accuracy, continuity and reliability, is performed to evaluate the contribution of multi-GNSS fusion to precise positioning, especially in constrained environments (e.g., urban canyons, open pits).

Results
In order to assess the precise positioning performance with current multi-constellation GNSS, we analyzed 100 days' data of MGEX and BETN 19 (BeiDou Experimental Tracking Network) networks from September 1 to December 9 (day of year 244 to 343) in 2013. Figure 1 shows the distribution of MGEX and BETN networks, which includes stations all over the world.
Taking the four-system station GMSD (Japan, 30.55uN, 131.01uE) as an example, the satellite visibilities for Galileo, BeiDou, GLONASS and GPS on September 3, 2013 are shown in Figure 2. It is obvious that a joint utilization of various GNSSs brings a significant improvement of satellite visibility. It also can be seen that the visibility is very different for the various orbital types. The BeiDou GEO satellites have the longest tracking periods, e.g., C01 , C04 are tracked for the whole day. The BeiDou IGSO satellites C06-C10 are tracked relatively shorter than GEOs with data gaps of about 4.5 h for one day. The tracking periods of MEO satellites are the shortest and about 6 , 10 h for one day. The special constellation of BeiDou, including five GEOs and five IGSOs, guarantees sufficient visible satellites in the Asia-Pacific area. The number of visible satellites for each GNSS and the total satellite number of four systems at GMSD on that day are shown in Figure 3. At each epoch, there are about 7 , 13 BeiDou, 5 , 10 GLONASS, and 7 , 14 GPS satellites, which ensures their autonomous positioning. The Galileo satellites in view vary between 0 , 3 and therefore this GNSS currently cannot perform autonomous positioning. The blue line shows, that the total satellite number of all four systems can even reach up to about 23 , 35 at each epoch. Such a redundancy of the satellite visibility significantly increases the reliability and robustness of precise positioning. Figure 4a, b and c show the sky plots (azimuth vs. elevation) of GPS, GLONASS, Galileo 1 BeiDou at GMSD, respectively. Compared with Figure 4d, the improvement of the spatial geometry is clearly visible when observations from the four navigation systems are used simultaneously. The PDOP values in single-system and four-system modes on that day are also calculated and shown in Figure 5. The GPS PDOP values vary between 1.5 and 3.5, while the GLONASS PDOP values are about 1.2 , 3. The BeiDou PDOP values are larger than GPS and GLONASS, about 1.8 , 4.5. The blue line of Figure 5 shows that the PDOP values in four-system mode are significantly reduced to below 1.0 in most of the time and the corresponding time series is very stable. Figure 6 shows the signal-to-noise ratio (SNR) and multipath combination (MPC 12 , which mainly contains code noise and multipath) values at GMSD to assess and compare the observation quality for the different satellite systems. As typical examples, the SNRs and MPCs of the satellites G02, R04 and E20 on the first frequency along with their variation of the elevations are shown in Figure 6a, b and c, respectively. The SNRs and MPCs of BeiDou GEO C02, IGSO C06  and MEO C13 are shown in Figure 6d, e and f, respectively. It can be clearly seen that both the SNR and MPC values are strongly correlated with the elevation variations. The higher the elevation angles are, the larger the SNR and the smaller the MPC values. Generally speaking, the SNR are higher than 40 and the MPC values are smaller than 0.5 m when the elevation angles are higher than 15u. We also find that the different orbital types show different variation characteristics for SNR and MPC. The IGSO satellites can be observed much longer and the corresponding SNR, MPC and elevation values change much slower than for the MEO satellites. Especially, the elevation angles of each GEO satellite stay around one specific value, the corresponding SNR and MPC values are the steadiest and nearly don't exhibit systematic variations.      GEO and MEO. It can be explained that the proportion of lowelevation periods for IGSO is lower than that of MEO as shown in Figure 6e, and the GEO satellites have the lowest low-elevation proportion as their elevations are always around some relatively high values within their service regions as shown in Figure 6d.
The multi-GNSS data from MGEX and BETN networks are processed in both single-and multi-system modes following the positioning algorithms presented in the Section 'Methods'. All the estimated station coordinates are compared with the IGS SINEX or weekly solution. Figure 8 shows static PPP solutions of single-, dualand four-system modes at the four-system station CUT0 in Australia, which is covered by the BeiDou service of the Asia-Pacific area. The left sub-figures show the single-system PPP results of GPS-only, BeiDou-only, GLONASS-only and Galileo-only, respectively. For the GPS-only solution, the positioning accuracy can be better than 1 dm after a convergence time of about 30 minutes. About 2 hours of convergence is required to ensure an accuracy of better than 5 cm in all three components. The mm accuracy can be achieved after the long convergence time of several hours. The convergence of GLONASS-only PPP is longer compared to GPS-only PPP, about 3 hours to achieve an accuracy of a few centimeters. Meanwhile, the GLONASS positioning accuracy after sufficient convergence time is also slightly worse compared to the GPS solution. The BeiDou-only PPP presents good performance in the horizontal components, few cm accuracy can be achieved within one hour. However, the vertical component is much more unstable than GPS and GLONASS. A Galileo-only PPP solution cannot be derived at this station as not enough satellites can be observed.
The combined GPS/BeiDou,GPS/GLONASS,GPS/Galileo,GPS 1 BeiDou 1 GLONASS 1 Galileo PPP solutions are shown in the right sub-figures. Obviously, the multi-GNSS combination significantly improves the PPP performance, compared to the left sub-figures of the single-system solutions. It can be clearly seen that the combined GPS/BeiDou and GPS/GLONASS solutions significantly shorten the convergence time and improve the position series compared to single-system PPP. The Galileo satellites also contribute to the combined GPS/Galileo PPP solution to some extent, although they are not sufficient for autonomous positioning. Especially, the combined GPS 1 BeiDou 1 GLONASS 1 Galileo PPP present the fastest convergence, the most stable position series and highest accuracy for all three components. It only takes several minutes to achieve an accuracy of better than 10 cm, less than 30 minutes to be better than 5 cm, and a few hours to reach mm level accuracy. Figure 9 shows static PPP solutions at another four-system station LMMF (14.59uN, 260.99uW, Martinique, Caribbean Sea), which is out of the Asia-Pacific region of BeiDou service. At this location, both BeiDou and Galileo cannot provide continuous precise positioning as a stand-alone system. The cm level accuracy for BeiDou-only or Galileo-only PPP is only available after a static observation of even longer than 16 hours. Importantly though, the addition of BeiDou or Galileo to the GPS-only processing, obviously shorten the convergence time and improve the position series. The combined GPS 1 BeiDou 1 GLONASS 1 Galileo PPP also shows the fastest convergence, the most stable position series and highest accuracy in all of the three components. Figure 10 shows the statistical results of the static PPP solutions with different session lengths, of 15 and 30 min, 1, 2, 4, 6, and 12 h. The root mean square (RMS) values are calculated from all static PPP solutions over all the selected days and stations (MGEX and BETN networks from September 1 to December 9 in 2013). The singlesystem PPP results are shown in the left sub-figures and it can be seen that the positioning accuracy is evidently improved along with the increase of the observational length. With the same session length, the GPS-only PPP can achieve the highest accuracy, while the BeiDou-only PPP has the worst performance. In addition, the positioning accuracy in north component is generally better than east and vertical components. The multi-system PPP solutions are shown in right sub-figures. Obviously, the multi-GNSS combination significantly improves the PPP performance, compared to the left subfigures of single-system solutions. Usually, single-system PPP requires two hours or longer observation to achieve a positioning accuracy of few centimeters. For multi-GNSS PPP, the accuracy of few centimeters is available within 30 minutes for all three components. With the same session length, multi-PPP accuracy is significantly better than single-system PPP. The north component is still the most accurate component in multi-PPP solutions (may be caused by the satellite constellation configuration). The RMS values of static PPP solutions with different session lengths in single-and four-system modes are listed in Table 1. Figure 11 shows the statistical RMS values of static PPP solutions under different cut-off elevation angles, ranging from 10u to 40u. The GPS-only PPP results are shown in the left sub-figures, while the four-system PPP solutions are shown in the right sub-figures. It can be clearly seen that the accuracy of GPS-only PPP decreases significantly when the cut-off elevation angle increases, especially for short session lengths. At the 40u cut-off elevation, the positioning accuracy with 30 min observations is reduced to about 18, 17, and 24 cm in east, north and vertical components, respectively. Even with 2 hours observations, the vertical accuracy is still worse than 10 cm. When the session length is 30 minutes or longer, it can be seen from the right sub-figures that the accuracy of multi-PPP is not obviously affected by high cut-off elevation and few centimeters are achievable in all the three components even with 40u cut-off elevation. Only the vertical accuracy is slightly reduced with the increase of the cut-off elevation angle. With only 15 min observations, the vertical accuracy decreases gradually as the cut-off elevation angle increases, but the horizontal accuracy is only slightly decreased even with 40u cut-off elevation. When the cut-off elevation is increased to 30u or 40u, a static observation of about four hours is required to obtain a positioning accuracy of few centimeters for GPS-only PPP. However, multi-system PPP only needs 30 min static observations to achieve few centimeters accuracy. This is important, since such high cut-off www.nature.com/scientificreports elevation capability will significantly increase the GNSS applicability in constrained environments, such as, e.g., in urban canyons, open pits or when serious low-elevation multipath or ionospheric scintillations are present.
All the Multi-GNSS data are also processed in kinematic mode, i.e. the station coordinates are estimated epoch-by-epoch without any constraints between the epochs. The kinematic PPP results of singleand four-system modes at station CUT0 are compared in Figure 12 as a typical example. The GPS-only, GLONASS-only and BeiDou-only position series are shown by the blue, yellow and green lines, respectively. The red line shows the four-system PPP results. It is clearly indicated that the position series of multi-PPP is significantly more stable compared to the single-system solutions with much smaller and fewer fluctuations. Several spikes, appearing in single-system solutions, can be removed when multi-GNSS observations are used together.
We calculated the RMS values of kinematic PPP solutions in single-and four-system modes for all the selected days and stations. The statistical results in east, north, and up components are shown in Table 2. For the single-system solutions, GPS has the best accuracy of 1.4, 1.2, and 4.4 cm in the east, north and up components, respectively. The horizontal accuracy of BeiDou-only kinematic PPP is comparable with that of GPS, but the BeiDou vertical accuracy is worse than GPS. The GLONASS-only PPP has the worst performance in all the three components with RMS values of 2.9, 1.8, and 5.9 cm. Compared to the multi-PPP results, one can see that the  We also analyzed the kinematic PPP performance in single-and multi-system modes under different cut-off elevation angles, ranging from 10u to 40u. The position series of the station CUT0 are shown in Figure 13 as a typical example. The GPS-only PPP results are shown by the red lines, while the four-system PPP solutions are shown by the blue lines. We can find that the accuracy and reliability of GPS-only PPP decreases dramatically as the cut-off elevation angle increases.
When the cut-off elevation is increased to 30u or 40u, the GPS-only PPP cannot provide continuous precise positioning, the PPP results are very unreliable and precise position estimates are frequently not available. However, the positioning accuracy of multi-system PPP is nearly not decreased and few centimeter are still achievable in horizontal components even with 40u cut-off elevation. The multi-PPP vertical accuracy decreases gradually as the cut-off elevation angle increases, but much better than that of single-system solutions.
Observation residuals, which mainly contain the observation noises, multipath, and other errors that are not fully modeled, can   also be used as an important index to assess the observation quality and positioning accuracy. Figure 14 shows the phase and code residuals of several typical satellites from different satellite systems and orbital types. It can be seen that the phase residuals are generally within 62.0 cm and the code residuals are generally within 63.0 m, except some low-elevation periods. The large errors of the code observations will not influence the positioning results significantly as the weighting of code observations is much smaller (usually about 1510,000) than that of phase observations. Figure 15 shows the RMS values of the phase and code residuals for different satellite systems and orbital types. The RMS value of GPS code residuals is about 1.6 m. The RMS value of GLONASS code residuals is about 1.3 m and smaller compared to GPS. This conclusion is different from the findings in some other publications, e.g., Cai and Gao (2013) 10 . These authors found that the GLONASS residuals are much larger compared to GPS. It can be explained by the fact that the GLONASS inter-frequency biases are well considered in our multisystem positioning model. In addition, the inter-frequency bias parameters can also absorb the systematic part of code errors to some extent. The RMS value of Galileo code residuals is smallest and about 1.2 m, while the RMS value of BeiDou MEO's code residuals is largest and about 1.7 m. For the different orbital types of BeiDou satellites, the RMS of GEO residuals is smallest and about 1.3 m, while the RMS of the IGSO residuals is about 1.5 m, which is between that of GEO and MEO. It is similar to the situation of MPCs and the possible reason is that the proportion of low-elevation periods for IGSO is lower than that of MEO and the GEO satellites have the lowest lowelevation proportion as their stable elevations within service regions. For the BeiDou phase observations, the RMS of the GEO residuals is also smallest and about 0.9 cm, while the RMS of MEO residuals is also largest and about 1.4 cm. We can also find that the RMS value of BeiDou MEO residuals is even slightly smaller than GPS residuals, which is about 1.7 cm. The GLONASS residuals are slightly larger than GPS ones and the corresponding RMS is about 1.9 cm. The Galileo residuals have the largest RMS of 2.5 cm, which can be caused by some modeling parameters that are not accurate enough (e.g. phase center offset and variation values).

Discussion
In order to take full advantage of the current available GNSSs, we developed a BeiDou 1 Galileo 1 GLONASS 1 GPS four-system positioning model to fully exploit all currently available GNSS observations. With the multi-GNSS data from the MGEX and BETN ground tracking networks, we first fully analyzed and compared the visibility, PDOP, SNR and MPC of the different constellations with real data. Precise positioning is then performed and the results show that the addition of BeiDou, Galileo and GLONASS systems to the standard GPS-only processing significantly shorten the convergence time and improve the positioning accuracy. Meanwhile, the position series of multi-PPP are much more stable than GPS-only solutions, with much smaller and fewer fluctuations. Some spikes, appeared in single-system solutions, can be easily solved when multi-GNSS observations are used together. The accuracy and reliability of GPS-only PPP decreases dramatically when the cut-off elevation angle increases, especially in kinematic applications. However, the accuracy of multi-PPP is not obviously affected by high cut-off elevation and few centimeters are achievable even with 40u cut-off elevation. In kinematic multi-PPP, the vertical accuracy decreases gradually as the cut-off elevation angle increases, but it is much better than that of single-system solutions.
The fusion of multiple GNSS can significantly increase the number of observed satellites, optimize the spatial geometry and improve convergence, accuracy, continuity and reliability of precise positioning. Especially, the high cut-off elevation capability of multi-GNSS will significantly increase its applicability in constrained environments, such as e.g. in urban canyons, open pits or when serious low-elevation multipath or ionospheric scintillations are present. In the sequential studies, the three-frequency observation data should be fully exploited to further improve the multi-GNSS performance. In addition, the multi-GNSS not only enhances precise positioning applications, but also offers an increased number of signals for GNSS based remote sensing as troposphere/ionosphere sounding with ground based and satellite based (radio occultation) techniques with numerous applications in atmosphere science and operational weather and space weather forecast 24,25 .

Methods
The GNSS observation equations for carrier phase L and pseudorange P respectively, can be expressed as following: where indices s, r , and j refer to the satellite, receiver, and carrier frequency, respectively; t s and t r are the clock biases of satellite and receiver; N s r,j is the integer ambiguity; b r,j and b j s are the receiver-and satellite-dependent uncalibrated phase delay 20 (UPD); l j is the wavelength; d r,j and d j s are the code biases of the receiver and the satellite; I s r,j is the ionospheric delay of the signal path at frequency j, the ionospheric delays at different frequencies can be expressed as Eq (3); Z r is the tropospheric zenith wet delay at the station r, the slant tropospheric delay consists of the dry and wet components and both can be expressed by their individual zenith delay and mapping function, the tropospheric delay is usually corrected for its dry component with an a priori model, while the residual wet part of the tropospheric delay is estimated from the observations; m s r is the wet mapping function; e s r,j and e s r,j denote the sum of measurement noise and multipath error for the pseudorange and carrier phase observations. Furthermore, r g denotes the geometric distance between the phase centers of the satellite and receiver antennas at the signal transmitting and receiving time, respectively. This means, that the phase center offsets and variations    www.nature.com/scientificreports and station displacements by tidal loading must be considered. Phase wind-up and relativistic delays must also be corrected according to the existing models 21 , although they are not included in the equations. Usually, the ionosphere-free linear combination is used for PPP to eliminate the ionospheric delays 21 . In this contribution, we use the raw carrier phase and pseudorange observations of Eqs. (1), (2), and (3) and estimate the slant ionospheric delay as unknown parameters 22 . The linearized equations for (1) and (2) can be respectively expressed as following, where l s r,j and p s r,j denote ''observed minus computed'' phase and pseudorange observables from satellite s to receiver r at the frequency j; u s r is the unit vector of the direction from receiver to satellite; r r denotes the vector of the receiver position increments relative to a priori position which is used for linearization. For the multiconstellation case, the combined BeiDou 1 Galileo 1 GLONASS 1 GPS observation model can be formulated as, where indices G, R, E and C refer to the GPS, GLONASS, GALILEO, and BeiDou satellites, respectively; R k denotes the GLONASS satellite with frequency factor k that are used for the computation of the carrier phase frequencies of the individual GLONASS satellites; d rG , d rR , d rE , and d rC denote the code biases of the receiver r for G, R, E and C, respectively. Because of the different frequencies and signal structure of the individual GNSS, the code bias values d rG , d rR , d rE , and d rC are different in one multi-GNSS receiver. The differences between them are usually called inter-system biases (ISB) for code observations. Similarly, the phase delays b rG , b rR , b rE and b rC are also different and their differences are inter-system biases for phase observations; As GLONASS satellites emit the signals on individual frequencies, it will also lead to frequency-dependent biases in the receivers. For the GLONASS satellites with different frequency factors, the receiver code bias d rR , as well as phase delay b rR , are different. Their differences are usually called inter-frequency biases (IFB). The inter-system and inter-frequency biases must be considered in a combined processing of multi-GNSS observations. The code bias parameters are setup for each system and each frequency of GLONASS. In order to eliminate the singularity between receiver clock and code bias parameters, the code bias for GPS satellites is set to zero. This means that all estimated biases of other systems are relative to the biases for the GPS satellites. These estimated biases can be interpreted as a relative calibration of ''BeiDou/Galileo with respect to GPS'' and ''each individual frequency used by a GLONASS satellite with respect to the GPS frequency''. It is worthwhile to notice that such a receiver internal bias is only relevant for processing the code data. When analyzing the phase measurements the corresponding phase ambiguity parameters will absorb the phase delays. They become only relevant if ambiguities between different GNSS are resolved to their integer values.
For precise positioning users, precise satellite orbit, clock and DCB products (e.g. from MGEX products) need to be applied. With these products, the positioning model then can be simplified as, In multi-GNSS PPP, the parameters to be estimated in the combined processing contain the receiver position increments r r , receiver clock bias t r , zenith tropospheric www.nature.com/scientificreports wet delay Z r , slant ionospheric delay I s r,1 , phase ambiguity N s r,j , and the system/frequency dependent code biases in the receiver end, i.e. d rE , d rC and d rR relative to the GPS biases d rG . The phase delays b r and b 5 will be absorbed by phase ambiguity parameters. The estimated parameters can be expressed as, In order to strengthen the solution, a priori knowledge of the ionospheric delays including the temporal correlation, spatial characteristics and external ionospheric model is also utilized to constrain the estimated ionospheric parameters 22,23 . These constraints, to be imposed on observations of a single station can be summarized as, where t is the current epoch and t 2 1 is the previous epoch; w t is a zero mean white noise with variance s 2 wt ; vI s r is the vertical ionospheric delay with a variance of s 2 vI ; f s r,IPP is the mapping function at the ionospheric pierce point (IPP); the coefficients a i describe the trend; dL and dB are the longitude and latitude difference between the IPP and the station location;Ĩ s r is the ionospheric delay obtained from external ionospheric model with a variance of s 2 I . Table 3 summarizes our multi-GNSS data processing strategy in detail. In our PPP model, all the observations from different GNSS (four systems) are processed together in one common least square estimator. The raw-observation model with ionospheric constraints of equation (12) is adopted to improve the PPP performance. The receiver positions are estimated in both static and kinematic modes. The tropospheric zenith wet delay Z r is described as a random walk process with noise of about 225 mm . ffiffiffiffiffiffiffiffiffi ffi hour p . The receiver clock is estimated epoch-wise as white noise. The carrier-phase ambiguities are estimated as constant for each arc. The ISB and IFB parameters are estimated as constant over time and GPS is selected as reference.