Research ArticleOriginal Article
Open Access

A multi-frequency method to improve the long-term estimation of GNSS clock corrections and phase biases

A. Rovira-Garcia, J. M. Juan, J. Sanz, G. González-Casado, J. Ventura-Traveset, L. Cacciapuoti and E. Schoenemann
NAVIGATION: Journal of the Institute of Navigation December 2021, 68 (4) 815-828; DOI: https://doi.org/10.1002/navi.453

Abstract

The space segment of the Global Navigation Satellite System (GNSS) is equipped with highly stable atomic clocks. In order to use these clocks as references, their time offsets must be estimated from ground measurements as accurately as possible. This work presents a multi-frequency and multi-constellation method for estimating satellite and receiver clock corrections, starting from unambiguous, uncombined, and undifferenced carrier-phase measurements. A byproduct of the estimation process is phase biases (i.e., the hardware delays of the carrier-phase measurements occurring at receivers and satellites). The stability and predictability of our clock estimates for receivers and satellites (GPS and Galileo) are compared with those obtained by the International GNSS Service (IGS), whereas the phase biases are assessed against two independent determinations involving combinations of carrier-phase measurements. We conclude that the method reduces day boundary discontinuities in the clock corrections, and that the estimated phase biases reproduce variabilities already observed by other authors.

KEYWORDS

1 INTRODUCTION

Time and frequency transfer is one of the main uses of the Global Navigation Satellite System (GNSS), besides geodetic or navigation applications. Stable ground clocks equipped with a passive hydrogen maser (PHM) or rubidium atomic frequency standard (RAFS) can be used as global time references. In such a way, the clocks associated with other ground receivers or satellites can be estimated and referred to any of these reference PHMs.

The International GNSS Service (IGS) has been routinely producing precise orbit and clock corrections for geodetic and timing applications in daily batches since 1998 (Beutler et al., 1999; Johnston et al., 2017). These estimations are performed simultaneously with other parameters such as orbit parameters, Earth rotation parameters, station coordinates, and tropospheric delays or hardware biases among others.

Due to the large number of parameters to be determined and their different stability on time, the estimation is usually done in a two-step process (Bock et al., 2009). In the first step, all parameters are computed with a low sampling rate (e.g., five minutes). Then, in a second step starting from the solutions for the low-varying parameters, satellite and receiver clock corrections are estimated with a higher sampling rate (e.g., 30 seconds, even 1 second), as seen in Bock et al. (2009) and Rovira-Garcia et al. (2016).

In such an estimation procedure, GNSS measurements are collected by a worldwide network of receivers located at permanent stations with coordinates accurate to the centimeter. This allows subtraction from the measurements’ well-known geodetic models, such as those used in the precise point positioning (PPP) technique (Malys & Jensen, 1990; Zumberge et al., 1997). In this manner, PPP accurately corrects antenna offsets, solid and ocean tide displacements, relativistic and gravitational effects, tropospheric propagation delays, and windup effect among other terms.

However, the ionospheric delay experienced by measurements cannot be modeled with enough accuracy. To eliminate first-order ionospheric delay, the standard procedure for estimating satellite and receiver clock corrections uses the dual-frequency ionospheric free (IF) combination of code pseudoranges 𝑃𝐼𝐹 and carrier-phase 𝐿𝐼𝐹 measurements (Subirana et al., 2013).

The obtained residuals Δ𝑃𝐼𝐹 and Δ𝐿𝐼𝐹 (i.e., measurements corrected from the PPP models) for any receiver-satellite pair can be used to estimate the satellite and receiver clock corrections as white noise processes and the receiver hardware biases between different constellations as a random walk process. As commented before, this usually occurs in the second step at a high throughput, starting from the orbits and clock corrections previously estimated in the first step.

Alternative to the use of IF combinations for estimating clock corrections, some authors like Schöenemann et al. (2011), Odijk et al. (2016), and Strasser et al. (2019) have obtained similar results using uncombined measurements, the so-called raw observation processing, in which ionospheric delays are estimated as white noise processes or through the estimation of additional parameters of an ionospheric model.

Since the early days of GPS, the importance of integer ambiguity resolution (IAR) has been recognized to enhance the accuracy of geodesy applications including clock determination (Blewitt, 1993). IAR reduces the number of parameters to be estimated and avoids the correlations between ambiguities and clocks. Thus, IAR is a standard procedure performed by different analysis centers (ACs) contributing to the IGS. In particular, the Double-Differenced (DD) method (Beutler et al., 1986) is applied to compute daily batches of precise orbits and clock corrections. Moreover, different authors showed the advantages of performing undifferenced IAR in PPP (Ge et al., 2008; Laurichese & Mercier, 2007).

However, as the processing is done in daily batches, IAR does not guarantee that the same integer values for ambiguities can be obtained every day. That is, in the daily estimates, arcs of continuous carrier-phase measurements spanning two consecutive days can present different values each day on their ambiguities and, consequently, in the clock estimates.

This problem is known as day boundary discontinuity (DBD), such as seen for instance in Ray and Senior (2003). If continuity is not ensured after performing the daily IAR process, DBDs in clock estimates can reach the decimeter level (i.e., several tenths of a ns).

The reason behind DBDs is that the noise and the instrumental delays present in the code pseudoranges can vary differently from those in the carrier-phase measurements. A solution to the DBDs proposed in Collins et al. (2010) was to define a different clock for pseudorange and carrier-phases (i.e., the decoupled clocks technique). DBDs on satellite clock solutions introduce artificial discontinuities that limit the accuracy of clock analysis (e.g., frequency stability) for periods longer than one day (Collins et al., 2010; Dach et al., 2006; Petit et al., 2015).

Recently, Rovira-Garcia et al. (2021) presented a methodology using only unambiguous carrier-phase measurements [as in Collins et al. (2010)]. The method guarantees the continuity of the unambiguous carrier-phase measurements over days, even when IAR is performed in daily batches and thus is able to reduce the DBDs between the clock estimates of such adjacent daily batches.

The aim of the present work is, on one hand, to extend the study of Rovira-Garcia et al. (2021) to hundreds of days and, on the other hand, to show that the methodology can be also applied when more than two frequencies are involved. Notice that handling more than two frequencies implies addressing different temporal evolution of the phase biases at different frequencies, as pointed out in Montenbruck et al. (2012) for the GPS satellites of Block IIF.

2 DATA SET

We have applied our methodology to more than 340 days in the year 2017, extending the proof of concept of 16 days of data presented in Juan et al. (2020) and Rovira-Garcia et al. (2021). In order to estimate the satellite and receiver clock corrections, together with the phase biases, we downloaded RINEX observation files containing code pseudoranges and carrier-phase measurements every 30 s.

The frequencies used by this method are L1, L2, and L5 for GPS, and E1 (L1), E5a (L5), and E5b (L7) for Galileo. This choice applies to each of the approximately 150 receivers located at permanent stations within the IGS network that were used in the study. The redundant number of receivers allows robust estimation of the phase biases, despite the fact that some regions (e.g., oceans) are not as well covered as continental areas.

In order to compare our estimates, we used RINEX clock files with precise clock solutions computed by different IGS ACs contributing to the Multi-GNSS Experiment (MGEX) (Montenbruck et al., 2017). Specifically, the selected clock products for the study are: the GRM product of the Centre National d’Études Spatiales (CNES), the GBM product of the GeoForschungsZentrum (GFZ), the COM product from the Center for Orbit Determination in Europe (CODE), and the repro4 product from the European Space Operations Center (ESOC). The first three products are publicly available at the IGS MGEX product archive, whereas repro4 is an adhoc process of ESOC in the context of the projects named the Gravitational Redshift Experiment with Eccentric Satellites (GREAT) funded by ESA (Delva et al., 2018; Herrmann et al., 2018; Rovira-Garcia et al., 2021). All ACs provided satellite clock corrections at an interval of 30 s, although COM initially contained estimates every 300 s. Eventually, after day 226 in 2017, COM provided determinations every 30 s as well.

We focused our study on GPS and Galileo clock estimates, which were available from all ACs except ESOC, which only provided estimates for Galileo satellites in its repro4. In this ESOC solution, the reference receiver for every day during the entire 2017 was the station BRUX, equipped with a very stable PHM atomic clock.

In contrast to ESOC, the other ACs typically re-aligned the clock estimates to the broadcast GPS time for the satellites as an underlying timescale (Ray et al., 2017). That is, first the ACs estimated all clock corrections by setting the value of one reference station (that could vary from day to day) to zero, and then, in a second step, the ACs subtracted from the clock estimates the mean value of all satellite clock determinations.

This approach is robust to variations of the reference station from one day to another, thus should produce continuous satellite clocks. However, when the satellite set changes, the mean value of the satellite clocks changes accordingly, in a similar way as the jumps in the differential code biases (DCBs) reported in Sanz et al. (2017). This issue, among others, can produce DBDs as shown in the results section.

The present study focuses on the aforementioned second step, in which satellite and receiver clock corrections are estimated starting from a previous determination of slow varying parameters. In this regard, our starting product is the repro3 solution; the penultimate version of the products computed by ESOC for the GREAT study. The main differences between repro3 and repro4 are that repro3 includes satellite orbits and clock corrections for GPS and Galileo and the reference station changes from day to day. The clock corrections examined in the results section correspond to the repro4 product.

Our study sought to evaluate the nominal performance of the clock corrections determined by our approach and those computed by the aforementioned four IGS ACs. Then, in order to perform a fair comparison, we scanned the products computed during the entirety of 2017. We looked for periods in which none of the five products presented a large discontinuity associated, for instance, to the chosen clock reference or any other potential problem associated with the clock determination methodology.

After the examination, we selected two common periods: days 191 through 230 and days 281 through 310. In these two periods, none of the five solutions presents an anomalous behavior and thus are representative of a nominal performance.

3 METHODOLOGY

We propose a method to estimate clock corrections and its associated inter-frequency biases using only carrier-phase measurements. As in Strasser et al. (2019), our approach uses uncombined carrier-phase measurements (instead of the standard IF combination of carrier phases) in which the integer part of each carrier-phase ambiguity has been fixed through a previous estimate of the phase biases. Moreover, our approach supports more than two frequencies in processing multiple constellations and guarantees the continuity of phase biases over days. In our methodology, the carrier-phase and pseudorange code measurements at a given frequency m between a receiver i and a satellite j can be written as:

Embedded Image 1

Embedded Image 2

where Embedded Image is the geometric distance between the receiver and satellite antennae, c is the speed of light, 𝑇𝑖 and 𝑇𝑗 are the receiver and satellite clock offsets, the factor Embedded Image converts from the total electron content unit (TECU) to meters of delay at the frequency 𝑓𝑚, Embedded Image and Embedded Image account for the slant tropospheric and ionospheric delays between the receiver and the satellite. The 𝐷𝑚𝑖 and 𝐷𝑚𝑗 are the DCB for the receiver and the satellite, whereas 𝛿𝑚𝑖 and 𝛿𝑚𝑗 are the phase instrumental delays, and Embedded Image is the integer carrier-phase ambiguities expressed in cycles and converted to meters by the corresponding wavelength 𝜆𝑚. Finally, 𝜖 and 𝜀 represent noise terms.

The key point of the methodology is to resolve the carrier-phase ambiguities Embedded Image in an undifferenced way. In order to guarantee the continuity of the unambiguous carrier-phase measurements over days, we have to do a rough estimation for the phase biases of both satellites and receivers at all frequencies and constellations. This previous process is described in detail in Rovira-Garcia et al. (2021) and Juan et al. (2020), allowing the estimation of the integer part of the carrier-phase ambiguities using standard PPP models.

Once the integer part of the ambiguities is removed to the raw undifferenced carrier-phase measurements in Equation (1), we then apply again geodetic models in order to obtain the phase residuals Δ𝐿 (i.e., undifferenced and ambiguity-resolved carrier-phase measurements corrected from effects modeled with PPP). Δ𝐿 at a given frequency m between a receiver i and a satellite j can be written as:

Embedded Image 3

where the instrumental delays 𝑘𝑚𝑖 +𝑘𝑚𝑗 have been redefined in Equation (3) as:

Embedded Image 4

Note that Equation (3) is only applied for carrier-phase measurements for which the ambiguities are fixed, otherwise, those measurements are discarded. This corresponds to an IAR success rate around a 90%, similar to Rovira-Garcia et al. (2021) in its 16-day study. The modeling error and other errors such as thermal noise and multipath of Δ𝐿 in Equation (3) are typically below one centimeter.

The method to solve Equation (3) uses the Kalman filter to estimate the clock corrections of satellites and receivers, the slant ionospheric delay between each receiver-satellite pair, and the phase biases of satellites and receivers. The estimates are produced at regular intervals of 30 s. Note that the value of the phase biases is strongly correlated with the ionospheric term in the right-hand side of Equation (3).

The process noise used in the filter to estimate the clock corrections and ionospheric delays is white noise, whereas the phase biases are estimated as a random walk. The reason not to estimate them later as a constant offset is to allow studying the evolution of the phases biases with time that can be related to variations in the satellite temperature or in the orbital position (Montenbruck et al., 2012).

In order to solve the rank deficiency in Equation (3), we have to define one reference frequency for each device (satellite or receiver) and one reference clock for referring the rest of the clock corrections. We have selected the receiver CEBR as the reference clock, as this station is equipped with a very stable PHM atomic clock.

For the satellites, the reference frequency can be imposed by assuming that the phase biases are the same for two of the frequencies involved in the processing. That is, for each GPS or Galileo satellite, we impose the constraint:

Embedded Image 5

From Equation (3), this is equivalent to take the IF combinations built with frequencies 𝑓1 and 𝑓2 for each GPS satellite, and 𝑓1 and 𝑓5 for each Galileo satellite, as the satellite clock reference, respectively.

In the same way, for each receiver, we have to define one reference frequency. The reference frequency can be defined linked to the IF combination of GPS frequencies 𝑓1 and 𝑓2, by assuming that the hardware biases are equal in the constraint:

Embedded Image 6

which is equivalent to consider the IF combination of these two carrier-phases as the reference for the receiver clock corrections. In this sense, when considering a triple-frequency processing per constellation, the method has to estimate for each receiver two different phase biases for GPS (𝑘1𝐺𝑃𝑆,𝑖 and 𝑘5𝐺𝑃𝑆,𝑖) and three phase biases for Galileo (𝑘1𝐺𝐴𝐿,𝑖 and 𝑘5𝐺𝐴𝐿,𝑖 and 𝑘7𝐺𝐴𝐿,𝑖) because, for the latter, we have not imposed any constraint relation between them.

Extending the same reasoning to the satellite phase biases, the method estimates two different phase biases for each GPS satellite (𝑘1𝐺𝑃𝑆,𝑗 and 𝑘5𝐺𝑃𝑆,𝑗) and two different phase biases for each Galileo satellite (𝑘1𝐺𝐴𝐿,𝑗 and 𝑘7𝐺𝐴𝐿,𝑗).

4 RESULTS

This section analyzes the estimates performed by the proposed method. First, we turn our attention to the input data used by the estimation. Second, we analyze the satellite and receiver clock results, specifically their predictability and stability. Finally, we analyze the satellite phase biases observing its temporal evolution.

4.1 Carrier-phase residuals

Figure 1 depicts an example of the carrier-phase residuals computed with Equation (3) for up to three frequencies and different Galileo and GPS satellites using their pseudorandom numbers (PRNs). In order to magnify the continuity of the residuals at the day change, we have applied a linear detrending. That is, for each receiver-satellite pair, we have subtracted from the residuals the value of a linear model adjusted by the whole data.

FIGURE 1
FIGURE 1

Linearly detrended carrier-phase residuals of different GPS satellites: G01 (Block IIF), G22 (Block IIR-B), and G32 (Block IIF); and Galileo satellites: E09 (GSAT209 on PHM) and E22 (GSAT204 on RAFS) at midnight from day 110 through 111 in 2017, for the f1 (top), f5 (middle), and f2 ∕f7 (bottom)

In this sense, we eliminate large trends (bias and drift) in the receiver and satellite clock corrections 𝑇𝑖 and 𝑇𝑗, but the evolution of all other terms in Equation (3) can still be observed in Figure 1, that is, mainly the clock corrections and ionosphere. The magnitude of the ionosphere variation is at the level of a few tens of decimeters and depends on the frequency as indicated at the top (f1), middle (f5), and bottom (f2) and (f7) panels.

Apart from the different stability of the depicted satellite clock corrections, the detrended residuals confirm that the daily IAR pre-process does not introduce any discontinuity at the midnight epoch. Because these residuals are the input data in Equation (3) for estimating the receiver and satellite clock corrections, and the remaining terms (ionospheric delay and phase biases) should exhibit a continuous behavior, the elimination of DBDs in the input data guarantees the inter-day continuity of the clock solutions.

Notice that in Figure 1 it is also possible to see the different stability of different satellite clocks. Indeed, the noisier clock occurs in the GPS satellite G22 (Block IIR-B), while the clocks of the Block IIF GPS satellites are more stable. It is even possible to distinguish between Galileo satellites which, according to the Notice Advisory to Galileo Users (NAGUS) published by the European GNSS Service Centre (EUSPA, 2017), during these days, E09 was operating with a PHM, whereas E22 had activated its RAFS.

4.2 Satellite clocks estimates

For two Galileo satellites operating with PHMs, Figure 2 illustrates the clock estimation performed by the proposed approach together with the products from the IGS ACs over days 191 through 230 and days 281 through 310. As an example, we selected two Galileo satellites operating with PHMs, but other Galileo satellites have similar performances. As in Figure 1, in order to better depict the continuity of the clock estimates over the day changes, we have detrended the clock solutions by subtracting a linear model computed for each AC and for each period of more than one week. Moreover, in order to enhance the plot visibility, the time series of each AC are shifted with one arbitrary value.

FIGURE 2
FIGURE 2

Satellite clock offsets for Galileo satellites E01 (GSAT210 on PHM) in the top row and E30 (GSAT206 on PHM) in the bottom row, estimated by the proposed approach (black), ESOC (grey), CODE (red), GFZ (blue), and CNES (green); the left column depicts days 191–199 within the first period in 2017, whereas the right column depicts days 281–289 within the second period in 2017

Figure 2 depicts different patterns of the detrended clock solutions linked to mismodeling the different effects in the GNSS signals. For instance, it is well-known that the radial component of the satellite position is quite correlated with the satellite clock corrections. Therefore, as it is stated in Delva et al. (2018) and Herrmann et al. (2018), an accurate estimate of the satellite clock correction requires accurate force models, otherwise, clock estimates would present patterns linked to the satellite orbits.

However, errors in the radial component are not the only source of error affecting satellite clock estimates. Indeed, as it is seen in previous references, other effects related to satellite position, such as geomagnetic field, satellite temperature, or even relativistic correction, can affect the satellite clock solution. In any case, these types of errors related with satellite position are expected to be continuous over time, so they are not the main cause of the DBDs.

Regarding DBDs, it can be seen that the clock corrections estimated by the proposed approach (labeled as gAGE), which relies on a careful IAR pre-process, are continuous regardless of the day change. The remaining discontinuities observed in our clock corrections (e.g., in day 284 for E01 or day 192 for E30) are caused by the starting products (precise orbit and clock corrections) used to model the carrier-phase measurements in Equation (3). In contrast, we can observe DBDs in the clock estimates of other ACs, which reach the decimeter level (i.e., the ns level).

A standard metric to evaluate the stability of clock corrections is through the Allan deviation (ADEV) as defined in Allan (1966). Figure 3 illustrates the ADEVs for time intervals up to 533,000 seconds, which correspond to 6.17 days, for the two previous Galileo satellites E01 and E30. Note that the solution computed by CODE is excluded from the first period, as the clock estimates were produced at a sampling interval of 300 seconds and therefore does not reflect a fair comparison with clock solutions every 30 seconds.

FIGURE 3
FIGURE 3

Allan deviations for Galileo satellites E01 (GSAT210 on PHM) in the top row and E30 (GSAT206 on PHM) in the bottom row, estimated by the proposed approach (black), ESOC (grey), CODE (red), GFZ (blue), and CNES (green); the left column corresponds to the entire first period analyzed in 2017, whereas the right column depicts the second period

Moreover, we have excluded days 293, 301, and 310 of the CNES solution due to an anomalous behavior affecting all the satellites simultaneously. We can see that our approach produces clock estimates that are similarly as stable as those computed by other IGS ACs. In this sense, the ADEV does not evidence the advantage of the DBD removal observed in Figure 2. This is because the DBD magnitude is approximately a few decimeters and only occurs once per day (every midnight in GPS time), whereas the ADEV computation includes many other points that belong to the same day (i.e., unaffected by DBDs). In some solutions, the fluctuations of the clock estimates are comparable to the DBDs. These two facts dilute the effect of the presence of DBDs in the ADEVs depicted in Figure 3.

In order to emphasize the advantage of the proposed method in reducing DBDs, we used one hour of clock solutions during the last hour of one day to fit a linear model using least squares. Then, this linear model was propagated one hour forward into the following day. Finally, we computed the differences between the linear model predictions and the actual clock estimates over the first hour of the following day. This procedure assessed the extrapolation errors, which is of primary interest for navigation applications, specifically, for high-accuracy service providers questing for an optimal refresh rate of satellite clock corrections.

Figure 4 illustrates the extrapolation errors, aggregating all GPS and all Galileo satellites during each period considered in 2017. The histogram of prediction errors obtained from the clock estimates determined by the proposed approach reached higher peaks centered around zero than those obtained by other ACs, evidencing the benefits of the DBD removal.

FIGURE 4
FIGURE 4

Histogram depicting the prediction errors after linear extrapolation for all GPS and Galileo satellites in the top and bottom panels, respectively; the left column depicts the first period considered in 2017, whereas the right column depicts the second period; the bin size of the histograms was 5 mm

The peaks of our approach for Galileo are higher than those of GPS, confirming the better predictability (i.e., stability) of Galileo clocks. The satellite clock estimates computed with three frequencies are similar to those previously estimated with two frequencies in Rovira-Garcia et al. (2021). The necessary estimation of the additional phase bias mitigated the redundancy of having an additional frequency.

The effect of the lower stability of GPS clocks can also be seen in the histograms of GPS clock estimates of other ACs: they look more similar between different ACs than those of Galileo, suggesting that the lower predictability of GPS clocks dominates the differences in estimation methodologies.

We can see that ESOC and GFZ improve the Galileo determinations from the first to the second period, but the reason is unclear. In contrast, our estimates do not substantially vary from the first to the second period maintaining a smaller distribution of prediction errors in comparison to other ACs, whose distributions are not centered in zero because the DBDs are not necessarily symmetrical.

4.3 Receiver clock estimates

Figure 5 depicts the results obtained for receiver clock corrections during eight days in each of the aforementioned periods. In order to better illustrate the continuity of the clock estimates over the midnight epoch, we selected two IGS stations, BRUX and MGUE. These two stations are equipped with stable PHM clocks and hence a high predictability can be expected. Moreover, all the ACs provide clock solutions for these stations during the selected periods.

FIGURE 5
FIGURE 5

Receiver clock offsets for permanent stations equipped with a PHM clock: BRUX (top row) and MGUE (bottom row), estimated by the proposed approach (black), ESOC (grey), CODE (red), GFZ (blue), and CNES (green); the left column depicts days 191–199 within the first period in 2017, whereas the right column depicts days 281–289 within the second period in 2017

As in the case of the satellites in Figure 2, we have detrended each time series of the clock estimates with a linear model. We can observe that the receiver clock corrections determined by the proposed methodology present a high continuity, because of the mitigation of the DBDs in input carrier-phase measurements.

The most stable solution is obtained by ESOC in the station BRUX, but as noted Section 2, BRUX has been chosen as reference in the daily solutions of ESOC. Apart from this case, the receiver clock estimates using the proposed approach seem to be more stable than the other clock solutions.

We can observe that our solution is slightly noisier than the one obtained by the ACs. The reasons that can explain this fact are: first, our receiver estimates were given every 30 s, whereas the ACs gave their receiver clock estimates every 300 s; and second, our approach used orbits and clock corrections externally computed by ESOC as starting products. Therefore, any difference on the models for computing these products with respect to our model could impact our clock solutions.

The solution of CODE presents a nearly pure linear drift for MGUE during the days 195, 197, 282, 284, and 287, and for BRUX during the day 285, which corresponds to these stations being selected as reference stations in their processing.

Figure 6 depicts the ADEV analysis for the two receiver clocks previously shown in Figure 5. In order to compute meaningful ADEVs, we selected days within the first and second period in which the receivers did not experience any reset. The clock results for the stations present lower ADEVs than those obtained for the satellites in Figure 3.

FIGURE 6
FIGURE 6

Allan deviations for permanent stations BRUX (top row) and MGUE (bottom row), estimated by the proposed approach (black), ESOC (grey), CODE (red), GFZ (blue), and CNES (green); every column represents a period in which the receivers did not experience any reset within the first period (left column) and second period (right column) studied in 2017

The highest stability (i.e., lowest values of ADEVs) was obtained by ESOC in station BRUX, but it is recalled that it corresponds to the reference station. Moreover, during certain days, other solutions used BRUX or MGUE as a reference station. We can observe that our clock estimates reach 10−15 in both stations, for long averaging times. However, as in the case of satellites, our ADEVs do not evidence the advantage of the DBD removal observed in Figure 5. This is a consequence of the large number of points that intervened in every plot, whereas decimeter-level DBDs only occur once per day.

In order to quantify the improvement of the receiver clock estimates using the new approach, we have applied the same procedure as the satellites to evaluate the DBD mitigation. That is, to compute the residuals of the first hour of the day with respect to a linear prediction fitted with the clock estimates of the last hour of the preceding day.

To enhance the visibility of the comparison, we have calculated the root mean square (RMS) of the residuals during the first hour of each day. We use those RMS values to evaluate the error of the linear extrapolation during the entirety of the two periods.

Figure 7 shows the cumulative distribution function (CDF) of the daily RMS of the residuals for the two stations previously depicted in Figure 5. We can observe that the prediction error of gAGE clock estimates is the lowest, with values at the level of 3 cm (approximately 100 ps).

FIGURE 7
FIGURE 7

CDF of the daily RMS values of the prediction errors after linear extrapolation for permanent stations equipped with a PHM clock: BRUX (top row) and MGUE (bottom row), estimated by the proposed approach (black), ESOC (grey), CODE (red), GFZ (blue), and CNES (green); the left column corresponds the first study period in 2017, whereas the right column depicts the second period

As in the case of satellites, these errors are not far from the nominal accuracy of the IGS final products (approximately 75 ps) and are lower than the errors reported in Ray and Senior (2003) with just GPS satellites.

In contrast, we can see that the clock estimates of the other ACs presented median errors at the level of one decimeter, except ESOC estimates for BRUX station, which was at the same level or better than the clocks determined by our approach. However, BRUX is the reference clock in its daily solutions.

4.4 Satellite phase biases

As commented in Section 3, the method estimates the satellite and receiver phase biases simultaneously with clock offsets and ionospheric delays. Because of the one-to-one correlation between the biases and the slant ionospheric delay, their absolute value is not well determined. However, the difference of phase biases (e.g., 𝑘2 −𝑘1, 𝑘5 −𝑘1, or 𝑘7 −𝑘5) can be accurately determined thanks to the elimination of the ionosphere. In what follows, we evaluate our phase biases with respect to two independent determinations.

First, we use the combination ionosphere- and geometry-free (IGF), 𝐿𝐼𝐺𝐹, which can be built as the difference of two IF combinations as seen in Li et al. (2020). For instance, using the GPS measurements and frequencies 𝑓1, 𝑓2, and 𝑓5, the IGF for a receiver i, and a satellite j is built as:

Embedded Image 7

where Δ𝑘 = 𝑘5 −𝑘1 stands for the phase bias difference between frequencies 𝑓5 and 𝑓1. Recall that because of the constraints imposed by Equations (5) and (6), the biases for 𝑓1 and 𝑓2 cancel for both receiver and satellites in the IF combination of reference. We can proceed in the same manner for Galileo satellites to obtain the differences of phase biases, for instance, for 𝑘7 −𝑘1.

Because the unambiguous 𝐿𝐼𝐺𝐹 values are accurately known, Equation (7) allows estimating the phase bias differences of satellites and receivers. In order to solve the rank deficiency of the problem, we have to constrain to zero one of the Δ𝑘 values in Equation (7), such as the value for the reference receiver CEBR.

In this way, the rest of Δ𝑘 can be determined in a similar way as the DCB estimation through an ionospheric model [see Sanz et al. (2017)]. During some epochs, it can sometimes occur that there are no observations in common throughout the satellites and stations with the reference. This can generate discontinuities in the clock correction estimates, and in some cases, can impede the estimation. An example of this can be seen in the bottom-left plot of Figure 8, where the red crosses toward the end of day 284.

FIGURE 8
FIGURE 8

Time series depicting differences of phase biases for GPS satellites of the Block IIF G08, G09, and G10 in the top-left, top-right, and bottom-left panel, respectively, whereas Galileo satellite E30 (GSAT206 on PHM) is depicted in the bottom-right panel

A second independent estimation of the phase biases can be performed using the Hatch-Melbourne-Wübenna (HMW) combination (Hatch, 1982). That is, the difference of the wide-lane (WL) combination of carrier-phases and the narrow-lane (NL) combination of code pseudoranges.

The HMW cancels out the terms related to 𝛼𝑚 in Equation (1) and Equation (2) and therefore, once having the ambiguities resolved, it can be written as:

Embedded Image 8

where 𝑓𝑚 and 𝑓𝑛 correspond to the frequencies 𝑓1 and 𝑓5 for GPS and 𝑓5 and 𝑓7 for Galileo. Note the participation of the code pseudorange in Equation (8), unlike Equations (3) and (7) that are based on phase-only measurements.

These comparisons provide us an additional assessment of the suitability of our methodology to handle more than two frequencies. Figure 8 depicts the differences of the estimated phase biases from Equation (3), together with the two independent determinations computed from Equations (7) and (8). In order to magnify this variability, the comparison is done for one Galileo satellite and three GPS satellites of Block IIF, which depict larger fluctuations in comparison with the biases of the Galileo satellite.

We can observe the agreement in the patterns of the differences of estimated phase biases, with those computed directly from measurements, for both the IGF combination and the HMW combination, the latter only biased by the DCB value as expressed by Equation (4). This agreement confirms that the estimation process is fully compatible for three frequencies.

Moreover, it also shows that estimating the hardware delays with a random walk, instead of a bias, is an adequate choice to accommodate periodic variations of Block IIF satellites. Specifically, the time variations of the biases in the third GPS frequency, 𝑓5, which were already pointed out in Montenbruck et al. (2012) to be linked to temperature changes during the eclipsing periods.

It is apparent that the Galileo satellite presents a higher stability with respect to time than the GPS satellites. This is thanks to the fact that Galileo satellites maintain the temperature of the onboard atomic clocks, leaving them highly stable. The drawback of estimating hardware delays such as random walk, instead of biases, is that the use of three frequencies does not strengthen the clock estimates with respect to the estimates using just two frequencies.

5 CONCLUSIONS

The present contribution shows that the methodology defined in Juan et al. (2020) and Rovira-Garcia et al. (2021) to estimate receiver and satellite clock corrections can be applied using three frequencies. The method requires unambiguous carrier-phases as input measurements. Avoiding pseudorange measurements in the estimation process largely mitigates discontinuities in the clock estimates in the day change.

Unlike other clock solutions, DBDs in our clock estimates are reduced to the level of few centimeters (approximately 100 ps), which is at the same level as the IGS final products (75 ps). This significant DBD reduction of our proposed methodology may benefit the statistical characterization of long-term phenomena correlated with the onboard satellite clocks. The stability of our satellite clock estimates has been measured with the ADEV computed over periods of tens of days, presenting similar stabilities as the clock solutions of other ACs, despite the latter containing DBDs of a few decimeters.

The approach considered here uses several carrier-phase measurements without the need of building IF combinations, which is the standard way for estimating clock corrections. Because the estimation process uses three frequencies, clocks offsets are estimated jointly with instrumental phase biases. These hardware delays are estimated as random walk processes, instead of computing them as constant parameters.

This treatment of the process noise allows us to take into account the dependency of the temperature of some of these biases already pointed out by some authors. Differences of our phase bias estimates have been compared with independent determinations built with ionosphere- and geometry-free combinations of GNSS observables. From these comparisons, we had an additional confirmation of the suitability of the methodology.

HOW TO CITE THIS ARTICLE

Rovira-Garcia A, Juan JM, Sanz J, et al. A multi-frequency method to improve the long-term estimation of GNSS clock corrections and phase biases. NAVIGATION. 2021;68:815–828. https://doi.org/10.1002/navi.453

ACKNOWLEDGMENTS

The present work was supported in part by the European Space Agency contract (REL-GAL) N.4000122402/17/NL/IB, by the project RTI2018-094295-B-I00 funded by the MCIN/AEI 10.13039/501100011033 which is co-founded by the FEDER programme, and by the Horizon 2020 Marie Skłodowska-Curie Individual Global Fellowship 797461 NAVSCIN. The authors acknowledge the use of data and products provided by the International GNSS Service.

Footnotes

  • Funding information

    H2020 Marie Skłodowska-Curie Actions, Grant/Award Number: 797461 NAVSCIN; Ministerio de Ciencia, Innovación y Universidades, Grant/Award Number: RTI2018-094295-B-I00; European Space Agency, Grant/Award Number: N.4000122402/17/NL/IB

  • Received March 16, 2021.
  • Revision received July 16, 2021.
  • Accepted September 25, 2021.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

REFERENCES

  1. ↵
    1. Allan, D. W.
    (1966). Statistics of atomic frequency standards. Proceedings of the IEEE, 54(2), 221–230. https://doi.org/10.1109/PROC.1966.4634
    CrossRef
  2. ↵
    1. Beutler, G.,
    2. Gurtner, W.,
    3. Rothacher, M.,
    4. Schildknecht, T., &
    5. Bauersima, I.
    (1986). Determination of GPS orbits using double difference carrier phase observations from regional networks. Bulletin Géodésique, 60(3), 205–220. https://doi.org/10.1007/BF02521052
  3. ↵
    1. Beutler, G.,
    2. Rothacher, M.,
    3. Schaer, S.,
    4. Springer, T.,
    5. Kouba, J., &
    6. Neilan, R.
    (1999). The International GPS Service (IGS): An interdisciplinary service in support of Earth sciences. Advances in Space Research, 23(4), 631–653. https://doi.org/10.1016/S0273-1177(99)00160-X
    CrossRef
  4. ↵
    1. Bock, H.,
    2. Dach, R.,
    3. Jäggi, A., &
    4. Beutler, G.
    (2009). High-rate GPS clock corrections from CODE: support of 1 Hz applications. Journal of Geodesy, 83(11), 1083–1094. https://doi.org/10.1007/s00190-009-0326-1
  5. ↵
    1. Blewitt, G.
    (1993). Advances in global positioning system technology for geodynamics investigations: 1978–1992. In D. E. Smith & D. L. Turcotte (Eds.), Contributions of space geodesy to geodynamics: Technology, Volume 25 (pp. 195–213).
  6. ↵
    1. Collins, P.,
    2. Bisnath, S.,
    3. Lahaye, F., &
    4. Hérroux, P.
    (2010). Undifferenced GPS ambiguity resolution using the decoupled clock model and ambiguity datum fixing. NAVIGATION, 57(2), 123–135. https://www.doi.org/10.1002/j.2161-4296.2010.tb01772.x
  7. ↵
    1. Dach, R.,
    2. Schildknecht, T.,
    3. Hugentobler, U.,
    4. Bernier, L.-G., &
    5. Dudle, G.
    (2006). Continuous geodetic time-transfer analysis methods. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 53(7), 1250–1259. https://doi.org/10.1109/tuffc.2006.1665073
    PubMed
  8. ↵
    1. Delva, P.,
    2. Puchades, N.,
    3. Schönemann, E.,
    4. Dilssner, F.,
    5. Courde, C.,
    6. Bertone, S.,
    7. Gonzalez, F.,
    8. Hees, A.,
    9. Poncin-Lafitte, C.,
    10. Meynadier, F.,
    11. Prieto-Cerdeira, R.,
    12. Sohet, B.,
    13. Ventura-Traveset, J., &
    14. Wolf, P.
    (2018). Gravitational redshift test using eccentric galileo satellites, Physical Review Letters, 121(23). https://doi.org/10.1103/PhysRevLett.121.231101
  9. ↵
    1. European Union for the Space Programme Agency
    (2017). Notice Advisory to Galileo Users (NAGU) 2017012. https://www.gsc-europa.eu/notice-advisory-to-galileo-users-nagu-2017012
  10. ↵
    1. Ge, M.,
    2. Gendt, G.,
    3. Rothacher, M.,
    4. Shi, C. &
    5. Liu, J.
    (2008). Resolution of GPS carrier-phase ambiguities in precise point positioning (PPP) with daily observations. Journal of Geodesy, 82(7), 389–399. https://doi.org/10.1007/s00190-007-0187-4
  11. ↵
    1. Hatch, R.
    (1982). The synergism of GPS code and carrier measurements. Proc. of the 3rd International Geodetic Symposium on Satellite Doppler Positioning, Las Cruces, NM, 1213–1231.
  12. ↵
    1. Herrmann, S.,
    2. Finke, F.,
    3. Lülf, M.,
    4. Kichakova, O.,
    5. Puetzfeld, D.,
    6. Knickmann, D.,
    7. List, M.,
    8. Rievers, B.,
    9. Giorgi, G.,
    10. Günther, C.,
    11. Dittus, H.,
    12. Prieto-Cerdeira, R.,
    13. Dilssner, F.,
    14. Gonzalez, F.,
    15. Schönemann, E.,
    16. Ventura-Traveset, J., &
    17. Lämmerzahl, C.
    (2018). Test of the gravitational redshift with Galileo satellites in an eccentric orbit. Physical Review Letters, 121(23). https://doi.org/10.1103/PhysRevLett.121.231102
  13. ↵
    1. Johnston, G.,
    2. Riddell, A., &
    3. Hausler, G.
    (2017). The International GNSS Service. In P. Teunissen & O. Montenbruck (Eds.), Springer handbook of global navigation satellite systems (pp. 967–982). Springer International Publishing. https://doi.org/10.1007/978-3-319-42928-1_33
  14. ↵
    1. Juan, J. M.,
    2. Sanz, J.,
    3. Rovira-Garcia, A.,
    4. González-Casado, G.,
    5. Ventura-Traveset, J.,
    6. Cacciapuoti, L. &
    7. Schoenemann, E.
    (2020). A new approach to improve satellite clock estimates, removing the Inter-day jumps. Proc. of the 51st Annual Precise Time and Time Interval Systems and Applications Meeting, San Diego, CA, 279–301. https://doi.org/10.33012/2020.17306
  15. ↵
    1. Laurichesse, F. &
    2. Mercier, D.
    (2007) Integer ambiguity resolution on undifferenced GPS phase measurements and its application to PPP. Proc. of the 20th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2007), Fort Worth, TX, 839–848. https://www.ion.org/publications/abstract.cfm?articleID=7584
  16. ↵
    1. Li, P.,
    2. Jiang, X.,
    3. Zhang, X.,
    4. Ge, M., &
    5. Schuh, H.
    (2020) GPS+Galileo+BeiDou precise point positioning with triple-frequency ambiguity resolution. GPS Solutions, 24. https://doi.org/10.1007/s10291-020-00992-1
  17. ↵
    1. Malys, S., &
    2. Jensen, P. A.
    (1990). Geodetic point positioning with GPS carrier beat phase data from the CAS UNO Experiment. Geophysical Research Letters, 17(5), 651–654, https://doi.org/10.1029/GL017i005p00651
  18. ↵
    1. Montenbruck, O.,
    2. Hugentobler, U.,
    3. Dach, R.,
    4. Steigenberger, P., &
    5. Hauschild, A.
    (2012). Apparent clock variations of the Block IIF-1 (SVN62) GPS satellite. GPS Solutions, 16(3), 303–313. https://doi.org/10.1007/s10291-011-0232-x
  19. ↵
    1. Montenbruck, O.,
    2. Steigenberger, P.,
    3. Prange, L.,
    4. Deng, Z.,
    5. Zhao, Q.,
    6. Perosanz, F.,
    7. Romero, I.,
    8. Noll, C.,
    9. Stürze, A.,
    10. Weber, G.,
    11. Schmid, R.,
    12. MacLeod, K., &
    13. Schaer, S.
    (2017). The Multi-GNSS Experiment (MGEX) of the International GNSS service (IGS) – Achievements, prospects and challenges. Advances in Space Research, 59(7), 1671–1697. https://doi.org/10.1016/j.asr.2017.01.011
  20. ↵
    1. Odijk, D.,
    2. Zhang, B.,
    3. Khodabandeh, A.,
    4. Odolinski, R., &
    5. Teunissen, P. J. G.
    (2016). On the estimability of parameters in undifferenced, uncombined GNSS network and PPP-RTK user models by means of S-system theory. Journal of Geodesy, 90(1), 15–44. https://doi.org/10.1007/s00190-015-0854-9
  21. ↵
    1. Petit, G.,
    2. Kanj, A.,
    3. Loyer, S.,
    4. Delporte, J.,
    5. Mercier, F., &
    6. Perosanz, F.
    (2015). 1 × 10−16 frequency transfer by GPS PPP with integer ambiguity resolution. Metrologia, 52(2), 301–309. https://doi.org/10.1088/0026-1394/52/2/301
  22. ↵
    1. Ray, J. &
    2. Senior, K.
    (2003). IGS/BIPM pilot project: GPS carrier phase for time/frequency transfer and timescale formation. Metrologia, 40(3), S270–S288. https://doi.org/10.1088/0026-1394/40/3/307
  23. ↵
    1. Ray, J.,
    2. Gurtner, W. &
    3. Coleman M. J.
    (2017). RINEX extensions to handle clock information. https://files.igs.org/pub/data/format/rinex_clock304.txt
  24. ↵
    1. Rovira-Garcia, A.,
    2. Juan, J. M.,
    3. Sanz, J.,
    4. González-Casado, G., &
    5. Bertran, E.
    (2016). Fast precise point positioning: A system to provide corrections for single and multi-frequency navigation. NAVIGATION, 63(3), 231–247. https://doi.org/10.1002/navi.148
  25. ↵
    1. Rovira-Garcia, A.,
    2. Juan, J. M.,
    3. Sanz, J.,
    4. González-Casado, G.,
    5. Ventura-Traveset, J.,
    6. Cacciapuoti, L., &
    7. Schoenemann, E.
    (2021). Removing day-boundary discontinuities on GNSS clock estimates: methodology and results. GPS Solutions, 25(2). https://doi.org/10.1007/s10291-021-01085-3
  26. ↵
    1. Sanz, J.,
    2. Juan, J. M.,
    3. Rovira-Garcia, A., &
    4. González-Casado, G.
    (2017). GPS differential code biases determination: methodology and analysis. GPS Solutions, 21(4), 1549–1561. https://doi.org/10.1007/s10291-017-0634-5
  27. ↵
    1. Schöenemann, E.,
    2. Becker, M., &
    3. Springer, T.
    (2011). A new approach for GNSS analysis in a multi-GNSS and multi-signal environment. Journal of Geodetic Science, 1(3), 204–214. https://doi.org/10.2478/v10156-010-0023-2
  28. ↵
    1. Strasser, S.,
    2. Mayer-Gürr, T., &
    3. Zehentner, N.
    (2019). Processing of GNSS constellations and ground station networks using the raw observation approach. Journal of Geodesy, 93(7), 1045–1057. https://doi.org/10.1007/s00190-018-1223-2
  29. ↵
    1. Subirana, J. S.,
    2. Zornoza, J. M. &
    3. Hernández-Pajarez M.
    (2013). GNSS Data Processing, Vol. I: Fundamentals and Algorithms. European Space Agency Communications. https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf
  30. ↵
    1. Zumberge, J. F.,
    2. Heflin, M. B.,
    3. Jefferson, D. C.,
    4. Watkins, M. M., &
    5. Webb, F. H.
    (1997). Precise point positioning for the efficient and robust analysis of GPS data from large networks. Journal of Geophysical Research: Solid Earth, 102(B3), 5005–5017. https://doi.org/10.1029/96JB03860
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools
A multi-frequency method to improve the long-term estimation of GNSS clock corrections and phase biases
A. Rovira-Garcia, J. M. Juan, J. Sanz, G. González-Casado, J. Ventura-Traveset, L. Cacciapuoti, E. Schoenemann
NAVIGATION: Journal of the Institute of Navigation Dec 2021, 68 (4) 815-828; DOI: 10.1002/navi.453
Twitter logo Facebook logo Mendeley logo
Bookmark this article

Jump to section

Keywords