Sage Journals: Discover world-class research

Abstract

A new method based on the innovation of Kalman filter and extrapolation-accumulation is proposed to detect weak frequency anomalies in atomic clocks with short detection time and high detection probability. In this method, the detection statistics of the innovation extrapolation method in several epochs are accumulated. It avoids the influence of outliers, increases the noncentrality parameter of chi-square distribution, and realizes more effective detection of weak frequency anomalies. The simulation results show that compared with the innovation method and the innovation extrapolation method, the new method has a higher detection probability for micro frequency anomalies, and the detection time is shortened. The new method is used to analyze the real data of the cesium atomic clock, and the results are consistent with the simulation results.

Keywords

Frequency anomaly atomic clock Kalman filter extrapolation-accumulation chi-square statistics

Introduction

Atomic clocks are an important part of the Global Navigation Satellite System (GNSS). They are used to generate standard time and frequency signals in both ground stations and satellites so that users can obtain accurate position and time information. It can be found that different numbers of hydrogen, cesium and rubidium atomic clocks have been equipped on the ground stations and satellites of GPS, GLONASS, Beidou, and Galileo.^1–4 Even though atomic clocks are ultra-stable oscillators, they are subject to temperature, magnetic fields, radiation and aging, which can cause frequency anomalies. The frequency anomalies of atomic clocks mainly include frequency jump and frequency drift jump. If these anomalies are not detected correctly and timely, it will result in thousands of meters of range error.⁵ In Figure 1, the stability of a cesium atomic clock is simulated when it is subjected to a small frequency jump and a frequency drift jump. Stability is deteriorating, and that is not acceptable. Therefore, atomic clock frequency anomaly detection is a key technique in GNSS and time-keeping.

Figure 1.

The influence of weak frequency anomalies on the stability of a cesium atomic clock.

In the past few years, some effective methods for detecting frequency anomalies have been proposed. In 2008, Riley⁶ proposed some simple frequency anomaly detection methods of BLKAVG, SEQAVG, and CUSUM. Nunzi and Carbone,⁷ Galleani,⁸ and Galleani and Tavella⁹ worked to detect the frequency anomalies with generalized likelihood ratio test (GLRT) and dynamic Allan variance (DAVAR), which could not only detect the anomalies but also determine the types of anomalies. In recent years, Galleani and Tavella¹⁰ and Song et al.¹¹ proposed innovation based on Kalman, an adaptive Kalman algorithm to detect the anomaly of atomic clock frequency, and both of them achieved good detection results. Since 2014, Huang et al.^12–14 have been committed to the detection technology of weak frequency anomalies of atomic clocks and proposed Kalman N-step prediction, parity transformation, and other methods, achieving more significant detection results. However, the detection time is also very important, and none of these papers mentioned the detection time of their algorithms. Higher detection probability and shorter detection time are our goals. We not only analyzed the detection probability but also analyzed the detection time for the first time.

The innovation method based on the Kalman filter has a weak ability to detect small frequency jump and frequency drift jump. When constructing detection statistics, the abnormal frequency value will “pollute” the state estimation of the filter. To improve the detection probability and shorten the detection time, a novel innovation extrapolation-accumulation method based on the Kalman filter is proposed to detect the small frequency jump and frequency drift jump of atomic clocks. In this method, the detection statistics of the innovation extrapolation method in several epochs are accumulated. It avoids the influence of outliers, increases the noncentrality parameter of chi-square distribution, and realizes more effective detection of weak frequency anomalies. The correctness, superiority, and timeliness of the method are verified by simulating the frequency anomaly data of a cesium atomic clock and the measured frequency anomaly data of a cesium atomic clock in our laboratory.

It is worth mentioning that the proposed method is more sensitive to weak frequency jump than other methods under the same conditions. This is good news for some laboratories where atomic clocks work with precision and stability. The new method is not only adapted to detecting anomalies in a local ensemble of clocks but also provides an effective solution to autonomous integrity monitoring of GNSS satellite clocks.

The paper is organized as follows. In section 2, the Kalman filter model of the atomic clock, the innovation extrapolation method and the chi-square detection of multi-epoch accumulation are described in detail. The superiority of the new method is proved mathematically strictly in theory. The implementation of the algorithm is also described. In section 3, the simulation of weak frequency jump and frequency drift jump is carried out. We give the optimal extrapolation number and the optimal accumulation number under certain conditions. In section 4, the new method is applied to the cesium atomic clock frequency data collected in our laboratory. Section 5 is the conclusion.

Innovation extrapolation-accumulation

The Kalman filter model

The three-state clock model is chosen here, which is expressed as follows:

X (k + 1) = Φ \cdot X (k) + ξ (k)

(1)

where $X (k) = [\begin{matrix} x_{1} (k) \\ x_{2} (k) \\ x_{3} (k) \end{matrix}]$ , $Φ = [\begin{matrix} 1 & T_{s} & 0.5 T_{s}^{2} \\ 0 & 1 & T_{s} \\ 0 & 0 & 1 \end{matrix}]$ , $ξ (k) = [\begin{matrix} ξ_{1} (k) \\ ξ_{2} (k) \\ ξ_{3} (k) \end{matrix}]$ . $x_{1} (k)$ is the phase deviation and $x_{2} (k)$ is the frequency deviation, which is defined as the derivative of the phase deviation $x_{1} (k)$ . $x_{3} (k)$ is the drift rate. $T_{s}$ represents the sampling interval. The process noise vector $ξ (k)$ has three components: $ξ_{1} (k)$ represents the random variation due to White Phase Noise (WPN), and similarly $ξ_{2} (k)$ accounts for White Frequency Noise (WFN), $ξ_{3} (k)$ is the random variation due to Random Walk Frequency Noise (RWFN), respectively.¹⁵ $ξ (k)$ is a multivariate Gaussian random vector with zero mean and a covariance of Q, which can be expressed as:

E (ξ (k)) = 0

(2)

\begin{matrix} E (ξ (k) \cdot ξ {(k)}^{T}) = \\ Q = [\begin{matrix} q_{1} T_{s} + \frac{q_{2} T_{s}^{3}}{3} + \frac{q_{3} T_{s}^{5}}{20} & \frac{q_{2} T_{s}^{2}}{2} + \frac{q_{3} T_{s}^{4}}{8} & \frac{q_{3} T_{s}^{3}}{6} \\ \frac{q_{2} T_{s}^{2}}{2} + \frac{q_{3} T_{s}^{4}}{8} & q_{2} T_{s} + \frac{q_{3} T_{s}^{3}}{3} & \frac{q_{3} T_{s}^{2}}{2} \\ \frac{q_{3} T_{s}^{3}}{6} & \frac{q_{3} T_{s}^{2}}{2} & q_{3} T_{s} \end{matrix}] \end{matrix}

(3)

where $q_{1}, q_{2}, q_{3}$ can be derived from the noise coefficients of WPN, WFN, and RWFN, respectively.^16,17

The frequency measurement model can be expressed as:

Z (k) = H \cdot X (k) + w (k)

(4)

where $H = [\begin{matrix} 0 & 1 & 0 \end{matrix}]$ , $Z (k)$ is the measured frequency deviation and $w (k)$ is White Gaussian Noise (WGN) with variance $R$ .

Innovation extrapolation method

The innovation extrapolation method is an improvement of the innovation method. The innovation extrapolation method is no longer based on state prediction at the previous moment and the observed value at the current moment to get state estimation. Instead, it uses the state estimation before the filter, avoiding the influence of abnormal observations and reducing the “pollution” of the state estimation. The Kalman filter equations are described as follows:

\bar{X} (k + 1) = Φ \hat{X} (k)

(5)

P (k + 1) = Φ \hat{P} (k) Φ' + Q (k)

(6)

\begin{matrix} \bar{K} (k + 1) \\ = P (k + 1) H' (HP (k + 1) H' + R (k + 1))^{- 1} \end{matrix}

(7)

\begin{matrix} \hat{X} (k + 1) = \bar{X} (k + 1) \\ + \bar{K} (k + 1) (Z (k + 1) - H \bar{X} (k + 1)) \end{matrix}

(8)

\hat{P} (k + 1) = P (k + 1) - \bar{K} (k + 1) HP (k + 1)

(9)

where $\bar{X} (k + 1)$ is state prediction vector and $\hat{X} (k + 1)$ is state estimation vector. $P (k + 1)$ represents error covariance matrix and $\hat{P} (k + 1)$ is the error covariance update matrix. $\bar{K} (k + 1)$ is the Kalman gain matrix.

Innovation is defined as:

r (k + 1) = Z (k + 1) - H \bar{X} (k + 1)

(10)

When there is no frequency anomaly, the innovation follows a Gaussian distribution with a mean of zero,¹⁸ which satisfies:

r (k + 1) : N (0, D (k + 1))

(11)

where $D (k + 1)$ is variance, expressed as follows:

D (k + 1) = HP (k + 1) H' + R (k + 1)

(12)

When the frequency anomaly occurs, the mean of the innovation is not zero, which satisfies:

r (k + 1) : N (μ, D (k + 1))

(13)

where $μ$ represents the size of the abnormal. When there is a small frequency jump or a frequency drift jump, if innovation detection method continues to be used, there will be:

\begin{matrix} r (k + 1) = Z (k + 1) - H \bar{X} (k + 1) \\ = Z (k + 1) - H Φ \hat{X} (k) \end{matrix}

(14)

$\hat{X} (k)$ contains the frequency anomaly value at time k, resulting in a small value of innovation, so that the anomaly cannot be found in time. If the state estimation before T of the filter is extrapolated to the current moment, the problem of untimely can be effectively avoided. The specific extrapolation process can be expressed as:

{\bar{X}}^{e} (k + 1) = Φ {\hat{X}}^{e} (k), {\hat{X}}^{e} (0) = \hat{X} (k - T)

(15)

P^{e} (k + 1) = Φ P^{e} (k) Φ' + Q (k), P^{e} (0) = \hat{P} (k - T)

(16)

where ${\hat{X}}^{e} (0)$ is an extrapolation initial value, $P^{e} (0)$ is the first matrix in the extrapolation of the error covariance matrix. The innovation and the variance of innovation generated by extrapolation are:

r (k - T + i) = Z (k - T + i) - H {\bar{X}}^{e} (k - T + i)

(17)

D (k - T + i) = H P^{e} (k - T + i) H' + R (k - T + i)

(18)

Where $1 \leq i \leq T$ . The innovation generated by the extrapolation process is no longer affected by the state prediction of the previous moment and does not include the frequency outliers. Therefore, the innovation is not “polluted.” When the weak frequency anomalies accumulate to a certain degree, the anomaly can be detected effectively and timely by accumulating multiple epoch statistics.

Chi-square detection of multi-epoch accumulation

According to the innovation extrapolation method, a detection statistic is adopted:

λ_{avg} = {r'}_{avg} D_{avg}^{- 1} r_{avg}

(19)

r_{avg} = (D_{avg}^{- 1})^{- 1} \sum_{i = 1}^{T} (D (k - T + i))^{- 1} r (k - T + i)

(20)

D_{avg} = (\sum_{i = 1}^{T} (D (k - T + i))^{- 1})^{- 1}

(21)

To better detect the small frequency jump and the frequency drift jump, a new detection statistic is constructed by accumulating the extrapolated statistic of N epochs:

T_{D} (i) = \sum_{i = k - N + 1}^{k} λ_{avg} (i)

(22)

where $λ_{avg} (i)$ is the i-th extrapolation detection statistic. According to the statistical distribution theory, $λ_{avg}$ follows a chi-square distribution with 1 degree of freedom when no frequency anomaly exists. It has a non-central chi-square distribution with 1 degree of freedom and noncentrality parameter $λ$ when a frequency anomaly occurs.

λ = μ' (D_{avg})^{- 1} μ

(23)

The new statistic $T_{D}$ follows a chi-square distribution with N degrees of freedom when no frequency anomaly exists. It has a non-central chi-square distribution with N degrees of freedom and noncentrality parameter $λ'$ when a frequency anomaly occurs.

λ' = \sum_{i = k - N + 1}^{k} μ' {(D_{avg})}^{- 1} μ \approx N μ' {(D_{avg})}^{- 1} μ

(24)

In a short time, we can approximately assume that the innovation covariance matrix remains constant. The ratio between the new detection statistic and the innovation extrapolation detection statistic is:

\frac{λ'}{λ} \approx \frac{N μ' {(D_{avg})}^{- 1} μ}{μ' {(D_{avg})}^{- 1} μ} = N

(25)

Compared with the innovation extrapolation method, the new detection statistic increases the decentralized parameter by N times. Therefore, the new detection statistic can significantly improve the detection effect of the small frequency jump and the frequency drift jump. In theory, the greater the cumulative epoch number is, the better the detection effect will be. However, with the increase of N, the innovation covariance matrix is no longer approximately unchanged, and the detection time will increase, which is not conducive to detecting anomalies in time.

The judgment criterion based on innovation extrapolation-accumulation is:

$H_{0}$ : There is no frequency anomaly

T_{D} \leq TH

(26)

$H_{1}$ : There is a frequency anomaly

T_{D} > TH

(27)

The detection threshold TH can be calculated by the given probability of false alarm (PFA).

\begin{matrix} P_{fa} = \int_{TH}^{+ \infty} P (T_{D} | H_{0}) d T_{D} \\ = \int_{TH}^{+ \infty} \frac{1}{2^{\frac{N}{2}} Γ (\frac{N}{2})} {T_{D}}^{\frac{N}{2} - 1} e^{- \frac{T_{D}}{2}} d T_{D} \end{matrix}

(28)

Then, the probability of detection is:

\begin{matrix} P_{D} = \int_{TH}^{+ \infty} P (T_{D} | H_{1}) d T_{D} \\ = \int_{TH}^{+ \infty} \frac{1}{2} {(\frac{T_{D}}{λ^{'}})}^{\frac{N - 2}{4}} e^{- \frac{T_{D} + λ^{'}}{2}} I_{\frac{N}{2} - 1} (\sqrt{λ^{'} T_{D}}) d T_{D} \end{matrix}

(29)

where $I_{a} (•)$ is the first modified Bessel function.

Implementation details

The detailed implementation steps of the proposed algorithm can be described as follows:

Step 1: State estimation vector $\hat{X} (k)$ and error covariance update matrix $\hat{P} (k)$ are obtained by Kalman filter;

Step 2: Select the number of extrapolation. It is assumed to be T. We can get $\hat{X} (k - T + 1)$ and $\hat{P} (k - T + 1)$ by the specific extrapolation process.

Step 3: The detection statistic $λ_{avg}$ is calculated by equations (19)–(21).

Step 4: Selective accumulation number. It is assumed to be N. The new detection statistic is obtained through equation (22).

Step 5: Each value of the detection statistic is compared with the detection threshold TH, which can be confirmed by equation (28) with a given probability of false alarm. If the value of the detection statistic exceeds the threshold, a frequency anomaly is detected.

The proposed algorithm is illustrated schematically in Figure 2.

Figure 2.

Implementation scheme of the innovation extrapolation-accumulation based on the Kalman filter.

Numerical simulation

To evaluate the performance of the innovation extrapolation-accumulation based on the Kalman filter in detecting small frequency jump and frequency drift jump, the following numerical simulation is carried out. According to the clock difference generation method provided by Xu et al.,¹⁹ the numerical simulation produced the clock difference data of 10,000 s duration for a cesium atom clock of a 5071 A, and its frequency data are shown in Figure 3(a). The sampling time is 1 s. The relevant noise parameters are:

{\begin{matrix} q_{1}^{2} = \frac{1}{3} {(3.59 \times 10^{- 11})}^{2} \\ q_{2}^{2} = {(7.92 \times 10^{- 12})}^{2} \\ q_{3}^{2} = \frac{3}{2} {(2.37 \times 10^{- 24})}^{2} \\ R = {(1 \times 10^{- 10})}^{2} \end{matrix}

(30)

Figure 3.

Frequency deviation under different anomalies of numerical simulation: (a) normal, (b) frequency jump, and (c) frequency drift jump.

Frequency jump

Frequency jump is a typical frequency anomaly of atomic clocks. It shows that the mean of frequency deviation suddenly increases or decreases at a certain moment. Therefore, it can be modeled as:

{\tilde{x}}_{2} (k) = {\begin{matrix} x_{2} (k), & k < n_{0} \\ x_{2} (k) + Δ f, & k \geq n_{0} \end{matrix}

(31)

where $n_{0}$ is the time instant at which the frequency anomaly occurs, and $Δ f$ is the size of the jump.

Figure 3(b) shows the frequency jump of size $Δ f = 3 \times 10^{- 11}$ occurring at 5000 s. The thin solid line is the measured frequency value, and the thick solid line represents the predicted frequency value, generated by the Kalman filter. It can be seen that when a frequency jump occurs, the predicted value can closely “track” the measured frequency value.

Optimal extrapolation and accumulation number

According to equations (19) to (22) and (29), the detection performance is related to the extrapolation number and the accumulation number. In Figure 4, the relationship between the accumulation number and the probability of detection (PD) under different extrapolation numbers are described in the case of given frequency jump is $3 \times 10^{- 11}$ and the probability of false alarm is 10⁻⁵. The relationship between the accumulation number and time of detection (TD) is plotted in Figure 5, and the other conditions are the same as Figure 4.

Figure 4.

The PD with a given PFA and frequency jump as a function of the accumulation number under different extrapolation numbers.

Figure 5.

The TD with a given PFA and frequency jump as a function of the accumulation number under different extrapolation numbers.

In Figure 4, it can be seen that under the same extrapolation number, the probability of detection increases with the increase of the accumulation number. Under the same accumulation number, the bigger the extrapolation number is, the higher the probability of detection is.

In Figure 5, we use different marks to correspond to the relationship between the time of detection and the accumulation number under different extrapolation numbers. It can be seen that with the same N, when T is greater than 100, the TD will be longer with the increase of T. When T is less than 100, the TD may be the same in the case of the same N and different T. For example, when T = 70, N = 30; T = 80, N = 30, and T = 90, N = 30, both of which have a detection time of 70 s. In Figure 5, when T is equal to 60, the relationship between different N and TD does not appear. That’s because it was missed detection.

It can be seen that the minimum TD is 70 s, respectively at T = 70, N = 30; T = 70, N = 30, and T = 90, N = 30. Therefore, we plotted the corresponding detection results in these three cases as shown in Figure 6. As can be seen from Figure 6, under the same N, the greater T is, the greater the detection result is, and the easier it is to detect. The detection threshold is 71.59. So we can get the optimal T and N, which are 90 and 30 in the case of a given frequency jump is $3 \times 10^{- 11}$ and the probability of false alarm is 10⁻⁵.

Figure 6.

Detection results of frequency jump under the same accumulation number and different extrapolation numbers.

Detection performance

To evaluate the new method, the detection probability of the innovation detection method, the innovation extrapolation method and the innovation extrapolation-accumulation method at given weak frequency jumps are analyzed. The probability of detection for values of PFA given in equation (32), under different frequency jumps with the sizes given in equation (33) are calculated, which are shown in Figure 7. The extrapolation number is 90 and the accumulation number is 30.

P_{fa} = 10^{- 7}, 5 \times 10^{- 7}, 10^{- 6}, 5 \times 10^{- 6}, 10^{- 5}, 10^{- 4}, 10^{- 3};

(32)

Δ f = 3 \times 10^{- 11}, 5 \times 10^{- 11};

(33)

As can be seen from Figure 7, for the same frequency jump value, the innovation extrapolation-accumulation method has the maximum detection probability, followed by the innovation extrapolation method, and the innovation detection method is the worst. The simulation results also presented as expected: for the same method, the bigger the frequency jump value, the easier the detection and the higher the detection probability. Meanwhile, with the increase of PFA, the PD also increases. Detailed statistical results are shown in Table 1.

Figure 7.

The receiver operating characteristic (ROC) curves obtained with the frequency jumps are given in (33).

Table 1.

The detection results of the three methods when $Δ f = 3 \times 10^{- 11}$ ( $T = 90, N = 30$ ).

	PFA	10⁻⁷	5 × 10⁻⁷	10⁻⁶	5 × 10⁻⁶	10⁻⁵	10⁻⁴	10⁻³
I	PD	2.45 × 10⁻⁷	1.13 × 10⁻⁶	2.19 × 10⁻⁶	1.01 × 10⁻⁵	1.95 × 10⁻⁵	1.72 × 10⁻⁴	1.51 × 10⁻³
	Time/s	/	/	/	/	/	/	/
II	PD	0.002	0.005	0.007	0.018	0.025	0.076	0.202
	Time/s	/	/	/	/	/	/	/
III	PD	0.999999	0.999999	0.999999	0.999999	0.999999	0.999999	0.999999
	Time/s	100	100	100	100	70	70	70

It can be seen from Table 1 that when the extrapolation number is equal to 90, the innovation extrapolation-accumulation method can greatly improve the detection probability only by accumulating the detection statistics of 30 epochs. At the same time, we see that both the innovation method and the innovation extrapolation method have a situation of missing alarm. This is because the abnormal detection statistics are too small to be below the detection threshold.

As far as detection time is concerned, the detection time of the innovation extrapolation-accumulation method is about 1 min. The innovation extrapolation-accumulation method can detect frequency anomalies in short detection time and with high detection probability.

Frequency drift jump

The small frequency drift jump is showed as that the mean value of the frequency drift suddenly starts to change from a certain moment and lasts for a certain time. Therefore, it can be modeled as:

{\tilde{x}}_{2} (k) = {\begin{matrix} x_{2} (k), & k < n_{0} \\ x_{2} (k) + Δ d (k - n_{0}), & k \geq n_{0} \end{matrix}

(34)

where $n_{0}$ is the time instant at which the frequency drift jump occurs, and $Δ d$ is the size of the jump.

Figure 3(c) shows the frequency drift jump with size $Δ d = 1 \times 10^{- 12}$ from 5000 s. The thin solid line is the measured frequency value, and the thick solid line represents the predicted frequency value, generated by the Kalman filter. As can be seen from the figure, when the frequency drift jump is small, there is a certain deviation between the predicted output value of the filter and the measured value of the frequency. However, as time goes on, the Kalman filter gradually converges, and the predicted output value of the filter “tracks” the measured value of the frequency well.

Optimal extrapolation and accumulation number

Similar to frequency jump, the relationship between the accumulation number and the probability of detection under different extrapolation numbers are described in Figure 8 when frequency drift jump $3 \times 10^{- 11} s^{- 1}$ is given and the PFA is 10⁻⁵.

Figure 8.

The PD with a given PFA and frequency drift jump as a function of the accumulation number under different extrapolation numbers.

Figure 8 shows a curvilinear trend similar to Figure 4. The relationship between N and TD is plotted in Figure 9, and the other conditions are the same as Figure 8. It can be seen that when N is less than 5, the detection time corresponding to different extrapolation numbers are very small, but combined with Figure 8, their detection probabilities are very low. As N increases, the detection time increases, but for different extrapolation numbers, the detection time is the same.

Figure 9.

The TD with a given PFA and frequency drift jump as a function of the accumulation number under different extrapolation numbers.

Suppose to ensure that the detection probability exceeds 0.995 and the detection time is as small as possible. Only T = 60, N = 20; T = 70, N = 20; T = 80, N = 20, and T = 100, N = 20 satisfy this assumption. Therefore, we plotted the detection results of these four cases as shown in Figure 10. In Figure 10, when frequency drift jump occurs, the detection results of anomaly will far exceed the threshold value. Even if T = 60, N = 20 corresponding detection results compared to the other three cases are smaller, but it is large enough. This is enough to detect anomalies. At the same time, the smaller T is, the lower the computational complexity is. So we get the optimal extrapolation number and the accumulation number, which are T = 60 and N = 20 in the case of given frequency drift jump is $3 \times 10^{- 11} s^{- 1}$ and the probability of false alarm is 10⁻⁵.

Figure 10.

Detection results of frequency drift jump under the accumulation number and different extrapolation numbers.

Detection performance

The simulation process is similar to the frequency jump. The probability of detection for values of PFA given in equation (35), under different frequency drift jumps with the sizes given in equation (36) are calculated, which are shown in Figure 11.

P_{fa} = 10^{- 7}, 5 \times 10^{- 7}, 10^{- 6}, 5 \times 10^{- 6}, 10^{- 5}, 10^{- 4}, 10^{- 3};

(35)

Δ d = 2 \times 10^{- 11}, 3 \times 10^{- 11} s^{- 1};

(36)

In Figure 11, when frequency drift jump is $3 \times 10^{- 11} s^{- 1}$ and T = 60, N = 20, the detection probability of both the innovation method and the innovation extrapolation method are extremely low, while the new method effectively improves the detection probability by more than 99%. Detailed statistical results are shown in Table 2. In Table 2, we can see that the new method greatly increases the detection probability on the premise of hardly increasing the detection time. The goal of high detection probability and short detection time is realized.

Figure 11.

The ROC curves obtained with the frequency drift jumps are given in (36).

Table 2.

The detection results of the three methods when $Δ d = 3 \times 10^{- 11} s^{- 1}$ ( $T = 60, N = 20$ ).

	PFA	10⁻⁷	5 × 10⁻⁷	10⁻⁶	5 × 10⁻⁶	10⁻⁵	10⁻⁴	10⁻³
I	PD	2.45 × 10⁻⁷	1.13 × 10⁻⁶	2.19 × 10⁻⁶	1.01 × 10⁻⁵	1.95 × 10⁻⁵	1.72 × 10⁻⁴	1.51 × 10⁻³
	Time/s	/	/	/	/	69	33	19
II	PD	5.66 × 10⁻⁴	0.002	0.002	0.006	0.010	0.034	0.111
	Time/s	19	19	18	17	17	16	15
III	PD	0.969653	0.984222	0.988481	0.994926	0.996587	0.999276	0.999905
	Time/s	20	20	20	20	20	20	20

“I” represents the innovation detection method, “II” is the innovation extrapolation detection method, and “III” represents the innovation extrapolation-accumulation detection method. “/” represents a situation of missed detection.

Figure 11 is drawn under the conditions of T = 60 and N = 20. It is the optimal extrapolation number and the accumulation number when $Δ d = 3 \times 10^{- 11} s^{- 1}$ . However, for $Δ d = 2 \times 10^{- 11} s^{- 1}$ , such a condition cannot provide the best detection effect, but compared with the other two methods, the improvement is also quite obvious.

Real data analysis

To obtain the measured data of the frequency jump of atomic clocks, an experimental platform was built. The experimental principle block diagram is shown in Figure 12(a), and part of the measured scenes are shown in Figure 12(b).

Figure 12.

Experimental principle block diagram and test environment: (a) principle block diagram and (b) test-site.

In the experiment, a cesium atomic clock 5071 was used as the atomic clock to be measured, and one output signal of 10 MHz was used as the input of a High-Resolution Phase and Frequency Offset Generator (HROG-10), and the other output signal was used as the input of a Multi-channel Measurement System (MMS). The standard 10 MHz frequency signal of NTSC to MMS as a reference, and 1PPS signal to HROG-10 for initial phase synchronization, and finally data acquisition through a computer. The HROG-10 is used to adjust the frequency offset.

Data were collected from September 1, 2020 to September 2, 2020, and weak frequency jumps of ±3 × 10⁻¹¹ were inserted through the HROG-10 on September 2. They correspond to 11015, 25895, 40535, and 54840 s. The duration is about 4 h. The sampling interval is 1 s. The whole test process was conducted in a constant temperature and humidity environment. During the test, the laboratory was closed and equipped with UPS emergency power supply to avoid interference of other factors on the test results.

Figure 13(a) shows the deviation of the normal atomic clock frequency on September 1, and Figure 13(b) shows the deviation on September 2. In the figure, “R” represents the measured frequency deviation value, and “P” represents the Kalman forecast value.

Figure 13.

Real data of frequency deviation: (a) normal and (b) frequency jump.

We use normal data to calculate the noise parameters of the atomic clock. The relevant noise parameters are:

{\begin{matrix} q_{1}^{2} = \frac{1}{3} {(1.55 \times 10^{- 10})}^{2} \\ q_{2}^{2} = {(5.95 \times 10^{- 12})}^{2} \\ q_{3}^{2} = \frac{3}{2} {(1.43 \times 10^{- 24})}^{2} \\ R = {(1 \times 10^{- 10})}^{2} \end{matrix}

(37)

So the covariance matrix of the Kalman filter can be initialized well, otherwise, it will cause a large error, and even the algorithm does not converge.

To highlight the advantages of the new method, all the three detection methods mentioned in this paper have been used to test the measured data. According to 3.1.1, for $Δ f = 3 \times 10^{- 13}$ , the optimal extrapolation number is 90 and the optimal accumulation number is 30. Therefore, based on this condition and PFA is equal to 10⁻⁵, the same analysis is carried out on the measured data. The detection results are shown in Figure 14.

Figure 14.

Detection results of real data: (a) innovation, (b) extrapolation, and (c) extrapolation-accumulation.

Figure 14(a) only draws the detection statistics of the innovation detection method. It can be seen from the figure that the detection statistics are very small and the detection threshold is 19.51. It is impossible to distinguish the abnormal detection statistics from the normal detection statistics and cannot identify the anomalies. Because of this, the detection results of the three methods are not drawn together. Otherwise, the detection results of the innovation method are too small to be covered up.

Figure 14(b) shows the detection statistics value of the innovation extrapolation method and its threshold value of 19.51, as well as the frequency deviation value of 10¹¹ times. The detection statistics corresponding to four anomalies are 7.92, 6.24, 6.78, and 7.20. Although there is a difference between the abnormal and the normal, the detection results of the abnormal do not reach the threshold value, so the anomaly identification is still unable to be achieved. This phenomenon is consistent with the results of our simulation.

Figure 14(c) is plotted in a similar way to Figure 14(b), showing the detection results of the innovation extrapolation-accumulation method. The value of the detection statistic increases at the time value corresponding to the frequency jumps. The detection statistics corresponding to four anomalies are 120.08, 120.67, 106.93, and 128.04. The detection threshold of the new method is 71.59, and the detection statistics at the abnormal place are far beyond the threshold, which can be detected well. The detection moments of the four anomalies are 11100, 25980, 40620, and 54840 s, and the detection time is 85 s. From the point of view of detection time, compared with the detection time obtained by simulation, it lags 15 s. This could be due to HROG-10. Because when we put in a frequency offset, the HROG-10 is relatively slow to adjust the frequency offset, rather than suddenly changing the frequency offset. This result is acceptable and understandable to us.

In Figure 14, you can see that the innovation extrapolation-accumulation based on the Kalman filter of the weak frequency anomaly shows sufficient superiority. If the frequency jumps were larger, the second method might also be able to pick up the anomaly, but the advantage of the new method is that it can be more sensitive to weak frequency anomalies as the detection values are larger. The innovation extrapolation-accumulation method can expand the detection statistics of frequency anomalies, while the detection statistics of no anomalies remain unchanged, widening the gap between the detection statistic results and making it easier to detect anomalies, which is consistent with the previous theoretical proof.

Conclusions

In this paper, a new method that innovation extrapolation-accumulation based on the Kalman filter is proposed for the detection of weak frequency anomalies in atomic clocks. The innovation extrapolation-accumulation detection method shows the advantages of high detection probability and short detection time whether it is small frequency jump or frequency drift jump. Both the simulation results and the measured data show that the innovation extrapolation-accumulation detection method is an effective method to detect the small frequency anomalies of atomic clocks, which can realize the real-time and high-performance detection of the small frequency anomalies of atomic clocks.

Our team is committed to the development of a seamless switcher for the main and backup atomic clocks, and it is expected that the new detection algorithm will be applied to it in the future.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors thank the ZFS Satellite Navigation Test System of China (Grant No. ZFS20001D) funding support.

ORCID iD

Ran Yan

References

Cernigliaro

Sesia

. Satellite clocks characterization and monitoring for global navigation satellite systems. In: 2011 XXXth URSI general assembly and scientific symposium, Istanbul, Turkey, 13–20 August 2011, pp.1–4. IEEE.

Wang

Rochat

. Backup hydrogen maser steering algorithm for Galileo precise timing facility. In: IEEE international frequency control symposium joint with the 21st European frequency and time forum, Geneva, Switzerland, 29 May–1 June 2007, pp.463–466. IEEE.

Lee

Kim

Lee

YJ.

Protecting signal integrity against atomic clock anomalies on board GNSS satellites. IEEE Trans Instrum Meas 2011; 60(7): 2738–2745.

Weiss

Beard

Shome

GPS signal integrity dependencies on atomic clocks. In: Proceedings of the 38th annual precise time and timing interval systems and applications meeting, Reston, VA, 7–9 December 2006.

Kaplan

Hegarty

CJ.

Understanding GPS: principles and applications. 2nd ed. Boston, MA: Artech House Inc, 2006.

Riley

WJ.

Algorithms for frequency jump detection. Metrologia 2008; 45(6): 154–161.

Nunzi

Carbone

Monitoring signal integrity of atomic clocks by means of the GLRT. Metrologia 2008; 45(6): S103.

Galleani

The dynamic Allan variance II: a fast computational algorithm. IEEE Trans Ultrason Ferroelectr Freq Control 2009; 57(1): 182–188.

Galleani

Tavella

The dynamic Allan variance. IEEE Trans Ultrason Ferroelectr Freq Control 2009; 56(3): 450–464.

10.

Galleani

Tavella

Detection of atomic clock frequency jumps with Kalman filter. IEEE Trans Ultrason Ferroelectr Freq Control 2012; 59(3): 504–509.

11.

Song

Dong

WJ.

Detecting an atomic clock frequency anomaly using an adaptive Kalman filter algorithm. Metrologia 2018; 55(3): 150.

12.

Huang

Gong

Zhu

, et al. Detection of weak frequency anomalies for atomic clocks with a Kalman filter. Metrologia 2014; 51(3): 183–190.

13.

Huang

Gong

Detection of weak frequency jumps for GNSS onboard clocks. IEEE Trans Ultrason Ferroelectr Freq Control 2014; 61(5): 747–755.

14.

Huang

Zhao

, et al. An integrity monitoring algorithm for GNSS satellite atomic clocks. J Navig 2019; 72(6): 1533–1549.

15.

Fan

Shi

XH.

An algorithm for the rubidium atomic clock control based on the Kalman filter. J Astronaut 2015; 36(1): 90–95.

16.

Diez

D’Angelo

Fernández

Clock errors simulation and characterisation. In: Proceedings of the 19th international technical meeting of the satellite division of the institute of navigation, Fort Worth, TX, 26–29 September 2006, pp.815–821.

17.

Stein

Filler

. Kalman filter analysis for real time applications of clocks and oscillators. In: Proceedings of the 42nd annual frequency control symposium, Baltimore, MD, 1–3 June 1988, pp.447–452. IEEE.

18.

Song

Zhang

, et al. RAIM method based on robust extended Kalman filter and extrapolation-accumulation. Syst Eng Electron 2017; 39(9): 2094–2099.

19.

Liu

, et al. A modified method for clock error simulation of high precision frequency source. J Astronaut 2017; 38(9): 998–1004.

The weak frequency anomaly detection method of atomic clocks based on Kalman filter and extrapolation-accumulation

Abstract

Keywords

Introduction

Innovation extrapolation-accumulation

The Kalman filter model

Innovation extrapolation method

Chi-square detection of multi-epoch accumulation

Implementation details

Numerical simulation

Frequency jump

Optimal extrapolation and accumulation number

Detection performance

Frequency drift jump

Optimal extrapolation and accumulation number

Detection performance

Real data analysis

Conclusions

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References