# Utilizing least squares variance component estimation to combine multi-GNSS clock offsets

Gustavo Mansur, Andreas Brack, Pierre Sakic, Benjamin Männel, Harald Schuh

► To cite this version:

Gustavo Mansur, Andreas Brack, Pierre Sakic, Benjamin Männel, Harald Schuh. Utilizing least squares variance component estimation to combine multi-GNSS clock offsets. GPS Solutions, 2024, 28 (2), pp.70. 10.1007/s10291-023-01604-4. hal-04873366

HAL Id: hal-04873366

https://hal.science/hal-04873366v1

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Utilizing least squares variance component estimation to combine multi-GNSS clock offsets

Gustavo Mansur<fnref n="1" />,<fnref n="2" /> · Andreas Brack<fnref n="1" /> · Pierre Sakic<fnref n="3" /> · Benjamin Männel<fnref n="1" /> · Harald Schuh<fnref n="1" />,<fnref n="2" />

Received: 9 October 2023 / Accepted: 20 December 2023 / Published online: 9 February 2024
© The Author(s) 2024

## Abstract

The International GNSS Service (IGS) provides combined satellite and station clock products, which are generated from the individual clock solutions produced by the analysis centers (ACs). Combinations for GPS and GLONASS are currently available, but there is still a lack of combined products for the new constellations such as Galileo, BeiDou, and QZSS. This study presents a combination framework based on least squares variance component estimation using the ACs' aligned clock solutions. We present the various alignments required to harmonize the solutions from the ACs, namely the radial correction derived from the differences of the associated orbits, the alignment of the AC clocks to compensate for different reference clocks within each AC solution, and the inter-system bias (ISB) alignment to correct for different AC ISB definitions when multiple constellations are used. The combination scheme is tested with IGS MGEX and repro3 products. The RMS computed between the combined product and the aligned ACs' solutions differ for each constellation, where the lowest values are obtained for Galileo and GPS with an average below 45 psec (13 mm) and reaching more than 150 psec (45 mm) for QZSS. The same behavior is repeated when the process is performed with the repro3 products. A clock and orbit combination validation is done using precise point positioning (PPP) that shows ionosphere-free phase residuals below 10 mm for all constellations, comparable with the AC solutions that are in the same level.

**Keywords** Multi-GNSS · Clock combination · Variance component estimation · Clock bias · Clock offset

## 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

<figure><img src="image_2.png" /></figure>

<footnote>Gustavo Mansur
mansur@gfz-potsdam.de</footnote>

<footnote n="1">¹ GFZ German Research Centre for Geosciences,
Telegrafenberg, 14473 Potsdam, Germany</footnote>

<footnote n="2">² Technische Universität Berlin, Str. des 17. Juni 135,
10623 Berlin, Germany</footnote>

<footnote n="3">³ CNRS, Institut de Physique du Globe de Paris, Université Paris Cité, 1 Rue Jussieu, 75005 Paris, France</footnote>
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} + dt_{k} - dt^{s} + m_{k}^{s} \tau_{k} + \mu_{f} t_{k}^{s} + d_{k,f} - d_{f}^{s} + \epsilon_{k,f}^{s} \\
\varphi_{k,f}^{s} &= \rho_{k}^{s} + dt_{k} - dt^{s} + m_{k}^{s} \tau_{k} - \mu_{f} t_{k}^{s} + \delta_{k,f} - \delta_{f}^{s} + a_{k,f}^{s} + \epsilon_{k,f}^{s}
\end{aligned} 
\quad (1) $$

with $\rho_k^s$ the receiver-satellite ranges, $dt_k$ the receiver clock offset, $dt^s$ the satellite clock offset, $\tau_k$ the tropospheric delay, $m_k^s$ the mapping function, $t_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,IF} &= \rho_k + dt_k - dt^s + m_k^s \tau_k + d_{k,IF} - d_{IF}^s + \epsilon_{k,IF}^s \\
\varphi_{k,IF} &= \rho_k + dt_k - dt^s + m_k^s \tau_k + \delta_{k,IF} - \delta_{IF}^s + a_{k,IF}^s + \epsilon_{k,IF}^s
\end{aligned} 
\quad (2) $$

with $(\cdot)_{IF} = \frac{1}{\mu_2 - \mu_1} (\mu_2(\cdot)_1 - \mu_1(\cdot)_2)$. 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}
dt_{k,IF} &= dt_k + d_{k,IF} \\
dt_{IF}^s &= dt^s + d_{IF}^s
\end{aligned} 
\quad (3) $$

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

$$ 
\begin{aligned}
p_{k,IF}^s &= \rho_k^s + dt_{k,IF} - dt_{IF}^s + m_k^s \tau_k + \epsilon_{k,IF}^s \\
\varphi_{k,IF}^s &= \rho_k^s + dt_{k,IF} - dt_{IF}^s + m_k^s \tau_k + (\delta_{k,IF} - d_{k,IF} - \delta_{IF}^s + d_{IF}^s + a_{k,IF}^s) + \epsilon_{k,IF}^s
\end{aligned} 
\quad (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, $dt_{1,IF}$, and the estimable clock parameters are:

$$ 
\begin{aligned}
\tilde{d}_{t_{k,IF}} &= dt_{k,IF} - dt_{1,IF} \quad k \neq 1 \\
\tilde{d}_{t_{IF}^s} &= dt_{IF}^s - dt_{1,IF}
\end{aligned}
\qquad (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. Choose a common reference station provided by all the ACs.

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{d}_{k,IF} &= \tilde{d}t_{k,IF} - \tilde{d}t_{l,IF} \\
&= dt_{k,IF} - dt_{1,IF} - dt_{l,IF} + dt_{1,IF} \\
&= dt_{k,IF} - dt_{l,IF} \\
\tilde{d}_{IF}^s &= \tilde{d}t_{IF}^s - \tilde{d}t_{l,IF} \\
&= dt_{IF}^s - dt_{1,IF} - dt_{l,IF} + dt_{1,IF} \\
&= dt_{IF}^s - dt_{l,IF}
\end{aligned}
\qquad (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

<figure><img src="image_3.png" /><figcaption>Fig. 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</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Comparison between GFZ and COD satellite clocks aligned to MIT</figcaption></figure>

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.

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).

## 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,IF}^{s*} &= \rho_{k}^{s*} + dr_{k,IF}^{s*} - dr_{IF}^{s*} + m_{r}^{s*} \tau_{r} + \epsilon_{k,IF}^{s*} \\
q_{k,IF}^{s*} &= \rho_{k}^{s*} + dr_{k,IF}^{s*} - dr_{IF}^{s*} + m_{k}^{s*} \tau_{k} + (\delta_{k,IF}^{s*} - d_{k,IF}^{s*} - \delta_{IF}^{s*} + d_{IF}^{s*} + a_{k,IF}^{s*}) + \epsilon_{k,IF}^{s*}
\end{aligned}
\quad (7) $$

with this formulation the index $s*$ as well as the receiver hardware biases $d_{k,f}^*$ and $\delta_{k,f}^*$ are system specific and the estimable receiver clocks $dr_{k,IF}^*$ are given as follows:

$$ dt_{k,IF}^{s*} = dt_{k} + d_{k,IF}^{s*} \quad (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 $dr_{k,IF}^G$, and time-constant ISB parameters $d_{k,IF}^{G*}$ for all other systems where:

$$ 
\begin{aligned}
dt_{k,IF}^{G} &= dt_{k} + d_{k,IF}^{G} \\
dt_{k,IF}^{s*} &= dt_{k,IF}^{G} + d_{k,IF}^{s*} - d_{k,IF}^{G} \quad * \neq G \\
& \underbrace{\quad \text{ISB} d_{k,IF}^{G*}}_{\quad}
\end{aligned}
\quad (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,IF}^{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,IF}^{G*} = d_{k,IF}^{G*} - d_{1,IF}^{G*} \quad k \neq 1 \quad (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,IF}^{G*}$ in addition to the clock of the reference receiver:

$$ \tilde{dt}_{IF}^{s*} = dt_{IF}^{s*} - dr_{1,IF}^{G*} - d_{1,IF}^{G*} \quad * \neq G \quad (11) $$

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

$$ 
\begin{aligned}
\tilde{dt}_{k,IF}^{G} &= dt_{k,IF}^{G} - dr_{1,IF}^{G} \quad k \neq 1 \\
\tilde{dt}_{IF}^{s*} &= dt_{IF}^{s*} - dr_{1,IF}^{G} \\
\tilde{d}_{k,IF}^{G*} &= d_{k,IF}^{G*} - d_{1,IF}^{G*} \quad k \neq 1, * \neq G \\
\tilde{dt}_{IF}^{s*} &= dt_{IF}^{s*} - dr_{1,IF}^{G} - d_{1,IF}^{G*} \quad * \neq G
\end{aligned}
\quad (12) $$

In the presented combination, our interest is the satellite clocks $\tilde{dt}_{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 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 ±2.5 ns.

## 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{(\hat{x}_c^s - \bar{x}_r^s)^T \bar{x}_r^s}{\|\bar{x}_r^s\| \|c\|} \quad (13) $$

where $\hat{x}_c^s$ describes the combined orbit solution, $\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
<figure><img src="image_5.png" /><figcaption>Fig. 3 MGEX clock difference between COD and GFZ before (top) and after (bottom) the ISB alignment for all satellites over all constellations</figcaption></figure>

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 ns, where one nanosecond is equivalent to a range difference of approximately 30 cm.

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 mm level. Bigger differences are found for the 15-min interpolated orbits, especially at the end of the day, reaching 10 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 mm, while the 15-min can reach 60 mm for certain days (Figs. 5, 6). Note that the high RMS values might be related to specific problematic epochs, like the day boundaries.

## Clock least squares variance component estimation

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

$$ \underline{dt}_r = [\underline{dt}_r(t_1), ..., \underline{dt}_r(t_P)]^T \in \mathbb{R}^{SP} \quad (14) $$

where $t_p, p \in \{1, ..., P\}$, the time at epoch $p$. The clock system model is then given as follows:

$$ \begin{bmatrix} \underline{dt}_1 \\ \vdots \\ \underline{dt}_R \end{bmatrix} = [\underline{1}_R \otimes I_{SP}] \cdot \underline{dt}_c + \begin{bmatrix} \underline{\theta}_1 \\ \vdots \\ \underline{\theta}_R \end{bmatrix} \quad (15) $$

where $\underline{\theta}_r \in \mathbb{R}^{SP}$ is an additive noise term. $\underline{1}_R$ is an $R$-vector of ones, $I_{SP}$ an identity matrix of size SP and $\underline{dt}_c$ is the unknown combined clock vector. We assume for the stochastic model that $\underline{\theta}_r \sim \mathcal{N}(\underline{0}, \sigma_r^2 I_{SP})$ and with the uncorrelated solutions between ACs we get the full noise vector $\underline{\theta}$ containing all $R$ AC contributions from (15),

$$ \underline{\theta} \approx \mathcal{N}(\underline{0}, \underline{Q}) \quad \text{with} \quad \underline{Q} = \sum_{r=1}^{R} \sigma_r^2 \underline{Q}_r \quad (16) $$

with $\underline{Q}_r = \text{diag}(c_r) \otimes I_{SP}$, where $c_r$ is the $r$th canonical unit vector in $\mathbb{R}^R$. The models presented in (15) and (16) can be solved for the unknown combined clock $\underline{dt}_c$ and the variance components $\sigma_r^2$ (see Mansur et al. 2022). The combined least squares clock solution is derived as follows:

$$ \underline{dt}_c = \frac{1}{\sum_{r=1}^{R} \frac{1}{\sigma_r^2}} \cdot \sum_{r=1}^{R} \frac{1}{\sigma_r^2} \cdot \underline{dt}_r \quad (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{dt}_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{dt}_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
<figure><img src="image_6.png" /><figcaption>**Fig. 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</figcaption></figure>

<figure><img src="image_7.png" /></figure>

valid example for the vector of misclosures contains the differences $\underline{dt}_r - \underline{dt}_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, \sigma_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{dt} \in \mathbb{R}^{S_g \times PR}$, which contains only the clock offset of a specific constellation so that the adapted system model is given as follows:

<figure><img src="image_8.png" /><figcaption>**Fig. 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</figcaption></figure>

<figure><img src="image_9.png" /></figure>
$$ \begin{bmatrix} \frac{1}{G}\frac{\mathrm{d}t}{\mathrm{d}t} \\ \vdots \\ \frac{1}{G}\frac{\mathrm{d}t}{\mathrm{d}t} \end{bmatrix} = \left[ \begin{bmatrix} \frac{1}{R} \otimes \mathbf{I}_{S_1 P} \\ \vdots \\ \frac{1}{R} \otimes \mathbf{I}_{S_G P} \end{bmatrix} \cdot \begin{bmatrix} \frac{1}{G}\frac{\mathrm{d}t_c}{\mathrm{d}t} \\ \vdots \\ \frac{1}{G}\frac{\mathrm{d}t_c}{\mathrm{d}t} \end{bmatrix} + \begin{bmatrix} \frac{1}{G}\frac{\theta}{\theta} \\ \vdots \\ \frac{1}{G}\frac{\theta}{\theta} \end{bmatrix} \right] \quad (18) $$

The covariance matrix $Q$ is then blockdiagonal where each element only depends on $g_r\sigma_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. 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. 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. 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. 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. Following the outlier step, the ISB alignment is applied as previously mentioned, computing an offset value between the GPS constellation and the others.

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. 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.

### 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 psec,

<table><caption>Table 1 MGEX clk ACs' contributions and constellations</caption><thead><tr><td>Analysis center (abbr.)</td><td>Constellations</td><td>Time interval</td></tr></thead><tbody><tr><td>GRG (Loyer et al. 2018)</td><td>G+R+E</td><td>30 s</td></tr><tr><td>COD (Prange et al. 2020)</td><td>G+R+E+C+J</td><td>30 s</td></tr><tr><td>GFZ (Deng et al. 2016)</td><td>G+R+E+C+J</td><td>30 s</td></tr><tr><td>JAX (Kawate et al. 2023)</td><td>G+R+J</td><td>30 s</td></tr><tr><td>SHA (Chen et al. 2012)</td><td>G+R+E+C</td><td>5 min.</td></tr><tr><td>WHU (Guo et al. 2016)</td><td>G+R+E+C+J</td><td>30 s</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Constellation specific RMS differences between the individual AC clocks and the combined solution for AC plus constellation weighting</figcaption></figure>

and an average RMS of 45 psec. SHA was excluded from
the average computation since this AC presents high RMS
values around 122 psec. The second best is GPS, where
most of the ACs' RMS is around 50 psec while GFZ shows
an RMS of around 100 psec. We note that SHA's RMS val-
ues are around 269 psec and JAX solution shows an RMS
of 860 psec before September 2019. After configuration
changes JAX's RMS values are below 50 psec. For Bei-
Dou, the COD and WHU agree with mean RMS of around
70 psec and 78 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 psec for BDS while reaching
an RMS of 120 psec in July 2020. Again, SHA is not shown
in the plot due to high RMS values of more than 1000 psec.
For GLONASS, the lowest RMS is found for GRG with
around 70 psec; while, COD has 72 psec, WHU 88 psec,
and GFZ regularly above 200 psec. Once more, SHA val-
ues exceed 1000 psec, and for JAX, the values are at the

<figure><img src="image_11.png" /><figcaption>**Fig. 8** Weights for MGEX AC based on constellation-specific weighting</figcaption></figure>

same level as SHA until September 2019. After this date,
the JAX values are comparable with the other ACs, reach-
ing an RMS of 98 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 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.

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. Unfor-
tunately, there is a lack of documentation about the MGEX
processing from most of the ACs, making it difficult to vali-
date 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; other-
wise, the process either does not converge or yields nega-
tive 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

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10** Repro3 constellation specific RMS differences between the individual AC clocks and the combined solution</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11** Repro3 clock constellation-specific weights</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12** Daily RMS of the Ionosphere-free phase residual per AC for all stations processed in 2020 (IGS=final IGS combination; CMG=our approach)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 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)</figcaption></figure>

## 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

<figure><img src="image_17.png" /></figure>
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 psec, where COD has the lowest RMS of 27 psec and GFZ the highest of 62 psec. For Galileo, the average RMS is 53 psec, and again COD has the lowest RMS with 35 psec. The constellation with the highest average RMS is GLONASS with 68 psec. Again, COD shows the lowest with an average value of 50 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.

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 psec for GPS, 97 psec for Galileo and GLONASS with 101 psec. The agreement between both solutions is around 0.1 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.

## 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 mm and the lowest value for GPS with 9.5 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 mm for each component, except for the COD, where the height component RMS is 13 mm, while for the combined solution is around 6 mm. For CAS1, the CMG solution has the lower RMS, especially for the East component, with 6 mm, while the other solutions are above 15 mm.

The results among the solutions are similar for the station SCOR, where larger differences occur in the East component with RMS of 8.4 mm (IGS), and 3.8 mm (COD). The RMS over all stations, shows similar values for all solutions, with values of around 12 mm for the North, 20 mm for East, and 22 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 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 psec for most of the ACs. GLONASS and BeiDou show RMS, on average, under 150 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 psec, and the RMS with the IGS ACC clocks is below 100 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.

**Author contributions** The study was designed by GM and AB. The data analysis was performed by GM, AB, PS and BM. The manuscript was written by GM and AB. All authors contributed to the discussion about the content and provided comments on the manuscript.

**Funding** Open Access funding enabled and organized by Projekt DEAL. Open Access funding enabled and organized by Projekt DEAL. The PhD of Gustavo Mansur is funded by the Deutscher Akademischer Austauschdienst (DAAD, German Academic Exchange Service. Grant Number: 57381412).

## Declarations

**Conflict of interest** The authors have no relevant financial or non-financial interests to disclose.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Banville S, Geng J, Loyer S, Schaer S, Springer T, Strasser S (2020) On the interoperability of IGS products for precise point positioning with ambiguity resolution. J Geod 94(1):1–15. https://doi.org/10.1007/s00190-019-01335-w

Chen J, Wu B, Hu X, Li H (2012) SHA: The GNSS Analysis Center at SHAO. In: Proceedings of China satellite navigation conference (CSNC) 2012, vol 243 Springer Heidelberg, p213–221. https://doi.org/10.1007/978-3-642-29175-3_19

Chen G, Guo J, Wei N, Li M, Zhao Q, Tao J (2022) Multi-GNSS clock combination with consideration of inconsistent nonlinear variation and satellite-specific bias. Earth Planet Space 74(1):1–17. https://doi.org/10.1186/s40623-022-01702-6

Chen G, Guo J, Geng T, Zhao Q (2023) Multi-GNSS orbit combination at Wuhan University: strategy and preliminary products. J Geod 97(5):41

Deng Z, Fritsche M, Nischan T, Bradke M (2016) Multi-GNSS ultra rapid orbit-, clock-and EOP-product series. https://doi.org/10.5880/GFZ.1.1.2016.003

Geng J, Mao S, Pan Y (2021) Pride ckcom: Repro3 clock combination. Presented at Tour de l'IGS 1st Stop: repro3 URL https://files.igs.org/pub/resource/pubs/workshop/2021/06-Geng.pdf

Gonzalez Martinez FJ (2014) Performance of new GNSS satellite clocks. KIT Sci Publ. https://doi.org/10.5445/KSP/1000036610

Guo J, Xu X, Zhao Q, Liu J (2016) Precise orbit determination for quad-constellation satellites at Wuhan University: strategy, result validation, and comparison. J Geod 90(2):143–159. https://doi.org/10.1007/s00190-015-0862-9

Iglewicz B, Hoaglin DC (1993) How to detect and handle outliers, vol 16. ASCQ Quality Press, ISBN, p 9780873892605

Kawate K et al (2023) MADOCA: Japanese precise orbit and clock determination tool for GNSS. Adv Space Res 71(10):3927–3950. https://doi.org/10.1016/j.asr.2023.01.060

Kouba J, Springer T (2001) New IGS station and satellite clock combination. GPS Solut 4(4):31–36. https://doi.org/10.1007/PL00012863

Li X, Zhu Y, Zheng K, Yuan Y, Liu G, Xiong Y (2020) Precise orbit and clock products of Galileo, BDS and QZSS from MGEX since 2018: comparison and PPP validation. Remote Sens 12(9):1415. https://doi.org/10.3390/rs12091415
Loyer S, Banville S, Geng J, Strasser S (2021) Exchanging satellite attitude quaternions for improved GNSS data processing consistency. Adv Space Res 68(5):2441–2452. https://doi.org/10.1016/j.asr.2021.04.049

Loyer S, Perosanz F, Versini L, Katsigianni G, Mercier F, Mezerette A (2018) CNES/CLS IGS Analysis center: recent activities. Presented at IGS workshop 2018, Wuhan, China

Mansur G, Sakic P, Männel B, Schuh H (2020) Multi-constellation GNSS orbit combination based on MGEX products. Adv Geosci 50:57–64. https://doi.org/10.5194/adgeo-50-57-2020

Mansur G, Sakic P, Brack A, Männel B, Schuh H (2022) Combination of GNSS orbits using least-squares variance component estimation. J Geod 96(11):92. https://doi.org/10.1007/s00190-022-01685-y

Montenbruck O, Steigenberger P, Prange L, Deng Z, Zhao Q, Perosanz F, Romero I, Noll C, Stürze A, Weber G, Schmid R, MacLeod K, Schaer S (2017) The Multi-GNSS Experiment (MGEX) of the International GNSS Service (IGS)—Achievements, prospects and challenges. Adv Space Res 59(7):1671–1697. https://doi.org/10.1016/j.asr.2017.01.011

Prange L, Arnold D, Dach R, Kalarus S, Schaer S, Stebler P, Villiger A, Jäggi A (2020) CODE product series for the IGS-MGEX project. https://doi.org/10.7892/boris.75882.3

Rebischung P (2021) Terrestrial frame solutions from the IGS third reprocessing. In: EGU general assembly conference abstracts, EGU21–2144 @@@

Sakic P, Mansur G, Viegas E, Männel B, Schuh H (2018) Towards a multi-constellation combination: improving the IGS orbit & clock combination software for MGEX products. Presented at IGS Workshop 2018, Wuhan, China

Sakic P, Mansur G, Männel B (2020) A prototype for a multi-GNSS orbit combination. In: Proceedings of the 2020 European Navigation Conference (ENC), Dresden, Germany, https://doi.org/10.23919/ENC48637.2020.9317316

Sakic P, Mansur G, Männel B, Brack A, Schuh H (2022) An experimental combination of IGS repro3 campaign’s orbit products using a variance component estimation strategy. In: International Association of Geodesy Symposia, Springer, Berlin, https://doi.org/10.1007/1345_2022_158

Schönemann E, Becker M, Springer T (2011) A new approach for GNSS analysis in a multi-GNSS and multi-signal environment. J Geod Sci 1(3):204–214. https://doi.org/10.2478/v10156-010-0023-2

Sośnica K, Zajdel R, Bury G, Bosy J, Moore M, Masoumi S (2020) Quality assessment of experimental IGS multi-GNSS combined orbits. GPS Solut 24(2):54. https://doi.org/10.1007/s10291-020-0965-5

Springer T, Beutler G (1993) Towards an official IGS orbit by combining the results of all IGS Processing Centers. In: proceedings of the 1993 IGS workshop, Bern, Switzerland, p242–249. https://doi.org/10.1007/BF00806733

Strasser S, Mayer-Gürr T, Zehentner N (2019) Processing of GNSS constellations and ground station networks using the raw observation approach. J Geod 93(7):1045–1057. https://doi.org/10.1007/s00190-018-1223-2

Teunissen PJG, Amiri-Simkooei AR (2008) Least-squares variance component estimation

Zajdel R, Masoumi S, Sośnica K, Gałdyn F, Strugarek D, Bury G (2023) Combination and SLR validation of IGS Repro3 orbits for ITRF2020. J Geod 97(10):87. https://doi.org/10.1007/s00190-023-01777-3

Zhou W, Cai H, Chen G, Jiao W, He Q, Yang Y (2022) Multi-GNSS combined orbit and clock solutions at IGMAS. Sens 22(2):457. https://doi.org/10.3390/s22020457

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_18.png" /></figure>

**Gustavo Mansur** received his Master's degree in Geodesy Science from the Federal University of Paraná (UFPR), Brazil, in 2017. In August 2023, he completed his Ph.D. at the Technische Universität Berlin (TU Berlin). He is currently working in the Geodesy Department of the GFZ German Research Centre for Geosciences, Potsdam, Germany, where he is working on Multi-GNSS orbit and clock combination.

<figure><img src="image_19.png" /></figure>

**Andreas Brack** received his Ph.D. degree in electrical and computer engineering on the topic of high-precision GNSS from the Technical University of Munich (TUM), Germany, in 2019. He is leading the real-time GNSS group at the GFZ German Research Centre for Geosciences, Potsdam, Germany, where he is working on precise GNSS positioning, orbit and clock determination, and atmospheric sensing.

<figure><img src="image_20.png" /></figure>

**Pierre Sakic** is specialized in space and marine geodesy. He obtained his Ph.D. in 2016 on seafloor precise positioning methods applied to tectonics deformation monitoring. Between 2017 and 2021, he was a research scientist at GFZ, Potsdam, Germany, where he was in charge of the operational GNSS processing within the framework of the Galileo Geodetic Service Provider. He is now a research engineer at the Institut de Physique du Globe de Paris, where he manages the geodetic data

handling and processing of the Volcanological and Seismological Observatories.
<figure><img src="image_21.png" /></figure>

**Benjamin Männel** received his Ph.D. in Geodesy from the Institute of Geodesy and Photogrammetry at ETH Zurich, Switzerland, in 2016. He is leading the IGS Analysis Center at the GFZ German Research Centre for Geosciences, Potsdam, Germany. His main research interests are the combination of ground and space-based GNSS observations and the impact of surface loading corrections on geodetic products.

<figure><img src="image_22.png" /></figure>

**Harald Schuh** is Director of the Geodesy Department at the GFZ German Research Centre for Geosciences, Potsdam, Germany, and professor of Satellite Geodesy at Technische Universität Berlin, Germany. He has engaged in space geodetic research for more than 40 years. He was Chair of the IVS (2007–2013), President of the IAU Commission 19 “Rotation of the Earth” (2009–2012), and President of the IAG (2015–2019).

<figure><img src="image_23.png" /></figure>