Performance of Galileo satellite products determined from multi-frequency measurements

Each Galileo satellite provides coherent navigation signals in four distinct frequency bands. International GNSS Service (IGS) analysis centers (ACs) typically determine Galileo satellite products based on the E1/E5a dual-frequency measurements due to the software limitation and the limited tracking capability of other signals in the early time. The goal of this contribution is to evaluate the quality of Galileo satellite products determined by using different dual-frequency (E1/E5a, E1/E5b, E1/E5, E1/E6) and multi-frequency (E1/E5a/E5b/E5/E6) measurements based on different sizes of ground networks. The performance of signal noise, the consistency of frequency-specific satellite phase center offsets and the stability of satellite phase biases are assessed in advance to confirm preconditions for multi-frequency processing. Orbit results from different dual-frequency measurements show that orbit precision determined from E1/E6 is clearly worse (about 35%) than that from other dual-frequency solutions. In view of a similar E1, E5a, E5b and E6 measurement quality, the degraded E1/E6 orbit performance is mainly attributed to the unfavorable noise amplification in the respective ionosphere-free linear combination. The advantage of using multi-frequency measurements over dual-frequency for precise orbit determination is clearly visible when using small networks. For instance, the ambiguity fixing rate is 80% for the multi-frequency solution while it is less than 40% for the dual-frequency solution if 150 s data sampling is employed in a 15-station network. Higher fixing rates result in better (more than 30%) satellite orbits and more robust satellite clock and phase bias products. In general, satellite phase bias products determined from a 20-station (or more) network are precise enough to conduct precise point positioning with ambiguity resolution (PPP-AR) applications. Multi-frequency kinematic PPP-AR solutions always show 5–10% precision improvement compared to those computed from dual-frequency observations.

Pilot Project) metadata SINEX file (Steigenberger and Montenbruck 2022b). As confirmed in Steigenberger et al. (2018), the impact of the antenna thrust on Galileo satellite orbits can be larger than 2 cm in the radial direction for the FOC (Full Operational Capability) satellites. In addition to the officially published metadata, radiator emission power estimated and navigation antenna thermal properties of Galileo satellites are published in Sidorov et al. (2020), Bhattarai et al. (2022), Dilssner et al. (2022) and Duan and Hugentobler (2022). With all this metadata, Galileo satellite orbits can be determined with a precision of better than 2 cm in the radial component.
The coherent Galileo navigation signals are transmitted in four frequency bands: E1, E5a, E5b and E6 based on the same onboard master clock. One signal (e.g., E5a) may contain several components, consisting of at least one pair of pilot (Q) and data (I) channels. The E5a and E5b signals are part of the E5 signal in its full bandwidth, an AltBOC signal (ICD 2021;Lestarquit et al 2008). This AltBOC signal can also be tracked as a single signal, providing a very large signal bandwidth (around the 1191.795 MHz center frequency), and thus is less affected by multipath errors (Falcone et al. 2017;Prochniewicz and Grzymala 2021;Diessongo et al. 2014). All the Galileo frequency bands are part of the spectrum assigned by the International Telecommunication Union (ITU), for Radio Navigation Satellite Services (RNSS). E1, E5a, E5b are included in the spectrum for ARNS (Aeronautical Radio Navigation Services), allowing dedicated safety-critical applications. The E6 signal is transmitted in the frequency band 1260 to 1300 MHz, sharing with radar systems of the radio navigation and radiolocation service and thus could be affected by different types of radio-frequency interference (i.e., radar, amateur radio and other unidentified low-power sources).
In view of the community with the GPS L1 and L5 frequencies and the fact that only E1 and E5a Galileo signals were tracked by GNSS receivers in the early time, most of the Galileo satellite products within the IGS are traditionally based on the dual-frequency (E1/E5a) measurements and refer the estimated clock offsets to the ionosphere-free linear combination of those signals. However, with the continued modernization of GNSS receivers, more than 100 IGS stations are able to track all-frequency Galileo signals since the beginning of 2021. Satellite PCO/PCVs (phase center offset/variations) of all the Galileo carrier frequencies are physically calibrated, playing an important role for the independent determination of the GNSS scale . The corresponding PCO/PCVs are provided as part of the igsR3.atx antenna model ), but the number of the provided receiver antenna types is still limited. Geo++ (robot method) and University of Bonn (chamber method) have made great efforts to complement the calibration of the remaining receiver antenna types (Wübbena et al. 2019;Kersten et al. 2022). TUG (Technische Uni-versität Graz), one of the IGS ACs in repro-3, employs multi-frequency raw measurements in an uncombined mode to determine all the satellite products (Strasser et al. 2019). However, Galileo E6 signals were not involved due to the insufficient ground tracking capability during the IGS repro-3 time period. With multi-frequency satellite products, advantages in PPP/PPP-AR (Zumberge et al. 1997) and PPP-RTK (Real-time Kinematic) are presented in Li (2018), Li et al. (2017), Zhang et al. (2019Zhang et al. ( , 2020, Li et al. (2022a), Li et al. (2022b), Wang and Rothacher (2013), Geng et al. (2022) and Jiang et al. (2023). As an example, the convergence time in the multi-frequency kinematic PPP-AR is reduced by about 15% compared to the dual-frequency results.
The navigation messages of the Galileo system are updated every 10 min, resulting in a global mean RMS (root mean square) of SISRE (signal-in-space range error) of about 10-20 cm in real time . Overall, positioning errors at the few-decimeter (30-40 cm) level could be achieved in kinematic PPP solutions with Galileo broadcast ephemerides (Carlin et al. 2021). The Galileo High Accuracy Service (HAS), aiming at obtaining more precise global PPP results, is currently under development and once completed will transmit precise orbit, clock and bias products, as well as atmospheric data (regional level) through the Galileo E6-B signal to the public free of charge. The HASsupported constellations will be GPS and Galileo, based on L1/L2 and E1/E5a/E5b/E5/E6 signals, respectively. The preliminary performance of the HAS products for the Galileo system is evaluated in Fernandez-Hernandez et al. (2022) and Hauschild et al. (2022). Kinematic PPP results are shown to fulfill the targeted two-decimeter user horizontal accuracy. However, the phase bias products were not yet available in the initial experimental HAS products. In this context, the precision of satellite phase bias products based on a small ground network needs further assessment.
With the above background, the present contribution aims at evaluating the performance of satellite products determined from different dual-frequency (E1/E5a, E1/E5b, E1/E5, E1/E6) and multi-frequency (E1/E5a/E5b/E5/E6) measurements based on different sizes of ground networks. The methodology and processing strategy are described in Sect. 2. An analysis of the quality and consistency of measurements on individual frequencies is provided in Sect. 3. Data and standards applied in the processing are introduced in Sect. 4. Detailed analyses of Galileo satellite products, i.e., orbit quality, precision of clock offsets and phase biases, and performance in kinematic PPP-AR are presented in Sect. 5. Finally, a summary and conclusions are presented based on the discussed results.

Methodology
In this section, we describe the procedures applied in the dualfrequency and multi-frequency processing. Pseudorange (P) and carrier phase (L) measurements between one receiver (subscript r ) and one Galileo satellite (superscript s) on frequency l are expressed as Here ρ s r denotes the geometric distance between the antenna reference points of the satellite and receiver, while the related PCO/PCV corrections are expressed by δ s r ,l . dt r and dt s denote the receiver and satellite clock offsets, d and b represent the pseudorange and phase biases, c denotes the speed of light in vacuum. The first-order ionospheric delay (I s r / f 2 l ) affects pseudorange and carrier phase with opposite signs and varies with the inverse square of the signal frequency f l . For Galileo, frequency 1, 5, 6, 7, 8 denote E1, E5a, E6, E5b, E5, respectively, in accord with the frequency band numbers employed in the Receiver Independent Exchange (RINEX) format (Romero 2021). The neutral tropospheric delay (T s r ) is the same for all the frequencies observed by a given station. λ l denotes the wavelength of frequency l and N s r ,l represents the phase ambiguity. Other corrections, e.g., the phase windup correction (Wu et al. 1993) are modeled and are not shown in the equation. e s r ,l and ε s r ,l contain unmodeled error sources of pseudorange and phase observations, most notably, measurement noise, multipath and residual tropospheric delay modeling error.

Dual-frequency procedure
The first-order ionospheric delay is dispersive and is usually eliminated by forming the ionosphere-free (IF) linear combination of two frequencies. For instance, IGS ACs determine precise Galileo satellite products based on E1/E5a dualfrequency observations (Laurichesse et al. 2009;Schaer et al. 2021;Katsigianni et al 2019;Deng 2022;Geng et al. 2019). To support PPP-AR applications, ACs may provide users directly with satellite wide-lane and narrow-lane biases or incorporate narrow-lane biases into satellite clock offsets (e.g., Center National d'Etudes Spatiales/Collecte Localization Satellites, CNES/CLS products, Loyer et al. 2012). Alternatively, ACs may provide users with observablespecific signal biases (OSB) on each code and phase signal (e.g., Center for Orbit Determination in Europe, CODE prod-ucts, Villiger et al. 2019;Schaer et al. 2021). Theoretically, there should be no difference between these two types of products and both are proven to perform very well in the single-receiver ambiguity resolution (Teunissen and Khodabandeh 2015;Banville et al. 2020;Montenbruck et al. 2021;Mao et al. 2021;Fernández et al. 2022).
Same as CODE, we determine Galileo satellite orbit, clock offset and OSB products. The IF linear combination of code and phase equations for E1 and E5a frequency are expressed as The measurement noise in this IF linear combination is amplified by about a factor of 2.6 compared to that of a singlefrequency measurement. The IF combined phase ambiguity cannot be expressed in the form λ IF N IF where N IF is an integer ambiguity with a reasonable wavelength. The classical ambiguity resolution approach for Eq. (2) is to introduce the integer wide-lane ambiguity into the IF ambiguity and then fix the deduced narrow-lane ambiguity (Ge et al. 2005). The overall processing scheme applied in our study is illustrated in Fig. 1. First, we determine OSB values for all the code signals. In the IGS ground network, some receiver types (e.g., Septentrio) track only the pilot (Q) component while some (e.g., Trimble) track only the combination of pilot and data (X) components. Sleewagen and Clemente (2018) show that differences of phase biases for pilot and data tracking are confined to a few milli-cycles for all current pilot-data signals. As such, it is reasonable to ignore also the possible phase bias difference between Q and X tracking of the Galileo signals. However, differences of code biases between Q and X for Galileo IOV (In-Orbit Validation) satellites could cause a difference of up to 0.05 wide-lane cycles in the HMW (Hatch-Melbourne-Wübbena: Hatch 1983;Melbourne 1985;Wübbena 1985) linear combination . Therefore, we adopt the pilot-only tracking observations (as identified by RINEX observation codes C1C and C5Q, Romero 2021) as the datum in the clock estimation and assume that code biases of satellites (d s IF ) and receivers (d r ,IF ) in the IF linear combination for these two signals are equal to zero. Satellite differential code biases (DCB) of C1C-C5Q are first estimated by forming the geometry-free pseudorange linear combination and applying a global ionosphere model (Montenbruck et al. 2014;Wang et al. 2016Wang et al. , 2020   Adding the clock datum condition that d s IF of C1C and C5Q is equal to zero, OSB products of C1C and C5Q are solved for each satellite. Then, we compute OSB products for C1X and C5X signals as part of the clock estimation mixing receivers tracking Q and X in the network. Satellite orbits in this step can be fixed to the a priori solutions, i.e., predicted orbits or float-ambiguity solutions. PCOs are applied both in the DCB and the clock estimation. Second, the HMW linear combination is used to resolve wide-lane (WL) ambiguities HMW(P s r ,E1 + OSB P,E1 , P s r ,E5a + OSB P,E5a , L s r ,E1 , L s r ,E5a ) where P s r ,E1 , P s r ,E5a , L s r ,E1 , L s r ,E5a represent code and phase measurements on frequency E1 and E5a. The determined code OSB products (in step 1) OSB P,E1 (OSB C1C or OSB C1X ) and OSB P,E5a (OSB C5Q or OSB C5X ) are applied in advance to correct the measured pseudoranges prior to forming the HMW linear combination. Daily constant WL biases of all the satellites (b s WL ) and receivers (b WL,r ) are estimated while fixing wide-lane ambiguities to integers, employing a zero mean condition of all satellite HMW biases. The estimated satellite WL biases can be assumed to contain only phase biases of E1 (OSB L1C =OSB L1X ) and E5a (OSB L5Q =OSB L5X ) as code OSBs are corrected in advance. This assumption is consistent with clock parameters and is one of the conditions to solve phase biases on each frequency. Third, the integer wide-lane ambiguity is introduced into the IF ambiguity to compute integer-ambiguity satellite orbits. The IF phase equation in Eq.
(2) is expressed as Heredt r anddt s contain both the clock offset and the phase bias of the receiver and satellite, respectively. N s r ,1 denotes the narrow-lane ambiguity with a wavelength of about 11 cm. With only phase measurements in this step, we need to select a reference ambiguity for each receiver and fix it directly to the nearest integer value to cope with the singularity between all the ambiguities and the receiver clock offsets. To cancel out potential receiver-specific errors in the undifferenced float ambiguities, we form pass-wise satellite-satellite single differences between float narrow-lane ambiguities of the same receiver and resolve single-differenced phase ambiguities. Clock parameters are pre-eliminated every epoch and we do not recover them in this step. Once all the narrow-lane ambiguities are fixed, daily normal equations are stacked to produce 3-day-arc solutions, comprising satellite orbits and station-related parameters (i.e., station coordinates and tropospheric delays).
In step 4, clock offsets are computed while fixing satellite orbits and station-related parameters that are determined in step 3 as known, leading to Hereρ denotes the geometric distance introducing satellite orbits and station-related parameters as known. dt r and dt s contain the IF code biases (C1C/C5Q) and this defines the datum of clock offsets. If C1X/C5X signals in the ground network are used, the corresponding satellite code OSB products have to be used to make all the signals consistent with the definition of the clock datum. Rigorously speaking, the phase clock offsetsdt r anddt s differ from those of the code model, but the difference is within one narrow-lane cycle and can typically be ignored considering the noise and multipath error of pseudorange measurements. Thus, we estimate common clock offsets for code and phase measurements. In step 5, satellite phase biases are estimated while fixing satellite orbits, clock offsets (dt s ) and station related parameters as known, leading to Here b s IF denotes satellite narrow-lane phase biases, which are estimated while fixing narrow-lane ambiguities to integers ( N s r ,1 ). Adding the condition that satellite WL biases (b s WL ) are known, satellite phase biases (OSB) are solved. We need to mention that the estimated satellite phase biases and clock offsets need to be used consistently in PPP-AR.
In step 6, the data sampling is changed from 150 to 30 s and PPP-AR processing of each station is performed in parallel to fix more ambiguities, because ambiguity resolution always benefits from the use of higher sampling data. For this, satellite clock products are densified to 30 s interval (Bock et al. 2009). If new phase ambiguities are fixed, steps 4-6 are iterated to compute a new set of satellite products.

Multi-frequency procedure
In our multi-frequency ambiguity resolution, we introduce all the integer phase ambiguities of E1, E5a that are fixed in the IF linear combination as known to reduce the correlation between parameters. All the other phase ambiguities in the uncombined phase equation can then be reasonably resolved. Only in very rare cases when no ambiguities are fixed for one satellite pass in the IF linear combination, correlations between slant ionospheric parameters and phase ambiguities still exist. To cope with these cases, we employ additionally a relative constraint on the change of the ionospheric slant delay with time. The detailed processing procedures are given in Fig. 2.
First, we compute code OSBs of C1C and C5Q in the same way as in the dual-frequency procedure. Then, we employ the same clock datum (C1C, C5Q) and compute code OSBs of C1X, C5X and also of signals of other frequencies as part of the clock estimation using uncombined multi-frequency measurements (Eq. 1).
Second, we introduce integer phase ambiguities of E1, E5a into the uncombined multi-frequency phase equation and determine ambiguity-fixed satellite orbits. The multifrequency phase observation model is given by where N s r ,1 , N s r ,5 denote the integer ambiguities as obtained from the wide-and narrow-lane ambiguities (N WL = N 1 − N 5 ; N NL = N 1 + N 5 ) that were previously fixed using the HMW and IF combination of E1 and E5a observations. Clock parametersdt r anddt s include phase biases of E1/E5a. Satellite and receiver phase biases b r ,6,7,8 , b s 6,7,8 will be shown to be mostly constant (at the phase measurement noise level) over one day in Sects. 3 and 4. Therefore, we do not estimate explicit receiver phase biases for frequency 6, 7 and 8, which will be fully absorbed by the corresponding float ambiguity parameters, e.g.,Ñ s r ,6 includes N s r ,6 and b r ,6 . By forming single differences between float ambiguities of the same frequency and receiver, receiver phase biases are canceled out and phase ambiguities can be fixed to integers. As only phase measurements are used in this step to determine satellite orbits (3-day-arc) and station-related parameters, we do not need to recover the pre-eliminated clock parameters.
With the determined satellite orbits and station-related parameters, we estimate satellite clock offsets and phase biases sequentially following similar procedures as in the dual-frequency processing. Steps 3-5 in Fig. 2 are iterated until no new ambiguities can be fixed. We would like to mention that multi-frequency biases determined in the uncombined mode can be used to fix phase ambiguities in any linear combinations (i.e., HMW and narrow-lane).

Signal and measurement performance
Leaving aside the availability of pilot and data channels, which may be tracked in either a combined or a standalone mode, the Galileo system offers a total of five open signals in different frequency bands. As such, at total of five distinct measurements are commonly made available by geodetic receivers for use in dual-and multi-frequency orbit determination, point positioning and timing. Prior to using multi-frequency Galileo observations for those applications, we evaluate and compare the quality of pseudorange and carrier measurements on the individual frequencies.
Considering the signal-in-space itself, the E5a, E5b, and E6 signals are transmitted with a power and antenna gain pattern that ensures a common total received minimum signal power of −157.25 dBW (ICD 2021) including pilot and data components. For the E1 open service signal a 2 dB lower value applies, which would translate into moderately (≈ 20%) increased thermal receiver tracking noise under otherwise equal conditions (Ward 2017). Vice versa, tracking of the entire E5 signal benefits from the combined E5a and E5b power and offers a 30% reduced noise level compared to E5a or E5b tracking with identical correlator and tracking loop settings. Ranging codes in the E1, E6, and E5a/b signals exhibit frequencies of 1.023 MHz, 5.115 MHz, and 10.23 MHz, respectively, which corresponds to chip lengths between about 30 and 300 m. However, in contrast to the binary or quadrature phase shift keying (BPSK/QPSK) modulation used in E6 and E5a/b, a binary offset carrier modulation (BOC) is used in E1, which helps to reduce the code tracking noise and multipath sensitivity despite the longer chip length (Betz 2016).
Despite obvious differences in the properties of the transmitted signals, the practical performance of actual GNSS measurements appears mostly driven by the specific characteristics of the user equipment (antenna, receiver) and the local environment (multipath, interference). For a practical assessment, we therefore made use of code and carrier phase observations from different sites to compare the measurement quality for the five different frequencies. All sites are equipped with PolaRx5 receivers that constitute the majority of IGS stations in the networks used in our subsequent processing, but are connected to different antennas (Javad RingAnt-DM, Septentrio PolaNt Choke Ring B3E6, Leica AR25).
For all three antenna types, carrier-to-noise density (C/N 0 ) ratios of E5a, E5b and E6 observations agree within about 1 dB (Fig. 3), which suggests a similar amount of carrier phase tracking noise for the respective measurements. As expected, the combined E5(a+b) AltBOC tracking results in a 3 dB higher C/N 0 , while the E1 measurements show a roughly 2 dB lower C/N 0 than E5a/E5b/E6 as a result of different receiver antenna gain patterns. The analysis of code measurement errors is based on a ionosphere-corrected code-carrier difference (Rocken and Meertens 1992;Kee and Parkinson 1994) that is commonly known as "multipath" combination and provides a measure of the combined receiver noise and multipath for pseudorange measurements on individual signals. Results shown in Fig. 4 demonstrate a good overall consistency of the measurement quality among the various signals. However, the wideband E5 AltBOC signal stands out with a 30-50% lower error budget that relates to the reduced thermal noise and the superior multipath resistance. E1 measurements, on the other hand, exhibit a slightly increased error budget and a higher site dependency, which reflects the scatter in the E1 gains of the employed antennas as well as an increased multipath susceptibility.
For the assessment of carrier phase noise observations, we make use of 1 Hz high-rate phase measurements ϕ and form time/frequency differences All values in mm of observations at epochs t i and t i−1 and signals a and b. These double-differences are free of geometry and clock contributions and represent the sum of the between-epoch and between-frequency difference of the ionospheric path delays and the corresponding double-difference of carrier phase noise ε. For a sufficiently quiet ionosphere, the change of path delays within a 1 s interval can be neglected relative to the contribution of phase noise. Furthermore, we may assume uncorrelated measurement errors across frequencies and epochs in view of independent front-ends and phase-locked loop (PLL) bandwidths that notably exceed the sampling rate. As a result, the variance of the epoch/frequency double difference equates to twice the sum of the carrier phase noise variances for the individual signals considered in the frequency difference. Using a 3cornered hat approach, the noise standard deviations for three signals a, b, c are finally obtained from as well as the corresponding equations with cyclically permuted indices. Example data for a single station at different elevations (and thus C/N 0 values) are collated in Table 1. As in the analysis of pseudorange measurements, the results indicate a similar quality of carrier observations from all five signals with a slightly reduced noise of the E5 AltBOC observations.
Despite a comparable quality of the individual multifrequency observations, care must be taken that the use of a dual-frequency combination for eliminating first-order ionospheric path delays will result in different levels of noise amplification. Considering combinations of observations with frequency f a in the lower L-band (i.e., E1) and f b in the upper L-band (i.e., E5a, E5, E5b, or E6), a wider separation of the two frequencies obviously results in a smaller noise and vice versa. With α values in the range of 1.3-−1.4, the IF linear combinations of E1 with E5a, E5, and E5b exhibit a roughly 2.7 times higher noise than individual single-frequency measurements, while the E1/E6 combination exhibits a 3.5 times higher noise.
Further to the assessment of carrier phase noise, we evaluate the consistency of the various IF carrier phase combinations, which is of practical relevance for the compatibility of clock and bias products in multi-GNSS solutions. As discussed in Simsky (2006) and Montenbruck et al. (2012a), the difference represents a geometry-and ionosphere-free triple-frequency linear combination that is nominally constant when correcting the observations for antenna patterns and wind-up. As such, the combination is well suited to assess temporal variations in carrier phase biases that were, for example, identified in GPS Block IIF satellites and related to sunincident-dependent temperature variations (Montenbruck et al. 2012b). For Galileo, we specifically evaluated the consistency of the E1/E5b, E1/E5, and E1/E6 combinations with that of E1/E5a, which is most widely used for generation of Galileo orbit and clock products within the IGS. A globally-distributed 15-station network with adequate depthof-coverage was used to estimate epoch-wise values of the triple-carrier combinations along with site-and pass-specific carrier phase biases over a continuous 5-day arc. Other than for GPS IIF satellites, orbit-periodic contributions are essentially absent and the respective amplitudes of firstto fourth-order harmonics are confined to 1-3 mm in all cases. On the other hand, obvious drifts may be recognized for selected satellites (Fig. 5). These are most pronounced for the early set of in-orbit validation (IOV) satellites (i.e., space vehicle numbers SVN E101 to E103) and may reach a magnitude of roughly 25 mm/d for the E1/E6-E1/E5a difference. For full operational capability (FOC) satellites, in contrast, the drift does not exceed a magnitude of 5 mm/d with the exception of E202 and E217. Likewise, drifts of triple-frequency combinations of E1, E5a, and either E5b or E5 amount to less than a few mm on all satellites over one day. Given the fact that drifts between the various IF carrier phase combinations are observed in the analysis of multi-station data and exhibit a clear satellite-dependence, it appears plausible to suspect small inconsistencies in the onboard signal generation (including modulation and up-conversion) as the root cause of the observed inter-carrier drifts. These are obviously most pronounced when comparing E1, E5a, and E6 signals, while the other triple-frequency combinations show Fig. 5 Bar chart of satellite-specific drifts of the ionosphere-free E1/E6, E1/E5b, and E1/E5 carrier phase combinations relative to the E1/E5a combination a better compatibility due to the fact that E5, E5a, and E5b are generated as a compound wide-band signal.
For dual-frequency processing, linear drifts between carriers of individual signals are essentially absorbed in the receiver clock estimate and therefore hardly-discernible for common users. However, inter-carrier-phase inconsistencies will inevitably show up in triple-or multi-signal processing where they require due consideration in terms of timevarying inter-signal biases.

Station networks and processing standards
Six ground networks (with 90, 50, 40, 30, 20, 15 stations) are employed (as shown in Fig. 6) to assess the performance of Galileo satellite products determined by different observations and networks. All the stations are able to track all-frequency Galileo measurements. Q and X signals are mixed in the networks. All the stations are included in the IGS14 frame with coordinates (at the reference epoch) and velocities given as known. The analysis period covers 30 days from 2021-01-01 to 2021-01-30. PCO/PCVs of receiver antennas are calibrated in the igsR3.atx file. The 90-station network is used to generate a set of reference solutions for E1/E5a dual-frequency measurements. Other smaller networks (≤ 50) are used to determine satellite products (e.g., satellite phase biases) based on different dual-frequency (E1/E5a, E1/E5b, E1/E5, E1/E6) and multifrequency (E1/E5a/E5b/E5/E6) measurements. Table 2 summarizes important modeling options and details of the estimated parameters. A modified version of the Bernese GNSS Software 5.3 ) is used to perform all the computations. The data sampling interval is set to 150 s in the network processing and 30 s in the PPP-AR processing. The elevation cutoff is set to 5 • and an elevation-dependent weighting is employed. The a priori physical macro-model  and the 7-parameter ECOM2 model (Arnold et al. 2015) are used to describe solar radiation pressure and thermal thrust forces. Earth radiation (Rodriguez-Solano et al. 2012) is modeled by using the published Galileo metadata (https://www.euspa. europa.eu/). Antenna thrust effect is considered based on the transmit power measured by DLR ). Station-related parameters, satellite orbits, Earth rotation, and geocenter parameters are based on 3-day-arc solutions, while clock offsets and biases are obtained in 1-day-arc solutions. Phase ambiguities are fixed to integers using the Bernese SIGMA method .

Consistency of frequency-specific PCOs for Galileo satellites
Satellite-and frequency-specific PCOs obtained from chamber calibrations have been published by the EUSPA (EUSPA 2022) and were adopted in the igsR3.atx antenna model. In order to evaluate their consistency, satellite antenna PCOs have been estimated with the NAPEOS software (Springer 2009) for the E1/E5a, E1/E5b, and E1/E6 IF linear combinations. A dedicated global network of 148 stations was used and the terrestrial scale was fixed to ITRF2020 (Altamimi et al. 2022). Further details on the PCO estimation are given in Steigenberger and Montenbruck (2022a). The mean PCO differences between the chamber calibrations and the estimates from the time interval 1 July until 31 December 2021 as well as their standard deviations are given in Table 3. The x-and y-components of the estimated PCOs agree on the millimeter level for the different IF linear combinations. The x-component shows a systematic offset of about 2 cm w.r.t. the calibrated PCOs whereas the y-component agrees almost perfectly with them. The z-component of the estimated PCOs differs by −11 to −16 cm from the chamber calibrations. The overall systematic bias translates into a small-scale inconsistency between the ITRF and station positions obtained with the chamber calibrations ). The largest z-PCO differences are present between E1/E5b and E1/E6 and amount to 4.7 cm. The major part of this difference generates a constant bias that can be absorbed by the different satellite clock offset estimates for different linear combinations or by the satellite bias estimates of the multi-frequency processing. The boresight-dependent effect scales with 1 − cos θ  where θ is the maximum nadir angle at the satellite. With θ GAL ≈ 12 • , this effect is only 2 %. For the PCO differences in Table 3, a maximum value of less than 1 mm is obtained, which is negligible in view of the phase measurement noise. Therefore, it is reasonable to make use of the unmodified EUSPA PCO values for a consistent multi-frequency processing of Galileo observations.

Stability of phase biases
The stability of phase biases has to be known in advance in the multi-frequency processing to determine the sampling interval for which satellite and receiver phase biases can be treated as constant values. For instance, phase biases of GPS Block IIF satellites in the L5 frequency show apparent temperaturedependent variations (Montenbruck et al. 2012a) and cannot be determined as a daily constant bias in the multi-frequency processing (i.e., triple-frequency or more). Similar works for GPS and Multi-GNSS of two frequencies are analyzed in Li et al. (2022c), Håkansson (2017) and Odijk et al. (2017).  Figure 7 illustrates the phase clock differences between different dual-frequency solutions in one day. All the differences are constant at the level of the phase measurement noise. Only phase clock differences of E1/E5a-E1/E6 for the IOV satellites show a clear drift over time. This is in line with what we have observed in Fig. 5. Table 4 shows the mean STD (standard deviation) of phase clock differences over 30 daily solutions for Galileo satellites (IOV and FOC) and different receiver antenna types. STD values from E1/E5a-E1/E5b and E1/E5a-E1/E5 are, in general, smaller than 3.5 mm while those from E1/E5a-E1/E6 can exceed 6 mm for some receivers. As discussed already in Sect. 3, this is attributed to the larger noise amplification factor of E1/E6 (and also E1/E5a-E1/E6). Phase clock differences of E1/E5a-E1/E6 for the IOV satellites show a clear shift, resulting in a clearly larger STD value than that of the FOC satellites. Considering the largest STD value (about 6 mm) to the wavelength (about 20 cm) of each phase ambiguity, we estimate daily constant phase biases for all the satellites and receivers on each phase signal.

Results and discussion
With all the data, settings and conclusions in Sect. 3 and 4, we evaluate the quality and performance of Galileo satellite orbit, clock offset and daily bias products that are determined by different measurements and networks. All the scenarios are described in Table 5. The 90-station network is first used to determine a set of reference solutions for the E1/E5a IF linear combination. Next, the 15-to 50-station networks are used to assess the performance of satellite products determined from different networks and frequencies. Time periods are the same as in Sect. 4.

Orbit quality
Ambiguity fixing rates from different observations and networks are taken as the first indicator to show the quality of the determined satellite products. Figure 8 illustrates the WL, NL, and uncombined ambiguity fixing rates from different solutions. The left figure shows fixing rates from the first iteration after the network solution (i.e., step 3 in Fig. 1), the right figure shows fixing rates from the last iteration after the PPP-AR solution (i.e., step 6 in Fig. 1). Fixing rates of the wide-lane ambiguity are stable for different networks because the HMW linear combination is geometryand ionosphere-free and depends only on observations of a single station. WL fixing rates are around 87% in the first iteration and about 95% in the last iteration due to the use of 30 s sampling data. NL fixing rates from dual-frequency solutions in the left figure (first iteration) drop from more than 70% for the 20-station network to less than 40% for the 15-station network. Compared to E1/E5(a,b), the E1/E6 combination not only suffers from increased noise amplification in the IF combination but also from a shorter NL wavelength (10.5 cm).The ambiguity fixing rate for uncombined multifrequency processing remains around 80% because 2.5 times more measurements are employed. This is the same as shown in the right figure if the higher sampling 30 s data are used in the last iteration. This indicates that our ambiguity resolution method always gets a higher fixing rate if higher sampling and multi-frequency data are employed.
Orbit precision is evaluated by the orbit misclosures (m) between two consecutive 3-day arcs.  where r s i+1 denotes the orbit position vector of satellite s of day i +1. Figure 9 illustrates the 3D (position distance) mean RMS of orbit misclosures determined using different observations and networks. For float-ambiguity processing, only the E1/E5a and E1/E6 results are displayed in the figure as float-ambiguity solutions of E1/E5a, E1/E5b and E1/E5 are similar. The RMS of the E1/E6 solutions is about 10% larger than that from E1/E5a because of the larger noise amplification in the E1/E6 IF linear combination. We also find that phase residuals determined from E1/E6 are about 15% larger. Ambiguity fixed orbits in Fig. 9 are the first-iteration results, corresponding to the fixing rates in the left Fig. 8. In general, the precision of ambiguity-fixed orbits is about two times better than float-ambiguity orbits if 20 (or more) tracking stations are used. The multi-frequency orbit solution shows a clear improvement of about 30% with respect to the E1/E5a dual-frequency solution for the 15-station network due to the higher ambiguity fixing rate.
As fixing rates in the first and the last iteration exhibit clear differences, we iterated the determination of satellite orbits after the PPP-AR step. We find that there is almost no difference in orbit quality between the first and the last iteration if 20 (or more) tracking stations are employed. However, noticeable improvements could be observed for the 15-station network solution. Figure 10 shows 3D RMS of orbit misclosures from the 15-station network solution in the last iteration. Overall, an improvement of about 20% is observed in the last iteration compared to the orbit misclosures in the first iteration.
Satellite laser ranging (SLR) measurements collected by the International Laser Ranging Service (ILRS, Pearlman et al. 2019) are used as an external reference to assess Galileo satellite orbits. For the computation of SLR residuals, stations recommended by Bury et al. (2021) are employed. Station coordinates are fixed to the SLRF2014 and ocean tidal loading is considered with FES2004. Statistics of all the integer-ambiguity orbit solutions are shown in Table 6, outliers exceeding ± 25 cm are rejected. The number of nor-

Clock and phase bias quality
The precision of satellite clock offsets is assessed by comparing clock products determined from different scenarios to the reference solution (90-station network solution). A mean difference of all the satellites every epoch is subtracted from all the clock differences to remove the datum difference. This is partly because that the E1/E5a solutions use the same frequency as that in the reference solution and observations are partially overlapped (the reference network contains all the stations in other networks). The precision of the E1/E6 clock solutions is about two times worse due to the higher noise amplification of E1/E6 and the drifts of clock differences (E1/E5a-E1/E6) for IOV satellites (which increase the STD value of E1/E6 by 5-10%). We would like to mention that the precision of different satellite clock products in Fig. 11 only represents the consistency of individual solutions with respect to the E1/E5a reference solution. However, the precision of clock offsets from different networks of the − 0.9 ± 3.1 − 0.8 ± 3.4 − 0.9 ± 3.2 − 0.8 ± 4.1 − 0.9 ± 3.0 20 − 0.9 ± 3.1 − 0.8 ± 3.0 − 0.9 ± 2.9 − 0.9 ± 3.2 − 0.9 ± 2.9 30 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 3.1 − 0.9 ± 2.9 40 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 3.0 − 0.9 ± 2.9 50 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 2.9 − 0.9 ± 2.9 90 − 0.9 ± 2.9 ----  Fig. 11 Mean STD of clock differences between different solutions and the reference solution same frequencies can still be used to assess the relationships between clock precision and the size of the network. The precision of satellite E1 and E5a bias products is evaluated by comparing individual bias solutions (WL and NL) to the reference solution (90-station network bias solution). The comparison of satellite wide-lane biases is easily possible since they are independent of satellite orbits and clock offsets. However, a direct comparison of satellite narrow-lane biases needs more considerations as these satellite narrowlane biases are correlated with satellite orbits, clocks and integer wide-lane ambiguities. To that end, we compute the total difference of daily orbit differences (in the radial direction), mean clock differences, integer wide-lane differences and IF phase bias differences. A reference satellite is selected and is subtracted (total differences) from other satellites to remove the datum difference. The total difference should in theory be equal to unknown integer cycles of the narrow-lane ambiguity. The fractional part of the total difference is a representation of the difference of the narrow-lane bias. More details of the comparison method for satellite narrow-lane biases are given in . Figure 12 shows the ratio of WL and NL bias differences that are smaller than 0.05 WL or NL cycles for different net-

PPP-AR applications
PPP results based on different measurements are analyzed as an indirect assessment of the quality of satellite products. For this, we compute kinematic PPP/PPP-AR (KPP/KPP-AR) solutions of 6 stations that are not included in any of the networks in Fig. 6. The stations are equipped with different receiver and antenna types and also globally distributed. PCO/PCVs of receiver antennas are calibrated in the igsR3.atx file. All the station-related model parameters (i.e., elevation cutoff and ocean loading) are the same as introduced in Sect. 4. Figure 13 shows kinematic positioning results compared to the static PPP-AR solution for one station in North, East and Up components. E1/E5a dual-frequency and E1/E5a/E5b/E5/E6 multi-frequency KPP/KPP-AR solutions based on satellite products from the 20-station network are taken as an example. E1/E5a dual-frequency KPP-AR results show clear improvements of about 35% with respect to the E1/E5a KPP solutions in the East component (similar as shown by Bertiger et al. 2010), the use of multi-frequency observations (KPP-AR) shows further improvement of about 7%. Also, multi-frequency positioning results are more robust if the number of dual-frequency observations is low at some epochs. The station-averaged STD values of KPP/KPP-AR solutions using different observations and satellite products are shown in Table 7. KPP-AR results using satellite products from the 15-station network are not as good as other solutions. For results based on satellite products from other larger networks, multi-frequency KPP-AR positioning results always show further improvement of 5-10% over the corresponding dual-frequency KPP-AR positioning results.

Summary and conclusion
The Galileo system has been providing initial operational global positioning and timing services since 2017. From the beginning of 2021, more than 100 globally distributed ground tracking stations are able to track Galileo all-frequency signals. To assess the quality of satellite products determined by different dual-frequency and multi-frequency measurements, and to suggest a minimum number of tracking stations ensuring a precise determination of phase bias products, we determine Galileo satellite products based on different measurements and networks.
We first assess signal and measurement performance of different frequencies for Galileo. C/N 0 ratios of E5a, E5b and E6 observations agree within about 1 dB, the combined E5(a+b) AltBOC tracking results in a 3 dB higher C/N0, while the E1 measurements show a 2-5 dB lower C/N 0 than E5a/E5b/E6 as a result of different receiver antenna gain patterns. The wideband E5 AltBOC signal stands out with a instead of 30-50% lower code measurement error budget that relates to the reduced thermal noise and the superior multipath resistance. E1 measurements, on the other hand, exhibit a slightly increased error budget and a higher site dependency. In view of a similar overall performance of measurements for all individual signals, the increased orbit misclosures (35%) in E1/E6 solutions are mainly attributed to the larger noise amplification factor in the respective IF linear combination.
Then, we confirm that the manufacturer-provided, frequencyspecific satellite PCOs and PCVs may be used in the multi-frequency processing to determine a common and self-consistent set of satellite products. The stability of frequency-specific phase biases is evaluated by both forming a triple-frequency linear combination and computing phase clock differences determined by two IF linear combinations. Obvious drifts of up to 25 mm/d for the E1/E5a-E1/E6 differences are observed for the IOV satellites. Other triplefrequency combinations and phase clock differences show a better compatibility due to the fact that E5, E5a and E5b are generated as a compound wide-band signal. Considering the wavelength (about 20 cm) of each phase ambiguity, we may nevertheless estimate daily constant phase biases for all the satellites on each frequency.
Ambiguity resolution results from different observations and networks demonstrate that ambiguity fixing rates always benefit from the use of more observations (e.g., high sampling or multi-frequency data), particularly for the 15-station network solution. For instance, narrow-lane fixing rates are less than 40% for all the dual-frequency solutions and increase to more than 80% after the PPP-AR step. This is proven to improve satellite orbit precision by about 20%. The use of multi-frequency observations in the 15-station network  The precision of clock offsets and biases depends on the number of tracking stations. A larger network always provides more precise and robust solutions. For instance, the STD of satellite clock differences from the 50-station network solution relative to those of a 90-stations network is 13 ps and increases to 47 ps for the 15-station network solution. Accordingly, phase biases from the 15-station network are not well determined. The ratio of narrow-lane bias differences (compared to the reference solution) within 0.05 narrow-lane cycles is about 50%. By using 2.5 times more observations in the multi-frequency processing, the ratio improves to about 70%. When using 20 (or more) tracking stations, the ratios are higher than 80%.
Kinematic PPP/PPP-AR results indicate that phase bias products from the 20-station network are precise enough to perform single-receiver ambiguity resolution. KPP-AR positioning solutions show an improvement of 30-40% compared to float-ambiguity solutions in the East component. Multi-frequency KPP-AR results show further improvement of 5-10% with respect to dual-frequency KPP-AR solutions.
Author Contributions BD and UH proposed the general idea of this contribution. BD adopted the software, did the computation and wrote the draft version. OM analyzed the signal and measurement performance and wrote Sect. 3. PS computed Galileo satellite PCOs of different frequencies and wrote the corresponding part in Sect. 4. UH and OM commented regularly on the results and gave suggestions. UH, OM and PS reviewed and revised the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.
Data availability GNSS and SLR data used are publicly available from International GNSS Service and International Satellite Laser Ranging Service of respective data centers.

Conflict of interest
The authors declare that they have no conflict of interest.
Code availability Calculations are based on the Bernese GNSS Software (license available, University of Bern) and its specific modifications by the authors (not public). PCO estimates were computed with the NAPEOS Software (license available, ESA).
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://creativecomm ons.org/licenses/by/4.0/.