Abstract
Over the past years, the International GNSS Service (IGS) has been putting efforts into extending its service by setting up and running the Multi-GNSS experiment and pilot project (MGEX). Several MGEX analysis centers (ACs) contribute by providing solutions containing not only GPS and GLONASS but also Galileo, BeiDou, and QZSS. As the current IGS combination software can only handle the orbits of one constellation at a time, it requires substantial modifications to obtain a consistent MGEX orbit product. In this contribution, we present a least-squares framework for a Multi-GNSS orbit combination, where the weights used to combine the ACs’ orbits are determined by least-squares variance component estimation. We introduce and compare two weighting strategies, where either AC-specific weights or AC and constellation-specific weights are used. An automated Z-score test is implemented yielding a common set of core satellites that are used to determine the weights. Both strategies are tested using MGEX orbit solutions for a period of two and a half years. They yield similar results with an agreement with the ACs’ orbits at the one centimeter level for GPS and up to a few centimeters for the other constellations. The 3D-RMS is generally slightly better with the AC and constellation weighting. A comparison of our combination approach with the official IGS combination using three years of GPS and GLONASS orbits shows an agreement of better than 5 mm and 12 mm for GPS and GLONASS, respectively, while the agreement of the official IGS combination with the ACs’ GPS solutions is only around \(15\,\textrm{mm}\). An external validation using satellite laser ranging shows that the mean residuals of our combined products are around \(-3\,\textrm{mm}\) for Galileo, \(6\,\textrm{mm}\) for GLONASS, \(-8\,\textrm{mm}\) for BeiDou, and \(-31\,\textrm{mm}\) for QZSS.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
1 Introduction
The International GNSS Service (IGS, Johnston et al. 2017) publishes on a regular basis GPS and GLONASS orbit and clock offset products with an accuracy of around \(2\,\textrm{cm}\) and \(75\,\textrm{ps}\) \(\textrm{RMS}\), respectively. These final products result from a combination using the orbits and clock offsets determined by the IGS Analysis Centers (ACs) as input. The motivation for a combination remains the same since the early years of the IGS: to smooth the errors and fill the possible gaps made by ACs during their independent processing, and thus provide a reliable, accurate, and unique product to the users (Springer and Beutler 1993; Beutler et al. 1995). Such combination is also a suitable way to have a reference for the comparison and validation of the ACs’ individual results (e.g., Griffiths 2019).
The method to perform the combination was developed by Springer and Beutler (1993) and described in detail by Kouba et al. (1995). Some updates were made over the years mainly to improve the clock combination and the alignment with the current ITRF (Kouba and Springer 2001). The complete history of changes in the combination strategy is described in the header of the combination summary files.
Within the IGS, the Multi-GNSS Experiment and Pilot Project (MGEX, Montenbruck et al. 2017) is working on the extension of its services regarding new systems. This project has gained attention since more receivers are capable of tracking new signals, and more constellations were declared operational in the last years, e.g., Galileo in 2016 (ESA 2016) or BeiDou in 2018 (CSNO 2018). Several ACs have recently been extending their products and provide the new constellations in their solutions.
However, so far no official combination for the new constellations is available. The IGS Analysis Center Coordinator (ACC) proposes an experimental combined MGEX solutionFootnote 1 since September 2019 based on an adaptation of the legacy strategy (Sośnica et al. 2020). The GFZ IGS team has been putting efforts to study and provide a dedicated and entirely consistent orbit and clock combination using the solutions provided by the ACs contributing to the MGEX project. Previous studies were made adapting the IGS combination software (Fritsche 2016; Mansur et al. 2020; Sakic et al. 2020), and have shown that significant modifications are needed to combine MGEX products, including an improved weighting strategy adapted to the different constellations. In the current software, the combination for GPS and GLONASS is done completely independent for each constellation, which is not ideal in terms of interoperability. A simple adaptation of the existing software may therefore not be suited to handle the different constellations with different quality of their orbit solutions.
To this end, we propose a unified multi-GNSS orbit combination approach based on variance component estimation (VCE). VCE has been extensively used in different applications in geodesy over the years. Some noticeable examples can be found in Koch and Kusche (2002), where VCE is used for geopotential regularization, and in Kusche (2003) within the weight determination for GOCE (Gravity Field and Steady-State Ocean Circulation Explorer) gravity recovery solutions. In the context of combining Terrestrial Reference Frames (TRF), Bähr et al. (2007) presents a comparison between three VCE methods, and more recently Kobel et al. (2019) implemented a VCE approach to combine the orbits of the Low Earth Orbit (LEO) satellite Sentinel-3. In this contribution, we make use of least-squares variance component estimation (LSVCE), originally developed by Teunissen (1988). The method is derived and described in detail in Amiri-Simkooei (2007) and Teunissen and Amiri-Simkooei (2008). Applications of LSVCE for GPS data processing are provided in Amiri-Simkooei et al. (2009) and Amiri-Simkooei et al. (2013). The method is flexible and relies on the well-known least-squares theory. We present the application of LSVCE to an orbit combination, where we study two different models: (1) Considering AC-specific orbit variances and (2) introducing AC plus constellation-specific orbit variances. The final combined orbits are derived as a weighted mean of the Helmert-transformed orbits of the ACS, where the weights are the normalized inverse estimated variance components, so that the sum of the weights equals one.
The overall structure of this manuscript contains five sections. Section 2 provides a brief explanation of the steps implemented to combine the orbits, whereas in Sect. 3 the mathematical least-squares framework of our orbit combination is introduced. Section 4 summarizes the experiments, compares the two methods, and shows a validation of the combined orbits using satellite laser ranging (SLR) and precise point positioning (PPP). A summary and conclusions are given in Sect. 5.
2 Combination steps
The task of the orbit combination is to transform the orbit products of different ACs to a single combined orbit product. In short, the current IGS combination procedure is as follows. It starts with a simple mean of the orbit solutions from all ACs. Then, seven parameter Helmert transformations between the individual orbits and the mean orbit are performed, and the ACs’ orbits are aligned to the mean via these transformations. A weight is derived for each AC from the differences of these transformed orbits and the mean orbit, and a new weighted mean orbit is computed. The same steps are repeated once, starting with the weighted mean orbit, finally yielding the combined orbit product.
In this contribution, we present a combination approach that is based on the least-squares principle. The combined orbits are then again a weighted mean of the orbits of the ACs, where the weights are proportional to the inverse variances of the ACs’ solutions. These unknown variances are estimated by means of LSVCE (Amiri-Simkooei 2007).
The orbit combination process is implemented in Python within the framework of the GeodeZYX toolbox (Sakic et al. 2019). The five main combination steps listed below are illustrated in Fig. 1. Detailed mathematical derivations are provided in Sect. 3.
-
(1)
We start with a pre-check to identify and exclude satellites that appear in only one of the ACs’ solutions or satellites with missing epochs. A first rough exclusion procedure is performed, in which the satellite coordinates are compared between ACs. If one AC provides an orbit for a satellite that is very far from the other ACs (e.g., threshold set at \(500\,\textrm{m}\)), this solution is already excluded at this stage. A mean orbit is computed using the pre-checked orbits.
-
(2)
The second step consists of a Helmert transformation between this mean orbit and the ACs’ solutions, which provides the seven parameters per AC that are used to align the orbits.
-
(3)
In the third step, a set of core satellites is defined. The idea is to select a subset of the orbits that are used for the VCE. The set of core satellites ensures that the satellites for computing the variance components are the same for all ACs and that outliers are detected and excluded to get proper estimates of the weights. The outlier detection is done by applying the Modified Z-Score (Iglewicz and Hoaglin 1993) to the radial, along-track, and cross-track (RAC) components of each set of AC coordinates for all satellites. If one of the components is identified as an outlier for a satellite, this satellite is excluded from the set of core satellites. The detection threshold is chosen less strict for BeiDou and QZSS since the quality of the ACs’ orbits for these systems is still not on the level of GPS, GLONASS, and Galileo.
-
(4)
The fourth step is the estimation of the variance components as described in Sect. 3.3, using only the above defined set of core satellites.
-
(5)
With the variance components determined, the final step is the orbit combination. As described in Eq. (12) in Sect. 3.2, the inverse variances are normalized and then applied as weights to the ACs’ orbits. The orbit combination is also provided for satellites that are not in the set of core satellites, where, if a solution is missing or excluded for one or more ACs for a specific satellite, the normalization is based only on the variances of the ACs that are included.
This process is repeated as mentioned at the end of Sect. 3.1, starting with the Helmert transformation that is computed with respect to the combined orbits of the previous iteration. The iteration is terminated once the 3D-RMS difference of the orbits of two subsequent iterations is below a threshold of \(1 \,\textrm{mm}\).
During the processing, we can define three groups of satellite orbits. (1) The first group contains the ACs’ orbits of the core satellites that pass the outlier detection and are common among all ACs. They are the ones defining the weights of the combination. With constellation-specific weighting (cf. Sect. 3.4), the set of core satellites is defined per constellation, so that ACs which provide only some constellations can still be considered. ACs with too few common satellites can be excluded in order to have enough core satellites. (2) The second group is used to compute the Helmert transformations. It includes the orbits of the core satellites of group (1), but also of the satellites that are not available at all ACs. Satellites that are excluded as outliers for some of the ACs when defining the set of core satellites are also included for the remaining ACs. (3) The third group of satellite orbits is the largest group, which is used to compute the combined orbits. In addition to the second group, it also contains orbits of satellites that are detected as outliers for all ACs. These satellites are included in order to keep the combined solution as complete as possible.
3 Mathematical framework
In this section, the mathematics for the combination framework of Sect. 2 are derived. In short, it is the solution to the problem of a joint least-squares estimation of the combined orbits and variance components for the input orbits.
3.1 System model
The input of the orbit combination is the orbits of the ACs. Let the Earth-fixed coordinates of satellite \(s \in \{ 1,\dots , S\}\) provided by AC \(r \in \{1, ..., R\}\) at the time t be given by \({\underline{x}}_{r}^{s}(t) \in {\mathbb {R}}^{3}\). We define
and
with \(t_p\), \(p \in \{1, ..., P\}\), the time at epoch p. The unknown combined orbit vector \({\underline{x}}_{\textrm{c}} \in {\mathbb {R}}^{3SP}\) is defined according to \({\underline{x}}_{r}\). The relation between the ACs’ orbits and the combined orbits can be described by means of a linearized Helmert transformation (Watson 2006) as
where the AC-specific Helmert transformation matrices \({H}_{r}\) are in the form of
with \(\theta _x\), \(\theta _y\), \(\theta _z\), and \(\mu \) the rotation and scale parameters, and \(\varDelta {\underline{x}}_r\) the translation parameters, that align the ACs’ orbits and the combined orbits, and \({\underline{\varepsilon }}_r^s(t) \in {\mathbb {R}}^{3}\) an additive noise term. Analogous to Eq. (2), \({\underline{\varepsilon }}_r \in {\mathbb {R}}^{3SP}\) is composed of all \({\underline{\varepsilon }}_r^s(t_p)\). In order to facilitate our derivations, we apply the following invertible transformation to the input orbits \({\underline{x}}_{r}^{s}\left( t \right) \):
Note that all information is kept through this transformation, so that the transformed orbits \(\bar{{\underline{x}}}_{r}^{s}\left( t \right) \) can equivalently be used instead of the original orbits. The final system model can now simply be written as
where \({\underline{1}}_{R}\) is an R-vector of ones, \(I_{3SP}\) an identity matrix of size 3SP, \(\otimes \) the Kronecker product, and \(\bar{{\underline{\varepsilon }}}_r\) are the transformed versions of \({\underline{\varepsilon }}_r\).
For the stochastic model, we assume that
i.e., there is an AC-specific variance \(\sigma _r^2\) that is common for all coordinates. This assumption is a simplification, since, for instance, the positions of each satellite are derived from the same set of orbital parameters within each AC, and the orbits of different satellites are connected through the receiver clock parameters in a network solution. Including such information would, however, require the covariance matrices of the ACs’ orbit estimation, which are not accessible. Assuming uncorrelated solutions between ACs, we therefore have
where \(\sum _{r=1}^{R} \sigma _{r}^2 Q_{r}\) is positive definite, \(\sigma _r^2\) are the variance components, and \(Q_{r} = \textrm{diag}({\underline{c}}_r) \otimes I_{3SP}\) is a known sparse block matrix, with \({\underline{c}}_r\) the rth canonical unit vector in \({\mathbb {R}}^R\). These AC-specific variance components \(\sigma _r^2\) lead to AC-specific weights for the orbit combination, cf. Sect. 3.2. An extension allowing for AC plus constellation-specific weights is presented in Sect. 3.4.
We can now solve the system model as described in (6) and (8) for the unknown combined orbits \({\underline{x}}_{\textrm{c}}\) and the variance components \(\sigma _r^2\), \(r\in \{1\dots ,R\}\). The estimation follows the LSVCE framework presented in Amiri-Simkooei (2007) and Teunissen and Amiri-Simkooei (2008). It is noted that for the computation of the Helmert transformation parameters in Eqs. (3)–(5), the unknown combined orbits \({\underline{x}}_{\textrm{c}}\) are required so that the above estimation problem is solved iteratively, where the combined orbits can be initialized with a simple mean and then updated in each iteration.
3.2 Combined orbit solution
We recall from Eqs. (6)–(8) that our system is of the form \(\bar{{\underline{x}}}=A{\underline{x}}_c + \bar{{\underline{\varepsilon }}}\), so that the weighted least-squares solution \(\hat{{\underline{x}}}_{c}\) for the combined orbit is given by
Using the mixed-product property of the Kronecker product and that the transposition is distributive on the Kronecker product (properties recalled by Teunissen and Amiri-Simkooei 2008 and proven by, e.g., Rao 1973), we have
and
With Eqs. (10) and (11), the combined orbit solution from Eq. (9) can be rewritten as
which shows that the combined orbits are indeed a weighted mean of the ACs’ orbits. The estimation of the variance components \(\sigma _r^2\) is shown in the next section.
3.3 Variance component solution
Following the LSVCE theory presented by Amiri-Simkooei (2007) and Teunissen and Amiri-Simkooei (2008), let B be a matrix whose columns span the null space of \(A^\textrm{T}\) so that \(A^\textrm{T}B = 0\). It follows that the random vector of misclosures \({\underline{t}}=B^{T} \bar{{\underline{x}}}\) is zero mean and uncorrelated with \(\hat{{\underline{x}}}_\textrm{c}\). There is a one-to-one correspondence between \(\left[ \hat{{\underline{x}}}_\textrm{c}^\textrm{T}, {\underline{t}}^\textrm{T} \right] ^\textrm{T}\) and \(\underline{\bar{x}}\), cf. Eq. (3) in Teunissen and Amiri-Simkooei (2008), so that instead of \(\underline{\bar{x}}\) we can equivalently use \(\hat{{\underline{x}}}_\textrm{c}\) and \({\underline{t}}\) to estimate the unknowns. As \(\hat{{\underline{x}}}_\textrm{c}\) is already the least-squares estimate of the parameter vector \({\underline{x}}_\textrm{c}\) and does not leave any redundancy, and as it is uncorrelated with \({\underline{t}}\), the variance components \(\sigma ^2_r\) are estimated only from \({\underline{t}}\). With \({\underline{t}}\) zero mean, the covariance matrix of \({\underline{t}}\) is given by \(\textrm{E}\left\{ {\underline{t}}{\underline{t}}^\textrm{T} \right\} = B^\textrm{T} Q B = \sum _{r=1}^R \sigma _r^2 B^\textrm{T}Q_r B\), which can be considered as a matrix observation equation for estimating \(\sigma _r^2\). Using the half-vectorization operator \(\textrm{vh}(\cdot )\) on \({\underline{t}}{\underline{t}}^\textrm{T}\) that transforms a matrix into a vector whose components are the \(\frac{n(n+1)}{2}\) entries on and below the main diagonal (Turkington 2005), thereby eliminating the duplicate entries of a symmetric matrix, a standard linear observation model is obtained as
with \({\underline{\sigma }}=\left[ \sigma _1^2,\dots ,\sigma _R^2 \right] ^\textrm{T}\) and \(A_\textrm{vh} = \left[ \textrm{vh}\left( B^\textrm{T}Q_1 B\right) ,\dots ,\textrm{vh} \left( B^\textrm{T} Q_R B \right) \right] \). The solution of the variance components \(\hat{{\underline{\sigma }}}\) defining the weights in Eq. (12) is given by
where \(W_{\textrm{vh}}\) is a weight matrix. With \(N={A}_{\textrm{vh}}^{\textrm{T}} W_{\textrm{vh}} {A}_{\textrm{vh}}\) and \({\underline{l}}={A}_{\textrm{vh}}^{\textrm{T}} W_{\textrm{vh}} \textrm{vh}\left( {\underline{t}} {\underline{t}}^{\textrm{T}}\right) \), we can write
Using the inverse covariance matrix \(Q_\textrm{vh}\) of \(\textrm{vh}\left( {\underline{t}} {\underline{t}}^{\textrm{T}}\right) \) as weight matrix, the elements of the normal-matrix N and the vector \({\underline{l}}\) are derived as (cf. Eqns. 4.85 and 4.86 in Amiri-Simkooei 2007)
with \(Q_{t}=B^{\textrm{T}} Q B\) and \({\text {tr}}(\cdot )\) the trace operator. Since \(Q_t\) depends on the variance components \({\underline{\sigma }}\), (15) is again solved iteratively within each iteration of solving Eq. (6), where \(Q_t\) is computed with the estimates of the variance components \(\sigma _r^2\) from the previous iteration.
In order to solve Eq. (15), we therefore need to find a matrix B, whose columns are a basis of the null space of \(A^\textrm{T}\) such that \(A^\textrm{T}B = 0\), as already mentioned. There are infinitely many choices for B, all of which will lead to the same result \(\hat{{\underline{\sigma }}}\). Our goal is to find a matrix B that leads to an efficient implementation avoiding large scale matrices that could be obtained for daily combinations with many satellites and ACs. With \(A=\left[ {\underline{1}}_{R} \otimes I_{3SP} \right] \), one solution for B is given by
with the differencing operator
The property \(A^\textrm{T}B = 0\) follows from \(D^\textrm{T}{\underline{1}}_R = {\underline{0}}\) and the mixed-product property of the Kronecker product. The vector of misclosures \({\underline{t}}=B^{\textrm{T}} \bar{{\underline{x}}}\) in Eq. (17) contains the between AC differenced orbits \({\underline{x}}_r - {\underline{x}}_1\), \(r \in \lbrace 2,\dots ,R \rbrace \), and \(\hat{{\underline{\sigma }}}\) is a function of only those differences.
With \(Q_r = \textrm{diag}({\underline{c}}_r) \otimes I_{3SP}\), we have
and matrix \(Q_t = B^\textrm{T} Q B\) required in Eqs. (16) and (17) can be rewritten as
Consequently, Eqs. (16) and (17) are given by
The blocks \({D}^{\textrm{T}} {\text {diag}}\left( {\underline{c}}_{r} \right) {D}\) used in this representation are only of dimension \(R-1 \times R-1\) and therefore suitable for an efficient implementation. Using Eqs. (22) and (23), the variance components \({\underline{\sigma }}\) are estimated, cf. Eq. (15), and used to compute the combined orbits, cf. Eq. (12).
3.4 Introducing constellation-specific weights
So far, we assumed that there is one variance component for each AC, see Eqs. (7) and (8), which translates into a weight per AC in the orbit combination in Eq. (12). However, this assumption might be too simple for a multi-GNSS solution, in which the orbits of different constellations might be of different quality. We can therefore extend the above model by introducing a variance component \(_g\sigma _{r}^2\), \(g\in \{ 1,\dots , G\}\), for each AC and constellation, with G the number of constellations. We define \(_g\bar{{\underline{x}}} \in {\mathbb {R}}^{3S_gPR}\) as a vector similar to \(\bar{{\underline{x}}}\) in Eq. (6) that only contains all coordinates of a specific constellation with \(S_g\) satellites. The system model then reads
where \(_g{\underline{x}}_\textrm{c}\) and \(_g\bar{{\underline{\varepsilon }}}\) are the constellation-specific versions of \({\underline{x}}_\textrm{c}\) and \(\bar{{\underline{\varepsilon }}}\). Assuming no correlations between different constellations, which is a natural consequence of the assumption in Sect. 3.1 that no correlations are present between the coordinates of different satellites, the covariance matrix Q of the noise vector \(\left[ _1\bar{{\underline{\varepsilon }}}^\textrm{T}, \dots , _G\bar{{\underline{\varepsilon }}}^\textrm{T} \right] ^\textrm{T}\) is blockdiagonal
where \(_gQ\) only depends on \(_g\sigma _r^2\), with \(r\in \{1,\dots ,R\}\). Therefore, the sets of parameters \(_g{\underline{x}}_\textrm{c}\) and the variance components \(_g\sigma _r^2\) are mutually exclusive between the constellations, and the measurements \(_g\bar{{\underline{x}}}\) are uncorrelated so that the joint orbit combination and LSVCE can be computed separately for each constellation. It is important to note that this is only true within each iteration of the combination process and that the constellations are still connected through the common Helmert transformation in Eq. (3).
4 Results
The two weighting approaches described in Sect. 3 are tested using products that are provided to the IGS by the ACs. We divide these products into two categories:
-
LegacyFootnote 2: GPS and GLONASS products provided by the ACs within the regular IGS operational processing, for which the IGS provides a separate combination (Kouba 2009; Johnston et al. 2017)
-
MGEXFootnote 3: products provided by the ACs to the Multi-GNSS Experiment and Pilot Project (Montenbruck et al. 2017)
Table 1 shows the ACs that contribute to the two processing lines and the constellations contained in the solutions. Some ACs provide both legacy and MGEX solutions. Except for GRG, these solutions are based on different processing setups and are therefore not identical. More information is provided in Mansur et al. (2020).
Two periods of processing are used for our analysis. For MGEX, we use orbit products between GPS week 1982 (January 2018) and GPS week 2116 (July 2020). GPS week 1930 (January 2017) until GPS week 2116 (August 2020) is considered for the legacy products.
4.1 Differences between individual AC orbit solutions and combined orbits
The results from the MGEX orbit combination are compared to the ACs’ individual solutions. Figure 2 shows the time series of the RMS orbit differences with the AC plus constellation weighting approach. The RMS values are computed based on the formula in Kouba et al. (1995). The smallest RMS values are obtained for GPS, and they are better than \(25\,\textrm{mm}\) for all ACs. QZSS shows the highest RMS values, which can reach around \(80\,\textrm{mm}\). This behavior is also noted by Sakic et al. (2020) who used an adaptation of the IGS legacy software to combine the MGEX orbits. Similar results are provided by the ACC where another adaptation of the legacy software is usedFootnote 4. Moreover, for QZSS, especially in 2018, there are days without a combined solution due to missing data or high discrepancy between the satellite coordinates provided by the ACs yielding solutions with high RMS values. These discrepancies can be related to the incomplete solar radiation models applied by the ACs to this constellation (Montenbruck et al. 2013). The days without a solution or with unreasonably high RMS values are not presented. One can note the gaps in Fig. 2. The RMS analysis shows similar characteristics for some ACs for certain constellations, see, for example, the RMS values of GLONASS for GRG and CODE, or for GFZ and WHU from 2018 until early 2019. TUM has the most discrepant RMS values. Similar results are obtained with the AC-only weighting approach, where the largest difference is obtained for QZSS with RMS values of around \(90\,\textrm{mm}\) instead of \(80\,\textrm{mm}\) (not shown).
Figure 3 shows the average of the daily RMS results per AC over the entire test period for both approaches. It is noted that TUM is excluded from the AC-only weighting approach. As described in Sect. 2, our algorithm defines a list of core satellites (satellite group 1), which are satellites provided by the majority of the ACs. Since TUM products do not contain GPS and GLONASS, this AC is excluded from the combination. However, a comparison of the TUM orbits with the combination is still provided. Comparing both strategies, only slight differences at the millimeter level are noted. The most significant differences are observed for QZSS with 13, 11, and \(8\,\textrm{mm}\) for CODE, WHU, and TUM, respectively. The differences are less than \(5\,\textrm{mm}\) for GPS, followed by GLONASS, BeiDou, and Galileo.
Given that the Helmert transformation connects the orbits across constellations, as shown in Eq. (3), we ensure small differences between the seven parameters obtained by our two strategies. On average, there are only differences in the millimeter level. These differences are related to the weighting scheme since, within both approaches, distinct mean orbits are computed and are used in the Helmert transformation.
4.2 Weights
Figures 4 and 5 show the weights used in the combination for the two approaches, and we can observe changes and trends of the ACs’ weights over time. For the AC plus constellation weighting in Fig. 4, GPS has the most consistent AC weights over the test period, followed by GLONASS and Galileo, while the weights for BeiDou and QZSS vary between the ACs. We might attribute this behavior to the number of satellites in the newer system, the different orbit types, or changes in the ACs’ processing routines. We note that in our processing we only consider the MEO satellites of the BeiDou constellation.
For the GPS constellation, WHU has the largest weight, followed by GFZ, with an average of \(28\%\) and \(21\%\), respectively. We observe that since May 2020, CODE’s weight has increased by \(5\%\). Similarly, for GLONASS the most dominating weights are obtained for GFZ and WHU, and again CODE’s weight increased at the end of the test period. For Galileo, we can see an increase in the GRG weight at the end of 2018. We relate this to the time when GRG activated ambiguity fixing in their solution (Loyer et al. 2018). In addition, for this constellation, CODE is the most stable AC in the years 2018 and 2019, but with a larger change in 2020, reducing the average weight of around \(40\%\) to \(25\%\). As mentioned before, BeiDou and QZSS have overall larger variations in the weights. For the Chinese system, the values are more or less evenly distributed among the ACs at the beginning of 2018. From April 2018, SHA had a big increase until January of 2019. In 2020, GFZ reaches an average weight of \(50\%\) while WHU and SHA’s weights decreased. From June 2020, CODE’s BeiDou orbits are excluded from the combination. As mentioned before, our algorithm defines a set of core satellites commonly available at the majority of the ACs, on which the computation of the weights is based. At this period, other ACs also provide BeiDou-3 solutions, which allows for a bigger set of core satellites without CODE; therefore, the CODE BeiDou-2 only solutions were not used. For QZSS, CODE got the major weight, with on average of more than \(30\%\), followed by GFZ. Until July 2018, TUM had a significant weight, decreasing after this period while WHU’s weights are more stable in 2019 and diminishing at the beginning of 2020. Generally, an increase of the weight of an AC goes along with a reduction of the RMS, see Figure 2, for example, for GRG and the Galileo constellation in the end of 2018.
The AC-only weights shown in Fig. 5 reveal that the values are more evenly distributed over the time of processing, but with clear advantages for GFZ and WHU. There is an increase of around \(5\%\) for CODE and GRG, especially in the first half of 2020. However, one has to keep in mind that the weights are defined mainly based on GPS, GLONASS, and Galileo, since the satellites of these constellations are most represented in the core set. For both strategies, the number of iterations needed to achieve the convergence is generally small. For the ACs only weighting, where satellites of different constellations are used in the core, it takes no longer than four iterations on average, whereas for the AC plus constellation weighting on average three iterations are needed for GPS, GLONASS and Galileo, and four for BeiDou and QZSS. The weight differences are small after the second iteration, varying around \(1\%\).
4.3 Validation
To provide a validation of our combination, three different validation steps have been performed. Firstly, we combine the orbit products provided within the IGS legacy chain using the AC-only weighting approach and compare the results to the official IGS combination. Further, the MGEX-based combined orbit solutions with both weighting strategies are validated using the SLR technique, and finally the consistency of the combined orbits is verified by looking at the phase residuals of a PPP analysis.
4.3.1 Combination of IGS legacy products
Since only the constellations GPS or GLONASS are processed, we choose the AC-only weighting approach. Figure 6 shows the GPS-only RMS orbit differences between our combined solution, the ACs’ orbits, and also the legacy IGS rapid and final orbits as reference. Similar to the ACC’s official IGS results (Villiger and Dach 2020), an agreement of around \(15\,\textrm{mm}\) is shown in our solution. Compared to the IGS rapid and final orbits, the RMS is \(6\,\textrm{mm}\) and \(4\,\textrm{mm}\), respectively, showing a very good agreement. Regarding the weights, CODE, on average, has the highest weight of \(15\%\), followed by NGS with \(12\%\) and JPL with \(11\%\). In addition, we observe that MIT’s influence increased in 2020 when their weight raised from around \(6\%\) in 2019 to \(15\%\) in 2020, which can be related to the new release of the GAMIT/GLOBK 10.71, that follows the recommendations for the IGS’s repro3 (MIT 2020). For GLONASS, the RMS difference to the IGS final solution is around \(12\,\textrm{mm}\). The weights are more evenly distributed at the beginning of 2019, when ESA’s weight decreased from an average of \(55\%\) to \(30\%\), while CODE and GRG show an increase of \(10\%\) each.
4.3.2 SLR-based validation
SLR measurements to the GNSS satellites provided by the International Laser Ranging Service (ILRS) establish an independent external validation of their orbits (Pearlman et al. 2019). It is a tracking technique that depends on short laser pulses not susceptible to ionospheric delays, ambiguities, and other effects like tropospheric delays caused by water vapor that introduce errors in the observations. The SLR residuals are directly obtained from the difference between the geometric distance that comes from the GNSS orbits and the ITRF station coordinates and the SLR measurements, after applying proper corrections (Sośnica et al. 2015; Dach et al. 2019). Currently, the satellites of Galileo, QZSS, GLONASS, and BeiDou are equipped with retro-reflectors and tracked by the ILRS. The number of SLR observations depends on the constellation and the ILRS priority list, with GLONASS and Galileo having in general more measurements than BeiDou and QZSS (Sośnica et al. 2020).
The SLR residuals are computed using GFZ’s EPOS.P8 software (Gendt et al. 1998; Uhlemann et al. 2015) and confirm the similarity between both strategies, as the differences in the SLR residuals between both combination strategies are at the millimeter level. We show the results per space vehicle number (SVN) from the AC plus constellation weighting approach in Fig. 7. The constellation with the smallest standard deviation (STD) of the SLR residuals is Galileo with an average STD of \(20 \pm 6\,\textrm{mm}\), followed by GLONASS with \(28 \pm 8 \,\textrm{mm}\), whereas \(68 \pm 35\,\textrm{mm}\) and \(131 \pm 34 \,\textrm{mm}\) are obtained for BeiDou and QZSS, respectively.
In general, the offsets of Galileo are at the millimeter level, and the highest values are for the in-orbit validation (IOV) satellite group with negative offsets. However, the satellites with the minimum STD of \(14\,\textrm{mm}\) are E102 (IOV) and E209 from the full operational capability (FOC) satellite group. The highest values of the STD are obtained for the satellites E201 and E202 with eccentric orbits with 45 and \(35\,\textrm{mm}\). For GLONASS, large positive offsets are observed for R855, R856, R858 of the block M+ with 36, 27, and \(33\,\textrm{mm}\), and for block M satellite R857 with \(39\,\textrm{mm}\). This behavior can also be noted in the preliminary GFZ repro3 results presented in Männel et al. (2020a) and the analysis of the experimental ACC combination shown in Sośnica et al. (2020). Meanwhile, the lowest value of the standard deviation for GLONASS is obtained for R802 (K1B) with \(15\,\textrm{mm}\) and the highest for R730 (M) with \(47\,\textrm{mm}\). Unlike the other systems, for BeiDou, the smallest offsets and STDs are observed using the AC-only weighting approach. Using this strategy, the lowest offsets are for the satellites from BeiDou-3 with \(7\,\textrm{mm}\) for C207 (MEO) and \(1\,\textrm{mm}\) for C208 (MEO). Comparing the STD for these two satellites, the values are \(80\,\textrm{mm}\) and \(125\,\textrm{mm}\) for C207 and C208 when using the AC plus constellation weighting, whereas they are \(49\,\textrm{mm}\) and \(34\,\textrm{mm}\) when using the AC-only weighting. The QZSS constellation has the most significant differences between both strategies. For all the satellites of this constellation, the AC plus constellation weighting approach results in smaller offset and STD values. The minimum STD is observed for J001 (IQ) with \(98\,\textrm{mm}\) and the maximum for J003 (IIG) with \(178\,\textrm{mm}\). We note that for J002 (IIQ), we have the biggest difference of the STDs between both strategies with \(41\,\textrm{mm}\). The smallest offset is observed for J002 with \(23\,\textrm{mm}\), and the only satellite with negative offset values is J001.
In Fig. 8, we compare the mean value of the SLR residuals of our combination (CMG) with the ACs’ orbits. This evaluation is done using the same set of satellites for all ACs during the same period, which is why the results can differ significantly from Fig. 7, in particular for QZSS. The constellation with the lowest residuals is Galileo, where the CMG mean value is \(-3\,\textrm{mm}\) with an STD of \(6\,\textrm{mm}\) and the AC with the smallest value is WHU with a mean of \(1\,\textrm{mm}\) and an STD of \(6\,\textrm{mm}\). For Galileo, the STD values are similar for all ACs, where only SHA has a bigger absolute value of the mean with around \(-34\,\textrm{mm}\). For GLONASS, our combination has a mean of \(6\,\textrm{mm}\) and an STD of \(14\,\textrm{mm}\), and SHA has the lowest absolute value of the mean with \(1\,\textrm{mm}\) and an STD of \(15\,\textrm{mm}\). For BeiDou, the best AC is WHU with \(-8\,\textrm{mm}\) mean and \(19\,\textrm{mm}\) STD. The highest absolute values of the mean SLR residuals are obtained for QZSS, where the combination has a mean of \(-31\,\textrm{mm}\) and an STD of \(12\,\textrm{mm}\). For this constellation, JAX has the lowest values with a mean of \(-3\,\textrm{mm}\) and an STD of \(38\,\textrm{mm}\). We note that the STDs of the ACs and the combination are similar. Moreover, the combination does not have the lowest residuals for any of the constellations, which is what one would expect as the weighted mean of the ACs’ orbits should average out the orbital errors. Due to potential systematic errors, it is, however, possible that the orbits of one AC are more accurate than the combined orbits. The purpose of the combination should therefore also be seen in terms of providing the most complete and stable solutions and not necessarily always the solution with the highest accuracy.
4.3.3 PPP validation
The primary purpose of the combined GNSS orbits lies in the field of high precision GNSS applications. A main technique in this context is PPP (Zumberge et al. 1997), which can provide sub-centimeter level positioning information and is used in various fields in academia and industry (Bisnath and Gao 2009). Most important, the PPP quality relies on the availability and consistency of precise satellite orbit and clock solutions.
We perform a PPP analysis using GFZ’s EPOS.P8 software. In order to provide a fair comparison, the same constellations are used for all ACs. We therefore focus on GPS, Galileo, and GLONASS so that five ACs can be compared. Daily coordinates are computed using the ACs’ orbit and clock products. The combined orbits are used together with combined clocks that are also based on the above described VCE procedure.
The average daily RMS residuals of the ionosphere-free phase combination using ten IGS stations during GPS week 2096 (March 2020) are shown in Fig. 9 for each constellation. The averaged residuals of COD, GRG, SHA, and WHU are all below \(10\,\textrm{mm}\). The GLONASS residuals are slightly larger than the ones of GPS and Galileo, which agrees with the SLR analysis (cf. Fig. 8). Residuals derived from the GFZ products are clearly larger due to problems with their inter-system clock biases (solved on day of the year 161 of 2021). The GFZ clock products are therefore excluded from the clock combination. The phase residuals with the combined CMG products are on the same level as the individual ACs with values of less than \(10\,\textrm{mm}\).
Furthermore, a PPP processing is performed using our combined solution from Sect. 4.3.1 and the IGS final products, for GPS only. The residuals of the ionosphere-free phase combination in both cases are on the same level with an RMS value of \(8\,\textrm{mm}\). In addition, Fig. 10 shows the comparison of the RMS coordinate differences between the PPP results using IGS and CMG and the final IGS combined station coordinates for four stations, where we can see that the two solutions have similar values for the three components.
The PPP analysis demonstrates that the combined products are suitable for high precision GNSS applications and that the assumption of constant inter-system biases is not violated with the constellation-specific weights.
5 Conclusions
In this study, we present and compare two approaches developed to combine GNSS orbits, where AC or AC and constellation-specific weights are used. Section 3 describes the assumed functional and stochastic models, based on which a least-squares framework is formulated that jointly estimates the combined orbits and the weights used in the combination. Aiming for a reduced computational load, Eqs. (5) and (8) are derived to be easily implemented. Among the models’ derivations, one important aspect is to find a basis matrix B for the null space of \(A^\textrm{T}\). Given the large dimensions of the matrices when considering daily solutions with many satellites and several ACs, a direct computation may require heavy usage of the computer memory. Therefore, simplifications are needed for operational usage. With the Helmert transformation applied in Eq. (5), a simple and structured form of the system matrix A is found, enabling a closed form solution for B. The structure of this matrix B again allows for an efficient implementation of the least-squares estimator of the variance components, which define the weights in the orbit combination, in which only matrices of small dimension have to be handled, see Eqs. (16) and (17). The method is initially derived for an AC-only weighting. It was shown that the AC plus constellation weighting scheme can effectively be reduced to the methods of the AC-only scheme. In this AC plus constellation weighting strategy, the weights and combined orbits can be solved individually for each constellation and are only connected by a common Helmert transformation in each iteration. Another important topic of the proposed combination approach is the definition of a set of common core satellites, which are used to compute the weights of the combination. We propose a method based on the Modified Z-Score test that only selects satellites without outliers in the RAC components.
Both combination methods are tested using the MGEX ACs’ orbits for a test period from January 2018 until August 2020. The analysis shows that using AC-only or AC plus constellation weighting leads mostly to similar results where a significant RMS difference can be found only for QZSS. As expected, GPS has the smallest RMS with less than \(25\,\textrm{mm}\), while GLONASS and Galileo are below \(50\,\textrm{mm}\). For BeiDou and QZSS, the RMS reaches \(80\,\textrm{mm}\). Concerning the weights determined by the procedures, the AC plus constellation weighting scheme shows interesting results where one can notice the different behavior of the ACs for each system over time. For example, in the Galileo weights, there are significant changes for GRGS, CODE, and GFZ. Not just for the European system but likewise for BeiDou and QZSS, significant weight variations are observed, especially for WHU, GFZ, and CODE.
A validation using products from the legacy IGS chain is performed. The RMS difference between the GPS-only combined orbit and the official IGS final and rapid orbits show an agreement of better than \(4\,\textrm{mm}\) and \(6\,\textrm{mm}\). These small differences for the standard GPS-only case confirm the validity of the new combination method. As an external validation, the orbits generated by the two approaches are tested using SLR, assessing if the combined orbits can adequately represent the satellites’ positions. Within the SLR validation, Galileo showed the smallest offsets, while QZSS has more discrepant values. When comparing the SLR residuals of each individual AC and the combination, the result shows that the combination adapts the dispersion of the ACs’ solutions, and, therefore, is not the solution with the lowest residuals. A PPP validation using the combined orbits results in sub-centimeter GNSS phase residuals, which demonstrates the consistency of the combined orbit product and its suitability for high precision GNSS applications.
The results support the application of LSVCE to determine the weights for an orbit combination. The implementation of the presented algorithm is sufficiently efficient to generate products on an operational basis, with a preliminary computational time of around ten minutes for daily solutionsFootnote 5. Judging the two strategies, even if the results are mostly similar, the AC plus constellation weighting reveals more information about the behavior of the ACs’ solutions through the distribution of the weights for each constellation, so that, for instance, problems and changes in the processing of each AC are easier visible. With more constellations becoming fully operational, several ACs are adapting their software to handle the new systems. It is therefore important to verify that an adaptation done for one system does not affect the others. Generally, from the experiments, one can see that the MGEX solutions submitted to the IGS are getting more stable over the years, and the ACs are putting efforts into providing products, including all constellations available. So far, the alignment with the ITRF is not implemented in the algorithm since there are no SINEX files provided by the ACs. A future addition will include the reference frame alignment approach from Mansur et al. (2020) in the combination. The application of the VCE-based orbit combination to the third reprocessing campaign of the IGS (repro3) is presented in Sakic et al. (2022). Summarizing, the presented combination approach was shown to be suitable to provide reliable and accurate combined multi-GNSS orbits on a regular basis, which are an important product for GNSS users.
Data availability
The GNSS products and observation data provided by the International GNSS Service (IGS) are available at https://cddis.nasa.gov/archive/gnss/. SLR observation data provided by the International Laser Ranging Service (ILRS) are available at https://cddis.nasa.gov/archive/slr/. The combined orbit products described in this study can be obtained directly from the authors upon request.
Notes
On a single computer equipped with a CPU Intel i5-4595 (4) @ 3.700GHz.
References
Amiri-Simkooei AR (2007) Least-squares variance component estimation - theory and GPS applications. PhD thesis, Delft University of Technology
Amiri-Simkooei AR, Teunissen PJG, Tiberius CCJM (2009) Application of least-squares variance component estimation to GPS observables. J Surv Eng 135(4):149–160. https://doi.org/10.1061/(ASCE)0733-9453(2009)135:4(149)
Amiri-Simkooei AR, Zangeneh-Nejad F, Asgari J (2013) Least-squares variance component estimation applied to GPS geometry-based observation model. J Surv Eng 139(4):176–187. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000107
Bähr H, Altamimi Z, Heck B (2007) Variance component estimation for combination of terrestrial reference frames. Tech. rep., Karlsruher Institut für Technologie (KIT), https://doi.org/10.5445/KSP/1000007363
Beutler G, Kouba J, Springer T (1995) Combining the orbits of the IGS analysis centers. Bull Géod 69(4):200–222. https://doi.org/10.1007/BF00806733
Bisnath S, Gao Y (2009) Current state of precise point positioning and future prospects and limitations. Observing our changing earth. Springer, New York, pp 615–623. https://doi.org/10.1007/PL00012883
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, Berlin, Heidelberg. pp 213–221. https://doi.org/10.1007/978-3-642-29175-3_19
CSNO (2018) The BDS-3 preliminary system is completed to provide global services. http://en.beidou.gov.cn/WHATSNEWS/201812/t20181227_16837.html
Dach R, Sušnik A, Grahsl A, Villiger A, Schaer S, Arnold D, Prange L, Jäggi A (2019) Improving GLONASS orbit quality by re-estimating satellite antenna offsets. Adv Space Res 63(12):3835–3847. https://doi.org/10.1016/j.asr.2019.02.031
Dach R, Schaer S, Arnold D, Kalarus M, Prange L, Stebler P, Villiger A, Jaeggi A (2020). CODE final product series for the IGS. https://doi.org/10.7892/boris.75876.4
Deng Z, Nischan T, Bradke M (2017) Multi-GNSS rapid orbit-. Clock- & EOP-product series GFZ data services. https://doi.org/10.5880/GFZ.1.1.2017.002
ESA (2016) Galileo begins serving the globe. http://www.esa.int/Our_Activities/Navigation/Galileo_begins_serving_the_globe
Fritsche M (2016) Multi-GNSS orbit and clock combination: preliminary results. In: EGU general assembly conference abstracts, Vienna, Austria, EPSC2016-7823
Gendt G, Dick G, Söhne W (1998) GFZ analysis center of IGS–annual report. IGS 1996 Ann Rep 169–181
Griffiths J (2019) Combined orbits and clocks from IGS second reprocessing. J Geodesy 93(2):177–195. https://doi.org/10.1007/s00190-018-1149-8
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 Geodesy 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. American Society for Quality Control, Milwaukee
Johnston G, Riddell A, Hausler G (2017) The international GNSS service. Springer handbook of global navigation satellite systems. Springer, New York, pp 967–982. https://doi.org/10.1007/978-3-319-42928-1_33
Katsigianni G, Loyer S, Perosanz F, Mercier F, Zajdel R, Sośnica K (2019) Improving Galileo orbit determination using zero-difference ambiguity fixing in a Multi-GNSS processing. Adv Space Res 63(9):2952–2963. https://doi.org/10.1016/j.asr.2018.08.035
Kobel C, Arnold D, Jäggi A (2019) Combination of precise orbit solutions for sentinel-3A using variance component estimation. Adv Geosci 50:27–37. https://doi.org/10.5194/adgeo-50-27-2019
Koch KR, Kusche J (2002) Regularization of geopotential determination from satellite data by variance components. J Geodesy 76(5):259–268. https://doi.org/10.1007/s00190-002-0245-x
Kouba J (2009) A Guide to using international GNSS service (IGS) products. Tech. rep, geodetic survey division natural resources canada ottawa
Kouba J, Springer T (2001) New IGS station and satellite clock combination. GPS Solut 4(4):31–36. https://doi.org/10.1007/PL00012863
Kouba J, Mireault Y, Lahaye F (1995) 1994 IGS orbit/clock combination and evaluation. In: IGS annual report 1994. Jet Propulsion Laboratory publication. pp 70–94
Kusche J (2003) Noise variance estimation and optimal weight determination for GOCE gravity recovery. Adv Geosci 1:81–85. https://doi.org/10.5194/adgeo-1-81-2003
Loyer S, Perosanz F, Versini L, Katsigianni G, Mercier F, Mezerette A (2018) CNES/CLS IGS analysis center: recent activities. In: Proceedings of IGS workshop 2018, Wuhan, China
Männel B, Brandt A, Sakic P, Brack A, Nischan T (2020a) Status of IGS reprocessing activities at GFZ. International association of geodesy symposia. Springer, Berlin
Männel B, Brandt A, Nischan T, Brack A, Sakic P, Bradke M (2020) GFZ final product series for the International GNSS Service (IGS). GFZ Data Services. https://doi.org/10.5880/GFZ.1.1.2020.002
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
MIT (2020) Release notes for gamit/globk 10.71. Tech. rep., Massachusetts Institute of Technology
Montenbruck O, Steigenberger P, Kirchner G (2013) GNSS satellite orbit validation using satellite laser ranging. In: Proceedings of ILRS workshop, Fujiyoshida, Japan. pp. 1–6
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
Pearlman MR, Noll CE, Pavlis EC, Lemoine FG, Combrink L, Degnan JJ, Kirchner G, Schreiber U (2019) The ILRS: approaching 20 years and planning for the future. J Geodesy 93(11):2161–2180. https://doi.org/10.1007/s00190-019-01241-1
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
Prange L, Villiger A, Sidorov D, Schaer S, Beutler G, Dach R, Jäggi A (2020) Overview of CODE’s MGEX solution with the focus on Galileo. Adv Space Res 66(12):2786–2798. https://doi.org/10.1016/j.asr.2020.04.038
Rao CR (1973) Linear statistical inference and its applications, 2nd edn. Wiley, New Jersey. https://doi.org/10.1002/9780470316436
Sakic P, Mansur G, Kitpracha C, Ballu V (2019) The geodeZYX toolbox: a versatile Python 3 toolbox for geodetic-oriented purposes. V. 4.0. GFZ data services https://doi.org/10.5880/GFZ.1.1.2019.002
Sakic P, Mansur G, Männel B (2020) A prototype for a Multi-GNSS orbit combination. In: Proceedings of the European navigation conference (ENC) 2020. pp 1–11. 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. Springer, Berlin, Heidelberg. pp. 1–9. https://doi.org/10.1007/1345_2022_158
Selmke I, Duan B, Hugentobler U (2018) Status of the TUM MGEX orbit and clock products. In: Proceedings of IGS workshop 2018, Wuhan, China
Sośnica K, Thaller D, Dach R, Steigenberger P, Beutler G, Arnold D, Jäggi A (2015) Satellite laser ranging to GPS and GLONASS. J Geodesy 89(7):725–743. https://doi.org/10.1007/s00190-015-0810-8
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 TA, 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. pp 242–249
Teunissen PJG (1988) Towards a least-squares framework for adjusting and testing of both functional and stochastic models, vol 26. Delft University of Technology, reprint 2004
Teunissen PJG, Amiri-Simkooei AR (2008) Least-squares variance component estimation. J Geodesy 82(2):65–82. https://doi.org/10.1007/s00190-007-0157-x
Turkington DA (2005) Matrix calculus and zero-one matrices: statistical and econometric applications. Cambridge University Press, Cambridge
Uhlemann M, Gendt G, Ramatschi M, Deng Z (2015) GFZ global multi-GNSS network and data processing results. International association of geodesy symposia, vol 143. Springer, Cham, pp 673–679. https://doi.org/10.1007/1345_2015_120
Villiger A, Dach R (2020) International GNSS service technical report 2019 (IGS Annual Report). Tech. rep., IGS Central Bureau and University of Bern. https://doi.org/10.7892/boris.144003
Watson GA (2006) Computing helmert transformations. J Comput Appl Math 197(2):387–394. https://doi.org/10.1016/j.cam.2005.06.047
Zumberge JF, Heflin MB, Jefferson DC, Watkins MM, Webb FH (1997) Precise point positioning for the efficient and robust analysis of GPS data from large networks. J Geophys Res Solid Earth 102(B3):5005–5017. https://doi.org/10.1029/96JB03860
Acknowledgements
We would like to thank the International Laser Ranging Service (ILRS), the International GNSS Service (IGS), and in particularly the MGEX Analysis Centers for providing their products.
Funding
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).
Author information
Authors and Affiliations
Contributions
GM, PS, AB, and BM designed the study. AB and GM developed the mathematical framework. GM and PS implemented the algorithm. GM, PS, AB, and BM wrote the paper. All the authors joined the research discussion and gave their feedback on the manuscript.
Corresponding author
Rights and permissions
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/.
About this article
Cite this article
Mansur, G., Sakic, P., Brack, A. et al. Combination of GNSS orbits using least-squares variance component estimation. J Geod 96, 92 (2022). https://doi.org/10.1007/s00190-022-01685-y
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00190-022-01685-y