Two methods to determine scale-independent GPS PCOs and GNSS-based terrestrial scale: comparison and cross-check

The GPS satellite transmitter antenna phase center offsets (PCOs) can be estimated in a global adjustment by constraining the ground station coordinates to the current International Terrestrial Reference Frame (ITRF). Therefore, the derived PCO values rest on the terrestrial scale parameter of the frame. Consequently, the PCO values transfer this scale to any subsequent GNSS solution. A method to derive scale-independent PCOs without introducing the terrestrial scale of the frame is the prerequisite to derive an independent GNSS scale factor that can contribute to the datum definition of the next ITRF realization. By fixing the Galileo satellite transmitter antenna PCOs to the ground calibrated values from the released metadata, the GPS satellite PCOs in the z-direction (z-PCO) and a GNSS-based terrestrial scale parameter can be determined in GPS + Galileo processing. An alternative method is based on the gravitational constraint on low earth orbiters (LEOs) in the integrated processing of GPS and LEOs. We determine the GPS z-PCO and the GNSS-based scale using both methods by including the current constellation of Galileo and the three LEOs of the Swarm mission. For the first time, direct comparison and cross-check of the two methods are performed. They provide mean GPS z-PCO corrections of -186±25\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- 186 \pm 25$$\end{document} mm and -221±37\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- 221 \pm 37$$\end{document} mm with respect to the IGS values and +1.55±0.22\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+ 1.55 \pm 0.22$$\end{document} ppb (parts per billion) and +1.72±0.31\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+ 1.72 \pm 0.31$$\end{document} in the terrestrial scale with respect to the IGS14 reference frame. The results of both methods agree with each other with only small differences. Due to the larger number of Galileo observations, the Galileo-PCO-fixed method leads to more precise and stable results. In the joint processing of GPS + Galileo + Swarm in which both methods are applied, the constraint on Galileo dominates the results. We discuss and analyze how fixing either the Galileo transmitter antenna z-PCO or the Swarm receiver antenna z-PCO in the combined GPS + Galileo + Swarm processing propagates to the respective freely estimated z-PCO of Swarm and Galileo.


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 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., PCO GPS ). Receiving antennas are indicated by subscripts (e.g., PCO LEO ). The estimated PCO differences in the z-direction with respect to the a priori values are indicated by z-ΔPCO , e.g., z-ΔPCO GPS and z-ΔPCO 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 PCO GPS in the z-direction (z-PCO 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).

Method I: Galileo with calibrated antenna offsets
The basic idea of this method is to separate z-PCO GPS and terrestrial scale by adding Galileo (GAL) observations. With the PCO 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-PCO 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. (2018Villiger et al. ( , 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-PCO GPS from igs14.atx and z-PCO GAL from GSA) a terrestrial scale bias of 1.16 ± 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 Galileobased 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-ΔPCO GPS have been computed by Villiger et al. (2020). He reported a system-wise change of −150 mm and −221 mm for the z-PCO 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-PCO GPS becomes possible. However, there are three limitations. Firstly, there are not enough space-based observations to solve for all z-PCO 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-PCO 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-PCO LEO can significantly bias the derived z-PCO GPS .

Validation
In general, the validation of z-PCO GPS is challenging as the phase center offsets cannot be observed by space geodetic techniques. However, we can evaluate the different z-PCO 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 PCO GAL or PCO LEO ) allows, secondly, to assess the agreement between the z-PCO 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 ( Δ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-ΔPCO GPS . In cases 2 and 3, z-ΔPCO GPS are estimated by fixing either PCO LEO or PCO GAL . Moreover, we combine GPS and Table 1 Six cases for the estimation of the GPS z-PCO and GNSSbased terrestrial scale deriving The column named "Satellites" shows the satellites included in the processing. "Fixed" means that the corresponding satellite PCOs are fixed to the a priori values. The last column shows the estimated corrections of satellite z-PCO in the processing. The name of each case is based on the included satellites and the estimated z-ΔPCO Galileo ground-based observations and space-based GPS observations in cases 4, 5, and 6. In case 4, both PCO GAL and PCO LEO are fixed to estimate the z-ΔPCO GPS . In cases 5 and 6, we determine z-ΔPCO LEO or z-ΔPCO GAL jointly with z-ΔPCO GPS while only fixing PCO GAL or PCO 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-ΔPCO . For example, GEL-GE means that GPS (G), Galileo (E), and Swarm (L) satellites are all included in the processing and z-ΔPCO GPS (-G) and z-ΔPCO GAL (-E) are estimated, while only the PCOs of Swarm satellites are fixed.

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.

Swarm orbit quality
For the gravitational constraint strategy, we included the three spacecraft of the Swarm mission, which is a minisatellite 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 ( ± 19.5 mm). With similar numbers of observations, the SLR residuals (with variation) of Swarm-A and Swarm-C are 3.7 ± 25.1 mm and 2.4 ± 25.0 mm. Compared with previous studies, the orbit quality of our solution is at a comparable level.

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 PCO GPS , the GSA calibrated PCO GAL and the ground calibrated PCO 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-PCO 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.

Results
We analyze the results from three aspects. Firstly, considering the relationship that 130-mm error in GPS z-PCO GPS leads to one ppb terrestrial scale (Zhu et al. 2003), we discuss both the estimated z-ΔPCO 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-ΔPCO GPS and scale. The variation of the estimated daily z-ΔPCO values, the formal error of z-Δ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-ΔPCO estimated by fixing only the PCO GAL in GEL-GL will be analyzed. The impact of fixed PCO GAL on the estimation of z-ΔPCO LEO is shown. At last, mainly based on GEL_GE, the z-ΔPCO GAL estimated by fixing only the PCO 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-1PCO GPS and terrestrial scale
In Fig. 5, the satellite-specific z-ΔPCO 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-ΔPCO 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.
In the case G-G, the estimated z-ΔPCO 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-ΔPCO 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-ΔPCO 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-PCO 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-ΔPCO GPS values have differences of about 30 mm for all satellites. The empirical STD of z-ΔPCO GPS , the formal error of z-ΔPCO 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 Initial epoch state vector and five solar radiation pressure parameters (initial orbital elements are generated from broadcast ephemeris) LEOs orbits Initial epoch state vector; piece-wise empirical force (90-min interval) and atmosphere drag (four hours interval) parameters (initial orbital elements are generated from official products) GPS, Galileo, and LEOs PCOs satellite-specific daily solution of ionosphere-free combined PCOs in the z-direction; fixed to a priori values or freely estimated depending on the case Earth rotation Rotation pole coordinates and UT1 for 24-h intervals, piece-wise linear modeling Tropospheric delay For each ground station; piece-wise constant zenith delays for 1-h intervals; piece-wise constant horizontal gradients for 4-h intervals Satellite and receiver clocks Epoch-wise; pre-eliminated 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-ΔPCO GPS values in GEL-GE are collectively larger than that of GL-G by about 10 mm. The empirical STD of z-ΔPCO GPS , the formal error of z-ΔPCO 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-ΔPCO 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-ΔPCO 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-ΔPCO 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-ΔPCO 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.
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-Δ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-ΔPCO 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 PCO 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 PCO 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-ΔPCO for both GPS and LEOs. In GEL-GE, with only three PCO-fixed LEOs, the correlation between z-Δ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-ΔPCO 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-ΔPCO GPS and 0.1 ppb larger empirical STD of the scale. Moreover, the correlation coefficient between z-ΔPCO 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-ΔPCO 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-ΔPCO GPS derived by GE-G is Fig. 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-ΔPCO , for example, GEL-GE means that GPS, Galileo, and Swarm satellites are included and z-ΔPCO GPS and z-ΔPCO GAL are estimated Table 4 The estimated z-ΔPCO 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 esti-mated 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 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-1PCO LEO
In the case GEL-GL, z-ΔPCO GPS and z-ΔPCO 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-ΔPCO LEO . The mean values and the empirical STD are −2.2 ± 2.5 mm, −2.6 ± 2.1 mm, and −1.1 ± 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-ΔPCO LEO on those days. It also affects all the time series of z-ΔPCO GPS , z-ΔPCO 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.
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-ΔPCO GPS and z-ΔPCO 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-ΔPCO 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-ΔPCO 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-ΔPCO LEO to the values offered by ESA. The time series of z-ΔPCO 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-ΔPCO 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.

Estimated z-1PCO GAL
In the case GEL-GE, z-ΔPCO GAL and z-ΔPCO 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-ΔPCO GAL with respect to the GSA values. Since the z-ΔPCO 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-ΔPCO GAL in GEL-GE are small. The satellitespecific z-ΔPCO GAL values are presented as blue bars 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-ΔPCO GPS directly by onboard GPS observations; therefore, the variations of z-ΔPCO GPS are smaller. Evaluating the impact on the whole Galileo constellation, the empirical STD of z-ΔPCO GAL averaged over all satellites is only 10 mm larger than that of GPS. The formal errors of z-ΔPCO GAL and the empirical STD of the formal errors are only about 6 mm and 2 mm larger than those of the estimated z-ΔPCO GPS . In general, the gravitational constraint on the LEOs acts on the estimation of z-ΔPCO GAL . We also see some unexpected phenomena on satellite E102. We expected a systematic change of the estimated z-ΔPCO GAL due to the scale change, but the absolute value of the estimated z-ΔPCO GAL of E102 is much larger than all the other satellites. For example, the z-ΔPCO 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).
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).
We also plot the time series of system-wise z-ΔPCO GAL and z-ΔPCO GPS in Fig. 12. As explained above, the variation of the z-ΔPCO GAL time series is slightly larger (10 mm) than that of z-ΔPCO GPS . However, due to the fixed LEO PCOs and the same ground stations, they agree with each other.

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, Fig. 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 Fig. 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 Fig. 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 we found a −186 ± 25 mm z-PCO correction with respect to the IGS values, and a +1.55 ± 0.22 ppb terrestrial scale with respect to the IGS14. The results of the gravitational constraint method are −221 ± 37 mm for the z-PCO and +1.72 ± 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-ΔPCO GPS and scale, the formal error of z-ΔPCO 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 Swarmbased solution leads to better results.
The joint estimation of z-ΔPCO GPS and z-ΔPCO LEO by only fixing PCO GAL showed that the z-ΔPCO GPS and the derived scale factor are very close to the solutions derived by the case including GPS and Galileo and fixing the PCO GAL . Consequently, the constraint from Galileo is very strong and is nearly unaffected by including LEOs. The z-ΔPCO 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-ΔPCO LEO time series. This is also visible in the time series of z-ΔPCO GPS , z-ΔPCO 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-ΔPCO GAL estimated by only fixing PCO 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-ΔPCO 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-ΔPCO GAL are both worse than those of the simultaneously estimated z-ΔPCO 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. Acknowledgment Wen Huang is financially supported by the Chinese Scholarship Council. The authors want to thank ESA, GSA, and IGS for offering the data and products.
Funding Open Access funding enabled and organized by Projekt DEAL.
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://creat iveco mmons .org/licen ses/by/4.0/.