1 Introduction

Pseudorange measurements of the Global Navigation Satellite Systems (GNSSs) suffer from the delays in signal generation, transmission and reception chains in both satellite and receiver sides (Montenbruck and Hauschild 2013). Pseudorange or code delays strongly depend on the signal frequency and signal modulation type (Sleewagen and Clemente 2018). Different bias characteristics manifest for signals from different GNSSs using distinct signal modulations, even if they are modulated on the same carrier frequencies (e.g., GPS and Galileo L1/E1 and L5/E5a signals). As receiver manufacturers might employ different options for the optimal tracking of GNSS signals, the variation of pseudorange biases also shows dependency on the receiver model. Depending on the receiver front-end characteristic and the correlator design, different biases may be present for different receiver hardware settings, resulting in the inconsistent pseudorange biases for a GNSS network with mixed types of receivers (Hauschild and Montenbruck 2016). The inter-receiver comparisons have demonstrated that the variation of such pseudorange biases can be significantly reinforced in case of the activation of multipath mitigation settings (Hauschild and Montenbruck 2016). Even though the receivers are configured with their standard options, the internal thermal stabilization of hardware may also cause short-term bias variations (Zhang and Teunissen 2015), or the immediate bias changes in case of the re-starts of the receivers (Wanninger 2011). In addition to the aforementioned frequency-dependent pseudorange delays, GNSS pseudorange measurements are also subject to the frequency- and elevation-dependent group delay variations (GDV) at the transmission and reception part of the antennas (Wanninger et al. 2017). The problem of pseudorange biases remains a challenge to the estimation of bias-free ionospheric total electron contents (TECs) from GNSS data (Li et al. 2015; Wang et al. 2016a, 2019a), Melbourne–Wübbena (MW) linear combination for widelane ambiguity resolution (Loyer et al. 2012), and other pseudorange-based positioning and timing applications (Montenbruck and Hauschild 2013).

Consideration and proper handling of pseudorange biases are of crucial importance in GNSS ionospheric analysis. The estimation of slant TECs (STECs) from the geometry-free combination of dual-frequency GNSS data needs the calibration of differential code biases (DCBs) in satellite and receiver parts to derive the bias-free STEC estimates (Sardón and Zarraoa 1997). Aside from the timing or broadcast group delays (TGD/BGD) provided by the GNSS providers themselves (Wilson et al. 1999; Wang et al. 2019b), a small set of DCB products have been defined and generated using the legacy GPS and GLONASS signals, with the use of different global ionospheric mapping techniques by individual Ionosphere Analysis Centers (IACs, Hernández-Pajares et al. 2009) of the International GNSS Service (IGS, Dow et al. 2009). The modernization of GPS and the construction of Galileo and BeiDou lead to a multi-constellation and multi-frequency GNSS environment. It becomes difficult to maintain the multi-GNSS DCBs covering all possible signal pairs. As an example, diverse signal types are presently supported by the receiver independent exchange (RINEX, IGS RINEX WG and RTCM-SC104 2018) format v3.04 for each individual GNSS constellation, which means the possibility of forming a large number of signal pairs and associated DCBs. Depending on the availability of observation data from a given tracking station network, different sets of DCBs might be estimated and maintained by different ACs, e.g., the multi-GNSS DCB products generated by German Aerospace Center (DLR, Montenbruck et al. 2014a) and Chinese Academy of Sciences (CAS, Wang et al. 2016b) from the multi-GNSS experiment (MGEX, Montenbruck et al. 2017) network of the IGS. To derive the required DCBs, which are not generated by the providers, a linear combination of different sets of DCBs can be formed with the use of a computationally lean algorithm as reported in Montenbruck and Hauschild (2013). Ideally, the reconstruction of DCBs from the chaining of multiple DCBs needs the definition of an independent set of DCBs covering all relevant signals of each satellite navigation system. Aside from the difficulties in establishing a non-redundant set of DCBs, the reconstructed DCBs might differ from the direct DCB estimates since different types of receivers might be involved in the computation of associated DCBs.

Another solution for code bias handling is to parametrize the biases of individual observations in their undifferenced forms. As pseudorange absolute code biases are not directly estimable, specific datums need to be established for the transition from differential to pseudo-absolute code biases (Villiger et al. 2019). Pseudo-absolute code biases have been defined in the standard of Radio Technical Commission for Maritime Service (RTCM-SC 2016) and formatted in RTCM v3 State Space Representations (SSRs). Some real-time ACs of the IGS, like the Centre National d’Etudes Spatiales (CNES, Laurichesse 2011), have distributed their own estimated code bias corrections through the Network Transportation of RTCM by Internet Protocol (NTRIP). In addition to differential signal biases (DSBs), the observable-specific signal biases (OSBs) are also defined in the new bias SINEX format 1.00 (Schaer 2016). The OSB parameterization is more straight forward, since the derived pseudo-absolute code biases can be directly corrected in the undifferenced observation equations. Multi-GNSS OSB products have been generated at the Center for Orbit Determination in Europe (CODE, Dach et al. 2009) by the combination of their clock and ionosphere analysis since mid-2016 (Villiger et al. 2019), and a first analysis of GPS, GLONASS, Galileo and BeiDou-2 satellite OSBs covering an 18-month period was presented in the same literature.

The purpose of the paper is to report a method for GPS and GLONASS OSB computation as part of the local ionospheric analysis and examine the impact of different networks (IGS vs. MGEX) and receiver groups on the estimation results. The proposed method can be easily extended for the computation of Galileo, BeiDou and QZSS OSBs as well (Wang et al. 2018). The OSB parameterization and estimation in ionosphere analysis is first presented, followed by the description of an automatic bias (AutoBIA) realignment method to handle the discontinuities in OSB series. With the use of observation data from the IGS and MGEX stations, GPS and GLONASS OSBs from the distinct networks are estimated and analyzed during a 2-year period (2017–2018), and the effects of different groups of receivers on the resulting OSB estimates are also discussed. Summary and conclusions are finally given.

2 Method

Pseudorange absolute code biases are not directly estimable, but only the biases in differenced forms can be estimated in the linear combination of associated observations. In the ionosphere-free (IF) linear combination for clock analysis, the bias parameters are lumped into the clock offsets, which are impossible to separate them from each other. In the geometry-free (GF) linear combination for ionosphere analysis, the satellite-plus-receiver (SPR) GF biases are directly estimable in case the two signals are modulated on the same frequency; otherwise, the ionospheric delay needs to be removed by a priori ionospheric model or jointly estimated with SPR GF biases (Montenbruck et al. 2014a). To change from the differential to observable-specific biases, the establishment of specific conventions is required to solve the rank deficient problem in the parameterization and estimation of OSBs (Villiger et al. 2019). This section starts with the discussion of OSB parameterization in ionosphere analysis, followed by a detailed description of a local ionospheric modeling technique for OSB estimation and an AutoBIA method for OSB realignment, respectively.

2.1 OSB parameterization in ionosphere analysis

Ignoring the inconsistent signal distortions caused by the different settings in receiver’s front-end bandwidth and employed filters (Hauschild and Montenbruck 2016), the practical approach in network-based bias estimation is to separate the pseudorange biases of code division multiple access (CDMA) signals into a receiver-dependent part and a satellite-dependent part (Montenbruck et al. 2014a). Code biases in satellite and receiver sides are commonly assumed to be constant through the data arc (e.g., one day), although the internal thermal stabilization of receiver hardware might cause short-term bias variations in certain receivers (Zhang et al. 2019). As for the bias estimation of GLONASS which uses frequency division multiple access (FDMA) to distinguish signals from different satellites, while the parameterization of satellite-receiver pair biases has been employed by CODE and Jet Propulsion Laboratory (JPL) in their global ionospheric map (GIM) computation (Dach et al. 2009; Vergados et al. 2016), the generation of satellite- and receiver-specific GLONASS biases is still commonly applied in global ionosphere analysis (Roma-Dollase et al. 2017). The estimated ionospheric parameters are affected by the simple separation of GLONASS biases into satellite-independent and receiver-independent parts considering the influences of GLONASS inter-channel biases (Wanninger 2011).

The parameterization of OSBs in ionosphere analysis starts from the observation equation for pseudorange measurement \( P_{ix}^{s} \) of a satellite s with the tracking mode/channel \( x \) on the frequency \( f_{i} \)

$$ P_{ix}^{s} = \rho_{{}}^{s} + \alpha_{i} I_{{}}^{s} + T_{{}}^{s} + c\left( {\delta t - \delta t^{s} } \right) + \underbrace {{b_{ix} + b_{ix}^{s} }}_{{{\text{SPR}}_{ix} }} + \varepsilon_{ix}^{s} $$
(1)

which uses the brief notation \( \rho_{{}}^{s} \) for the geometrical distance between satellite and receiver, \( I_{{}}^{s} \) for the ionospheric slant TEC in TEC unit (TECU) multiplied by the frequency-dependent factor \( \alpha_{i} = {{40.3} \mathord{\left/ {\vphantom {{40.3} {f_{i}^{2} }}} \right. \kern-0pt} {f_{i}^{2} }} \), i for the considered frequency brand, \( T_{{}}^{s} \) for the tropospheric delay, \( c \) for the speed of light, \( \delta t \) and \( \delta t^{s} \) for the receiver and satellite clock offsets, respectively, \( {\text{SPR}}_{ix} \) for the satellite-plus-receiver code bias, which is the sum of absolute code biases in satellite (\( b_{ix}^{s} \)) and receiver (\( b_{ix} \)) parts, and \( \varepsilon_{ix}^{s} \) for the residual error comprising multipath and thermal errors.

Forming differences of observations between tracking modes \( x \) and \( y \) on the two distinct frequencies \( f_{i}^{{}} \) and \( f_{j}^{{}} \) eliminates the satellite-receiver ranges, satellite and receiver clock offset as well as the tropospheric delay (i.e., the geometry part), which leads to

$$ \Delta P_{ij,xy}^{s} = \alpha_{ij} I^{s} + \underbrace {{\overbrace {{b_{ix} - b_{jy} }}^{{{\text{DCB}}_{ij,xy} }} + \overbrace {{b_{ix}^{s} - b_{jy}^{s} }}^{{{\text{DCB}}_{ij,xy}^{s} }}}}_{{{\text{SPR}}_{ij,xy} }} + \Delta \varepsilon_{ij,xy}^{s} $$
(2)

where \( \alpha_{ij} \) equals to \( 40.3\left( {f_{i}^{ - 2} - f_{j}^{ - 2} } \right) \),\( b_{ix} \),\( b_{jy} \),\( b_{ix}^{s} \) and \( b_{jy}^{s} \) are the receiver- and satellite-related code biases,\( \Delta \varepsilon_{ij,xy}^{s} \) denotes the differential multipath and thermal errors, \( {\text{DCB}}_{ij,xy} \) and \( {\text{DCB}}_{ij,xy}^{s} \) are the associated inter-frequency biases in receiver and satellite sides, respectively.

The creation of GF linear combination of different observations on the same frequency also cancels out the ionospheric part

$$ \Delta P_{i,xy}^{s} = \underbrace {{\overbrace {{b_{ix} - b_{iy} }}^{{{\text{DCB}}_{i,xy} }} + \overbrace {{b_{ix}^{s} - b_{iy}^{s} }}^{{{\text{DCB}}_{i,xy}^{s} }}}}_{{{\text{SPR}}_{i,xy} }} + \Delta \varepsilon_{i,xy}^{s} $$
(3)

and, aside from the differential multipath and noise errors, only the code biases in receiver and satellite parts remain. The associated receiver and satellite intra-frequency biases are denoted as \( {\text{DCB}}_{i,xy} \) and \( {\text{DCB}}_{i,xy}^{s} \), respectively.

It is obvious that the intra-frequency SPR biases \( {\text{SPR}}_{i,xy} \) are directly estimable referring to Eq. (3), with a zero-mean assumption of thermal noises and multipath errors over the data arc (Montenbruck et al. 2014a). The inter-frequency SPR biases \( {\text{SPR}}_{ij,xy} \) can be determined following Eq. (2), with the use of a priori ionospheric information or the estimation of local/global ionospheric parameters (Vergados et al. 2016; Wang et al. 2016b). To generate the satellite- and receiver-specific biases, a zero-constellation-mean constraint is set by the IGS convention, assuming a comparable stability of satellite code biases within each individual constellation (Li et al. 2012). In case of OSB estimation, additional conditions need to be employed to eliminate the over parameterization and associated rank deficient problem. As an example, the separation of satellite OSBs \( b_{ix}^{s} \), \( b_{iy}^{s} \) and \( b_{jy}^{s} \) illustrated in Eqs. (2) and (3) needs the definition of one additional constraint between the arbitrary two OSB parameters. To ensure the derived OSBs compatible with the IGS clock product, the constraint between the two signals, which is conventionally employed in the generation of satellite clock offsets, is applied (Villiger et al. 2019).

$$ \frac{{f_{i}^{2} }}{{f_{i}^{2} - f_{j}^{2} }}b_{ix}^{s} - \frac{{f_{j}^{2} }}{{f_{i}^{2} - f_{j}^{2} }}b_{jy}^{s} = 0 $$
(4)

where \( b_{ix}^{s} \) and \( b_{jy}^{s} \) denote the code biases of the two reference signals CIX and CJY, identified by their RINEX v3.04 names. For GPS, the two signals are C1W and C2W, and for GLONASS, which are C1P and C2P according to the IGS clock convention (Montenbruck et al. 2014b). Obviously, an identical constraint is required for the generation of receiver OSBs.

The OSB parameterization in ionosphere analysis is summarized as follows. One bias parameter is first set for the individual pseudorange measurement. The priority list of pseudorange observations is then defined considering the availability of observation data from the involved GNSS network, or simply including all signals supported by the RINEX v3.04 standard. The pre-defined pseudorange priority lists are employed to form the GF linear combinations and to generate the associated SPR differential biases following Eqs. (2) and (3). For global VTEC modeling, the spherical harmonic (SH) expansion is commonly applied under an ionospheric single-layer assumption. To cover all possible signal pairs, the global ionospheric parameters generated from the initial step of the ionospheric analysis can be used as a priori information for the subsequent estimation of SPR differential biases (Villiger et al. 2019). For local VTEC modeling, the generalized trigonometric series (GTS, Yuan and Ou 2004) function can be employed for the joint estimation of local ionospheric parameters and all relevant SPR biases in differenced forms at individual stations (Li et al. 2012; Wang et al. 2016b). The derived SPR differential biases serve as the pseudo-observations for the subsequent OSB estimation (Wang et al. 2018), which are comparable to the bias normal equations (NEQs) without any fixed datum definitions as reported in Villiger et al. (2019).

Aside from the IF constraints (identical to the IGS clock conventions) applied in satellite and receiver parts within individual constellations, zero-mean conditions are also required for each type of satellite biases in the generation of satellite- and receiver-specific OSBs. To determine the satellite-receiver pair OSBs in case of GLONASS, only the aforementioned IF constraints need to be set for each GLONASS satellite. Since the constraints are applied after the computation of all relevant SPR differential biases, it allows the determination of adequate datum definitions after analyzing the contents of the involved bias types. The daily OSB solution can be easily estimated from the one-day pseudo-observations of SPR biases in differenced forms by the least-squares. The combination of SPR differential biases from different days leads to a mean OSB estimate across the covered time period.

2.2 Local ionosphere modeling for OSB estimation

GPS L1, L2 and L5 pseudorange biases in dual-frequency GF linear combinations for OSB estimation are illustrated in Fig. 1. The reference observation is first defined for pseudorange measurements on the individual frequency, i.e., C1W for L1, C2W for L2 and C5Q for L5 frequency, respectively. As the estimation of intra-frequency SPR biases is free of ionospheric modeling errors, GF combinations are formed between each pseudorange measurement and the reference observation within the respective frequency band. For the estimation of inter-frequency SPR biases, the reference frequency, i.e., L1 frequency, is selected. The pre-defined reference observations are then used to form GF combinations between L2/L5 and L1 reference frequency. A similar pre-defined observation priority list can be set for GLONASS and easily extended to other GNSSs like Galileo, BeiDou and QZSS. The resulting SPR differential biases shall cover all relevant code observations of each constellation.

Fig. 1
figure 1

Illustration of GPS pseudorange biases in GF linear combinations for OSB estimation

Full size image

Other than the intra-frequency SPR biases, the estimation of inter-frequency SPR biases requires the proper handling of ionospheric delay errors. Different from the global ionospheric modeling using SH expansions employed by CODE (Villiger et al. 2019), a modified GTS (MGTS) function is presented here for the joint estimation of local ionospheric activities and inter-frequency SPR biases based on an ionospheric thin-layer assumption. The VTEC \( V_{{}}^{s} \) in MGTS is modeled as the sum of a two-dimensional polynomial of the latitude \( \varphi_{d} \) and longitude \( \lambda_{d} \) plus a finite Fourier series expansion depending on the local time t

$$ \left\{ {\begin{array}{*{20}l} {V_{{}}^{s} \left( {\varphi_{d} ,\lambda_{d} ,h} \right) = \sum\limits_{n = 0}^{{n_{\text{max} } }} {\sum\limits_{m = 0}^{{m_{\text{max} } }} {\left\{ {E_{n,m} \cdot \varphi_{d}^{n} \cdot \lambda_{d}^{m} } \right\}} } + \sum\limits_{k = 0}^{{k_{\text{max} } }} {\left\{ {C_{k} \cos \left( {k \cdot h} \right) + S_{k} \sin \left( {k \cdot h} \right)} \right\}} } \hfill \\ {h = {{2\pi \left( {t - 14} \right)} \mathord{\left/ {\vphantom {{2\pi \left( {t - 14} \right)} T}} \right. \kern-0pt} T},\;\;T = 24\;{\text{h}}} \hfill \\ \end{array} } \right. $$
(5)

where n, m and k are the degrees of polynomial and Fourier series expansions, with the corresponding maximum degrees nmax, mmax and kmax, respectively, and En,m, Ck and Sk denote the model coefficients to be estimated at the individual station.

In the original GTS function, the local VTEC is described as a function of the geographic latitude and local time (Yuan and Ou 2004; Wang et al. 2016b), whereas in the proposed MGTS, the spherical cap coordinate system (Haines 1988; Liu et al. 2018) is used instead of the geographic coordinate system. The latitude \( \varphi_{d} \) and longitude \( \lambda_{d} \) in Eq. (5) are computed as follows

$$ \left\{ {\begin{array}{*{20}l} {\varphi_{d} = ds \cdot \cos \varphi_{c} ,\quad \lambda_{d} = ds \cdot \sin \varphi_{c} } \hfill \\ {ds = \left( {R + H} \right) \cdot \lambda_{c} } \hfill \\ {\varphi_{c} = \arccos \left\{ {\sin \varphi_{0} \cdot \sin \varphi + \cos \varphi_{0} \cdot \cos \varphi \cdot \cos \left( {\lambda - \lambda_{0} } \right)} \right\}} \hfill \\ {\lambda_{c} = \arcsin \left\{ {{{\sin \left( {\lambda - \lambda_{0} } \right) \cdot \cos \varphi_{c} } \mathord{\left/ {\vphantom {{\sin \left( {\lambda - \lambda_{0} } \right) \cdot \cos \varphi_{c} } {\sin \varphi_{c} }}} \right. \kern-0pt} {\sin \varphi_{c} }}} \right\}} \hfill \\ \end{array} } \right. $$
(6)

which uses the notations \( (\varphi_{c} ,\lambda_{c} ) \) for the latitude and longitude of the ionospheric pierce point (IPP) in spherical cap coordinate system, \( (\varphi_{0} ,\lambda_{0} ) \) and \( (\varphi_{i} ,\lambda_{i} ) \) for the geographic latitude and longitude of the station and the IPP, respectively, R for the mean radius of the Earth (6378 km) and H for the altitude of the assumed single-layer ionosphere (450 km).

A stochastic model depending on satellite elevation (e) and local time (t), as shown in Eq. (7), is employed in MGTS

$$ \left\{ {\begin{array}{*{20}l} {w\left( {e,t} \right) = w_{s} \left( e \right) \cdot w_{i} \left( t \right)} \hfill \\ {w_{s} \left( e \right) = {1 \mathord{\left/ {\vphantom {1 {\left( {1 + \cos^{2} e} \right)}}} \right. \kern-0pt} {\left( {1 + \cos^{2} e} \right)}}} \hfill \\ {w_{i} \left( t \right) = 1 + \frac{1}{ 2}\left( {\cos \left( {\frac{t - 2}{12}\pi } \right) - 1} \right)} \hfill \\ \end{array} } \right. $$
(7)

With the use of a satellite-elevation and local-time-dependent stochastic model, the effects of observation noises and day-night differences of the ionospheric variability are expected to be reduced.

2.3 OSB realignment method

As a zero-constellation-mean condition is conventionally set to generate the satellite- and receiver-specific OSBs, systematic offsets might exist between OSB estimates of different days in case of changes in satellite constellations. To allow for the direct comparison between different OSB series, the realignment of OSBs to a common reference datum is required. One practice to remove the systematic offset is to align the bias solutions to a common fixed set of satellites across the entire period. However, such a realignment method is not flexible to changes in satellite constellations. In case, a new satellite is incorporated or an old satellite is excluded from the zero-mean condition, the bias series needs to be realigned again, thus not applicable for the routine computation of continuous OSB series. Following Sanz et al. (2017), an AutoBIA realignment method was presented in Wang et al. (2019b) to handle the discontinuities in multi-GNSS DCB series. The idea behind AutoBIA is to align the bias series by automatically removing the bias datum differences with a common set of reference satellites. Here, AutoBIA method is extended for the routine realignment of OSB parameters as well as the generation of continuous OSB series with common bias datums.

For a given group (g) of OSB series (Bg) from one provider, a common set of \( N_{g}^{\text{svn}} \) satellites between the current day and previous Nd days is first determined based on the satellite space vehicle number (SVN). For each individual satellite OSB series, the difference between OSB values of the ith day and mean value of the entire Nd days is first checked. The outlier, which is empirically defined as 0.5 ns, is removed from the OSB series and the number of common satellites is then updated to \( N_{g}^{\text{s}} \). For all involved groups of OSB series from different providers, a common set of \( N_{{}}^{s} \) satellites is finally derived on the basis of the individual \( N_{g}^{s} \;\;\left( {g = 1 \ldots G} \right) \), where the notation G denotes the maximum number of the involved OSB groups. For the considered OSB series Bg, the reference datum offset \( \delta_{g} \left( n \right) \) is calculated as the difference between OSB mean values of the current day (\( \bar{B}_{g,c} \)) and previous Nd days (\( \bar{B}_{g,p} \))

$$ \left\{ {\begin{array}{*{20}l} {\delta_{g} \left( n \right) = \bar{B}_{g,c} \left( n \right) - \bar{B}_{g,p} \left( n \right)} \hfill \\ {\bar{B}_{g,c} \left( n \right) = \frac{1}{{N^{s} }}\sum\limits_{j = 1}^{{N^{s} }} {B_{g}^{j} \left( n \right)} } \hfill \\ {\bar{B}_{g,p} \left( n \right) = \frac{1}{{N_{d} }}\sum\limits_{i = 1}^{{N_{d} }} {\frac{1}{{N^{s} }}\sum\limits_{j = 1}^{{N^{s} }} {B_{g}^{j} \left( {n - i} \right)} } } \hfill \\ \end{array} } \right. $$
(8)

where n denotes an arbitrarily given day, Nd is the selected time span which is fixed as 7 days in our analysis, Ns is the number of reference satellites determined by all of the involved groups of OSB series.

The aligned OSB value of a given satellite j, which refers to the reference datum \( \bar{B}_{g,p} \left( n \right) \), is calculated as follows

$$ \bar{B}_{g}^{j} \left( n \right) = B_{g}^{j} \left( n \right) - \delta_{g} \left( n \right) $$
(9)

and, the comparison of the aligned OSBs between two different groups g (i.e., \( \bar{B}_{g}^{j} \)) and k (i.e., \( \bar{B}_{k}^{j} \)) is easily achieved by a conversion of the reference datum between \( \bar{B}_{g,p} \) and \( \bar{B}_{k,p} \), as they refer to a common set of reference satellites.

The systematic offsets between different sets of OSB series are eliminated using the aforementioned AutoBIA method. The common satellites are determined by their SVNs, and the potential bias jumps or drops are detected and removed using an outlier check procedure. Obviously, changes in satellite constellations (e.g., the launch of new satellites or decommission of old satellites) will not affect the continuity of the entire OSB series, thus suitable for the analysis of bias repeatability within individual OSB series as well as the comparison of OSB products from different providers.

3 Data sets

GPS and GLONASS observation data from the IGS and MGEX networks are used to estimate the respective OSBs and examine the consistency between GPS/GLONASS OSB solutions from different networks and different types of receivers. The analysis was presented during a 24-month period from January 2017 to December 2018. An overview of GPS L1/L2/L5 and GLONASS L1/L2 OSBs generated from the IGS and MGEX stations is summarized in Table 1. As the IGS observation data are commonly provided in RINEX v2.x format, they are first mapped into RINEX v3.04 observation codes. For GPS, the pseudorange C1, P1 and P2 are mapped to C1C, C1W and C2W observation codes, respectively. Considering that diverse observation types are defined on GPS L2 frequency in RINEX v3.04 standard, the legacy C2 is simply converted to C2C observable. The cross-correlation (C–C) GPS receivers, whose C2 observable is actually C1 + (P2 − P1), are not included in our analysis. As for GLONASS, the relation between RINEX v2.x and v3.04 observation codes is C1 for C1C, P1 for C1P, C2 for C2C and P2 for C2P. In case of MGEX observation data, no modern GNSS receivers are currently capable of GPS C2C signal, but only C2W, C2S, C2L and C2X signals are trackable on the L2 frequency. The remaining observation types on GPS L1 frequency and GLONASS L1/L2 frequency bands are identical for the MGEX and IGS receivers. As for OSBs on GPS L5 frequency, only the estimation of C5Q and C5X OSBs is supported by the MGEX receivers.

Table 1 Overview of GPS and GLONASS OSBs generated from the IGS and MGEX networks
Full size table

GPS and GLONASS pseudorange observables on each individual frequency, as shown in Table 1, are given in their priorities. For GPS, C1W and C2W are selected as the reference observables on L1 and L2 frequencies, respectively, and for GLONASS, the reference observables are C1P and C2P. The GF linear combinations are first formed between each involved pseudorange measurement and the reference observable within individual frequencies to estimate the intra-frequency SPR biases. The aforementioned reference observables are then used to form the GF combinations between L2/L5 and L1 frequency, and to generate the inter-frequency SPR biases. Since a limited number of MGEX receivers support the tracking of GPS C1W signal, the C1C signal is selected as an additional reference observation to form the corresponding GF linear combinations between L2/L5 and L1 frequency. The MGTS function is employed for local ionospheric activity modeling in both IGS and MGEX OSB analysis. To generate the satellite- and receiver-specific OSBs from all relevant intra- and inter-frequency SPR biases, the IF constrains and zero-mean conditions are applied to the reference observables in satellite/receiver sides and each individual satellite OSB type, respectively.

As depicted in Fig. 2, the number of IGS stations used for our GPS/GLONASS OSB analysis varies between 400 and 500 during the 2-year test period. The number of MGEX stations is less than that of IGS stations, which changes from about 200 in early 2017 to roughly 280 by the end of 2018. An apparent increase in MGEX stations is noticed since the fourth quarter of 2018. Note that around 80% of the MGEX stations also belong to the IGS network by the end of 2018. The AutoBIA method is applied for the realignment of GPS and GLONASS OSB series, to allow for the direct comparisons between IGS (denoted by CASIGSRAP) and MGEX (denoted by CASMGXRAP) OSB solutions. The comparison of the original and realigned C1W OSBs for one GPS satellite G066 (identified by its SVN) is depicted in Fig. 3, in which the top and bottom plots correspond to the IGS and MGEX solutions, respectively. The variation of the realigned OSBs is less scatted compared with the original ones, but there exists an obvious drop during days 118–161 of the year 2017. To interpret this unexpected discontinuity, the experimental OSB product from CODE is included for comparison. As shown in Fig. 3 (the pink line), the bias deviations in CODE OSB series reach to 0.13 and 0.38 ns at the two epochs, which keep in proper accord with the Notice Advisory to Navstar Users (NANU) messages as explained in Villiger et al. (2019). Moreover, the stability of IGS and MGEX OSB estimates is at a comparable level for the selected SVN G066, with a standard deviation (STD) value of 0.126 and 0.133 ns, respectively, across the entire period.

Fig. 2
figure 2

The number of IGS and MGEX stations used for GPS/GLONASS OSB estimation during 2017–2018

Full size image
Fig. 3
figure 3

The original and realigned C1W OSB series of GPS SVN G066 during a two-year period. Top and bottom plots correspond to the IGS and MGEX solutions, respectively. The pink line denotes GPS C1W OSB estimates from CODE

Full size image

4 Results and discussion

Based on the realigned OSB solutions over the years 2017–2018, the analysis of GPS and GLONASS satellite OSBs from IGS and MGEX networks is first shown in this section, followed by the verification of different groups of receivers on the resulting OSB estimates. Receiver OSB results are not discussed in this study.

4.1 GPS OSB solutions

One Block IIF satellite, SVN G066, is first selected as a representative to check the long-term variability of GPS satellite OSB estimates from the independent IGS and MGEX networks. Time series of satellite OSBs within each calendar month and their standard deviations are depicted in Fig. 4. The C2W shows a large scatter compared to C1C and C1W OSBs. For the IGS solution, the mean STD values across the test period are 0.074, 0.071 and 0.112 ns for C1C, C1W and C2W OSBs, respectively. For the MGEX solution, the values are 0.075 ns for C1C, 0.070 ns for C1W and 0.113 ns for C2W, respectively. The inferior repeatability of C2W OSB can be partly attributed to the IF constraint between C1W and C2W observables given in Eq. (4), which results in the large variation range of C2W OSBs. A steep decrease is found in C1C, C1W and C2W OSB series in May 2017. The discontinuity is also noticed in the experimental OSB series of CODE, which is more likely caused by the in-orbit changes of G066 satellite vehicle (SV) as documented in the NANU messages. A positive deviation on the magnitude of 0.1–0.2 ns exists between MGEX and IGS OSB estimates. Considering a comparable stability of the two OSB products, the differences can be largely attributed to the receivers of different models/brands used in the generation of respective OSBs, rather than just the overall large number of the IGS stations.

Fig. 4
figure 4

GPS satellite (G066) OSB estimates of each calendar month from IGS and MGEX solutions during 2017–2018. Top, middle and bottom plots correspond to C1C, C1W and C2W OSBs, respectively, and the length of vertical lines of individual satellites denotes the respective standard deviations

Full size image

Figure 5 depicts the satellite L1/L2/L5 OSBs of the whole GPS constellation generated from the MGEX network. The OSBs have been aligned across the entire period, and their weekly mean values are presented. Individual satellites are labeled by their SVNs and sorted by Block types on the horizontal axis. The overall OSB range is about 40 ns across the entire GPS constellation, and the variation of OSBs shows dependence on Block types. For GPS L1/L2 OSBs, a significant deviation is found between Block IIF and other Blocks. The difference is about 10 ns for a given signal within Block IIF, and the scatter is at the same level in case of Blocks IIA, IIR-A and IIR-M. The range of OSBs within Block IIR-B is significantly larger than that within IIR-A and IIR-M. As for OSBs of the new civil signals on L2 frequency, i.e., C2S, C2L and C2X, they exhibit similar bias patterns. The OSBs on the L5 frequency, i.e., C5Q and C5X, are only supported by Block IIF satellites. The overall range of C5Q/C5X OSBs is within 20 ns, and the two biases agree well with each other.

Fig. 5
figure 5

Weekly aligned GPS satellite OSBs from the MGEX solution (2017–2018). Individual satellites are labeled by their SVNs for unambiguous identifications

Full size image

The weekly STDs of GPS C1W, C2W and C5Q OSBs are calculated for the individual satellite and plotted in Fig. 6, to check the stability of GPS L1/L2/L5 OSB estimates generated from the MGEX network. Except for G061, OSB estimates of Block IIA, IIR-A and IIR-B satellites reveal overall superior repeatability compared with the other two Block types. In comparison of the stability of satellite OSB estimates within Blocks IIR-M and IIF, there exist notably large STD values for some IIR-M (G048 to G050) and IIF (G067 to G069) satellites. While the OSB estimates of individual satellites show stable characteristics during most of the test period, pronounced scatters are still noticed for each satellite, which indicates the possibility of the poor stability of satellite OSB estimates in certain cases. In the use of zero-mean conditions to separate the satellite and receiver biases, it is assumed that the bias stability of all involved satellites is at a comparable level. In case the satellites with unstable bias characteristics (presumably as a result of the in-orbit status of SVs) are involved in the separation of satellite and receiver OSBs, OSB estimates of the remaining satellites might be contaminated. Considering the different levels of the variability of satellite code biases, an iterative method was reported in Li et al. (2012) to generate a reasonable set of satellites with stable bias variabilities as reference datums. Although the number of selected reference satellites might differ from the whole constellation, the resulting bias estimates do not affect the practical applications. This strategy will be further verified for the quality improvement of the derived satellite OSBs.

Fig. 6
figure 6

Weekly STDs of GPS C1W, C2W and C5Q satellite OSBs from the MGEX solution during the two-year test period

Full size image

The stability of GPS satellite OSBs generated from the IGS and MGEX networks is summarized in Table 2, in which the result of C2C OSB is missing in MGEX solution since the signal is not supported by modern GNSS receivers. The STD averages are 0.062, 0.060, 0.095 and 0.101 ns for C1C, C1W, C2W and C2C OSBs of the IGS solution, respectively, and the corresponding values of the MGEX-based OSB estimates are 0.068 ns for C1C, 0.064 ns for C1W and 0.102 ns for C2W. The STD of satellite OSBs is at the same level for the new L2 civil signals, which is 0.106 ns for C2X/C2X and 0.107 ns for C2L. The STD of OSBs at the L5 frequency is slightly larger compared with that of L1/L2 OSBs, which is 0.122 and 0.117 ns for C5Q and C5X, respectively. In view of the superior repeatability of L1/L2 OSBs generated from the IGS network, the IGS-based OSB estimates are used as references to analyze the consistency between MGEX and IGS OSB estimates.

Table 2 Mean STD of GPS and GLONASS satellite OSB estimates during 2017–2018 (unit: ns)
Full size table

The distribution of GPS satellite OSB differences between MGEX and IGS solutions during the two-year period is depicted in Fig. 7. There is no pronounced systematic deviation between the two OSB estimates since OSB differences are almost centered near zero. A slightly worse agreement is found for C2W OSB in comparison with C1C and C1W OSBs. This can be verified in the analysis of absolute OSB differences within the range of 0.3 ns. The percentage of absolute OSB differences within 0.3 ns is about 99% for C1C and C1W, which notably drops to around 93% for C2W. The decreased consistency of C2W OSBs between the two products is related at first, to the large variation range of C2W OSBs resulting from the IF constraint between C1W and C2W observations, and second, to the different data sets for L1/L2 inter-frequency SPR bias estimation in IGS and MGEX solutions. Due to the poor coverage of simultaneous tracking of GPS C1W and C2W measurements by the MGEX receivers, both C1C and C1W signals are used to form the GF combinations with respect to the C2W in MGEX OSB computation. As for the IGS OSB solution, only the GF combination between C1W and C2W signals is formed. The differences in the involved receivers (e.g., the numbers, distributions and types) introduce different levels of ionospheric modeling residuals, which might, propagate into the resulting OSB errors.

Fig. 7
figure 7

Histogram of the distribution of GPS satellite OSB differences between IGS and MGEX estimation results during 2017–2018

Full size image

For each individual satellite, the root-mean-square (RMS) of GPS OSB differences between IGS and MGEX estimates is calculated based on the daily solution during the 2-year period. As depicted in Fig. 8, GPS OSB RMS differences show no obvious dependence on the Block types. The C2W OSB exhibits the worst consistency between the two products, which mainly varies within the range of 0.3 ns. RMS differences of C1C and C1W OSB estimates are at the same level, which are below 0.2 ns. Overall, the OSB RMS difference between MGEX and IGS estimation results is 0.122 ns for C1C, 0.106 ns for C1W and 0.174 ns for C2W, respectively. It indicates that the consistency between GPS satellite OSB estimates from different global networks (MGEX vs. IGS) with mixed receivers can reach up to the level of 0.1–0.2 ns with the use of identical estimation method. Aside from the stabilization of hardware characteristics in the satellite part, the global coverage of IGS/MGEX receivers also improves the repeatability of GPS satellite OSB estimates.

Fig. 8
figure 8

GPS OSB RMS differences between IGS and MGEX solutions for individual satellites during the 2-year period

Full size image

Since the OSB computation is based on a local ionospheric modeling method, the estimation results might be affected by the different levels of solar activities. To examine the sensitivity of the method to changes of solar conditions, the IGS and MGEX OSB solutions were generated during a three-month period (August to October) of the years 2014 and 2017, respectively. The comparison of IGS- and MGEX-based GPS L1/L2 OSB estimates during high and low solar conditions is summarized in Table 3. For the IGS solution, the STD of GPS satellite OSBs varies within the range of 0.1–0.15 ns across the entire 3 months in 2014, which drops to 0.07–0.1 ns for the test period in 2017. For the MGEX solution, the STD of GPS OSB estimates is on the order of 0.13–0.2 ns during the period of high solar activity, and 0.07–0.11 ns during the period of low solar activity. The stability of MGEX OSB solution is notably larger during high solar conditions in comparison with the IGS solution. Note that about 90 and 220 stations are involved in the MGEX OSB estimation during the test period of 2014 and 2017, respectively. The station numbers are almost the same, 450 in 2014 and 470 in 2017, in the IGS OSB estimation. The inferior stability of MGEX solutions in 2014 also relates to the notably small network of MGEX receivers. Overall, the precision of GPS L1/L2 satellite OSBs during high solar conditions is about 1.5 times worse than that during low solar conditions for both MGEX and IGS solutions. The OSB RMS differences between the two products are on the level of 0.23–0.35 and 0.07–0.12 ns during high and low levels of solar activity, respectively.

Table 3 Comparison of GPS satellite OSBs estimated from the IGS and MGEX networks during high and low levels of solar activity (unit: ns)
Full size table

4.2 GLONASS OSB solutions

Like the bias analysis of GPS, the weekly aligned OSB estimates of GLONASS satellites generated from the IGS network is first depicted in Fig. 9, where the length of individual colors denotes the range of OSB variations. The OSB range is roughly 40 ns across the GLONASS constellation. One GLONASS-M satellite with slot number 6 (i.e., SVN 733) shows the largest bias variation during the test period. As for each individual GLONASS satellite, the OSB range is notably larger compared to that of GPS, which can be partly attributed to the influences of GLONASS inter-channel biases (Wanninger 2011). Different from the GPS constellation, the receiver biases of GLONASS strongly depend on the frequency channel because of the FDMA modulation employed by the GLONASS. In the analysis of satellite-receiver pair biases performed at CODE for both GPS and GLONASS (see Villiger et al. 2019), it was demonstrated that the scatter of GLONASS OSBs (1.0–2.9 ns) is roughly by a factor of two larger than GPS OSBs (0.5–1.6 ns). Since the satellite-receiver pair bias is separated into a satellite-dependent part and a receiver-dependent part in the OSB estimation of GLONASS, a high scatter of GLONASS OSBs can be expected in both satellite and receiver sides.

Fig. 9
figure 9

Weekly aligned GLONASS satellite OSB estimates from the IGS solution (2017–2018). Individual satellites are labeled by their SVNs and followed by the frequency channel numbers in the brackets

Full size image

As shown in Table 2, the stability of GLONASS satellite OSBs generated using IGS stations is 0.092 ns for C1C, 0.096 ns for C1P, 0.137 ns for C2P and 0.134 ns for C2C, respectively, across the two-year period. When calculating from the MGEX network, they are 0.118, 0.105, 0.156 and 0.153 ns for C1C, C1P, C2P and C2C OSBs, respectively. Obviously, GLONASS OBS estimates are much scatted compared with GPS OSBs. The consistency between GLONASS satellite OSBs calculated from the independent IGS and MGEX networks is investigated by analyzing the distribution of OSB differences within each individual bias bin of 0.04 ns. As depicted in Fig. 10, the slightly large range of OSB differences indicates the inferior agreement between the IGS- and MGEX-based OSB estimates of GLONASS, although no systematical offset is observed. In the analysis of the percentage of absolute OSB differences within the range of 0.5 ns, the biases on the L2 frequency band show worse consistency than those on the L1 frequency band. Similar to GPS, this behavior can be partly explained by the IF constraint between GLONASS C1P and C2P observations, which leads to the large range of OSB variations on the L2 frequency.

Fig. 10
figure 10

Histogram of the distribution of GLONASS satellite OSB differences between IGS and MGEX solutions during the period 2017–2018

Full size image

The RMS differences of GLONASS satellite OSBs between IGS and MGEX solutions are depicted in Fig. 11. The bias RMS differences of C1C and C1P are generally below 0.4 ns, which vary within the range of 0.6 ns for C2C and C2P OSBs. It is also interesting to note that there is a correlation of the bias RMS differences with frequency number, in particular for satellites located close to the upper edge of the frequency channel. The RMS difference between MGEX and IGS OSB solutions is 0.269 ns for C1C, 0.203 ns for C1P, 0.364 ns for C2C and 0.328 ns for C2P across the entire test period. Overall, the consistency between IGS- and MGEX-based satellite OSBs is 0.2–0.4 ns for GLONASS, which is roughly by a factor of two worse than GPS. Obviously, the bias estimates will be systematically affected by the differences in the involved tracking networks and employed receivers, although the same estimation method is used. This deviation behavior would be more pronounced for the OSB estimates of GLONASS compared with GPS.

Fig. 11
figure 11

GLONASS OSB RMS differences between IGS and MGEX solutions for each individual satellite during the two-year period

Full size image

4.3 Influences of different receiver groups on OSB estimation

The inter-receiver comparisons have demonstrated that the separation of pseudorange biases into a satellite-dependent and receiver-dependent part results in the inconsistent biases between receivers with different correlator types and front-end designs (Hauschild and Montenbruck 2016). Reference station networks commonly consist of geodetic receivers using specific correlators and front-end bandwidths, which may differ between different manufacturers or even between different generations of receivers from the same provider. As such, these receivers will generate pseudorange biases containing satellite-dependent biases, which are consistent within the specific group of receivers.

The impacts of different groups of receivers on the resulting OSB estimates are investigated in this section. The receivers were restricted to Septentrio, Javad and Trimble brands for both IGS and MGEX networks, but no detailed classification on antenna types was made. An estimation of GPS and GLONASS OSBs from each group of receivers was performed from June to August 2018. Among the selected IGS stations, 77 Javad, 61 Septentrio and 132 Trimble receivers are included. For the MGEX receivers, the numbers are 56 for Javad, 46 for Septentrio and 81 for Trimble.

The experimental results are presented in terms of the RMS differences between satellite OSB estimates from the independent receiver groups, which are summarized in Tables 4 and 5, corresponding to the IGS and MGEX results, respectively. Within the IGS network, the comparison between Septentrio and Javad/Trimble receivers shows similar bias RMS differences, which are on the order of 0.29–0.49 ns for GPS and 0.85–1.53 ns for GLONASS. The bias RMS difference between Javad and Trimble receivers is slightly larger, which increases to 0.45–0.81 and 1.34–1.67 ns for GPS and GLONASS, respectively. Ideally, different sets of biases need to be estimated for distinct respective groups of receivers with different receiver models (Hauschild and Montenbruck 2016). In case, the bias corrections generated from certain receiver group are applied to a different set of receivers, which exhibits distinct bias characteristics, a decrease in the performance can be expected.

Table 4 RMS differences of satellite OSBs generated from the separate types of receivers included in the IGS network (unit: ns)
Full size table
Table 5 RMS differences of satellite OSBs generated from the separate types of receivers included in the MGEX network (unit: ns)
Full size table

A similar systematic deviation between OSBs generated from different groups of receivers can be also observed in MGEX-based solutions. The comparison between Septentrio and Javad receivers presents the smallest bias inconsistency, which varies within the range of 0.23–0.38 ns in case of GPS OSBs. The bias RMS difference increases significantly to the order of 0.43–0.86 ns in the comparison of Septentrio/Javad and Trimble receivers. As for GLONASS, the bias RMS difference is at a comparable level between the three-receiver groups, which varies mainly in the range of 1.0–1.55 ns. The effect of different types of receivers is more pronounced on the estimated OSBs of GLONASS. Overall, the bias consistency of GLONASS is roughly by a factor of 2–3 worse than GPS in both IGS and MGEX analysis.

When comparing OSB estimates between the identical receiver groups, the OSB RMS differences are significantly smaller in comparison with the biases generated from different groups of receivers. Take the Trimble receiver based OSB estimates for example, the bias RMS differences between the IGS and MGEX solutions are 0.12–0.2 and 0.11–0.21 ns for GPS and GLONASS, respectively. We further restrict the IGS/MGEX Trimble receivers to the same number (70 stations) for the respective OSB computation. It is found that the bias RMS differences drop to 0.04–0.1 ns for GPS and at the same magnitude (0.1–0.22 ns) for GLONASS. The comparison between Septentrio-only, Javad-only, Trimble-only solutions and the combined solutions also reveals an interesting result that, the bias RMS differences are slightly smaller in this case than those estimated from individual receiver groups for both GPS and GLONASS. As the biases are commonly estimated in local/global ionospheric analysis under the ionospheric thin-shell assumption, they are significantly affected by thin-shell errors using fixed altitude (Themens et al. 2015). For a detailed analysis of different groups of receivers on the resulting OSB estimates, the mitigation of ionospheric modeling errors as well as the effects of different antenna types is worthy of further investigation.

5 Summary and conclusions

The maintenance of DCB parameters covering all possible signal pairs is becoming infeasible in a multi-constellation and multi-frequency GNSS environment. As a result, the observable-specific signal bias is proposed for the transition from differential to pseudo-absolute code biases in support of both real-time and post-processing applications. In this paper, the estimation of OSB parameters is achieved in local ionosphere analysis with the use of a modified GTS function for the station-specific ionospheric VTEC modeling. As an experimental step, GPS L1/L2/L5 and GLONASS L1/L2 OSBs are computed using observation data from the independent IGS and MGEX network.

The analysis of the generated GPS and GLONASS satellite OSBs is performed during a two-year period (2017–2018). To allow for the direct comparison of OSB solutions from different networks, an AutoBIA method is applied for the proper handling of discontinuities in the derived OSB series. The stability of satellite OSB estimates is on the order of 0.06–0.1 ns across the GPS constellation. Since GLONASS OSB is separated into a satellite-dependent part and a receiver-dependent part in our analysis, the value slightly increases to 0.09–0.15 ns due to the effects of GLONASS inter-channel biases. In the comparison of OSB estimates from the separate global network (IGS vs. MGEX) with mixed types of receivers, the RMS difference is at the level of 0.1–0.2 ns for GPS, which shows better bias consistency by a factor of two compared with GLONASS. GPS and GLONASS OSBs generated from three different groups of receivers of the IGS/MGEX network are also analyzed during a three-month period of 2018. The maximum RMS difference between different receiver groups reaches up to the magnitude of 0.6–0.9 and 1.4–1.7 ns for GPS and GLONASS, respectively. The OSB estimation can be systematically affected by the bias inconsistencies between receivers of different models/brands, which is more pronounced on the estimated OSBs of GLONASS.

The receiver groups were roughly classified by the manufactures in the analysis of different types of receivers on the estimation of OSBs. Since different antenna models may also result in different signal distortions, the effects of different antenna types on the resulting OSBs will be investigated using identical receivers. The multi-GNSS OSBs, which also cover Galileo, BeiDou and QZSS signals, are now routinely generated at CAS using the same estimation method. The bias characteristics between the new signals of new GNSS constellations will be analyzed and presented in future studies.