Resilient timekeeping algorithm with multi-observation fusion Kalman filter

The timescales incorporated into the Primary Frequency Standard (PFS) exhibit excellent stability and accuracy. However, during the dead time of PFS, the reliability of the timescale can be compromised. To address this issue, a resilient timekeeping algorithm with a Multi-observation Fusion Kalman Filter (MFKF) is proposed. This algorithm fuses the frequency measurements from hydrogen masers with various reference frequency standards, including PFS and commercial cesium beam atomic clocks. The simulation results show that the time deviation and instability of the timescale generated by MFKF are improved compared to those with Kalman filtering. The experimental results demonstrate that even within 70 days of PFS dead time the resilient timescale generated by MFKF can operate reliably. Furthermore, it is theoretically proven that MFKF produces a smaller post-covariance than that with single-observation Kalman filtering.

By incorporating the Cs fountain clock (Primary Frequency Standard, PFS) into the generation of timescale, its autonomy, accuracy, stability, and other characteristics can be easily improved (Levine, 1997;Wang et al., 2020Wang et al., , 2022)).Moreover, the emergence of optical clocks (Barbiero et al., 2022;Poli et al., 2013) as the next generation of PFS promotes the research on timekeeping methods for incorporating PFS (Formichella et al., 2021;Galleani et al., 2019Galleani et al., , 2020;;Grebing et al., 2016;Hachisu et al., 2018;Yao et al., 2017;Yao et al., 2019).In general, the frequency of the Hydrogen Maser (HM) is estimated using PFS and steered to generate a timescale.However, PFS is a laboratory instrument and usually cannot operate continuously for an extended period like commercial clocks.Therefore, maintaining the performance of the timescale when the PFS is not operational presents an urgent issue.
Multiple frequency standards are used as redundant backups to ensure a stable time scale during PFS failures.For example, Sytèmes de Référence Temps Espace of Laboratoire National de Métrologie et d'Essais (LNE-SYRTE) (Rovera et al., 2016) steers HMs using four PFS to generate Coordinated Universal Time (UTC) in France generated at Observatoire de Paris (OP).Furthermore, Physikalisch-Technische Bundesanstalt (PTB) (Bauch et al., 2012) and Italian National Institute of Metrological Research (INRiM) (Galleani et al., 2020) employ frequency standards like PFS and commercial Cs clocks to steer HMs, resulting in the generation of continuous and reliable timescales UTC (PTB) and Italian time scale UTC (IT).Although these algorithms utilize these standards to generate reliable timescales, they lack the integrated performance superiority seen with ensemble timescale algorithms involving multiple atomic clocks.Unfortunately, continuous operation over a long period is difficult for most PFS, which the traditional ensemble timescale algorithms require.Therefore, we propose a resilient timekeeping algorithm utilizing a Multi-observation Fusion Kalman Filter (MFKF).This algorithm resiliently fuses the HM frequency measurements using various standards such as PFS and commercial Cs beam clocks.The filtering results are then used to steer the HM and generate timescale TA_MFKF.TA_MFKF can utilize the measurements with other frequency standards, such as commercial Cs beam clocks, during PFS outages to ensure reliable and continuous operation.
Additionally, this algorithm that fuses multiple frequency standards exhibits a superior performance compared to the one utilizing only a single frequency standard.The simulation results show that TA_MFKF exhibits smaller time deviations and better stability than the Kalman Filter (KF) with a single observation.In addition, the time difference between the rapid realization of the Coordinated Universal Time (UTCr) (Petit et al., 2014) and TA_MFKF is indirectly obtained by utilizing the optical fiber link between the Beijing Satellite Navigation Center (BSNC) and National Institute of Metrology (NIM) (Wang et al., 2020(Wang et al., , 2022)).The comparison results show that the TA_MFKF can maintain reliable operation within 70 days of PFS's dead time.
Section "Measurement system" presents the frequency measurement systems utilized for Cs fountain clocks, Cs beam atomic clocks, and HMs.Section "Theory" introduces the theory of timescale based on the MFKF algorithm and theoretically proves that MFKF improves the performance of the fused timescale compared to the KF with a single measurement source.Section "Simulation results and Discussions" simulates two PFSs and one HM to compare the performance of MFKF and KF during PFS regular operation.In the fifth section, one PFS, five commercial Cs beam clocks, and one HM of BSNC are resiliently integrated using MFKF to conduct resilient timescale and then compared with single-observation KF.Finally, the article is concluded with prospects.

Measurement system
Cs fountain clocks cannot generate time signals directly like commercial atomic clocks.Instead, they require a frequency signal from another atomic clock to serve as the frequency reference, which is then processed to estimate the frequency compared to that atomic clock.Therefore, a measurement system has been devised to integrate Cs fountain clocks into the conventional atomic clock comparison measurement system.
Figure 1 displays the measurement principle and schematic diagram of the HM, Cs fountain clock, and Cs beam clocks.The HM functions as the master clock for the Cs fountain clock, providing a 5 MHz frequency signal as its reference.Each day, a frequency deviation arises between the Cs fountain clock and the HM.Moreover, HM's 1 PPS (1 Pulse Per Second) signal is utilized as an input in the time interval counter, enabling the measurement of clock differences with other atomic clocks.The clock deviation signal is processed daily through software.All frequency comparison data and time comparison data are stored in a database.At a fixed time, the data processing software automatically removes the frequency shift of the Cs fountain clock, determining the frequency deviation between the HM and the Cs fountain.Simultaneously, the data processing software utilizes the Least Squares (LS) method to fit the phase comparison data on the previous day, deriving the daily frequency deviation between the HM and other atomic clocks.
The data processing software automatically executes initial data processing at a set time every day, effectively eliminating the frequency offset of the Cs fountain clock.Following that, the software calculates the discrete frequency deviation between the master clock H0 and Cs fountain: where f HM and f CsF1 are the frequency of the master clock and Cs fountain clock, respectively.In calculating the frequency difference between Cs beam clocks ( Cs − i ) and HM, the LS method is applied to the time difference data to estimate the frequency deviations between the master clock H0 and Cs beam clocks, denoted as f i : The frequency differences ( f CsF1 -HM and f i ) of Cs fountain and Cs beam clocks with respect to HM are utilized as the data source of MFKF.

Theory
Most Cs fountain clocks cannot operate continuously for a long period of time like commercial atomic clocks and have different measurement mechanisms, which makes it hard to compare their time difference with the atomic clock ensemble directly.As a result, it is difficult for the (1) traditional ensemble clock timescale algorithms to incorporate Cs fountains into the generation of time scales.By utilizing MFKF, it is possible to combine Cs fountain clocks with other frequency standards and achieve resilient timekeeping.The frequency estimation and prediction of the HM are first calculated by using MFKF, and then the frequency of the HM is steered by using the predicted value to generate the timescale.The scheme of the resilient timekeeping algorithm based on MFKF is shown in Fig. 2.
In this section, we begin by introducing the concept of multi-observation fusion Kalman filtering (Gan & Harris, 2001), followed by the development of a dynamic model for HMs and an observation model for frequency standards such as Cs fountain clocks and Cs beam clocks.Based on these models, we present an innovative and resilient timekeeping algorithm with the MFKF method.Finally, the estimation error is analyzed.

Multi-observation fusion Kalman filter
The model of MFKF (Gan & Harris, 2001) describes the tracking of a target by N sensors: where t represents the time index, x(t) represents the state vector, A(t) denotes the state transition matrix, B(t) denotes the input control matrix, u(t) denotes the input vector, y j (t) and C j (t) are measurement vectors and measurement matrices, respectively, and v(t) and w j (t) are noise vectors of dynamic target and measurement with the covariance matrices Q(t) and R j (t) , respectively.In this article, the dynamic target is a HM, while the sensors are Cs fountain clocks, Cs beam clocks, and other frequency standards. (3)

CsF1
Multi-channel interval counter The difference between MFKF and the classical Kalman filtering is that MFKF's measurement equation is extended from one to N. MFKF enables the integration of N measurement models into a single unified measurement model (Gan & Harris, 2001): where w(t) ∼ N (0, R(t)) .y(t) , C(t) , and w(t) are the fusion of N measurement equations: The Kalman filter can obtain an optimal estimate for the model described in (3) and ( 5).The solution is given by Eqs. ( 9) to ( 13).
where the superscript '^' denotes the estimation of that variable, and P(t|t) is the post-covariance matrix.(Gan & Harris, 2001).
In the usual Kalman ensemble timescale algorithm (Galleani & Tavella, 2010), the state vector, x(t) , repre- sents the estimated clock difference of other clocks relative to the reference clock, while the measurement vector represents the actual clock difference measurements obtained with the measuring instruments.Our proposed algorithm uses the state vector to represent the frequency and the drift of the HM itself, and the measurement ( 5) vector is the result of the frequency standard's estimation of the HM.Thus, to implement this algorithm, it is imperative to develop a dynamic model (Eq.3) for the HM and a measurement model (Eq.4) for the frequency standards.

Dynamic model
In the dynamic model, Eq. (3), the state vector is where f (t) and d(t) are the frequency and drift of the HM, respectively.Then, the state transition matrix is where τ 0 is the sampling period, which is the time differ- ence of t and t − 1 .The noise covariance matrix of atomic clocks (Galleani & Tavella, 2010;Zucca & Tavella) is written as: where q f and q d are the intensity of the atomic clock's White Frequency Modulation (WFM) noise and Random Walk Frequency Modulation (RWFM) noise, respectively.They can be calculated from the Allan variance, and the formula (Galleani & Tavella, 2010;Zucca & Tavella) is written as: Since the dynamic model in this paper are HMs, the noise covariance matrix of the clock can be simplified to:

Measurement model
In the measurement model, Eq. (4), the measurement vector represents the frequency of the HM measured by the frequency standards.Therefore, the measurement matrices C j (t) in measurement Eq. ( 4) is written as: ( 14) At this point, the vector R j (t) degenerates into a sca- lar R j (t) .The sampling time interval in the measurement system is one day, so when sampling over this period, the noise in usual frequency standards, such as Cs fountain clocks and Cs beam clocks, is dominated by WFM.Therefore, the measurement covariance R j (t) is given by: where q j is the WFM intensity of the frequency stand- ards.In Cs beam atomic clocks, q j is also known as q f and can be calculated using Eq. ( 17).
The dominant noise type in Cs fountain clocks is WFM.Moreover, the stability of Cs fountain is written as (Santarelli et al., 1999): For the meanings of each symbol in ( 21), refer to (Wang et al., 2020).Therefore, RWFM of the PFS can be ignored, so the formula of q j related to Allan variance is written as: q j of PFS can be expressed by:

Resilient timekeeping algorithm based on MFKF
MFKF is a model that measures one dynamic target using multiple sensors, and N frequency standards, such as PFS and commercial Cs clocks, measure the frequency of the HM.Therefore, the MFKF model for this resilient timekeeping algorithm is as follows: where v(t) ∼ N (0, Q(t)) , w j (t) ∼ N (0, R j (t)) .The spe- cific meanings of each expression are explained in detail in Sections "Dynamic Model and Measurement Model".
To avoid the impact of specific frequency standard outages on fusion results, a resilient factor ρ j (t) is added to the MFKF, and then the fusion equations are: where where ρ j (t) = 1 when the j-th frequency standard is avail- able and ρ j (t) = 0 when it is unavailable.The incorpo- ration of resilient factors enables the adaptive update of the weights of frequency standards and observed noise variance in the circumstances where specific frequency standards become unavailable, thereby reducing the impact of interruptions to some frequency standards on the frequency estimation of the HM.
By utilizing the five solution equations to, we can obtain the estimated and predicted frequency of the HM through equations and.By steering the HM with the predicted frequency, the resilient timescale TA_ MFKF can ultimately be obtained.

Estimation error analysis
This section gives the post-covariance of the estimates provided by MFKF and compares it with that of KF with a single measurement source.The post-covariance matrices of KF and MFKF have the same formula expressed as Eq. ( 13).Let the post-covariance matrix of timescale provided by MFKF be P(t|t) , and that pro- vided by Kalman filter based on only one of the observation sources be P i (t|t) .Substituting (10) into (13), we get Therefore, the estimated covariance matrix is correlated with the previous round and iteratively propagated.Assuming that the covariance matrix of the estimated timescale at the previous round by using MFKF is the same as that of the classical Kalman filter, namely ( 26) the estimated error of the two algorithms can be compared.Moreover, suppose that MFKF and Kalman filter observe the same object, so Therefore, the difference between P(t|t) and P i (t|t) is Substitute ( 11) into (34), then where Therefore, the magnitude of P(t|t) and P i (t|t) is only related to R(t) and R i (t).
Due to Eq. ( 8), R(t) < R i (t) , thus P(t|t) is smaller than P i (t|t) , there- fore, it is proved that MFKF has a smaller post-covariance than KF with a single measurement source.

Simulation results and discussions
This section uses simulation experiments to demonstrate that the timescale generated by the MFKF algorithm using multiple PFS outperforms that generated by the KF algorithm using one PFS.Section 3 has theoretically proven that the MFKF algorithm has a smaller post-covariance than the KF algorithm.However, it is difficult to verify in our laboratory as we only have one Cs fountain clock and several commercial Cs beam atomic (31) clocks.The noise of the Cs beam atomic clocks is an order of magnitude higher than that of the Cs fountain clock, and therefore their fusion with the Cs fountain clock does not significantly improve performance, as shown in the next section.To address this limitation, two Cs fountain clocks with a similar performance and one HM were simulated, and the MFKF and KF algorithms were used to generate timescales, respectively.The time differences and stability analysis were compared with an ideal reference.
Based on the performance of Cs fountain clocks and HMs in our laboratory, the noise parameters of two Cs fountain clocks and one HM were determined and presented in Table 1.
Using Stable32 (Riley, 2020), several clocks' frequency and time difference data were simulated relative to an ideal reference.The simulation had a sampling interval of one day and lasted for 500 days.Figure 3 presents the stability analysis of the simulated clocks compared to the ideal noiseless reference.
By subtracting the data of the HM and the Cs fountain, we obtained their differences, namely HM1_CsF1 and HM1_CsF2.KF was then used to filter both HM1_CsF1 and HM1_CsF2 separately, and the filtering result steered the HM.This process led to the creation of two timescales, TA_KF_SIM_1 and TA_KF_SIM_2.With MFKF, HM1_CsF1 and HM1_CsF2 data were fused and filtered, generating TA_MFKF_SIM by steering the HM with the filter result.Figures 4 and 5 present the three timescales' time differences and stability relative to the ideal reference, respectively.In Fig. 4, the maximum time difference of 12.0 ns between TA_KF_SIM_1 and TA_KF_SIM_2 is observed on the 490th day.The corresponding average frequency deviation ( f d ) of TA_KF_SIM_1 and TA_KF_SIM_2 within 490 days is 2.83 × 10 -16 .Based on the WFM of CsF1_SIM and CsF2_SIM, their stability for the averaging times of 500 days is calculated to be 1.08 × 10 -16 ( σ 1 ) and 1.11 × 10 -16 ( σ 2 ), respectively.The stability of the deviation between CsF1_SIM and CsF2_SIM is 1.55 × 10 -16 ( σ c ).The ratio of f d to σ c is 1.83, indicating that the 12.0 ns time deviation of the two falls within the normal range.
The Root Mean Square (RMS) of TA_MFKF_SIM, TA_KF_SIM_1, and TA_KF_SIM_2 is 0.84, 2.99, and 2.67 ns, respectively.As shown in Fig. 5, the stability of TA_MFKF_SIM is either superior or comparable to TA_KF_SIM_1 and TA_KF_SIM_2 at nearly all averaging times.As indicated by the simulation results, the MFKF algorithm exhibits a better performance in terms of time deviation and stability compared to the traditional KF algorithm.
Therefore, by utilizing MFKF to fuse multiple frequency standards, a better timescale is obtained compared to that using a single frequency standard alone.

Experiment results and analysis
In this section, the Cs fountain clock CsF1 (Shi et al., 2019;Wang et al., 2020Wang et al., , 2022) ) and five commercial Cs beam clocks in BSNC were integrated using MFKF to predict the HM's frequency, and then steer the HM to generate the resilient timescale TA_MFKF.The time difference between UTCr and TA_MFKF was indirectly obtained by the optical fiber link between BSNC and NIM.The time stability of the optical fiber link is less than 6.0 ps/s and 0.9 ps/100 s. (Liang et al., 2015).
The experimental data comprised the frequency comparison data obtained with Cs fountain clock CsF1, five commercial Cs beam clocks, and one HM for the period from MJD58284 to MJD58448 (165 days).The observed HM's daily stability was 8.5 × 10 -16 , and its frequency drift was 4.21 × 10 -17 d −1 .The noise parameters of frequency standards CsF1 and five Cs beam clocks are shown in Table 2. Furthermore, the frequency of Cs beam clocks was calibrated using CsF1.
To evaluate the performance of TA_MFKF in the event of a failure of CsF1, a data interruption of CsF1 was simulated during the period MJD58354-58424 (70 days).MFKF was used to resiliently generate the timescale TA_ MFKF using the data from both Cs beam clocks and the interrupted CsF1.Furthermore, another timescale, TA_ PFS, was generated by the Kalman filter using the data   from a continuous CsF1.Both TA_MFKF and TA_PFS were compared and evaluated for their performance during CsF1 interruption and its normal operation.
Figure 6 shows the measured and filtered values obtained with the two methods.The red star represents the fusion measurement result of one discontinuous Cs fountain clock and five commercial Cs clocks based on MFKF, and the red curve is its filter result.The blue dots represent the measurement result of the continuous Cs fountain clock, and the blue curve is its filter result based on the Kalman filter.Measurement_MFKF represents the HM's fusion measurement with reference to CsF1 and five Cs beam clocks, including a simulated 70-day CsF1 interruption.Measurement_PFS represents the measurement of the HM with reference to continuous CsF1.Estimate_MFKF and Estimate_PFS represent the filtering results obtained with MFKF and KF, respectively.The frequency difference between Estimate_MFKF and Esti-mate_PFS is shown in Fig. 7.
Two timescales, TA_MFKF and TA_PFS, were generated by steering the HM using its frequency estimates, expressed as Estimate_MFKF and Estimate_PFS, respectively.Figure 8 illustrates the indirect comparison of the time difference between UTCr and the two timescales.This comparison is made by utilizing an optical fiber link to UTC(NIM) and referencing the downloaded time difference data between UTCr and UTC(NIM) from the BIPM FTP server.TA_MFKF and TA_PFS were generated using the interrupted and continuous data of Cs Fountain, respectively.
Figures 6 and 7 show that when the Cs fountain clock is available (MJD58284-58353), Estimate_MFKF is nearly identical to Estimate_PFS, with a maximum frequency difference of only 3.61 × 10 -16 .This is because the measurement noise of commercial Cs clocks is approximately ten times higher than that of Cs fountains.As a result, the weight assigned to commercial Cs clocks tends toward zero, while the weight assigned to the Cs fountain clock approaches 100% in Eq. during observation fusion.Consequently, when the Cs fountain clock operates normally, the measurements and estimates obtained from the fusion observation and the Cs fountain clock are nearly identical.
Figures 7 and 8 show that during the 70-day period when the Cs fountain clock is unavailable (MJD58354-58424), the frequency difference between Estimate_ MFKF and Estimate_PFS remains within ± 1 × 10 -14 , while the time deviation between TA_MFKF and UTCr within ± 5 ns.When the Cs fountain clock is in regular operation, the frequency difference between Estimate_ MFKF and Estimate_PFS immediately is reduced to 1.96 × 10 -15 , and their maximum frequency difference is reduced to 3.78 × 10 -16 after five days.
The classic weighted average algorithm introduces a prediction term to prevent phase and frequency jumps when some clocks fail.However, the MFKF approach does not require additional operations and avoids jumps.Due to that, Kalman filtering performs a weighted average of the measured value and the predicted value, where the Kalman gain determines the weight.As shown in Eq. ( 11), even if the clock with the best performance (the Cs fountain clock in this experiment) experiences an interruption, the covariance of observation noise will increase, which causes the Kalman gain to decrease.As The Cs clock ensemble's role was demonstrated by conducting an experiment during the 70-day interruption period of the Cs fountain clock, where the HM was allowed to run freely without being steered by the Cs clock ensemble.The time scale generated by this experiment is referred to as TA_Free.Figure 9 shows that there is a time deviation of more than 1100 ns between TA_ Free and UTCr during the 70-day interruption of CsF1.The reason is that there is the frequency drift and bias (usually at the level of 1 × 10 -13 ) of the HM, which causes a huge time deviation when it runs freely and without being steered by the Cs clock ensemble, even though the short-term stability of the HM is good.Therefore, during the interruption of the Cs fountain clock, if the proposed resilient timescale algorithm is applied to fuse the Cs clock ensemble and steer the HM, the performance of the time scale can effectively be improved.
Furthermore, the stability of both TA_MFKF and TA_Free relative to UTCr were computed separately and depicted in Fig. 10.The red and green curves represent the stability of TA_MFKF and TA_Free.
As shown in Fig. 10, TA_MFKF exhibits much better stability than TA_Free.When the Cs fountain clock is interrupted, TA_Free relies on a free-running HM, which has frequency drift and bias.As a result, TA_Free does not receive the correction from other frequency standards, leading to frequency drift and deviation.On the other hand, during the Cs fountain clock interruption period, TA_MFKF utilized five Cs beam atomic clocks as frequency standards, which were fused to steer the HM and correct its frequency drift and deviation.This difference in the fusion of frequency standards explains why TA_MFKF demonstrates much better stability than TA_Free.When the Cs fountain is interrupted, the predicted data from the previous day can also be utilized as a fixed reference to steer the HM during the interruption period, resulting in a time scale TA_Prediction.As shown in Fig. 11, TA_Prediction follows the trend before the interruption, while the HM experiences a small frequency deviation.Due to the lack of fused frequency standards to steer the HM, the resulting time deviation (9.7 ns) is greater than that of TA_MFKF and TA_PFS.
Therefore, the experiments verify the TA_MFKF's resilient timekeeping ability, whether the Cs fountain clock is in regular operation or not.

Conclusion
A resilient timekeeping algorithm with a multi-observation fusion Kalman filter is proposed.It integrates the frequency measurements of a HM with several frequency standards, including one discontinuous Cs fountain clock and five commercial Cs clocks, to generate the timescale (TA_MFKF).During the 70-day interruption of the fountain clock, the time deviation between TA_MFKF and UTCr is maintained within ± 5 ns.In contrast to the traditional Kalman filtering algorithms, which fail when reference clock data is interrupted, the MFKF algorithm can resiliently adjust the observation noise covariance matrix and fuse the measurement data from other continuous frequency standards, thus maintaining reliable operation of the timescale.
Furthermore, the simulations have shown that the timescale generated by MFKF is more stable than that of KF due to the resilient fusion of multiple PFS, which reduces the impact of measurement noise on HM frequency estimates.Moreover, externally compared data can also be incorporated into the resilient timekeeping algorithm as a measurement reference source for the institutions that do not emphasize autonomy.

Fig. 3
Fig. 3 Stability of simulated atomic clocks relative to the ideal reference

Fig. 4
Fig. 4 Time deviation with reference to the ideal reference of TA_ MFKF_SIM, TA_KF_SIM_1, and TA_KF_SIM_2, respectively

Fig. 6 Fig. 7
Fig. 6 Measurement and estimation of the frequency difference between frequency standards and hydrogen maser with MFKF and KF, respectively

Table 1
Noise parameters of simulated clocks a HM1_SIM's initial frequency offset is 2.0 × 10 -13 b White Phase Modulation noise

Table 2
Noise parameters of frequency standards