A variant of raw observation approach for BDS/GNSS precise point positioning with fast integer ambiguity resolution

The Precise Point Positioning (PPP) technique uses a single Global Navigation Satellite System (GNSS) receiver to collect carrier-phase and code observations and perform centimeter-accuracy positioning together with the precise satellite orbit and clock corrections provided. According to the observations used, there are basically two approaches, namely, the ionosphere-free combination approach and the raw observation approach. The former eliminates the ionosphere effects in the observation domain, while the latter estimates the ionosphere effects using uncombined and undifferenced observations, i.e., so-called raw observations. These traditional techniques do not fix carrier-phase ambiguities to integers, if the additional corrections of satellite hardware biases are not provided to the users. To derive the corrections of hardware biases in network side, the ionosphere-free combination operation is often used to obtain the ionosphere-free ambiguities from the L1 and L2 ones produced even with the raw observation approach in earlier studies. This contribution introduces a variant of the raw observation approach that does not use any ionosphere-free (or narrow-lane) combination operator to derive satellite hardware bias and compute PPP ambiguity float and fixed solution. The reparameterization and the manipulation of design matrix coefficients are described. A computational procedure is developed to derive the satellite hardware biases on WL and L1 directly. The PPP ambiguity-fixed solutions are obtained also directly with WL/L1 integer ambiguity resolutions. The proposed method is applied to process the data of a GNSS network covering a large part of China. We produce the satellite biases of BeiDou, GPS and Galileo. The results demonstrate that both accuracy and convergence are significantly improved with integer ambiguity resolution. The BeiDou contributions on accuracy and convergence are also assessed. It is disclosed for the first time that BeiDou only ambiguity-fixed solutions achieve the similar accuracy with that of GPS/Galileo combined, at least in mainland China. The numerical analysis demonstrates that the best solutions are achieved by GPS/Galileo/BeiDou solutions. The accuracy in horizontal components is better than 6 mm, and in the height component better than 20 mm (one sigma). The mean convergence time for reliable ambiguity-fixing is about 1.37 min with 0.12 min standard deviation among stations without using ionosphere corrections and the third frequency measurements. The contribution of BDS is numerically highlighted.


Introduction
with precise orbit and clock corrections without using direct and explicit reference sites (Kouba and Héroux 2001). The PPP technique started with the Ionosphere-Free (IF) combination approach and then developed with the raw observation approach (Wübbena et al. 2005;Teunissen et al. 2010;Li et al. 2013;Chen and Zhao 2014;Zhang et al. 2019). The raw observation approach directly uses uncombined and undifferenced observations, therefore the ionosphere delays are estimated. While the IF combination approach eliminates the ionosphere delay in the observation domain.
Integer Ambiguity Resolution (IAR) is a key for PPP technique to improve the accuracy and shorten the convergence time. The ambiguity term estimated with the IF approach is a combination of Narrow-Lane (NL) and Wide-Lane (WL) ambiguity (Ge et al. 2008), while the raw observation approach produces L1 and L2 ambiguity parameters simultaneously together with the ionosphere term. To fix the ambiguities of a single receiver, additional corrections, i.e., satellite hardware biases (or so-called uncalibrated phase delays, or fractional cycle biases), are needed to be calibrated in a reference network. This is often achieved through merging the ambiguity terms derived from individual sites of the network into a single set of corrections that can be applied to the entire network. The IF combination approach uses two processes (Ge et al. 2008;Laurichesse et al. 2009;Collins et al. 2010;Geng et al. 2012). One process is to compute the geometry-free (GF) Melbourne-Wübbena combinations (Melbourne 1985;Wübbena 1985) to obtain WL satellite phase biases and integers of those ambiguities. The other is to carry out the Geometry-Based (GB) positioning to compute the IF ambiguity terms and further separate NL biases with the involvement of the integers of GF WL ambiguities. The raw observation approach uses only one process that produces GB L1 and L2 ambiguities simultaneously, which are then used to form GB WL ambiguities and further produce WL satellite phase biases. Due to the fact that the ambiguity terms are highly correlated with the ionosphere terms, an IF combination of L1 and L2 ambiguities is used to form GB NL ambiguities in this approach (Li et al. 2013;Gu et al. 2015). The same as the IF approach, the NL satellite hardware biases are then obtained in this step. The difference is that one estimates the NL ambiguities using the IF combination observations, while the other forms the NL ambiguities through the estimates of L1 and L2 ambiguities. In Teunissen et al. (2010) and Zhang et al. (2019), the satellite hardware biases are estimated together with satellite clock biases, ionosphere delays, and L1 and L2 ambiguities.
This contribution develops a variant of the raw observation approach. A reparameterization is applied to the design matrix coefficients and the parameter domain, while the observations, as input, are undifferenced and uncombined. A computational procedure is developed at the network side to estimate WL and L1 satellite hardware biases as corrections for the IAR purpose, without using the IF combination operation (Zhang et al. 2019;Li et al. 2020). The PPP IAR solutions are therefore obtained at user side with incorporating of these corrections.
As is well known that there are currently four constellations of GNSSs, namely American GPS, Russian GLO-NASS, Chinese BeiDou and European Galileo. Among them, GPS, BeiDou and Galileo systems are based on the same Code-Division Multiple Access (CDMA) technique. The presented PPP-IAR method can be directly employed for these three systems. The GLONASS system is not considered in this contribution for a concern of loss of focus. It is noticed that the European Galileo system has been achieving its nearly full capability by early 2020 with 22 medium earth orbit satellites in normal operation (Falcone et al. 2017) (also https://www.gsc-europa.eu/ system-service-status/constellation-information accessed on August 11, 2021). Another two satellites (i.e. E14 and E18) were launched into incorrect orbit, but they were moved to usable orbit later in 2014 and 2015. The navigation messages of these two satellites are available since August 2016. The two satellites are fully usable when the precise orbit and clock products are available (Paziewski et al. 2018).
In the meanwhile, the BeiDou system has been emerged since the BDS-3 (BeiDou Phase 3) primary system was announced to provide global services on December 27, 2018 (Yang et al. 2017(Yang et al. , 2020. The constellation of BeiDou system currently includes 7 Geostationary Earth Orbit (GEO) satellites, 10 Inclined GeoSynchronous Orbit (IGSO) satellites, and 27 Medium Earth Orbit (MEO) satellites in normal and healthy operation (http:// www.csno-tarc.cn/en/system/health accessed on August 11, 2021) in combining with BDS-2 (BeiDou Phase 2) satellites. Earlier studies estimated the fractional cycle bias for GPS/BDS-2/Galileo based on the international GNSS monitoring network and assessed the performance of PPP-IAR using the IF combination approach Hu et al. 2020) and using the raw observation approach (Wang et al. 2021). However, the BDS-3 was not finalized at that time and the performance with its combination with BDS-2 was not evaluated yet. This provides an opportunity to further investigate the contribution of BDS-3 together with BDS-2 satellites. In addition, it also makes sense to evaluate the contribution of the Galileo to the accuracy and convergence of PPP IAR solutions. As most of GNSS users still use the dual-frequency receivers, the focus of this contribution is on evaluation of the performance of dual-frequency measurements, though these three GNSS constellations provide signals of three or five frequencies.
The proposed method is applied to process the data of a GNSS network covering a large part China. The spacing of the sites is about 160 -280 km, which is much longer than that of traditional network RTK (Real-Time Kinematic). The WUM (Wuhan University Multi-GNSS) precise orbit and clock products are used for hardware bias estimation and PPP positioning. The L1 and WL satellite biases of BeiDou, GPS, and Galileo are produced and applied to user sites. The ambiguity-float and -fixed solutions are produced to evaluate the contribution of Galileo and BeiDou system. The paper is organized as follows. "Section Methodology" introduces the algorithm of the raw observation approach and its variant. The hardware bias derivation and PPP-IAR methodology are also described. The experimental network and the used orbit and clock products are detailed in "Section Experiment and result analysis". The hardware biases are demonstrated afterwards. The positioning performance in terms of precision and convergence time are analyzed in the section as well. The conclusions and outlook are provided in "Section Conclusions and outlook".

Methodology
This section describes the methodology of the proposed PPP algorithms, hardware bias derivation, and PPP integer ambiguity fixing. We start with the original raw observation approach on the basis of Teunissen et al. (2010), Li et al. (2013), Odijk (2017) and Zhang et al. (2019).

Raw observation approach
The observation equation for code p and carrier phase φ measurement tracked by a receiver at location r from a satellite s (s = 1, 2, . . . , m) of a constellation g, indicating BDS, GPS and Galileo system, at frequency i (i = 1, 2) can be expressed in meters as follows (Kleusberg and Teunissen 1996;Hauschild 2017): where E{·} is the expectation operator. ρ s,g r is the geometry distance between the receiver and the satellite. T r is the zenith tropospheric delay with its corresponding mapping function m s,g r , while δt r and δt s,g are receiver and satellite clock bias in meters, respectively. I s,g r the ionosphere delays with a factor of γ and d s,g φ i are phase hardware biases on receiver and satellite, respectively. It is assumed that these hardware biases are relatively stable in time.
To obtain precise point positioning, precise orbit and clock products are applied. The raw observation approach directly uses the uncombined and undifferenced measurements. The receiver and satellite code biases are not estimable in a single receiver, a reparameterization is therefore needed, which implies an application of the S-basis transformation (Teunissen 1985;Odijk et al. 2016). Taking the fact that satellite clock product containing part of satellite code biases into account, the raw observation equation of GB model can be formed as follows (Zhang et al. 2012;Odijk 2017): where δp s i,r and δφ s i,r are Observed Minus Computed (OMC) values, which are calculated from raw observations and a priori coordinates as well as precise satellite orbit and clock. µ s,g r is the unit vector of LOS (Line-Of-Sight) between the receiver and the satellite, while ∆x r the receiver positioning increments with respect to the a priori coordinates. The receiver clock parameter in Eq.
(2) is formed as: and the ionosphere term is parameterized as: The ambiguity terms containing part of code hardware delays in Eq. (2) are given as follows: It should be noticed that the ambiguity terms have no integer nature as they are not separated with the receiver and satellite hardware biases. For the detail of derivation of Eqs. (3)-(5), one can refer to Odijk (2017). In order to achieve integer ambiguity-fixed solution, the (3) t g r = δt r + γ ,g 2 γ ,g satellite biases have to be computed in a regional or global network.
The stochastic model for the observations of one satellite is given as follows: where e s,g r is the elevation of the satellite, σ (e s,g r ) p and σ (e s,g r ) φ are a priori formal errors for code and phase measurements, respectively. The elevation-dependent data weighting is applied. Assume that they are the same for L1 and L2 frequency. Their correlations between code and phase, as well as between L1 and L2 measurements, are assumed sufficiently small and therefore are not considered.
Compared with the IF combination approach, the raw observation approach estimates the ionospheric parameters instead of eliminating them. This gives an opportunity to constrain the ionospheric parameters when their corrections are provided, resulting in stronger geometry strength and potentially faster convergence. However, the raw observation approach produces almost the same positions as those with the IF combination approach, if the ionosphere parameters are not constrained and the stochastic model and the process noise in parameters are properly considered the similar in both cases.

A variant of raw observation approach
The raw observation equation shown in Eq. (2) can be further manipulated only for phase observation of the second frequency: where the WL ambiguity term is defined as follows: which includes a WL integer ambiguity plus a receiver and satellite wide-lane biases. The term ,g 2Ñ s,g 1,r was added as an additional parameter, while its negative counterpart is combined with the previous term ,g 2Ñ s,g 2,r . The widelane integer ambiguity is the difference between L1 and L2 integer one, as defined in Allison (1991), Chen et al. (2018), Liu et al. (2021): Note that the WL ambiguity has the L2 wavelength as coefficient, rather than the WL wavelength in Ge et al. (2008), Laurichesse et al. (2009), Collins et al. (2010. The WL receiver hardware bias is parameterized as a (6) Q y s,g r = diag{σ (e s,g combination of code and phase receiver hardware biases on L1 and L2: and the WL satellite hardware bias is parameterized as a combination of code and phase satellite hardware biases on L1 and L2: The L1 ambiguity term in Eq. (7) includes the integer plus a receiver and satellite biases.
where the receiver bias is formed as: similarly the satellite bias is formed as follows: The observation equations for all observed satellites of one constellation in the matrix-vector form is expressed as: where the code OMC vector P g = δp I m is a unit matrix and e 2m = [1 1 · · · 1] T . ⊗ is the Kronecker product operator. The stochastic model is expressed as: where the code and phase part are, respectively, denoted as and The variant of the approach is equivalent to the original raw observation approach. The only difference is that the WL ambiguity term replaces the L2 ambiguity term via a slight manipulation of the elements of the design matrix. All other parameters are the same as what they were in the original raw observation approach. In other words, the variant produces the equivalent positioning as the original raw observation approach, if the same stochastic model and computation mode are used.

Hardware bias estimation
In order to produce satellite hardware biases corrections, the positions of the reference sits are considered known, the position increments ∆x r and their elements in design matrix, i.e., unit vector, are vanished. The terms b s,g wl in Eq. (11) and b s,g 1 in Eq. (14) are therefore the quantities to be extracted from the ambiguities obtained at the individual sites. The advantage of the manipulation is directly estimating the WL and L1 ambiguities and their variance-covariances, without forming WL and IF (or NL) combination ambiguity after estimating the L1 and L2 ambiguities as done in Li et al. (2013) and Gu et al. (2015). It is also possible to directly estimate the single-differenced WL/L1 ambiguity term (to eliminate the receiver hardware biases) through a further reparameterization shown as follows: the last satellite m is assumed as the reference satellite, e.g. with the highest elevation, s = 1, 2 . . . , m − 1 in sequel unless otherwise specified.
Through a least-square estimation, the single-differenced WL and L1 ambiguity terms of one satellite are obtained together with their variance-covariance When the formal error σˆÑ sm,g wl,r is smaller than a threshold (e.g., 0.2), the fractional part of single-differenced WL ambiguities is therefore determined as follows: where · the ambiguity fixing operator, e.g., a rounding operator. Assume that the ambiguities are correctly rounded in the network, a sm,g wl,r is therefore taken as an input to estimate the satellite WL phase biases. The functional model for estimating the satellite WL phase biases can be formed as follows: where a g wl,r = a 1m,g wl,r a 2m,g wl,r · · · a m−1m,g wl,r T with r = 1, 2, . . . , n being the reference sites of the network. The WL bias of a (reference) satellite is assigned to an as a datum, adding into Eq. (23) in order to solve the rank defect issue.
Equation (23) assumes that all reference sites are tracking the same satellites, and the same reference satellite is chosen. However, this is not in the reality. Even though, it is straight forward to handle. We just need to replace 1 to 0 for the corresponding elements in design matrix when one satellite does not appear in one site. Each individual site can select its own reference satellite. We just need fill in -1 to the corresponding column when a different satellite is selected as the reference satellite for that site. More generic theory on how to handle the connectivity of ambiguities in the network can be found in Khodabandeh and Teunissen (2019).
The estimated satellite phase biases b s,g wl and b m,g wl from the network are used to correct all float WL ambiguities for all sites as the following: The WL integer ambiguities ǎ sm,g wl,r are therefore determined on the basis of the corrected float ambiguities and its formal error σˆÑ sm,g wl,r . Ideally, this can be rigorously achieved with the all WL ambiguity terms at all individual sites with their full variance-covariances using the LAMBDA (Least-square AMBiguity Decorrelation Adjustment) method to search for the integer ambiguity candidates (Teunissen 1995). However, the system complexity and computation burden might exponentially increase for the network solution.
The differences of the WL float and fixed ambiguities are used to correct the L1 ambiguities estimated previously together with their variance-covariance elements: The variance of the L1 ambiguity is also improved in this case.
This step is different from that done in Li et al. (2013) and Gu et al. (2015), where the IF or NL combination operator is applied to form a NL ambiguity term. This contribution considers the correlation of L1 and WL ambiguity terms. The real values of L1 estimates as well as their variances are virtually corrected by their correlation with the WL estimates after the WL satellite phase biases are solved. This is actually an extension and application of so-called integer bootstrapping technique (Teunissen 2001) to the estimation of satellite phase biases.
The formal error σˆÑ sm,g 1|wl,r is also used to determine the fractional part of single-difference L1 ambiguities.
The value a sm,g 1|wl,r is taken as an input to estimate the satellite L1 phase biases. The similar functional model as Eq. (23) is described in the following: T . Again, an arbitrary datum b m,g 1,0 has to be added in Eq. (28) in order to solve the rank defect issue. The estimates of L1 satellite phase biases b s,g 1 and b m,g 1 , together with estimated WL satellite phase biases in previous step, are ready to broadcast for users to have integer ambiguity resolution.

PPP IAR solutions
The parameters of the functional model in the user side can be divided into two groups, as follows: The single-differenced ambiguity terms Ñ sm,g 1,r and Ñ sm,g wl,r are included in the group of x 2 , while the ambiguity terms of the reference satellite Ñ m,g 1,r and Ñ m,g wl,r are included in the group of x 1 with other remaining parameters. The corresponding stochastic model is the same as Eq.16. The least-squares estimation is carried out and the ambiguity float solution is therefore obtained as: The single-differenced ambiguity term x 2 can be directly corrected with the WL and L1 satellite hardware biases. This recovers the integer nature of the ambiguities. In this case, the integer searching and fixing scheme, e.g. the LAMBDA method, is triggered. The procedure is applied to the corrected x 2 and Qx 2 ,x 2 . Assume that they are fixed to their integer ambiguities x 2 . The ambiguity fixed solution of the remaining parameters is then given by The advantage of fixing WL and L1 ambiguities, compared with fixing of L1 and L2 ambiguities, is less computational efforts. The WL part includes already an explicit decorrelation that might have to be derived with L1 and L2 ambiguities (Chen and Zhao 2014;Zhu et al. 2021). However, This decorrelation might be only a start between the frequencies. The correlation among satellites has to be reduced for both cases. It is also possible to first fix the WL ambiguities and then fixing the L1 ambiguities through the integer bootstrapping, as we did for the network side. However, the highest success rate of ambiguity resolution can be achieved in theory by resolving the full set of ambiguities simultaneously through the integer least-squares estimator (Teunissen 1999).

Experiment and result analysis
In this section, the experimental data and processing strategy are introduced. The satellite hardware biases produced with the proposed method are first evaluated. The positioning performance is then analyzed in terms of both precision and convergence time.

Experimental measurements
In order to test the methodology and processing strategy, we collected and processed GNSS measurements with 2 s sampling interval at 54 selected sites covering large part of China (mainly in North, South and East China) from Wuhan University's experimental network. These sites are equipped with Unicore UB4B0-MAX geodetic receivers with Dywell GNSS miniaturized antenna. The receiver UB4B0 series support to track multi-frequency signals of BDS, GPS, GLONASS and Galileo navigation systems. The antenna applicable signal bands are GPS L1/L2, GLONASS L1/L2, BDS B1/B2/B3 as well as Galileo E1/E5. The spacing among the sites are about 160 -280 KM. We selected 32 sites as reference stations to derive satellite hardware biases and the rest of 22 sites as users to evaluate the performance of positioning. The site distribution can be seen in Fig. 1. The observations of DOY 127, 2021 were used in the experiment.

WUM orbit and clock products
Orbit and clock are primary elements to compute PPP solutions. In this contribution, the orbit and clock products are produced by Wuhan University as an IGS (International GNSS Service) analysis center. The computation strategy used for the WUM (Wuhan Multi-GNSS) products was presented in Guo et al. (2016). Software Position And Navigation Data Analyst (PANDA) is applied for data analysis (Liu and Ge 2003). Generally, three steps are used. Firstly, the GPS and GLONASS data with 300 s interval from IGS network are processed to estimate the orbits, clocks, Earth Rotation Parameters (ERP), station coordinates and other parameters. The clocks with 30 s interval are determined further by using the epoch-differenced phase observations. In the second step, the GPS-only PPP is applied to determine the station coordinates, epoch-wise receiver clock offsets and Zenith Troposphere Delay (ZTD) parameters with two hour interval, and others for each station with the daily data from IGS and iGMAS networks by fixing the orbits, 30 s interval clocks, and ERP determined previously. Finally, Galileo, BDS, and QZSS orbit and 30 s interval clock are determined with the data from the same stations. In this step, the station positions, ZTD, and receiver clocks obtained in the previous GPS-only PPP are introduced as known parameters. As stated in Guo et al. (2016), we follow the recommendations of the 2nd IGS reprocess campaign (see http://acc.igs.org/ reprocess2.html) for the specific measurement mode, reference frame, orbit model, and parameters to be estimated. The external and independent assessments from the point of user side indicates the WUM BDS-2 product has the comparable quality as other analysis centers produced (Liu et al. 2019). The analysis by Steigenberger and Montenbruck (2020) shows the good quality of WUM products and the best consistency between WUM and GBM Galileo and BDS orbit products among all MGEX

Hardware bias products
BDS, GPS and Galileo measurements from those 32 sites in Fig. 1 are used to derive satellite hardware biases. The coordinates of these reference stations are fixed as the "ground truth", which are calculated from one day's data in the static mode. The orbit and clock are interpolated for each epoch (2 s of data interval is used) to produce the OMC values. The design matrix is formed as shown in Eq. (19). A preprocessing procedure is applied in single satellite basis to detect large cycle slips (Liu 2016). A quality control procedure is applied to reject the outliers and small cycle slips on the basis of the geometry-based residuals and their variance-covariance (Teunissen 1998). These two procedures are crucial to obtain satisfactory results in this study. The WL and L1 ambiguities are estimated for individual sites through a forward computation. In parallel, those ambiguities are synchronized and used to estimate satellite hardware biases as described in Eqs. (23) and (28). Figure 2 displays the time series of hardware biases for BDS GEO/IGSO satellites. They are shifted by an integer cycle within [− 3 3] to distinguish them from each another, while their variations are kept relatively visible (the same has been done for all the time series of hardware biases in sequel figures). The GEO satellites are always tracked in the region and their time series of the biases are complete during the entire day. For IGSO satellites, there are no bias presented for about 3.5-5 h, as they are not sufficiently tracked (above certain elevations) at all sites in those periods. It is obvious that the WL biases are much more stable than the L1 ones. Most of these satellites' WL biases vary about 0.2 cycles, while the L1 biases vary about 0.2-0.7 cycles, except for C05 having about 1 cycle variations. At the beginning of their presents in the network, the variations seem slightly larger. In general, the longer a satellite is tracked in the network, the more stable the hardware biases are over the period. The larger bias variation may affect the ambiguity fixing for the satellites at the beginning, but does not affect the positioning results when there are sufficient number of MEO satellites with ambiguity-fixing. It should be also mentioned that the hardware biases do not need so stable as a constant. The biases absorb some LOS orbit errors (Chen et al. 2018;Liu et al. 2021), particularly for regional networks, the LOSs of individual satellites for users in the network are similar as the reference sites. Therefore the corrections have an advantage more than just correct the hardware biases, but also help users to remove the residual errors of applied precise orbits. Therefore, the regional satellite biases may be better in improving the convergence speed and precision than global ones. The similar observations were obtained in a GPS only case (Wang et al. 2018). Figure 3 demonstrates the time series of hardware biases for BeiDou MEO satellites. There are 27 MEO satellites having the corresponding biases estimated. Most of MEO satellites are BDS-3 except for C11, C12 and C14 being BDS-2 satellites. C59 and C60 are tracked by the receivers of the network, but their orbit and clock products are not available yet, and their hardware biases are not calculated. Again, the WL biases are much more stable than the L1 ones. Most of these satellites' WL biases vary about 0.1 cycles, while the L1 biases vary about 0.1-0.3 cycles, except for three BDS-2 satellites (WL biases of these satellites indicated in lightblue, red and green curves change smoothly within 0.3 cycles while the change range of L1 counterparts is about 1.0 cycles), they vary more significantly than the rest of BDS-3. It is likely that the BDS-2 orbits and clocks are not as good as that of the BDS-3 satellites, therefore, part of orbit and clock errors are absorbed into the hardware biases. The advantage of BDS constellation in our experimental network is that there are generally more than 20 satellites with hardware biases provided at all time. Figure 4 displays the time series of hardware biases for GPS satellites. There are 31 MEO satellites with satellite hardware biases provided. The variation behavior of GPS WL biases is similar as that of BDS-3 MEO satellites, while the GPS L1 biases are slightly more stable than that of BDS-3 satellites. It is noticed that the hardware biases availability of BDS-3 MEO satellites are, in general, about one hour longer than the GPS satellites, probably due to their different orbit periods. This potentially gives users more ambiguity-fixed satellites at least in this regional network. Figure 5 demonstrates the time series of hardware biases for Galileo satellites. There are 24 MEO satellites having the corresponding biases estimated in this study (including E14 and E18). Compared with GPS and BDS MEO satellites, the hardware biases (both WL and L1) of Galileo satellites have better stability. This could be two reasons: one is that the hardware biases of Galileo satellites are actually more stable. The other is that Galileo

Fig. 3 Time series of hardware WL (top) and L1 biases (bottom) for BeiDou MEO satellites
orbit and clock may also be more accurate than GPS and BDS MEO satellites, therefore, less orbit and clock residuals are absorbed in their biases, a similar behavior observed in Li et al. (2018), Li et al. (2020) as well. It is also noticed that the hardware bias availability of most Galileo satellites is about one hour longer than the GPS satellites and similar as the BDS MEO satellites. The Galileo has the same advantage as the BDS constellation in this aspect.

Positioning performance
The positioning performance is investigated in terms of their precision and convergence with the hardware biases being applied.

Precision evaluation
The raw dual-frequency measurements of GPS, Galileo and BeiDou satellites for 22 user receivers are processed with WUM precise orbit and clock products. Besides the kinematic coordinates, the parameters to be estimated are tropospheric zenith delay, slant ionospheric delay, receiver clock errors of the used constellations, and WL and L1 ambiguities. The satellite elevation mask is set to degree 5 and the elevation dependent weighting is applied. The used correction models are referred to Kouba and Héroux (2001). The standard PPP ambiguityfloat solutions are calculated without using hardware bias corrections. While ambiguity-fixed solutions are only obtained when the hardware biases are applied to correct the single-differenced ambiguity term as indicated in Eq. (29) and the sufficient number of ambiguities are successfully fixed. The modified LAMBDA (Chang et al. 2005) is applied to fix the ambiguities into their integers with a Partial Ambiguity Resolution (PAR) scheme based on Cao et al. (2007), Wang and Feng (2013), Wang and Verhagen (2014), Li and Zhang (2015). The processing scheme attempts to fix ambiguities as many and fast as possible. Three constellations' measurements are separately or combined employed to compute PPP/PPP-IAR solutions for comparing the performance of three constellations. The positioning performance is investigated with the ambiguity-float and ambiguity-fixed solutions. Figure 6 demonstrates the standard deviations of GPS only PPP and PPP-IAR kinematic solutions for each individual user site for an entire day (the position biases are negligible as their ground truth are calculated using the same data). The positions of the first 20 minutes of ambiguity-float solutions are not included in the statistics (the same for next two figures). The overall standard deviations of ambiguity-float solutions for all these sites are   Figure 7 displays the standard deviations of BeiDou only ambiguity-float and -fixed PPP kinematic solutions. The overall performance of ambiguity-float solutions is slightly worse than that of GPS solutions. The average standard deviations are 3.11, 3.14 and 5.31 cm for North, East and Height components, respectively. However, the ambiguity-fixed solutions are 0.68, 0.73 and 2.39 cm for three components, which are significantly better than that of GPS counterparts. The improvements with respect to the float solutions are 78.3%, 76.8% and 55.1%, respectively for three components, which is rather impressive. It is interestingly noticed that the BDS float solutions have almost equal precision for two horizontal components, while GPS or Galileo ambiguity-float solutions often have much better precision in North than in East. This is obviously the contribution of GEO and IGSO satellites to such a phenomenon in favor of BeiDou system. According to the statistics, there are averagely 86.6% ambiguities are fixed to their integers in each epoch, slightly higher than that of GPS. The average number of the used Bei-Dou satellites is 21.6 in the entire day, about two and half times of that in the GPS solutions. This means that in each epoch there are averagely about 36 single-differenced ambiguities fixed, which significantly contributes to the improvement of the positioning accuracy. Figure 8 plots the standard deviations of Galileo only PPP and PPP IAR kinematic solutions. The performance of Galileo seems the worst among three systems. There are a couple of sites which have relatively larger standard deviations (e.g. site 4596, site 4682 and hb06). The overall standard deviations of ambiguity-float solutions for all 22 sites are 1.96, 3.88 and 5.54 cm for North, East and Height components, respectively, while the ambiguityfixed solutions are 0.99, 1.17 and 3.36 cm. The improvements are 49.9%, 69.9% and 39.4%, respectively for three components. Our statistics shows that there are in average 83.5% ambiguities are fixed to their integers in each epoch. The average number of the used Galileo satellites

Fig. 5 Time series of hardware WL (top) and L1 biases (bottom) for Galileo satellites
is 7.8, about 1.5 satellite less than in the GPS case. The situation will be improved with the future deployment of Galileo satellites. We also calculate the GPS/Galileo and GPS/Galileo/ BeiDou combination solutions. The statistics of these solutions are included in Tab. 1. The precision of these combination float solutions is significantly improved compared with that of single constellation solutions. It is noticed that BeiDou only PPP IAR solutions have the similar accuracy with GPS/Galileo PPP IAR solutions. Both types of solutions have 7 mm standard deviations in horizontal components and better than 25 mm in height component. The improvement room becomes smaller when the BeiDou system joins GPS and Galileo for the float solutions, which seem hardly improved compared with GPS/Galileo float solutions (the horizontal precision seems even slightly worse, while the height component is still obviously better). This reminds us that more efforts should be made to improve particularly the stochastic model for BeiDou measurements to improve the performance of the horizontal components for float solutions. The ambiguity-fixed GPS/Galileo/BeiDou combination solutions are obviously improved compared with GPS/Galileo ambiguity-fixed solutions. The accuracy for horizontal components is better than 6 mm, and for the height component better than 20 mm. This is rather promising.

Convergence assessment
This section is to assess the convergence time of ambiguity-fixed solutions. To do so, all the parameters in the PPP estimation are completely reset every two hours to simulate the loss of signal tracking in receivers. This means the data are divided into 2-hour sessions and processed independently. The assessment is to account the time used for the first reliable ambiguity fixing for all the 2-hour sessions at all user sites. The criteria of reliable ambiguity fixing is that (1) sufficient number of the ambiguities (e.g., eight ambiguities per constellation) are fixed to their integers, (2) the horizontal and vertical positioning accuracies are respectively better than 5 cm and 10 cm, and they are stable for at least 10 min. Figure 9 displays the time series of GPS/Galileo/BeiDou combination PPP and PPP-IAR solutions for Site 4843. The site is located in the Southern China, therefore, the ionosphere is relatively more active than most of other sites. The results of this site can be a representative to demonstrate the convergence performance of three-system combination solutions. The result looks rather promising. Some of the float solutions are converged quite fast, while the  In order to display the details of the convergence and ambiguity-fixing, Fig. 10 zooms in the time series of 30 min after one of the resets (at 6:00 h as an example). The GPS only, and GPS/Galileo solutions are also displayed together with the GPS/Galileo/BeiDou ones for a comparison. The PPP-IAR solutions are not available at the beginning of the reset, when there are no sufficient number of the ambiguities being fixed. As can be seen, the reliable PPP-IAR solutions start at about 9 min after the reset for GPS only case, about 4 min for GPS/Galileo, and 1.5 min for GPS/Galileo/BeiDou case. The PDOP values are above 1.5 for GPS only, below 1.5 but larger than 1.0 for GPS/Galileo solutions. While the PDOP values are well below 1.0 for GPS/Galileo/BeiDou solutions. The smaller value the PDOP is, the larger the geometric strength will be. This leads to faster convergence. Figure 11 zooms in the time series of 30 min after one of the reset(at 4:00 h as an example) for BeiDou only, Bei-Dou/GPS and BeiDou/Galileo. As can be seen, the reliable PPP-IAR solutions start at about 7 min after the reset for BeiDou only case, which is earlier than that of GPS only solutions. Both BeiDou/GPS and BeiDou/Galileo IAR solutions start about 4 min, comparable to the GPS/ Galileo IAR solutions. The PDOP values are larger than 1.0, but smaller than 1.5 for BeiDou only solutions. This is better than the GPS only solutions (higher than 1.5). It is understandable that the number of BeiDou satellites is many more than that of GPS satellites. The BeiDou/GPS and BeiDou/Galileo PDOP values are close each other. In this case, their convergence performance is also similar. Figure 12 demonstrates the average convergence times and their corresponding variabilities at all sites for GPS only, GPS/Galileo and GPS/Galileo/BeiDou solutions, respectively. As can be seen that the site average values of convergence time for GPS only solutions are between 10 and 12 min (mean is 11.06 min), rather similar for all sites except for hb17, which is longer than 12 min. The average-variability is about 1.54 min in standard deviation. In other words, the GPS only solution needs about 14 min to achieve reliable ambiguity solutions at 95% confidence level. This is widely achieved in earlier studies. For the GPS/Galileo combination solutions, the ambiguity-fixed solutions are achieved in 4.15 min with 0.46 min standard deviation. This is certainly much faster than the GPS only solutions. Regarding to the GPS/Galileo/BeiDou solutions, the convergence time is about 1.37 min with 0.12 min standard deviation, which tells that 95% resets can be converged at 1.6 min, i.e., within 100 seconds. This is nearly real-time RTK solutions that needs additional atmosphere corrections. The contribution of BeiDou in this stage is remarkable (Fig. 13). In order to further assess the potential of BeiDou system, we computed BeiDou-only and BeiDou/Galileo and BeiDou/GPS solutions separately. Figure 13 demonstrates the average convergence times and their corresponding variations. The average convergence time for BeiDou-only solutions are between 6 -8 min (mean is 6.08 min). The average variation is about 0.49 min. This is significantly better than GPS-only solutions. The Bei-Dou-only solutions require only about 7 min to achieve integer ambiguity resolutions at 95% confidence level, rather than 14 min in the case of GPS only solutions. For the BeiDou/GPS combination solutions, the ambiguityfixed solutions are achieved in 3.40 min with 0.29 min standard deviation. This is also better than GPS/Galileo solutions. Regarding to the BeiDou/Galileo solutions, the

Conclusions and outlook
This contribution develops a variant of the original raw observation approach. The approach estimates the WL and L1 ambiguities via a reparameterization process, instead of L1/L2 or IF ambiguities. In the network side, the WL and L1 hardware biases are estimated, rather than WL/NL ones. These WL and L1 hardware bias corrections are broadcast to users, which are used to correct the WL and L1 ambiguities in user side and directly fix them to their corresponding integers. There is no any IF combination used in the above processing steps as adapted in the earlier studies. In order to validate the approach and assess the BeiDou contribution in mainland China, we processed the data from a regional network. The WL and L1 satellite biases of GPS, Galileo and BeiDou satellites are produced with the data from 32 sites with spacing of 160 -280 KM. The other 22 sites are taken as users to evaluate the precision of PPP and PPP-IAR solutions and the convergence time of the ambiguity-fixing solutions. The following main observations and conclusions can be obtained:    Sites Fig. 13 Average convergence time and its variability of BeiDou only (blue), BeiDou/GPS (green), BeiDou/Galileo (red) solutions eter precision can be achieved within 100 seconds, which is nearly real-time RTK. The contribution of BeiDou system is significant.
We have concluded the positioning performance of GPS, Galileo and BeiDou in terms of achieved precision and convergence time. In particular, the performance of Bei-Dou system and the capacities of the combined GPS/Galileo/BeiDou to achieve fast centimeter-level positioning are numerically highlighted for a regional network. We believe that future research can be considered to achieve even faster ambiguity-fixing in real-time by combining GLONASS and QZSS system, as well as three or five frequency measurements. In addition, the research could be extended to assess the BeiDou contribution to the positioning performance globally, in particular, beyond Asia-Pacific region where fewer GEO and IGSO satellites are tracked.