## Abstract

Onboard time synchronization is an important requirement for a wide range of low Earth orbit (LEO) missions such as altimetry or communication services, and extends to future position, navigation, and timing (PNT) services in LEO. For GNSS-based time synchronization, continuous knowledge about the satellite’s position is required and, eventually, the quality of the position solution defines the timing precision attainable through GNSS measurements. Previous research has shown that real-time GNSS orbit determination of LEO satellites can achieve decimeter-level accuracy.

This paper characterizes the performance of GNSS-based real-time clock synchronization in LEO using the satellite Sentinel-6A as a real-world case study. The satellite’s ultra-stable oscillator (USO) and triple-frequency GPS/Galileo receiver provide measurements for a navigation filter representative of real-time onboard processing. Continuous evaluation of actual flight data over 14 days shows that a 3D orbit root-mean-square (RMS) error of 11 cm and a 0.9-ns clock standard deviation can be achieved.

## 1 INTRODUCTION

Next to positioning and navigation, global navigation satellite systems (GNSSs) are widely used for precise timing. In terrestrial applications, time synchronization is crucial for a wide range of critical infrastructure and services, including communication networks, power grids, and financial transactions. Here, GNSS offers an affordable way to achieve time synchronization with respect to Universal Time Coordinated (UTC) with accuracies up to the ten-nanosecond level.

Timing is also an important requirement in many satellite missions, like the time-tagging of altimetry measurements, payload synchronization of satellites flying in a formation such as bistatic synthetic-aperture radar satellites, and synchronization of constellation-wide networks. The majority of today’s active LEO satellites operate in the areas of remote sensing, telecommunications, and broadband services, where precise time information is mandatory. Furthermore, the feasibility of position, navigation, and timing (PNT) services from LEO constellations has been highlighted in recent years. Advantages and disadvantages of LEOs compared to conventional navigation from MEOs, such as a potentially higher signal strength or an increased number of satellites needed for global coverage, have been widely addressed in scientific publications like Reid et al. (2018). With the announcement and installation of multiple LEO mega constellations like Starlink or OneWeb, the number of LEO satellites is increasing faster than ever before. Such constellations provide new opportunities for transmitting navigation signals synchronized to a GNSS for augmentation purposes.

Other than in terrestrial timing applications, which benefit from the static position of the receiving GNSS antenna, spaceborne time synchronization is intimately tied to the problem of real-time orbit determination. Knowledge of the instantaneous position, particularly its radial component, directly affects the estimated receiver clock offset and, thus, determines the achievable time synchronization performance.

GNSS has been well proven and established for orbit determination in LEO for many years. It has been shown previously that dual-frequency, multi-constellation GNSS orbit determination in LEO using a carrier-phase approach with a reduced dynamics orbit model and broadcast ephemerides can achieve standard deviations of around one decimeter without external augmentation data (Hauschild & Montenbruck, 2021).

Using a similar algorithm, this study characterizes the achievable accuracy and precision of real-time clock determination on a LEO satellite equipped with a GNSS receiver and an external oscillator. Sentinel-6A, the latest Earth observation satellite of Europe’s Copernicus program, is used as a case study. The satellite provides a suitable architecture with two redundant GPS/Galileo receivers connected to an external oscillator, similar to a potential LEO satellite navigation payload. While the oscillator is running freely, the GNSS measurements are time-tagged with the current oscillator beat count. Furthermore, raw code and numerically-controlled oscillator (NCO) phase data are accessible, which allow us to incorporate the oscillator in the observable calculation. This renders Sentinel-6A an interesting candidate to evaluate the potential performance of orbit and clock determination using an online estimation process, which is comparable to real-time onboard processing.

The paper first provides a general overview of Sentinel-6A, including the GNSS receiver and available data as well as the external oscillator. Afterwards, the methodology for orbit and clock determination is described, along with the calculation of observables from raw correlation data and oscillator beats, the definition and generation of different timescales, and the online processing architecture using an extended Kalman filter (EKF). The section concludes with the procedure to calculate a precise reference orbit and clock product for performance analysis. The results of the simulated onboard estimation are presented and compared against the post-processed reference solution to evaluate the achieved orbit and clock precision and accuracy. Furthermore, different biases affecting the clock solution as well as radiation effects on the oscillator and clock model limitations are discussed.

## 2 SENTINEL-6A

Sentinel-6A was launched on November 21, 2020, in an approximately circular low Earth orbit with a nominal altitude of 1,336 km and an inclination of 66° (see Table 1) with the mission to observe the global mean sea level rising (Donlon et al., 2021). Orbit determination was performed using GNSS as well as the Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS; Auriol & Tourain, 2010; Willis et al., 2010) system, and was verified through satellite laser ranging. The french DORIS system consists of a worldwide network with around 60 ground stations and is used for precise orbit determination (POD) of LEO satellites equipped with a DORIS receiver. In the DORIS system, the ground stations act as transmitters and emit signals on two different frequencies at approximately 400 MHz and 2 GHz to cope with ionospheric effects. Those signals are used by DORIS receivers in LEOs to measure Doppler shifts and ranges, allowing for centimeter-level orbit determination. A key component of the DORIS instrument is an ultra-stable oscillator (USO), which is characterized briefly below.

For GNSS reception, Sentinel-6A carries two redundant precise orbit determination receivers (PODRIX) developed by RUAG Space, which are able to track GPS and Galileo signals on the L1/E1, L2, and L5/E5a bands. Similar to Sentinel-3A and Sentinel-3B, the GNSS receivers are connected to the DORIS USO and, thus, allow for the estimation of the USO’s phase and frequency errors in terms of clock bias and drift using GNSS measurements.

### 2.1 DORIS Ultra-Stable Oscillator

The DORIS USO, an oven-controlled crystal oscillator (OCXO), is a core component of the DORIS receiver system and serves as time and frequency reference for DORIS measurements. Its very high short-term stability provides an approximately constant frequency for Doppler measurements during ground station passes, which can last for roughly 1,000 seconds for the Sentinel-6A orbit altitude. For larger timescales, the oscillator frequency has a very small and approximately linear drift in the order of nanoseconds per second compared to timescales such as UTC. Figure 1 shows the estimated phase offset and fractional frequency offset with respect to GPS time for the DORIS USO over a period of 14 days. The clock phase offset is measured relative to the initial epoch and shows an approximately linear behavior. Over the 14-day interval, the clock offset accumulated to roughly 2 milliseconds, which reflects a mean frequency offset of 1.6 · 10^{−9}. The oscillator frequency was not fully constant, and shows an additional drift of 0.25 · 10^{−9} over the two weeks.

However, radiation impacts the USO’s stability with rapid frequency changes occurring over the South Atlantic Anomaly (SAA; Fiandrini et al., 2004). In this region, Earth’s magnetic field is considerably weaker than elsewhere, allowing the Van-Allan belt to reach lower altitudes and, consequently, increase radiation exposure to LEO satellites. The SAA impact on a DORIS USO was first observed on Jason-1 (Willis et al., 2003), which led to pre-launch irradiation of radiation hardening of the oscillators used in the next-generation DGXX DORIS instruments (Mercier et al., 2010). Previous studies have used GNSS measurements in post-processing to observe the SAA impacts on the Sentinel-3A and Sentinel-3B missions (Jalabert & Mercier, 2018; Štěpánek et al., 2020). The radiation effects on the oscillator become visible when looking at the fractional frequency offset over a single day as shown in the bottom plot of Figure 1. The SAA transitions cause the frequency to suddenly decrease at a high rate before recovering to its nominal drift.

### 2.2 GNSS Data

The PODRIX receivers used on Sentinel-6A support the tracking of GPS and Galileo satellites on the L1/E1, L2, and L5/E5a frequencies. In the case of Galileo, the receiver tracks E1-C with BOC(1,1) modulation and E5a-Q pilot signals, which are available on all Galileo satellites. For GPS on the other hand, the PODRIX on Sentinel-6A is configured to track different signals depending on the satellite block. For Block IIR GPS satellites, the measurements of the legacy signals L1 C/A and L1/L2 P(Y) (semi-codeless) are provided, while L1 C/A and L2-CL measurements are provided for the newer Block IIR-M/IIF and GPS III satellites.

The GNSS data used for the experiment are directly extracted from telemetry records of the PODRIX receiver on Sentinel-6A. For increased flexibility in post-processing, the PODRIX receivers’ output correlation results for each tracked signal at the measurement epoch. This includes the code phase (i.e., the signal transmit time) as well as the local NCO phase and carrier-to-noise density ratio (*C* / *N*_{0}). The availability of raw code and NCO phase values provide the flexibility to form pseudorange and carrier-phase observations relative to arbitrary receiver timescales. Aside from the coarse estimate of GNSS time obtained from the receivers’ internal navigation solution, more precise receiver timescales may thus be generated as part of an orbit and clock determination process.

The PODRIX maintains an internal receiver timescale called Instrument Measurement Time (IMT), represented by a discrete counter value *b*. The IMT counts full cycles of the DORIS USO relative to receiver boot time and is obtained from a hardware counter in the Advanced GPS/GLONASS and Galileo Application-Specific Integrated Circuit (AGGA) of the PODRIX. It refers to a core frequency *f*_{core} = 8/3 · *f*_{osc}, which is derived from the native 10-MHz USO frequency *f*_{osc} by phase-coherent upconversion. For each GNSS measurement, the corresponding IMT count is provided as part of the receiver telemetry. This count can later be associated with the GNSS time of the measurement epoch as derived in the receiver’s navigation solution or a more refined orbit and time determination process. The onboard time synchronization is, thus, accomplished through the relation between individual oscillator clock beats and the GNSS system timescale at individual measurement epochs.

## 3 METHODOLOGY

This section describes the algorithm for online orbit and clock determination using an EKF, including observable generation, applied clock model, and data sources. The different biases that need to be considered are introduced and the applied timescales are defined. Finally, the calculation of the precise reference solution is briefly described.

### 3.1 Observable Calculation

For each measurement epoch, the Sentinel-6A telemetry records provide a GNSS navigation solution calculated by the PODRIX receiver, the current IMT beat count, as well as code and NCO phases for each tracked signal. The navigation solution is only used at the first measurement epoch to initialize the clock model and EKF state vector, while phase and IMT data are used to continuously form the GNSS observables at each epoch.

The pseudorange is obtained from the time of signal transmission, represented by the code phase *ϕ*_{code} in units of time with (Won & Pany, 2017):

1

where *c* is the speed of light in vacuum and mod(*x, y*) represents the modulo operator. The time *t* is an arbitrary timescale realized by the receiver and given in seconds at signal reception.

The carrier phase *ϕ* for each tracked signal is calculated with (Won & Pany, 2017):

2

Here *ϕ*_{NCO} is the phase of an NCO, which is used to compensate the Doppler shift and intermediate frequency *f*_{IF} of the incoming signal after down-conversion, while *τ* denotes the time elapsed since the start of tracking (Petovello & O’Driscoll, 2010). *ϕ*_{NCO} constitutes the native measurement provided by the AGGA chip, from which the common GNSS carrier-phase observation needs to be obtained by a PODRIX user.

Due to the varying sun incident angle and payload activities, the front-end temperature of the receiver exhibits variations in the range of about 16°C to 21 °C over the evaluated period. The temperature variations cause signal-dependent differential code biases with meter-level amplitudes, which have been characterized in a ground-based calibration. Together with the epoch-wise front-end temperature provided in the receiver telemetry, the calculated ranges are corrected accordingly.

### 3.2 Receiver Clock Model

For the simulated real-time processing performed in this study, a receiver clock model was introduced that aimed to align the internal receiver time estimate *t*_{rcv} as closely as possible with the GPS system time. Compared to a free-running receiver clock, a GPS-aligned timescale facilitates the synchronization of observations to integer seconds of GPS or UTC time as well as the generation of a corresponding pulse-per-second signal. An overview of the timescales used in this study is given in Table 2.

As described earlier, each measurement epoch is tagged with an IMT value representing the current discrete counter value *b* of the upconverted oscillator frequency in the AGGA chip of the PODRIX. Therefore, the elapsed time between two consecutive epochs *i*−1 and *i* can be described with:

3

where *f*_{core} is the core frequency of the AGGA chip’s hardware counter. However, *τ _{i}* is affected by the oscillator’s drift between the two epochs with respect to GPS time. To compensate for this timing error, a complementary first-order receiver clock offset polynomial is introduced. The polynomial’s constant term is defined as:

4

and describes the residual offset between the receiver and GPS timescale by the speed of light *c*, while the linear coefficient:

5

is the oscillator’s fractional frequency offset. The resulting polynomial can be combined with Equation (3) to describe the receiver’s time propagation between the two epochs with:

6

While *τ _{i}* is calculated directly from the IMT count at the measurement epoch, the clock model coefficients need to be estimated by the Kalman filter.

Changes of the receiver timescale directly affect the measurement generation. For the discussion of clock model estimation, we assume that an EKF measurement update at epoch *i*−1 has just been performed. At this point, the a-posteriori coefficients and with covariances and have been obtained by the Kalman filter based on the pseudorange and carrier-phase observations **ρ**_{i−1} and *ϕ*_{i−1}. The observables were formed with a receiver time predicted prior to the filter’s measurement update stage. The residual time prediction error is compensated by applying the estimated clock offset with:

7

to obtain the receiver’s best estimate of the GPS time at epoch *i*−1. The observables are consequently adjusted to the new receiver time estimate with:

8

These corrected observables can be used to calculate a precise solution in post-processing, yielding the remaining clock error which the EKF is unable to compensate in real time and, thus, determining the estimation performance.

When advancing to the next epoch *i*, the clock coefficients are propagated with:

9

during the time update stage, leading to an a-priori prediction for the new receiver time at *t _{i}* based on Equation (6) with:

10

This timestamp is used to form the new pseudorange * ρ_{i}* and carrier phase

*from the corresponding raw code and NCO phase as shown in Equations (1) and (2). Proceeding to the measurement update stage, new clock coefficients and are obtained, and the clock and observables update process for epoch*

**ϕ**_{i}*i*can be performed with Equations (7) and (8) in preparation for the next filter epoch

*i*+ 1.

### 3.3 Data Processing

Data processing is done with RTNAV (Hauschild & Montenbruck, 2021; Montenbruck & Ramos-Bosch, 2008), a multi-GNSS-capable software for LEO satellite navigation in a simulated real-time mode. It is based on an extended Kalman filter (EKF), using pseudoranges and carrier phases as measurement input and a reduced-dynamics orbit model (Wu et al., 1991) for propagation of the satellite’s state vector. While designed for playback processing of pre-recorded observations and ephemeris data, RTNAV operates in a purely sequential manner and makes exclusive use of auxiliary data available in an onboard environment. It is, thus, fully representative of an actual real-time navigation processor in terms of the achieved orbit and clock determination performance. Hauschild and Montenbruck (2021) demonstrated the algorithm’s capability to reach decimeter-level accuracy for real-time orbit determination using multi-GNSS measurements and broadcast ephemerides.

For this study, RTNAV was modified to enable replay and processing of raw Sentinel-6A telemetry files in binary format and real-time epoch-wise observable calculation. Therefore, the overall software creates a realistic simulation environment for a real-time onboard navigation system which is tightly integrated with a GNSS receiver. The general schematics of this process is shown in Figure 2, while algorithms and models are summarized in Table 3.

The EKF state vector is defined as:

11

where * r* and

*are the satellite’s position and velocity in the International Celestial Reference Frame (ICRF),*

**v***C*and

_{R}*C*are the dimensionless radiation and drag coefficients, and

_{D}

**a**_{emp}represents the estimated empirical accelerations that are not incorporated in the force model. The vector

*contains carrier-phase ambiguities which are estimated for each receiver channel. Furthermore, the filter state includes the receiver clock offset*

**N***a*

_{0}and the fractional frequency offset

*a*

_{1}, both scaled by the speed of light and referred to GPST. The inter-system bias ISB

_{GAL}is further elaborated in the next section.

In common navigation filters, where only position and velocity are of concern rather than timing, the clock offset is often estimated as a free parameter. This means its covariance is reset to an arbitrarily large value prior to the measurement update to allow the filter to estimate a new clock offset at every epoch solely on the basis of the observables. A free clock offset estimation is generally adequate for low-grade oscillators, but causes an increased sensitivity of the resulting estimate to broadcast ephemeris errors. For high-precision oscillators, in contrast, realistic process noise values for the predicted clock states may be applied based on the oscillator stability (Weinbach & Schön, 2011), describing the relative uncertainty of the predicted clock offset between two consecutive estimation epochs. Using the oscillator’s Allan deviation (Allan, 1966), the clock estimation can be constrained with (Wang & Rothacher, 2013):

12

where σ_{ADEV,τ} is the Allan deviation at sample time *τ _{s}*, scaled to distance units with the speed of light

*c*. By way of example, an Allan deviation of 10

^{−12}at

*τ*= 30 s translates into a 1-cm uncertainty of the predicted clock state, which would typically be less than the clock estimation uncertainty introduced by broadcast ephemeris errors. The applicability of clock model constraints for the DORIS USO is further discussed in Section 4.3.

_{s}In preparation for the measurement update stage, ionosphere-free combinations are formed from dual-frequency observations. The resulting signal combinations used as measurements in the EKF are listed in Table 4 along with the corresponding identifiers of the Receiver-INdependent EXchange format (RINEX; Romero, 2020). The ionosphere-free code and phase ranges are analyzed prior to the measurement update stage by calculating the residuals relative to the expected geometric ranges. With this adaptive data editing, outliers above pre-defined thresholds are detected and can be excluded from the measurement update to improve the estimation results.

#### 3.3.1 Time References and Bias Handling

The usage of dual-frequency observations from both GPS and Galileo requires careful consideration of different biases to obtain a consistent timescale based on all available observations. All GNSSs realize their own timescale based on a specific set of signals. Galileo relates the FNAV broadcast time to the E1/E5a signal combination, while GPS uses L1/L2 P(Y) as clock reference signals (Bar-Sever, 2020). Hence the satellites broadcast their clock offsets with respect to their corresponding reference signals. As different codes and frequencies cause individual instrumental delays in the satellite hardware, each of the other available GNSS signals has a bias relative to its respective reference signal. Therefore, when using such a signal, satellite-specific differential code biases (DCBs) need to be applied to correct those differences and align the signal with the system clock reference.

As seen in Table 4, the GPS Block IIR and Galileo measurements used for this study are already consistent with the clock reference signals and, therefore, do not need any further correction. However, the GPS L1 C/A and L2-CL signals need to be corrected with the respective DCBs relative to the L1/L2 P(Y) signals. GPS provides a timing group delay (TGD) as well as additional inter-signal correction (ISC) parameters for satellites of Block IIR-M and later generations. The ISC* _{x}* describes the instrumental delay difference between L1 P(Y) and the corresponding signal

*x*for a certain satellite, while the TGD describes the difference between the L1 P(Y) signal and the L1/L2 P(Y) combination (Tetewsky et al., 2009). As soon as any signal other than L1/L2 P(Y) is used, the corresponding ISCs and TGD are needed even for dual-frequency processing to remove any signal-dependent DCBs. In accordance with a real-time environments, the algorithm uses broadcast navigation messages as a data source and extracts the required TGD and ISC parameters from the GPS CNAV messages.

While the DCBs allow for the alignment of any signal of a GNSS with respective reference signals of the same GNSS, differences between individual GNSS need to be considered as well. Under nominal conditions, GPS system time (GPST) and Galileo system time (GST) differ at the nanosecond level with slow variations. In addition to timescale offsets, constellation-dependent biases in the receiver hardware need to be considered. The combined effect of system time differences and receiver-specific biases is commonly compensated by choosing one GNSS as time reference and estimating an inter-system bias (ISB; Hauschild, 2017), also called user Galileo-to-GPS Time Offset (GGTO; Defraigne et al., 2021), relative to this system for all remaining GNSS.

For this study, GPST was chosen as the reference time and an ISB for Galileo relative to GPS was estimated at each epoch so that the relation between Galileo and GPS clock offsets can be described as:

13

where *dt _{i}* denotes the receiver clock offset specific to observations of constellation

*i*and ISB

_{GAL}is the ISB of Galileo relative to GPS. The ISB-corrected clock offset is applied in the measurement model for all Galileo satellites.

As an alternative to estimating an ISB, the modernized navigation messages of most GNSSs include X-to-Y timing offset (XYTO) information, which describes the time difference between constellations X and Y. For Galileo, the GGTO is transmitted via INAV and FNAV messages (Galileo ICD, 2021). It describes the predicted GST-GPST difference and is updated on a daily basis. A GGTO correction polynomial (albeit with an opposite sign convention) is also defined in the modernized civil navigation messages of GPS (GPS ICD, 2020), but not yet part of the current pre-operational CNAV transmission.

Neglecting receiver-specific biases as well as the errors of the GGTO determination and prediction, itself, the GGTO as broadcast in the Galileo navigation messages matches the negative value of the ISB as defined in Equation (13). The officially stated accuracy of Galileo’s GGTO values ranges from 2.2 ns to 4.6 ns during the first quarter of 2021 (GSA, 2021). Sesia et al. (2021) showed that the broadcast GGTO varied approximately ± 5 ns over the course of one month in 2018. The GGTO provided in the INAV and FNAV messages is represented by a first-order polynomial, allowing the user to obtain an epoch-wise value for time correction relative to GPS. Note that any receiver-dependent biases are neglected if only the GGTO is applied to the observables instead of estimating an ISB. Furthermore, consecutive sets of the daily updated GGTO parameters can have differences of several nanoseconds and cause discontinuities at such transitions.

Defraigne et al. (2021) assessed the impact of the GGTO on the navigation solution and came to the conclusion that, for high-precision GNSS receivers, a better performance could be achieved when estimating the ISB rather than fixing it to a known value. In the present study, a sufficient number of GNSS satellites from both constellations was always available to enable precise estimation of the ISB from dual-frequency code and phase observations within the navigation filter. In view of the daily GGTO discontinuities and the presence of remaining receiver-specific biases at the few-nanosecond level, incorporation of the broadcast GGTO provides no practical benefits in the present context. For further information, a comparison of estimated ISB time series with the broadcast GGTOs is given in Section 4.2.1.

### 3.4 Reference Orbit and Clock Solution

The estimation results for orbit and clock data obtained from online processing with RTNAV need to be compared against a precise reference solution to analyze the estimation performance. For this purpose, the observables calculated with the predicted GPS time from Equation (10) and corrected with the estimated offset after the measurement update as shown in Equation (8) are logged to a RINEX file during the execution of the real-time navigation filter. These corrected observables and the associated receiver time represent the best real-time estimate of the GPS time at each measurement epoch. By performing a POD using these observables in combination with precise GNSS orbit and clock products, the error between the estimated GPS time and the timescale of the precise GNSS clock product can be determined.

The post-processing was performed with a batch least-squares algorithm, a more complex reduced force model, and precise orbit, clock, and bias products from the Center for Orbit Determination in Europe (CODE; Dach et al., 2020; Schaer et al., 2021). The CODE clock products are based on observations of a global monitoring network and are closely aligned with the GPS broadcast timescale. More specifically, epoch-wise receiver clock offsets of the individual receivers from GPS broadcast time are first determined through single point positioning with broadcast ephemerides and later refined through carrier-phase observations in the orbit determination and time synchronization process (Dach et al., 2015).

In this way, an intermediate timescale representing the ensemble mean of all station-wise GPS broadcast time realizations is obtained. After selecting a time reference station equipped with a highly stable hydrogen maser, the intermediate reference timescale was adjusted in such a way that the clock offset of the reference station was rigorously represented by a first-order polynomial in time. The final CODE timescale, thus, reflects the physical short-term behavior of the reference station but is closely aligned to the GPS broadcast timescale over a day.

Details of the algorithms for precise orbit and clock offset determination of Sentinel-6A are discussed in Montenbruck et al. (2021). Even though the reference POD solution is based on the same GNSS observations as the playback real-time solution, it benefits from the use of more accurate GNSS orbit and clock products and the simultaneous adjustment of all observations in a least-squares process. Based on validation with satellite laser ranging measurements, a 1D RMS position accuracy at the 1-cm level may be inferred, which outperforms the real-time navigation solution by roughly one order of magnitude.

## 4 RESULTS AND DISCUSSION

For the evaluation of the previously described online processing algorithm, a continuous 14-day arc of telemetry data from Sentinel-6A was used. It contained all PODRIX data from January 17, 2021, to January 30, 2021. Raw measurements were decimated to 30-s sampling for the present study in accordance with the desired update rate of the navigation filter. The telemetry records exhibited an 18-min data gap on January 25 starting at 02:29. The EKF was not reinitialized and solely relied on propagation during this period.

The results provided by the EKF include the estimated orbit, DORIS USO drift, as well as the clock offset, which depicts the clock model correction term. The timescale investigated here is based on the corrected time after measurement update as presented in Equation (7).

This section starts with a brief presentation of the orbit determination results compared to a precise reference orbit. Thereafter, the estimated timescale and its error with respect to precise clock products is presented and discussed in terms of precision and accuracy, focusing on the influence of different biases on the final result. For the clock model, an analysis regarding its constraints and limitations is presented, particularly under the influence of the SAA. Finally, a brief outlook of UTC synchronization is given.

### 4.1 Orbit Estimation

The estimated real-time orbit of Sentinel-6A was compared to the precise post-processed reference solution along the satellite’s radial-transverse-normal axes and is shown in Figure 3. After a short convergence time, the error in all three axes remained well below 25 cm, with standard deviations of 39 mm, 67 mm, and 80 mm for radial, along-track, and cross-track axes. Errors in the radial and cross-track directions exhibited an obvious once-per-revolution periodicity, which is characteristic for satellites in near-circular orbits and can be understood by the respective orbit dynamics (Colombo, 1989). In the along-track direction, which was less tightly constrained by the dynamic model, more rapid variations of the position error could be observed, but the overall errors were of a similar magnitude as for the other components. The data gap of around 17 minutes on January 25 did not cause visible outliers, demonstrating the algorithm’s capability to cope even with larger data outages without significantly deteriorating the orbit determination performance.

### 4.2 Time Synchronization

The timestamps derived from the oscillator ticks and corrected by the residual clock offset as shown in Equation (7) represent the best estimate of the GPST at each measurement epoch and form a continuous timescale over the whole evaluation period between January 17 to January 30, 2021. For the reference solution, daily CODE “_M” orbit, clock, and bias products with a 30-s clock sampling rate were used. Thus, the clock error analyzed in the following sections describes the difference between the approximate GNSS timescale *t*_{rcv} realized by the EKF based on GPS/Galileo broadcast ephemerides and the inherent timescale of the CODE clock products *t*_{CODE}.

This section covers the problem of combining measurements from both GPS and Galileo and how filter settings affect the final timescale estimate, followed by an analysis of the clock estimation error in terms of stability and accuracy.

#### 4.2.1 Combination of Different Timescales

By tracking signals from both GPS and Galileo, two different system timescales were observed that needed to be combined into a single receiver timescale. The system time differences, as well as further system-specific receiver biases, were considered by introducing an ISB as an additional filter state. As described in Section 3, the ISB lumps together system time differences and GNSS/signal-specific receiver biases. The actual time difference between GPST and GST is known to be at the few-nanosecond level with slow variations (Sesia et al., 2021). Furthermore, post-processing showed that the physical biases of the PODRIX receiver between Galileo and GPS observations were almost constant over one day (Montenbruck et al., 2021).

Figure 4 shows the estimation error for two clock solutions referring to GPST and GST, respectively. Both solutions used otherwise identical filters and incorporated all available measurements as listed in Table 4, while the ISB was propagated as a kinematic parameter with an arbitrarily large process noise.

The GPST estimate shows an approximately constant phase offset of −4.1 ns over the 14-day period. Systematic variations with peak-to-peak amplitudes of 2–4 ns may have been observed between consecutive midnight epochs, which exhibited an almost parabolic shape centered around the noon epoch. These variations are assumed to represent actual differences in the realization of the broadcast timescale of GPS. While their exact cause could not be fully identified from the present data, they may have originated from either the clock steering process in the GPS control segment, or small errors in the prediction of quadratic satellite clock model terms.

It is instructive to note that similar variations were earlier reported for GPST-UTC(NIST) differences determined with a GPS receiver at the National Institute of Standards and Technology (Levine, 2020). Compared to the GPS-based receiver clock offset, no obvious quadratic variations across a day could be recognized in the Galileo-based clock estimate. This suggests that the observed features of the GPS-based Sentinel-6A real-time clock solution are related to the realization of the GPS broadcast timescale, which ultimately limits the stability of GPST access for broadcast ephemeris users.

The GST-based solution, on the other hand, shows long-term variations, reflecting the differences between GPS and Galileo system timescales. These variations would not show up when comparing against a precise GST-referenced GNSS clock product. However, no such product is currently available to the public and, therefore, the following combined clock solutions discussed in this paper all refer to GPST.

Dual-constellation processing requires careful tuning of the ISB process noise to exploit the full set of multi-GNSS observations in the realization of a GPST-aligned receiver timescale. While the daily variations need to be reflected by the ISB, short-term noise up to one orbit revolution of Sentinel-6A are unrelated to system time differences. By accounting for the actual ISB stability in the navigation filter, the overall impact of observation noise and measurement model errors on the clock offset estimation can be reduced.

As the state model assumes the ISB to be constant between consecutive epochs, the process noise needed to reflect the expected slow ISB variations. A suitable value for the ISB process noise was determined with a parametric search at σ_{ISB} = 6.0 mm over 30 s, with the results shown in the second row of Figure 5. With this setting, the clock solution would benefit from the Galileo measurements and achieve a lower standard deviation of 0.87 ns compared to the kinematic results with 1.04 ns. Reducing the process noise further to 0.5 mm already induces slight long-term deviations from GST. Such a setting is, therefore, already assumed to be over-constrained, despite having a lower clock error standard deviation of 0.70 ns over the evaluation period. In general, lowering the ISB process noise increases the weight of the Galileo measurements, with the extreme case of zero process noise rendering the estimated timescale to be a hybrid solution which effectively averages the observed GPST and GST.

While the receiver cannot directly estimate its own system-related inter-signal biases, they can be roughly characterized by comparing the estimated ISB to the broadcast GGTO. The right-hand side plots in Figure 5 show the broadcast GGTO with an inverted sign to match the ISB definition. During the evaluation period, the GGTO polynomial parameters were updated once a day at approximately 12:30 pm. Those transitions can be seen clearly creating jumps of up to 2 nanoseconds. The estimated ISB tends to follow the GGTO with an average offset of approximately 1 ns. Although the results indicate very low time variations for the receiver biases, the 14-day evaluation period was too short to draw conclusions on whether this good match between ISB and GGTO would be permanent.

#### 4.2.2 Timescale Stability and Accuracy

The stability of the real-time clock solution *t*_{rcv} with respect to the CODE reference timescale can be visualized with the time deviation (TDEV), which is based on the modified Allan deviation (Allan, 1966). It indicates time stability with respect to a reference clock for different sample intervals *τ*.

The previously identified ISB process noise setting of σ_{ISB} = 6.0 mm was found to be a good compromise between noise reduction and low influence of GST on the GPST-based receiver timescale estimation, with the resulting TDEV shown in Figure 6. A first peak was observed at approximately 3,400 s and a minimum at 6,800 s. These correspond to a half and a full orbital revolution of Sentinel-6A and were caused by the already mentioned once-per-revolution variations visible in both orbit and clock solutions. The daily quadratic variations caused by limitations in the broadcast GPST realization are also reflected in the TDEV, with clearly visible local maxima at 12 h and 36 h, and minima at 24 h and 48 h. Overall, the TDEV indicates a standard deviation of the clock solution below the nanosecond level at all sample intervals, with the long-term stability improving as expected due to continuous time estimation with GNSS measurements.

In contrast to stability or precision, the accuracy represented by the −4.1 ns phase offset was a result of external factors rather than the estimation algorithm. As explained earlier, a GNSS user has only indirect access to the respective GNSS timescales via broadcast ephemerides and clock corrections, and may suffer from a lacking calibration of receiver-related biases. This also applies to the International GNSS Service (IGS; Johnston et al., 2017), to which the CODE products contribute.

As described in Section 3.3.1, the various GNSS timescales all refer to a specific signal or pair of signals, which in the case of GPS was L1/L2 P(Y) and in the case of Galileo was E1/E5a. While the latter signals were directly observed in this study, the GPS L1 C/A and L2C signals need to be aligned with the clock reference signals using the respective TGD and ISC parameters. Consequently, the GPST estimation in this study was a hybrid solution which referred to the P-code signals via a combination of civil signals and broadcast DCB corrections. Providers of precise products, on the other hand, usually estimate the DCBs from observations of all available signals along with their orbit and clock products. The reference solution in this study used CODE bias products, which are published in Bias Solution INdependent EXchange Format (Bias-SINEX). Therefore, different DCBs were used for the real-time and post-processed solution. While the DCBs had no significant effect on the position solution and clock standard deviation, they applied a certain bias on the clock solution (i.e., an approximately constant clock phase offset).

Therefore, the mean clock error of −4.1 ns as seen in Figure 5 is not a direct quality measure for the clock solution, but rather a result of different biases, observations, and reference times employed in the real-time and precise products. The major part of the mean error can be traced back to two different sources.

First, the CODE clock products showed an almost constant bias of approximately −2 ns compared to the GPS broadcast clock offsets during the evaluation period, indicating that the CODE timescale realization reproduced the average frequency of GPST, but not necessarily the absolute phase. The second part was caused by the different bias products, namely broadcast TGDs and ISCs for the real-time estimation and the CODE DCBs for post-processing, which also show an average difference of around −2 ns over the 14 days. The two errors sum up to approximately −4 ns, closely matching the observed mean clock error and, thus, providing indirect evidence that the PODRIX receiver exhibits a factory calibration with an accuracy at the few-nanosecond level.

Eventually, it is important to emphasize that the onboard time as well as CODE clock products were only estimates of the real and only indirectly accessible GPST. Both timescale realizations took different approaches regarding tracked signals and auxiliary information such as DCBs, ultimately leading to differences in the form of an overall bias or phase offset. If a user wants to avoid such a bias, a calibration of the whole receiver and processing chain with respect to the target reference timescale would be required.

### 4.3 Clock Model Constraints and Impact of the South Atlantic Anomaly

A specific characteristic of LEOs is the presence of the SAA region, where satellites are exposed to increased radiation doses. The fractional frequency offset in Figure 1 already showed that despite the DORIS USO’s improved radiation hardening (Auriol & Tourain, 2010), the SAA passes still cause rapid changes in the oscillator’s frequency.

In view of its high frequency stability, the DORIS USO’s drift is assumed constant between two consecutive measurement epochs of 30 s, which is also reflected in the EKF’s system model in Equation (9). Similar to the ISB propagation which was discussed in Section 4.2.1, the actual changes in the clock drift or frequency needed to be compensated for by choosing an appropriate process noise for the clock model.

Constraining the clock according to its physical properties as described in Wang and Rothacher (2013) and shown in Equation (12) can improve the estimation performance if the oscillator stability is sufficiently high. In practice, this would imply a TDEV of the oscillator lower than the TDEV of the unconstrained clock estimate shown in Figure 6.

However, Weinbach and Schön (2011) noted that OCXOs typically do not perform well enough to effectively apply clock constraints due to flicker frequency noise. Several tests with the experimental setup used for this study confirm this statement, as clock constraints could not significantly improve the overall estimation performance. The carrier-phase residuals in Figure 7 show that over-constraining leads to increased errors in the observation model, with the strongest impact during SAA transitions, consequently degrading both orbit and clock solutions.

Therefore, we conclude that despite the excellent performance of the DORIS USO, deviations from a linear clock phase variation are still too large in comparison with the carrier-phase precision to benefit from clock constraints. To better cope with the SAA effects, more sophisticated clock models could be used to incorporate the rapid frequency changes. In the case of Sentinel-6A’s DORIS USO during the evaluation window, a linear frequency increase could be assumed for the nominal case, and the SAA spikes would be modeled based on the satellite’s current position and a physical model (Lemoine & Capdeville, 2006).

While the clock constraints could not improve the estimation of the clock offset, itself, they enabled estimation of the clock rate and, thus, improved the prediction of the modeled clock offset in the time update step of the Kalman filter. This is illustrated in Figure 8, which shows the distribution of clock offset corrections obtained in the measurement update steps. Case 1 corresponds essentially to a free estimation, while Case 2 applies moderate clock constraints which still allow the clock model to cope with frequency spikes during SAA transitions. Case 3 finally represents an over-constrained filter in which confidence in a linear clock variation is no longer in accord with the actual changes of the clock drift.

When moderately constraining the clock model, as in Case 2, the distribution narrows, meaning that the a-priori prediction from Equation (10) would become more precise, while the final orbit and clock solution would remain similar to Case 1. However, over-constraining the clock model, as in Case 3, significantly degrades the quality of both orbit and clock solution.

The time prediction between consecutive 30-s measurement epochs range in the order of tens-of-picoseconds, and is mainly a result of the high-frequency stability of the DORIS USO. When using less stable oscillators, the prediction error is expected to increase, as integrating the clock drift between two measurement epochs cannot accurately describe the actual behavior of the oscillator.

### 4.4 Synchronization to UTC

The algorithm presented here provides a precise referencing of the individual oscillator clock beats to the GPS timescale as realized through GNSS observations and broadcast ephemerides. However, some applications might require a synchronization to UTC. All GNSSs maintain their own timescales, which differ from UTC mostly in the order of nanoseconds, neglecting the known integer leap second offset since the respective first epoch. According to GPS ICD 200L (GPS ICD, 2020), the GPST is kept within a 1*σ* standard deviation of 20 ns with respect to UTC. The GST, on the other hand, is guaranteed to stay within 30 ns around UTC, and showed an average difference from UTC to be around 4.3 ns for the first quarter of 2021 (GSA, 2021).

Predicted GNSS-UTC offset parameters are broadcast via the system’s respective navigation messages to allow a user to derive UTC based on the obtained GNSS time. Sesia et al. (2021) examined the quality of broadcast GNSS-UTC time offset predictions and the feasibility of using UTC as a common timescale for multi-GNSS processing. They specifically highlighted the fact that broadcast GNSS-UTC offsets of various constellations refer to different UTC realizations. Therefore, additional external information is required to reference all GNSS observations to a common UTC timescale. Using such information as well as a calibrated receiver, consistency could be achieved at the few-ns level for GPS and Galileo, while BeiDou and GLONASS would exhibit deviations of up to 10 ns in the examined 1-month time interval.

Based on these results, UTC synchronization of a LEO-based oscillator using multi-GNSS observations appears to be a potential alternative to the use of a single-GNSS timescale, but requires additional information or a better harmonization of the GNSS-specific reference UTC realizations.

## 5 SUMMARY AND CONCLUSION

GPS time synchronization of a LEO satellite clock with sub-nanosecond precision was demonstrated using multi-GNSS observations collected aboard the Sentinel-6A satellite. The spacecraft was equipped with a geodetic-grade PODRIX GNSS receiver connected to a USO. Dual-frequency code and carrier-phase observations of GPS and Galileo enabled a ground-based POD with few-cm accuracy. In contrast to the GPS-disciplined clocks used for terrestrial time synchronization applications, LEO satellites are fast-moving objects and accurate knowledge of their position is a prerequisite for GNSS-based clock synchronization. We made use of actual Sentinel-6A GNSS flight data to study joint orbit determination and time synchronization of the attached USO in a playback processing that emulated the function of a real-time navigation system embedded in the GNSS receiver. The processing was based on an EKF with a reduced-dynamic orbit model considering various perturbations.

Measurements processed in the filter include pseudorange and carrier phase of civil signals and P(Y)-code of GPS L1/L2, as well as Galileo Open Service E1/E5a signals. The achievable real-time positioning accuracy depended critically on the respective broadcast ephemerides’ quality and reached a 1-decimeter 3D RMS level in the combined GPS/Galileo processing. While use of Galileo is crucial from a performance point of view, the addition of GPS offers a larger overall number of visible satellites. Furthermore, GPS provides a practical benefit for real-time applications by broadcasting Earth orientation parameters (EOPs) for use in reference system transformations, and finally, offers access to the GPS system time.

The USO time synchronization was accomplished by associating the oscillator beat count at the instant of the GNSS observations with the corresponding estimate of GPS time obtained in the navigation filter. A linear model was used to describe the variation of the oscillator-based receiver timescale between consecutive measurements and an ISB with properly tuned process noise adjusted to combine observations from multiple constellations. The filtered clock offset and rate of a chosen reference constellation were used to continuously align the receiver timescale to the best available GNSS time estimate. The resulting receiver time was, furthermore, used to form pseudorange and carrier-phase observations from the raw code and NCO correlator measurements provided by the PODRIX receiver. These can subsequently be used in a POD process to determine the residual error of the receiver timescale relative to precise GNSS clock products as provided by the IGS and its analysis centers.

Based on the processing of a continuous 2-week data arc, a 0.9-ns precision of the receiver time alignment relative to the GPS-based timescale realized by a precise clock product of the CODE analysis center was obtained. The achieved precision was partly limited by apparent quadratic variations of the GPS-related receiver timescale between consecutive midnight epochs. These are, likewise, present in the estimated Galileo-GPS ISB and must be attributed to inherent GPS broadcast timescale variations. Small inaccuracies in the broadcast values of the second-order clock offset coefficients might offer a possible explanation, but further investigation would be required to fully identify the root cause of the observed quadratic variations in the GPS-based receiver timescale. The estimated ISB was consistent with the broadcast GGTO to within about 1 ns, but did not exhibit linearization errors and discontinuities. This indicates that the PODRIX receiver is largely free of systematic inter-system signal biases and suggests that the estimated ISB represents a good measure of the actual broadcast system time differences between GPS and Galileo.

The mean difference between the estimated receiver timescale and the timescale of the CODE reference clock product amounts to approximately −4 ns, which dominates the overall accuracy and highlights the problem of systematic offsets in individual timescale realizations. The PODRIX was not intended to serve as a timing receiver, and the CODE clock product timescale was not designed to ensure a well-defined offset from a GNSS system timescale or UTC, but only aimed to stay within a few nanoseconds of GPS time. Further work will, therefore, be required to establish an absolute synchronization between LEO satellite clocks and a GNSS or UTC timescale.

In this context, the generation of precise GNSS clock products closely aligned to a GNSS system timescale by the scientific community is strongly encouraged. In view of apparent sub-daily variations in the realization of the GPS broadcast timescale, Galileo offers itself as a new constellation with potentially more precise access to the actual GNSS system time through its superior broadcast ephemerides. Use of Galileo as a primary reference for precise clock products should, therefore, be considered as an alternative to the present GPS alignment.

Overall, the present study demonstrated that LEO satellite clocks could be synchronized in real time with a common timescale with nanosecond precision using onboard GNSS receivers and a proper navigation filter for real-time orbit determination. This is of interest for, e.g., bistatic synthetic aperture radar (SAR) satellite formations that require coordinated operation of their instruments, but also the proposed use of LEO satellite constellations as GNSS augmentation systems. Availability of GNSS-disciplined onboard clocks with properly calibrated receivers would enable the dissemination of synchronized ranging signals with representative user range errors at the half-meter level without a need for costly atomic clock standards or excessive ground infrastructure for orbit and time synchronization. This would allow for seamless integration of LEO navigation satellites into existing GNSS architectures and offer improved satellite visibility, availability, and dilution of precision even in challenging environments.

## HOW TO CITE THIS ARTICLE

Kunzi, F., & Montenbruck, O. (2022) Precise onboard time synchronization for LEO satellites. *NAVIGATION, 69*(3). https://doi.org/10.33012/navi.531

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.