1. INTRODUCTION
Atomic clocks are of vital importance to Global Navigation Satellite Systems (GNSS), and the performance of the atomic clocks ensures GNSS users have precise position and time information (Lee et al., Reference Lee, Kim and Lee2011). Most of the studies into navigation system integrity monitoring are related to satellite integrity and not specifically to that of the embedded clock, while clock anomalies are the most prevalent source of space segment anomalies and the most common source of major service anomalies. These anomalies can result in thousands of metres of range error (Kaplan and Hegarty, Reference Kaplan and Hegarty2006). Rapid detection of anomalies is therefore of great importance.
Several detection methods have been proposed to detect and identify clock anomalies (Riley, Reference Riley2008; Nunzi et al., Reference Nunzi, Galleani, Tavella and Carbone2007; Nunzi and Carbone, Reference Nunzi and Carbone2008; Galleani, Reference Galleani2008a; Reference Galleani2010; Galleani and Tavella, Reference Galleani and Tavella2008; Reference Galleani and Tavella2009; Sesia et al., Reference Sesia, Galleani and Tavella2011). The Dynamic Allan Variance (DAVAR), a representation of the time-varying stability of an atomic clock, is used to detect and identify clock anomalies and a fast detector by means of the Generalised Likelihood Ratio Test (GLRT) is adopted. These two techniques are complementary, since the GLRT allows a fast detection of the fault while the DAVAR gives an identification of the fault type. Another useful method based on the prediction interval approach using the one-step predicted residue of the Kalman filter was proposed in Galleani and Tavella (Reference Galleani and Tavella2012) and Lee et al. (Reference Lee, Kim, Jeong and Lee2011). This method detects clock anomalies quickly but demands filtering computation every measurement period. It would be advantageous to provide new ways to detect atomic clock anomalies with tolerable computational complexity.
In this paper, an algorithm for detection of two typical anomalies of atomic clocks, namely, the phase jump and the frequency jump, is developed. A detection statistic based on parity transformation, derived as Receiver Autonomous Integrity Monitoring (RAIM) (Brown, Reference Brown1992; Lee, Reference Lee1992; Kelly, Reference Kelly1996), is constructed to detect clock anomalies in real time. The concept of the optimal accumulation number is provided, and a method to find out an optimal accumulation number is given.
The paper is organised as follows. In Section 2 an extended measurement model is achieved by the parity relation-based fault detection approach using temporal redundancy. In Section 3 an integrity monitoring algorithm based on an extended measurement model is proposed, and the option of the optimal accumulation number is also discussed. Finally, in Section 4 the algorithm is applied to simulated data and its validity is evaluated for two typical anomalies.
2. THE EXTENDED MEASUREMENT MODEL
The two-state clock model is chosen here, which is expressed as follows:
$${\bi x}\lpar k+1\rpar ={\bf A}\cdot {\bi x}\lpar k\rpar +{\bf B}\cdot \xi \lpar k\rpar $$
where 
${\bf A}=\left[\matrix{1 & {T_{s} } \cr 0 & 1 \cr}\right]\comma \; \xi \lpar k\rpar =\left[\matrix{\xi_{1} \lpar k\rpar \cr \xi_{2} \lpar k\rpar \cr }\right]\comma \; {\bi x}\lpar k\rpar =\left[\matrix{x_{1} \lpar k\rpar \cr x_{2} \lpar k\rpar \cr }\right]\comma \; {\bf B}={\bf I}_{2\times 2}$. I2×2 is the 2-by-2 identity matrix. x 1(k) is the phase deviation and x 2(k) is only a component of the frequency deviation, which is defined as the derivative of the phase deviation x 1(k). T s represents the sampling interval. The process noise vector ξ (k) has two components: ξ 1 (k) represents the random variation due to White Frequency Noise (WFN), and similarly ξ 2 (k) accounts for the random variation due to Random Walk Frequency Noise (RWFN), respectively (Tavella, Reference Tavella2008; Galleani, Reference Galleani2008b). ξ (k) is a multivariate Gaussian random vector with zero mean and a covariance of Q, which can be expressed as:
$${\rm E}\lpar \xi \lpar k\rpar \rpar =0$$
$${\rm E}\lpar \xi \lpar k\rpar \cdot \xi \lpar k\rpar ^{\rm T}\rpar ={\bf Q}=\left[\matrix{{q_1 T_s +q_2 \displaystyle{{T_s^3 }\over{3}}} & {q_2 \displaystyle{{T_s^2 }\over{2}}} \cr {q_2 \displaystyle{{T_s^2 }\over{2}}} & {q_2 T_s } \cr } \right]$$where q 1 and q 2 are noise parameters of the covariance matrix Q. The parameter q 1 drives the white frequency noise, corresponding to a random walk on the time deviation, while q 2 drives the frequency random walk, corresponding to an integrated random walk on the time deviation (Tavella, Reference Tavella2008; Galleani, Reference Galleani2008a; Reference Galleani2008b).
The clock phase measurement model can be expressed as:
$$z\lpar k\rpar ={\bf H}\cdot {\bi x}\lpar k\rpar +w\lpar k\rpar $$
where 
${\bf H}^{\rm T}=\left[\matrix{1 \cr 0 \cr }\right]$, z(k) is the measured clock phase deviation and w(k) is White Gaussian Noise (WGN) with variance σ 2.
An extended measurement model is achieved by the parity relation-based fault detection approach using the temporal redundancy according to Ding et al. (Reference Ding, Guo and Jeinsch1999). The extended measurement model is the direct extension of the original measurement model given as Equation (4). It is similar to the batch approach, which is applied to a sequence of measurements and system dynamics over a finite time window. The extended measurement model generated using this approach is expressed by:
$${\bf Z}\lpar k\rpar ={\bf H}_{{\bf A}} \cdot {\bi x}\lpar k-{\rm N}+1\rpar +{\bf F}\cdot {\brGamma} \lpar k\rpar +{\bf G}\cdot {\bf W}\lpar k\rpar $$where N is the order of the parity relation.
The corresponding coefficient matrix of the extended measurement model is given as:
$${\bf Z}\lpar k\rpar = \left[\matrix{{z\lpar k-N+1\rpar } \cr {z\lpar k-N+2\rpar } \cr \vdots \cr {z\lpar k\rpar } \cr }\right]$$
$${\brGamma} \lpar k\rpar = \left[\matrix{{\xi \lpar k-N+1\rpar } \cr {\xi \lpar k-N+2\rpar } \cr \vdots \cr {\xi \lpar k\rpar } \cr }\right]=\left[\matrix{ \xi_1 \lpar k-N+1\rpar \cr \xi_2 \lpar k-N+1\rpar \cr \xi_1 \lpar k-N+2\rpar \cr \xi_2 \lpar k-N+2\rpar \cr \vdots \cr \xi_1 \lpar k\rpar \cr \xi_2 \lpar k\rpar \cr }\right]$$
$${\bf W}\lpar k\rpar = \left[\matrix{ {w\lpar k-N+1\rpar } \cr {w\lpar k-N+2\rpar } \cr \vdots \cr {w\lpar k\rpar } \cr }\right]$$
$${\bf H}_{\bf A} =\left[\matrix{{\bf H} \cr {\bf HA} \cr \vdots \cr {\bf HA}^{{\rm N}-1}}\right]$$
$${\bf F} = \left[\matrix{0 & 0 & \cdots & 0 & 0 \cr {\bf HB} & 0 & \cdots & 0 & 0 \cr {\bf HAB} & {\bf HB} & \cdots & 0 & 0 \cr \vdots & \vdots & \vdots & \vdots & \vdots \cr {\bf HA}^{{\rm N}-2}{\bf B} & {\bf HA}^{{\rm N}-3}{\bf B} & \cdots & {\bf HB} & 0 \cr }\right]$$
$${\bf G} = {\bf I}_{{\rm N}\times {\rm N}} $$The extended measurement equation can be rewritten as:
$${\bf Z}\lpar k\rpar ={\bf H}_{\bf A} \cdot {\bi x}\lpar k-N+1\rpar +{\brPsi} \lpar k\rpar $$where
$${\brPsi} \lpar k\rpar =\left[\matrix{ {w\lpar k-N+1\rpar } \cr {\bf HB}\cdot \xi \lpar k-N+1\rpar +w\lpar k-N+2\rpar \cr \vdots \cr \sum\limits_{i=0}^{{\rm N}-1} {\bf HA}^i{{\bf B}}\cdot \xi \lpar k-i\rpar +w\lpar k\rpar }\right]\comma \;$$which follows a generalised Gaussian distribution. This can be expressed as follows:
$${\brPsi} \lpar k\rpar \sim N\lpar {\bf 0}\comma \; {\bf C}\rpar $$The covariance matrix C is expressed as:
$${\bf C}=\left[\matrix{{\sigma_{1\comma 1}^2 } & \cdots & {\sigma_{1\comma {\rm N}}^2 } \cr \vdots & \ddots & \vdots \cr {\sigma_{{\rm N}\comma 1}^2 } & \cdots & {\sigma_{{\rm N}\comma {\rm N}}^2 }}\right]$$Elements of matrix C are given as follows (the calculation process is provided in detail in the Appendix):
$$\sigma _{{\rm n}\comma {\rm n}}^2 = \left\{\matrix{\sigma^2\comma \; \hfill & {\rm n}=1 \cr \sum\limits_{i=0}^{{\rm n}-2} {\bf HA}^{i} {\bf BQB}^{\rm T}\lpar {\bf A}^i\rpar^{\rm T}{\bf H}^{\rm T}+\sigma^2\comma \; & {\rm n}>1}\right. $$
$$\sigma_{\rm m\comma n}^2 =0\comma \; \quad {\rm m}=1\hbox{ or }{\rm n}=1 $$
$$\sigma_{\rm m\comma n}^2 = \left\{\matrix{{\bf H}\left(\sum\limits_{\mathop {i = 1}\limits_{m \ge 2} }^{m-2} {\bf A}^{{\rm n}-{\rm m}+i}{\bf BQB}^{\rm T}\lpar {\bf A}^{\rm T}\rpar^i\right){\bf H}^{\rm T}\comma \; & {\rm m}< {\rm n} \cr {\bf H}\left(\sum\limits_{\mathop {i = 1}\limits_{n \ge 2} }^{n-2} {\bf A}^{{\rm m}-{\rm n}+i}{\bf BQB}^{\rm T}\lpar {\bf A}^{\rm T}\rpar^i\right){\bf H}^{\rm T}\comma \; & {\rm m}>{\rm n}}\right. $$As elements of the extended measurement vector are correlated, it is a good choice to pre-whiten the noise vector. For any C that is positive definite, it can be shown that C−1 exists and is also positive definite. Consequently, it may be factored as C−1 = DTD, where D is called the pre-whitening matrix, which produces WGN (Kay, Reference Kay1998). The pre-whitening matrix D is achieved by Cholesky factorisation.
The pre-whitened extended measurement model can be expressed as:
$${\bf DZ}\lpar k\rpar ={\bf DH}_{\bf A} \cdot {\bi x}\lpar k-{\bf N}+1\rpar +{\bf D}{\brPsi} \lpar k\rpar $$Equation (18) can be also rewritten as:
$$\tilde{{\bf Z}}\lpar k\rpar =\tilde{{\bf H}}\cdot {\bi x}\lpar k-N+1\rpar +\tilde{{\brPsi}}\lpar k\rpar $$
where 
$\tilde{{\bf Z}}\lpar k\rpar ={\bf DZ}\lpar k\rpar \comma \; \tilde{{\bf H}}={\bf DH}_{\bf A}\comma \; \tilde{{\brPsi}}\lpar k\rpar ={\bf D}{\brPsi}\lpar k\rpar $.
Therefore, the extended measurement model for an atomic clock is built. It can be seen from the extended measurement model that the new measurements not only include the satellite clock observations, but also the dynamic information of the satellite clock. This is the highlight of the extended measurement model.
By adding a third component, an atomic clock model considering the frequency drift is built, which can be used to model rubidium and maser clocks. In particular, the vector and matrices used in Equation (1) become:
$${\bf A}=\left[\matrix{1 & {T_s } & {\displaystyle{1 \over 2}T_s^2 } \cr 0 & 1 & {T_s } \cr 0 & 0 & 1 \cr }\right]\comma \; \xi \lpar k\rpar =\left[\matrix{\xi_1 \lpar k\rpar \cr \xi_2 \lpar k\rpar \cr \xi_3 \lpar k\rpar \cr }\right]\comma \; {\bi x}\lpar k\rpar =\left[\matrix{x_1 \lpar k\rpar \cr x_2 \lpar k\rpar \cr x_3 \lpar k\rpar \cr }\right]$$while the matrix Q defined in Equation (3) is given by:
$${\bf Q}=\left[\matrix{{q_1 T_s +q_2 \displaystyle{{T_s^3 }\over{3}}+q_3 \displaystyle{{T_s^5 }\over{20}}} & {q_2 \displaystyle{{T_s^2 }\over{2}}+q_3 \displaystyle{{T_s^4 }\over{8}}} & {q_3 \displaystyle{{T_s^3 }\over{6}}} \cr {q_2 \displaystyle{{T_s^2 }\over{2}}+q_3 \displaystyle{{T_s^4 }\over{8}}} & {q_2 T_s +q_3 \displaystyle{{T_s^3 }\over{3}}} & {q_3 \displaystyle{{T_s^2 }\over{2}}} \cr \noalign{\vskip6pt} {q_3 \displaystyle{{T_s^3 }\over{6}}} & {q_3 \displaystyle{{T_s^2 }\over{2}}} & {q_3 T_s } \cr } \right]$$3. THE INTEGRITY MONITORING ALGORITHM
3.1. Detection
Detection of a clock anomaly is achieved by parity transformation similar to RAIM (Brown, Reference Brown1992; Lee, Reference Lee1992; Kelly, Reference Kelly1996).
Performing a transformation on 
$\tilde{{\bf Z}}\lpar k\rpar $, the parity vector p can be found and this is expressed as:
$${\bf p}={\bf P}\tilde{{\bf Z}}\lpar k\rpar $$
where the parity transformation matrix P is defined as an (N-2)-by-N matrix, which can be obtained by QR factorisation of the 
$\tilde{{\bf H}}$ matrix. The rows of P are mutually orthogonal, unity in magnitude, and mutually orthogonal to the columns of 
$\tilde{{\bf H}}$.
If there is no fault, then the parity vector p can be given as follows:
$${\bf p}={\bf P}\cdot \tilde {{\bf Z}}\lpar k\rpar ={\bf P}\cdot \tilde {{\brPsi} }\lpar k\rpar $$ 
$\tilde {{\brPsi} }\lpar k\rpar $ is WGN with variance 1 according to Equations (18) and (19), and the rows of P are mutually orthogonal and unity in magnitude, thus 
${\bf P}\cdot \tilde {{\brPsi} }\lpar k\rpar $ follows a N-2 degrees of freedom WGN distribution with variance 1.
In the presence of faults, the measurement model has the form as follows:
$$\tilde {{\bf Z}}\lpar k\rpar =\tilde {{\bf H}}\cdot {\bi x}\lpar k-N+1\rpar +{\bf L}\lpar k\rpar +\tilde {{\brPsi}}\lpar k\rpar $$
where 
${\bf L}\lpar k\rpar ={\bf D}\cdot \left(\matrix{{{\rm L}_{1} } & {{\rm L}_{2} } & \cdots & {{\rm L}_{\rm N}}}\right)^{\rm T}$ and Li is the anomalous value of clock phase deviation.
The parity vector p with clock faults can be given as:
$${\bf p}={\bf P}\cdot {\bf L}\lpar k\rpar +{\bf P}\cdot \tilde {{\brPsi} }\lpar k\rpar $$which follows a N-2 degrees of freedom WGN distribution with variance 1, and non-centrality parameter P · L(k).
In order to detect clock anomalies, a detection statistic is adopted:
$${\rm FD}={\bf p}^{\rm T}\cdot {\bf p} $$It can be concluded from Equations (23) and (25) that the proposed detection statistic FD follows a chi-square distribution with N-2 degrees of freedom when no clock anomaly exists. It has a non-central chi-square distribution with N-2 degrees of freedom and non-centrality parameter λ (function of the fault vector P · L(k)) when a clock anomaly occurs.
In this paper, H1 and H0 are defined as the statistical hypotheses of clock anomalies existing or not existing, including phase and frequency jumps, respectively.
Assuming the probability of false alarm is P fa, and the probability of detection is P d, the detection threshold can be confirmed by false alarm as follows. Details of calculation can be found in Kay (Reference Kay1998).
$$\eqalign{{\rm P}_{fa} & =\Pr \lcub {\rm FD}>\gamma \semicolon \; H_0 \rcub \cr & =\vint_\gamma^\infty \left(2^{\displaystyle{{{\rm N}-2}\over{2}}}\Gamma \left(\displaystyle{{{\rm N}-2}\over{2}}\right)\right)^{-1} {\rm FD}^{\displaystyle{{{\rm N}-2}\over{2}}-1}e^{-\displaystyle{{{\rm FD}}\over{2}}} d\lpar {\rm FD}\rpar } $$The probability of detection is given as:
$$\eqalign{{\rm P}_d &=\Pr \lcub {\rm FD}>\gamma \semicolon \; H_1 \rcub \cr & = \vint_\gamma ^\infty \displaystyle{{1}\over{2}}\left(\displaystyle{{{\rm FD}}\over{\lambda}}\right)^{\displaystyle{{{\rm N}-2}\over{4}}-\displaystyle{{1}\over{2}}} e^{-\displaystyle{{{\rm F}D+\lambda }\over{2}}}I_{\displaystyle{{{\rm N}-2}\over{2}}-1} \lpar \sqrt {{\rm FD}\lambda}\rpar d\lpar {\rm FD}\rpar } $$where the non-centrality parameter λ equals:
$$\lambda =\lpar {\bf P}\cdot {\bf L}\lpar k\rpar \rpar ^{\rm T}\cdot\, \lpar {\bf P}\cdot {\bf L}\lpar k\rpar \rpar $$Performance levels of the detector can be obtained for a given probability of false alarm according to Equations (27) and (28).
It can be seen that detection performance of the detector is decided by the non-centrality parameter λ, which is related to values of fault vector L(k), the accumulation number N and the parity transformation matrix P.
The key idea of the proposed detection method is to detect the discrepancy between the normal and abnormal measurements in a finite window: if the discrepancy denoted by FD is larger than the threshold, an anomaly is detected.
3.2. The optimal accumulation number
Since the accumulation number means the degree of freedom for the distribution of the proposed detection statistic, it is necessary to show the relationship between the accumulation number and the detection statistic and to study the influence of the accumulation number on the detection performance. On one hand, increasing the accumulation number leads to an increase of anomalous observations that can be used and so gives an improvement of the performance and on the other hand, more observations mean more measurement noise in the prediction residues, which might lead to a decrement of performance when the accumulation number increases. Consequently, an optimal accumulation number might exist.
The optimal accumulation number is defined in this paper as the accumulation number from among all of those for which the target of detecting a given weak anomaly is met or exceeded, for which the accumulation number is a minimum. It is the objective of this section to find out the possible optimal accumulation number for a given anomaly.
To find out an optimal accumulation number in terms of the detection performance, a method is provided, whose detailed implementation steps can be described as follows:
Step 1: Set target probabilities of detection and false alarm, P d and P fa, respectively, and confirm the theoretical non-centrality parameter λ with accumulation number n according to Equations (27) and (28). The accumulation number n is first initialised to 1.
Step 2: Calculate the actual non-centrality parameter
$\tilde {\lambda }$ with accumulation number n for the given anomaly according to Equation (29).Step 3: Compare the theoretical non-centrality parameter λ against the actual non-centrality parameter
$\tilde {\lambda }$: if the difference is bigger than a given tolerance, end; otherwise set n = n + 1 and go back to Step 2.
The optimal accumulation number can be also confirmed by simulation, which is an effective way to avoid heavy computation.
A situation that the actual non-centrality parameter 
$\tilde {\lambda }$ never exceeds the theoretical λ might happen. Target detection performance could not be achieved in this situation due to the limitation of the algorithm. A reasonable treatment is to choose the accumulation number, from among all of the possible numbers, for which the probability of detection is best with a given probability of false alarm.
Set target probabilities of detection and false alarm, P d = 0 · 99 and P fa = 10−5 respectively, and the theoretical non-centrality parameter λ under different accumulation numbers are calculated according to Equations (27) and (28), which are shown in Figure 1. The theoretical non-centrality parameter λ can be used as a reference value, and if the actual non-centrality parameter exceeds it, then the optimal accumulation number is found.
3.3. Implementation details
The detailed implementation steps of the proposed algorithm can be described as follows:
Step 1: Select N, the order of the parity relation, and then construct the measurement vector Z(k).
Step 2: Compute the pre-whitening matrix D and parity transformation matrix P.
Step 3: Compute values of the proposed detection statistic FD, which can be achieved by FD = |p|2 = |P · D · Z(k)|2.
Step 4: Each value of the detection statistic is compared with the detection threshold, which can be confirmed by Equation (27) with a given probability of false alarm. If the value of the detection statistic exceeds the threshold, a clock anomaly is detected.
The proposed integrity monitoring algorithm is illustrated schematically in Figure 2.
Figure 2. Implementation scheme of the new method.
4. NUMERICAL SIMULATIONS
To evaluate the proposed method, numerical simulation of phase and frequency jump anomalies for an atomic clock with implementation structure as shown in Figure 2 were performed. The case of the two-state model was considered, whose noise parameters for a sampling time T s = 1 s are given by:
$$q_1 =\lpar 5\times 10^{-12}\rpar ^2 $$
$$q_2 = \lpar 1\times 10^{-16}\rpar ^2 $$The initial value of the clock state vector is given as:
$$\left[{\matrix{ x_1 \lpar 0\rpar \cr x_2 \lpar 0\rpar \cr }} \right]=\left[{\matrix{ 0 \cr 0 \cr }} \right]$$The variance value of measurement noise for a sampling time T s = 1 s is given as:
$$\sigma ^2=\lpar 5\times 10^{-12}\rpar ^2 $$A time-invariant clock was simulated, that is a clock whose parameters given in Equations (30), (31), (32) and (33) do not change over time. In Figure 3, an illustration of a measured clock phase deviation sequence with length 1,000 seconds is shown.
Two typical anomalies of atomic clocks are discussed in the following, including phase jump and frequency jump. Monte-Carlo (M-C) simulation is also conducted to verify the performance of the proposed method.
4.1. Phase jump
A typical clock anomaly is the phase jump. It corresponds to a sudden variation in the mean value of the normalised phase deviation and consequently can be modelled as:
$$\tilde {x}\lpar k\rpar = \left\{\matrix{x\lpar k\rpar \comma \; \hfill & k < n_0 \cr x\lpar k\rpar +\Delta x\comma \; & k\ge n_0}\right. $$where n 0 is the time instant at which the phase jump occurs, and Δ x is its size.
4.1.1. Validity of the method
A phase jump anomaly with a jump size of 6 × 10−11 s is simulated. The time instant at which the phase jump occurs is at 500 seconds. The corresponding proposed test statistic computed with the order of the parity relation N = 6 is given in Figure 4. Assuming the probability of false alarm is constant with a value of 10−6, the detection threshold can be confirmed by Equation (27), which is shown as the red line in Figure 4. As can be seen from the figure, that phase jump is detected in real time.
Figure 4. Values of the proposed test statistic for phase jump.
4.1.2. Detection performance
To verify the detection performance of the new method, M-C simulation is adopted. The order of the parity relation is assumed to be six. Figure 5 shows the results of 10,000 repeated M-C simulations. The probability of detection for values of P fa given in Equation (35), under different phase jumps with the sizes given in Equation (36) are calculated, which are also shown in Figure 5.
Figure 5. Receiver Operating Characteristic (ROC) curves for M-C simulations and theoretical calculation.
$$P_{fa} = 10^{-6}\comma \; 10^{-5}\comma \; 10^{-4}\comma \; 10^{-3} $$
$$\Delta x = 5\comma \; 6\comma \; 7\times 10^{-11}\, {\rm S} $$The curves shown in Figure 5 are known as Receiver Operating Characteristics (ROC) (Kay, Reference Kay1998). It can be seen that the M-C simulation results are highly consistent with results of the calculation formula, which illustrates the correctness of the calculation formula for detection probability. Moreover, for a given value of phase jump, the detection probability increases with the probability of false alarm. As expected, as the size of the phase jumps increases, the detector performance improves.
4.1.3. Optimal accumulation number
Numerical simulation is adopted to verify the relationship between detection performance and accumulation number. It is assumed that the size of the phase jump is 6 × 10−11 s and the target probabilities of detection and false alarm are 0·95 and 10−6 respectively. The detection probability values with given value of P fa and size of phase jump as a function of the accumulation number are shown in Figure 6.
Figure 6. Detection probability values with given value of Pfa and size of phase jump as a function of the accumulation number.
It can be seen that detection performance improves at first as the accumulation number increases. It becomes worse instead of better when the accumulation number exceeds a specific number, and the maximum probability of detection is close to 0·81, less than the target. That is to say a situation in that the actual non-centrality parameter never exceeds the theoretical one happens. Therefore, the optimal accumulation number should be chosen as the one that gives the best probability of detection for a given probability of false alarm. As can be seen, the optimal accumulation number is 11 in the given situation.
4.2. Frequency jump
A type of frequency anomaly known as frequency jump is now considered. This corresponds to a sudden variation in the mean value of the normalised frequency, and can be modelled as:
$$\tilde {y}\lpar k\rpar = \left\{\matrix{y\lpar k\rpar \comma \; \hfill & k< n_0 \cr y\lpar k\rpar +\Delta f\comma \; & k\ge n_0}\right. $$where n 0 is the time instant at which the frequency jump occurs and Δ f is its size.
A frequency jump anomaly with jump size 2 × 10−11 Hz is simulated. The time instant at which the frequency jump occurs is at 500 seconds. The corresponding proposed test statistic computed with the order of the parity relation N = 6 is given in Figure 7. It is assumed that the probability of false alarm is constant with a value of 10−6 and the detection threshold can be confirmed by Equation (27). As can be seen from Figure 7, that frequency jump is effectively detected.
Figure 7. Values of the proposed test statistic for frequency jump.
Numerical simulation is adopted in order to verify the relationship between detection performance and the accumulation number. It is assumed that the size of the frequency jump is 2 × 10−11 Hz, and the target probabilities of detection and false alarm are 0·95 and 10−6 respectively. Results of 10,000 repeated M-C simulations are provided. The detection probability values with given value of P fa and size of frequency jump as a function of the accumulation number are shown in Figure 8.
Figure 8. Detection probability values with given value of Pfa and size of frequency jump as a function of the accumulation number.
It can be seen that the detection performance improves with the accumulation number. Moreover, when the accumulation number exceeds 20, the probability of detection is above 0·95. The optimal accumulation number should be chosen to ensure the target probability of detection is achieved. As can be seen, the favourable accumulation number is 21 in the given situation.
4.3. Performance comparison
Detection performances of the proposed method and the Kalman filter method proposed in Lee et al. (Reference Lee, Kim, Jeong and Lee2011), which is a useful clock anomaly detection method, are simulated. Values of Pd as a function of the sizes of phase and frequency jumps, with a constant Pfa value, are shown in Figures 9 and 10 respectively.
Figure 9. Phase jump detection performances of the proposed and Kalman filter methods.
Figure 10. Frequency jump detection performances of the proposed and Kalman filter methods.
The accumulation number is 8, and the probability of false alarm is 10−6. It is observed that the phase jump detection performance of the new method is poorer than that of the Kalman filter method, while its frequency jump detection performance is better than that of the Kalman filter method. This is related to the number of measurements each method used.
The detector of the Kalman filter method only used one-step prediction (that means only one measurement is used), which is sensitive to phase jump as the phase jump is directly reflected by the phase deviation measurements, while it is less sensitive to the frequency jump because the frequency jump would be reflected slowly by the phase deviation measurements.
In contrast, the new method uses more than one measurement and will be helpful in detecting frequency jumps. However, as more measurements mean more measurement noise, the performance of phase jump detection will be not as good as one-step prediction.
5. REAL DATA ANALYSIS
The frequency deviation of the GPS PRN 16 satellite clock is illustrated in Figure 11(a), for the period 29 August – 14 November 2010, obtained from the International GNSS Service (IGS) database. The sampling interval is 15 minutes. Two obvious frequency jumps are recognised.
Figure 11. Detection Results of Real Data.
Figure 11(b) and 11(c) show the corresponding values of the test statistics of the proposed method and the Kalman filter method. It can be seen that:
(1) The value of the test statistic increases at the time value corresponding to the frequency jump. This result demonstrates that the method can be used to detect clock anomalies.
(2) Both the new approach and the Kalman filter approach are effective in detecting frequency jumps, while the new approach is more sensitive as the detection values are larger.
6. CONCLUSIONS
An algorithm for atomic clock integrity monitoring is proposed. Numerical simulation of two frequent atomic clock anomalies, namely, the phase jump and the frequency jump, were performed, which verified the validity of the proposed algorithm for detecting clock anomalies. M-C simulation was also adopted to verify detection performance of the proposed method.
Such a method is simple to implement with tolerable computational complexity, and it can detect two main anomalies of atomic clocks and could provide an effective solution to GNSS satellite clock autonomous integrity monitoring.
APPENDIX: EXPRESSION OF THE COVARIANCE MATRIX C
The covariance matrix of Ψ(k) is expressed as
$$\eqalign{\textbf{C}=E\{\boldsymbol{\Psi} (k)\cdot \boldsymbol{\Psi} (k)^{\rm T}\} \cr &=\left[\vphantom{\begin{matrix}{E\left\{w(k-N+1)\cdot \left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)}+w(k)\right)^{\rm T}\right\}} \cr \vdots \cr {E\left\{\left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)}+w(k)\right)\cdot \left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi(k-i)} +w(k)\right)^{\rm T}\right\}}\end{matrix}}\begin{matrix} {E\{w(k-N+1)\cdot w(k-N+1)^{\rm T}\}} & \cdots \cr \vdots & \ddots \cr {E\left\{\left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)}+w(k)\right)\cdot w(k-N+1)^{\rm T}\right\}} & \cdots\end{matrix}\right.\notag\\&\quad\left.\begin{matrix}{E\left\{w(k-N+1)\cdot\left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)}+w(k)\right)^{\rm T}\right\}} \cr \vdots \cr {E\left\{\left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)}+w(k)\right)\cdot \left(\displaystyle\sum\limits_{i=0}^{{\rm N}-1} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi(k-i)} +w(k)\right)^{\rm T}\right\}}\end{matrix}\right]}$$The elements can be calculated as follows:
$$\sigma _{\rm m,n}^2 =\begin{cases} E\{w(k-N+1)\cdot w(k-N+1)^{\rm T}\}=\sigma ^2, & {\rm m}={\rm n}=1 \cr E\left\{\left(\displaystyle\sum\limits_{i=0}^{{\rm n}-2} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)} +w(k)\right)\cdot \left(\displaystyle\sum\limits_{i=0}^{{\rm n}-2} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)} +w(k)\right)^{\rm T}\right\} & \cr \qquad\quad =\displaystyle\sum\limits_{i=0}^{{\rm n}-2} {\textbf{HA}^i} \textbf{BQB}^{\rm T}(A^i)^{\rm T}\textbf{H}^{\rm T}+\sigma^2, & {\rm m}={\rm n}>1 \cr 0, & {\rm m}=1\hbox{ or }{\rm n}=1 \cr E\left\{\left(\displaystyle\sum\limits_{i=0}^{{\rm m}-2} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)} +w(k)\right)\cdot \left(\displaystyle\sum\limits_{i=0}^{{\rm n}-2} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)} +w(k)\right)^{\rm T}\right\}\\ \qquad\quad =\textbf{H}\left(\displaystyle\sum\limits_{\substack{i=1 \cr m\ge 2}}^{{\rm m}-2} \textbf{A}^{{\rm n}-{\rm m}+i}\textbf{BQB}^{\rm T}(\textbf{A}^{\rm T})^i\right)\textbf{H}^{\rm T}, & {\rm m}<{\rm n} \cr E\left\{\left(\displaystyle\sum\limits_{i=0}^{{\rm n}-2} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)} +w(k)\right)\cdot \left(\displaystyle\sum\limits_{i=0}^{{\rm m}-2} {\textbf{H}\textbf{A}^i{\textbf{B}}\cdot \xi (k-i)} +w(k)\right)^{\rm T}\right\}\\ \qquad \quad =\textbf{H}\left(\displaystyle\sum\limits_{\substack{i=1 \cr {\rm n}\ge 2}}^{{\rm n}-2} {\textbf{A}^{{\rm m}-{\rm n}+i}\textbf{B}} \textbf{QB}^{\rm T}(\textbf{A}^{\rm T})^i\right)\textbf{H}^{\rm T}, & {\rm m}>{\rm n}\end{cases}$$Thus, elements of matrix C can be expressed as follows, which are shown as Equations (15) to (17) in the paper:
$$\sigma _{\rm n\comma n}^2 =\left\{\matrix{{\sigma}^2\comma \; \hfill & {\rm n}=1 \cr \sum\limits_{i=0}^{{\rm n}-2} {{\bf HA}^i} {\bf BQB}^{\rm T}\lpar {\bf A}^i\rpar^{\rm T}{\bf H}^{\rm T}+{\sigma}^2\comma \; & {\rm n}>1}\right. $$
$$\sigma_{\rm m\comma n}^2 =0\comma \; \quad {\rm m}=1\hbox{ or }{\rm n}=1 $$
$$\sigma _{\rm m\comma n}^2 = \left\{\matrix{{\bf H}\left(\sum\limits_{\mathop {i = 1}\limits_{m \ge 2} }^{m-2} {{\bf A}^{{\rm n}-{\rm m}+i}{\bf B}} {\bf QB}^{\rm T}\lpar {\bf A}^{\rm T}\rpar^i\right){\bf H}^{\rm T}\comma \; & {\rm m}< {\rm n} \cr {\bf H}\left(\sum\limits_{\mathop {i = 1}\limits_{n \ge 2} }^{n-2} {{\bf A}^{{\rm m}-{\rm n}+i}{\bf B}} {\bf QB}^{\rm T}\lpar {\bf A}^{\rm T}\rpar^i\right){\bf H}^{\rm T}\comma \; & {\rm m}>{\rm n}}\right. $$
