1 Introduction

Precise point positioning (PPP) is a high-precision positioning technique using a stand-alone receiver (Zumberge et al. 1997) and has been widely used for scientific research and civilian applications. As is well known, PPP ambiguity resolution (PPP-AR) (Gabor 1999; Ge et al. 2008; Laurichesse et al. 2008; Collins 2008) plays a key role in the development and it is becoming more and more reliable due to the sustained improvement of the precise satellite orbit and clock products (An et al. 2020). The principle of PPP-AR is that the uncalibrated phase delays (UPD) (Blewitt 1998) at both the receiver and satellite side are temporally stable, and thus, the satellite UPDs can be estimated precisely from a reference network and applied to a stand-alone receiver to recover the integer feature of the ambiguities (Gabor 1999).

Usually, the ionosphere-free (IF) combination is used in PPP processing in order to eliminate the first-order ionospheric delays; however, its ambiguity has to be separated into wide-lane (WL) and narrow-lane (NL) ambiguities for fixing (Dong et al. 1989; Blewitt et al. 1998). Therefore, it is also obvious that an IF ambiguity can be fixed only if both its WL and NL components are fixable. The WL ambiguity fixing is carried out first as the prerequisite of NL fixing because WL ambiguity can be easily estimated and fixed due to its long wavelength. It must be pointed out that the UPD is sensitive to the pseudorange and phase measurement noise in different frequency observation linear combinations (Duong et al. 2019; Xiao et al. 2019). The accuracy of UPDs is crucial for PPP-AR, as even a small bias in UPDs will contaminate the integer feature directly and result in wrong or failed ambiguity fixing (Li et al. 2018). All in all, both WL and NL UPDs are essential and must be estimated in advance in order to enable the PPP-AR at the user-side.

While being provided directly to users, the NL UPD can also be expressed implicitly, for example, merged into satellite clock corrections as decoupled or integer clocks (Laurichesse et al. 2008; Collins 2008). All the three methods are theoretically the same but with different expressions (Geng et al. 2010; Teunissen and Khodabandeh 2015), so we work on the direct way to estimate and provide the WL and NL UPDs in this study.

In practice, WL ambiguity is estimated directly from the Melbourne–Wübbena (MW) observation, which is a linear combination of the pseudorange and phase observations (Melbourne 1985; Wübbena 1985), while NL ambiguity is derived from the estimated IF ambiguity. Thanks to its long wavelength property and the continuous improvement of pseudorange accuracy, WL ambiguity can be well estimated and thus fixed reliably. Therefore, a lot of efforts have been made on improving the accuracy of estimating NL ambiguity, for example, in the way of reducing modeling errors in IF ambiguity estimation (Ge et al. 2008; Laurichesse et al. 2008; Collins 2008). However, the pseudorange observations could be biased with respect to different carrier frequencies, using different tracking methods and the various receiver realization of manufacturers. For example, it has been demonstrated that the differential code biases (DCB) exist among pseudoranges of different frequencies and code types (Lanyi and Roth 1988; Langley 1996), which prevents the use of mixed pseudorange types in WL fixing and UPD estimation (Ge et al. 2008). Fortunately, statistical studies showed that P1C1 DCBs can be estimated and corrected accordingly (Gao et al. 2001), and then, the remaining biases are negligible for WL UPD estimation with different types of receivers. In general, it is usually assumed that the UPDs at the receiver side can be removed by forming inter-satellite single difference (SD), so their characteristics at the receiver side are often overlooked (Gabor et al. 1999).

Along with the development of real-time precise positioning services, the contribution of pseudorange observation to PPP convergence is well recognized and more studies are focusing on the pseudorange quality (Fan and Ma 2014; Seepersad and Bisnath 2015), especially on the biases related to the code tracking performance depending on the specific satellites or configurations of the receivers (Hauschild and Montenbruck 2016; Liu et al. 2017). Hauschild and Montenbruck (2016) assessed the pseudorange difference in different satellites through zero-baseline tests and demonstrated its relating to the receiver type. Besides, the pseudorange biases can reach up to 20 cm when using the data collected from different receivers connected to a common antenna splitter. Based on the essential relationship between pseudorange bias and the non-ideal channel, Liu et al. (2019) found that different characteristics of the transmitting channel of a receiver also contribute to the distortion of pseudorange bias. Besides, the observable-specific signal bias is proposed for the pseudorange biases (Wang et al. 2020), but further investigation still needs to be carried out to take into account the different receivers.

Although such satellite-specified pseudorange biases are detected, neither their effect on WL UPDs estimation nor consequently their impact on PPP-AR has been investigated or evaluated so far, probably because such biases are simply considered to be small and negligible compared to the pseudorange noise and WL wavelength. However, such a bias of 20 cm will shift the corresponding WL ambiguity by about 0.25 cycles and makes its reliable fixing difficult and even impossible. Therefore, it is worthful to investigate the existence of such biases and to evaluate their impact on PPP-AR with the increasing number of receiver types in use. On the other hand, the huge amount of GNSS data provided by International GNSS Service (IGS) networks and other national/regional Continuously Operating Reference Station (CORS) networks makes such a study possible and reliable.

In this study, we focus on the WL UPDs estimation from a large number of GNSS stations with various receivers. Statistical analysis of the WL fractional part of the float ambiguity (FPA) and UPD are carried out carefully in order to investigate possible satellite-specified biases among different receiver types, respectively. Furthermore, the estimation of such receiver-type-dependent satellite-specified WL UPDs biases and their application as corrections to both reference stations and user stations are studied in order to reduce their impact on both UPD estimation and user-side ambiguity fixing, respectively. Finally, PPP-AR with and without such calibration is performed and compared to demonstrate the impact of this receiver-dependent deviation.

2 PPP-AR and UPD estimation model

We start with the basic GPS observational equations and then give a detailed description of the UPD estimation method as well as the K-Means clustering model (Krishna and Narasimha 1999) for classifying the fractional parts.

2.1 PPP using ionosphere-free observations

For satellite \(s\) observed by receiver \(r\), the IF observation equation of pseudorange \(P\) and carrier phase \(L\) observations can be expressed as,

$$ P_{{r,{\text{IF}}}}^{s} = \rho _{r}^{s} + c\left( {{\text{d}}t_{r} - {\text{d}}t^{s} } \right) + T_{r}^{s} + b_{{r,{\text{IF}}}} + b_{{{\text{IF}}}}^{s} + {\epsilon }_{{r,{\text{IF}}}}^{s}$$
(1)
$$ L_{{r,{\text{IF}}}}^{s} = \rho_{r}^{s} + c\left( {{\text{d}}t_{r} - {\text{d}}t^{s} } \right) + T_{r}^{s} + \lambda_{{{\text{IF}}}} \cdot \overline{N}_{{r,{\text{IF}}}}^{s} + \varepsilon_{{r,{\text{IF}}}}^{s} $$
(2)

where \({\rho }_{r}^{s}\) is the geometric distance, \(c\) is the speed of light, \({dt}_{r}\) and \({dt}^{s}\) are the clock offsets of receiver and satellite, respectively; \({T}_{r}^{s}\) is the slant troposphere delay; \({b}_{IF}^{s}\) and \({b}_{r,IF}\) are the code hardware and signal processing delays from the satellite signal generator to the satellite antenna and that from receiver antenna to the signal correlator in the receiver, respectively; \({\lambda }_{IF}\) is the wavelength of the IF carrier phase; \({\overline{N} }_{r,IF}^{s}\) is IF PPP ambiguity including the receiver and satellite UPDs; and \({\epsilon }_{r,IF}^{s}\) and \({\varepsilon }_{r,IF}^{s}\) are the measurement noise of pseudorange and carrier phase. Other items, such as the phase center offset (PCO) and variation (PCV) (Schmid et al. 2005), phase wind-up (Wu et al. 1993), relativistic effect, and tidal loading, are assumed to be precisely corrected with their corresponding models (Petit and Luzum 2010).

2.2 Estimation of UPD

For PPP-AR solution, \(\overline{N}_{{r,{\text{IF}}}}^{s}\) is usually decomposed into the following combination of integer WL \(N_{{r,{\text{WL}}}}^{s}\) and float NL \(\overline{N}_{{r,{\text{NL}}}}^{s}\) ambiguities for ambiguity fixing (Dong et al. 1989; Blewitt et al. 1998).

$$ \overline{N}_{{r,{\text{IF}}}}^{s} = \frac{{f_{1} }}{{f_{1} + f_{2} }}\overline{N}_{{r,{\text{NL}}}}^{s} - \frac{{f_{1} f_{2} }}{{f_{1}^{2} - f_{2}^{2} }}N_{{r,{\text{WL}}}}^{s} $$
(3)

where \({f}_{1}\) and \({f}_{2}\) are the carrier frequency of \({L}_{1}\) and \({L}_{2}\), respectively. It should be mentioned that here only integer WL ambiguity is introduced for easy realization, as its fractional part can be absorbed in the NL UPDs.

Using the MW combination at each epoch, we can obtain an observation of the WL ambiguity for epoch i:

$$ N_{{r,{\text{WL}}}}^{s, i} = \left[ {\left( {f_{1} L_{r,1}^{s,i} - f_{2} L_{r,2}^{s,i} } \right)/\left( {f_{1} - f_{2} } \right) - \left( {f_{1} \overline{P}_{r,1}^{s,i} + f_{2} \overline{P}_{r,2}^{s,i} } \right)/\left( {f_{1} + f_{2} } \right)} \right]/\lambda_{{{\text{WL}}}} $$
(4)

where \({\overline{P} }_{r,j}^{s}\) are the pseudorange observations which are corrected by satellite DCB. \({\lambda }_{WL}\) is the wavelength of WL. The estimation is calculated by taking the temporal average of the MW combination observables as,

$$ \overline{N}_{{r,{\text{WL}}}}^{s} = \langle{N_{{r,{\text{WL}}}}^{s, i}}\rangle = N_{{r,{\text{WL}}}}^{s} + d_{{r,{\text{WL}}}} - d_{{{\text{WL}}}}^{s} $$
(5)

where \(\overline{N}_{{r,{\text{WL}}}}^{s}\) is the WL float ambiguity, \(d_{{r,{\text{WL}}}}\) and \(d_{{{\text{WL}}}}^{s}\) are WL UPDs in receiver and satellite side, respectively. With all the float WL ambiguities of a network, the UPDs can be estimated at each epoch. The details will be discussed later on.

In case that both receiver- and satellite-side WL UPDs are already known, the WL ambiguity fixing can be attempted, and as soon as its fixing is successful, the integer \({N}_{r,WL}^{s}\) is introduced to Eqs. (3) to get the NL float ambiguity,

$$ \overline{N}_{{r,{\text{NL}}}}^{s} = \left[ {\lambda_{{{\text{IF}}}} \overline{N}_{{r,{\text{IF}}}}^{s} - cf_{2} N_{{r,{\text{WL}}}}^{s} /\left( {f_{1}^{2} - f_{2}^{2} } \right)} \right]/\lambda_{{{\text{NL}}}} = N_{{r,{\text{NL}}}}^{s} + d_{{r,{\text{NL}}}} - d_{{{\text{NL}}}}^{s} $$
(6)

where \(N_{{r,{\text{NL}}}}^{s}\) is the NL integer ambiguity, \(\lambda_{{{\text{NL}}}}\) is wavelength of NL, \(d_{{r,{\text{NL}}}}\) and \(d_{{{\text{NL}}}}^{s}\) are NL UPD in receiver and satellite side, respectively. As soon as we have all the float ambiguities at a single epoch, the UPDs can be estimated by an iterative procedure (Li et al. 2016; Zhang et al. 2020) as follows.

Each of the float ambiguity can be expressed by the integer ambiguity \({N}_{r}^{s}\) and the related UPDs as,

$$ \overline{N}_{r}^{s} = N_{r}^{s} + d_{r} - d^{s} $$
(7)

If we have an approximate value for each UPD, for example, taken from the fractional part of ambiguity at a single station assuming the receiver UPD is zero, we can try to get the integer ambiguity value according to its estimates from Eq. (7). The difference between the floating ambiguity and the integer at the current epoch for all station-satellite pairs can be combined to form a system of equations as,

$$ \left[ {\begin{array}{*{20}c} {B_{1}^{1} } \\ \vdots \\ {B_{1}^{s} } \\ {B_{2}^{1} } \\ \vdots \\ {B_{2}^{s} } \\ \vdots \\ {B_{r}^{1} } \\ \vdots \\ {B_{r}^{s} } \end{array} } \right] = \left[ {\begin{array}{*{20}c} 1 & 0 & \cdots & 0 & -1 & \cdots &0 \\ &&\vdots &&&\ddots& \\ 1 & 0 & \cdots & 0& 0 & \cdots & -1 \\ 0 & 1 & \cdots &0 &-1 & \cdots & 0\\ &&\vdots &&&\ddots& \\ 0 & 1 & \cdots &0 &0 & \cdots & -1\\ \cdots &\cdots &\cdots &\cdots &\cdots &\cdots &\cdots \\ 0 & 0 & \cdots &1&-1&\cdots& 0\\&&\vdots &&&\ddots& \\ 0&0&\cdots&1&0&\cdots& -1\\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} d_{1} \\ d_{2} \\ \vdots \\ d_{r} \\ d^{1} \\ \vdots \\ d^{s} \end{array} } \right]$$
(8)

where \(B_{r}^{s}\) denote the FPA between station \(r\) and satellite \(s\) taken from the related float ambiguity, and \(d_{r}\) and \(d^{s}\) on the right side express the UPDs of the receiver and satellite side, respectively. Since the receiver- and the satellite-side UPDs are linearly related in the equation, the matrix is rank deficient and the rank deficiency is 1. Therefore, a datum should be imposed, such as fixing one UPD as a reference or adding the constraint that the sum of all the satellite UPDs is zero so that the parameters can be solved properly (Li et al. 2016). This approach can be applied to estimate both WL and NL UPDs and in both post-processing and real-time application.

Usually, on the user-side, after applying the satellite UPD products to WL and NL ambiguities, the integer nature of IF ambiguity can be recovered in positioning. Then, the inter-satellite SD ambiguities are aimed to be fixed. It should be noted that the undifferenced or SD ambiguities are estimable if and only if they are in double-differenced (DD) form and that PPP single-receiver ambiguity resolution should be seen as a special case of network DD ambiguity resolution (Khodabandeh and Teunissen 2019). Denoting the reference satellite as \(s_{0}\), the SD float ambiguity can be expressed as,

$$ \Delta \overline{N}_{r}^{{s,s_{0} }} = \overline{N}_{r}^{s} - \overline{N}_{r}^{{s_{0} }} = \Delta N_{r}^{{s,s_{0} }} - \Delta d^{{s,s_{0} }} $$
(9)

where \(\Delta d^{{s,s_{0} }}\) and \(\Delta N_{r}^{{s,s_{0} }}\) refer to SD UPD and integer part between satellites, and \(\Delta d^{{s,s_{0} }}\) is provided in UPD products. After both the WL and NL ambiguities are fixed, the integer IF ambiguity can be introduced into positioning to obtain the PPP-AR solution.

2.3 UPD clustering

The pseudoranges of a single receiver could be biased differently for each satellite, i.e., the satellite-specified pseudorange biases, and the WL ambiguities could be biased accordingly. It is obvious that they will be canceled by forming DD if the satellite-specified biases are the same for all receivers, i.e., they can be absorbed by satellite- and/or receiver-side UPDs. Therefore, only the inconsistency of satellite-specified pseudorange biases between different receiver types will bias WL ambiguities differently and thus contaminate UPD estimates. These inconsistent biases can be revealed in the DD FPAs or UPDs where common impacts at both receiver and satellite side are already removed (Hauschild and Montenbruck 2016) as follows:

$$ d_{r,q}^{s,p} = d_{r}^{s,p} - d_{q}^{s,p} = \left( {d_{r}^{s} - d_{r}^{p} } \right) - \left( {d_{q}^{s} - d_{q}^{p} } \right) $$
(10)

where \(d_{r,q}^{s,p}\) refers to the DD of UPDs between satellite \(s\) and \(p\) and receivers \(r\) and \(q\). The DD UPDs should be very close to zero if the UPDs estimated from different receivers are consistent with each other. Under the normal distribution assumption for UPD uncertainty, a significant inconsistency can be detected as outliers under the 3σ principle.

As soon as a large set of UPDs samplings is available, we can classify them into groups according to their agreement. This group is also called as a cluster, a region in which the density of samples is locally higher than in other regions (Krishna and Narasimha 1999). Significant differences among different groups would indicate that they have inconsistent UPDs.

There are a number of mathematical methods for automatic classification or clustering, for example, the K-Means algorithm. The K-Means algorithm aims to choose the best clusters that minimize the sum square error of the distances between a sample and its cluster centroid. The problem can be expressed mathematically as (Aristidis et al. 2003; Hansen et al. 2005),

$$ P\left( {S\left\{ {S_{1} , \cdots ,S_{k} } \right\}} \right) = \min \sum\limits_{{i = 1}}^{k} {\sum\limits_{{j = 1}}^{{n_{{S_{i} }} }} {\left\| {w_{{{\text{ij}}}} - \bar{w}_{i} } \right\|} } $$
(11)

where the \(S\) is the set of total samples to be classified, \(\left\{ {S_{1} , \cdots ,S_{k} } \right\}\) are the k found which meet minimization criteria. For each cluster \(S_{i}\), there are \(n_{{S_{i} }}\) samples, \(w_{{{\text{ij}}}}\) is the \(j\) th sample of the \(i\) th cluster, and \(\overline{w}_{i}\) is the centroid of the \(i\) th cluster which is defined as,

$$ \overline{w}_{i} = \frac{1}{{n_{{S_{i} }} }}\mathop \sum \limits_{j = 1}^{{n_{{S_{i} }} }} w_{{{\text{ij}}}} ,\quad i = 1, \ldots ,k $$
(12)

The algorithm of clustering consists of looping of the following steps.

  • Step 1 Select the data as initial partition of the whole data set (Peña et al. 1999) for the groups \(\left\{ {S_{1} , \cdots ,S_{k} } \right\}\). In this study, the clustering is applied on a satellite base. We set the initial groups according to the receiver types and sort the corresponding FPAs into the groups.

  • Step 2 For each group, calculate its centroids \(\overline{w}_{i}\) from all members \(w_{i,j}\) by Eq. (12) and the standard deviation (STD) of the residuals; take three-time averaged group STD as the threshold for testing the group significance. Two groups are merged if their distance is shorter than the threshold. This will be repeated until no further groups can be merged.

  • Step 3 For each member \(w_{i,u}\) in group \(S_{u}\) reassign it to the closest group \(S_{v}\) if \(w_{i,u} - \overline{{w_{v} }} \le w_{i,u} - \overline{{w_{t} }}\) for all \(t = 1, \ldots ,k\).

  • Step 4 Recalculate centroids for all the clusters after all the member in \(S_{u}\) are checked and then repeat Step 3 for all the groups

  • Step 5 If there is no new group reassignment in Step 3 and 4, the processing stops. Otherwise, repeat Step 2 to 5.

As shown in the above-mentioned steps, the K-Means algorithm has provided an initial partition of the deviation, and the average value of these initial groups is calculated. Then, the instances of the database are relocated to the group represented by the nearest centroid in an attempt to reduce the sum square error. If an instance in the relocation step changes its group membership, then the centroids of the groups and the sum square error should be recomputed. This process is repeated until convergence, that is to say, until the sum square error cannot be further reduced which means no instance changes its group membership.

3 Experimental validation

In the following, we introduce the data and network information and data processing strategy used in the experimental validation.

3.1 Data set

A collection of 1560 regionally distributed reference stations from the National Geodetic Survey (NGS) of USA permanent network during a period of 30-day in January, 2019 was used in this experimental validation. To ensure the stability and precision of UPD estimation, we only estimate UPD for the receivers type which is employed by over 40 stations (Wang et al. 2018), meaning that if one type of receiver has less than 40 stations then it is not used for the UPD estimation. Among them, the majority receivers are from Trimble with 7 different types, accounting for about three-quarters of the total. The second one is Leica installed at about 200 stations and then Septentrio at 72 stations. The distributions of stations are shown in Figs. 1 and 2, with the different colors representing different receiver hardware and firmware versions and the related number of stations are summarized in Table 1.

Fig. 1
figure 1

Distribution of GNSS stations of the USA NGS GNSS network (as of January 2019). Red circles denote the RCV_SEP, green circles denote the RCV_LCA3, and violet circles denote the RCV_LCA6

Full size image
Fig. 2
figure 2

Distribution of different types of Trimble receivers in the USA NGS network (as of January 2019). Red circles denote the RCV_TRMS_32, green circles denote the RCV_TRMS_30, violet circles denote the RCV_TRMS_25, blue circles denote the RCV_TRM9_33, dark slate blue circles denote the RCV_TRM9_30, teal circles denote the RCV_TRM9_25, and dimgray circles denote the RCV_TRM5

Full size image
Table 1 Receiver types and the corresponding number used in this study
Full size table

As shown in Figs. 1 and 2, the numbers of RCV_TRM9 and RCV_TRMS receivers are the largest among all the types. Each type of receiver is enough to estimate UPD on the American region except the RCV_LCA6, which only covers the north-east region. During the test period, all GPS satellites are in good condition except the satellite G04, which is unusable for most days, therefore is excluded in the processing.

3.2 Data processing

For all the stations, PPP is performed to estimate IF ambiguities and WL ambiguities for each station on a daily base. The details of the processing strategies are listed in Table 2.

Table 2 Processing strategies for UPD estimation of the service-side and PPP-AR on the user-side
Full size table

For the UPD estimation, the station-wise static PPP processing is performed, where the station coordinate is estimated as daily constant or tightly constrained to the precise a priori value. The WL ambiguities are calculated based on the MW combinations simultaneously and both WL and IF ambiguities are output for the estimation of WL and then NL UPD as described above and by Ge et al. (2008) and Li et al. (2016). For the user-side positioning, the kinematic mode can be applied and integer ambiguity fixing will be implemented using the estimated UPDs where WL ambiguities are usually fixed by rounding to the closest integer while NL ambiguities fixed using the LAMBDA method (Teunissen 1995).

In this study, we analyze the deviation in the fractional part of daily WL ambiguities from all the stations and the deviations in the UPDs estimated for each type of receivers. Furthermore, their impact on the positioning will also be evaluated.

4 Analysis of WL deviations

After the statistical analysis of the observations used, the temporal stability of the WL FPAs of all stations of about 30 daily solutions is first demonstrated. Then, the averaged FPAs are used in the investigation on consistency among different receivers and in the classification as well. Moreover, the receiver-type-dependent inconsistency of the estimated satellite WL UPDs is compared.

4.1 Observations statistics

The number of ambiguities contributing to the UPD estimation of a single satellite, which is also the number of all continuously observed data arcs of the satellite, is a measure of observation redundancy. As there might be outliers that are excluded in the UPD estimation, the data usage rate, i.e., the percentage of the number of ambiguities used, must also be considered for quality control. For a direct illustration, the statistics of the number of ambiguities and the usage rate in UPD estimation over the test periods (30 days) are shown in Fig. 3 for RCV_TRM9_22 and RCV_TRM9_30, respectively.

Fig. 3
figure 3

The statistics of the observation number (a) and the usage rate (b) over the 30 days test period for receiver type RCV_TRM9_22 (left) and RCV_TRM9_30 (right)

Full size image

In Fig. 3, the boxplot is used to show the statistics where the box contains a 50% confidence interval, while the top and bottom lines contain a 97% confidence interval; the line in the box is the median value. As shown in Fig. 3, for all the satellites the number of data arcs is more than 200, and the data usage rate is above 95%, which means they have a large number of input observations for parameter estimation. For the G11, G18, and G23 satellites, the observation number fluctuation is large than others because some of the poor quality data are eliminated during the data preprocessing.

4.2 Analysis of WL FPAs

As the WL float ambiguities are determined by the observations without any modeling, they are therefore highly correlated with the receiver's hardware characteristics. Each ambiguity or its FPA can be considered as an observation of the corresponding receiver and satellite UPDs, so FPAs can be applied to investigate the UPD consistency among receivers and receiver types as well.

For a station-satellite pair, there may be several ambiguities in a daily solution; their FPAs should agree with each other very well if the continuous observations are long enough (at least 30 min). In principle, the daily FPAs are very stable as reported (Li et al. 2017, 2018). We show here three typical examples of FPA with the good, moderate, and poor agreement in Fig. 4.

Fig. 4
figure 4

Daily FPA estimates for all GPS satellites at three stations: a ACSX with good stability, b CACC with moderate stability, and c P192 with poor stability. The x-axis is the DOY in 2019

Full size image

As shown in Fig. 4c, at the station P192 with a poor stability, there are few satellites with large fluctuation in the daily FPAs, e.g., G09, G13. The other satellites, however, all have stable performance during the 30-day period. All the satellites at the stations with the good (ACSX) and moderate (CACC) agreement have very stable variation, with the biases smaller than 0.05 cycles. Therefore, the averaged FPAs over 30 days can be applied for the investigation of receiver-type-dependent UPD biases.

For each station-satellite pair, there is only one averaged FPA over the test period, and for each satellite, the maximum number is the number of stations that have enough observations to this satellite. For each satellite, the distribution of the FPAs gives a measure of their agreement or disagreement among receivers. We count the FPA number for each type of receivers, i.e., distribution of each receiver type, and put these of all receiver types together to depict their agreement within a receiver type and among different receiver types. Figure 5 shows three typical cases with almost no, small, and significant differences among different receiver types. It should be noted that the G01 satellite is selected as the reference to remove the receiver UPDs which could vary in each daily solution.

Fig. 5
figure 5

Distribution of satellite WL FPAs of receiver types with almost no (a), small (b), and significant (c) differences among different receiver types. The histogram of each receiver type is in a different color. The centralization of each histogram represents the agreement among the receivers and the distance between the histogram of different colors shows bias among receiver types

Full size image

The histogram of each receiver type is in a different color. The centralization of each histogram represents the agreement among the receivers and the distance between the histogram of different colors shows bias among receiver types. From Fig. 5a, for G10 satellite all the FPAs agree with each other within ± 0.05 cycle, i.e., the various types of receivers show high consistency. While for G09 satellite, RCV_TRMS_33 and RCV_TRM9_22 have the deviation of −0.21 cycles from the others, the G32 satellite shows a rather large difference in five groups.

In general, the FPA differences among different receiver types are very significant and vary from satellite to satellite. Therefore, the automatic classification algorithm K-Means described in Sect. 2.3 is applied to obtain the estimates of the differences in a rather convenient way.

The wavelength of the GPS WL ambiguity is 86 cm, and the deviation effect of 0.05 cycles is less than 5 cm. We set the threshold as 0.05 cycles in the classification algorithm, i.e., when the three-time STD is less than 0.05 cycles, this threshold is used to avoid potential incorrect outlier detection. The result is shown in Table 3 with the classified groups and their FPA differences referring to the receiver type of RCV_TRMS_32 (f).

Table 3 Result of K-Means automatic classification of fractional parts of WL ambiguities of receiver types
Full size table

Table 3 shows that most of the receiver types have very good agreement with each other with few exceptions. For the G07, G10, and G25 satellites, all the receiver types have the most consistent characteristics. On the other hand, for the G03, G05, G12, G15, and G32 satellites, there are some receiver types with a large deviation with respect to the reference type RCV_TRMS_32(f). We have checked the observations and estimation procedure for those satellites and have not found anything unusual compared to other satellites. Moreover, the results of the classification also show that RCV_TRMS_30 (e) and RCV_TRMS_32 (f) are classified in the same group for all satellites, indicating that these two types of receivers have very similar performance in terms of WL FPAs. However, RCV_TRM9_22 (a) and RCV_TRM9_30 (b) have the largest difference compared to the reference: RCV_TRMS_32 (f).

The deviations of the WL FPAs among different types of receivers will have negative impact on not only the estimation of satellite UPDs but also PPP ambiguity fixing at the user-side. When a mixture of different receiver types is involved to provide UPD service, the satellite UPDs could be contaminated. More importantly, PPP ambiguity fixing will be degraded if the UPDs are estimated from different types of receivers with significant UPD deviations.

4.3 Satellite WL UPD deviations

We select the RCV_TRMS_32 as reference type because it is the reference type in FPA clustering and compare its UPDs estimates with the other six types of Trimble receivers. The daily UPDs of each receiver type are estimated, i.e., 30 days of UPD solutions in total. The characteristics of the satellite UPD differences of the six receiver types with respect to the reference are shown in the boxplot in Fig. 6 with the average and the STD.

Fig. 6
figure 6

WL UPD differences of different receiver types with respect to the RCV_TERMS_32 type, a RCV_TRM9_22, b RCV_TRMS_30 c RCV_TRM9_33 d RCV_TRMS_25 e RCV_TRMS_30 f RCV_TRM5. For each satellite, the 50% and 95% are within the box and upper and bottom lines, respectively. The red markers show those with a bias larger than 0.1 cycles

Full size image

Figure 6 shows that RCV_TRMS_30 has the most agreement with the reference receiver RCV_TRMS_32, with a bias of ± 0.03 cycles for all satellites. On the other hand, some satellites in the other three types of RCV_TRM9 series and as well as RCV_TRM5 show a large UPD deviation. For the RCV_TRM9, the GPS Block IIF satellites (G03, G06, G08, G09, G25), and Block IIR satellites (G11, G12, G13, G14, G19, and G23) have deviations larger than 0.1 cycles; while for RCV_TRM5, the Block IIF (G06, G08, G09) and Block IIR (G11, G19) satellites show large bias up to 0.2 cycles. These indicate the existence of large biases among different versions of receivers from the same brand. For RCV_TRMS_25 UPDs, the G11, G13, G19, and G29 of IIR type satellites are close to 0.1 cycles. The RCV_TRM9_30 has the worst performance, which has eleven satellites with STD larger than 0.1 cycles, and among them, the STDs of G03, G06, G08, and G09 satellites are even larger than 0.2 cycles. Nonetheless, the WL UPDs deviations, between RCV_TRMS_32 and each Trimble type, including G01, G05, G07, G10, G17, G18, G28, and G31 satellites, remain within ± 0.1 cycles.

A similar phenomenon is also reported in the investigation of FPA characteristics in Sect. 4.2. For example, RCV_TRMS_30 and RCV_TRMS_32 agree very well and the FPA deviation of G09 satellite reaches 0.25 cycles for the RCV_TRM9 types. The UPD differences of the receiver types from other manufacturers are also calculated in the same way with respect to RCV_TRMS_32, as illustrated in Fig. 7. RCV_ALL means UPD estimated from all the receivers despite the receiver-type-dependent biases.

Fig. 7
figure 7

UPD deviation of four no-Trimble receiver types with respect to RCV_TRMS_32. a RCV_SEP b RCV_LCA3 c RCV_LCA6 d RCV_ALL

Full size image

As we can see from Fig. 7, the three receivers show significant deviation with respect to the reference receiver: RCV_TRMS_32. Compared with the mixture of all receivers, i.e., RCV_ALL, the G09 satellite has a slightly larger deviation over 0.1 cycles; while for other satellites, the deviation is within 0.1 cycles. The major reason is that most of the receivers are from the reference type RCV_TRMS_32 and RCV_TRMS_30, and most satellites do not have a receiver-type-dependent bias. The G02, G03, and G19 satellites of RCV_SEP also have significant deviation. However, for RCV_LCA, these two types (RCV_LCA3 and RCV_LCA6) show a similar distribution. The satellite UPD deviations show the same features with FPAs distribution similarity with Table 3. In the type of RCV_LCA6, only G29 and G30 satellites have the deviations of about 0.1 and −0.1 cycles in UPDs and the same approximate 0.08 and −0.09 cycles in FPAs, respectively. Besides, the three types of RCV_TRM9 receivers all have similar large deviations in FPAs and UPDs. It also means that the deviations in FPAs could affect the UPD estimation.

4.4 Deviation correction

To uniformly compensate the WL UPD deviations among different types of receivers for both UPD estimation and for user-side PPP-AR, the satellite-specified WL UPD deviations can be estimated by the classification algorithm or by simply taking the average of the 30 daily UPD estimates depicted in Fig. 8.

Fig. 8
figure 8

Receiver-type-dependent WL UPD biases with respect to RCV_TRMS_32 by taking an average of 30 daily UPD estimates

Full size image

RCV_TRM9_30 has the largest deviation in the pseudorange of ± 0.26 cycles, while RCV_TRMS has the smallest deviation in the pseudorange of ± 0.04 cycles. Besides, for satellite G03, G06, G08, G09, and G19, the UPD deviations among different receiver types could reach 0.2 cycles.

Since it is hard to find out the mechanism of such differences and to precisely characterize receiver WL UPD deviation by a model, therefore, the deviations among different receiver types are calculated and compensated to avoid the trouble caused by inconsistent UPDs at the service and user side. In this study, we calibrate the WL deviations using the model proposed in this paper that classifies according to FPAs. To verify the effectiveness of the calibration after classification, the UPD residuals are calculated for each receiver type, and that of RCV_TRM9_22 and RCV_SEP types are shown as examples in Fig. 9.

Fig. 9
figure 9

WL UPD residuals after the calibration of receiver-type-dependent biases

Full size image

Compared to the subplots in Fig. 9 with the corresponding subplot for RCV_TRM9_22 in Fig. 6 and RCV_SEP in Fig. 7, we can see clearly that the deviations after compensation are reduced to within ± 0.03 cycles, very close to zero and its STDs have similar fluctuations as before. Besides, comparing the box diagram in Fig. 9 with Figs. 6 and 7, we can find that the 30-day residuals of all satellites are reduced and kept very close to 0 without deviations.

After applying the calibration in Table 3, the WL UPDs of all the receiver types can agree within ± 0.03 cycles as depicted in Fig. 10. It is obvious that the remaining receiver-type-dependent WL UPD bias can be ignored for WL ambiguity fixing. As the corrections are derived with respect to a common receiver type, the corrected UPDs, therefore, refer to a unique reference no matter whether the reference receiver type is involved in the processing or not. Of course, corresponding corrections must be applied to the user-side to use the corrected WL UPDs.

Fig. 10
figure 10

The WL UPD deviations after being calibrated

Full size image

5 PPP-AR result

To validate the reliability and accuracy of the corrected WL UPDs, the performance of PPP-AR using calibrated and uncalibrated UPDs is compared in terms of the ambiguity fixing rate, the time to first fix (TTFF), and the positioning accuracy. The TTFF indicates the time when the ambiguity is first successfully fixed continuously from the positioning (Feng and Wang 2008). The processing strategies of PPP-AR are listed in Table 2. Furthermore, the satellite with the highest elevation angle and continuous observations of more than 5 min would be selected as the reference satellite to form the SD ambiguities. The station coordinate from average of PPP daily static solution is used as a reference.

All 10 types of UPDs are utilized for kinematic PPP-AR for the test stations to statistically assess the impact of WL UPD biases on the PPP-AR performance. The ambiguity fixing rate, which is defined as the percentage of the number of epochs with fixed ambiguities, and positioning accuracy, are calculated and shown in Table 4 for PPP-AR using receiver types without UPD differences, or with rather small UPD differences for a few satellites, and in Table 5 with large UPD deviations for a few satellites.

Table 4 PPP-AR result using receiver types without or with small WL UPD deviation
Full size table
Table 5 PPP-AR result using receiver types with large deviation for a few satellites
Full size table

As shown in Table 4, the results of those four types of receivers do not show any significant difference, between the solutions using the WL UPD from the same receiver types and the others, because of the rather small biases in UPDs as discussed above. However, in Table 5, the difference between the solution using UPD from the same receiver type and that mixed is relatively obvious, for example, the fixing rate for the former one is about 97% and about 92% for the latter one, similar for the positioning accuracy. This is clearly the effect of the deviations in WL UPDs for different satellites which leads to poor ambiguity fixing.

Although the majority of stations have similar positioning performance, as we can see from Table 4 and even Table 5, PPP uses the UPD estimated from the same receiver types or the types with small deviation, or for fewer satellites with slightly large deviations, the result of positioning does not show any significant difference. However, there are still some test stations that suffer from significant degradation in positioning performance because of using mismatched UPDs. A set of stations equipped with RCV_TRMS_32 receiver is chosen to compare and analyze the performance using the other five types of UPDs with large deviations and as well as their calibrated ones. Among the results of five types of UPD, the UPD named calibrated indicates that the calibration is used and the results of RCV_TRMS_32 are used as the reference uniformly. The ambiguity fixing rate in this paper is defined as the ratio of the number of fixed epochs to the number of total epochs (Feng and Wang 2008; Li et al 2017). This fixing rate is a commonly used and easy-to-implement index to indicate the availability of the AR solution. We count the number of fixing epochs and do not identify whether the ambiguities are correctly fixed or not because it is difficult to check the correctness of the fixed ambiguity, especially in real-time positioning (Feng and Jokinen 2017). The statistics are summarized in Table 6.

Table 6 Statistics of the PPP-AR results using UPDs with large deviations
Full size table

In Table 6, using RCV_TRM9_22 UPDs for positioning, there are significant low fixing rates and precision divergence, but these positioning results show normal after calibration on WL UPDs. Using the calibrated RCV_TRM9_22 UPD, the accuracy of N, E, and U is improved from 4.9, 7.9, and 10.1 cm to 1.5, 1.6, and 4.0 cm, while those of the RCV_SEP are improved from 3.8, 5.2, and 8.7 cm to 1.4, 1.7, and 3.5 cm, respectively. In these two receiver types, more than 13 satellites were calibrated with an average deviation of about 0.2 cycles. The accuracy and fixing percentage of RCV_TRM9_30, RCV_TRM9_33, and RCV_TRM5 is also improved by about 50%. It is indicated that the deviation in WL UPD makes the ambiguity fixing abnormal, and affects the convergence time and fixing rate and finally results into large positioning errors. However, after calibrating the WL UPD deviations, the fixing rate, TTFF and the accuracy of the PPP-AR have been improved to the same level of that using the UPDs of the same receiver type.

Table 6 shows two stations with obvious impact as examples for discussion, namely CSST and BSRY both equipped with the RCV_TRMS_32 receiver. The positioning differences with respect to the ground true on DOY 6 2019 are shown in Fig. 11, where three PPP-AR solutions using UPDs corresponding to RCV_TRM9_30, RCV_TRM9_30 calibrated, and RCV_TRMS_32, respectively.

Fig. 11
figure 11

Positioning biases of PPP-AR solutions at stations CSST (left) and BSRY (right), from top to bottom is N, E, U, respectively. Both stations equipped with RCV_TRMS_32 receiver using UPDs of RCV_TRM9_30 (red) and RCV_TRM9_30 calibrated (blue) and, RCV_TRMS_32 as reference (green)

Full size image

From Fig. 11, the differences of the PPP-AR solution with RCV_TRM9_30 UPD are rather large up to 10 cm and 20 cm in horizontal for CSST and BSRY, respectively, and even larger in vertical. To investigate the reason, Fig. 12 shows the number of fixed ambiguities at each epoch for both stations. It can be seen that for the RCV_TRMS_32 receivers, the fixing rate using RCV_TRM9_30 UPDs is significantly lower than using the UPD from the same receiver type and the calibrated one. After 2500 epochs, the number of fixed satellites using RCV_TRM9_30 UPD at the BSRY station increased to 7, indicating that the ambiguity fix was again successfully performed. However, it can be seen from Fig. 11 that there is a significant deviation of almost 15 cm in the positioning bias after 2500 epochs, indicating an incorrect fixing at this stage. The accuracy in terms of the RMS of positioning differences in N, E, and U is improved from 2.5, 4.1, 5.8 cm to 1.5, 1.3, 4.3 cm, and from 11.2, 4.2, 12.2 cm to 1.5, 1.7, 4.8 cm for CSST, and BRSY stations, respectively, with the calibrated UPDs.

Fig. 12
figure 12

The number of fixed satellites for CSST (left) and BSRY (right)

Full size image

Furthermore, the WL and NL UPDs differences between RCV_TRM9_22 and RCV_TRMS_32 are plotted in Fig. 13. There are four satellites with deviations larger than 0.2 cycles and seven satellites larger than 0.1 cycles in WL UPDs, whereas the NL UPDs also exist the fluctuations. This might be caused by using the mismatch WL UPDs, which results in degraded WL ambiguity fixing and consequently different NL ambiguities involved in the NL UPD estimation. Meanwhile, Fig. 14 gives the differences of carrier and pseudorange residuals of PPP-AR using RCV_TRM9_22 and RCV_TRMS_32 UPDs on CSST and BSRY station, respectively.

Fig. 13
figure 13

The WL (left) and NL (right) UPDs difference between RCV_TRM9_22 and RCV_TRMS_32

Full size image
Fig. 14
figure 14

Differences of pseudorange (left) and phase (right) residuals between PPP-AR solutions using RCV_TRM9_22 and RCV_TRMS_32 UPD at station CSST (top) and BSRY (bottom)

Full size image

It can be seen from the residual differences that both the carrier and the pseudorange exhibit diverging characteristics when the number of fixed satellites is decreased. Around 500 and 1500 epochs, at 4:10 and 12:30, the position differences of the solution with uncalibrated UPDs begin to diverge, at CSST and BRSY station, respectively. However, with the calibrated UPDs the convergence and stability of the solution can be retained in this stage. Therefore, with UPDs with a large deviation, the positioning performance could be abnormal, and the effectively calibrated UPDs can significantly improve the positioning accuracy and fixing rate.

6 Conclusion and remarks

In this study, we present a data processing strategy to investigate the inconsistency of WL UPDs estimated from different receiver types caused by the related pseudorange biases and an approach to estimate and calibrate the UPD differences to improve the WL ambiguity fixing in PPP.

An experimental validation has been carried out with about 1560 stations of 10 receiver types from the NGS GNSS permanent network to evaluate the characteristics of the WL UPD deviations among different receiver types and their impact on WL UPD estimation and consequently the effect on PPP ambiguity solution.

The experiment demonstrates that WL UPDs estimated from different receiver types could be different remarkably; for example, there are large deviations up to 0.3 cycles between the same Trimble receiver with different firmware versions. It is confirmed that such UPD biases would affect satellite UPD estimation in case of receivers of mixed types are involved and also lead to possible wrong PPP-AR if inconsistent UPDs are implemented.

Fortunately, the deviations are very stable and therefore can be estimated very precisely and applied as calibration, so that the calibrated UPDs can have the same accuracy as the UPDs estimated from the same receiver type. With the calibrated UPDs, the PPP-AR accuracy can be improved by 50%, and the TTFF, as well as the fixing rate, also has been improved.