P-step Kalman-like unbiased FIR algorithm for detecting anomalies of atomic clocks

Abstract

An extended measurement model anomaly detection algorithm for atomic clocks is proposed based on the p-step Kalman-like iterative Unbiased finite impulse response (UFIR) algorithm. The novel of the proposed algorithm is that it breaks the disadvantage of initializing the noise matrix in the traditional Kalman filter improved algorithms. By accumulating the prediction residuals generated by the process of p-step Kalman-like iterative, the detection statistics are constructed to realize the weaker phase and frequency jump detection of atomic clocks. Simulation results show the effectiveness of the proposed algorithm, and Monte-Carlo (M-C) experiment is used to verify the theoretical calculation of detection probability. We also analyze the measured data of a cesium atomic clock in our timekeeping laboratory, and the results show that the new method is more sensitive to weak frequency jump than the improved algorithm based on the Kalman filter.

Introduction

Atomic clocks have been widely used in many fields such as wireless communication, satellite navigation, aerospace, astronomical observation, and time-frequency measurement, providing stable and reliable high-precision time-frequency reference, and becoming one of the important devices to ensure the stable operation of the system. Perhaps the most famous application of atomic clocks is satellite navigation. In Global Navigation Satellite Systems (GNSS), such as GPS, GLONASS, BeiDou, or Galileo, atomic clocks are used in both onboard and ground stations, which make it possible for users to obtain accurate location and time information. ¹ However, it is subject to various factors, especially in the space station. Abnormalities in atomic clocks will seriously damage the stability and reliability of the system’s output signals, and ultimately the users will suffer. Although the time error up to 10 ns will lead to a positioning error of more than 3 m, so it is very important and worthy of research to quickly and accurately identify atomic clock anomalies.

Anomaly detection of atomic clocks is a continuous research direction for scientists. Riley ² proposed three simple frequency jump detection techniques, they are BLKAVG, SEQAVG, and CUSUM. Subsequently, the Generalized Likelihood Ratio Test (GLRT) was gradually widely used in the detection of atomic clock frequency anomalies.^3–7 However, as a statistical method, it could not accurately obtain the relationship between the degree of anomalies and the parameters to be set. Dynamic Allan Variance (DAVAR) and spectrogram are also used to monitor the stability of atomic clocks. ⁸ INRIM has developed clock analysis tools for monitoring atomic clocks in space, including the DAVAR and Dynamic Hadamard Deviation (DHDEV). ⁹ Since Galleani and Tavella proposed using the Kalman filter to detect the abnormal frequency of atomic clocks in 2012, improved algorithms based on the Kalman filter have been continuously proposed to detect abnormal frequency of atomic clocks,^10–14 including one-step prediction, N-step prediction, adaptive Kalman filter, etc. Huang et al. ¹⁵ proposed an algorithm for GNSS satellite atomic clock integrity monitoring based on the parity transformation.

However, some improved algorithms based on Kalman can detect atomic clock anomalies better, but they rely on the initialization of the noise matrix and error matrix. Once the state noise matrix cannot be initialized accurately, the Kalman algorithm will introduce large errors and even lead to non-convergence of the algorithm, which is unacceptable to engineers.^16,17

In this paper, a new detection algorithm is proposed. Its innovation is that it can ignore the initialization of the noise matrix and detect weak atomic clock anomalies. This is a point that all previous algorithms based on the Kalman filter cannot achieve. It does not depend on the initialization of the noise matrix on the premise of high detection probability and real-time detection, which is beneficial to the realization of engineering. The p-step UFIR prediction algorithm has been applied in GPS time synchronization and timekeeping technology as early as 2006, achieving higher accuracy than Kalman prediction.^18–20 This is the first time that it has been extended to anomaly detection in atomic clocks. The p-step Kalman-like UFIR detection algorithm can effectively postpone the influence of observed outliers, increase the difference between the predicted value and observed value at an abnormal time, and thus increase the anomaly detection statistics, which is more conducive to detection. Please refer to Appendix B for a detailed description. The detection statistics are constructed by using the residual generated in the process of the Kalman-like iteration. The detection statistics in several epochs are accumulated to increases the non-centrality parameter of the chi-square distribution. The results show that the method can effectively detect the phase or frequency jumps of the atomic clock by analyzing the abnormal data of the simulated cesium atomic clock and measured the cesium atomic clock in our laboratory.

The paper is organized as follows. In Section 2, the extended model of the p-step FIR filter is described in detail. In Section 3, an integrity detecting algorithm based on p-step Kalman-like UFIR is proposed and the option of the horizon is also discussed. In Section 4, the theoretical simulation results under different phase and frequency jumps of a cesium atomic clock are given and verified by the M-C test. In Section 5, we apply the new method and the improved detection algorithm based on the Kalman filter to the analysis of the measured cesium atomic clock anomaly data. The conclusion is in Section 6.

The extended model of the p-step FIR filter

The clock time-invariant state and observation equations can be respectively expressed as:

$${\mathbf{x}}_n=\mathbf{A}·{\mathbf{x}}_{n-1}+{\mathbf{w}}_n$$

(1)

$$z_n=\mathbf{H}·{\mathbf{x}}_n+v_n$$

(2)

where ${\mathbf{x}}_n=\left[\begin{array}{c}{x}_{1n}\hfill \\ \multicolumn{1}{c}{x_{2n}}\\ \multicolumn{1}{c}{x_{3n}}\end{array}\right]$ , $\mathbf{A}=\left[\begin{array}{ccc}1\hfill & T_s\hfill & 0.5T_s^2\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & T_s\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$ , ${\mathbf{w}}_n=\left[\begin{array}{c}{w}_{1n}\hfill \\ \multicolumn{1}{c}{w_{2n}}\\ \multicolumn{1}{c}{w_{3n}}\end{array}\right]$ , $\mathbf{H}=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]$ . $x_{1n}$ is the phase deviation and $x_{2n}$ is the frequency deviation. $x_{3n}$ is the drift rate. $T_s$ represents the sampling interval. The process noise vector ${\mathbf{w}}_n$ has three components: $w_{1n}$ represents the random variation due to White Phase Noise (WPN), and $w_{2n}$ accounts for White Frequency Noise (WFN), $w_{3n}$ is the random variation due to Random Walk Frequency Noise (RWFN), respectively. ²¹ ${\mathbf{w}}_n$ is an approximately multivariate Gaussian random vector with zero mean and a covariance of Q. $z_n$ is the measured data and $v_n$ is White Gaussian Noise (WGN) of variance $R$ .

The p-step FIR prediction algorithm estimates the clock state at n + p (p > 0) with N measurements on a horizon from m = n − N + 1 to n, as shown in Figure 1.

Figure 1.

Clock prediction problem solved with p-step FIR structure.

Therefore, (1) and (2) need to be extended on [m, n]. That can be done if to use the recursively computed solutions and write

$${\mathbf{x}}_m={\mathbf{x}}_{m-1}+{\mathbf{w}}_m$$

(3)

$${\mathbf{x}}_{m+1}=\mathbf{A}{\mathbf{x}}_m+{\mathbf{w}}_{m+1}=\mathbf{A}{\mathbf{x}}_m+\mathbf{A}{\mathbf{w}}_m+{\mathbf{w}}_{m+1}$$

(4)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{x}}_{m+2}=\mathbf{A}{\mathbf{x}}_{m+1}+{\mathbf{w}}_{m+2}={\mathbf{A}}^2{\mathbf{x}}_m+{\mathbf{A}}^2{\mathbf{w}}_m}\\ +\mathbf{A}{\mathbf{w}}_{m+1}+{\mathbf{w}}_{m+2}\hfill \end{array}$$

(5)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{x}}_n=\mathbf{A}{\mathbf{x}}_{n-1}+{\mathbf{w}}_n={\mathbf{A}}^{n-m}{\mathbf{x}}_m+{\mathbf{A}}^{n-m}{\mathbf{w}}_m}\\ +{\mathbf{A}}^{n-m-1}{\mathbf{w}}_{m+1}+\cdots +{\mathbf{w}}_n\hfill \end{array}$$

(6)

At the initial horizon point, equation (3) implies that the initial state must be known in advance or estimated optimally. ²² Equation (7) can be obtained through the recursively computed solutions from equations (3) to (6).

$${\mathbf{X}}_{n,m}={\mathbf{A}}_{n-m}·{\mathbf{x}}_m+{\mathbf{B}}_{n-m}·{\mathbf{W}}_{n,m}$$

(7)

where the extended vectors are ${\mathbf{X}}_{n,m}={[{\mathbf{x}}_n{\mathbf{x}}_{n-1}\dots {\mathbf{x}}_m]}^T$ , ${\mathbf{A}}_i={\left[\begin{array}{cccc}\left({\mathbf{A}}^i\right)\hfill & ···\hfill & \mathbf{A}\hfill & \mathbf{I}\hfill \end{array}\right]}^T$ , ${\mathbf{W}}_{n,m}={\left[\begin{array}{cccc}{\mathbf{w}}_n\hfill & {\mathbf{w}}_{n-1}\hfill & ···\hfill & {\mathbf{w}}_m\hfill \end{array}\right]}^T$ .

Similarly, extend the observation equation (2) and write

$$z_m=\mathbf{H}·{\mathbf{x}}_m+v_m=\mathbf{H}·{\mathbf{x}}_m+\mathbf{H}·{\mathbf{w}}_m+v_m$$

(8)

$$\begin{array}{c}\multicolumn{1}{c}{z_{m+1}=\mathbf{H}·{\mathbf{x}}_{m+1}+v_{m+1}=\mathbf{HA}{\mathbf{x}}_m+\mathbf{HA}{\mathbf{w}}_m}\\ +\mathbf{H}{\mathbf{w}}_{m+1}+v_{m+1}\hfill \end{array}$$

(9)

$$\begin{array}{c}\multicolumn{1}{c}{\begin{array}{c}\multicolumn{1}{c}{z_n=\mathbf{H}{\mathbf{x}}_n+v_n=\mathbf{H}{\mathbf{A}}^{n-m}{\mathbf{x}}_m+\mathbf{H}{\mathbf{A}}^{n-m}{\mathbf{w}}_m}\\ +\mathbf{H}{\mathbf{A}}^{n-m-1}{\mathbf{w}}_{m+1}+\cdots +\mathbf{H}{\mathbf{w}}_n+v_n\hfill \end{array}}\end{array}$$

(10)

Therefore, the extended observation vector ${\mathbf{Z}}_{n,m}$ is shown in equation (11).

$${\mathbf{Z}}_{n,m}={\mathbf{H}}_{n-m}·{\mathbf{x}}_m+{\mathbf{C}}_{n-m}·{\mathbf{W}}_{n,m}+{\mathbf{V}}_{n,m}$$

(11)

where the extended vectors are ${\mathbf{Z}}_{n,m}={[z_n{z}_{n-1}\dots z_m]}^T$ , ${\mathbf{H}}_i={\left[\begin{array}{cccc}\mathbf{H}{\mathbf{A}}^i\hfill & ···\hfill & \mathbf{HA}\hfill & \mathbf{H}\hfill \end{array}\right]}^T$ , ${\mathbf{V}}_{n,m}={\left[\begin{array}{cccc}{v}_n\hfill & v_{n-1}\hfill & ···\hfill & v_m\hfill \end{array}\right]}^T$ . ${\mathbf{B}}_{n-m}$ and ${\mathbf{C}}_{n-m}$ are both coefficient matrices produced by extension.

$${\mathbf{B}}_{n-m}=\left[\begin{array}{ccccc}\mathbf{I}\hfill & \mathbf{A}\hfill & \cdots \hfill & {\mathbf{A}}^{n-m-1}\hfill & {\mathbf{A}}^{n-m}\hfill \\ \multicolumn{1}{c}{0}& \mathbf{I}\hfill & \cdots \hfill & {\mathbf{A}}^{n-m-2}\hfill & {\mathbf{A}}^{n-m-1}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 0\hfill & \mathbf{I}\hfill \end{array}\right]$$

(12)

$${\mathbf{C}}_{n-m}=\left[\begin{array}{ccccc}\mathbf{H}\hfill & \mathbf{HA}\hfill & \cdots \hfill & \mathbf{H}{\mathbf{A}}^{n-m-1}\hfill & \mathbf{H}{\mathbf{A}}^{n-m}\hfill \\ \multicolumn{1}{c}{0}& \mathbf{H}\hfill & \cdots \hfill & \mathbf{H}{\mathbf{A}}^{n-m-2}\hfill & \mathbf{H}{\mathbf{A}}^{n-m-1}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & 0\hfill & \mathbf{H}\hfill \end{array}\right]$$

(13)

where I stands for identity matrix.

Assuming that the p-step prediction is ${\tilde{\mathbf{x}}}_{n+p|n}$ , and N measurements on a horizon are used as the input of the filter, then ${\tilde{\mathbf{x}}}_{n+p|n}$ can be expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{{\tilde{\mathbf{x}}}_{n+p|n}={\mathbf{G}}_p{\mathbf{Z}}_{n,m}={\mathbf{G}}_p{\mathbf{H}}_{n-m}{\mathbf{x}}_m+{\mathbf{G}}_p{\mathbf{C}}_{n-m}{\mathbf{W}}_{n,m}}\\ +{\mathbf{G}}_p{\mathbf{V}}_{n,m}\hfill \end{array}$$

(14)

where ${\mathbf{G}}_p$ is the gain of the FIR.

UFIR filter means to ensure the unbiased estimation, as shown in equation (15), but cannot guarantee the minimum mean square error (MSE).^23,24

$$\mathbb{E}\left\{{\tilde{\mathbf{x}}}_{n+p|n}\right\}=\mathbb{E}\left\{{\overline{\mathbf{x}}}_{n+p|n}\right\}=\mathbb{E}\left\{{\mathbf{x}}_{n+p|n}\right\}$$

(15)

It has been shown in Ibarra-Manzano ¹⁶ that the unbiasedness constraint for linear FIR filters is equal to the deadbeat constraint associated with noiseless both system and measurement.

The estimated value of p-step UFIR can be expressed as follow:

$${\overline{\mathbf{x}}}_{n+p|n}={\mathbf{A}}^{n-m+p}{\left({\mathbf{H}}_{n-m}^T{\mathbf{H}}_{n-m}\right)}^{-1}{\mathbf{H}}_{n-m}^T{\mathbf{Z}}_{n,m}$$

(16)

The specific derivation process can be referred to the references.^25,26

The integrity detecting algorithm

Although Section 2 presents a batch form of p-step UFIR estimators based on convolution rules, when N is large, all matrix and vector dimensions are large, which will cause computational problems. For fast computation, the unbiased estimator for batch processing can be expressed in the form of a Kalman-like recursion. ²⁷

P-step Kalman-like UFIR algorithm

The principle of p-step Kalman-like UFIR iterative is shown in Figure 2. One can observe that the iterative computation starts at m + v with v = K and finishes at v = N − 1. The true estimate is taken at each v = N − 1 and the procedure is repeated recursively.

Figure 2.

Operational principles of the p-step Kalman-like UFIR.

The predicted vector of the atomic clock state ${\overline{\mathbf{x}}}_{s+p}$ can be calculated by using the initial vector ${\mathbf{x}}_m$ using equation (16), which can be used as the initial vector of the whole iteration process. As time goes on, new observations $z_{m+v}$ are constantly injected and new state predictions ${\overline{\mathbf{x}}}_{m+v+p|m+v}$ are updated based on the previous moment. It is worth noting here that the observed lags by p steps each time, thus effectively postponing the “pollution” of outliers, which is beneficial to the construction of anomaly detection statistics. Please refer to Appendix B for a detailed description.

The detailed algorithm process is described as follows^16,27 :

Number of atomic clock states K = 3, The horizon N as input. The algorithm starts at n = N.

Set the intermediate variables

$$m=n-N+1$$

(17)

$$s=m+K-1$$

(18)

Set auxiliary matrices

$$\mathbf{\Xi }={\mathbf{A}}^T{\mathbf{H}}^T\mathbf{HA}$$

(19)

$${\mathbf{F}}_{K-1}={\mathbf{A}}^{K-1}{\left({\mathbf{H}}_{K-1}^T{\mathbf{H}}_{K-1}\right)}^{-1}{\left({\mathbf{A}}^{K-1}\right)}^T$$

(20)

According to the observations in the range of [m,s], the predicted vector ${\overline{\mathbf{x}}}_{s+p|s}$ is calculated as the initial vector of iteration.

$${\overline{\mathbf{x}}}_{s+p|s}={\mathbf{A}}^{K-1+p}{\left({\mathbf{H}}_{K-1}^T{\mathbf{H}}_{K-1}\right)}^{-1}{\mathbf{H}}_{K-1}^T{\mathbf{Z}}_{s,m}$$

(21)

The iterative starts with v = K and ends with v = N − 1.

$${\mathbf{F}}_v=\mathbf{A}{\mathbf{F}}_{v-1}{\mathbf{A}}^T-\mathbf{A}{\mathbf{F}}_{v-1}\times {(\mathbf{I}+\mathbf{\Xi }{\mathbf{F}}_{v-1})}^{-1}\mathbf{\Xi }{\mathbf{F}}_{v-1}{\mathbf{A}}^T$$

(22)

$$\begin{array}{c}{\overline{\mathbf{x}}}_{m+v+p|m+v}=\mathbf{A}{\overline{\mathbf{x}}}_{m+v+p-1|m+v-1}\hfill \\ +{\mathbf{A}}^p{\mathbf{F}}_v{\mathbf{H}}^T(z_{m+v}-\mathbf{H}{\mathbf{A}}^{1-p}{\overline{\mathbf{x}}}_{m+v+p-1|m+v-1})\hfill \end{array}$$

(23)

Note that (23) is the Kalman-like estimate with the bias correction gain ${\mathbf{G}}_{m+v,p}={\mathbf{A}}^p{\mathbf{F}}_v{\mathbf{H}}^T$ that, unlike the Kalman gain, does not depend on noise and initial errors.

The p-step UFIR estimation error can be evaluated using the following expression

$${\mathbf{P}}_{n+p}=\mathbb{E}\left[({\mathbf{x}}_{n+p}-{\overline{\mathbf{x}}}_{n+p|n}){({\mathbf{x}}_{n+p}-{\overline{\mathbf{x}}}_{n+p|n})}^T\right]$$

(24)

According to equations (6) and (24), the error matrix of the p-step Kalman-like UFIR algorithm can be deduced as follows:

$$\begin{array}{c}{\mathbf{P}}_{n+p}^{\mathrm{UB}}=\hfill \\ \{\begin{array}{cc}(\mathbf{I}-{\mathbf{G}}_n\mathbf{H})(\mathbf{A}{\mathbf{P}}_{n-1}{A}^T+\mathbf{Q}){(\mathbf{I}-{\mathbf{G}}_n\mathbf{H})}^T+{\mathbf{G}}_{n}R{\mathbf{G}}_n^T\hfill & p=0\hfill \\ \multicolumn{1}{c}{\begin{array}{c}(\mathbf{A}-{\mathbf{G}}_{n+p}\mathbf{H}{\mathbf{A}}^{1-p}){\mathbf{P}}_{n+p-1}(\mathbf{A}-{\mathbf{G}}_{n+p}\mathbf{H}{\mathbf{A}}^{1-p}{)}^T\hfill \\ \multicolumn{1}{c}{+\mathbf{Q}+{\mathbf{G}}_{n+p}R{\mathbf{G}}_{n+p}^T}\end{array}}& p>0\hfill \end{array}\hfill \end{array}$$

(25)

${\mathbf{P}}_{n+p}^{\mathrm{UB}}$ is the upper bound (UB) on MSE because the process noise covariances are accumulated at each iteration. The MSE lower bound (LB) can be found without the system noise. ²³ The ${\mathbf{P}}_{n+p}^{\mathrm{LB}}$ can be found by iterating:

$${\mathbf{P}}_{n+p}^{LB}=(\mathbf{A}-{\mathbf{G}}_{n+p}\mathbf{H}{\mathbf{A}}^{1-p}){\mathbf{P}}_{n+p-1}{(\mathbf{A}-{\mathbf{G}}_{n+p}\mathbf{H}{\mathbf{A}}^{1-p})}^T+{\mathbf{G}}_{n+p}R{\mathbf{G}}_{n+p}^T$$

(26)

So far, we also get the Kalman-like iteration form of the error matrix, including upper bound and lower bound.

The problem with (25) arises instantly when the noise behavior is unknown as it usually is. We explore an effective method for detecting anomalies of atomic clocks when the system noise matrix is unknown. We use the lower bound of the error matrix as the error matrix in the whole detection algorithm to further realize the anomaly detection of the atomic clock.

Innovative construction of p-step UFIR abnormal detection statistics

The predicted residual is expressed as:

$$r_{n+p}=z_{n+p}-\mathbf{H}{\overline{\mathbf{x}}}_{n+p|n}$$

(27)

When there is no anomaly, the residual nearly follows a Gaussian distribution with a mean of zero, which satisfies:

$$r_{n+p}~\mathrm{N}(0,D_{n+p})$$

(28)

where $D_{n+p}={\mathbf{HP}}_{n+p}^{\mathrm{LB}}\mathbf{H}\prime +R$ .

When the anomaly occurs, the mean of the residual will not be equal to zero, which satisfies:

$$r_{n+p}~\mathrm{N}(\mu ,D_{n+p})$$

(29)

where $\mu$ represents the size of the abnormal.

For the iterative process, we can construct:

$$r_{v-K+1}=z_{m+v+p}-\mathbf{H}{\overline{\mathbf{x}}}_{m+v+p|m+v}$$

(30)

$$D_{v-K+1}={\mathbf{HP}}_{m+v+p}^{\mathrm{LB}}\mathbf{H}\prime +R$$

(31)

where $K\le v\le N-1$ .

According to the residuals generated by the iterative process, detection statistics are adopted:

$${\lambda }_{\mathrm{avg}}=r_{\mathrm{avg}}\prime D_{\mathrm{avg}}^{-1}{r}_{\mathrm{avg}}$$

(32)

$$r_{\mathrm{avg}}=(D_{\mathrm{avg}}^{-1}{)}^{-1}\sum _{v=K}^{N-1}{D}_{v-K+1}^{-1}{r}_{v-K+1}$$

(33)

$$D_{\mathrm{avg}}={(\sum _{v=K}^{N-1}(D_{v-K+1}^{-1})}^{-1}$$

(34)

To better detect weak phase or frequency jump, the new detection statistics can be accumulated by M epochs from the value obtained in equation (32):

$$T_D=\sum _{i=1}^M{\lambda }_{\mathrm{avg}}(i)$$

(35)

where ${\lambda }_{\mathrm{avg}}(i)$ is the detection statistic constructed by the iterative process. According to the statistical distribution theory, when the observed value is not abnormal, ${\lambda }_{\mathrm{avg}}$ obeys the chi-square distribution with 1 degree of freedom. When a jump occurs, it follows a noncentral chi-square distribution with 1 degree of freedom. The non-central parameter is $\lambda$ .

$$\lambda =\mu \prime D_{\mathrm{avg}}^{-1}\mu$$

(36)

When the observed value is not abnormal, the new statistic $T_D$ obeys the chi-square distribution with M degree of freedom. When a jump occurs, it follows a noncentral chi-square distribution with M degree of freedom. The non-central parameter is $\lambda \prime$ .

$${\lambda }^{\prime }=\sum _{i=1}^M\lambda =\sum _{i=1}^M\mu \prime D_{\mathrm{avg}}^{-1}\mu \approx M\mu \prime D_{\mathrm{avg}}^{-1}\mu$$

(37)

where M stands for accumulated epochs, $\mu$ is the size of the abnormal, and $D_{\mathrm{avg}}$ is the covariance matrix calculated in the iterative process. By observing equations (36) and (37), it can be found that equation (37) is the accumulation of M detection statistics of equation (36). If M is not very large, we can assume that $D_{\mathrm{avg}}$ remains constant. The ratio of the two non-central parameters is:

$$\frac{\lambda \prime }{\lambda }\approx \frac{M\mu \prime D_{\mathrm{avg}}^{-1}\mu }{\mu \prime D_{\mathrm{avg}}^{-1}\mu }=M$$

(38)

It can be seen from equation (38) that the new detection statistics accumulated by M epochs are about M times larger than that in equation (36). Therefore, compared with the previous detection statistics, the new detection statistics can increase the statistical value of abnormal places and keep the statistic value of no abnormal places, which is beneficial for anomaly detection.

The judgment criterion based on the p-step Kalman-like UFIR algorithm is:

$H_0$ : There is no anomaly

$$T_D\le \mathrm{TH}$$

(39)

$H_1$ : There is an anomaly

$$T_D>\mathrm{TH}$$

(40)

The detection threshold TH can be calculated by the given probability of false alarm. ²⁸

$${\mathrm{P}}_{fa}={\int }_{TH}^{+\infty }\mathrm{P}(T_D|H_0)d{T}_D={\int }_{TH}^{+\infty }\frac{1}{2^{\frac{M}{2}}\mathrm{\Gamma }(\frac{M}{2})}{T_D}^{\frac{M}{2}-1}{e}^{-\frac{T_D}{2}}d{T}_D$$

(41)

Then, the probability of detection is:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{P}}_D={\int }_{\mathrm{TH}}^{+\infty }\mathrm{P}(T_D|H_1)d{T}_D}\\ ={\int }_{\mathrm{TH}}^{+\infty }\frac{1}{2}{\left(\frac{T_D}{\lambda \prime }\right)}^{\frac{M-2}{4}}{e}^{-\frac{T_D+\lambda \prime }{2}}{\mathrm{I}}_{\frac{M}{2}-1}\hfill \\ \left(\sqrt{\lambda \prime T_D}\right)d{T}_D\hfill \end{array}$$

(42)

where ${\mathrm{I}}_a(•)$ is the first modified Bessel function.

The horizon

The error matrix of the p-step Kalman-like UFIR algorithm is related to the length of the horizon N and the number of prediction steps p. Shmaliy and Ibarra-Manzano ²⁹ mentioned that if N is fixed, the prediction error increases with the increase of p. If p is fixed, the prediction error decreases as N increases.

Given ${\mathbf{P}}_{n+p}^{\mathrm{UB}}$ , by (25), the optimal horizon N_opt can be estimated by minimizing the trace of ${\mathbf{P}}_{n+p}^{\mathrm{UB}}$ as^23,27:

$${\stackrel{⌢}{N}}_{opt}=\arg {\min }_N\left\{tr\left({\mathbf{P}}_n^{UB}(N)\right)\right\}$$

(43)

Note that equation (43) is applicable if $R$ and noise covariance $\mathbf{Q}$ are both known. In the case that the state noise matrix Q is unknown, the horizon N can be obtained through the bias correction gain ${\mathbf{G}}_{m+v,p}={\mathbf{A}}^p{\mathbf{F}}_v{\mathbf{H}}^T$ in the iterative process. ³⁰ In other words, when equation (44) is satisfied, the horizon N can be obtained. $\Vert {\mathbf{G}}_{i,p}{\Vert }_2$ represents the two-norms of the bias correction gain matrix.

The horizon N obtained by this method may not guarantee the minimum estimation error, but for the UFIR algorithm, it is low sensitive to the value of horizon. ³¹

$$\{\begin{array}{c}\mathrm{\Delta }{\Vert {\mathbf{G}}_{m+v,p}\Vert }_2=|{\Vert {\mathbf{G}}_{m+v+1,p}\Vert }_2-{\Vert {\mathbf{G}}_{m+v,p}\Vert }_2|\to 0\hfill \\ \multicolumn{1}{c}{{\Vert {\mathbf{G}}_{m+v,p}\Vert }_2\to {\Vert {\mathbf{G}}_{\infty ,p}\Vert }_2}\end{array}$$

(44)

The horizon N is solved for the phase and frequency anomalies, respectively. By comparing the situation that the increment of the bias correction gain changes with the increase of the number of iterations N under different predicted steps, as shown in Figure 3.

Figure 3.

The incremental gain Δ||G_i,p||₂ given by (44) as functions of the number of iterations for the phase and frequency detection: (a) phase jump and (b) frequency jump.

As the phase fluctuates greatly in a short time, the predicted step number p should not be too large, as shown in Figure 3(a). Therefore, for phase jump detection, the horizon N = 50.

However, for frequency detection, if the number of predicted steps is too small, the anomaly cannot be distinguished. Therefore, the number of predicted steps should be larger, as shown in Figure 3(b). Therefore, for frequency detection, the horizon N = 100.

In Appendix A, we will explain that there is little difference between the horizon N obtained by satisfying equation (44) and the optimal horizon N_opt for UFIR Kalman-like atomic clock anomaly detection algorithms.

Numerical simulation

To evaluate the performance of the p-step Kalman-like UFIR algorithm in detecting small phase jump and frequency jump, the following numerical simulations are performed. According to the simulation generation method of atomic clock difference provided in the references,^32,33 we simulated a cesium atomic clock of 5071A with a clock difference of 10,000 s. The sampling time is 1 s.

Phase jump

Phase jump is one of the typical phase anomalies of atomic clocks. It shows that the mean value of the phase increases or decreases suddenly at a certain moment and lasts for a while. Therefore, it can be modeled as:

$${\tilde{x}}_1(k)=\{\begin{array}{cc}{x}_1(k),\hfill & k<n_0\hfill \\ \multicolumn{1}{c}{x_1(k)+\mathrm{\Delta }x,}& k\ge n_0\hfill \end{array}$$

(45)

where $n_0$ is the time instant at which the phase anomaly occurs, and $\mathrm{\Delta }x$ is the size of the jump.

Figure 4 shows the phase jump of size $\mathrm{\Delta }x=3\times {10}^{-10}$ occurring at 5000 s.

Figure 4.

Time error of a time-invariant clock.

Parameter analysis

From equations (32)–(35) and (42), it can be seen that the detection probability is related to the predicted steps and the accumulation numbers. Therefore, given that the probability of false alarm is 10⁻⁵ and the accumulated number M = 1, the relationship between the probability of detection (P_D) and the predicted steps is shown in Figure 5. The relationship between the time of detection and the predicted steps is shown in Figure 6. The time of abnormal occurrence is 5000 s, and the detection time represents the time when abnormal detection can be carried out normally. For example, if the detection time is 27 s, it means that an abnormity is detected in 5027 s and a warning is given.

Figure 5.

The relationship between detection probability and the number of prediction steps when M = 1.

Figure 6.

The relationship between detection time and the number of prediction steps when M = 1.

In Figure 5, the detection probability decreases gradually with the increase of the number of predicted steps, indicating that for phase detection, the number of predicted step p should not be too large.

In Figure 6, due to the small number of predicted steps, the anomaly detection statistic fails to reach the detection threshold, and the detection time is 0, indicating that no anomaly is detected. When the number of prediction steps is too large, the overall prediction residual becomes larger, and the abnormal detection statistic is submerged in the non-abnormal detection statistics, which also leads to the detection failure, and the detection time is 0.

Combined with Figures 5 and 6, it can be concluded that under the current conditions, it is appropriate to take the prediction step p = 10. At the same time, it can be seen that when p = 10, the detection probability is close to 1, and there is no need to accumulate epochs. At this time, M = 1 and the detection time is 27 s.

Detection performance

To make the theoretical simulation detection results more convincing, we adopted the M-C test to verify the theoretical detection performance. According to 3.3 and 4.1.1, the horizon is 50, the number of predicted steps is 10, and the number of the accumulated epoch is 1. For a set of different phase jumps, as shown in equation (47), the theoretical detection probabilities under a given set of different P _fa , as shown in equation (46), are analyzed. Meanwhile, the statistical results of 100 M-C tests are also shown in Figure 7.

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

(46)

$$\mathrm{\Delta }x=0.8\times {10}^{-10},1\times {10}^{-10},2.0\times {10}^{-10};$$

(47)

Figure 7.

Detection probabilities with given different phase jumps in equation (47).

It can be seen from Figure 7 that for a given set of phase jumps, the detection probability increases with the increase of the probability of false alarm. As the size of the phase jump increases, the detection probability also increases, which is consistent with the expected results. At the same time, it can be seen that the statistical results of the M-C test are in good agreement with the theoretical calculation results, which proves the correctness of the theoretical simulation results. Detailed statistical results are presented in Table 1.

Table 1.

The detection results under different phase jumps ( $N=50,p=10,M=1$ ).

${\mathrm{P}}_{\mathrm{fa}}$	10−7	5 × 10−7	10−6	5 × 10−6	10−5	10−4	10−3
$\mathrm{\Delta }x=0.8\times {10}^{-10}$	0.0610	0.1063	0.1332	0.2163	0.2620	0.4560	0.6878
$\mathrm{\Delta }x=1.0\times {10}^{-10}$	0.1458	0.2254	0.2678	0.3849	0.4423	0.6486	0.8369
$\mathrm{\Delta }x=2.0\times {10}^{-10}$	0.9994	0.9998	0.9999	0.9999	0.9999	0.9999	0.9999

Frequency jump

The frequency jump is one of the typical frequency anomalies of atomic clocks. It shows that the mean value of the normalized frequency deviation increases or decreases suddenly at a certain moment and lasts for a while. It can be modeled as:

$${\tilde{x}}_2(k)=\{\begin{array}{cc}{x}_2(k),\hfill & k<n_0\hfill \\ \multicolumn{1}{c}{x_2(k)+\mathrm{\Delta }f,}& k\ge n_0\hfill \end{array}$$

(48)

where $n_0$ is the time instant at which the frequency anomaly occurs, and $\mathrm{\Delta }f$ is the size of the jump.

Figure 8 shows the normalized frequency deviation jump of size $\mathrm{\Delta }f=3\times {10}^{-11}$ occurring at 5000 s.

Figure 8.

The normalized frequency deviation of a time-invariant clock.

Parameter analysis

Similar to phase jump, given that the probability of false alarm is 10^-5 and the relationship between detection probability and the number of accumulation under different prediction steps is shown in Figure 9. The relationship between detection time and the number of accumulation under different prediction steps is shown in Figure 10.

Figure 9.

The relationship between detection probability and the number of accumulation under different prediction steps.

Figure 10.

The relationship between detection time and the number of accumulation under different prediction steps.

As can be seen from Figure 9, under the same accumulation number, the larger the predicted step size is, the lower the detection probability is. Because the larger the number of predicted steps, the greater the estimation error. Under the same number of predicted steps, the detection probability increases with the increase of accumulation number. If the multi-epoch accumulation is not carried out like the phase detection, the detection probability is very low and no anomaly can be detected.

In Figure 10, if the number of predicted steps is too small, the abnormal detection statistic is too small to reach the detection threshold, and it is impossible to detect the anomaly by accumulating the epoch, just as p = 20 and the detection time is 0, which means that the detection failed. Therefore, for frequency anomaly detection, the number of predicted steps should be slightly larger. To make the detection time as small as possible and maintain a high detection probability, the predicted step is set as 70, the accumulated number is set as 25, and the detection time is 50 s at this time.

Detection performance

The simulation process is similar to the phase jump. For a set of different frequency jumps, as shown in equation (49), the theoretical detection probabilities under a given set of different P _fa , as shown in equation (46), are analyzed. Meanwhile, the statistical results of 100 M-C tests are also shown in Figure 11.

$$\mathrm{\Delta }f=1\times {10}^{-11},1.5\times {10}^{-11},2\times {10}^{-11};$$

(49)

Figure 11.

Detection probabilities with given different frequency jumps in equation (49).

It can be seen that the detection probability calculated by the theoretical formula is close to the statistical probability of 100 M-C experiments, with the same trend, which once again proves the correctness of the detection probability formula calculated by the theoretical formula. Detailed statistical results are presented in Table 2.

Table 2.

The detection results under different frequency jumps ( $N=100,p=70,M=25$ ).

${\mathrm{P}}_{\mathrm{fa}}$	10−7	5 × 10−7	10−6	5 × 10−6	10−5	10−4	10−3
$\mathrm{\Delta }f=1.0\times {10}^{-11}$	0.0018	0.0044	0.0064	0.0146	0.0207	0.0624	0.1692
$\mathrm{\Delta }f=1.5\times {10}^{-11}$	0.1291	0.1973	0.2339	0.3365	0.3880	0.5817	0.7788
$\mathrm{\Delta }f=2.0\times {10}^{-11}$	0.7728	0.8469	0.8742	0.9259	0.9431	0.9805	0.9955

Real data analysis

To verify the usability of the new method, we set up an experimental platform, and the experimental principle block diagram and part of the tested scenes are shown in Figure 12. The measured data of a cesium atomic clock of 5071A are obtained.

Figure 12.

Experimental principle block diagram and part of the tested scenes.

The whole experimental platform is set up in a time-keeping laboratory, and a cesium atomic clock of 5071A is used as the atomic clock to be measured. One of its outputs 10 MHz signal as input to a High-Resolution Phase and Frequency Offset Generator(HROG-10), as a branch to be tested. The other output signal is direct as the input of the Multi-channel Measurement System (MMS). According to the user manual of MMS, its Allan deviation is less than 2.5e−13 at 1 s. Therefore, it is feasible to use it as an accurate measuring device. VCH-1003M is the main clock of the whole timekeeping system. Its 10 MHz output signal serves as the reference signal of MMS, and the 1PPS signal serves as the initial phase synchronization signal of MMS. The Allan deviation analysis of the reference clock and the atomic clock to be measured is given in Figure 13.

Figure 13.

The Allan deviation of the reference and the atomic clock to be measured.

The HROG-10 is a device that can be manually set with different phase or frequency jumps, and the whole experiment is to configure different phase or frequency jumps through it.

Data were collected on December 9, 2020, and different frequency jumps $\mathrm{\Delta }f=2\times {10}^{-11},1\times {10}^{-11},3\times {10}^{-11}$ , phase jump $\mathrm{\Delta }x=3\times {10}^{-10}$ were respectively set through the HROG-10, as shown in Figure 14. Data is collected every second. The whole test process was carried out under constant temperature and humidity. During the test, to ensure that the measured data will not be interfered with by other factors, we try our best to close the laboratory and provide high-power UPS emergency power to ensure the continuity of data.

Figure 14.

Real data of the normalized frequency deviation with phase and frequency anomalies.

In Figure 14, the blue line describes the abnormal normalized frequency deviation data collected after HROG-10, the red line represents the estimated value of UFIR, and the yellow line is the estimated value using the Kalman extrapolation-accumulation method. ³⁴ As can be seen from the figure, at the time of frequency jump, the estimation error of UFIR is greater than that of Kalman, which is conducive to the construction of anomaly detection statistics. Figure 15 shows the root mean square error (RMSE) of the two methods. The statistical RMSE of the new method is 8.169e−12, and the statistical RMSE of the Kalman method is 6.487e−12. It can be seen that the consistency of the two methods is very good, and the new method is slightly inferior to the Kalman method. Because Kalman algorithm guarantees the minimum root mean square error, but UFIR cannot.

Figure 15.

Consistency of the proposed method and Kalman method estimates of real data.

The performance of the new method and the Kalman filter method proposed in Yan et al. ³⁴ are both used for the analysis of measured data.

In Figure 16, the blue line represents the detection statistics calculated by equation (35) of the new method, while the red line represents the detection statistics calculated by the improved the Kalman filter method. For comparison, the measured normalized frequency deviation of the atomic clock multiplied by 5 × 10⁻¹¹ is superimposed on the detection statistics, as described in the purple line.

Figure 16.

Detection results of real data by the proposed method and Kalman.

The yellow dotted line corresponds to the detection threshold. According to equation (41), the detection threshold is related to P _fa and accumulation number M, while the detection threshold of measured data is calculated at P_fa = 10⁻⁵ and M = 25.

To compare the effectiveness of the detection, the detection statistics of the normal data without anomalies are plotted in Figure 17. It can be seen that:

Figure 17.

Detection results of normal real data by the proposed method and Kalman filter.

1. The value of detection statistics increases when the frequency jumps or the phase jumps and exceeds the detection threshold. The results show that the new method can be used for atomic clock anomaly detection.

2. Both the new method and the improved the Kalman filter method are effective in detecting phase or frequency jumps, but the new method is more sensitive to weak jumps, which is represented by larger detection statistics.

3. The proposed method can detect both weak frequency and phase jump simultaneously, while the Kalman filter method cannot.

4. The detection statistics of the new method are larger than that of the Kalman filter method at every time, but the normal time is not greater than the detection threshold, and the abnormal time is far more than the detection threshold. This is beneficial for anomaly detection.

Conclusions

This paper presents a new anomaly detection algorithm for the atomic clock. P-step Kalman-like UFIR algorithm breaks through the shortcomings of the Kalman filter method which relies on accurate initialization of state noise matrix and error matrix. This advantage is very useful for real-world applications. Although it does not guarantee the minimum estimation error as Kalman does, it is very practical. Through the analysis of the simulation data and the measured data, it is proved that the new method can detect the anomaly of the atomic clock accurately and timely.

The correctness of the proposed method is also verified by the M-C simulation experiment.