Real-Time Estimation of LEO Satellite Clocks Using a Combination of Low- and High-Frequency Clock Solutions

2.1 Low-Frequency LEO Satellite Clock Estimation and Prediction

The reduced-dynamic BLS POD process is often used to produce high-precision LEO satellite orbits and clocks, owing to its higher accuracy and stability compared with kinematic or dynamic POD methods (Bock et al., 2011; Mao et al., 2021; Wang et al., 2023a). By using onboard GNSS dual-frequency ionosphere-free (IF) carrier phases and code observations and integrating existing dynamic models, a series of dynamic parameters (X_D), the float IF ambiguities (N_IF), and the estimable LEO clock (∆t_r), we obtain an estimate for epoch t_i as follows:

$E(\Delta p_{r,IF}^s(t_i\left)\right)={\left(A_{r,D}^s\right(t_i\left)\right)}^T{X}_D+c\times \Delta t_r\left(t_i\right)$ 1

$E(\Delta {\phi }_{r,IF}^s(t_i\left)\right)={\left(A_{r,D}^s\right(t_i\left)\right)}^T{X}_D+c\times \Delta t_r\left(t_i\right)+{\lambda }_{IF}\times N_{IF}$ 2

where $\Delta p_{r,IF}^s$ and $\Delta {\phi }_{r,IF}^s$ denote the observed-minus-computed (O-C) term of the IF combined carrier code and phase observations, respectively. X_D represents the estimable dynamic parameters, including the six orbital Keplerian parameters at the initial time, the solar radiation pressure (SRP) parameters, and a series of stochastic accelerations. More details have been given by Dach et al. (2015). $A_{r,D}^s$ contains the partial derivatives of the observations with respect to X_D. c is the speed of light, λ_IF is the wavelength of the IF combination, and E(∙) represents the expectation operator. The subscript r and the superscript s denote the LEO (as the user here) and GNSS satellites, respectively. The state vectors of the BLS ${\hat{X}}_{BLS}$ can be expressed as follows:

${\hat{X}}_{BLS}={(\Delta K_1,\dots ,\Delta K_6,\Delta D_1,\dots ,\Delta D_q,\Delta Q_1,\dots ,\Delta Q_p,N_{r,IF}^1,\dots ,N_{r,IF}^z,\Delta t_{1,r},\dots ,\Delta t_{u,r})}^T$ 3

where ∆K_i (i = 1, ⋯ ,6) denotes the six orbital Keplerian parameters at the initial time, ∆D_q denotes the SRP parameters, and q denotes the number of SRP parameters. ∆Q_p represents the stochastic piecewise constant accelerations in the radial (R), along-track (S), and cross-track (W) directions, and p denotes the number of stochastic accelerations. z and u represent the number of float IF ambiguities and estimated LEO clock, respectively.

With the high-precision LEO satellite clocks determined in near-real time within a BLS reduced-dynamic POD process, the clock prediction is achieved by using an m-degree polynomial and k harmonic series for a period of up to 15 min, considering various periodic effects within the clock estimates as follows (Wang & EI-Mowafy, 2022):

$\widehat{ClK}(t_p-t_0)={\hat{a}}_0+{\hat{a}}_1(t_p-t_0)+\cdots +{\hat{a}}_m{(t_p-t_0)}^m+{\sum }_{i=1}^k{\hat{A}}_i\sin \left(\frac{2\pi }{{\hat{T}}_1}\right(t_p-t_0)+{\hat{\phi }}_i)$4

where $\widehat{ClK}$ denotes the predicted clocks and t₀ and t_p denote the initial and prediction epochs. ${\hat{a}}_j(j=1,\cdots ,m)$ represents the m polynomial fitting coefficients (here, m is 1). k is the number of periodic terms (considered here as 5), and ${\hat{T}}_l$, ${\hat{A}}_i$, and ${\hat{\phi }}_i$ represent the period, amplitude, and phase of the periodic terms, respectively.

2.2 High-Frequency LEO Satellite Clock Estimation with Constrained Predicted Clocks

High-frequency LEO satellite clocks are typically estimated in real time by using dual-frequency GNSS IF carrier phase and code observations, employing methods such as the Kalman filter or sequential least-squares (SLS) adjustment. These observations can be expressed as follows:

$E(\Delta p_{r,IF}^s(t_k\left)\right)={\left(A_{r,K}^s\right(t_k\left)\right)}^T{X}_K\left(t_k\right)+c\times \Delta t_r\left(t_k\right)$ 5

$E(\Delta {\phi }_{r,IF}^s(t_k\left)\right)={\left(A_{r,K}^s\right(t_k\left)\right)}^T{X}_K\left(t_k\right)+c\times \Delta t_r\left(t_k\right)+{\lambda }_{IF}\times N_{IF}\left(t_k\right)$ 6

where X_K denotes the three-dimensional orbital correction in the Earth-centered Earth-fixed system. $A_{r,K}^s$ contains the partial derivatives of the observations with respect to X_K. Temporal constraints are placed on N_IF, expressed as follows:

$N_{IF}\left(t_k\right)={\hat{N}}_{IF}\left(t_{k-1}\right)$7

$Q_N\left(t_k\right)=Q_{\hat{N}}\left(t_{k-1}\right)$8

The ambiguities to be calculated in this epoch N_IF (t_k) are constrained to those estimated in the previous epoch ${\hat{N}}_{IF}\left(t_{k-1}\right)$ when a cycle slip does not occur. Q_N (t_k) and $Q_{\hat{N}}\left(t_{k-1}\right)$ are the variance-covariance matrices of the estimated ambiguities in the current and previous epoch, respectively.

In addition to the temporal constraints mentioned above, the predicted clocks $(\Delta {\hat{\tilde{t}}}_r(t_k\left)\right)$ over a certain prediction period based on near-real-time BLS reduced-dynamic POD (see Section 2.1) are introduced and used for constraining the high-frequency clock estimates. To constrain a predicted clock, the estimated clocks ∆t_r are set equal to the predicted clocks $\Delta {\hat{\tilde{t}}}_r$ as pseudomeasurements associated with a proper variance, expressed as follows:

$\Delta {\hat{\tilde{t}}}_r\left(t_k\right)=\Delta t_r\left(t_k\right)$ 9

The covariance matrix Q_t (t_k) for the constraint presented in Equation (9) affects the influences of the introduced predicted clocks, where a small variance indicates a strong constraint and a large variance implies a weak constraint. Thus, the variances must be defined based on the precision of the predicted clocks σt_p according to their prediction time t_p. However, because the introduced predicted clock errors are not white noise but biases, the variance of the constraint in Equation (9) may need to be increased to achieve a better clock solution. Consequently, different variances are tested for the constraint in Equation (9) for the predicted clocks introduced from each tested prediction window. This step will be discussed in detail in Section 3.2. The state vector for the Kalman filter or SLS adjustment $\hat{X}\left(t_i\right)$ can be expressed for epoch t_i:

$\hat{X}\left(t_i\right)={\left(X_K\right(t_i),\Delta t_r(t_i),N_{r,IF}^1,\dots ,N_{r,IF}^w)}^T$10

where w denotes the number of float IF ambiguities in epoch t_i.

2.3 Data Processing Procedure and Test Setup

A flowchart of the data processing procedure for the proposed approach for combining low- and high-frequency clock solutions is presented in Figure 2. The procedure consists of two parts: the first part focuses on data and products, and the second part is related to data processing. First, the required real-time GNSS products, onboard GNSS observations, and other necessary model information are collected. Then, errors, such as relativistic errors, phase center offsets (PCOs), and phase center variations (PCVs), can be corrected according to the models provided by Kouba and Héroux (2001). Subsequently, the short-term predicted LEO satellite clocks from a low-frequency BLS solution are used to update the O-C terms with appropriate constraints on the LEO satellite clocks in a high-frequency solution. Subsequently, the Kalman-filter-based processing method (high-frequency solution) is utilized to estimate all unknown parameters, including the real-time LEO satellite orbits, clocks, and float ambiguities. The LEO satellite clocks and orbits at the initial epoch are estimated via the single point position approach.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 2

Flowchart of the data processing procedure for LEO satellite clock estimation using a combination of low- and high-frequency clock solutions; DCB: differential code bias; EOP: Earth orientation parameter

The timeline for applying the combined low- and high-frequency clock solutions is shown in Figure 3. Different sessions of the BLS POD are initiated every 3 min, with the processing starting time for each round shifted by 3 min. This approach ensures the continuity of predicted clocks within the same prediction window (e.g., 0–3 min) across consecutive processing rounds, facilitating their continuous integration into the high-frequency clock solution. Depending on the processing time T_pro of the low-frequency solution, the predicted clocks of the prediction window [T_pro, T_pro + ∆T] can be introduced into the high-frequency clock estimation, where ∆T denotes the time interval between two subsequent sessions of the low-frequency clock estimation, which is 3 min in this case, as mentioned before. In this study, different values of T_pro (0, 3, 6, 9, and 12 min) are tested, leading to predicted clock constraints for prediction periods of 0–3, 3–6, 6–9, 9–12, and 12–15 min.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 3

Timeline of the combined low- and high-frequency LEO satellite clock determination

To validate the proposed approach, GPS dual-frequency phase and code observations tracked onboard Sentinel-3B with an orbital altitude of approximately 810 km (Li et al., 2022) were used. The details of the processing strategies are given in Table 1. As mentioned earlier, BLS sessions shifted by 3 min each were performed using data from day of year (DOY) 226 to 231, 2018. The estimated LEO satellite clocks were predicted for 15 min in each session, generating clocks within prediction periods of 0–3, 3–6, 6–9, 9–12, and 12–15 min (see Table 2). Real-time GPS satellite products were provided by the National Centre for Space Studies (CNES) (Laurichesse et al., 2013) in France. To avoid the stability problem encountered for real-time time references, the real-time GNSS clocks (here, the CNES products) are re-referenced to a stable GPS satellite clock (Wang et al., 2024). Thus, the time reference stability in this study is expected to be better than that of the LEO clock estimates.

For filter-based high-frequency LEO satellite clock estimation in real time, predicted clocks from different prediction periods (see Table 2) are introduced and constrained with different variances, which will be introduced in Section 3. Note that the duration of the prediction periods was determined based on the processing time of the BLS reduced-dynamic POD. For processors operating at 3.0 GHz or higher, the processing time amounts to approximately 3 min. Therefore, five prediction windows of 3 min each were set to cover a total prediction range of 3–15 min (see Table 2). Detailed processing information for the high-frequency solutions is presented in Table 2. Ten strategies for different constraining variances (K1–K10) were designed for the introduced clocks of each prediction window for comparison and analysis purposes (see Table 3). K1 denotes the case in which no external clocks are introduced. K2 to K10 refer to scenarios with decreasing constraint strength.

3.1 Precision of Predicted LEO Satellite Clocks

The precision of the introduced clocks for five different prediction periods, i.e., 0–3, 3–6, 6–9, 9–12, and 12–15 min, was assessed. Reference clocks were obtained by post-processing the LEO satellite clocks using the final GNSS products from the Center for Orbit Determination (CODE) in Europe (Dach et al., 2013). The clock errors in this study are the real-time processing errors, with expectation values equal to zero, which are obtained via a traditional double-differencing procedure, i.e., first removing the reference differences of the CNES and CODE clock products and then comparing the re-referenced real-time LEO satellite clocks with the final clocks.

Figure 4 shows the connected predicted clock errors for different prediction windows mentioned above on DOY 226, 2018. The precision of the predicted clocks degrades with increasing prediction time, giving STDs of 0.12, 0.20, 0.32, 0.43, and 0.56 ns, respectively, on the test day. Table 4 lists the average STDs of predicted clock errors for different prediction periods from DOY 226 to 231 in 2018. The spikes are related to the large biases in the predicted clock errors, particularly for increasing prediction time.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 4

Predicted clock errors for prediction periods of 0–3, 3–6, 6–9, 9–12, and 12–15 min on DOY 226, 2018

The performances of the predicted clocks were also assessed using the MDEVs. As shown in Figure 5, the short-term stability decreases with increasing prediction time, particularly for an averaging time of 1000 s or less. MDEV values for averaging times of 500, 1000, and 2000 s are listed in Table 5 for different prediction periods.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 5

MDEVs of predicted clock errors for prediction periods of 0–3, 3–6, 6–9, 9–12, and 12–15 min on DOY 226, 2018

For future dedicated LEO PNT broadcasting GNSS-like navigation signals, onboard LEO satellite clocks are likely to be good clocks. Examples include the USOs on CentiSpace (Chang et al., 2025) and even atomic clocks incorporated in other LEO navigation constellations. For LEO satellites with low-stability clocks, the prediction precision of low-frequency clocks is degraded, and the results in this section are not appliable.

3.2 Real-Time Clock Performance with the Incorporation of BLS-Based Predicted Clocks

As mentioned before, predicted clocks of different prediction periods are constrained in the real-time filter-based clock determination. Filter-based clock errors applying different constraints for the introduced clocks are shown for different prediction windows for DOY 226, 2018, in Figure 6. In each subplot, the scenario presenting the best solution is indicated. Taking the prediction period of 0–3 min (top left) as an example, a long convergence time and some spikes can be observed for K1 (no constraint) owing to poor geometry, yet a significantly shortened convergence time and improved stability can be observed for K5, exhibiting a 55.3% improvement in precision compared with K1.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 6

Filter-based clock errors applying different constraining strategies, with the introduction of BLS-based predicted clocks for 0–3 min (top left), 3–6 min (top right), 6–9 min (middle left), 9–12 min (middle right), and 12–15 min (bottom) for DOY 226, 2018

Figure 7 shows the average improvements of the real-time clocks applying different strategies. Cases that do not show improvement are not illustrated in this plot. For the prediction period of 0–3 (purple bars) and 3–6 min (green bars), introducing the predicted clocks applying all of the tested constraining strategies (K2-K10) results in significant improvement, particularly for scenarios K5–K7. For a prediction window of 6–9 min (blue bars), only the weak constraints of K5–K10 improve the real-time LEO clock performance. For prediction periods longer than 9 min (red and yellow bars), one must apply weak clock constraints to obtain improvement.

The STDs and percentage improvements for all strategies are listed in Tables 6 and 7. When introducing clocks from the five prediction windows (0–3, 3–6, 6–9, 9–12, and 12–15 min), the best constraining strategies are shown to be K6, K8, K8, K9, and K9, leading to improvements of 46.38%, 39.33%, 31.77%, 25.62%, and 18.35%, respectively. The STD of the real-time clock errors is reduced from 0.24 ns without a constraint to 0.13, 0.15, 0.16, 0.18, and 0.20 ns for the five prediction windows applying the corresponding best constraining strategies. It can be seen that as the prediction time of the introduced clock increases, the proper constraints to be applied become weaker, and the improvement is smaller. This trend suggests that the time used for low-frequency BLS processing (T_pro) (see Figure 2) is a main factor in improving the real-time clocks.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 7

Average improvement of real-time filter-based clocks applying different constraints compared with K1 (no constraint) for DOY 226 to 231, 2018

3.3 Discussion of SISRE

In real-time PNT service, the orbits and clocks are typically used in a combined form; consequently, the SISRE is of higher concern than the individual precisions. When LEO observations are used for PNT, the LEO satellite SISRE can be expressed as follows (Heng et al., 2011):

$SISRE=\sqrt{{({\omega }_R\cdot \Delta R-c\cdot dclk)}^2{\omega }_{SW}^2(\Delta A^2+\Delta C^2)}$11

where ω_R and ω_sw are projection coefficients, which depend on the orbital heights (Montenbruck et al., 2018). dclk represents the LEO satellite clock errors with the daily mean value removed to address the clock precision (STD) instead of the RMS. ∆R, ∆A, and ∆C correspond to the LEO satellite orbital errors in the radial, along-track, and cross-track components, respectively. In this subsection, to assess changes in the real-time orbit/clock correlations, the LEO satellite orbits estimated together with the filter-based LEO satellite clocks are used for calculating the SISREs. The precise orbits of the Sentinel-3B satellite are provided by the European Space Operations Center (Montenbruck et al., 2021). In the following, the RMS of the SISRE in Equation (11) is directly denoted as the SISRE.

The SISRE without a clock constraint is 0.09 m, and the average improvements in the SISRE for different constraint strategies are illustrated in Figure 8. For a prediction window of 0–3 min, introducing predicted clocks applying all of the tested strategies leads to improvements compared with K1 (no constraint), particularly for K2–K7. For each strategy, the improvement with respect to K1 decreases with increasing prediction time. When introducing predicted clocks with a prediction time greater than 9 min (red and yellow bars), only slight improvements can be achieved by applying weak constraints (K8–K10). As shown by the SISREs and improvements in Tables 8 and 9, the best constraining strategies for the five prediction windows are K4, K7, K8, K9, and K9, bringing improvements of 34.73%, 26.48%, 18.58%, 13.55%, and 8.93%, respectively. The SISREs are accordingly reduced from 0.09 m to 0.06, 0.07, 0.07, 0.08, and 0.08 m, respectively.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 8

Average improvements in the SISRE when introducing predicted clocks applying different constraining strategies for five different prediction periods from DOY 226 to 231, 2018

3.4 Introducing Filtering-Based Predicted Clocks in Real-Time Clock Estimation

As shown in Section 3.3, shortening the prediction time of the introduced clocks is essential for improving the real-time clock precision. Filter-based clock processing is more efficient than BLS clock processing and can also be predicted and introduced again into the real-time filter-based clock determination. In this subsection, as an extension of the previous results, the real-time clock performances are compared when predicted clocks of 6–9 min based on BLS clock determination and those of 0–3 min based on filter-based clock determination are introduced. The usage of different prediction times for the introduced clocks reflects the different efficiencies of the BLS and filter-based processing. The latter case faces shorter delays and thus requires a shorter prediction time to be introduced into real-time clock processing.

Figure 9 shows the clock errors (left) and MDEVs (right) that arise when the two types of predicted clocks are introduced. The short-term noise based on BLS adjustment (blue lines) is higher than that of the filter-based candidate (red lines), but larger systematic effects can be observed in the filter-based cases (red lines), leading to larger MDEVs in the mid- to long-term.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 9

Predicted clock errors (left) and MDEV (right) for a prediction window of 6–9 min based on BLS adjustment and for a prediction window of 0–3 min based on a Kalman filter for DOY 226, 2018

To further demonstrate the effects that arise when the above two types of predicted clocks are introduced, Figure 10 shows real-time clock errors without a constraint (green line), errors for Kalman-filter-based predicted clocks with the K6 constraint (red line), and errors for BLS predicted clocks with the K8 constraint (blue line). While the BLS-predicted clocks help with the stability of the real-time clocks, introducing Kalman-filter-based predicted clocks does not help reduce systematic effects contained in the real-time clocks. In contrast, the stability becomes worse in this case, owing to the prediction loss of the introduced Kalman-filter-based predicted clocks.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 10

Real-time clock errors when predicted clocks based on different approaches are introduced for DOY 226, 2018

The STDs of the real-time clock errors agree with the conclusions above (see Figure 11). For instance, introducing filter-based predicted clocks does not improve the filter-based clock estimation itself; instead, it degrades the solution precision.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 11

STDs of clock errors when predicted clocks based on two different approaches are introduced for DOY 226, 2018