The CODE ambiguity-fixed clock and phase bias analysis products: generation, properties, and performance

Abstract

The generation and use of GNSS analysis products that allow—particularly for the needs of single-receiver applications—precise point positioning with ambiguity resolution (PPP-AR) are becoming more and more popular. A general uncertainty concerns the question on how the necessary phase bias information should be provided to the PPP-AR user. Until now, each AR-enabling clock/bias representation method had its own practice to provide the necessary bias information. We have generalized the observable-specific signal bias (OSB) representation, as introduced in Villiger (J Geod 93:1487–1500, 2019) originally exclusively for pseudorange measurements, to carrier phase measurements. The existing common clock (CC) approach has been extended in a way that OSBs allowing for flexible signal and frequency handling between multiple GNSS become possible. Advantages of the proposed OSB-based PPP-AR approach are: GNSS biases can be provided in a consistent way for phase and code measurements and it is capable of multi-GNSS and suitable for standardization. This new, extended PPP-AR approach has been implemented by the Center for Orbit Determination in Europe (CODE). CODE clock products that adhere to the integer-cycle property have been submitted to the International GNSS Service (IGS) since mid of 2018 for three analysis lines: Rapid, Final, and MGEX (Multi-GNSS Extension). Ambiguity fixing is performed not only for GPS but also for Galileo. The integer-cycle property of between-satellite clock differences is of fundamental importance when comparing satellite clock estimates among various analysis lines, or at day boundaries. Both kinds of comparisons could be exploited at a very high level of consistency. Any retrieved comparison essentially indicated a standard deviation for between-satellite clocks from CODE of the order of 5 ps (1.5 mm in range). Finally, the integer-cycle property that may be recovered between the CODE Final clock and the accompanying bias product of consecutive daily sessions (using clock estimates additionally provided for the second midnight epoch) allows us to deduce GPS satellite clock and phase bias information that is consistent and continuous with respect to carrier phase observation data over two, three, or, in principle, yet more days. Phase-based clock densification from initially estimated integer-cycle-conform clock corrections at intervals of 300 s to 30 s (5 s in case of our Final clock product) is a matter of particular interest. Based on direct product comparisons and GRACE K-band ranging (KBR) data analysis, the quality of accordingly densified clock corrections could be confirmed to be on a level similar to that of “anchor” (300 s) clock corrections.

1 Introduction

Precise point positioning (PPP) is a well-established tool to enable single-receiver users to compute their receiver positions with an accuracy conforming to geodetic network standards (Zumberge et al. 1997; Kouba and Héroux 2001). More and more analysis centers of the International GNSS Service (IGS, Johnston et al. 2017) started to provide precise satellite orbit and “integer,” or ambiguity-fixed clock products, usually along with information on phase biases, potentially enabling integer ambiguity resolution (IAR or briefly AR) for PPP users. The methods that have been developed over the last decade for the representation of clock and appropriate bias information, however, are manifold. We would like to refer here primarily to the contribution by Teunissen and Khodabandeh (2015), where a comprehensive overview and comparison of documented IAR-enabling representation methods are given.

Integer-ambiguity, or simply “integer” PPP (IPPP) is a synonym commonly used for IAR-enabled PPP. Within the IGS, PPP-AR (PPP with AR) is used as generic term addressing this particular refinement of PPP and its underlying analysis results. An essential feature of PPP-AR is that it is a relative technique where the single receiver is processed relative to the global network solution providing the satellite clock corrections. In the context of PPP-AR this means that the single-receiver users’ integer phase ambiguities correspond in fact to double-differenced ambiguities (Teunissen and Khodabandeh 2015). PPP-AR, therefore, has to be considered an undifferenced analysis strategy that competes directly against traditional ambiguity-fixed double-difference analysis. This is one reason why PPP-AR becomes a key tool for single-receiver users that are processing GNSS observation data originating from remote receivers.

A general uncertainty concerns the question on how the necessary phase bias information should be provided to the PPP-AR user. Until now, each AR-enabling clock/bias representation method had its own practice to provide the necessary bias information. We consider the observable-specific signal bias (OSB) representation as a suitable starting point for a generally applicable bias parameterization for PPP-AR. That is our motivation to proposing to use the OSB bias representation, as introduced in Villiger et al. (2019) originally exclusively for pseudorange measurements, concurrently for all involved types of GNSS measurements, also including carrier phase measurements. This finally enabled us to further extend the existing common clock (CC) approach in a way that OSBs allowing for flexible signal and frequency handling between multiple GNSS becomes possible. Further key benefits of the proposed PPP-AR bias description are: It is capable of multi-GNSS, easy to use, suitable for standardization, and it is incidentally relying on a GNSS bias data format officially approved by the IGS. A particular PPP-AR bias conditioning should finally become distinct when carrier phase as well as pseudorange bias data is considered simultaneously in the bias data interface as offered, e. g., by Bias-SINEX (Schaer 2018). Besides, the bias sign convention is a crucial question of detail we would like to address.

The new, extended PPP-AR approach could be implemented by the Center for Orbit Determination in Europe (CODE, Dach et al. 2009), a joint venture of the Astronomical Institute of the University of Bern (AIUB, Bern, Switzerland), the Swiss Federal Office of Topography (swisstopo, Wabern, Switzerland), the Federal Agency for Cartography and Geodesy (BKG, Frankfurt am Main, Germany), and the Institut für Astronomische und Physikalische Geodäsie of the Technische Universität München (IAPG/TUM, Munich, Germany). CODE clock products that adhere to the integer-cycle property have been submitted to the International GNSS Service (IGS) since mid of 2018 for three analysis lines (Schaer et al. 2018): Rapid, Final, and MGEX (Multi-GNSS Extension, Montenbruck et al. 2017). The fact that single-receiver ambiguity resolution is being performed not only for GPS but also for Galileo can be seen as another milestone in terms of IGS analysis. In addition, the IAR-specific procedures in the various clock analysis lines of CODE are largely isolated and kept independent. This allows us to directly observe the exclusive impact of successful between-satellite ambiguity fixing on CODE clock products in the operational comparisons and combinations of IGS legacy clock products.

For the GPS constellation, integer-cycle clock results that are coming from three different analysis lines are at our disposal. We took this opportunity to recover their integer property and thus to examine their consistency among each other. An important step for the IGS toward the PPP-AR domain was made by Banville et al. (2020), where a preliminary combination of PPP-AR products from six IGS analysis centers over a 1-week period could be successfully achieved. This effort may be seen as a rigorous generalization of those comparisons we accomplished among the available analysis lines from our side. It should be noted that CODE contributed with its IGS Final analysis products to the mentioned IGS PPP-AR experiment. This paper shall help to make another step toward standardization in the PPP-AR domain by consistently relying on standards for OSB-based PPP-AR (as they should also apply to Bias-SINEX). Such a standardization has to be considered beneficial especially for all PPP-AR users, which thereby should get a common interface for standardized, easily applicable GNSS bias data that, moreover, may be expected to be consistent to existing IGS clock analysis products. The same should be valid for an IGS-combined integer clock/bias product that might eventually be strived for (Banville et al. 2020).

The availability of consistent and continuous orbit and clock analysis products that are covering session lengths of 30 h or longer is a request that has been expressed for a long time from the analysis community of LEO-originated GNSS data (e. g., Bock et al. 2007). Our new IGS Final analysis products actually allow to straightforwardly generate continuous integer-cycle clock information over 72 h on the basis of the clock and phase bias information of consecutive daily sessions.

The choice of the sampling interval for satellite clock estimates is an issue that gets rarely addressed in the literature, even though it is an utmost crucial one when generating high-quality GNSS clock corrections (Bock et al. 2009). That is why clock densification, particularly with regard to integer-cycle products, is an issue we will focus on in this contribution too.

The paper is organized as follows. Section 2 describes the relevant observation equations and the transformational links between OSB and widelane/narrowlane fractional phase biases. Section 3 describes the algorithm we are proposing and highlights the necessary steps for generation of ambiguity-fixed clocks and consistent phase bias information. Section 4 gives an overview of CODE’s improved IGS/MGEX clock analysis products that adhere to the integer-cycle property. Section 5 gives the basic equations for accessing the integer-cycle property of satellite clock estimates; it includes various validations and numerous comparisons, where this particular property could be accessed successfully; besides of that, Sect. 5 demonstrates that we are able to straightforwardly generate continuous integer clock information over three consecutive daily sessions and it finally addresses the question whether densification of clock corrections (by analyzing time differences of carrier phase observation data) is harmful for an ambiguity-fixed clock product. Section 6 summarizes the results in terms of method and implementation at CODE. A short supplementary guideline on how to use Bias-SINEX for the proposed OSB-based PPP-AR application is included as an appendix.

2 OSB code and phase bias representation for common clock (CC) approach

A primary precondition for all further considerations is that receiver clock ($t_k$) and satellite clock ($t^i$) corrections are assumed to be common for carrier phase and pseudorange code signals. Let us generally use B to describe a particular GNSS measurement bias. The bias is related to units of time. Typically, it needs to be estimated. In the following, bias parameters are all treated at a pseudo-absolute, undifferenced level, which, incidentally, is in complete agreement with the situation regarding clock offsets. Original dual-frequency GNSS measurements of carrier phase (L) and pseudorange code (P) signals expressed in units of length from station k to satellite i can be written as (Teunissen and Kleusberg 1998)

$$\begin{aligned} {L_1}_k^i&= \tilde{\rho }_k^i\! -\!\frac{I_k^i}{\nu _1^2}\! + c\left( t_k-t^i\right) \! +\!\lambda _1 {N_1}_k^i+c\left( {B_{L_1}}_k\!+\!B_{L_1}^i\right) \nonumber \\ {L_2}_k^i&= \tilde{\rho }_k^i\! -\!\frac{I_k^i}{\nu _2^2}\! + c\left( t_k-t^i\right) \! +\!\lambda _2 {N_2}_k^i+c\left( {B_{L_2}}_k\!+\!B_{L_2}^i\right) \nonumber \\ {P_1}_k^i&= \tilde{\rho }_k^i\! +\!\frac{I_k^i}{\nu _1^2}\! + c\left( t_k-t^i\right) \! +c\left( {B_{P_1}}_k\!+\!B_{P_1}^i\right) \nonumber \\ {P_2}_k^i&= \tilde{\rho }_k^i\! +\!\frac{I_k^i}{\nu _2^2}\! + c\left( t_k-t^i\right) \!+c\left( {B_{P_2}}_k\!+\!B_{P_2}^i\right) \,,\nonumber \\ \end{aligned}$$

(1)

where $\tilde{\rho }_k^i$ contains the geometrical distance between station k and satellite i at signal reception time $t_k$ and signal emission time $t^i$, respectively, including troposphere delay and corrections due to receiver and satellite antenna phase center calibrations; $I_k^i/\nu ^2$ denotes the first-order ionosphere delay, namely for the signal frequencies $\nu _1$ and $\nu _2$; c is the speed of light in vacuum; $t_k$ and $t^i$ are the receiver and satellite clock offsets, respectively; $\lambda _1=c/\nu _1$ and $\lambda _2=c/\nu _2$ are the wavelengths associated with the two frequencies; ${N_1}_k^i$ and ${N_2}_k^i$ correspond to undifferenced integer carrier phase ambiguities for the first and the second frequency. We want to emphasize that frequency numbers 1 and 2 are just a counting index and simply stand for two of the signals of the respective GNSS. Last but not least, each observation equation includes a specific signal bias term, $B_k+B^i$, consisting of a receiver component, $B_k$, and a satellite component, $B^i$, each expressed in units of time and each finally associated with a particular observation type ($L_1$, $L_2$, $P_1$, $P_2$). All signal bias terms introduced in Eq. (1) are to be interpreted as receiver and satellite hardware delays. Those concerning the phase signals are sometimes also called “the non-integer parts of ambiguities.” Other effects, such as, e. g., measurement noise, multipath, phase wind-up, are ignored for clarity.

The basic set of undifferenced observation equations, as introduced in Eq. (1) for the purpose of bias treatment related to PPP-AR, is specifically applicable to GPS and Galileo. In fact, it might be adopted to any GNSS following a CDMA (code-division multiple access) channel access method. With the objective of promoting Eq. (1) for the use of OSB-based PPP-AR, it should be noted that our observation equations are quite similar to those adopted by Geng et al. (2019). However, in contrast to their equations, we tried to establish a more intuitive, self-explanatory nomenclature for the multitude of bias terms. Apart from this, a deciding difference affects the sign for all satellite bias terms in Eq. (1): $+c\,B_{L_1}^i$, $+c\,B_{L_2}^i$, $+c\,B_{P_1}^i$, $+c\,B_{P_2}^i$. Specifically, these bias terms are vitally important for the OSB-based PPP-AR application. There are apparently different conventions for the sign of biases in the literature. In this article we refer to that of Eq. (1), which incidentally is also that of Bias-SINEX (Schaer 2018):

$$\begin{aligned} +\text {Bias} = \text {Observed} - \text {Computed}\,, \end{aligned}$$

(2)

or, in agreement with the arrangement in Eq. (1),

$$\begin{aligned} \text {Observed} = \text {Computed} + \text {Bias}\,. \end{aligned}$$

(3)

In addition, it should generally apply that satellite and receiver bias values should be considered with the same sign rule:

$$\begin{aligned} B_{\text { total}} = B_{\text { satellite}} + B_{\text { receiver}}\,. \end{aligned}$$

(4)

2.1 Melbourne–Wübbena linear combination

Single-difference ambiguities ($N_k^{ij} = N_k^i-N_k^j$) between pairs of satellites (i, j) have to be considered to eliminate receiver bias terms for single-receiver ambiguity resolution, The observation equation for the Melbourne–Wübbena (MW) linear combination (LC) of phase and code measurements (Melbourne 1985; Wübbena 1985), accordingly formed for equivalent single differences, enables, in a first step, ambiguity resolution (AR) for between-satellite widelane ambiguities, subsequently substituted with ${N_W}_k^{ij}={N_1}_k^{ij}-{N_2}_k^{ij}$:

$$\begin{aligned} \text {MW}_k^{ij} =\,&\kappa _{L_1} {L_1}_k^{ij} + \kappa _{L_2} {L_2}_k^{ij} + \kappa _{P_1} {P_1}_k^{ij} + \kappa _{P_2} {P_2}_k^{ij}\nonumber \\ =\,&\lambda _W \left( {N_1}_k^{ij}\!-\!{N_2}_k^{ij}\right) + c \left( \kappa _{L_1} B_{L_1}^{ij}\!+\!\kappa _{L_2} B_{L_2}^{ij}\right) \nonumber \\&+ c \left( \kappa _{P_1} B_{P_1}^{ij}\!+\!\kappa _{P_2} B_{P_2}^{ij}\right) \nonumber \\ =\,&\lambda _W {N_W}_k^{ij} + c\,B_{L_W}^{ij} + c\,B_{P_W}^{ij}\nonumber \\ =\,&\lambda _W \left( {N_W}_k^{ij} + b_{L_W}^{ij}\right) + c\,B_{P_W}^{ij}\,, \end{aligned}$$

(5)

where $\kappa _{L_1} = +\frac{\nu _1}{\nu _1-\nu _2}$, $\kappa _{L_2} = -\frac{\nu _2}{\nu _1-\nu _2}$, $\kappa _{P_1} = -\frac{\nu _1}{\nu _1+\nu _2}$, $\kappa _{P_2} = -\frac{\nu _2}{\nu _1+\nu _2}$ are the prefactors of the MW LC for phase (L) and code (P) measurements (approx. $+4.53$, $-3.53$, $-0.56$, $-0.44$ for GPS); $\lambda _W = \frac{c}{\nu _1-\nu _2}$ is the widelane wavelength (approx. 86.2 cm for GPS); $B_{L_W}^{ij}$, or $b_{L_W}^{ij}$ describe the between-satellite widelane fractional cycle bias (FCB) coming from the phase measurements ($L_1$ and $L_2$); $B_{P_W}^{ij}$ describes the between-satellite differential widelane bias induced by the code measurements ($P_1$ and $P_2$). Both $B_{L_W}^{ij}$ and $B_{P_W}^{ij}$ are expressed in units of time, whereas $b_{L_W}^{ij}$ (in lowercase) is kept in units of widelane ($\lambda _W$) cycles, given that $c\,B_{L_W}^{ij} = \lambda _W b_{L_W}^{ij}$. It is worth mentioning that the characteristics of GPS satellite-satellite widelane fractional biases ($b_{L_W}^{ij}$) was the subject of investigations already long time ago (see, e. g., Gabor and Nerem 2004).

In either case (including the double-difference case), widelane AR—when relying on the MW LC—is subject to differential code biases (DCBs). This is indicated in Eq. (5) with the code bias term $B_{P_W}^{ij} = B_{P_W}^i-B_{P_W}^j$ or, written in extenso,

$$\begin{aligned} B_{P_W}^{ij} = \kappa _{P_1} \left( B_{P_1}^{i}-B_{P_1}^{j}\right) + \kappa _{P_2} \left( B_{P_2}^{i}-B_{P_2}^{j}\right) \,. \end{aligned}$$

(6)

In our approach, GNSS code bias values are considered at a pseudo-absolute and, as a result, undifferenced level. They are directly accessible in the form of $B_{P_1}^i$ and $B_{P_2}^i$, for each satellite i and, in particular, for each code observation type involved. Such a bias parameter representation is commonly called OSB (observable-specific signal bias) representation (for more details, see Villiger et al. 2019). It can be easily verified that errors (or uncertainties) for GPS-related $B_{P_1}$ and $B_{P_2}$ should not exceed about $\pm 0.51$ ns and $\pm 0.66$ ns, respectively, to keep a DCB-induced widelane fractional bias error below a tenth ($\pm 0.10$) of a widelane cycle.

2.2 Relevance of single-differenced receiver phase and code biases to widelane AR

Receiver biases, ${B_{L_1}}_k$, ${B_{L_2}}_k$, ${B_{P_1}}_k$, ${B_{P_2}}_k$, were specified in Eq. (1) intrinsically for completeness. These biases are assumed to get eliminated after between-satellite differencing. This basic assumption is generally satisfied. There are cases, however, where it is no longer appropriate to comply with this general constraint. Such cases basically may happen as soon as signal observable types get involved that are different among a particular pair of observed satellites (and thus have to be considered potentially “unequal” in terms of their bias characteristics). This means that the bias elimination assumption may become wrong when phase or code signal observable types get mixed among each other:

$$\begin{aligned} \begin{aligned} B_{L_Wk}^{(ij)}=\,&B_{L_Wk}^{(i)}-B_{L_Wk}^{(j)}\ne 0\\ B_{P_Wk}^{(ij)}=\,&B_{P_Wk}^{(i)}-B_{P_Wk}^{(j)}\ne 0\,. \end{aligned} \end{aligned}$$

(7)

Note that a superscript in parentheses, e. g. $^{(i)}$, is used to distinguish between receiver-induced biases assigned to different satellites (i and j). In practice, corresponding receiver biases are commonly unknown and thus unavailable to a typical PPP-AR user. However, cases are conceivable in which the simplification initially made is no longer possible.

• For phase observables, single-difference widelane phase biases may become relevant as soon as a phase observation difference is affected by a quarter-cycle phase bias (e. g., concerning GPS L2C versus L2W). Such a 0.25-cycle bias gets directly transferred into the MW observation (leading to a 0.25-cycle bias also there).

• For code observables, single-difference widelane code biases may become relevant as soon as we get confronted with a mixture of different code observation types (e. g., C2C versus C2W for GPS). The resulting bias in the MW observation is a function of the diverse receiver code bias components that are present. The threshold values (we previously estimated) might help to decide whether widelane AR shall be attempted for a specific code observation difference that has to be assessed accordingly as “inconsistent.”

To be on the safe side, the PPP-AR user is well advised to omit widelane AR for specific single-difference combinations in any case of doubt (as we highlighted separately in the latter paragraph). It is obvious that omitting widelane AR for a particular widelane ambiguity (${N_W}_k^{ij}$) will make integer AR for the respective narrowlane ambiguity (${N_1}_k^{ij}$) impossible.

2.3 Ionosphere-free linear combination

The between-satellite single-difference observation equation for the ionosphere-free LC ($L_0$) of carrier phase measurements, ${L_0}_k^{ij}$, may be formed based on the LC factors, $\kappa _1$ and $\kappa _2$. As long as widelane integers ${N_W}_k^{ij}$ are substituted, with ${N_1}_k^{ij}-{N_W}_k^{ij}$, for the ambiguity term ${N_2}_k^{ij}$ belonging to the second frequency (typically the one with the higher noise level), this specific phase observation is called narrowlane. Provided that the underlying phase biases $B_{L_1}^{ij}$ and $B_{L_2}^{ij}$ for the satellite pair (i, j) are conditioned consistently, the ambiguity term ${N_1}_k^{ij}$, now addressed narrowlane, is of integer nature and, therefore, should enable integer AR. The narrowlane phase observation equation to be used for the purpose of single-receiver narrowlane AR may be formulated, in conformity with ${L_0}_k^{ij}={L_0}_k^i-{L_0}_k^j$, as

$$\begin{aligned} \begin{aligned} {L_0}_k^{ij} =&\kappa _1 {L_1}_k^{ij} + \kappa _2 {L_2}_k^{ij} = \tilde{\rho }_k^{ij} - c\,t^{ij} + \lambda _N {N_1}_k^{ij}\\&- \kappa _2 \lambda _2 {N_W}_k^{ij} + c \left( \kappa _1 B_{L_1}^{ij}\!+\!\kappa _2 B_{L_2}^{ij}\right) \\ =&\tilde{\rho }_k^{ij} - c\,t^{ij} + \lambda _N \left( {N_1}_k^{ij} - \kappa _{L_2} {N_W}_k^{ij}\right) + c\, B_{L_N}^{ij}\\ =&\tilde{\rho }_k^{ij} - c\,t^{ij} + \lambda _N \left( {N_1}_k^{ij} + b_{L_N}^{ij} - \kappa _{L_2} {N_W}_k^{ij}\right) \,,\\ \end{aligned} \end{aligned}$$

(8)

where $\kappa _1 = +\frac{\nu _1^2}{\nu _1^2-\nu _2^2}$ and $\kappa _2 = -\frac{\nu _2^2}{\nu _1^2-\nu _2^2}$ are the factors to form the ionosphere-free LC (approx. $+2.55$ and $-1.55$ for GPS); $\tilde{\rho }_k^{ij}=\tilde{\rho }_k^i-\tilde{\rho }_k^j$ is the differenced geometrical term (including troposphere delay); $t^{ij}=t^i-t^j$ represents the differenced satellite clock offset, being entirely consistent to $-c\,t^{ij}=+c\,\left( t^j-t^i\right)$ (cf. Eq. (1)); $\lambda _N = \frac{c}{\nu _1+\nu _2}$ is the resulting narrowlane wavelength (approx. 107.0 mm for GPS); $\kappa _{L_2} = -\frac{\nu _2}{\nu _1-\nu _2}$ (approx. $-3.53$ for GPS) is a factor we already used in Eq. (5); $B_{L_N}^{ij}$ and $b_{L_N}$, respectively, describe the between-satellite narrowlane fractional cycle bias (FCB); $B_{L_N}^{ij}$ is expressed in units of time, whereas $b_{L_N}^{ij}$ (in lowercase) is kept in units of narrowlane ($\lambda _N$) cycles, given that $c\,B_{L_N}^{ij} = \lambda _N b_{L_N}^{ij}$.

After successful ambiguity resolution with respect to widelane (${N_W}_k^{ij}$) and subsequently narrowlane (${N_1}_k^{ij}$) ambiguities, the ionosphere-free LC ($L_0$) of ambiguity-fixed phase measurements can be analyzed (taking into account all resolved integers with respect to the original carrier phase ambiguities ${N_1}_k^{ij}$ and ${N_2}_k^{ij} = {N_1}_k^{ij} - {N_W}_k^{ij}$):

$$\begin{aligned} {L_0}_k^{ij} = \tilde{\rho }_k^{ij}-c\,t^{ij} +\kappa _1\lambda _1 {N_1}_k^{ij}+\kappa _2 \lambda _2 {N_2}_k^{ij} + c\,B_{L_N}^{ij}\,, \end{aligned}$$

(9)

with

$$\begin{aligned} B_{L_N}^{ij} = \kappa _1\left( B_{L_1}^i-B_{L_1}^j\right) +\kappa _2\left( B_{L_2}^i-B_{L_2}^j\right) \,, \end{aligned}$$

(10)

where $B_{L_N}^{ij}$, $B_{L_1}^{i}$, $B_{L_2}^{i}$ are the underlying phase biases, expressed in units of time, associated with the narrowlane, $L_1$, and $L_2$ observations. It should be emphasized that Eq. (9) is equivalent to Eq. (8). The only difference concerns the fact that Eq. (9) is formulated on the basis of the original $({N_1}_k^{ij},{N_2}_k^{ij})$ ambiguity domain, whereas Eq. (8) is referred to the narrowlane/widelane $({N_1}_k^{ij},{N_W}_k^{ij})$ ambiguity domain.

The between-satellite single-difference observation equation for the ionosphere-free LC ($P_0$) of pseudorange code measurements, ${P_0}_k^{ij}$, finally reads as

$$\begin{aligned} {P_0}_k^{ij} = \kappa _1 {P_1}_k^{ij}+\kappa _2 {P_2}_k^{ij} = \tilde{\rho }_k^{ij}-c\,t^{ij}+c\,B_{P_0}^{ij}\,, \end{aligned}$$

(11)

with $B_{P_0}^{ij} = \kappa _1 B_{P_1}^{ij} + \kappa _2 B_{P_2}^{ij}$, representing code biases conforming with the ionosphere-free LC ($P_0$) and the two original pseudorange code signals ($P_1$ and $P_2$).

We have to be aware of the given fact that $B_{P_0}^{ij}$ bias parameters play an essential role for definition of both the bias datum and the clock datum. Within the IGS, GPS satellite clock analysis results are commonly assumed to be consistent with the ionosphere-free LC from GPS C1W and C2W code measurements. For $P_1^i$ and $P_2^i$ code bias results following Villiger et al. (2019), the $P_0$-specific bias term is controlled so that it completely vanishes (i. e., $B_{P_0}^{i} = 0$) as long as reference signal types of a particular GNSS are considered. This implies that, when relying on our PPP-AR clock/bias (de)composition, Eq. (11) could actually be reduced for consideration of all chosen reference signal types (in particular for GPS C1W and C2W code measurements) to simply

$$\begin{aligned} {P_0}_k^{ij} = \tilde{\rho }_k^{ij} - c\,t^{ij}\,. \end{aligned}$$

(12)

2.4 Integer clocks

The term integer recovery clock (IRC) was originally introduced by the CNES/CLS analysis group (Laurichesse et al. 2009; Loyer et al. 2012). “Integer clocks” or, to be more specific, integer-cycle satellite-to-satellite clock differences, subsequently denoted with $t_N^i-t_N^j=t_N^{ij}$, are constructed to permit narrowlane AR without the need of applying any further phase bias information. By re-ordering the existent terms of Eq. (9), the mechanism with respect to the narrowlane phase bias term, $B_{L_N}^{ij}$, may be easily understood:

$$\begin{aligned} {L_0}_k^{ij} = \tilde{\rho }_k^{ij} \underbrace{-c\,t^{ij} + c\,B_{L_N}^{ij}}_{\textstyle -c\,t_N^{ij}} + \kappa _1\lambda _1 {N_1}_k^{ij} + \kappa _2\lambda _2 {N_2}_k^{ij}\,. \end{aligned}$$

(13)

Using the substituted clock-determined term of Eq. (13), the relationship between common clock difference, $t^i-t^j=t^{ij}$, narrowlane phase bias, $B_{L_N}^{ij}$, and integer clock difference, $t_N^i-t_N^j=t_N^{ij}$, becomes clear. Moreover, it can be noticed that integer clocks, $t_N^i$, for satellite i may be deduced directly from $t^i$ and $B_{L_N}^i$ at zero-difference level:

$$\begin{aligned} t_N^i = t^i - B_{L_N}^i = t^i - \left( \kappa _1 B_{L_1}^i + \kappa _2 B_{L_2}^i\right) \,. \end{aligned}$$

(14)

Existing datum biases (due to the absolute alignment) concerning $t^i$ and $B_{L_N}^i$ may be expected to cancel out after between-satellite differencing of accordingly differenced integer clock information ($t_N^i-t_N^j=t_N^{ij}$).

Figure 1 visually compares the common clock/OSB (CC-OSB) representation we propose in this paper (and implemented at CODE) with the integer recovery clock (IRC) representation, specifically in the (appropriately signed) range, clock, and bias domain.

Fig. 1

Common clock/OSB (CC-OSB) representation compared to the integer recovery clock (IRC) representation (in the range/clock/bias domain). Code-type observations are marked in blue, phase-type observations, assumed to be ambiguity-aligned, in red; specific OSB biases are indicated in green

CC-OSB clock/bias information may be consistently converted to any sort of (code and phase) observations, whereas IRC clock information exhibits strong consistency to phase observations only (therefore often also called “phase clocks”). For the CC-OSB clock/bias (de)composition, the situation is illustrated for three sorts of GPS observations:

• C1W/C2W code observations, unbiased with respect to the raw (or reference) clock information ($t^{ij}$);

• C1C/C2W code observations, referring to a redefined clock offset laterally shifted by a correction due to C1W−C1C DCB, which in turn corresponds to $B_{P_0}^{ij}$ of Eq. (11);

• L1W/L2W phase observations, referring to a redefined clock offset laterally shifted by a narrowlane fractional cycle bias (NL FCB), which in turn corresponds to $B_{L_N}^{ij}$ of Eq. (10).

It is important to note that narrowlane phase observations, which are shown with red arrows in Fig. 1, are assumed to be ambiguity-aligned, or consequently unambiguous. Satellite clock information that has an integer property, or that is alternatively referred to as phase clock information, may essentially be assumed to be consistent (apart from an arbitrary integer number of narrowlane cycles) among different PPP-AR clock representation methods. IRC clock correction values, according to their declaration, refer directly to (uncorrected) narrowlane phase observations. In the example in Fig. 1, CC-OSB phase clock and IRC clock correction values nominally differ by one narrowlane cycle. However, it is also apparent in Fig. 1 that the IRC approach is unable to provide a clear, well-defined reference to code observations.

2.5 Transformation between OSB and widelane/narrowlane fractional phase biases

The transformation of widelane/narrowlane fractional phase bias parameters $[B_{L_W}^{ij},B_{L_N}^{ij}]$ into OSB phase bias parameters $[B_{L_1}^{ij},B_{L_2}^{ij}]$ has to be considered a decisive factor, specifically for two reasons:

• 1. For a PPP-AR clock/bias product provider, the bias transformation from the $[B_{L_W}^{ij},B_{L_N}^{ij}]$ back into the $[B_{L_1}^{ij},B_{L_2}^{ij}]$ domain must be accomplished. $[B_{L_1}^{ij},B_{L_2}^{ij}]$ are expected to be generally applicable to nominal phase observations $L_1$ and $L_2$.

• 2. For a PPP-AR user of a correspondingly generated clock/bias product, an essential advantage is that the resulting signs for $[B_{L_1}^{ij},B_{L_2}^{ij}]$ do not depend on a particular sign convention internally followed by a specific PPP-AR product provider (as imposed, e. g., by an opposed sign for the widelane ambiguity term ${N_2}_k^{ij}-{N_1}_k^{ij}$). $[B_{L_1}^{i},B_{L_2}^{i}]$ are strictly expected to be unambiguously conform to the basic set of observation equations we initially specified in Eq. (1).

As corresponding bias arithmetics may be applied not only to single-differenced but also to undifferenced phase bias parameters, it is formulated here generally for bias parameters at zero-difference level, namely for $[B_{L_W}^{i},B_{L_N}^{i}]$ and for $[B_{L_1}^{i},B_{L_2}^{i}]$. Looking at Eqs. (5) and (8) and taking over the MW and $L_0$ LC factors, the OSB-to-widelane/narrowlane bias transformation for undifferenced satellite phase biases $B^i$ (again expressed in units of time) may be summarized with

$$\begin{aligned} \left[ \begin{array}{l}B_{L_W}^{i}\\ B_{L_N}^{i}\end{array}\right] =&\left[ \begin{array}{ll}\kappa _{L_1} &{} \kappa _{L_2}\\ \kappa _1 &{} \kappa _2\end{array}\right] \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] \nonumber \\ =&\left[ \begin{array}{ll}+\frac{\nu _1}{\nu _1-\nu _2} &{} -\frac{\nu _2}{\nu _1-\nu _2}\\ +\frac{\nu _1^2}{\nu _1^2-\nu _2^2} &{} -\frac{\nu _2^2}{\nu _1^2-\nu _2^2}\end{array}\right] \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] \nonumber \\ \approx&\left[ \begin{array}{ll}+4.53 &{} -3.53\\ +2.55 &{} -1.55\end{array}\right] \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] \,. \end{aligned}$$

(15)

By inverting the coefficient matrix of Eq. (15), the widelane/narrowlane-to-OSB bias transformation finally results in

$$\begin{aligned} \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] =&\left[ \begin{array}{ll}-\frac{\nu _2}{\nu _1} &{} +\frac{\nu _1+\nu _2}{\nu _1}\\ -\frac{\nu _1}{\nu _2} &{} +\frac{\nu _1+\nu _2}{\nu _2}\end{array}\right] \left[ \begin{array}{l}B_{L_W}^{i}\\ B_{L_N}^{i}\end{array}\right] \nonumber \\ \approx&\left[ \begin{array}{ll}-0.78 &{} +1.78\\ -1.28 &{} +2.28\end{array}\right] \left[ \begin{array}{l}B_{L_W}^{i}\\ B_{L_N}^{i}\end{array}\right] \,. \end{aligned}$$

(16)

It is important to notice that the numerical examples in Eq. (15) and (16) reflect the GPS ordinary case.

3 Generation of phase biases and ambiguity-fixed clocks

3.1 Introductory remarks

IGS/MGEX data analyses at CODE are generally accomplished using the latest development version of the Bernese GNSS Software (Dach et al. 2015a). For the present contribution, this version is specifically able to

• solve for OSB code bias parameters,

• correct for given OSB code and phase bias values,

• determine a set of fractional phase biases based on ambiguity-float(ing) single-receiver (PPP) analysis results for a significant number of ground stations,

• resolve single-difference between-satellite phase ambiguities in PPP analysis mode,

• reintroduce single-difference between-satellite phase integers in PPP analysis mode.

The first two capabilities could already be acquired as part of the developments achieved in Villiger et al. (2019). Correcting for given OSB phase bias values, as it is included in the second bullet, turned out to be just a generalization regarding OSB handling. The algorithm dedicated to solving the task listed in the third bullet will be introduced shortly in Sect. 3.3. Single-receiver AR and the necessary accounting of integer phase ambiguities could be implemented as an extension of the existing double-difference AR.

3.2 Basic principles

We are assuming perfectly known geometry (especially regarding orbital information) from a double-difference network analysis solution with AR. Figure 2 gives a schematic diagram of those processing sequences that are critical for GNSS satellite clock estimation.

Fig. 2

Generation scheme for GNSS satellite clock estimation at CODE. The original part is indicated in blue, the added part—exclusively dedicated to the IAR task—is indicated in red

These sequences include

• determination of OSB satellite code biases $[B_{P_1}^i,B_{P_2}^i]$,

• GNSS satellite clock estimation of $t^i$ at intervals of 300 sec, conducted ambiguity-fixed, if admissible, else conduced ambiguity-float (specifically for the initial iteration 0),

• all required analysis steps for IAR and provision of OSB satellite phase biases $[B_{L_1}^i,B_{L_2}^i]$,

• phase-based densification of clocks from 300 sec to 30 sec sampling (Bock et al. 2009).

At CODE, GNSS satellite clocks and phase biases are generated completely independent for each daily (24 h) session. Bias parameters of all types are generally treated constant over a daily (24 h) session. Furthermore, phase and code observation data is always considered simultaneously (specifically for clock estimation and for all narrowlane-specific analysis steps). The processing loop, consisting of clock estimation and phase bias determination and IAR, actually could be conducted repeatedly. We will make use of this particular feature for a specific experiment in Sect. 3.4.

The processing part contained in the red box of Fig. 2 may be considered isolated and, therefore, could be segregated into in a separate processing procedure. This procedure includes the following processing steps:

• 1. widelane phase bias determination,

• 2. widelane integer AR,

• 3. narrowlane phase bias determination,

• 4. narrowlane integer AR,

• 5. phase bias conversion from $[B_{L_W}^i,N_{L_N}^i]$ to $[B_{L_1}^i,N_{L_2}^i]$, according to Eq. (16),

• 6. generation of various statistics and a summary related to IAR.

For widelane and narrowlane integer AR, we apply what we call a sigma-dependent AR strategy (Dach et al. 2015a), which checks explicitly all possible combinations of satellite-to-satellite ambiguities and finally prioritizes those ambiguity parameters with the smallest formal error. During that partial AR, inversion of the (remaining) normal equation matrix is performed sequentially after having a certain number of accordingly prioritized ambiguities resolved and fixed. The method used is basically equivalent to the double-difference AR method. The CPU time of our PPP-AR procedure typically takes a couple of minutes only (having the 4 times n PPP analysis subprocesses concurrently distributed to a significant number of computing nodes, where n is the number of stations analyzed). This portion of processing time is almost negligible compared to the one needed for the actual GNSS clock estimation.

3.3 Algorithm for retrieval of fractional phase biases

Retrieval of fractional phase biases, in particular of

• widelane fractional phase biases, $b_{L_W}^{ij}$ of Eq. (5), and

• narrowlane fractional phase biases, $b_{L_N}^{ij}$ of Eq. (8),

has to be seen as an important key task. We developed a dedicated algorithm that works on the basis of estimated float (real-valued) ambiguity parameters from n station-specific PPP solutions. These ambiguity parameters correspond to either ${N_W}_k^{i}$ of Eq. (5) or ${N_1}_k^{i}$ of Eq. (8). Code bias information, represented by $B_{P_W}^{ij}$ in Eq. (5), is assumed to be invariant in this context (especially when analyzing the MW LC).

Satellite-to-satellite differences, $N_{k}^{ij}=N_{k}^{i}-N_{k}^{j}$, should be formed for all possible combinations (i, j). These differences are then reduced by means of median statistics over two stages. Median values for real-valued $N_{k}^{ij}$ are determined for each individual station k first, then median values for $N^{ij}$ are derived from the initially determined global set of real-valued ambiguity differences $N_{k}^{ij}$. The latter preparation step yields a set of satellite-to-satellite fractional phase bias values ($-0.5\,\text {cyc}\le b_{L_{W\!/\!N}}^{ij}<+0.5\,\text {cyc}$) for all possible satellite pairs belonging to a particular GNSS. Note that $b_{L_{W\!/\!N}}^{ij}$ is used in the following to symbolize either a widelane ($b_{L_{W}}^{ij}$) or a narrowlane ($b_{L_{N}}^{ij}$) phase bias. Such a set of preliminary phase bias estimates (or better pseudo-observations) is finally fed into a weighted least-squares adjustment (with a priori variances, $\sigma ^2$, assumed to be proportional to 1/n, where n is the number of individual phase bias samples that actually contributed to each satellite combination (i, j)). Until here, the phase bias datum at zero-difference level has been chosen arbitrarily (e. g., by setting the first element of $b_{L_{W\!/\!N}}^i$ to zero).

A finishing step, a refinement with regard to bias datum, can be reached by finding a common shift ($b_0$) that is minimizing the L1 norm for all elements of $b_{L_{W\!/\!N}}^i$ (for a particular GNSS consisting of n satellites):

$$\begin{aligned} \sum _{i=1}^{n} \Bigl |\text {frac}(b_{L_{W\!/\!N}}^i+b_0)-0.5\,\Bigr |= \text {min}\,. \end{aligned}$$

(17)

$0 \le b_0 < 1$ defines the search range to be followed. Based on the shift value for $b_0$ obtained from Eq. (17), the final estimates for the widelane/narrowlane phase biases ($b_{L_{W\!/\!N}}^i$) can be provided as

$$\begin{aligned} b_{L_{W\!/\!N}}^i := \text {frac}(b_{L_{W\!/\!N}}^i+b_0)-0.5\,. \end{aligned}$$

(18)

Resulting elements of $b_{L_{W\!/\!N}}^i$ should thus be conditioned so that $|b_{L_{W\!/\!N}}^{ij}|$ become minimal across all possible satellite pair combinations (i, j). Moreover, a particular $b_{L_{W\!/\!N}}^{ij}$ (or alternatively $B_{L_{W\!/\!N}}^{ij}$ after converting it into units of cycles) should never exceed $\pm 1$ cycle. It can be summarized that the arithmetic operations proposed in Eqs. (17) and (18) ensure to keep the phase bias values applied later as a basis as small as possible. In principle, each involved widelane/narrowlane phase bias parameter could be changed by any integer. For a possible direct comparison of phase bias products, this is not only a cosmetic matter but provides a desirable pre-alignment (especially also between successive analysis sessions).

There are some additional rules and assumptions we generally comply with:

• Phase (fractional) biases, as we consider them, are assumed to be hardware-dependent and can be interpreted as “non-integer parts of the ambiguities.”

• We generally assume that all phase observations on a frequency are aligned (independent of their phase modulation). This basic assumption is common to that made for double-difference AR.

• Although a set of distinct (widelane) phase biases is generated for each differing code-type combination, all accordingly gathered sets are finally fed into one common least-squares adjustment (using their appropriate weights, $1/\sigma ^2 \sim n$). Note that code type combinations being linearly independent should be automatically detected and finally treated independently (e. g. in case of Galileo).

• The empirical determination of widelane ($b_{L_W}^{i}$) and narrowlane ($b_{L_N}^{i}$) phase biases follows exactly the same basic principle.

• Phase biases are generated completely independent for each daily session. These are determined, in accordance with each analysis session, with a temporal resolution of 24 h.

• The presented method does not require initial information (or any starter values) regarding phase biases. These are re-estimated in each iteration.

• The CC-OSB approach is applicable to any GNSS following CDMA. Phase bias determination, as it is specified here, is carried out separately for each GNSS that shall be prepared for single-receiver IAR.

• Finally, an important feature for the results presented below concerns the underlying GNSS orbit information. This is taken unchanged from our operational double-difference analyses and is not further adjusted or enhanced (in which the ambiguities were resolved using an ensemble of complementary double-difference AR methods).

The algorithm described was implemented in the development version of the Bernese GNSS Software. The CPU time of the corresponding auxiliary tool is usually a few seconds. It is interesting to add that we could see the potential of our algorithm for detecting of quarter-cycle phase biases.

3.4 Some specific experiments

We carried out a series of specific experiments in order to validate the developed algorithm. We assessed how crucial specific processing features are, namely

• 1. iterating across narrowlane AR,

• 2. reduced clock sampling for narrowlane AR,

• 3. IRC-like narrowlane bias representation.

3.4.1 Experiment 1: Iterating across narrowlane AR

In order to assess how relevant iterations are for the IAR-relevant processing loop (cf. Fig. 2), consisting of

• clock estimation and

• phase bias determination and IAR (internally substituted by the autonomous PPP-AR procedure),

we conducted a corresponding test analysis where we performed this loop up to 10 times. We could recognize that alterations in any descriptive parameter became marginal very rapidly after few iterations. For PPP-AR, the most prominent break could generally be observed between iteration 1 (still based on ambiguity-float clocks) and iteration 2 (for the first time ever based on ambiguity-fixed clocks). Specifically, the post-fit rms of unit weight became significantly smaller for the narrowlane AR PPP solutions. The standard deviation of narrowlane ambiguity fractionals also became considerably reduced as of iteration 2, whereas the total number of resolved narrowlane ambiguities got increased by only a few percent after iterating. This increase in percentage of resolved ambiguities is actually in line with the reduced rms of unit weight afterward. In the end, by crosschecking the solution vectors of resolved ambiguity integers across iterations (in particular between iterations 1 and 2), no differences could be revealed for in common resolved ambiguities. This first experiment therefore showed that our ambiguity-float clocks of iteration “0” (while consistently applied in conjunction with the underlying ambiguity-fixed orbits adopted from a double-difference analysis and with adequately determined phase biases) already allow for a reliable narrowlane AR. In the interest of speeding up, performing the ambiguity resolution loop for an operational processing just once is quite sufficient.

3.4.2 Experiment 2: Reduced clock sampling for narrowlane AR

The processing time required for clock densification (embedded in the clock generation scheme in Fig. 2) is substantial. For this reason, we asked ourselves whether satellite clocks at intervals of 30 sec are indispensable for the PPP-AR procedure (dashed red line in Fig. 2) or whether we could proceed with a narrowlane AR done on the basis of decimated (300-sec) clocks (solid red line in Fig. 2). It is obvious that by conducting clock densification just once (instead of twice) clock analysis could be accelerated to a notable extent. Please note that widelane AR (being independent from orbit and clock information) is done only once, generally using 30-sec observation data. Primarily due to the fact that formal errors of narrowlane ambiguity parameters (and consequently corresponding confidence intervals for statistical testing associated with AR) become larger due to reduced data sampling (here by a factor of $\sqrt{10}$ compared to the “regular” case), the success rate of resolved narrowlane ambiguities becomes accordingly slightly smaller for the “decimated,” but ultimately “faster” case. Once more, the correctness of the solution vector of resolved ambiguity integers turned out to remain unaffected by this short cut. This means that satellite clocks with limited sampling are sufficient for single-receiver narrowlane AR. This allows us to follow the solid red line in Fig. 2 for any clock analysis line.

3.4.3 Experiment 3: IRC-like narrowlane bias representation

This experiment was dedicated to investigate numerical differences between a clock product following our clock/bias representation (CC-OSB) and a clock product following the integer recovery clock (IRC) strategy as realized by Laurichesse et al. (2009), Loyer et al. (2012). We could manage such an experiment easily by just disregarding the (fractional) narrowlane bias term, $c\,B_{L_N}^{ij}$, in Eq. (9) for computation of ambiguity-fixed clock parameters ($t^{i}$). In our case, this is achieved by redefining $B_{L_N}^{i}:=0$. Note that narrowlane AR was completed in the usual manner beforehand. It is obvious that resulting IRC-like satellite clock parameters thus de facto correspond to integer clocks, $t_N^{i}$, we introduced in the context of Eq. (13). It is worth mentioning that, in our case, (fractional) narrowlane bias values ($b_{L_N}^{i}$ as well as $B_{L_N}^{i}$) may be expected in our case to be minimized according to Eq. (17). Numerical differences for satellite-to-satellite clock differences, $t_N^{ij}$, between CC-OSB results (reduced according to Eq. (14)) and IRC-like results (directly computed) turned out to be a level of several femtoseconds only. This extremely high consistency confirmed us that these two methods must rely on the same basic principle and, therefore, must be considered equivalent, at least in the sense of “phase clocks.”

3.5 GNSS clock generation scheme approved for use

The results of experiment 1 we summarized in Sect. 3.4.1 showed us that we can accomplish the generation of ambiguity-fixed GNSS clock corrections and associated phase biases by complying with the following sequence of processing steps:

• 1. GNSS satellite clock estimation, called “zero floated” iteration: using phase and pseudorange code tracking data, with floating phase ambiguities;

• 2. phase bias determination and single-receiver integer fixing (first for widelane, then for narrowlane LC), using the global tracking data sample from step 1;

• 3. GNSS satellite clock estimation, called “first fixed” iteration: with mostly fixed phase ambiguities, taking into account the OSB phase bias values determined for $L_1$ and $L_2$ and the resolved phase ambiguities for $L_1$ and $L_2$ as obtained jointly in step 2;

• 4. one, or two successive clock densification steps, using phase tracking data only.

Referring to Fig. 2, we consider it sufficient to execute the IAR task loop (shown there in red) once, particularly with regard to the narrowlane LC. This procedure applies in particular if the orbital information is already available with sufficient accuracy (as appropriately assumed at the end of Sect. 3.3). Let us note for completeness that in case further iterations are desired, steps 2 (exclusively for the narrowlane LC) and 3 would have to be repeated.

3.6 Advantages of the CC-OSB method

We could verify in our experiment 3 that the CC-OSB and the IRC method are practically identical in terms of “phase clocks.” Nevertheless, we consider the following advantages of the CC-OSB method (against the well-established original IRC method) to be decisive:

• GNSS biases can be provided in a consistent way for phase and code observation data.

• There is no supplementary sign dependency, of signs which may depend on providers’ analysis conventions (such as, e. g., the WSB sign convention).

• It is possible to use the same satellite clock corrections with other carrier phase signals that are not the ones used to generate them.

• It is capable of multi-GNSS, easy to use, suitable for standardization, and incidentally relying on a GNSS bias data format (Bias-SINEX) officially approved by the IGS.

• Associated GNSS clocks can be consistently used with phase and code observation data (thus specifically allowing for code-supported narrowlane AR).

• Individual satellites, for which (ambiguity-float and thus non-integer) clocks are still provided, might be excluded from AR in clearly justifiable exceptional cases (by following a “no phase bias no AR” rule).

• CC-OSB could, in principle, serve different phase alignments (hypothetically assuming alignments being dependent on specific phase observation data types). Such a generalized phase bias information could be directly transferred to the PPP-AR user.

Due to the more comprehensive GNSS phase bias information that must be dealt with, a little more work is required, especially for the extra correction of the $L_1$ and $L_2$ (or correspondingly converted narrowlane) phase biases, but this makes the CC-OSB method much more flexible. This additional work, once the method is properly implemented in software, is likely to take little additional CPU time without making much greater demands on other computer resources. The fact that “phase clocks” or “integer clocks” are no longer directly accessible can be seen as a disadvantage. However, integer clocks can be recovered using a simple relation (cf. Eq. (14)), which, moreover, is directly applicable to undifferenced clock and phase bias information.

4 Overview of CODE clock analysis products

CODE clock products that adhere to the integer-cycle property have been submitted to the IGS since mid of 2018 (Schaer et al. 2018), namely for the following three analysis lines:

• CODE (Early and Final) Rapid, starting with 3 June 2018 (start of GPS week 2004), considering GPS/GLONASS and, as of 22 September 2019 (start of GPS week 2072), also incorporating Galileo (Dach et al. 2020b);

• CODE Final, starting with 3 June 2018 (start of GPS week 2004), considering GPS/GLONASS (Dach et al. 2020a);

• CODE MGEX (COM), starting with 17 June 2018 (start of GPS week 2006), incorporating a total of five GNSS, specifically GPS/GLONASS/Galileo/BeiDou/QZSS (Prange et al. 2020a, b).

Single-receiver ambiguity fixing could be established not only for GPS but also for Galileo. Figure 3 gives an overview of these CODE ambiguity-fixed clock analysis products together with their main characteristics.

Fig. 3

CODE’s ambiguity-fixed clock analysis products and their main characteristics: GNSS considered (black) and additionally IAR-enabled (red), clock sampling interval, elevation mask angle, number of stations jointly used in analysis, start of production and submission to the IGS (GPS week number)

IGS Rapid analysis products should be available until 17:00 UT of the next day. In our case, a first update (“Early Rapid”) may be expected already around 09:00 UT. In addition, a second update for the CODE Rapid orbit and clock products (internally called “Final Rapid”) is provided, ultimately targeting the nominal deadline (17:00 UT) for IGS Rapid products. (Dach et al. 2015b). For the latter update, the “middle day” part of a long-arc orbit analysis is considered. The long-arc orbit strategy as pursued at CODE (Beutler et al. 1996; Lutz et al. 2016) is illustrated in Fig. 4, here with a specific focus on the evolution in time.

Fig. 4

Long-arc orbit strategy pursued at CODE

In this figure, a particular Rapid orbit solution (resulting from first update: “last day” part of 3-day orbit arcs) is indicated in orange; a particular Final-Rapid orbit solution (resulting from second update: “middle day” part of 2.5-day orbit arcs) is indicated in light red. In terms of orbits, the latter update takes place automatically with every CODE Ultra-rapid analysis update done for the past fraction of the next (current) day. Figure 4 finally illustrates the situation concerning CODE Final orbits which again correspond to respective middle-day slices of 3-day long-arc orbit solutions. In contrast to CODE’s Rapid and MGEX clock products that comprise 30-second clock corrections, the CODE Final clock product (consistently generated in conjunction with such middle-day orbit slices) comprises also high-rate clocks at dense intervals of 5 seconds. It includes, in addition, phase-derived satellite clock estimates for the second midnight epoch (accordingly marked with “24:00 epoch” entries in Fig. 3). Note that this extra epoch is currently available only internally. It is not yet included in those clock product files submitted to the IGS data centers (due to the current policy, where provision of overlapping, or redundant analysis results is not accepted).

5 Results

5.1 GNSS clock/bias product validation using PPP

The repeatability of station position results into North, East, and vertical (up) components is a reliable quality measure for validation of GNSS analysis products. Figure 5 shows corresponding daily repeatability results as we could accumulate from traditional PPP, where phase ambiguities are treated floating, and finally from PPP with satellite-to-satellite ambiguity resolution (PPP-AR).

Fig. 5

Validation of CODE Final (top) and CODE MGEX (bottom) integer-property clock/bias products by means of traditional (ambiguity-float) PPP (blue) and PPP with ambiguity resolution (PPP-AR) (red), respectively. Median standard deviation of daily station coordinate determinations using single-GNSS PPP is compared here

It is important mentioning that we generally apply partial AR or, to be more precise, a sequentially and repeatedly executed AR strategy, which achieves a bootstrapping effect by invoking it multiple times. The repeatability metrics were computed on the basis of 1 month of data for two specific analysis lines (CODE Final and CODE MGEX). The upper panel of Fig. 5 shows the situation concerning the CODE Final clock/bias product, evaluating the quality for the GPS constellation with 32 active satellites in September 2018. 295 stations out of a total of 337 processed were taken into account for the GPS-only PPP performance statistics. Only those observing stations were taken into consideration that contributed frequently (at least 27 of 30 days). The lower panel finally shows the situation concerning the CODE MGEX clock/bias product, evaluating separately the quality for the Galileo constellation with 24 active satellites in August 2019. 92 stations out of 140 were taken into account here. Note that PPP/PPP-AR analyses were restricted here to Galileo observation data only to get an individual quality assessment. Similar to the case of double-difference AR, PPP single-difference AR improves principally the position accuracy in the East direction, where we could observe a significant average improvement (by almost a factor of 2). GPS-based PPP-AR is still superior. However, Galileo-based PPP-AR is already capable to provide station position results at a competitive accuracy level. The differences in performance can essentially be explained by different availability in terms of the number of satellites (32 GPS and 24 Galileo) and the number of tracking stations ( $300$ used for GPS, 140 for Galileo). The validation results presented in Fig. 5 are encouraging, making it interesting to simultaneously analyze observation data from two (or more) IAR-enabled GNSS constellations.

5.2 GNSS phase bias arithmetics needed for clock/bias product comparison

The interested reader should recognize at some point that

• widelane fractional phase bias parameters, $B_{L_W}^{ij}$,

• narrowlane fractional phase bias parameters, $B_{L_N}^{ij}$,

• integer-property clock offset parameters, $t^{ij}$,

are connected with each other by transformational links as defined by the observation equations of Sect. 2. This implies that any correction applied to any parameter of the above listing can be transferred directly to the other parameters. In order to allow us to constitute product comparisons, we will derive in the following equivalence relations between a pivot product (subsequently labeled product A) and a sample product (subsequently labeled product B).

5.3 GPS and Galileo widelane phase bias results

To enable a direct comparison of phase bias products, let us rewrite Eq. (5) for Melbourne–Wübbena (MW) observations of any receiver k in the world, with respect to two competing products, a (pivot) product A and a (sample) product B:

$$\begin{aligned} \text {MW}_k^{ij} =\,&\left( \lambda _W {N_W}_k^{ij}\!+\!c\,B_{L_W}^{ij}\!+\!c\,B_{P_W}^{ij}\right) _{\text { A}}\nonumber \\ =\,&\left( \lambda _W ({N_W}_k^{ij}\!+\!\delta {N_W^{ij}}_{\text {B}})\!+\!c\,B_{L_W}^{ij}\!+\!c\,B_{P_W}^{ij}\right) _{\text { B}}\,. \end{aligned}$$

(19)

In the thus obtained equivalence equation, the widelane ambiguity parameter ${N_W}_k^{ij}$ is suggested to be common for the two products (A and B). For the latter product (B), the widelane ambiguity is further assumed to be shifted by an arbitrary integer number, denoted with $\delta {N_W^{ij}}_{\text {B}}$, where ij indicates for a specific pair of satellites. Note that all remaining parameters included in each bracket expression (each labeled with a subscript “Product X” or, if necessary, briefly with “X”) have to be interpreted specific to a particular product. Differencing the two MW observation expressions of Eq. (19) and finally solving for $\lambda _W \delta {N_W^{ij}}_{\text {B}}$ yields

$$\begin{aligned} \lambda _W (\delta {N_W^{ij}}_{\text {B}}+r_{W}^{ij})&= \left( c\,B_{L_W}^{ij}\!+\!c\,B_{P_W}^{ij}\right) _{\text { A}}\nonumber \\&\qquad - \left( c\,B_{L_W}^{ij}\!+\!c\,B_{P_W}^{ij}\right) _{\text { B}}\,. \end{aligned}$$

(20)

$r_W^{ij}$ was introduced to tolerate a (widelane phase bias) residual error in units of widelane cycles. Such residual errors ($r_W^{ij}$) are expected to not exceeding a small fraction of a cycle. Note that the term $\lambda _W {N_W}_k^{ij}$, which was accepted to be common for the two products, cancels out after differencing Eq. (19).

5.3.1 CODE internal validation

We used Eq. (20) for investigation of the reproducibility of GPS and Galileo widelane phase bias results (generated independently for each day). We compared corresponding daily results of a specific GNSS for day i (Product A) with those for day $i+n$ (Product B), where n can be considered in this context as a time lag. Figure 6 shows the standard deviation of $r_W^{ij}$, evaluated according to Eq. (20), as a function of varying time lag (n), specifically for lags of 1, 2, 3, 4, 5, 10, 15, 30 days.

Fig. 6

Reproducibility study in respect of GPS (top) and Galileo (bottom) widelane fractional cycle bias (FCB) results produced at CODE. 0.01 cyc level is indicated with an extra line (total scale is 0.05 cyc)

For GPS, the reproducibility is evaluated at a level of about 0.005 cyc for small time lags. It is finally increasing with increasing time lag (to around 0.03 cyc for a lag of 30 days). For Galileo, the reproducibility level is a little bit higher compared to that of GPS (at about 0.01 cyc), but almost independent of the time lag. This generally confirms that Galileo widelane phase biases are highly stable in time. However, this must also be put into perspective right away because phase bias discontinuity events seem to occur much more frequently for Galileo. Whereas, for GPS, the reproducibility as represented by the red curve (1 day lag) usually gets minimal, the light blue curve (associated with a lag of 10 days) tends to be lowest in case of Galileo. This is an interesting detail since the Galileo satellites repeat their ground track patterns every 10 days. Because the corresponding period is 1 day for GPS, this might be a contribution to the better performance of daily estimated biases for GPS.

5.3.2 Verification against external CNES/CLS WSB tabular data

We verified our GPS widelane phase bias results taking well-established CNES/CLS WSB (widelane satellite bias) tabular data into account (Laurichesse et al. 2009; Loyer et al. 2012). In order to do so, we had to bring their daily tabular values in a form being compatible to our formalism. Let us allocate a correspondingly expanded WSB bias representation to the second bracket expression of Eq. (20):

$$\begin{aligned} \lambda _W (\delta {N_W^{ij}}_{\text {B}}\!+\!r_{W}^{ij})&= \left( c\,B_{L_W}^{ij}\!+\!c\,B_{P_W}^{ij}\right) _{\text { A}}\nonumber \\&\qquad -\Big (\lambda _W (-\text {WSB}^{ij})\!+\!c\,\Delta B_{P_W}^{ij}\Big )_{\text { B}}\,. \end{aligned}$$

(21)

For the CNES/CLS GPS WSB product, the code-locked and thus DCB-dependent bias term, $\Delta B_{P_W}^{ij}$, is supposed to vanish ($\Delta B_{P_W}^{ij}=0$) in case of GPS C1W and C2W observation data. Accordingly, GPS C1W−C1C DCB values should be used in case of GPS C1C and C2W: $\Delta B_{P_W}^{ij} = \kappa _{P_1} (-\text {DCB}^{ij})$. Residual fractional parts ($r_W^{ij}$) between two synchronized sets of widelane phase biases can be evaluated as

$$\begin{aligned} r_{W}^{ij}&= \text {frac}\left( \frac{c}{\lambda _W}\left( B_{L_W}^{ij} + B_{P_W}^{ij}\right) \right) _{\text { A}}\nonumber \\&\qquad - \text {frac}\left( -\text {WSB}^{ij} + \frac{c}{\lambda _W} \Delta B_{P_W}^{ij}\right) _{\text { B}}\,, \end{aligned}$$

(22)

where $B_{L_W}^{ij} + B_{P_W}^{ij}= \kappa _{L_1} B_{L_1}^{ij}\!+\!\kappa _{L_2} B_{L_2}^{ij}\!+\!\kappa _{P_1} B_{P_1}^{ij}\!+\!\kappa _{P_2} B_{P_2}^{ij}$, finally reminding us of Eq. (5).

A comparison of CODE’s GPS widelane phase bias values with CNES/CLS GPS WSB daily tabular values confirmed an agreement generally at a level of 0.01 widelane cycles (or 0.03 ns) when referring to GPS C1W and C2W observation data (specifically for the CODE Final bias product).

5.4 Integer-property clock results

In order to constitute comparisons of integer-property clock products, equivalence relations are finally required for narrowlane observations of a particular GNSS. Let us take over the widelane ambiguity parameter (${N_W}_k^{ij}+\delta {N_W^{ij}}_{\text {B}}$) for the (sample) product B as it was introduced in Sect. 5.3, specifically in Eqs. (19) and (20). The following equivalence equation may be formulated based on Eq. (8) for hypothetical narrowlane observations as collected from a virtual receiver (R), again for a (pivot) product A and a (sample) product B:

$$\begin{aligned} {L_0}_k^{ij} =&\left( \rho _k^{ij}\!-\!c\,t^{ij} \!+\!\lambda _N ({N_1}_k^{ij}\!-\!\kappa _{L_2} {N_W}_k^{ij}) \!+\!c\, B_{L_N}^{ij}\right) _{\text { A}}\nonumber \\ =&\left( \rho _k^{ij}\!-\!c\,t^{ij} \!+\!\lambda _N ({N_1}_k^{ij}\!+\!\delta {N_1^{ij}}_{\text {B}}\right. \nonumber \\&\left. \!-\!\kappa _{L_2} ({N_W}_k^{ij}\!+\!\delta {N_W^{ij}}_{\text {B}})) \!+\!c\, B_{L_N}^{ij}\right) _{\text { B}}\,. \end{aligned}$$

(23)

Differencing the two narrowlane observation expressions of Eq. (23) and finally solving for $\lambda _N \delta {N_1^{ij}}_{\text {B}}$ yields

$$\begin{aligned} \lambda _N \delta {N_1^{ij}}_{\text {B}}+cR_{N}^{ij}&=\left( \rho _k^{ij}\!-\!c\,t^{ij}\right. \left. +c\,B_{L_N}^{ij}\right) _{\text { A}}\nonumber \\&\quad - \left( \rho _k^{ij}\!-\!c\,t^{ij}\right. -\lambda _N \kappa _{L_2} \delta {N_W^{ij}}_{\text {B}} \left. +c\,B_{L_N}^{ij}\right) _{\text { B}}&\,.\nonumber \\ \end{aligned}$$

(24)

$R_N^{ij}$ was introduced to describe a (clock) residual error in units of time. The prefactor $\kappa _{L_2}=-\frac{\nu _2}{\nu _1+\nu _2}$ was already introduced in Eq. (5). Note that the terms $\lambda _N {N_1}_k^{ij}$ and $-\lambda _N\kappa _{L_2}{N_W}_k^{ij}$, which were accepted to be common for the two products, cancel out after differencing Eq. (23).

An essential property of GNSS clock parameters ($t^{i}$ and corresponding single differences $t^{ij}$) is that supplementary consideration of pseudorange code measurements using Eq. (11) should strongly dictate their quasi-absolute alignment (at a few-cm, or sub-$\lambda _N$ accuracy level). This means that a self-adjusting integer-cycle clock compensation going with “$\text {round}(\lambda _N\kappa _{L_2} \delta {N_W^{ij}}_{\text {B}})$” must be presupposed. As a result of this “automatic” clock compensation, $\delta {N_1^{ij}}_{\text {B}}$ of Eq. (24) is in practice at most a single-cycle flip, i. e. $\delta {N_1^{ij}}_{\text {B}}=\{-1,0,+1\}$.

Ideally, widelane AR was done consistently for both products (A and B). If this is the case, the additional bias component ($-\lambda _N\kappa _{L_2} \delta N_{W_B}^{ij}$) vanishes from Eq. (24) (as $\delta N_{W_B}^{ij}=0$). Another clarification of Eq. (24) can be achieved by using integer phase clock parameters ($t_N^{ij}$) we introduced in Eqs. (13) and (14). The resulting equivalence relation then simply reads as

$$\begin{aligned} \lambda _N \delta {N_1^{ij}}_{\text {B}}\!+\!c\,R_{N}^{ij}&= \left( \rho _k^{ij}\!-\!c\,t_N^{ij} \right) _{\text { Product A}}\nonumber \\&\qquad - \left( \rho _k^{ij}\!-\!c\,t_N^{ij} \right) _{\text { Product B}}\,. \end{aligned}$$

(25)

This highly simplified equation clearly illustrates the mechanism that will play when comparing two integer-property clock products, or, to be more specific, when comparing a satellite-to-satellite clock difference of product A with a satellite-to-satellite clock difference of product B. Eq. (25) just includes an integer ambiguity ($\delta {N_1^{ij}}_{\text {B}}$), a residual (clock) error ($R_N^{ij}$), two satellite-to-satellite clock differences ($t_N^{ij}$), and finally two satellite-to-satellite geocentric distance differences ($\rho _k^{ij}$).

For PPP-AR on the one hand, consistency between $\rho _k^{ij}$ and $t_N^{ij}$ is absolutely important. For computation of residual clock errors ($R_N^{ij}$) on the other hand, an appropriate so-called radial orbit correction should be made in any case:

$$\begin{aligned} \delta \rho _k^{ij}= \left( \rho _k^{ij}\right) _{\text { Product A}}- \left( \rho _k^{ij}\right) _{\text { Product B}}\,. \end{aligned}$$

(26)

Note that $\rho _k^{i}$ corresponds to the length of the geocentric position vector pointing to satellite i and $\rho _k^{ij}$ describes a corresponding single difference ($\rho _k^{i}-\rho _k^{j}$).

5.4.1 Comparison of integer clocks among various CODE analysis lines

For comparison of integer (or phase) satellite clocks among various analysis lines, corresponding clock products of two selected analysis lines are assigned to product A and product B. We have been able to make appropriate comparisons for GPS clock results successfully among all three IAR-enabled clock analysis lines from CODE, specifically CODE Final versus CODE Rapid, CODE Final versus CODE MGEX, CODE Rapid versus CODE MGEX. Exemplary comparison results are shown in Fig. 7 for two specific comparisons between clock products stemming from CODE Final and CODE (Final-) Rapid.

Fig. 7

Standard deviation of GPS satellite clock differences between CODE (Final-) Rapid and CODE Final clock/bias analysis products, computed for daily sessions of 2 June 2018 (ambiguity-float) and 13 September 2018 (ambiguity-fixed results)

Satellite-to-satellite clock differences ($R_N^{ij}$) were computed according to Eq. (24) we established in Sect. 5.4. Radial orbit corrections were applied according to Eq. (26). Both panels of Fig. 7 show resulting $R_N^{ij}$ every 300 seconds for a particular daily session, color-coded for each involved GPS satellite (given with PRN number). The plotting scale for Fig. 7 was chosen according to $\pm 1$ narrowlane cycles ($\pm \lambda _N/c\approx \pm 356.8$ ps).

The lower panel of Fig. 7 shows the situation for a pair of ambiguity-fixed (or integer-property) clock products as generated for 13 September 2018 (day 256 of 2018). It should be mentioned that integer-cycle corrections were necessary for both $\delta {N_1^{ij}}_{\text {B}}$ and $\delta {N_W^{ij}}_{\text {B}}$ terms of Eq. (24). Clock differences for one satellite, PRN/SVN G04/G049, turned out to be unreliable. This particular satellite was marked unusable/unhealthy by the GPS system operators during that time. As a consequence of this, it was poorly/sparsely observed by the IGS tracking network. Therefore, we disregarded G04 for computation of overall alignment and statistics. For this daily comparison (between CODE Final and CODE Rapid), the standard deviation (STD) of $R_N^{ij}$ is 7.3 ps (for day 256 of 2018). This value is also representative for other days. If radial orbit corrections are not taken into account, the STD value reaches 23.8 ps.

The upper panel of Fig. 7 was added for the purpose of direct comparison—in order to show the benefit of our method (prior to activation of IAR on 3 June 2018) for a pair of ambiguity-float clock products as generated the last time for 2 June 2018 (day 153 of 2018). In this special case, the $\delta {N_1^{ij}}_{\text {B}}$ term of Eq. (24) finally was “misused” for compensation of a real-valued correction specific to each satellite involved (further setting $\delta {N_W^{ij}}_{\text {B}}=0$).

5.4.2 Regular comparisons and combinations made by IGS analysis center coordinator

IGS orbit and clock product comparisons and combinations are made by the IGS analysis center coordinator (ACC) on a regular basis (http://acc.igs.org). Figure 8 shows snapshots of corresponding comparison plots for the Rapid (top) and Final (bottom) analysis line of the IGS.

Fig. 8

IGS clock comparison for Rapid (top) and Final (bottom) analysis line, specifically the Std Dev statistics that is appropriate for (phase-based) PPP applications (courtesy of http://acc.igs.org)

These comparison plots specifically show the Std Dev statistics that is appropriate for (phase-based) PPP applications. Std Dev values are declared to be computed by removing a separate bias for each satellite/station clock. Note that is exactly what we did in Sect. 5.4.1 for generation of the upper panel of Fig. 7.

The point in time at which we have activated single-receiver IAR in both clock analyses (specifically GPS week 2004) is marked with an arrow. A very significant drop in Std Dev could be noticed for the CODE Rapid as well as the CODE Final clock product. Note that the conspicuous drop with respect to JPL’s Final clock product was by coincidence (as it was again included in the IGS combination starting with wk 2004). In our case, the reduction in Std Dev can be roughly quantified by 50%. Interesting is that a similar ratio can be seen in the designated formal errors of clock parameters between an ambiguity-fixed and an ambiguity-float solution (“fixed” divided by “float” formal errors). In late 2019 and early 2020, comparison/combination of IGS legacy clock products made by the IGS ACC manifests a standard deviation of about 7–8 ps STD for the CODE Rapid and about 6–7 ps STD for the CODE Final clock product.

Figure 9 shows a snapshot of a comparison plot where no separate (satellite-specific) biases are removed.

Fig. 9

IGS clock comparison for Rapid analysis line, specifically the RMS statistics that is suitable for pseudorange-based applications or when accessing the satellite time scale (courtesy of http://acc.igs.org)

This kind of comparison, called RMS statistics, is suitable for pseudorange-based applications or when accessing the satellite time scale. The RMS value for CODE clock products is typically around 75–100 ps STD (or 2–3 cm STD in range). Essential for us is that we are able to condition our clock products in a way so that both quality assessments are well fulfilled, specifically that one related to carrier phase measurements (as assessed, e. g., by Fig. 8) as well as that one related to pseudorange code measurements (as assessed, e. g., by Fig. 9).

5.4.3 Comparison of integer clocks at day boundaries

For comparison of integer (or phase) satellite clocks at day boundaries, corresponding clock results of a particular analysis line that contain a common midnight epoch are assigned to product A (day i) and product B (day $i+1$). It should be obvious that the same principles now apply as they were applied for the comparison between two analysis lines as described previously in Sect. 5.4.1. For comparing satellite clocks, an AR scheme effectively consisting of

• 1. widelane AR with respect to $\delta {N_W^{ij}}_{\text {B}}$ of Eq. (20), taking into account GNSS phase and code biases, and

• 2. narrowlane AR with respect to $\delta {N_1^{ij}}_{\text {B}}$ of Eq. (24), taking into account GNSS clocks and phase biases,

is necessary. GNSS bias information is taken from Bias-SINEX, GNSS clock information from Clock-RINEX. In order to do so, GNSS clock information generated for consecutive days must cover at least one common epoch.

5.4.4 GPS clocks

The CODE Final clock product (see Fig. 3) contains GPS clock corrections not only for the first midnight epoch (00:00:00) but also for the second midnight epoch (24:00:00). Figure 10 shows day-boundary comparison results for a total of common midnight epochs covered by GPS week 2018.

Fig. 10

Standard deviation of GPS satellite clock differences (or misclosures) at day boundaries for the CODE IGS Final clock/bias analysis product, computed for the 1-week period of 9–15 September 2018 (GPS week 2018)

Once again, satellite-to-satellite clock differences ($R_N^{ij}$) were computed according to Eq. (24) and radial orbit corrections were applied according to Eq. (26). The plotting scale for Fig. 10 was chosen according to $\pm 1$ narrowlane cycles ($\pm \lambda _N/c\approx \pm 356.8$ ps).

The lower panel of Fig. 10 shows the situation for integer-recovered (or phase) GPS clocks ($t_N^{ij}$) as obtained from CODE Final clock/bias solutions for 9–15 September 2018 (days 252–258, 2018). It should be mentioned that integer-cycle corrections were necessary only for the $\delta {N_1^{ij}}_{\text {B}}$ term of Eq. (24). For PRN/SVN G04/G049 (marked unusable and thus sparsely observed during that time), clock differences turned out to be unreliable once more. There are rare cases where the expected consistency is not achieved for all satellites (where satellite clocks at 24:00:00 and at 00:00:00 do not match in all cases). These exceptional cases need further investigation. We disregarded clear outliers (principally due to G04) for computation of overall alignment and statistics. For this weekly comparison of the CODE Final clock/bias product at day boundaries, the standard deviation of $R_N^{ij}$ is 8.3 ps (15.3 ps with disregarded radial orbit corrections).

The upper panel of Fig. 10 shows misclosures of raw (or nominal C1W/C2W pseudorange) GPS clocks ($t^{ij}$) as obtained from CODE Final clock solutions for direct comparison. Day-boundary misclosures for raw (pseudorange-leveled) GPS clocks are typically at a level better than 100 ps STD. By the way, this particular quality assessment of CODE’s GPS (raw) clocks is rather consistent to those results suggested by the IGS ACC in Fig. 9. The observed clock misclosures never exceeded $\pm 200$ ps (approximately half of a narrowlane cycle) for GPS satellites marked usable.

5.4.5 Galileo clocks

Figure 11 illustrates the scheme intended for linear extrapolation of Galileo clocks to a common epoch.

Fig. 11

Galileo clock extrapolation scheme for comparison at day boundaries

The proposed linear extrapolation is made only on the basis of the two closest epochs, either an extrapolation (of day i) forward to 24:00:00 (left) or an extrapolation (of day $i+1$) backward to 23:59:30 (right). In accordance with Fig. 11, the “overlap clocks” for satellite i can be derived as

$$\begin{aligned} \begin{aligned} t_{(24:00:00)}^{i}&= 2\,t_{(23:59:30)}^{i} \!-\! t_{(23:59:00)}^{i}\\ t_{(-0:00:30)}^{i}&= 2\,t_{(00:00:00)}^{i} \!-\! t_{(00:00:30)}^{i}\,. \end{aligned} \end{aligned}$$

(27)

Galileo clock estimates that originate from the CODE MGEX clock analysis, finally predicted according to Eq. (27), revealed satellite-to-satellite clock misclosures at a level of 12 ps STD (for extrapolation forward as well as backward), which also attests the high stability of Galileo passive hydrogen clock masers. This also means that (observed) clocks for the second midnight epoch are actually not needed for Galileo satellites for the purpose of determination of integer narrowlane-style ambiguities $\delta {N_1^{ij}}_{\text {B}}$ in Eq. (24).

5.4.6 GPS and Galileo clock misclosure levels in comparison

Let us briefly compare the clock misclosure levels we obtained at 7–8 ps STD for GPS and at around 12 ps STD for Galileo. How is the difference in performance to be classified? For the following consideration, let us generally assume that satellite clock estimates ($t^{i}$) are all of a certain accuracy, $\sigma _{t}$, and further that they are uncorrelated in time. In case of GPS, where we had “overlap clocks” available, clock differences (for satellite i) at day boundaries could be calculated simply as

$$\begin{aligned} \begin{aligned} \delta t_{(24:00:00)}^{i}&= t_{(24:00:00)}^{i} \!-\! t_{(00:00:00)}^{i}\,. \end{aligned} \end{aligned}$$

(28)

In case of Galileo, where a clock prediction had to be made, we calculated corresponding clock differences as

$$\begin{aligned} \begin{aligned} \delta t_{(00:00:00)}^{i}&= \left( 2\,t_{(23:59:30)}^{i} \!-\! t_{(23:59:00)}^{i}\right) \!-\! t_{(00:00:00)}^{i}\\ \delta t_{(23:59:30)}^{i}&= t_{(23:59:30)}^{i} \!-\! \left( 2\,t_{(00:00:00)}^{i} \!-\! t_{(00:00:30)}^{i}\right) \,. \end{aligned} \end{aligned}$$

(29)

Error propagation finally tells us that the expected accuracy of satellite clock differences may be estimated as $\sigma _{\delta t}=\sqrt{2}\,\sigma _{t}$ in the first case (GPS) and as $\sigma _{\delta t}=\sqrt{6}\,\sigma _{t}$ in the second case (Galileo). It is therefore remarkable that, conversely, $\sigma _{t}$ can be attested to be at the same level of roughly 5 ps STD, namely for both GPS and Galileo! This is all the more astonishing because in the case of Galileo one of the clock estimates was predicted (therefore to be expected noisier) and not observed. Short-term clock prediction—as illustrated in Fig. 11 and described in Eq. (27)—therefore seems to work almost perfectly for Galileo clocks. This in turn confirms the great performance of the Galileo satellite clocks.

5.5 Generation of continuous integer-property clock information over 72 h

The availability of consistent and continuous orbit and clock analysis products that are covering session lengths of 30 h or longer is a request that has been expressed for a long time from the analysis community of LEO-originated GNSS data (Bock et al. 2007). We could demonstrate in Sect. 5.4.3 that satellite clock misclosures at day boundaries (after resolving of potential narrowlane ambiguity flips) were obtained at a noise level that is consistent to that one attested for CODE satellite clock determinations. Such a connection of clock and bias analysis products can be done not only for day i and day $i+1$ but, on top of that, also for a third day (day $i-1$). As a result of this, the CODE Final analysis products actually allow to straightforwardly generate continuous integer clock information over 72 h (or, in principle, over a multiple of 24 h) on the basis of the clock and phase bias information of three consecutive daily sessions. This outstanding capability is also supported by the long-arc orbit modeling strategy that is pursued at CODE (Beutler et al. 1996; Lutz et al. 2016) since the beginning of the IGS service (see also Fig. 4).

The following three rules should be followed in order to get a consistent and continuous orbit/clock/bias product representation when considering a series of corresponding daily products:

• 1. Validity period of each individual daily product shall be from 00:00:00 until 23:59:59 (or closing epoch).

• 2. Information referring to the second midnight epoch (24:00:00), if available, is for interpolation within a particular day or ultimately for connecting of specific products (or alternatively for corresponding product validation purposes).

• 3. Clock (C) and bias (B) products of the two added days ($i-1$ and $i+1$) have to be properly aligned.

Let us discuss in the following in particular the latter rule. Resolution of integer narrowlane-style ambiguities $\delta {N_1^{ij}}_{\text {B}}$ in Eq. (24) seems almost trivial in view of (clock-induced) residual fractional parts of the order of only 0.02 narrowlane cycles.

5.5.1 Clock product realignment (strategy “C”)

We would like to point out two possible strategies for clock or bias product concatenation and realignment. To make our considerations as simple as possible, changes in the satellite constellation shall be prohibited. Strategy “C” denotes a realignment of the clock product. Here, the complete bias information of the middle day (i) is adopted unchanged for the previous day ($i-1$) and the subsequent day ($i+1$). This can be achieved just by extending the validity period for the middle-day (day i) bias product. All satellite clock values ($t^{i}$) for day $i-1$ and day $i+1$ are to be corrected as

$$\begin{aligned} t^{ij} := t^{ij} - \frac{\lambda _N}{c} \delta {N_1^{ij}}_{\text {B}}\,, \end{aligned}$$

(30)

or, written in undifferenced notation, as

$$\begin{aligned} t^{i} := t^{i} - \frac{\lambda _N}{c} \delta {N_1^{i}}_{\text {B}}\,, \end{aligned}$$

(31)

where $\lambda _N/c \approx 356.8$ ps in case of GPS. Eqs. (30) and (31) are obviously also valid for integer satellite clocks, $t_N^{i}$ (see Eq. (25)). Strictly speaking, clock corrections induced by integer widelane ambiguity flips ($\delta {N_W^{ij}}_{\text {B}}$) could be appropriately handled as well (cf. Eq. (24)). This is why strategy “C” has to be considered as a more general one (compared to strategy “B”).

5.5.2 Bias product realignment (strategy “B”)

Strategy “B” denotes a realignment of the bias product. This strategy is only reasonably applicable if no widelane ambiguity flips ($\delta {N_W^{ij}}_{\text {B}}=0$) are present. However, it is rather smart as just incremental corrections ($\delta B_{L_N}^{i}$) are necessary responding to narrowlane phase biases ($B_{L_N}^{i}$). The bias realignment we propose can be split into two computation steps. First, incremental corrections conforming to nominal phase data ($L_1$ and $L_2$) have to be computed using Eq. (16):

$$\begin{aligned} \left[ \begin{array}{l}\delta B_{L_1}^{i}\\ \delta B_{L_2}^{i}\end{array}\right] = \left[ \begin{array}{c}+\frac{\nu _1+\nu _2}{\nu _1}\\ +\frac{\nu _1+\nu _2}{\nu _2}\end{array}\right] \delta B_{L_N}^{i}\,, \text {with}~\delta B_{L_N}^{i} \!=\! \frac{\lambda _N}{c} \delta {N_1^{i}}_{\text {B}}\,. \end{aligned}$$

(32)

Secondly, the resulting incremental phase bias corrections ($\delta B_{L_1}^{i}$ and $\delta B_{L_2}^{i}$) have to be applied for realignment (redefinition) of $B_{L_1}^{i}$ and $B_{L_2}^{i}$, individually for all satellites i of a particular GNSS:

$$\begin{aligned} \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] := \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] + \left[ \begin{array}{l}\delta B_{L_1}^{i}\\ \delta B_{L_2}^{i}\end{array}\right] \,. \end{aligned}$$

(33)

By following strategy “B,” just a comparably small number of phase bias records is to be modified, specifically in Bias-SINEX daily files referring to day $i-1$ and day $i+1$, respectively. Clock product files, which are considerably bigger in size, would remain unchanged for strategy “B.”

We consider it a possible solution that in future we will be able to provide correspondingly supplemented Bias-SINEX files (still labeled with day i), which could contain consistently realigned bias information from the previous day (day $i-1$) and the following day (day $i+1$). In principle, such a supplemented Bias-SINEX file version could be regarded as a replacement for the current one (of pure 24 h files).

5.6 Densification of clocks

We have noticed in Sect. 3.2 that the actual GNSS satellite clock estimation (Fig. 2) is done at CODE—as for many other analysis centers—for principal epochs every 5 minutes (from 00:00 up to 23:55). Assuming that the geometrical, in particular orbital, information is given, an efficient algorithm developed at CODE comes into play where epoch-differenced phase observations of an arbitrarily selectable, global tracking (sub)network are analyzed for clock densification. Simultaneous consideration of corresponding phase epoch-differences allows on the one hand a reliable phase data screening and finally a phase-specific determination of clock epoch-differences, which are then fed to an interpolation between given “anchor clocks.” In this way, the mentioned algorithm allows to generate high-rate clock corrections within reasonably short time and with sufficient accuracy (Bock et al. 2009). Figure 12 illustrates two stages of clock densification as considered at CODE (in detail for the three last 5-minute intervals of a daily session).

Fig. 12

Clock densification scheme for two stages of densification, specifically for 300 sec to 30 sec (generally accomplished at CODE) and for 30 sec to 5 sec (only performed as part the CODE Final clock analysis)

The first stage of clock densification from 300 sec to 30 sec (subsequently covering a period of 00:00:00–23:59:30) is accomplished in each CODE clock analysis line. The second stage from 30 sec to 5 sec (subsequently covering a period 00:00:00–24:00:00) is only performed as part of the CODE Final clock analysis. The ultra-high-rate (5-sec) clock densification can be realized by using 1-Hz observation data as collected from real-time data streams (and finally decimated to 5-sec spacing data). Doubts whether such densified clocks would meet the high demands on the clocks are entirely justified because forming of temporal phase differences removes any phase ambiguity and the correlations of these newly formed observations are also neglected. Considering this background, we have carried out three different validations to find out how crucial densification of clocks (based on time differences of high-rate phase data) is for an ambiguity-fixed clock product.

5.6.1 Validation by comparing high-rate clocks of two analysis lines

Figure 13 shows a corresponding validation by comparing clock results from two analysis lines, namely CODE Final versus CODE (Final-) Rapid.

Fig. 13

1-day validation of densified clocks by comparing CODE Final versus CODE (Final-) Rapid clock results as computed for day 256, 2018. The standard deviation of GPS satellite clock differences resulting for principal clocks (every 300 s) is indicated in red and for densified clocks (every 30 s) in blue. The average quality level is indicated with a dotted line

We should mention that results of a similar comparison were already presented in Sect. 5.4.1, specifically in the lower panel of Fig. 7. However, the results shown in Fig. 13 now include not only “anchor” clocks from a rigorous clock estimation as computed for principal epochs every 300 s but also densified clocks every 30 s. Standard deviation values were computed accumulated for 10 specific bins, each shifted by a multiple of 30 seconds (as additionally indicated at the bottom of Fig. 12). There is no sign of deterioration for densified clocks (blue) compared to “anchor” clocks (red). This validation may seem a bit vague because the two analysis lines may be similar with regard to the procedure used and the data processed.

5.6.2 Validation by comparing clocks at day boundaries

Having a closer look at Fig. 12, it should become obvious that clock densification turns out to be generally an interpolation for all intervals from 1 up to $n-1$, with the exception of the last interval n, which turns out to be an extrapolation. Corresponding comparisons of clocks at day boundaries, which turn out to be along the way also an indirect validation of clock densification, were already made in Sect. 5.4.3, specifically in the lower panel of Fig. 10. The error-compliant consistency of the clocks at the day boundaries we could obtain there, therefore, also allows the conclusion that accordingly “extrapolated” (and thus also accordingly “interpolated”) satellite clock corrections do not suffer any loss of accuracy (at least over densification intervals not exceeding 5 minutes).

5.6.3 Validation using GRACE K-band ranging data

A direct comparison with GRACE K-band ranging (KBR) measurements was aimed for to assess the quality of the high-rate clocks of the CODE Final product. To achieve this, a reprocessing was accomplished on the basis of an adapted processing procedure that is fully conform to that of the operational CODE Final clock analysis. This dedicated reprocessing effort, made for a complete month (April) of 2007, has also shown that our advanced clock estimation procedure is well capable to cope with older data.

The resulting 1-month validation using GRACE KBR data is shown in Fig. 14.

Fig. 14

1-month validation of densified high-rate clocks using GRACE K-band ranging (KBR) data, computed for April (day 091–120) 2007 on the basis of a reprocessed time series of an integer-property clock/bias product conform to the CODE Final analysis line. GRACE KBR residuals’ results resulting from PPP kinematic solutions treated ambiguity-float (top) and ambiguity-fixed (bottom) are compared concurrently. The standard deviation resulting for principal clocks (every 300 s) is indicated in red and for densified clocks (every 10 s) in blue. The average quality level is indicated with a dotted line for each of the two PPP solution modes

The underlying positions for GRACE-A and GRACE-B were determined in purely kinematic mode (i. e., independently for each position epoch). Standard deviation values were computed accumulated (over 1 month of GRACE KBR data) for 30 specific bins, each shifted by a multiple of 10 seconds. The upper panel of Fig. 14 shows the situation as it appears for PPP solutions treated ambiguity-float, the lower panel finally shows the situation as it appears for PPP-AR solutions (with fixed ambiguities). The gain in GRACE KBR residuals’ data scattering after ambiguity fixing, from about 18.2 mm down to 3.8 mm STD, is comparably pronounced for the LEO-based kinematic PPP analysis scenario as it is investigated here.

Based on a dedicated analysis of GRACE KBR data, as shown in Fig. 14, the quality of accordingly densified clock corrections may be assessed to be on a level similar to that of “anchor” clock corrections (sampled every 5 minutes). Note that “humps” appear in Figs. 13 and 14 in the event of a deficiency. However, such humps are not visible. Another conclusion we may draw from Sect. 5.6 is that the use of any high-rate (preferably ambiguity-fixed) clock product for densification of another poorly sampled integer-property clock product must be seen as a feasible solution. This applies to densification of GPS clocks within 5-minute intervals, but especially for high-rate densification (within 30-second intervals). Furthermore, we have learned in Sect. 5.4.3 that such a densification is not needed for Galileo clocks (data at intervals of 30 seconds provided).

6 Summary and conclusions

6.1 Method

We have introduced a new, innovative approach for describing satellite bias parameters for PPP with ambiguity resolution (PPP-AR). The existing common clock approach, in particular “CC-1” as introduced in Teunissen and Khodabandeh (2015), has been extended in a way that OSBs (observable-specific signal biases) allowing for flexible signal and frequency handling between multiple GNSS becomes possible. Strictly speaking, the OSB approach, which was initially applied exclusively to pseudorange measurements (Villiger et al. 2019), has been appropriately generalized to also include carrier phase measurements with resolved ambiguities. Our experiences and especially those of first test users show that the usage of a consistently prepared OSB product should be straightforward for OSB-based PPP-AR. This is particularly useful as both phase and pseudorange bias values can be delivered in a consistent manner—and ultimately to multiple PPP-AR-enabled GNSS concurrently. Another key benefit is that it is suitable for standardization. The exchange of OSB values is incidentally also supported by a GNSS bias data format (Bias-SINEX) officially approved by the IGS.

With this contribution, we are providing a thorough description of the CC-OSB clock/bias methodology and guidance on how the related OSB-based PPP-AR application works. Against this background, we also have addressed the sign issue with regard to associated (satellite) OSB values. Essential is that the same sign convention should hold for bias values related to pseudorange measurements and for bias values related to phase measurements (actually independently of whether they are specific to satellites or specific to receivers). The set of undifferenced observation equations we initially specified in Eq. (1) for an arbitrarily chosen set of dual-frequency GNSS observations ($L_1$, $L_2$, $P_1$, $P_2$) represents the basis for the entire mathematical formalism. For a PPP-AR user, Eq. (1) is, in essence, the only relevant formula, as it unambiguously tells how the required bias information (as suggested directly for the nominal observations, i. e., $L_1$, $L_2$, $P_1$, $P_2$) has to be applied. The basic assumption that all phase observations on a frequency are aligned (independent of their phase modulation) is common to that made for double-difference AR. It should be noted that this specific assumption is dictated by the PPP-AR product provider to the PPP-AR user (via accordingly conditioned phase bias information).

One of the numerical analysis experiments we carried out (cf. Sect. 3.4.3) was particularly important to us. We were able to verify that the proposed CC-OSB and the well-established integer-recovery clock (IRC) method are de facto identical in terms of “integer phase clocks.” In the appendix (Sect. A.3), we will use the IRC method as an example to show that an OSB-conform bias description could be reached for any other IGS clock/bias analysis product that is relying on a “common clock” assumption.

6.2 Implementation at CODE

The new, extended PPP-AR approach could be implemented and thoroughly tested at CODE. CODE clock products that adhere to the integer-cycle property have been submitted to data centers of the IGS since 3 June 2018 (since 17 June 2018 in case of MGEX). All of our IGS clock analysis lines have been upgraded accordingly: (Early- and Final-) Rapid, Final, and MGEX. They cover up to five GNSS (GPS/GLONASS/alileo/BeiDou/QZSS). Single-receiver ambiguity fixing (considering satellite differencing) could be incorporated successfully not only for GPS but also for Galileo. This applies in particular to our contribution to MGEX and, after a corresponding enhancement with Galileo to a triple GNSS, also to our IGS Rapid solutions (since 22 September 2019). Galileo PPP-AR in an operational mode, as originally announced in Schaer et al. (2018), is an achievement that can be seen as another milestone in terms of IGS analysis. An interesting detail regarding Galileo AR is that we could recognize an effect of the ground track repetition period of Galileo satellites (of 10 days) in our widelane phase bias results generated daily for Galileo.

The new CODE clock products reveal a notably improved quality and, in the end, allow for PPP-AR (with up to two AR-enabled GNSS). Each clock product (Clock-RINEX) should to be used in conjunction with the accompanying phase and pseudorange bias product (Bias-SINEX) in order to achieve best possible performance. They are conditioned in such a way that maximum consistency can be ensured for ambiguity-float, ambiguity-fixed, and pseudorange-supported PPP applications. Certainly worth mentioning is the fact that we were able to show that a temporal resolution of 24 h is sufficient for phase bias parameters when relying on the CC-OSB approach. This is not entirely surprising, however, because in the IRC approach the total fraction of narrowlane phase biases ends up in the epoch-estimated satellite clock corrections, and in our approach, there is probably still a small portion which is still capable of absorbing small fluctuations in $L_1$ and $L_2$ phase biases.

CODE clock products are generally provided with a high data rate (30 seconds or higher). The CODE Final clock product covers even a higher data rate (5 seconds) and also provides data for the second midnight (24:00) epoch. Corresponding 24:00-epoch data, however, is not yet included in the CODE Final clock product files submitted to the IGS data centers. The IGS should open up the possibility to include the second midnight epoch (00:00 and 24:00) in orbit and clock product submissions. We believe that allowing this additional epoch would actually be very beneficial (specifically considering the needs of PPP-AR), and the rules for use should be simple and clear (cf. Sect. 5.5).

The integer-cycle property of between-satellite clock differences is of fundamental importance when comparing satellite clock estimates among various analysis lines, or at day boundaries. We have repeatedly shown that the integer property of CC-OSB clock and bias products can be easily recovered. The basic equation describing the relation between “common” (or CC-compatible) satellite clock corrections and “integer” (or IRC-compatible) satellite clock corrections is simple and even works with undifferenced clock and OSB bias values (cf. Eq. (14)). Both kinds of comparisons could be exploited at a very high level of consistency, specifically at 7–8 ps standard deviation for GPS. For Galileo, satellite clock estimates predicted for the second midnight epochs revealed between-satellite clock misclosures at a level of 12 ps, which also attests the high stability of Galileo passive hydrogen clock masers. Any retrieved comparison essentially indicated a standard deviation for between-satellite clocks from CODE of the order of 5 ps (1.5 mm in range). Our quality assessment is by the way in line with an STD estimate (4.6 ps) as reported in Banville et al. (2020) for a specific example day, where (partly experimental) integer-property clock and bias products from a total of six IGS analysis centers were compared and finally combined—ultimately keeping the integer nature of the resulting GPS clock and associated bias products. It should be noted that we contributed with our operational IGS Final product line to this promising IGS PPP-AR case study.

The integer-cycle property that can be recovered between the CODE Final clock and the accompanying bias product of consecutive daily sessions (using clock estimates additionally provided for the second midnight epoch) allows us to deduce GPS satellite clock and phase bias information that is consistent and continuous with respect to carrier phase observation data over two, three, or, in principle, yet more days. Linking of the clock/bias products over consecutive days is simplified by the long-arc orbit modeling strategy as pursued at CODE. Consideration of corrections of the radial orbit for the purpose of a reliable narrowlane AR (of residual cycles) is in fact no longer necessary when linking consecutive daily clock solutions from CODE.

We have presented three different validations, which all confirm that densification of clocks (based on time differences of high-rate phase data) is not harmful and thus may be approved for an ambiguity-fixed clock product. Another conclusion we could draw is that the use of any high-rate (preferably ambiguity-fixed) clock product for densification of another poorly sampled integer-property clock product must be seen as a feasible solution. This applies to densification of GPS clocks within 5-minute intervals, but especially for an additional densification within 30-second intervals. Also, we have learned that such a further densification is not needed for Galileo clocks. An interesting detail is that the short-term prediction method we used to generate 24:00-epoch Galileo clocks has finally revealed results that can be viewed as almost perfect in terms of the error budget.

A direct comparison with GRACE K-band ranging (KBR) measurements was aimed for to assess the quality of the high-rate clocks of the CODE Final product. To achieve this, a reprocessing was accomplished on the basis of an adapted processing procedure that is fully conform to that of the operational CODE Final clock analysis. This dedicated reprocessing effort, made for a complete month (April) of 2007, has also shown that our advanced clock estimation procedure is well capable to cope with older data. The gain in GRACE KBR residuals’ data scattering after ambiguity fixing, from about 18.2 mm down to 3.8 mm STD, is comparably pronounced for the LEO-based kinematic PPP analysis scenario as it was investigated here.

The CODE analysis products as described in this article are meanwhile consistently applied for LEO POD analyses at AIUB (Mao et al. 2021). Last but not least, we would like to point out that the CODE analysis products together with the accompanying bias product have been in operational use since 6 May 2020 within the Copernicus POD Services (GMV-GMESPOD-MEM-0036 v1.1), where numerous independent software packages are in use besides the Bernese GNSS software.

6.3 Outlook and further developments

The core modules for the PPP-AR enhancement at CODE were designed and implemented to cope with multiple GNSS. First tests have been carried out already in 2018 involving other GNSS, specifically BeiDou and QZSS, to verify the technical implementation. These first tests were very promising and, among other things, revealed an interesting detail regarding QZSS. At that time, QZSS consisted of only one satellite. Due to this circumstance, we were reminded that the described methods essentially also work in the case of a single-satellite constellation. This means that the indices “ij” that we used to designate a pair of satellites or, more generalized, a satellite-to-satellite difference, actually denote different satellite paths or different ambiguity arcs. On the other hand, this means that such a “single difference” may as well refer to one particular satellite only. It was reassuring for us to see that this was ensured in our test scenario with QZSS. Regarding multi-GNSS, it is important to emphasize that we consider generally one set of epoch-specific receiver clock parameters (strictly following the “CC” idea), but in combination with receiver biases (which are treated constant over 24 h) to intercept systematic differences among two or more GNSS involved (see, e. g., Villiger et al. 2019).

Appendices

Appendix

Usage of Bias-SINEX for the OSB-based PPP-AR application

Is the Bias-SINEX data format (Schaer 2018) suitable to serve the PPP-AR application? The answer is clearly yes. Phase as well as code biases required for PPP-AR may be appropriately described for any clock analysis product conforming to a common clock (CC) assumption. Let us establish Eq. (1) as a basis for the PPP-AR terminology to be used in consideration of Bias-SINEX. In view of Eq. (1), the use of the observable-specific signal bias (OSB) representation is compelling for both phase and code bias information. To have all associated bias values generally provided in units of time (typically in nanoseconds) shall be a strong recommendation on top of that. Figure 15 shows a Bias-SINEX V1.00 file excerpt with data content supporting OSB-based PPP-AR. It gives an impression on how the required OSB code and phase bias data records (shown here for three arbitrarily selected satellites of the GPS constellation) may look like. We would like to emphasize that the rules listed below actually only refer to the OSB bias description we are proposing for the PPP-AR application. These rules are therefore actually independent of a specific data exchange format.

Fig. 15

Bias-SINEX V1.00 file excerpt, illustrating OSB (observable-specific signal bias) entries for all involved pseudorange code (blue) and carrier phase (red) observable types, finally enabling PPP-AR when being considered simultaneously in conjunction with a consistent IAR-enabled clock product

1.1 PPP-AR user side

The interested PPP-AR user has to comply with the following rules:

• 1. GNSS orbits and IAR-enabled clocks (Clock-RINEX) have to be used in conjunction with accompanying biases (Bias-SINEX).

• 2. Such a (typically daily) Bias-SINEX file should contain OSB data records for

  • pseudorange code and

  • carrier phase observations

  that are consistent with the associated (ambiguity-fixed) clock product and, furthermore, have to be applied to each involved observation, specifically in the sense

  $$\begin{aligned} O_{\text { measured}} = O_{\text { computed}} + B_{\text { satellite}}\,, \end{aligned}$$

  (34)

  or alternatively, depending on the users’ preference,

  $$\begin{aligned} O_{\text { measured}} - B_{\text { satellite}} = O_{\text { computed}}\,, \end{aligned}$$

  (35)

  where O stands for observation and B for bias. OSB data records are generally labeled with 3-figure observation codes (e. g., C1W, C2W, L1W, L2W, etc. for GPS) as defined in the RINEX3 standard (Gurtner and Estey 2018). Note that only satellite bias ($B^i$) corrections are relevant for PPP-AR.

• 3. Following this convention, the users’ PPP software should finally be able to recover the integer nature of single-receiver between-satellite phase ambiguities automatically for both

  • widelane ambiguities using the MW LC and subsequently for

  • narrowlane ambiguities using the ionosphere-free LC (when reintroducing/fixing the previously resolved widelane integers).

  Ideally, narrowlane AR may be performed not only in phase-only but also consistently in pseudorange-supported mode. That’s all.

Satellites for which phase bias data is unavailable should be suspended for AR.

1.2 PPP-AR product provider side

Eq. (1) shall also be the basis for a PPP-AR product provider. He is asked to provide all required code and phase biases in a pseudo-absolute form, specifically by converting his phase biases consistently to nominal phase observations ($L_1$ and $L_2$) according to Eqs. (15) and (16). This becomes in particular important when GNSS with more than two frequencies are considered (which goes beyond the scope of this paper). Code biases for reference observation types (e. g., C1W and C2W for GPS) may be provided consistently with zero-values: $[B_{P_1}^i,B_{P_2}^i]=[0,0]$. Biases for other code observation types may be directly substituted with adequate DCB information, e. g., with $B_{P_1}^i=-\text {DCB}^i$, for GPS C1C, where $\text {DCB}^i$ represents a C1W−C1C DCB for satellite i. The provision of consistently generated orbits and (ambiguity-fixed) clocks is generally presupposed.

“No phase bias, no (PPP-)AR.” This specific rule actually allows a PPP-AR product provider to exclude selected satellites from AR and, finally, to report that (by refusal to provide individual OSB phase bias records) to PPP-AR users so that those could also omit affected satellites for their AR actions. However, this should only be permitted in clearly justifiable exceptional cases.

1.3 Special case of integer-recovery clock products

Let us provide an example of the well-established integer-recovery clock (IRC) representation on how it can be transferred into an OSB representation being conform to Bias-SINEX. CNES/CLS GPS widelane satellite biases (WSBs) can be converted according to Eq. (16), with $B_{L_N}^{i} = 0$ (due to the IRC assumption), into phase OSBs allowed in Bias-SINEX:

$$\begin{aligned} \left[ \begin{array}{l}B_{L_1}^{i}\\ B_{L_2}^{i}\end{array}\right] = \left[ \begin{array}{c}-\frac{\nu _2}{\nu _1}\\ -\frac{\nu _1}{\nu _2}\end{array}\right] B_{L_W}^{i} = \left[ \begin{array}{c}-\frac{\nu _2}{\nu _1}\\ -\frac{\nu _1}{\nu _2}\end{array}\right] \left( -\text {WSB}^{i}\right) \,. \end{aligned}$$

(36)

Note that the negative sign comes from a different widelane ambiguity definition (${N_2}_k^{i}-{N_1}_k^{i}$). Code OSBs ($B_{P_1}^i$ and $B_{P_2}^i$) can be substituted with: $B_{\text { C1W}}^i=0$, $B_{\text { C2W}}^i=0$, and accordingly $B_{\text { C1C}}^i=-B_{\text { C1W-C1C}}^i$.