Introduction

In October 2017, the European GNSS Agency (GSA) released a comprehensive set of satellite metadata for the Galileo FOC satellites. The available data set includes spacecraft properties, optical surface characteristics, the attitude law, and the phase center positions of the transmitter antenna with respect to the satellite centers of mass (PCOs) as well as azimuth- and nadir-dependent phase variations (PVs). Together with similar information released for the IOV satellites in December 2016, for the first time, this information became available for a whole GNSS. Since then, several studies have discussed resulting improvements in the geodetic analysis (Bury et al. 2019; Katsigianni et al. 2019; Zajdel et al. 2019; Li et al. 2019). In particular, the information about the PCOs and the PVs allows improved processing and new investigations, which are discussed within this study.

In order to realize the Terrestrial Reference System (TRS), which is the basis for all geodetic measurements on the earth, the geodetic datum has to be defined, i.e., origin, orientation, and scale have to be specified. Theoretically, GNSS can provide a terrestrial scale thanks to (1) centimeter-level accurate satellite orbits (Männel 2016) and (2) the precision of the GNSS phase measurements (observation error less than 2 mm). However, to link both orbit and observation, information is required about the transmitting point (reference for the observation) with respect to the satellite center of mass (reference for the orbit). Obviously, an unconsidered offset in the radial direction (i.e., in the z-direction of the spacecraft body-fixed frame) will shift the determined station heights and bias the eventually estimated terrestrial scale parameter. Unfortunately, the position of the transmitting point is usually not disclosed by the GNSS providers. For some recently launched satellites, ground calibrated PCOs are now provided, e.g., for Galileo, BeiDou-3, QZSS, and GPS III. For most currently and formerly operational satellites, however, the PCOs and the PVs have to be determined in global adjustments (Schmid et al. 2007). Over the past years, several PCO and PVs sets have been estimated for the different constellations, for example, by Steigenberger et al. (2016) and Schmid et al. (2016). Due to the high correlation between station height, troposphere delay, and the offsets of transmitting and tracking antennas, accurate calibrations of the tracking antennas are a prerequisite for estimating the transmitting antenna offsets. The corresponding robot calibrations are provided in the International GNSS Service (IGS) antenna exchange format (ANTEX). Moreover, thanks to a recent effort by Geo +  + , signal-specific (including Galileo frequencies) and multi-GNSS calibrations are available for many receiver antennas used within the IGS tracking network, for example, in the ANTEX file for IGS repro3 igsR3_2057.atx provided by Villiger (2019). In addition, the terrestrial scale had to be fixed, for example, to the current ITRF solution, to avoid a poorly conditioned normal equation system with less precise estimates. Consequently, the derived transmitter offsets and any further derived geodetic products are not independent of this ITRF scale anymore (Haines et al. 2015).

By fixing the transmitter antenna patterns of Galileo satellites to the ground calibrated values, a GNSS-based terrestrial scale becomes achievable. However, with the first operational Galileo satellites launched in 2012, a corresponding Galileo-only solution could cover only the most recent years (i.e., from 2017 onward). To process a long-time solution and determine the terrestrial scale back in time, the PCOs, which are independent of the terrestrial scale and derived by other techniques such as very long baseline interferometry (VLBI) and satellite laser ranging (SLR), are still required for GPS and GLONASS. We present two different approaches to re-adjust these offsets. It will cross-check and compare both approaches and present the resulting terrestrial scale values for an exemplary period (first half of 2019).

To improve readability, we will use the following naming convention. PCOs describe the offset between the center of mass and mean transmitting point onboard the spacecraft and the offset between the antenna reference point and mean transmitting center for receiving antennas. Deviation from the mean transmitting or receiving point is described by PVs, which are nadir and azimuth or elevation and azimuth dependent, respectively. Transmitter phase centers are identified by the satellite system in a superscript (e.g., \({\text{PCO}}^{{{\text{GPS}}}}\)). Receiving antennas are indicated by subscripts (e.g., \({\text{PCO}}_{{{\text{LEO}}}}\)). The estimated PCO differences in the z-direction with respect to the a priori values are indicated by z-\(\Delta {\text{PCO}}\), e.g., z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\). The manuscript is structured as follows. Following this introduction, the two PCO determination approaches are introduced in more detail and the validation and comparison scheme is explained. Subsequently, the processing period, the selection of ground stations, the quality check of the Swarm orbits, and the details about the processing strategy are introduced. The results are presented and discussed in the section afterward. The summary and our conclusion are given in the last section.

Methods for phase center offset estimation

This section describes two methods used to derive \({\text{PCO}}^{{{\text{GPS}}}}\) in the z-direction (z-\({\text{PCO}}^{{{\text{GPS}}}}\)) without fixing the terrestrial scale. Also, the validation procedure used later and the definition of the assessed cases are described. Both approaches rely on additional observations, either ground Galileo observations or GPS observations onboard low earth orbiters (LEOs). Figure 1 presents the basic setup consisting of ground stations, GPS and Galileo satellites, and LEOs. The satellite orbits are constrained by celestial mechanics. Ground-based and space-based observations connect the antenna phase center of different transmitters and receivers. The estimated coordinates of the ground station network have a scale factor with respect to the a priori terrestrial frame, in our case IGS14 (Rebischung and Schmid 2016).

Fig. 1
figure 1

Schematic diagram of the two methods of determining scale-independent GPS z-PCO

Method I: Galileo with calibrated antenna offsets

The basic idea of this method is to separate z-\({\text{PCO}}^{{{\text{GPS}}}}\) and terrestrial scale by adding Galileo (GAL) observations. With the \({\text{PCO}}^{{{\text{GAL}}}}\) fixed to the calibrated values provided by the GSA, a reliable scale-independent network solution is achieved. As GPS and Galileo are observed by the same stations whose coordinates are now estimated scale independently from the underlying reference frame, also the z-\({\text{PCO}}^{{{\text{GPS}}}}\) can be estimated scale independently. This method will fail if there is any systematic bias between independently estimated station coordinates for GPS and Galileo. Villiger et al. (2018, 2019) reported translational biases of several mm when applying the L1 and L2 PVs of GPS to the Galileo E1 and E5 signals. With the signal-specific antenna corrections provided by Geo +  +, this systematic discrepancy should not occur anymore. This assumption was tested by processing GPS and Galileo solutions independently in the framework of the next IGS reprocessing campaign (repro3). However, due to the different satellite PCOs used (z-\({\text{PCO}}^{{{\text{GPS}}}}\) from igs14.atx and z-\({\text{PCO}}^{{{\text{GAL}}}}\) from GSA) a terrestrial scale bias of 1.16 \(\pm\) 0.27 ppb (part per billion) was observed in the GFZ submission (Männel et al. 2020). When taking this terrestrial scale into account, GPS and Galileo-based coordinates agree on the level of a few millimeters in the height component. By fixing the antenna offsets of Galileo to the calibrated values, z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) have been computed by Villiger et al. (2020). He reported a system-wise change of \(- 150\) mm and \(- 221\) mm for the z-\({\text{PCO}}^{{{\text{GPS}}}}\) by using robot calibrations and chamber calibrations for ground stations, respectively. The re-adjusted PCOs have been updated in the IGS repro3 ANTEX file (igsR3_2077.atx) and will be used in the IGS repro3 processing.

Method II: gravitational constraint

The orbits of the LEOs are scale independent as their radial position is constrained by orbital dynamics (so-called gravitational constraint). Therefore, the estimation of scale-independent z-\({\text{PCO}}^{{{\text{GPS}}}}\) becomes possible. However, there are three limitations. Firstly, there are not enough space-based observations to solve for all z-\({\text{PCO}}^{{{\text{GPS}}}}\), LEO orbits, GPS orbits, etc. Therefore, ground- and space-based observations have to be combined. This approach is known as an integrated or one-step approach and has been studied for the past 15 years. It was already used to determine z-\({\text{PCO}}^{{{\text{GPS}}}}\) by Haines et al. (2015) and Männel (2016). To transfer the scale constraint offered by the space-based observations requires a fully consistent estimation of GNSS satellite orbits and clocks which link ground- and space-based observations. The second limitation is the availability and quality of the space-based observations. And thirdly, an error in the a priori calibrated z-\({\text{PCO}}_{{{\text{LEO}}}}\) can significantly bias the derived z-\({\text{PCO}}^{{{\text{GPS}}}}\).

Validation

In general, the validation of z-\({\text{PCO}}^{{{\text{GPS}}}}\) is challenging as the phase center offsets cannot be observed by space geodetic techniques. However, we can evaluate the different z-\({\text{PCO}}^{{{\text{GPS}}}}\) estimated by both methods by comparison and cross-check. First of all, scale independence can be mathematically assessed by comparing the correlation between scale and phase center parameters. Using only ground-based GPS observations results in a large correlation coefficient of the two parameters (Schmid et al. 2007). Using both methods with different observations and constraints (on \({\text{PCO}}^{{{\text{GAL}}}}\) or \({\text{PCO}}_{{{\text{LEO}}}}\)) allows, secondly, to assess the agreement between the z-\({\text{PCO}}^{{{\text{GPS}}}}\) estimates. For this purpose, we designed six cases that are listed in Table 1. Different combinations of included satellites, PCO constraints, and estimated satellite z-PCO increments with respect to the a priori values (\(\Delta {\text{PCO}}\)) are selected for the different cases. Since our focus is on the satellite z-PCO and the terrestrial scale, the satellite PCOs in x- and y-directions are always kept fixed. In case 1, as a reference case, we want to show the problem of a high correlation between the terrestrial scale and the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\). In cases 2 and 3, z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) are estimated by fixing either \({\text{PCO}}_{{{\text{LEO}}}}\) or \({\text{PCO}}^{{{\text{GAL}}}}\). Moreover, we combine GPS and Galileo ground-based observations and space-based GPS observations in cases 4, 5, and 6. In case 4, both \({\text{PCO}}^{{{\text{GAL}}}}\) and \({\text{PCO}}_{{{\text{LEO}}}}\) are fixed to estimate the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\). In cases 5 and 6, we determine z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) or z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) jointly with z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) while only fixing \({\text{PCO}}^{{{\text{GAL}}}}\) or \({\text{PCO}}_{{{\text{LEO}}}}\), respectively. These two cases allow the ultimate cross-check with the known Galileo and LEO offsets. However, it is debatable whether the gravitational constraint can be transferred from the GPS space-based observations to the Galileo satellites or reversely. (Unfortunately, space-based observations are available only for GPS.) This question will be discussed in the section of the results. To improve the readability, we name the six cases based on the included satellites and the estimated z-\(\Delta {\text{PCO}}\). For example, GEL-GE means that GPS (G), Galileo (E), and Swarm (L) satellites are all included in the processing and z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\)(-G) and z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\)(-E) are estimated, while only the PCOs of Swarm satellites are fixed.

Table 1 Six cases for the estimation of the GPS z-PCO and GNSS-based terrestrial scale deriving

Processing period and ground station selection

We selected day 1 to 180 of 2019 as our processing period. During this period, the Galileo constellation already had 24 satellites in operation. All selected ground stations are tracking both GPS and Galileo, and the network is globally and evenly distributed. As a prerequisite for the terrestrial scale realization, the stations should have accurate coordinates that are offered within the IGS products (i.e., in IGS14 reference frame). There are 68 to 94 stations (only a few days less than 75) that are selected for different days. The majority of the stations for each daily processing have Galileo-calibrated receiver antennas (Fig. 2), and for the others, the GPS L2 calibrations are applied for E5a. We used only stations that provide observations in RINEX3 format to the IGS data centers. The station number increases around DOY (day of the year) 87 because more stations started to offer RINEX3 observations from that day onward. Figure 3 presents the selected 75 stations for the processing of DOY 1 as an example.

Fig. 2
figure 2

Number of stations selected for each day

Fig. 3
figure 3

Distribution of the 75 stations selected for January 1, 2019. Stations with Galileo antenna calibrations are marked in blue. Red denotes the stations using phase center corrections of GPS for the Galileo signals

Swarm orbit quality

For the gravitational constraint strategy, we included the three spacecraft of the Swarm mission, which is a mini-satellite constellation mission to survey the geomagnetic field (Friis-Christensen et al. 2008). The three satellites (Swarm-A, B, and C) are operated in two different orbit configurations. Swarm-A and Swarm-C are flying at a mean altitude of 450 km and in an 87.4° inclined orbital plane, while Swarm-B has a higher inclination of 88° and a larger mean altitude of 530 km. To check the quality of the LEO observation data and to verify our orbit determination, we did a Swarm-only reduced-dynamic POD by using IGS final orbit and 30-s clock products. The data sampling rate is 30-s and the arc length is 24-h. The determined orbits are validated by comparing with an external solution and by SLR observation residual validation. The daily orbit is compared with the official precise orbit products which are offered by the European Space Agency (ESA, Olsen 2019) in the along-track, cross-track, and radial directions. The orbit RMS values averaged over 180 days are presented in Table 2. In general, the orbits of the three Swarm satellites are determined in similar accuracy with about 30-mm RMS in the along-track direction, 14-mm RMS in the cross-track direction, and 24-mm RMS in the radial direction. We used all SLR observations of the high-quality Yarragadee station in Australia during the 180 days to validate the Swarm orbits. The statistical results are also listed in Table 4, and all the epoch-wise solutions are shown in Fig. 4. With the most observations, the SLR residuals of Swarm-B have the largest mean (4.2 mm) and the smallest variation (\(\pm\) 19.5 mm). With similar numbers of observations, the SLR residuals (with variation) of Swarm-A and Swarm-C are 3.7 \(\pm\) 25.1 mm and 2.4 \(\pm\) 25.0 mm. Compared with previous studies, the orbit quality of our solution is at a comparable level.

Table 2 The validation of the orbits of Swarm satellites
Fig. 4
figure 4

SLR observation residuals for the Swarm-A, B, and C POD solution for all passes of the Yarragadee station in Australia for 180 days. The gaps are caused by missing SLR observations

Processing strategy

We use the software PANDA (Liu and Ge 2003) for the processing. We performed the one-step method (Montenbruck and Gill 2000) to estimate the orbits of GPS, Galileo, and Swarm satellites, the PCOs in the z-direction of different satellites, and the other parameters simultaneously. Detailed information on the orbit modeling, processing configuration, metadata, and estimated parameters is listed in Table 3. As the initial release of the antenna correction file for IGS Repro3, igsR3_2057.atx includes IGS estimated \({\text{PCO}}^{{{\text{GPS}}}}\) (Schmid et al. 2016), the GSA calibrated \({\text{PCO}}^{{{\text{GAL}}}}\) and the ground calibrated \({\text{PCO}}_{{{\text{station}}}}\) for multi-GNSS. Depending on the cases in Table 1, the satellite PCOs in the z-direction are fixed to their a priori values or estimated freely. To derive scale-independent z-\({\text{PCO}}^{{{\text{GPS}}}}\) and a GNSS-based terrestrial scale, the coordinates of the ground stations are constrained only by applying no-net-rotation condition. The scale of the ground-tracking network is not constrained in this study, and the scale is derived by Helmert transformation between the estimated coordinates of the ground network and the a priori coordinates (i.e., IGS14). Montenbruck et al. (2018) reported in-flight calibrated PVs for the Swarm satellites of up to 25 mm. As these in-flight calibration results might not be independent of the scale provided by VLBI and SLR, we decided not to apply them in this study.

Table 3 Processing details (orbit modeling, processing configurations, metadata, and estimated parameters)

Results

We analyze the results from three aspects. Firstly, considering the relationship that 130-mm error in GPS z-\({\text{PCO}}^{{{\text{GPS}}}}\) leads to one ppb terrestrial scale (Zhu et al. 2003), we discuss both the estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the derived terrestrial scale with respect to IGS14. The further comparisons and the estimation quality analysis are based on the daily estimates, the formal error of the estimates, and the correlation coefficient of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and scale. The variation of the estimated daily z-\(\Delta {\text{PCO}}\) values, the formal error of z-\(\Delta {\text{PCO}}\), and the derived scale between the processed days are shown by the empirical standard deviation (STD) of their time series with respect to the mean. Both satellite-specific results and the results averaged over satellites (system-wise) are discussed in detail. Secondly, the z-\(\Delta {\text{PCO}}\) estimated by fixing only the \({\text{PCO}}^{{{\text{GAL}}}}\) in GEL-GL will be analyzed. The impact of fixed \({\text{PCO}}^{{{\text{GAL}}}}\) on the estimation of z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) is shown. At last, mainly based on GEL_GE, the z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) estimated by fixing only the \({\text{PCO}}_{{{\text{LEO}}}}\) is analyzed. The effect of transferring the gravitational constraint directly to GPS and indirectly to Galileo via GPS satellites and ground stations is discussed.

Estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and terrestrial scale

In Fig. 5, the satellite-specific z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) with respect to the IGS values and averaged over the 180 processed days are shown as blue bars. The vertical lines denote the empirical STD for each time series. The last bar in each plot provides the mean value over all satellites; correspondingly, the empirical STD of the constellation-wise value is smaller than that of the satellite-specific values. The formal errors of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and their empirical STD are presented as green bars. Due to the evenly distributed ground network and satellite constellation, the formal errors are quite similar within one case. There is no obvious block-specific phenomenon visible. Although the z-PCOs of GPS satellites in the same block are similar, the z-PCOs corrections are similar for all satellites in every case.

Fig. 5
figure 5

Estimated GPS z-PCO differences with respect to IGS values (blue) and their formal errors (green). The name of each case is based on the included satellites and the estimated z-\(\Delta {\text{PCO}}\), for example, GEL-GE means that GPS, Galileo, and Swarm satellites are included and z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) are estimated. Each bar denotes the solution averaged over 180 processing days. The vertical lines denote the empirical standard deviation of the time series with respect to the mean. The x-label presents the space vehicle number of the satellites, and the satellites are sorted by block-wise

In the case G-G, the estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) values are smaller than 100 mm, but with large empirical STD (100 to 130 mm), formal error (about 46 mm), and empirical STD of formal error (about 22 mm) among all cases. The reason for this is the high correlation between the estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the terrestrial scale. Slight changes in any inputs of the estimation (e.g., the ground station network) lead to very different solutions; therefore, the precision of the estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the scale is low.

In the other five cases, the PCOs of either the Galileo or the Swarm satellites or both are fixed. Consequently, the precision of the z-PCO estimates is improved. In general, the results of the five cases show collective shifts of z-\({\text{PCO}}^{{{\text{GPS}}}}\) with respect to the IGS values, and the satellite-specific values have a good agreement among the five cases. Comparing the results based on the gravitational constraint (GL-G) and on Galileo (GE-G), the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) values have differences of about 30 mm for all satellites. The empirical STD of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\), the formal error of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\), and the empirical STD of the formal error of GL-G are 12 mm, 5 mm, and 3 mm larger than those of GE-G, respectively. That means the precision of the LEO-PCO-fixed case is slightly lower than that of the Galileo-PCO-fixed cases. It is explained by the stronger constraint transferred by many more observations from 24 Galileo satellites compared to the three Swarm satellites, which is verified later in this section.

In the GEL-G case, the PCOs of both Galileo satellites and Swarm satellites are fixed, but the results are quite similar to GE-G. Similar results are obtained in GEL-GL in which only Galileo PCOs are fixed. This demonstrates that the Galileo satellites are dominating the results due to the larger number of observations. In GL-G and GEL-GE, only the PCOs of Swarm satellites are fixed. However, the result differences between GL-G and GEL-GE are larger than the result differences between GE-G, GEL-G, and GEL-GL. The z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) values in GEL-GE are collectively larger than that of GL-G by about 10 mm. The empirical STD of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\), the formal error of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\), and the empirical STD of formal error are all smaller in GEL-GE than in GL-G. The differences between GL-G and GEL-GE are caused by including Galileo.

The time series of daily system-wise (averaged over satellites) z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the corresponding terrestrial scale are shown in Fig. 6 for G-G and in Fig. 7 for the other five cases. The corresponding mean values and the standard deviations of all the time series are presented in Table 4. Comparing the upper (z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\)) and the lower (scale factor) plots, we can see the relationship between the two parameters. The variation of the time series in G-G is quite large (103 mm empirical STD for z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and 0.823 ppb empirical STD for terrestrial scale). The solutions of the Galileo-PCO-fixed solutions (GE-G, GEL-G, and GEL-GL) are very similar. The time series of GL-G and GEL-GE have larger variation and -20 to -40 mm differences in mean values of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) than those of the Galileo-PCO-fixed solutions. By including Galileo satellites, GEL-GE is more stable and closer to the Galileo-PCO-fixed solutions than GL-G.

Fig. 6
figure 6

Time series of the differences between the estimated GPS z-PCO and IGS values averaged over satellites (upper) and the corresponding terrestrial scale with respect to IGS14 (lower) in the case including GPS only without fixing the scale

Fig. 7
figure 7

Time series of the estimated GPS z-PCO differences with respect to IGS values averaged over satellites (upper) and the corresponding terrestrial scale with respect to IGS14 (lower) for the five cases. The name of each case is based on the included satellites and the estimated z-\(\Delta {\text{PCO}}\), for example, GEL-GE means that GPS, Galileo, and Swarm satellites are included and z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) are estimated

Table 4 The estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) averaged over satellites and processed days and the scale factor with respect to IGS14 averaged over the processed days. The empirical standard deviations of the time series (STD) are also given. The correlation coefficients of the estimated satellite PCOs in the z-direction and the terrestrial scale. The values are averaged over satellites and processed days. Zero values are due to the fixing to the a priori values. The dash means not included

The impact of the 20 additional stations after DOY 89 (Fig. 2) on the estimates is not visible. Only a slight decrease of the formal errors is observed in the analysis.

The quality of the estimation in the different cases is also reflected in the correlation coefficients of the estimated z-\(\Delta {\text{PCO}}\) and the terrestrial scale. The corresponding correlation coefficients averaged over satellites and days are presented in Table 4. Overall, the coefficients are very stable with variations smaller than 0.01. G-G shows the largest correlation coefficient of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and terrestrial scale (0.87) which agrees to the analysis mentioned above. The correlation coefficient can be reduced effectively by introducing LEOs or by processing together with Galileo and fixing \({\text{PCO}}^{{{\text{GAL}}}}\). Derived by different numbers of observations, the Galileo-PCO-fixed case GE-G is more effective than the LEOs-PCO-fixed case GL-G in de-correlation (reduction of 0.74 versus 0.35). Due to the stronger impact of Galileo on transferring the constraint compared to Swarm, the correlation coefficient nearly does not change after fixing \({\text{PCO}}_{{{\text{LEO}}}}\) additionally (GEL-G and GEL-GL). The correlation coefficients in GEL-GL show that the fixed Galileo PCOs can separate the derived terrestrial scale and the estimated z-\(\Delta {\text{PCO}}\) for both GPS and LEOs. In GEL-GE, with only three PCO-fixed LEOs, the correlation between z-\(\Delta {\text{PCO}}\) and terrestrial scale is identical for GPS and Galileo satellites (0.56) and is close to that of GL-G (0.52).

To investigate the impact of the numbers of Galileo and Swarm satellites on the estimation, GE-G is processed again by only including three Galileo satellites (E101, E210, and E212) in three different orbital planes (GE-G*). The statistic of the solution for GE-G* is presented in Table 4. With fewer Galileo satellites, the results of GE-G* are different from those of GE-G with a system-wise difference of 25 mm for the estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\). This is caused by the weaker geometry and fewer observations of the three Galileo satellites compared to the full system. Without the advantage of more satellites, the precision of GE-G* becomes lower than that of the GL-G with 13-mm larger empirical STD of the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and 0.1 ppb larger empirical STD of the scale. Moreover, the correlation coefficient between z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the terrestrial scale increases from 0.13 (GE-G) to 0.54 (GE-G*), which exceeds that of GL-G by 0.02. In a summary, due to the much faster geometry change, including three Swarm satellites gives more precise z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) than including three Galileo satellites.

Besides the internal comparison and cross-check between the different cases, we also compared our results with other studies. The system-wise z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) derived by GE-G is between the robot-calibration-based solution and the chamber-calibration-based solution in Villiger et al. (2020). We compared the estimated GNSS-based scale with the scale determined by the VLBI and SLR. As reported by Altamimi et al. (2016), the scale factors determined by VLBI and SLR with respect to the ITRF 2014 are about + 0.77 ppb and -0.77 ppb, respectively. The GNSS-based scales derived by GL-G and GE-G cases are + 1.72 ppb and + 1.55 ppb, respectively. Therefore, the GNSS-based scale derived by both methods agree with each other well and agree with the VLBI-based scale better than the SLR-based scale does. After removing the systematic errors in SLR data by Luceri et al. (2019), the scale derived by SLR is about + 1 ppb toward ITRF 2014 scale. Therefore, the scales determined by GNSS in this study, by VLBI, and by SLR have an agreement within differences smaller than 1 ppb.

Estimated z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\)

In the case GEL-GL, z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) are estimated simultaneously by fixing the PCOs of all Galileo satellites to the GSA values. Figure 8 shows the time series of the estimated satellite-specific z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\). The mean values and the empirical STD are \(- 2.2 \pm 2.5\) mm, \(- 2.6 \pm 2.1\) mm, and \(- 1.1 \pm 2.4\) mm for Swarm-A, B, and C satellites, respectively. The plots of Swarm-A and Swarm-C are very similar, but that of Swarm-B is slightly different from them. This can be explained by their orbital configuration introduced in Section “Swarm orbit quality.” During DOY 55 to 57, the orbits of Swarm-B have a 10-mm larger RMS with respect to the official products than the other days, which might be caused by some unknown behavior of the spacecraft. This is assumed to cause the large deviation of the estimated z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) on those days. It also affects all the time series of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\), z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\), and the scale derived by the LEOs-PCO-fixed cases. Therefore, we concluded that orbit modeling quality has a large impact on the estimation.

Fig. 8
figure 8

Time series of the estimated Swarm-A, B, and C z-PCO differences with respect to the a priori values (ESA offered) in the case including GPS, Galileo, and Swarm satellites and fixing only the PCO of Galileo satellites

All of the three time series show a periodic behavior. The periodicity might be related to the draconitic period, i.e., the time between two passages of the satellite through its ascending node, of Swarm, Galileo, and GPS. The impact of the periodicity can also be observed from the time series variation of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) estimated in GL-G and GEL-GE in which only the PCOs of LEOs are fixed (Figs. 7 and 12). From the magnitude of the estimated z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\), the importance of the PCO accuracy of the LEOs can be realized. The scale factor between GL-G and GEL-GL has a 0.2 ppb (1.3 mm at the equator) difference, while the estimated z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) in GEL-GL is 1–2 mm with respect to the a priori values which are fixed in GL-G. Additionally, we processed an update for GL-G (GL-G*) using artificially modified Swarm PCOs by adding the estimated z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) to the values offered by ESA. The time series of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) estimated in GL-G* is presented in Fig. 9 (green). The curve of the updated case is systematically shifted from the curve of GL-G by 48 mm. However, the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) averaged over satellite and processed days is \(- 173\) mm which is much closer to the Galileo-PCO-fixed solution in GEL-GL (red) than that of GL-G (purple). This comparison shows the importance of the accuracy of the LEO PCOs again.

Fig. 9
figure 9

Time series of the estimated GPS z-PCO differences with respect to IGS values averaged over satellites in three cases. GL-G includes GPS and Swarm satellites and the PCOs of Swarm satellites are fixed. GEL-GL includes GPS, Galileo, and Swarm satellites, and only the PCOs of Galileo satellites are fixed. GL-G* is an update processing of GL-G by modifying the z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) artificially

Estimated z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\)

In the case GEL-GE, z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) and z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) are estimated simultaneously by fixing the PCOs of the three LEOs to a priori values and without constraint on the terrestrial scale. This allows us to discuss the estimated z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) with respect to the GSA values. Since the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the terrestrial scale derived in GEL-GE have small differences with respect to the solutions derived by the Galileo-PCO-fixed cases (GE-G, GEL-G, and GEL-GL), the estimated z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) in GEL-GE are small. The satellite-specific z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) values are presented as blue bars in Fig. 10. The z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) averaged over satellites is only \(- 21\) mm. The empirical STD values of the z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) time series are about 80 mm to 100 mm for different Galileo satellites, and they are about 30 mm to 50 mm larger than that of the GPS satellites shown in Fig. 5. The reason is that the gravitational constraint on LEOs is transferred only via the GPS satellites and the ground stations to Galileo (Fig. 1). Since the selected ground network changes day by day during the processed period, the impact on individual Galileo satellites will change. However, the constraints on the LEO orbits affect z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) directly by onboard GPS observations; therefore, the variations of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) are smaller. Evaluating the impact on the whole Galileo constellation, the empirical STD of z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) averaged over all satellites is only 10 mm larger than that of GPS. The formal errors of z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) and the empirical STD of the formal errors are only about 6 mm and 2 mm larger than those of the estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\). In general, the gravitational constraint on the LEOs acts on the estimation of z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\). We also see some unexpected phenomena on satellite E102. We expected a systematic change of the estimated z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) due to the scale change, but the absolute value of the estimated z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) of E102 is much larger than all the other satellites. For example, the z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) of E101 and E102 has a \(- 123\) mm difference. However, the a priori (GSA) z-PCOs of the two satellites have a \(+ 87\) mm difference. That means our estimated z-PCO values of the two satellites are much closer than their a priori values. This result agrees with the study by Steigenberger et al. (2016).

Fig. 10
figure 10

The estimated Galileo z-PCO differences with respect to igsR3_2057.atx (upper) and their formal errors (lower). Each bar denotes the solution averaged over 180 processing days. The thin errors bars denote the empirical standard deviation of the time series

For the four Galileo satellites E219, E220, E221, and E222, which were launched in July 2018, we found larger formal errors than for the other Galileo satellites. This is likely caused by fewer ground-based observations that were available during the first 43 processed days, as a part of the ground stations was not offering data for them. The numbers of observations are shown in Fig. 11. Moreover, the observations of satellites E222 and E219 reach a similar number as the other satellites in mid-2019. This is due to the limited capability of some ground receivers which only observe satellites with PRN smaller than 32 (Mozo 2018).

Fig. 11
figure 11

Number of ground-based observations that are used for the processing of Galileo satellites E102 (as a reference), E219, E220, E221, and E222 during the 180 processed days

We also plot the time series of system-wise z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) and z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) in Fig. 12. As explained above, the variation of the z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) time series is slightly larger (10 mm) than that of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\). However, due to the fixed LEO PCOs and the same ground stations, they agree with each other.

Fig. 12
figure 12

Time series of the estimated GPS and Galileo z-PCO differences with respect to igsR3_2057.atx averaged over satellites in the case including GPS, Galileo, and Swarm satellites and only fixing the PCOs of Swarm satellites

Summary and Conclusions

Using two different methods based on (1) the Galileo system with ground-calibrated antenna offsets and on (2) the gravitational constraint on LEO orbits, we determined the scale-independent GPS z-PCO and the corresponding GNSS-based terrestrial scale. Applying the first method, we found a \(- 186 \pm 25\) mm z-PCO correction with respect to the IGS values, and a \(+ 1.55 \pm 0.22\) ppb terrestrial scale with respect to the IGS14. The results of the gravitational constraint method are \(- 221 \pm 37\) mm for the z-PCO and \(+ 1.72 \pm 0.31\) ppb for the terrestrial scale. The solutions derived by the two independent methods with different observations and metadata agree well with each other. The Galileo-based solution agrees very well with the latest study by Villiger et al. (2020). Moreover, these two solutions also agree with the VLBI-based scale (+ 0.77 ppb) better than the SLR-based scale (-0.77) does. Compared with the updated SLR-based scale without systematic errors (Luceri et al. 2019), the scales determined by GNSS in this study, by VLBI, and by SLR agree with each other with differences smaller than 1 ppb.

Since Galileo offers many more observations which transfer the constraints than the Swarm constellation, Galileo dominated the results of the case in which the PCOs of both Galileo and LEOs are fixed. Based on the correlation coefficient of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and scale, the formal error of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the empirical STD of the time series, the precision and stability of the solution derived by the Galileo-PCO-fixed method is higher than that derived by the LEO-PCO-fixed method. This is mainly caused by the different number of satellites and observations from Galileo and Swarm. If Galileo is reduced from the full constellation to only three satellites, the better geometry of the Swarm-based solution leads to better results.

The joint estimation of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) by only fixing \({\text{PCO}}^{{{\text{GAL}}}}\) showed that the z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and the derived scale factor are very close to the solutions derived by the case including GPS and Galileo and fixing the \({\text{PCO}}^{{{\text{GAL}}}}\). Consequently, the constraint from Galileo is very strong and is nearly unaffected by including LEOs. The z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) precisely estimated at 1 to 2 mm with respect to the values offered by ESA. This shows the small difference between the two methods again. Moreover, the accuracy of the LEOs’ PCOs is very important for the gravitational constraint method. We realized some periodic variations in the z-\(\Delta {\text{PCO}}_{{{\text{LEO}}}}\) time series. This is also visible in the time series of z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\), z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) and scale derived by applying only the gravitational constraint and might be related to the draconitic period of GPS and Swarm constellations. Based on the unusual results of the Swarm-B satellite in three days, the importance of orbit modeling quality is shown.

The z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) estimated by only fixing \({\text{PCO}}_{{{\text{LEO}}}}\) in the GPS, Galileo, and Swarm joint processing differs on average by \(- 21\) mm from the GSA values. This difference corresponds to 0.13 ppb difference in the terrestrial scale. The estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\) and scale have slight differences from the results derived by the case, which includes only GPS and LEOs. The precision and stability of z-\(\Delta {\text{PCO}}^{{{\text{GAL}}}}\) are both worse than those of the simultaneously estimated z-\(\Delta {\text{PCO}}^{{{\text{GPS}}}}\). The gravitational constraint on the Swarm orbits is partially transferred to the Galileo satellites. This situation can be improved by including more LEOs and moreover, by including Galileo space-based observations, which might be available in the future.

A future study including a longer processing period (at least three years) and more LEOs will be performed to study the GNSS-based terrestrial scale in detail.