Introduction

As part of the International GNSS Service (IGS), the multi-GNSS Pilot Project, known under the acronym MGEX, is responsible for the integration of all multi-GNSS related activities, aiming to provide high-precision products for the scientific community as well as to engineering applications (Montenbruck et al. 2017). Within the MGEX, different Analysis Centers (ACs) contribute with products, including satellite orbits and clocks, which are then distributed to the public. Unlike the official IGS routine, where a GPS-based combination of the ACs solutions is performed and distributed, a multi-GNSS combination is still missing. The combined official IGS solutions are GPS only and based on the routines developed by Springer and Beutler (1993) and Kouba et al. (2001). In recent years, several groups have presented different combination processes and results, pointing out the advantages of having a multi-GNSS solution as well as the limitations of the current combination. Examples of a multi-GNSS orbit combination can be found in Sakic et al. (2018, 2020), where an adaptation of the existing software is presented for the combination of GPS, GLONASS, Galileo, BeiDou, and QZSS. The adapted algorithm handles all constellations as “one big constellation.” A similar adaptation is shown in the study of Sośnica et al. (2020), where the results from the experimental orbit combination performed by the IGS ACC (analysis center coordinator) are evaluated via satellite laser ranging (SLR), and in Zajdel et al. (2023) the repro3 orbits are also evaluated using SLR. Another adaptation of the legacy software with a focus on the alignment of the orbits to the international terrestrial reference frame (ITRF) is given in Mansur et al. (2020) and lately in Mansur et al. (2022) an orbit combination using least squares variance component estimation (LSVCE) is presented. In Chen et al. (2023) preliminary results for the combination algorithm developed at Wuhan University are given, with a weighting scheme based on the orbit residuals. A study about the clock combination is presented by Banville et al. (2020), where a clock combination takes into account the different clock and bias definitions used by the ACs. The clock and bias solutions are processed jointly, yielding in the final clock corrections and a set of observable-specific signal biases (OSB). Similar strategies for combining clock offsets and biases are also given by Geng et al. (2021) and Chen et al. (2022). Zhou et al. (2022) describe an orbit and clock combination where the weighting is based on the respective orbit and clock residuals.

The methods mentioned above for a multi-GNSS satellite clock combination generally require inter-system bias (ISB) products. ISBs define the differences of the estimable receiver clocks when using different constellations, but the chosen reference ISBs in the AC processing are also mapped to the satellite clocks of the respective constellations, so ISBs have to be considered when combining satellite clocks. ISB products are currently not provided by all ACs. Building a combination scheme relying on today’s set of available MGEX products impedes a proper ISB combination. We, therefore, derive a clock combination considering only the clock products provided by the ACs, meaning that there is no direct ISB combination or usage. The combined clocks are derived as a weighted mean of the aligned ACs clock solutions, where the weights applied are the normalized inverse of the estimated variance components. The computation of the variance components is based on the LSVCE theory (Teunissen and Amiri-Simkooei 2008).

We start by explaining the basics of clock estimation which is normally used by the ACs. Further, we discuss the necessary alignments to homogenize the clock products, including the alignment to a reference AC, the alignment applied to correct the ISBs between the constellations, and the radial orbit correction applied to the ACs clock solutions. In the next section, the mathematical formulation for the LSVCE is given, which is in line with the orbit combination derived in Mansur et al. (2022), and the combination workflow. The tests with the developed combination approach and results discussion are given first using MGEX products as input and later with repro3 solutions. Summary and conclusions are given in the last section.

Basics of GNSS clock estimation

We first recall the GNSS observation equations, which allow us to understand the actual definition of the clock parameters estimated by the ACs and the required alignments. Let the code and carrier-phase observations at the receiver \(k \in \{1,\dots , K\}\), satellite \(s \in \{1,\dots , S\}\) and frequency \(f \in \{1,\dots , F\}\) be defined for a single system as follows:

$$\begin{aligned} p_{k,f}^{s} = & \rho_{k}^{s} + {\text{d}}t_{k} - {\text{d}}t^{s} + m_{k}^{s} \tau_{k} + \mu_{f} \iota_{k}^{s} + d_{k,f} - d_{,f}^{s} + \epsilon_{k,f}^{s} \\ \varphi_{k,f}^{s} = & \rho_{k}^{s} + {\text{d}}t_{k} - {\text{d}}t^{s} + m_{k}^{s} \tau_{k} - \mu_{f} \iota_{k}^{s} + \delta_{k,f} - \delta_{,f}^{s} + a_{k,f}^{s} + \varepsilon_{k,f}^{s} \\ \end{aligned}$$
(1)

with \({\rho }_{k}^{s}\) the receiver-satellite ranges, \({\text{d}}{t}_{k}\) the receiver clock offset, \({\text{d}}{t}^{s}\) the satellite clock offset, \({\tau }_{k}\) the tropospheric delay, \({m}_{k}^{s}\) the mapping function, \({\iota }_{k}^{s}\) the first order ionospheric delay on the first frequency, \({\mu }_{f}\) frequency specific multiplier, \({d}_{k,f},{d}_{,f}^{s},{\delta }_{k,f},{\delta }_{,f}^{s}\) frequency-specific constant hardware biases, and \({a}_{k,f}^{s}\) the carrier-phase integer ambiguity. Further effects, such as phase wind-up or antenna offsets and variations, are assumed to be corrected by the ACs. Even though the antenna calibrations used by the ACs are identical within the IGS, their effect on the observations depends on the modeled attitude of the satellites, which might differ among the ACs, especially for low beta angles, see Loyer et al. (2021). These modeling differences are visible in the estimated satellite clocks, but can be corrected if all ACs provide attitude information about the satellites. Within this study, we do not apply such corrections.

Considering a network processing where most of the ACs use the ionosphere-free (IF) linear combination of two frequencies, the observation equations read as follows:

$$\begin{aligned} p_{{k,{\text{IF}}}}^{s} = & \rho_{k}^{s} + {\text{d}}t_{k} - {\text{d}}t^{s} + m_{k}^{s} \tau_{k} + d_{{k,{\text{IF}}}} - d_{{,{\text{IF}}}}^{s} + \epsilon_{k,IF}^{s} \\ \varphi_{{k,{\text{IF}}}}^{s} = & \rho_{k}^{s} + {\text{d}}t_{k} - {\text{d}}t^{s} + m_{k}^{s} \tau_{k} + \delta_{{k,{\text{IF}}}} - \delta_{,IF}^{s} + a_{{k,{\text{IF}}}}^{s} + \varepsilon_{k,IF}^{s} \\ \end{aligned}$$
(2)

with \({\left(\cdot \right)}_{\text{IF}}=\frac{1}{{\mu }_{2}-{\mu }_{1}}\left({\mu }_{2}{\left(\cdot \right)}_{1}-{\mu }_{1}{\left(\cdot \right)}_{2}\right)\). The pseudoranges \({\rho }_{k}^{s}\) are usually parameterized in terms of the receiver coordinates and orbit parameters. The subsequent system of equations would be rank-deficient so that absolute clock offsets and instrumental biases cannot be separated in the least squares adjustment. In other words, we cannot decouple the clock offsets from the frequency-specific hardware biases. With the ionosphere-free combination, a common solution to overcome the rank deficiencies is to define IF clocks as follows:

$$\begin{aligned} {\text{d}}t_{{k,{\text{IF}}}} = & {\text{d}}t_{k} + d_{{k,{\text{IF}}}} \\ {\text{d}}t_{{,{\text{IF}}}}^{s} = & {\text{d}}t^{s} + d_{{,{\text{IF}}}}^{s} \\ \end{aligned}$$
(3)

and therefore, equation (2) can be rewritten as follows:

$$\begin{aligned} p_{{k,{\text{IF}}}}^{s} = & \rho_{k}^{s} + {\text{d}}t_{{k,{\text{IF}}}} - {\text{d}}t_{{,{\text{IF}}}}^{s} + m_{k}^{s} \tau_{k} + \epsilon_{k,IF}^{s} \\ \varphi_{{k,{\text{IF}}}}^{s} = & \rho_{k}^{s} + {\text{d}}t_{{k,{\text{IF}}}} - {\text{d}}t_{{,{\text{IF}}}}^{s} + m_{k}^{s} \tau_{k} + (\delta_{{k,{\text{IF}}}} - d_{{k,{\text{IF}}}} - \delta_{{,{\text{IF}}}}^{s} + d_{{,{\text{IF}}}}^{s} + a_{{k,{\text{IF}}}}^{s} ) + \varepsilon_{k,IF}^{s} \\ \end{aligned}$$
(4)

where the term in the parenthesis in the bottom line is a time-constant ambiguity parameter. If relying on uncombined observations on an arbitrary number of frequencies (Schönemann et al. 2011), the resulting rank deficiencies have to be removed in a similar way, and the same IF clock parameters can be defined. This is the case for the TU Graz solution (Strasser et al. 2019). Even with the IF clock definition, the system of equations still suffers from a rank deficiency and has no unique solution since an arbitrary shift of all satellite clocks can simply be compensated by also shifting the receiver clocks by the same amount. To solve this issue, a reference clock is chosen by the AC, for example, \({\text{d}}{t}_{1,{\text{IF}}}\), and the estimable clock parameters are:

$$\begin{aligned} {\tilde{\text{d}}}t_{{k,{\text{IF}}}} = & {\text{d}}t_{{k,{\text{IF}}}} - {\text{d}}t_{{1,{\text{IF}}}} \;\;\;k \ne 1 \\ {\tilde{\text{d}}}t_{{\text{,IF}}}^{s} = & {\text{d}}t_{{,{\text{IF}}}}^{s} - {\text{d}}t_{{1,{\text{IF}}}} \\ \end{aligned}$$
(5)

Naturally, the ACs select different reference clocks so that the clock combination cannot be performed directly, as the clocks from different ACs represent different quantities.

Analysis centers product alignment

To remove the impact of different reference clocks, two approaches can be used:

  1. 1.

    Choose a common reference station provided by all the ACs.

  2. 2.

    Estimate offset and trend regarding a reference AC, assuming that the reference clocks have a linear behavior.

In the first approach, the AC clocks can be homogenized by shifting them to a common reference clock which is provided by all the ACs. This can be done on the product level even without the knowledge of the original reference clock by subtracting the clock of the new reference station \(l\) from all station and satellite clocks:

$$\begin{aligned} \tilde{\tilde{d}}_{{k,{\text{IF}}}} = & \;{\tilde{\text{d}}}t_{{k,{\text{IF}}}} - {\tilde{\text{d}}}t_{{l,{\text{IF}}}} \\ = & \;{\text{d}}t_{{k,{\text{IF}}}} - {\text{d}}t_{{1,{\text{IF}}}} - {\text{d}}t_{{l,{\text{IF}}}} + {\text{d}}t_{{1,{\text{IF}}}} \\ = & \;{\text{d}}t_{{k,{\text{IF}}}} - {\text{d}}t_{{l,{\text{IF}}}} \\ \tilde{\tilde{d}}_{{,{\text{IF}}}}^{s} = & \;{\tilde{\text{d}}}t_{{,{\text{IF}}}}^{s} - {\tilde{\text{d}}}t_{{l,{\text{IF}}}} \\ = & \;{\text{d}}t_{{,{\text{IF}}}}^{s} - {\text{d}}t_{{1,{\text{IF}}}} - {\text{d}}t_{{l,{\text{IF}}}} + {\text{d}}t_{{1,{\text{IF}}}} \\ = & \;{\text{d}}t_{{,{\text{IF}}}}^{s} - {\text{d}}t_{{l,{\text{IF}}}} \\ \end{aligned}$$
(6)

we can see that the reference clock is now \(l\) instead of 1.

This approach can be challenging as we often do not find any common station in all individual ACs solutions. Therefore, to overcome this issue the AC solutions are aligned to a reference AC by estimating AC-wise trend and offset based on the satellite clock offsets of each AC. This procedure is only an approximate alignment that is only valid if all reference clocks used in the AC solutions behave linearly over time. As an example, Fig. 1 shows the behavior of the station clock of NRC1 for GFZ and GODE for COD, where the reference clock for GFZ is GODE and the reference clock for COD is NRC1. That is, even though we seemingly compare the clocks of two different stations, we, in fact, compare the same quantity, namely the difference between the two clocks. Note that the sign of the COD values has been switched in the figure. In both plots, we can clearly see that at least one of the two involved clocks deviates from the assumption of linear behavior, but we cannot know which AC used a nonlinear clock as a reference.

Fig. 1
figure 1

Relative receiver clock offset between NRC1 and GODE stations for GFZ and COD, after trend and offset removed. Note that the sign in the bottom plot has been switched so that the plots are compatible

Full size image

For a better understanding, the GPS satellite clocks of GFZ and COD are given in Fig. 2, where the satellite clocks are aligned to the MIT solution by removing the common offset and trend. From (5), we remember that the behavior of the chosen reference clock is also absorbed by all satellite clocks. We observe the same pattern of the GFZ satellite clocks as for the receiver clock shown in Fig. 1, while for COD, there is no relevant trend or offset in the satellite clocks. As we expect the satellite clocks to behave linearly over time, we conclude that the reference station GODE that GFZ uses does not fit a linear behavior. As a consequence, the linear alignment, i.e., removing offset and trend with respect to a reference AC, is not valid for the GFZ solution. This problem can be solved by changing the reference clock of the GFZ solution to a suitable station, as explained in (6).

Fig. 2
figure 2

Comparison between GFZ and COD satellite clocks aligned to MIT

Full size image

Multi-constellation alignment

Another required alignment is related to the ISBs. Let us recall the code and carrier-phase observations (4), but now introduce the symbol \(*\) to mark different constellations (GPS, GLONASS, Galileo, BeiDou, QZSS)

$$\begin{aligned} p_{{k,{\text{IF}}}}^{{s_{*} }} = & \rho_{k}^{{s_{*} }} + {\text{d}}t_{{k,{\text{IF}}}}^{*} - {\text{d}}t_{{,{\text{IF}}}}^{{s_{*} }} + m_{r}^{{s_{*} }} \tau_{r} + \epsilon_{k,IF}^{{s_{*} }} \\ \varphi_{{k,{\text{IF}}}}^{{s_{*} }} = & \rho_{k}^{{s_{*} }} + {\text{d}}t_{{k,{\text{IF}}}}^{*} - {\text{d}}t_{{,{\text{IF}}}}^{{s_{*} }} + m_{k}^{{s_{*} }} \tau_{k} + \left( {\delta_{{k,{\text{IF}}}}^{*} - d_{{k,{\text{IF}}}}^{*} - \delta_{{,{\text{IF}}}}^{{s_{*} }} + d_{{,{\text{IF}}}}^{{s_{*} }} + a_{{k,{\text{IF}}}}^{{s_{*} }} } \right) + \varepsilon_{k,IF}^{{s_{*} }} \\ \end{aligned}$$
(7)

with this formulation the index \({s}_{*}\) as well as the receiver hardware biases \({d}_{k,f}^{*}\) and \({\updelta }_{k,f}^{*}\) are system specific and the estimable receiver clocks \({\text{d}}{t}_{k,{\text{IF}}}^{*}\) are given as follows:

$${\text{d}}t_{{k,{\text{IF}}}}^{*} = {\text{d}}t_{k} + d_{{k,{\text{IF}}}}^{*}$$
(8)

Since hardware biases are considered to be constant over one day, there is a constant offset between the receiver clocks of different constellations, so ACs usually estimate one receiver clock per epoch referring to the reference system, which is usually GPS \({\text{d}}{t}_{k,{\text{IF}}}^{\text{G}}\), and time-constant ISB parameters \({d}_{k,{\text{IF}}}^{{\text{G}}*}\) for all other systems where:

$$\begin{aligned} {\text{d}}t_{{k,{\text{IF}}}}^{G} = & {\text{d}}t_{k} + d_{{k,{\text{IF}}}}^{G} \\ {\text{d}}t_{{k,{\text{IF}}}}^{*} = & {\text{d}}t_{{k,{\text{IF}}}}^{G} + \underbrace {{d_{{k,{\text{IF}}}}^{*} - d_{{k,{\text{IF}}}}^{G} }}_{{{\text{ISB}} d_{{k,{\text{IF}} }}^{G*} }} \quad * \ne G \\ \end{aligned}$$
(9)

Similar to the single constellation framework, the resulting system of equations would be rank-deficient, as an additional rank deficiency is introduced due to the ISBs for each constellation other than the reference system. Similar to (5), where a reference clock was introduced, we can solve this issue by defining a reference ISB for each system, which is then absorbed by all station ISBs \({d}_{k,{\text{IF}}}^{{\text{G}}*}\) of the same system. If the ISBs of the first receiver are chosen as reference, the estimable ISB parameters are given by:

$${\tilde{d}}_{k,{\text{IF}}}^{{\text{G}}*}={d}_{k,{\text{IF}}}^{{\text{G}}*}-{d}_{1,{\text{IF}}}^{{\text{G}}*} \quad k\ne 1$$
(10)

From the observation Eq. (7), we can see that this shift is also absorbed by all satellite clocks of this system, which are also shifted by the reference ISB \({d}_{1,{\text{IF}}}^{{\text{G}}*}\) in addition to the clock of the reference receiver:

$$\tilde{{\text{d}}}{t}_{,{\text{IF}}}^{{s}_{*}}={\text{d}}{t}_{,{\text{IF}}}^{{s}_{*}}-{\text{d}}{t}_{1,{\text{IF}}}^{\text{G}}-{d}_{1,{\text{IF}}}^{{\text{G}}*} \quad *\ne {\text{G}}$$
(11)

In summary, the estimable clock and ISB parameters of a multi-GNSS solution of an AC are:

$$\begin{aligned} {\tilde{\text{d}}}t_{{k,{\text{IF}}}}^{G} = & {\text{d}}t_{{k,{\text{IF}}}}^{G} - {\text{d}}t_{{1,{\text{IF}}}}^{G} \quad k \ne 1 \\ {\tilde{\text{d}}}t_{{,{\text{IF}}}}^{{s_{G} }} = & {\text{d}}t_{{,{\text{IF}}}}^{{s_{G} }} - {\text{d}}t_{{1,{\text{IF}}}}^{G} \\ \tilde{d}_{{k,{\text{IF}}}}^{G*} = & d_{{k,{\text{IF}}}}^{G*} - d_{{1,{\text{IF}}}}^{G*} \quad k \ne 1,* \ne G \\ {\tilde{\text{d}}}t_{{,{\text{IF}}}}^{{s_{*} }} = & {\text{d}}t_{{,{\text{IF}}}}^{{s_{*} }} - {\text{d}}t_{{1,{\text{IF}}}}^{G} - d_{{1,{\text{IF}}}}^{G*} \quad * \ne G \\ \end{aligned}$$
(12)

In the presented combination, our interest is the satellite clocks \(\tilde{{\text{d}}}{t}_{,{\text{IF}}}^{{s}_{*}}\). All the alignments are performed before the AC clocks enter the combination process. As pointed out before, it is challenging to find a common reference station within all solutions, and the ISBs are not provided by all the MGEX ACs; therefore, equation (12) cannot always be applied directly. To overcome this issue, a common offset can directly be estimated for the satellite clocks of each constellation (not the reference system) with respect to a reference AC after the AC alignment of the previous section has been applied. To exemplify the effect of different reference ISBs, Fig. 3 shows the difference between the GFZ and COD satellite clocks for all constellations. In the top plot of Fig. 3, a clear but consistent offset between the constellations is noted, reaching more than \(20\hspace{0.17em}\mathrm{ ns}\) for BeiDou. For comparison, in the bottom plot of Fig. 3, the clock offsets are presented after the ISB alignment, where all constellations show similar clock differences over the constellations, with values between \(\pm 2.5\hspace{0.17em}\mathrm{ ns}\).

Fig. 3
figure 3

MGEX clock difference between COD and GFZ before (top) and after (bottom) the ISB alignment for all satellites over all constellations

Full size image

Radial orbit correction

In addition to a common clock and ISB reference, the clock offsets from the individual ACs should be corrected for the radial orbit differences in the AC solutions since the radial orbit component and the satellite clock offset are highly correlated (Martinez 2014). The correction \({\delta }_{r}^{s}\) applied to the satellite clock offsets consists of the radial orbit difference between the combined and the AC products:

$$\delta_{r}^{s} = \frac{{\left( {\hat{{\underline {x} }}_{{\text{c}}}^{s} - \underline{{\bar{x}}}_{r}^{s} } \right)^{\text{T}} \underline{{\bar{x}}}_{r}^{s} }}{{\left| {\left| {\underline{{\bar{x}}}_{r}^{s} } \right|} \right|c}}$$
(13)

where \({\hat{\underline{x}}}_{c}^{s}\) describes the combined orbit solution, \({\underline{\bar{x}}}_{r}^{s}\) the AC solution, and \(c\) the speed of light. The combined orbit solution used in this computation comes from a previous orbit combination as presented in Mansur et al. (2022).

Figure 4, for example, shows the GPS satellite G17 (SVN53, block IIRM), where the difference between clock offsets from GFZ and COD is computed, as well as the potential radial correction compensating the difference in radius between the GFZ and COD orbit products. We note that after applying the correction (blue line), the clocks are closer to a linear behavior (green line) than in the pure differences (orange line). The maximum corrections are in the order of \(0.1\hspace{0.17em}{\text{ns}},\) where one nanosecond is equivalent to a range difference of approximately 30 cm.

Fig. 4
figure 4

Satellite clock difference between GFZ and COD (satellite G17), radial correction

Full size image

Due to the already mentioned correlation, clock and orbit products should have equivalent epoch intervals. However, typically, orbits are given in 5 or 15-min intervals; while, clock products are given in 30-s or 5-min intervals. Consequently, for example, combining 5-min orbits and 30-s clock products requires an orbit interpolation and extrapolation for data beyond 23h 55min. We performed a test using MGEX orbits from DOY 235-341, 2022 to determine the difference between interpolated and original orbits. The original orbits were computed using GFZ’s GNSS processing software EPOS.P8 with a 30-s sampling interval and including the 24-h epoch. This interval was then decimated to create “artificial” intervals of 15 and 5 min, upon which an interpolation approach using a 10th-order Lagrange polynomial was applied to return it to a 30-s interval. A simple linear extrapolation was used at the end of the day boundary. Finally, the resulting orbits were then compared to the original ones. The 5-min orbit agrees at the \(1\hspace{0.17em}{\text{mm}}\) level. Bigger differences are found for the 15-min interpolated orbits, especially at the end of the day, reaching \(10\hspace{0.17em}{\text{mm}}\) for some satellites. In addition, the 3D-RMS of the geometric differences were computed, where we observe that particularly for satellites with eccentric orbits (E14 and E18) the 5-min interpolation has a better performance with RMS values below \(1\hspace{0.17em}{\text{mm}}\), while the 15-min can reach 60 \({\text{mm}}\) for certain days (Figs. 5, 6). Note that the high RMS values might be related to specific problematic epochs, like the day boundaries.

Fig. 5
figure 5

Daily RMS of orbit coordinates between 30 s and 15-min sampled orbits for the two Galileo eccentric satellites (E14 left, E18 right). The period used was in the year 2022

Full size image
Fig. 6
figure 6

Daily RMS of orbit coordinates between 30 s and 5-min sampled orbits for the two Galileo eccentric satellites (E14 left, E18 right). The period used was in the year 2022

Full size image

Clock least squares variance component estimation

With all the alignments done we can now define our combination system model. Let \(\underline{\text{d}} t_{r} \left( t \right) \in \mathbb{R}^{S}\) contain the aligned clock estimates for AC \(r \in \{1,\dots , R\}\) at time \(t\), so that:

$$\underline{\text{d}} t_{r} = \left[ {\underline{\text{d}} t_{r} \left( {t_{1} } \right), \ldots ,\underline{\text{d}} t_{r} \left( {t_{P} } \right)} \right]^{{\text{T}}} \in \mathbb{R}^{{{\text{SP}}}}$$
(14)

where \({t}_{p}\), \(p \in \{1,\dots , P\}\), the time at epoch \(p\). The clock system model is then given as follows:

$$\left[ {\begin{array}{*{20}c} {\underline{\text{d}} t_{1} } \\ \vdots \\ {\underline{\text{d}} t_{R} } \\ \end{array} } \right] = \left[ {\underline {1}_{R} \otimes {\varvec{I}}_{{\text{SP}}} } \right] \cdot \underline{\text{d}} t_{c} + \left[ {\begin{array}{*{20}c} {\underline {\vartheta }_{1} } \\ \vdots \\ {\underline {\vartheta }_{R} } \\ \end{array} } \right]$$
(15)

where \(\underline{\vartheta } _{r} \in \mathbb{R}^{{{\text{SP}}}}\) is an additive noise term. \({\underline{1}}_{R}\) is an \(R\)-vector of ones, \({{\varvec{I}}}_{{\text{SP}}}\) an identity matrix of size \({\text{SP}}\) and \(\underline{\text{d}} t_{c}\) is the unknown combined clock vector. We assume for the stochastic model that \(\underline{\vartheta } _{r} \sim {\mathcal{N}}\left( {\underline{0} ,\sigma _{r}^{2} \varvec{I}_{{{\text{SP}}}} } \right)\) and with the uncorrelated solutions between ACs we get the full noise vector \(\underline{\vartheta }\) containing all \(R\) AC contributions from (15),

$$\underline {\vartheta } \simeq {\mathcal{N}}\left( {\underline {0} ,{\varvec{Q}}} \right){\text{ with }}{\varvec{Q}} = \mathop \sum \limits_{r = 1}^{R} \sigma_{r}^{2} {\varvec{Q}}_{r}$$
(16)

with \({{\varvec{Q}}}_{r}={\text{diag}}\left({\underline{c}}_{r}\right)\otimes {{\varvec{I}}}_{{\text{SP}}}\), where \({\underline{c}}_{r}\) is the rth canonical unit vector in \({\mathbb{R}}^{R}\). The models presented in (15) and (16) can be solved for the unknown combined clock \(\underline{\text{d}}t_{c}\) and the variance components \({\sigma }_{r}^{2}\) (see Mansur et al. 2022). The combined least squares clock solution is derived as follows:

$$\underline{\hat{\text{d}}} t_{c} = \frac{1}{{\sum _{{r = 1}}^{R} \frac{1}{{\sigma _{r}^{2} }}}} \cdot \sum _{{r = 1}}^{R} \frac{1}{{\sigma _{r}^{2} }} \cdot \underline{\text{d}}t_{r}$$
(17)

and depends on the unknown AC-specific variance factors\({\sigma }_{r}^{2}\). Following the LSVCE theory in Teunissen and Amiri-Simkooei (2008), these are estimated only from the random vector of misclosures, which is uncorrelated with the clock estimates \(\underline{\hat{\text{d}}}t_{c}\) in (17) and of the same dimension as the redundancy in (15). The uncorrelated property implies that the joint least squares estimation of the combined clock \(\underline{\text{d}}t_{c}\) and the variance components \({\sigma }_{r}^{2}\) can be done consecutively by solving two separate least squares problems, starting with the variance components. For the given model (15), one valid example for the vector of misclosures contains the differences \(\underline{\text{d}}t_{r}\)−\(\underline{\text{d}}t_{1}\) of the AC clock solutions, meaning that the \(R\) variance components \({\sigma }_{r}\) are estimated in a least squares adjustment from the \(R-1\) vectors of clock differences. The solution is found iteratively as given in Eqs. (15), (22) and (23) in Mansur et al. (2022).

The previous formulation gives us AC-based weights considering all constellations at once. An extension for AC-constellation-specific weights can be derived by introducing additional variance components \({}_{g}\,{\upsigma }_{r}^{2}\),\(g\in \{ 1,\dots , G\}\), for each AC and constellation, with \(G\) as the number of constellations. Our input clocks are then defined as \({}_{g}\underline{\text{d}}t\in {\mathbb{R}}^{{S}_{g}PR}\), which contains only the clock offset of a specific constellation so that the adapted system model is given as follows:

$$\left[ {\begin{array}{*{20}c} {{}_{1}^{{}} \underline{\text{d}}t } \\ \vdots \\ {{}_{G}^{{}} \underline{\text{d}}t } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\underline {1}_{R} \otimes {\varvec{I}}_{{S_{1} P}} } & {} & {} \\ {} & \ddots & {} \\ {} & {} & {\underline {1}_{R} \otimes {\varvec{I}}_{{S_{G} P}} } \\ \end{array} } \right] \cdot \left[ {\begin{array}{*{20}c} {{}_{1}^{{}} \underline{\text{d}}t_{c} } \\ \vdots \\ {{}_{G}^{{}} \underline{\text{d}}t_{c} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {{}_{1}^{{}} \underline {\vartheta } } \\ \vdots \\ {{}_{G}^{{}} \underline {\vartheta } } \\ \end{array} } \right]$$
(18)

The covariance matrix \({\varvec{Q}}\) is then blockdiagonal where each element only depends on \({}_{g}{\upsigma }_{r}^{2}\) as we assume no correlations between different constellations, so that the combination and LSVCE can be done separately for each constellation.

Clock combination workflow

The clock combination can be divided into seven steps, summarized in the following:

  1. 1.

    The radial correction is applied to the ACs clock files (clk). The radial correction is based on the previous combined orbit, using also the LSVCE.

  2. 2.

    Initial check of the AC products, verifying if the ACs provides clock offsets for the entire day, as well as standardizing the sample rate (in our combination 5-min sampling rates are used). In addition, a simple detection per satellite is computed, where individual values are analyzed and, if needed, removed when detected as an outlier.

  3. 3.

    With the standardized data, the reference AC is selected based on a median trend computation; and then, AC-wise offsets and trends are computed. These offsets and trends are then applied in order to bring the clock offsets to a common reference.

  4. 4.

    A second outlier detection is performed based on a satellite-wise RMS computed for each AC. This is done to find problematic satellites or, in some worst cases, problematic ACs. RMS outlier values are detected through the modified Z-score method (Iglewicz et al. 1993).

  5. 5.

    Following the outlier step, the ISB alignment is applied as previously mentioned, computing an offset value between the GPS constellation and the others.

  6. 6.

    A final outlier control is done per satellite comparing the clock values among the ACs epoch-wise. This detection is again performed using the modified Z-score method. In the end, this outlier check allows to define a set of core satellites, which is the subset of satellites that will be used for the VCE. It contains only satellites that are available for all ACs (per constellation in case of constellation-specific weights). This outlier detection plays an important role in the processing, assuring that the VCE can be robustly computed.

  7. 7.

    Finally, the VCE is computed, and the combined clocks are computed as a weighted mean of the ACs products, where the weights are proportional to the inverse of the variances.

The main goal of the developed combination is to enable multi-GNSS PPP applications, which require a consistent set of orbits and clocks but not necessarily an alignment to a time global reference. We note, however, that a final alignment could be done, for instance, with respect to the broadcast clocks which are transmitted as part of the GNSS navigation messages, which is done in the official IGS combination, or with respect to a dedicated timing station, possibly even with calibrated ISBs. When performing this alignment, one must be careful not to violate the assumption of constant ISBs, for instance, by applying different offsets and trends for different constellations. To ensure complete clock products, satellites identified as outliers are still contained in the final result of the combination.

Combination of IGS MGEX products

The clock combination is performed using the individual ACs solutions submitted to MGEX. Table 1 summarizes the constellations and sampling intervals provided by each AC. The products are combined between GPS week 1982 (January 2018) and GPS week 2116 (July 2020). Since only three ACs provide all five constellations, we chose to use the AC-constellation weighting strategy, see (18), which does not require complete AC solutions covering all five constellations.

Table 1 MGEX clk ACs’ contributions and constellations
Full size table

MGEX differences between individual AC clock solutions and combined clocks

The comparison between the combined clock solution and the individual ACs is shown in Fig. 7. We note that for some constellations, JAX and SHA are not in the plot since their values are too dissonant from the main ACs levels. Galileo shows the best agreement between the ACs in the clock combination with ACs’ RMS below \(80\hspace{0.17em} {\text{psec}}\), and an average RMS of \(45 \hspace{0.17em}{\text{psec}}\). SHA was excluded from the average computation since this AC presents high RMS values around \(122\hspace{0.17em} {\text{psec}}\). The second best is GPS, where most of the ACs’ RMS is around \(50 \hspace{0.17em}{\text{psec}}\) while GFZ shows an RMS of around \(100\hspace{0.17em} {\text{psec}}\). We note that SHA’s RMS values are around \(269\hspace{0.17em} {\text{psec}}\) and JAX solution shows an RMS of \(860 \hspace{0.17em}{\text{psec}}\) before September 2019. After configuration changes JAX’s RMS values are below \(50 \hspace{0.17em}{\text{psec}}\). For BeiDou, the COD and WHU agree with mean RMS of around \(70 \hspace{0.17em}{\text{psec}}\) and \(78\hspace{0.17em} {\text{psec}}\), respectively, whereas from the middle of July 2019, we note an increase in the RMS for COD. GFZ has a mean RMS of \(214\hspace{0.17em}{\text{psec}}\) for BDS while reaching an RMS of \(120\hspace{0.17em} {\text{psec}}\) in July 2020. Again, SHA is not shown in the plot due to high RMS values of more than \(1000\hspace{0.17em} {\text{psec}}\). For GLONASS, the lowest RMS is found for GRG with around \(70 \hspace{0.17em}{\text{psec}}\); while, COD has \(72 \hspace{0.17em}{\text{psec}}\), WHU \(88\hspace{0.17em} {\text{psec}}\), and GFZ regularly above \(200 \hspace{0.17em}{\text{psec}}\). Once more, SHA values exceed \(1000 \hspace{0.17em}{\text{psec}}\), and for JAX, the values are at the same level as SHA until September 2019. After this date, the JAX values are comparable with the other ACs, reaching an RMS of \(98 \hspace{0.17em}{\text{psec}}\). The most problematic constellation is QZSS, where we noted already modeling issues for the satellite orbits (Mansur et al. 2022). The best RMS values are of around \(150\hspace{0.17em} {\text{psec}}\) for WHU and COD. GFZ and JAX present similar behavior with extremely high RMS values, and only in September 2019 did they reach a better level, getting to the same level as the other ACs.

Fig. 7
figure 7

Constellation specific RMS differences between the individual AC clocks and the combined solution for AC plus constellation weighting

Full size image

The results in our RMS comparison are in line with the analysis of individual ACs solutions presented by Li et al. (2020), especially regarding the behavior of SHA. Unfortunately, there is a lack of documentation about the MGEX processing from most of the ACs, making it difficult to validate RMS variations concerning processing changes.

One challenge of the clock combination relies on the computation of the core satellites which are used in the VCE. This means that a strict definition is required; otherwise, the process either does not converge or yields negative VCE values. This was strongly noted for the QZSS constellation, where the modeling from the ACs is still in constant improvement, but also for GLONASS, where the individual clocks need to be properly aligned with individual (satellite-wise) ISBs, related to the FDMA (frequency-division multiple access) technique. In general, COD shows the best RMS performance, followed by GRG and WHU. These three ACs are very similar in their behavior, which is also visible in the derived weights.

MGEX Clock weights

Figure 8 shows the weights computed in the LSVCE using the constellation-specific strategy. Overall, the weights presented are compatible with the RMS shown in the previous Sections. For GPS, SHA is downweighted on average to \(1\mathrm{\%}\) except in the first two months of the processing period, where it contributed with \(15\mathrm{\%}\). Also, with low weight, GFZ contributed, on average, \(6\%\) for the entire processing time. As mentioned in the RMS analysis, JAX became more compatible with the other ACs in September 2019, leading to an average weight of \(22\mathrm{\%}\). The main contributions for GPS are from WHU (\(20\mathrm{\%}\)), GRG (\(29\mathrm{\%}\)), and COD (\(37\mathrm{\%}\)). For the Galileo constellation, the weights are more evenly distributed, between \(22\%\) and \(27\%\) for the ACs, except for SHA, which contributes by only \(5\mathrm{\%}\). For GLONASS, GRG dominates with \(38\mathrm{\%},\) while the weight distribution improves after September 2019. SHA and GFZ contribute with less than \(1\mathrm{\%}\). Unlike the other constellations, BeiDou has two dominant ACs, WHU, with an average of 56%, and COD, with 48%, while the other two ACs have a small impact. GFZ gets higher weight values by the end of the processing test time up to (\(10\mathrm{\%}\)). Weights for QZSS are again mainly divided between WHU and COD. Equally, as for the previous constellation, JAX started to be more compatible with QZSS in September 2019. Subsequently, there is a better agreement between the four ACs, but still, GFZ contributes only \(10\mathrm{\%}\), while COD dominates with \(44\mathrm{\%}\).

Fig. 8
figure 8

Weights for MGEX AC based on constellation-specific weighting

Full size image

Overall, the weights are evenly distributed for constellations which are consistently modeled by the ACs, especially GPS and Galileo. For COD we note very stable weights for all MGEX constellations, indicating that this AC provides very stable solutions.

Combination of IGS legacy products

As a validation of our approach, the ACs clock products from the IGS legacy chain are used as input for the combination. The result is then compared with products from the official IGS combination after aligning them to the reference AC of the combination. Figure 9 shows the RMS comparing the combination to the individual solutions and the official IGS final and rapid products. On average, the IGS final combination has the smallest RMS with \(32 \hspace{0.17em}{\text{psec}}\text{,}\) while in the IGR solution, COD, and JPL have values around \(36\hspace{0.17em} {\text{psec}}\). The most incompatible ACs are GFZ with an average RMS of \(96 \hspace{0.17em}{\text{psec}}\) and MIT with \(65\hspace{0.17em} {\text{psec}}\). The GFZ solution uses a reference station with a nonlinear behavior, which yields an issue in the clock alignment of our combination process. As explained earlier, this issue can be removed by switching the reference clock of the GFZ solution before the alignment steps, and the average RMS is then reduced to \(53 \hspace{0.17em}{\text{psec}}\), which is similar to the other ACs. The major weights computed without this additional adjustment are on average, around \(27\mathrm{\%}\) for COD, followed by JPL with \(21\mathrm{\%}\) and \(14\mathrm{\%}\) for GRG. Overall, the resulting AC-specific RMS values are reasonable and compatible with the official IGS and IGR combinations, leading to the conclusion that our combination approach provides a fair clock combination.

Fig. 9
figure 9

RMS clock differences compared to the combined solution for the GPS only legacy processing. The IGS rapid (igr) and final (igs) solutions are included for validation

Full size image

Combination of IGS repro3 products

As a contribution to the new ITRF2020 realization, the IGS ACs reprocessed over more than 25 years of GNSS observations (Rebischung 2021). The reprocessing was performed within a common agreement on the application of the processing models as well as the inclusion of Galileo for some ACs. We combine the repro3 clocks from GPS week 2035 (January 2019) until 2138 (December 2020). The combined repro3 orbit used to perform the radial corrections and PPP tests is computed also using the LSVCE methodology; more information is in Sakic et al. (2022).

Repro3 differences between individual AC clock solutions and combined clocks

In general, for all constellations, the RMS values agree well among the ACs, especially for GPS, as shown in Fig. 10. For this constellation, the ACs have an average RMS of \(41 \hspace{0.17em}{\text{psec}}\), where COD has the lowest RMS of \(27 \hspace{0.17em}{\text{psec}}\) and GFZ the highest of \(62\hspace{0.17em} {\text{psec}}\). For Galileo, the average RMS is \(53 \hspace{0.17em}{\text{psec}}\), and again COD has the lowest RMS with \(35\hspace{0.17em} {\text{psec}}\). The constellation with the highest average RMS is GLONASS with \(68 \hspace{0.17em}{\text{psec}}\). Again, COD shows the lowest with an average value of \(50 \hspace{0.17em}{\text{psec}}\). The compatibility among the ACs in the repro3, when compared with the MGEX solutions, is better for the three constellations, which can be related to the agreed satellite modeling.

Fig. 10
figure 10

Repro3 constellation specific RMS differences between the individual AC clocks and the combined solution

Full size image

The IGS-ACC also provided combined clocks in September 2022, based on the combination workflow developed at Wuhan University (Geng et al. 2021). The comparison between both combinations shows an RMS of \(93 \hspace{0.17em}{\text{psec}}\) for GPS, \(97 \hspace{0.17em}{\text{psec}}\) for Galileo and GLONASS with \(101 \hspace{0.17em}{\text{psec}}\). The agreement between both solutions is around \(0.1\hspace{0.17em}{\text{ns}}\), noting that this agreement is the pure difference without any alignment between the solutions.

The weights used in the repro3 combination are shown in Fig. 11. In general, the weights are evenly distributed for all constellations, with COD being the AC with the highest percentage values, which is reflected in the low RMS values previously presented.

Fig. 11
figure 11

Repro3 clock constellation-specific weights

Full size image

Repro3 PPP processing

To validate orbits and clocks, a daily PPP processing for 50 stations was performed for the year 2020, using our combined orbit and clock solution (CMG) together with two other ACs (COD and GFZ) and the repro3 combination products provided by the IGS-ACC. Figure 12 shows the daily RMS residuals of the ionosphere-free phase linear combination for all the stations for each constellation. All the RMS values are similar, differing by less than one millimeter between the solutions and following the same pattern where GLONASS has the high residual with about \(11 \hspace{0.17em}{\text{mm}}\) and the lowest value for GPS with \(9.5\hspace{0.17em} {\text{mm}}\). Figure 13 shows the repeatability (RMS of station coordinate) for the three coordinate components based on a GPS only PPP. The processing was limited to GPS since this constellation is common in all the stations. The first three plots show the results for three different stations (AREQ, CAS1, SCOR) over more than 300 processed days, while the last plot shows the RMS averaged over all the stations. The results of the ACs solutions and both combinations are at the same level, with a few millimeters difference. For the station AREQ, we note that both combined products have a better performance than the other two ACs, mainly differing around \(1\hspace{0.17em}\mathrm{ mm}\) for each component, except for the COD, where the height component RMS is \(13 \hspace{0.17em}{\text{mm}}\), while for the combined solution is around \(6\hspace{0.17em} {\text{mm}}\). For CAS1, the CMG solution has the lower RMS, especially for the East component, with \(6\hspace{0.17em} \text{mm,}\) while the other solutions are above \(15\hspace{0.17em} {\text{mm}}\).

Fig. 12
figure 12

Daily RMS of the Ionosphere-free phase residual per AC for all stations processed in 2020 (IGS = final IGS combination; CMG = our approach)

Full size image
Fig. 13
figure 13

Station stability considering the coordinates repeatability derived by GPS only PPP processing during the year 2020. Note that the y-axis scale is different in the bottom plot (IGS = final IGS combination; CMG = our approach)

Full size image

The results among the solutions are similar for the station SCOR, where larger differences occur in the East component with RMS of \(8.4\hspace{0.17em} {\text{mm}}\) (IGS), and \(3.8 \hspace{0.17em}{\text{mm}}\) (COD). The RMS over all stations, shows similar values for all solutions, with values of around \(12\hspace{0.17em} {\text{mm}}\) for the North, \(20\hspace{0.17em}{\text{mm}}\) for East, and \(22 \hspace{0.17em}{\text{mm}}\) for the Up components.

Conclusions

In this study, we present a mathematical framework to combine clock solutions from different ACs in the multi-GNSS context, using the LSVCE to determine the individual weights to be applied in the actual combination. As the AC solutions make use of different clock and bias references, a series of alignments are required. The results show that the radial correction can reach, on average, \(1\hspace{0.17em}{\text{cm}}\) and, therefore is relevant to be applied, thus eliminating the orbital errors absorbed by the clocks. Further, a 5-min orbit is sufficient for the orbit interpolation, which is required to keep the consistency between clock and orbits. In addition, providing the 24-h epoch in the orbit solution would also be beneficial since the extrapolation at the end-of-day boundaries might not express the real behavior of the orbits, especially for satellites with eccentric orbits. The alignment used in the combination process selects one AC as a reference and aligns the others via drift and offset determination to this reference. This assumes that the reference clocks have linear behavior, which is true in most cases. Special cases where some AC does not provide a clock with linear behavior need to be handled by changing the reference clock to a clock with linear behavior within the AC solution. The next alignment applied refers to the ISBs present in multi-GNSS solutions. The ACs usually compute the ISBs using GPS as a reference, and a different ISB definition within each AC solution leads to a common offset of the satellite clocks for each constellation between different ACs. Not all ACs provide information about their ISBs, so we present a solution to compensate for these offsets through an alignment to a reference AC. Ideally, each AC would provide relevant information about the definition of their clock and bias parameters as well as their bias values so that this information could be used in the combination, as shown in Geng et al. (2021).

The combination process was first tested using MGEX ACs solutions for a period of 2 years, where the RMS for Galileo and GPS have similar values, both under \(80\hspace{0.17em}{\text{psec}}\) for most of the ACs. GLONASS and BeiDou show RMS, on average, under \(150\hspace{0.17em}{\text{psec}}\), and QZSS shows the highest RMS values. Further, the validation of the combination using the legacy products has shown similar results when compared with the official IGS final and rapid solution with RMS values of around 32 psec on average. A combination of the repro3 products was performed as well, showing improved consistency compared to the MGEX. Smaller RMS values for all constellations indicate good agreement between the ACs. The RMS is mostly below \(65 \hspace{0.17em}{\text{psec}}\), and the RMS with the IGS ACC clocks is below \(100 \hspace{0.17em}{\text{psec}}\). To validate orbits and clocks, a PPP processing was performed. The ionospheric-free phase residuals for the combined solution are on the same level as for the ACs and are compatible with results obtained from the IGS ACC combined solution. A comparison of coordinate repeatability showed that the combination is stable over time, and the coordinate RMS is at the same level for the four compared solutions.

Combining clock offsets requires a proper understanding of the alignments since in each AC, the clock offsets are computed with different reference clocks. Adding more constellations also brings a new challenge, where the biases need to be fully comprehended and properly treated. It would be beneficial for the combination if there was an agreement among the ACs to provide the daily bias products on a regular basis. Furthermore, based on the experiments, we can conclude that the combination process developed can achieve reliable results that are compatible with the existing ACs solutions. In addition, the combined solution can be used to assess the quality of the individual ACs and provide the necessary feedback. Finally, the combined solution provides a five-system solution since one of the goals of having a combination is to make the most complete product, aiming for the highest accuracy available to the users.