# Research on Self-Monitoring Method for Anomalies of Satellite Atomic Clock

Authors: Lei Feng; Guotong Li
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

Atomic clock is the core component of navigation satellite payload, playing a decisive role in the realization of positioning function. So the monitoring for anomalies of the satellite atomic clock is very important. In this paper, a complete autonomous monitoring method for the satellite clock is put forward, which is, respectively, based on Phase-Locked Loop (PLL) and statistical principle. Our methods focus on anomalies in satellite clock such as phase and frequency jumping, instantaneous deterioration, stability deterioration, and frequency drift-rate anomaly. Now, method based on PLL has been used successfully in China’s newest BeiDou navigation satellite.


---

## Body

## 1. Introduction

The most important function of navigation satellite is to support its users to acquire their position through the satellite signal, during which satellite time is one of the most important factors. Because of the changes in temperature, humidity, radiation, and the aging of the satellite clock, the physical and electric part of clock may both have problems, which will bring anomaly in clock signal, resulting in large error in the prediction of satellite time or even unpredictability, which may lead to disastrous consequence. So the anomaly monitoring of satellite clock is very important.So far, researchers have proposed schemes to monitor anomalies of clock, such as Interferometric Detection Method [1], Least Square (LS) Detection Method [2, 3], Generalized Likelihood Ratio Test (GLRT) [4–6], Kalman Filtering Method [7–10], and Dynamic Allan Variance (DAVAR) Method [5, 11–13]. Though their schemes have been proved to be effective for some (not all) anomalies, some extra work is still needed to realize Self-Monitoring.Under normal circumstance, ground station can evaluate the health condition and performance of the clock by continuously tracking satellite signals. But when the satellite flies beyond the ground station’s sight, or owing to some reasons, the satellite cannot contact with the ground station in a few hours or even days; the satellite needs to judge the status of the clock all by itself. Self-Monitoring for clock anomaly, which is in the absence of ground station, is that the satellite monitors its clock by itself to make a judgment on the satellite clock running state.The common anomalies of the satellite clock are signal loss, phase jumping, frequency jumping, instantaneous deterioration, stability, and frequency drift-rate deterioration.The contributions of this paper can be summarized as follows.In this paper, a set of Self-Monitoring algorithms is proposed to improve the reliability of satellite. Two methods are put forward to monitor satellite clock anomalies. The first method is based on PLL, and it can detect signal loss and phase and frequency jumping. Based on the measurement data from intercomparison among three clocks, Modified DAVAR is used to detect phase and frequency jumping and instantaneous deterioration; we use windowed overlapping Hadamard variance to evaluate clock stability in real time and the three-state Kalman filter to detect large drift rate.The method based on PLL has been proved effective and used in newest BeiDou satellite. And the other research on Self-Monitoring Method in this paper can be used in next generation navigation satellites.


## 2. Self-Monitoring Method for Anomaly of Satellite

Generally speaking, there are two methods to evaluate atomic clocks: (1) comparing the clock signal with standard reference whose stability is much better than the evaluated clock and (2) making intercomparison among three or more clocks whose stability is almost the same.Because there is no standard reference in the satellite and the performance of satellite clocks is similar, we make use of the second method to realize Self-Monitoring for anomalies. The schematic diagram is shown in Figure1.Figure 1
Global schematic diagram of Self-Monitoring Method.Firstly, we define thatΔ
t
1 is the time error of clock 1, which is the difference between clock time and the standard time. Δ
t
12 is the time difference between clock 1 and clock 2. As shown in Figure 1, three clocks are all powered up. Their 10 MHz signals act as the input of phase difference measurement module, through which we can get the time difference data Δ
t
12, Δ
t
13, and Δ
t
23 among them. The Signal Processing Module uses the time difference data Δ
t
12, Δ
t
13, and Δ
t
23 to evaluate the health state of three clocks with certain algorithm and then commands the master clock selector to choose suitable clock as the frequency and time source of the entire satellite.In this paper,Δ
t
12, Δ
t
13, and Δ
t
23 are used in Modified DAVAR to monitor phase and frequency jumping, used to evaluate the stability of three clocks, and used to monitor drift-rate anomaly.The detailed structure of Self-Monitoring Module based on PLL in Figure1 can be described in Figure 2.Figure 2
Schematic diagram of Self-Monitoring Method based on PLL.


### 2.1. Self-Monitoring Method Based on PLL

Figure2 shows the basic schematic diagram of this method. As a phase tracking system, PLL is used to adjust the phase of local signal to trace the reference signal. Voltage Controlled Oscillator (VCO) provides sampling clock and working clock for AD and FPGA, respectively. As the input of the Decision Module, the observed quantity of this method comes from the output of Phase Detector. Output of Decision Module will be sent to Signal Processing Module in Figure 1 to help choose the master clock. At the same time, 10 MHz signal from Figure 1 is sampled by AD in Figure 2. The working frequency of Phase Detector is 1000 Hz.Once phase or frequency jumping occurs, the output of Phase Detector in PLL will follow. In this section, the response of Phase Detector to these two anomalies will be derived.According to [14], assuming that the phase of reference signal of PLL is 2
π
f
r
t
+
φ
1
(
t
) and that of Direct Digital Synthesizer (DDS) output is 2
π
f
r
t
+
φ
2
(
t
), then we get(1)Kφ1s-φ2sFs1s=φ2s,where F
(
s
)
=
(
1
+
s
τ
2
)
/
s
τ
1 is the transfer function of the two-order ideal loop filter, 1
/
s is the normalized transfer function of the DDS, and K is the loop gain. From expression (1), we get(2)φes=ss+KFsφ1s,where φ
e
(
s
)
=
φ
1
(
s
)
-
φ
2
(
s
) is the phase difference between the reference signal and the local signal, and the error transfer function of the loop can be expressed as(3)Hes=ss+KFs=s2s2+2ξωns+ωn2,where ω
n
=
K
/
τ
1 is the undamped oscillation frequency and ξ
=
τ
2
/
2
K
/
τ
1 is the damped coefficient.In the following, the tracking property of Phase Detector for phase and frequency jumping will be deduced.


#### 2.1.1. Phase Jumping

Supposing that the phase jumping can be written asφ
1
(
t
)
=
Δ
φ
·
ε
(
t
), whose Laplace transform can be expressed as φ
1
(
s
)
=
Δ
φ
/
s, then the error response is(4)φes=s2s2+2ξωns+ωn2·Δφs=ss2+2ξωns+ωn2·Δφ.Through factorization, (4) is equivalent to(5)φes=As-s1+Bs-s2,where(6)s1=-ξωn+iωn1-ξ2,s2=-ξωn-iωn1-ξ2,(7)A=-Δφξ-i1-ξ22i1-ξ2,B=Δφξ+i1-ξ22i1-ξ2.Considering (6) and (7), the inverse Laplace transform of (5) can be expressed as(8)
φ
e
t
=
Δ
φ
e
-
ξ
ω
n
t
cos
⁡
ω
n
t
1
-
ξ
2
-
ξ
1
-
ξ
2
sin
⁡
ω
n
t
1
-
ξ
2
.From (8) we notice that the phase difference at t
=
0 reaches its peak value which has nothing to do with the loop parameters. Figure 3 is the simulation result, the PLL is locked at the beginning, and phase of reference signal jumps at t
=
500s which leads to obvious jumping in the output of Phase Detector. In Figure 3, the loop parameters of three PLLs are ω
n
=
1,4
,
10 and ξ
=
2
/
2 and the amplitude of phase jumping is equal to 1
/
1
0
8 period of reference signal. With different loop parameters, the relocking process is different. The narrower the loop bandwidth is, the slower the tracking will be.Figure 3
Response of Phase Detector to phase jumping during simulation.


#### 2.1.2. Frequency Jumping

Assuming that frequency jumping isφ
2
(
t
)
=
Δ
φ
t
·
ε
(
t
), whose Laplace transform is φ
2
(
s
)
=
Δ
φ
/
s
2, then the error response can be expressed as(9)φes=s2s2+2ξωns+ωn2·Δφs2=Δφs2+2ξωns+ωn2.After factorization,(10)A=Δω2iωn1-ξ2,B=-Δω2iωn1-ξ2.s
1 and s
2 can also be described by expression (6). According to (6) and (10), the inverse Laplace transform of (5) can be expressed as(11)φet=Δωωn1-ξ2e-ξωntsin⁡ωnt1-ξ2.As can be seen from (11), the maximum amplitude of the phase difference in tracking is inversely proportional to ω
n. Figure 4 is the simulation result, PLL is locked at the beginning, and frequency of reference signal jumps at t
=
500 s, which leads to obvious jumping in the output of Phase Detector. In Figure 4, the loop parameters are the same as in Figure 3 and frequency jumping is equal to 4
/
10
11 reference signal frequency. With different loop parameters, the relocking process is different. The narrower the loop bandwidth is, the slower the tracking will be; however the jumping is much more obvious.Figure 4
Response of Phase Detector to frequency jumping during simulation.


#### 2.1.3. Signal Loss

As shown in Figure5, assuming that the PLL has been locked and reference signal lost at t
=
150 s, the jumping amplitude of output of Phase Detector is far larger than threshold in the following section; then it turns to 0 immediately, which is easy to be detected.Figure 5
Response of Phase Detector to signal loss during simulation.


#### 2.1.4. Summary

It can be seen that, from Figures3, 4, and 5, phase jumping, frequency jumping, and signal loss will all lead to obvious jumping in Phase Detector output, which provides us with chance to monitor anomalies of satellite clock signal.


#### 2.1.5. Simulations and Detection Performance

In practice, Probability of False Alarm (PFA) and Detection Probability (PD) are usually used to evaluate the detection method. The basic principle in setting parameters (loop parameters and detection threshold) is to improve PD and minimize PFA at the same time.The loop parameters and detection threshold are mainly determined by clock noise level and required resolution. We usually use Allan variance (12) to calculate stability to evaluate the size of noise. And resolution is the minimum range of phase and frequency jumping that algorithm can distinguish.During the following simulations, we simulate 10000 realizations.Simulation 1. During the first simulation with MATLAB, ξ
=
2
/
2, ω
n
=
4, and we use two-order ideal loop filter. Assuming that the relative frequency deviation y
i (12) of clock signal follows Gauss distribution, whose Allan deviation can be expressed as 3
E
-
12
/
τ, the detection performance of the method for phase and frequency jumping is shown in Tables 1 and 2. They separately give the PD, PFA, and detection delay in 1
/
1
0
8 period phase jumping and 4
/
10
11 frequency jumping. Detection delay is defined as Δ
t
=
m
·
T, where T is the sampling interval and m is the number of sampling points that lasted from the moment anomaly occurred to the moment they are detected by the algorithm. So it is in fact determined by the output frequency of VCO in Figure 2.Table 1
Detection performance for phase jumping.
Loop parameter
Threshold
PD
PFA
Delay (m)
ξ
=
2
/
2
ω
n
=
4
8
E
-
9
1.0000
6.131
e
-
4
1
9
E
-
9
1.0000
1.193
e
-
4
1
10
E
-
9
0.9999
2.224
e
-
5
1
11
E
-
9
0.9978
2.412
e
-
6
1
12
E
-
9
0.9833
8.300
e
-
8
1Table 2
Detection performance for frequency jumping.
Loop parameter
Threshold
PD
PFA
Delay (m)
ξ
=
2
/
2
ω
n
=
4
8
E
-
9
1.0000
6.131
e
-
4
234
9
E
-
9
1.0000
1.193
e
-
4
296
10
E
-
9
1.0000
2.224
e
-
5
337
11
E
-
9
0.9933
2.412
e
-
6
380
12
E
-
9
0.9873
8.300
e
-
8
429Simulation 2. During the second simulation, detection threshold and resolution change with the stability of atomic clock. Tables 3 and 4 show it.Table 3
Detection performance of phase jumping for different clock stability.
Clock stability
Resolution
Threshold
PD
PFA
Delay (m)
3
E
-
12
/
τ
1
E
-
8
11
E
-
9
0.9978
2.412
e
-
6
1
3
E
-
11
/
τ
1
E
-
7
11
E
-
8
0.9972
6.290
e
-
7
1
3
E
-
10
/
τ
1
E
-
6
11
E
-
7
0.9955
4.370
e
-
7
1Table 4
Detection performance of frequency jumping for different clock stability.
Clock stability
Resolution
Threshold
PD
PFA
Delay (m)
3
E
-
12
/
τ
4
E
-
11
11
E
-
9
0.9933
2.412
e
-
6
380
3
E
-
11
/
τ
4
E
-
10
11
E
-
8
0.9947
6.290
e
-
7
384
3
E
-
10
/
τ
4
E
-
9
11
E
-
7
0.9978
4.370
e
-
7
381Analysis. When 0
<
ξ
<
1, PLL is called underdamped system, in which phase and frequency jumping will result in drastic oscillation. If ξ
>
1, PLL is overdamped and usually more stable and slow to anomaly. In practice, we often set ξ
=
0.707, which is an acceptable compromise between stability and response speed. From expression (11) and Figure 4, we notice that detection for frequency jumping will become difficult when ω
n is too large, and suppression for noise will also become weaker. Conversely, if ω
n is too small, on one hand, the locking process will become difficult, and detection delay becomes longer; on the other hand, the loop will be too sensitive, which leads the Decision Module to regard bottom noise as jumping by mistake frequently, which results in rising in PFA. During simulation, ω
n
=
4, which is also a compromise between PD and detection delay and can be adjusted as required.The noise level of atomic clock directly determines the detection resolution, which we can see from Tables3 and 4. The relationship between resolution and stability can be described as Re
⁡
s
(
p
)
≈
(
1
E
4
/
3
)
·
σ
(
1
) and Re
⁡
s
(
f
)
≈
(
4
E
1
/
3
)
·
σ
(
1
), while threshold can be set as T
h
r
≈
(
11
E
3
/
3
)
·
σ
(
1
). From Tables 1 and 2, we can see that PD will be more than 99
% and PFA less than 0.001
% with appropriate threshold for both phase and frequency jumping. Moreover, it should be noted that PD and PFA listed in the tables are for least phase and frequency jumping that the method can distinguish; the detection performance improves with jumping size increasing. In fact, before the method was used in the satellite, we have tested it in real circuit board for a long time and it works well. Detection delay depends on Phase Detecting frequency, which is 1000 Hz. Delay for phase jumping is 1 ms, and it is less than 0.5 s for frequency jumping.The method based on PLL can realize Self-Monitoring for phase jumping, frequency jumping, and signal loss. The computation complexity is low and costs little time to detect anomaly. But if we want to enhance its weak anomaly detection performance, we need to lower the working frequency of Phase Detector, which will lead to longer detection delay. In practice, we pay more attention to large frequency jumping in satellite atomic clock, which will obviously affect the positioning accuracy and our PLL method is designed for it.


### 2.2. Self-Monitoring Method Based on Statistics

#### 2.2.1. Allan Variance

We usually use Allan variance [15, 16] to evaluate the stability of atomic clock; it can be expressed as follows:(12)σy2τ=12M-1∑i=1M-1y¯i+1m-y¯im2,yi=xi+1-xiτ=fo-frfr,y¯im=1m∑j=1myj,where τ
=
m
τ
0 is the averaging time and M is the amount of y
¯
i
(
m
). What needs to be pointed out is that y
i is the relative frequency deviation. f
o is the instantaneous frequency, f
r is the nominal frequency, and x
i is the clock time error at the ith measuring instant.Power-law spectrum is used to analyze noise property in frequency domain:(13)Syf=h-2f-2+h-1f-1+h0+h1f+h2f2=12π2∑α=-2α=2hαfα,where S
y
(
f
) is spectrum density for relative frequency deviation and h
α is the amplitude corresponding to different noise type. The power-law spectrum model contains five kinds of noise (α
=
-
2
,
-
1,0
,
1,2); they are RW FM, Flicker FM, White FM, Flicker PM, and White PM. σ
y
2
(
τ
) is determined by these five kinds of noise:(14)σy2τ=3fhh22πτ2+1.038+3ln⁡2πfhτh12πτ2+h02τ+2ln⁡2h-1+2π2τh-23.The slope of Allan variance gives us a knowledge of noise distribution in different averaging time.After obtaining a sufficient number of measurement data, Allan variance can be used to calculate the stability of different averaging time. But when anomaly occurs, the results given by Allan variance may lose practical significance. As shown in Figure6(a), the frequency of clock signal jumped and then returned some time later. We cannot get correct judgment for noise distribution according to the computed result by (12) shown in Figure 6(b). Besides, we do not know the anomaly type and detection delay is also too long.Figure 6
Relative frequency deviation and its Allan deviation when anomaly occurs.
(a)
Measured value of relative frequency deviation
(b)
Allan deviation


#### 2.2.2. Dynamic Allan Variance (DAVAR)

As can be seen from Figure6(a), the main noise type does not change, but Figure 6(b) gives wrong judgment. Hence the conclusion is not consistent with the actual situation, and we cannot find the sign of frequency jumping either. Therefore, the traditional Allan variance cannot give believable information on such anomaly. In view of this, Galleani and Tavella put forward DAVAR, which can be expressed as (15) and can be used to evaluate the performance of clock in real time:(15)σy2n,k=12k2τ021N/k-1∑i=0N/k-2y¯i+1,k,n-y¯i,k,n2,y¯i,k,n=1k∑m=n-N+1+ikm=n-N+1+i+1kym.When we calculate DAVAR, a sliding window is used to cut the data. The window length isN, and σ
y
2
(
n
,
k
) will be updated when new measurement data y
n comes, so it can tell us health condition of clock in real time. τ
0 is the least measurement period, and k
τ
0 is the averaging time.


#### 2.2.3. Modified DAVAR

From expression (15), we know that DAVAR can be updated in real time, but in order to guarantee the reliability of long-term stability, N must be large enough, which will greatly reduce the detection probability of instantaneous anomaly. Because when frequency jumping occurs, …
,
y
i
-
1
-
y
i
-
2
≈
0
,
y
i
-
y
i
-
1
=
δ
,
y
i
+
1
-
y
i
≈
0
,
…, only one factor is not 0, DAVAR is not sensitive enough to weak frequency jumping. In this paper, we modify Dynamic Allan Variance to improve detection sensitivity for small frequency jumping:(16)σy2n=12m∑i=1my¯i+n-m-y¯i+n-2m2,(17)σy2n=σy2n-1+y¯n-y¯n-m2-y¯n-m-y¯n-2m22m.Expression (16) is the Modified DAVAR, and (17) is its iterative calculation method.


#### 2.2.4. Detection Performance of Modified DAVAR

In this section, we will firstly analyze and compare the detection performances of DAVAR and Modified DAVAR when facing phase and frequency jumping and then show that Modified DAVAR is also effective in detecting instantaneous stability deterioration.The monitoring method for phase and frequency jumping is based on statistics.Δ
t
12, Δ
t
13, and Δ
t
23 from Figure 1 will be used here. Assume that only one of the three clocks breaks down. If phase or frequency jumping occurs in clock 1, Δ
t
12, Δ
t
13 will be abnormal, while Δ
t
23 is still normal. Because the anomaly in Δ
t
12 is the same as that in Δ
t
13, we only need to analyze Δ
t
12.In Figure7(a), the amplitude of phase jumping is 12 times the standard deviation of the relative frequency deviation data y
i. In Figure 7(b), the frequency jumping is 4 times the standard deviation of y
i.Figure 7
Phase and frequency jumping.
(a)
Phase jumping
(b)
Frequency jumpingWe simulated 1000 sampling points and phase and frequency jumping occurred at the 500th point. We simulated one realization and saved the response data of DAVAR and Modified DAVAR to jumping at every time instant; then we repeated 10000 realizations in the same way. Of course the 1000 sampling points’ data is different in every realization. Then we got the average response at every time instant that is shown in Figures8 and 9.Figure 8
Comparison of two variances for phase jumping.
(a)
Average response of DAVAR to phase jumping
(b)
Average response of Modified DAVAR to phase jumpingFigure 9
Comparison of two variances for frequency jumping.
(a)
Average response of DAVAR to frequency jumping
(b)
Average response of Modified DAVAR to frequency jumpingBecause the peak value of response of DAVAR and Modified DAVAR to jumping determines if the jumping could be detected, we focus on the peak value in each realization. Assuming that, in one realization, the maximum value of response of DAVAR to frequency jumping isσ
i
,
m
2, the maximum value of response of Modified DAVAR to frequency jumping is υ
i
,
m
2, and then we saved σ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
). We studied the saved data to give the minimum value, maximum value, mean value, and standard deviation of σ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
) in Figures 8 and 9.To make the statistical result more clear, we list the statistical characteristic ofσ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
) in frequency jumping in Table 5. From Table 5 and Figure 9 we know that σ
i
,
m
2 is not big enough to be distinguished from the base noise, and Modified DAVAR is more sensitive to weak frequency jumping.Table 5
Statistical characteristic ofσ
i
,
m
2 and υ
i
,
m
2.
Statistical item
DAVAR
Modified DAVAR
Minimum value
1.60
e
-
23
3.09
e
-
23
Maximum value
5.59
e
-
23
1.49
e
-
22
Mean value
2.89
e
-
23
8.14
e
-
23
Standard deviation
5.01
e
-
24
1.70
e
-
23In Tables6 and 7, the detection performance of two kinds of variance for phase and frequency jumping is given. Tables 8 and 9 tell us the detection performance for different clock stability.Table 6
Detection performance for phase jumping.
Detection item
DAVAR
Modified DAVAR
Resolution
12
12
Threshold
4.5
E
-
23
4.5
E
-
23
PD
1.0000
0.9978
PFA
1.61
e
-
5
3.00
e
-
7
Detection delay (points)
1
1Table 7
Detection performance for phase jumping.
Detection item
DAVAR
Modified DAVAR
Resolution
12
4
Threshold
4.5
E
-
23
4.5
E
-
23
PD
0.9874
0.9927
PFA
1.61
e
-
5
3.00
e
-
7
Detection delay (points)
1
5Table 8
Detection performance of phase jumping for different clock stability.
σ
(
τ
0
)
Resolution
Threshold
PD
PFA
Delay (points)
3
E
-
12
12
4.5
E
-
23
0.9978
3.00
e
-
7
1
3
E
-
11
12
4.5
E
-
21
0.9961
5.00
e
-
7
1
3
E
-
10
12
4.5
E
-
19
0.9977
6.00
e
-
7
1Table 9
Detection performance of frequency jumping for different clock stability.
σ
(
τ
0
)
Resolution
Threshold
PD
PFA
Delay (points)
3
E
-
12
4
4.5
E
-
23
0.9927
3.00
e
-
7
5
3
E
-
11
4
4.5
E
-
21
0.9923
5.00
e
-
7
5
3
E
-
10
4
4.5
E
-
19
0.9934
6.00
e
-
7
5What should be pointed out is that the unit of resolution in phase and frequency jumping in Tables6–9 is the standard deviation of the relative frequency deviation data y
i, namely, σ
y. Detection delay is the number of sampling points from anomaly occurring to be detected.During the simulation, we simulate 10000 realizations and use same threshold for both DAVAR and Modified DAVAR. The window lengthN
=
10, and k
τ
0
=
τ
0
=
1 s. It should be noted that because k
τ
0
=
τ
0
=
1 s, it is precise enough for us to do the approximation to only consider WFM in the simulation. In fact, the Modified DAVAR is still effective in detecting phase and frequency jumping in the presence of other noise types.Table6 tells us that PD of Modified DAVAR is almost as good as DAVAR and PFA is lower. From Table 7, we know that Modified DAVAR is more sensitive to weak frequency jumping but detection delay is longer.From Tables8 and 9, we notice that the resolution of Modified DAVAR for phase and frequency jumping is the same for different clocks. What we need to do is only to reset the threshold according to expression (18):(18)Thr=5·σ21.In addition, the window lengthN is an important parameter for Modified DAVAR. The longer the window, the weaker the detection performance. However, if the window length is too short, PFA will rise and detection resolution will also deteriorate.Figure10 shows that the Modified DAVAR can also monitor instantaneous deterioration effectively.Figure 10
Instantaneous deterioration of clock.
(a)
Instantaneous deterioration
(b)
Modified DAVARModified DAVAR can be considered as a statistical tool; it is effective to detect phase and frequency jumping. Compared with PLL method, we need to measure the time error data firstly and then calculate the statistical characteristics of clock. Modified DAVAR can monitor weaker frequency jumping compared to PLL method, but PLL method is independent of a second standard reference and time-comparison device, which will give us more flexibility. Taking their respective characteristics into account, cooperation between them may be a good choice to improve the Self-Monitoring reliability.


#### 2.2.5. LS Method and Kalman Method on Detection for Frequency Jumping

In this section, we will introduce two existing methods, which are called LS method and Kalman filter method.LS Method. We may make use of LS algorithm to calculate the averaging frequency deviation y
¯ with M newest saved sampling points and then predict the time error Δ
t
^
k
+
1
=
Δ
t
k
+
y
¯
·
τ. After we can compare it with the real measurement Δ
t
k
+
1, if the difference ε
=
Δ
t
^
k
+
1
-
Δ
t
k
+
1 is beyond the threshold γ, we think that frequency jumping occurs.Kalman Filter Method. We can make use of Kalman filter to predict the next state y
^
n of clock and then compare it with the real measurement y
n. If the difference ε
=
y
^
n
-
y
n is larger than the configurable threshold, we think that anomaly occurs.Simulations. Because the sampling interval τ
=
1 s, we only consider WFM noise, whose standard deviation is σ
0
=
3
E
-
12. During the simulations for LS method, we choose M
=
20, while, for Kalman filter method, the state transition matrix Φ
=
1, observation matrix H
=
1, system error covariance matrix Q
=
(
1.0
×
10
-
13
)
2, and observation error covariance matrix R
=
(
3.0
×
10
-
12
)
2.Figures11 and 12 show the frequency jumping detection of the two methods.Figure 11
LS method in detecting frequency jumping.
(a)
ε
=
Δ
t
^
k
+
1
-
Δ
t
k
+
1
(b)
y
- and yFigure 12
Kalman filter method in detecting frequency jumping.
(a)
ε
=
y
^
n
-
y
n
(b)
y
^
n and y
nAfter we have done numerical simulations, we give Table10 to show the detection performance of the methods, in which PLL method, DAVAR method, and Modified DAVAR method are included.Table 10
Detection performance of frequency jumping for different clock stability.
Detection method
Resolution
PD
PFA
Delay
PLL method
4.0
E
-
11
0.9933
2.41
E
-
6
0.4 s
DAVAR method
3.6
E
-
11
0.9874
1.61
E
-
5
1 s
Modified DAVAR method
1.2
E
-
11
0.9927
3.00
E
-
7
5 s
LS method
1.5
E
-
11
0.9944
2.53
E
-
5
1 s
Kalman filter method
1.8
E
-
11
0.9926
2.60
E
-
6
1 sDiscussion. It should be noted firstly that the detection performance will be different with different parameters. We use the same noise level to test different methods to give Table 10, which can be a reference to show different characteristics of different methods.Different observation quantity is needed for different methods. In our opinion, DAVAR, Modified DAVAR, LS, and Kalman filter are all effective in detecting weak frequency jumping. They need a second standard reference and time-comparison device to get the time error measurementΔ
t
k to be used as observation quantity to run the algorithm. PLL method can realize Self-Monitoring for frequency jumping without standard source, which can give us more flexibility and has been used on BeiDou satellite. The computation complexity is also different among them. LS method, Kalman filter method, and PLL method should be of less computation quantity. The resolution of Modified DAVAR is the best but at the cost of longer delay, while the PLL method has the shortest delay but at the cost of worst resolution. Sometimes we may have to compromise between resolution and delay.


#### 2.2.6. Stability Evaluation of Satellite Clock

In this section, we use the well-known “three-cornered hat” approach [17–19] to evaluate the stability of atomic clocks.According to the references, we can getσ
1
2
(
τ
), σ
2
2
(
τ
), and σ
3
2
(
τ
) from Δ
t
12, Δ
t
13, and Δ
t
23 in Figure 1. For example, σ
1
2
(
τ
) can be expressed as(19)σ12τ=12σ122τ+σ132τ-σ232τ.Because Allan variance is convergent for the five kinds of noises at different averaging time, it is often used to evaluate the stability of clock. However Allan variance cannot rule out frequency drift. Especially when the drift is almost equal to Allan variance for certain averaging time, if we use Allan variance to calculateσ
12
2
(
τ
), σ
13
2
(
τ
), and σ
23
2
(
τ
) and then calculate σ
1
2
(
τ
), σ
2
2
(
τ
), and σ
3
2
(
τ
), the result cannot reflect the real situation of atomic clock. In order to avoid the influence of frequency drift, two-order difference for frequency data or three-order difference for phase data is needed, which is just the definition of Hadamard variance.In order to make full use of the measurement data and also track slow change of satellite clock timely, we use windowed overlapping Hadamard variance, as shown in (20); even though the additional overlapping differences are not all statistically independent, they nevertheless increase the number of degrees of freedom and thus improve the confidence in the estimation. Moreover, by using the latest data, the variance can evaluate the health condition of atomic clock in real time:(20)Hσy2n=16N-3kτ2·∑i=n-N+1n-3kxi+3k-3xi+2k+3xi+k-xi2,where N is the length of the window, namely, the amount of data used for each update, and b
τ
0 is the time interval between x
i and x
i
+
1, which is defined as sampling period. τ
0 is the measurement period, averaging time τ
=
m
τ
0, and k
=
m
/
b. We know that, from the characteristics of Hadamard variance, the longer the averaging time, the larger the amount of data needed.Expression (21) can be used as recursive algorithm to reduce computation complexity in the updating of H
σ
y
2
(
n
). Moreover, 6
(
N
-
3
m
)
τ
2 is a constant for each averaging time; division operation can be done only when necessary:(21)Hσy2n=Hσy2n-1+Δ1n-Δ2n6N-3mτ2,where Δ
1
n
=
x
n
-
3
x
n
-
m
+
3
x
n
-
2
m
-
x
n
-
3
m
2 and Δ
2
n
=
x
n
-
N
+
3
m
-
3
x
n
-
N
+
2
m
+
3
x
n
-
N
+
m
-
x
n
-
N
2.As opposed to the Allan variance, which makes use of a second difference, the Hadamard variance employs a third difference that leads to reduction in the degrees of freedom by one. The Hadamard variance requires more data to produce a single stability calculation, as compared to the Allan variance, given equal averaging timeτ. So it will be a better choice to use different statistical tool for different averaging time. When the averaging time is short, it is much more convenient to use an ADEV three-cornered hat method, while HDEV will be a better choice when the linear frequency drift is dominant.


#### 2.2.7. Detection for Frequency Drift-Rate Anomaly of Clock

According to [20–24], several drift-rate estimators are discussed and compared. We will firstly compare six different estimators by simulations in the following:(a)
Two points:z
^
1
=
(
y
(
n
)
-
y
(
1
)
)
/
n
-
1
τ
0
.
(b)
Two groups of points:z
^
2
=
(
2
/
n
τ
)
[
(
2
/
n
)
∑
i
=
n
/
2
+
1
n
y
(
i
)
-
(
2
/
n
)
∑
i
=
1
n
/
2
y
(
i
)
]
.
(c)
LS:z
^
3
=
(
6
/
n
n
2
-
1
τ
0
)
∑
i
=
1
n
(
2
i
-
n
-
1
)
y
(
i
)
.
(d)
Three points:z
^
4
=
(
x
2
n
+
1
-
2
x
1
+
n
+
x
1
)
/
(
n
τ
0
)
2
.
(e)
ω
4: z
^
5
=
(
6
/
N
3
τ
0
2
r
1
1
-
r
1
)
w
N
-
w
0
-
(
w
N
-
n
1
-
w
n
1
)
/
(
1
-
2
r
1
)
.
(f)
Kalman filter:X
k
=
Φ
X
k
-
1
+
ε
k, Z
k
=
H
X
k
+
ν
k,where X
k
=
x
k
,
y
k
,
z
k is the phase data, frequency deviation and drift rate, and observation matrix H
=
1,0
,
0. The state transition matrix is (22)ϕ=1ττ2201τ001and state error covariance matrix Q can be expressed as [25].To compare these six estimators, we generate simulation data by well-known Stable 32 software. The parameters of phase data are shown in Table11, and we consider three kinds of noise type which are WFM, FFM, and RWFM.Table 11
Parameters of generated phase data.
σ
0
σ
-
1
σ
-
2
Drift rate/day
Frequency deviation
1
E
-
12
1
E
-
14
1
E
-
16
1
E
-
12
1
E
-
12From Table12 and Figures 13 and 14, we know that ω
4 method and Kalman filter should be good choice for drift-rate estimation. To make sure that the simulation result is reliable, we make use of different simulation data to test the performance of the six estimators, and the result is just similar.Table 12
Statistical property of the estimation result of six estimators.
Estimator
a
b
c
d
e
f
E
(
z
)
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
9.97
E
-
13
σ
(
z
)
3.48
E
-
13
2.40
E
-
13
2.51
E
-
13
1.74
E
-
13
1.58
E
-
13
3.71
E
-
14Figure 13
Comparison of different drift-rate estimators.
(a)
Drift-rate estimation of different estimators
(b)
Partial enlarged view of (a)Figure 14
The drift-rate estimation of Kalman filter.In fact, no matter which method we choose to evaluate the drift rate, we must know the time error datax
k, which is equivalent to Δ
t in this section.When the satellite is within the sight of ground station, we can calculate the drift rate with certain estimator by comparing the satellite time with the timescale on the ground. But when the satellite cannot contact with the station, the only available data is the intercomparison dataΔ
t
12, Δ
t
13, and Δ
t
23.In order to evaluate the drift rate, there may be two methods. (1) Three-state Kalman filter can be used to evaluate the drift rate directly with the observation quantity Δ
t
12, Δ
t
13, and Δ
t
23, which is similar to the Kalman timescale algorithm on the ground. (2) Two steps are needed. The first step is to predict Δ
t
1, Δ
t
2, and Δ
t
3 from Δ
t
12, Δ
t
13, and Δ
t
23. And the second step is to evaluate drift rate with Δ
t
^
1, Δ
t
^
2, and Δ
t
^
3 (Δ
t
^
1, Δ
t
^
2, and Δ
t
^
3 are the prediction values of Δ
t
1, Δ
t
2, and Δ
t
3).In the following, we will compare these two methods. The three-state Kalman filter method is called method 1, and the combination of two-state Kalman filter withω
4 is called method 2.The basic Kalman fiter equations are as follows.System equation is(23)Xk=ΦXk-1+εk.Observation equation is(24)Zk=HXk+νk.For the two-state Kalman filter,(25)Xk=x1,k,y1,k,x2,k,y2,k,x3,k,y3,kT,Φ=diag⁡ϕ,ϕ,ϕ,ϕ=1τ01,Zk=Δt12,k,Δt13,kT,H=10-10001000-10.The system error covariance matrix Q can be expressed as [26–28].For the three-state Kalman filter,(26)Xk=x1,k,y1,k,z1,k,x2,k,y2,k,z2,k,x3,k,y3,k,z3,kT,Φ=diag⁡ϕ,ϕ,ϕ,ϕ=1ττ2201τ001,Zk=Δt12,k,Δt13,kT,H=100-100000100000-100,and we can get the system error covariance matrix Q according to [29].To compare these two methods, we make use of Stable 32 to generate simulation data. The parameters of the phase data are shown in Table13.Table 13
Parameter setting.
Parameters
σ
0
σ
-
1
σ
-
2
Drift rate/day
Frequency deviation
Clock 1
1
e
-
12
2
e
-
14
3
e
-
16
1
e
-
13
1
e
-
12
Clock 2
2
e
-
12
3
e
-
14
1
e
-
16
1
e
-
12
1
e
-
12
Clock 3
3
e
-
12
1
e
-
14
2
e
-
16
-
2
e
-
11
1
e
-
12Firstly, we will show the prediction performance of both two-state and three-state Kalman filter forΔ
t
1, Δ
t
2, and Δ
t
3; they are similar, which can be shown in Figure 15.Figure 15
Prediction performance of Kalman filter forΔ
t
1 and Δ
t
12.
(a)
Δ
t
12 and Δ
t
^
12
(b)
Δ
t
1 and Δ
t
^
1For method 2, after we have got the prediction ofΔ
t
1, Δ
t
2, and Δ
t
3, we will use ω
4 to evaluate the drift rate.Figure15(a) is the comparison of the clock difference Δ
t
12
(
Δ
t
12
=
Δ
t
1
-
Δ
t
2
) and its prediction Δ
t
^
12. Figure 15(b) compares the real time errors Δ
t
1 and Δ
t
^
1. From Figure 15, we know that the prediction error ε  (
ε
=
Δ
t
1
-
Δ
t
^
1
) will increase gradually, while Δ
t
^
12 is unbiased. In fact, this phenomenon is inevitable owing to the lack of absolute standard reference. No matter which method we choose, the prediction error will increase with time going. That also indicates that we cannot get precise drift rate, but it does not mean we can do nothing. Next we will prove that method 1 can give an alarm when the drift rate of one clock is much larger than the best one. It should be noted that when we make use of ω
4 method to evaluate drift rate for method 2, a sliding window is used to evaluate the drift rate in real time instead of batch processing.Figures16 and 17 show the drift-rate estimation of these two methods for three clocks. From the simulation result, we know that if the drift rate of one clock is much larger than the best clock, its drift-rate estimation is very close to its real value, and its computed estimation is much larger than the best clock, which enables method 1 to give an alarm. Besides, method 1 gives more precise and stable computation result compared to method 2.Figure 16
Drift-rate estimation of method 1.
(a)
Drift rate of clock 1
(b)
Drift rate of clock 2
(c)
Drift rate of clock 3Figure 17
Drift-rate estimation of method 2.
(a)
Drift rate of clock 1
(b)
Drift rate of clock 2
(c)
Drift rate of clock 3


## 2.1. Self-Monitoring Method Based on PLL

Figure2 shows the basic schematic diagram of this method. As a phase tracking system, PLL is used to adjust the phase of local signal to trace the reference signal. Voltage Controlled Oscillator (VCO) provides sampling clock and working clock for AD and FPGA, respectively. As the input of the Decision Module, the observed quantity of this method comes from the output of Phase Detector. Output of Decision Module will be sent to Signal Processing Module in Figure 1 to help choose the master clock. At the same time, 10 MHz signal from Figure 1 is sampled by AD in Figure 2. The working frequency of Phase Detector is 1000 Hz.Once phase or frequency jumping occurs, the output of Phase Detector in PLL will follow. In this section, the response of Phase Detector to these two anomalies will be derived.According to [14], assuming that the phase of reference signal of PLL is 2
π
f
r
t
+
φ
1
(
t
) and that of Direct Digital Synthesizer (DDS) output is 2
π
f
r
t
+
φ
2
(
t
), then we get(1)Kφ1s-φ2sFs1s=φ2s,where F
(
s
)
=
(
1
+
s
τ
2
)
/
s
τ
1 is the transfer function of the two-order ideal loop filter, 1
/
s is the normalized transfer function of the DDS, and K is the loop gain. From expression (1), we get(2)φes=ss+KFsφ1s,where φ
e
(
s
)
=
φ
1
(
s
)
-
φ
2
(
s
) is the phase difference between the reference signal and the local signal, and the error transfer function of the loop can be expressed as(3)Hes=ss+KFs=s2s2+2ξωns+ωn2,where ω
n
=
K
/
τ
1 is the undamped oscillation frequency and ξ
=
τ
2
/
2
K
/
τ
1 is the damped coefficient.In the following, the tracking property of Phase Detector for phase and frequency jumping will be deduced.


### 2.1.1. Phase Jumping

Supposing that the phase jumping can be written asφ
1
(
t
)
=
Δ
φ
·
ε
(
t
), whose Laplace transform can be expressed as φ
1
(
s
)
=
Δ
φ
/
s, then the error response is(4)φes=s2s2+2ξωns+ωn2·Δφs=ss2+2ξωns+ωn2·Δφ.Through factorization, (4) is equivalent to(5)φes=As-s1+Bs-s2,where(6)s1=-ξωn+iωn1-ξ2,s2=-ξωn-iωn1-ξ2,(7)A=-Δφξ-i1-ξ22i1-ξ2,B=Δφξ+i1-ξ22i1-ξ2.Considering (6) and (7), the inverse Laplace transform of (5) can be expressed as(8)
φ
e
t
=
Δ
φ
e
-
ξ
ω
n
t
cos
⁡
ω
n
t
1
-
ξ
2
-
ξ
1
-
ξ
2
sin
⁡
ω
n
t
1
-
ξ
2
.From (8) we notice that the phase difference at t
=
0 reaches its peak value which has nothing to do with the loop parameters. Figure 3 is the simulation result, the PLL is locked at the beginning, and phase of reference signal jumps at t
=
500s which leads to obvious jumping in the output of Phase Detector. In Figure 3, the loop parameters of three PLLs are ω
n
=
1,4
,
10 and ξ
=
2
/
2 and the amplitude of phase jumping is equal to 1
/
1
0
8 period of reference signal. With different loop parameters, the relocking process is different. The narrower the loop bandwidth is, the slower the tracking will be.Figure 3
Response of Phase Detector to phase jumping during simulation.


### 2.1.2. Frequency Jumping

Assuming that frequency jumping isφ
2
(
t
)
=
Δ
φ
t
·
ε
(
t
), whose Laplace transform is φ
2
(
s
)
=
Δ
φ
/
s
2, then the error response can be expressed as(9)φes=s2s2+2ξωns+ωn2·Δφs2=Δφs2+2ξωns+ωn2.After factorization,(10)A=Δω2iωn1-ξ2,B=-Δω2iωn1-ξ2.s
1 and s
2 can also be described by expression (6). According to (6) and (10), the inverse Laplace transform of (5) can be expressed as(11)φet=Δωωn1-ξ2e-ξωntsin⁡ωnt1-ξ2.As can be seen from (11), the maximum amplitude of the phase difference in tracking is inversely proportional to ω
n. Figure 4 is the simulation result, PLL is locked at the beginning, and frequency of reference signal jumps at t
=
500 s, which leads to obvious jumping in the output of Phase Detector. In Figure 4, the loop parameters are the same as in Figure 3 and frequency jumping is equal to 4
/
10
11 reference signal frequency. With different loop parameters, the relocking process is different. The narrower the loop bandwidth is, the slower the tracking will be; however the jumping is much more obvious.Figure 4
Response of Phase Detector to frequency jumping during simulation.


### 2.1.3. Signal Loss

As shown in Figure5, assuming that the PLL has been locked and reference signal lost at t
=
150 s, the jumping amplitude of output of Phase Detector is far larger than threshold in the following section; then it turns to 0 immediately, which is easy to be detected.Figure 5
Response of Phase Detector to signal loss during simulation.


### 2.1.4. Summary

It can be seen that, from Figures3, 4, and 5, phase jumping, frequency jumping, and signal loss will all lead to obvious jumping in Phase Detector output, which provides us with chance to monitor anomalies of satellite clock signal.


### 2.1.5. Simulations and Detection Performance

In practice, Probability of False Alarm (PFA) and Detection Probability (PD) are usually used to evaluate the detection method. The basic principle in setting parameters (loop parameters and detection threshold) is to improve PD and minimize PFA at the same time.The loop parameters and detection threshold are mainly determined by clock noise level and required resolution. We usually use Allan variance (12) to calculate stability to evaluate the size of noise. And resolution is the minimum range of phase and frequency jumping that algorithm can distinguish.During the following simulations, we simulate 10000 realizations.Simulation 1. During the first simulation with MATLAB, ξ
=
2
/
2, ω
n
=
4, and we use two-order ideal loop filter. Assuming that the relative frequency deviation y
i (12) of clock signal follows Gauss distribution, whose Allan deviation can be expressed as 3
E
-
12
/
τ, the detection performance of the method for phase and frequency jumping is shown in Tables 1 and 2. They separately give the PD, PFA, and detection delay in 1
/
1
0
8 period phase jumping and 4
/
10
11 frequency jumping. Detection delay is defined as Δ
t
=
m
·
T, where T is the sampling interval and m is the number of sampling points that lasted from the moment anomaly occurred to the moment they are detected by the algorithm. So it is in fact determined by the output frequency of VCO in Figure 2.Table 1
Detection performance for phase jumping.
Loop parameter
Threshold
PD
PFA
Delay (m)
ξ
=
2
/
2
ω
n
=
4
8
E
-
9
1.0000
6.131
e
-
4
1
9
E
-
9
1.0000
1.193
e
-
4
1
10
E
-
9
0.9999
2.224
e
-
5
1
11
E
-
9
0.9978
2.412
e
-
6
1
12
E
-
9
0.9833
8.300
e
-
8
1Table 2
Detection performance for frequency jumping.
Loop parameter
Threshold
PD
PFA
Delay (m)
ξ
=
2
/
2
ω
n
=
4
8
E
-
9
1.0000
6.131
e
-
4
234
9
E
-
9
1.0000
1.193
e
-
4
296
10
E
-
9
1.0000
2.224
e
-
5
337
11
E
-
9
0.9933
2.412
e
-
6
380
12
E
-
9
0.9873
8.300
e
-
8
429Simulation 2. During the second simulation, detection threshold and resolution change with the stability of atomic clock. Tables 3 and 4 show it.Table 3
Detection performance of phase jumping for different clock stability.
Clock stability
Resolution
Threshold
PD
PFA
Delay (m)
3
E
-
12
/
τ
1
E
-
8
11
E
-
9
0.9978
2.412
e
-
6
1
3
E
-
11
/
τ
1
E
-
7
11
E
-
8
0.9972
6.290
e
-
7
1
3
E
-
10
/
τ
1
E
-
6
11
E
-
7
0.9955
4.370
e
-
7
1Table 4
Detection performance of frequency jumping for different clock stability.
Clock stability
Resolution
Threshold
PD
PFA
Delay (m)
3
E
-
12
/
τ
4
E
-
11
11
E
-
9
0.9933
2.412
e
-
6
380
3
E
-
11
/
τ
4
E
-
10
11
E
-
8
0.9947
6.290
e
-
7
384
3
E
-
10
/
τ
4
E
-
9
11
E
-
7
0.9978
4.370
e
-
7
381Analysis. When 0
<
ξ
<
1, PLL is called underdamped system, in which phase and frequency jumping will result in drastic oscillation. If ξ
>
1, PLL is overdamped and usually more stable and slow to anomaly. In practice, we often set ξ
=
0.707, which is an acceptable compromise between stability and response speed. From expression (11) and Figure 4, we notice that detection for frequency jumping will become difficult when ω
n is too large, and suppression for noise will also become weaker. Conversely, if ω
n is too small, on one hand, the locking process will become difficult, and detection delay becomes longer; on the other hand, the loop will be too sensitive, which leads the Decision Module to regard bottom noise as jumping by mistake frequently, which results in rising in PFA. During simulation, ω
n
=
4, which is also a compromise between PD and detection delay and can be adjusted as required.The noise level of atomic clock directly determines the detection resolution, which we can see from Tables3 and 4. The relationship between resolution and stability can be described as Re
⁡
s
(
p
)
≈
(
1
E
4
/
3
)
·
σ
(
1
) and Re
⁡
s
(
f
)
≈
(
4
E
1
/
3
)
·
σ
(
1
), while threshold can be set as T
h
r
≈
(
11
E
3
/
3
)
·
σ
(
1
). From Tables 1 and 2, we can see that PD will be more than 99
% and PFA less than 0.001
% with appropriate threshold for both phase and frequency jumping. Moreover, it should be noted that PD and PFA listed in the tables are for least phase and frequency jumping that the method can distinguish; the detection performance improves with jumping size increasing. In fact, before the method was used in the satellite, we have tested it in real circuit board for a long time and it works well. Detection delay depends on Phase Detecting frequency, which is 1000 Hz. Delay for phase jumping is 1 ms, and it is less than 0.5 s for frequency jumping.The method based on PLL can realize Self-Monitoring for phase jumping, frequency jumping, and signal loss. The computation complexity is low and costs little time to detect anomaly. But if we want to enhance its weak anomaly detection performance, we need to lower the working frequency of Phase Detector, which will lead to longer detection delay. In practice, we pay more attention to large frequency jumping in satellite atomic clock, which will obviously affect the positioning accuracy and our PLL method is designed for it.


## 2.1.1. Phase Jumping

Supposing that the phase jumping can be written asφ
1
(
t
)
=
Δ
φ
·
ε
(
t
), whose Laplace transform can be expressed as φ
1
(
s
)
=
Δ
φ
/
s, then the error response is(4)φes=s2s2+2ξωns+ωn2·Δφs=ss2+2ξωns+ωn2·Δφ.Through factorization, (4) is equivalent to(5)φes=As-s1+Bs-s2,where(6)s1=-ξωn+iωn1-ξ2,s2=-ξωn-iωn1-ξ2,(7)A=-Δφξ-i1-ξ22i1-ξ2,B=Δφξ+i1-ξ22i1-ξ2.Considering (6) and (7), the inverse Laplace transform of (5) can be expressed as(8)
φ
e
t
=
Δ
φ
e
-
ξ
ω
n
t
cos
⁡
ω
n
t
1
-
ξ
2
-
ξ
1
-
ξ
2
sin
⁡
ω
n
t
1
-
ξ
2
.From (8) we notice that the phase difference at t
=
0 reaches its peak value which has nothing to do with the loop parameters. Figure 3 is the simulation result, the PLL is locked at the beginning, and phase of reference signal jumps at t
=
500s which leads to obvious jumping in the output of Phase Detector. In Figure 3, the loop parameters of three PLLs are ω
n
=
1,4
,
10 and ξ
=
2
/
2 and the amplitude of phase jumping is equal to 1
/
1
0
8 period of reference signal. With different loop parameters, the relocking process is different. The narrower the loop bandwidth is, the slower the tracking will be.Figure 3
Response of Phase Detector to phase jumping during simulation.


## 2.1.2. Frequency Jumping

Assuming that frequency jumping isφ
2
(
t
)
=
Δ
φ
t
·
ε
(
t
), whose Laplace transform is φ
2
(
s
)
=
Δ
φ
/
s
2, then the error response can be expressed as(9)φes=s2s2+2ξωns+ωn2·Δφs2=Δφs2+2ξωns+ωn2.After factorization,(10)A=Δω2iωn1-ξ2,B=-Δω2iωn1-ξ2.s
1 and s
2 can also be described by expression (6). According to (6) and (10), the inverse Laplace transform of (5) can be expressed as(11)φet=Δωωn1-ξ2e-ξωntsin⁡ωnt1-ξ2.As can be seen from (11), the maximum amplitude of the phase difference in tracking is inversely proportional to ω
n. Figure 4 is the simulation result, PLL is locked at the beginning, and frequency of reference signal jumps at t
=
500 s, which leads to obvious jumping in the output of Phase Detector. In Figure 4, the loop parameters are the same as in Figure 3 and frequency jumping is equal to 4
/
10
11 reference signal frequency. With different loop parameters, the relocking process is different. The narrower the loop bandwidth is, the slower the tracking will be; however the jumping is much more obvious.Figure 4
Response of Phase Detector to frequency jumping during simulation.


## 2.1.3. Signal Loss

As shown in Figure5, assuming that the PLL has been locked and reference signal lost at t
=
150 s, the jumping amplitude of output of Phase Detector is far larger than threshold in the following section; then it turns to 0 immediately, which is easy to be detected.Figure 5
Response of Phase Detector to signal loss during simulation.


## 2.1.4. Summary

It can be seen that, from Figures3, 4, and 5, phase jumping, frequency jumping, and signal loss will all lead to obvious jumping in Phase Detector output, which provides us with chance to monitor anomalies of satellite clock signal.


## 2.1.5. Simulations and Detection Performance

In practice, Probability of False Alarm (PFA) and Detection Probability (PD) are usually used to evaluate the detection method. The basic principle in setting parameters (loop parameters and detection threshold) is to improve PD and minimize PFA at the same time.The loop parameters and detection threshold are mainly determined by clock noise level and required resolution. We usually use Allan variance (12) to calculate stability to evaluate the size of noise. And resolution is the minimum range of phase and frequency jumping that algorithm can distinguish.During the following simulations, we simulate 10000 realizations.Simulation 1. During the first simulation with MATLAB, ξ
=
2
/
2, ω
n
=
4, and we use two-order ideal loop filter. Assuming that the relative frequency deviation y
i (12) of clock signal follows Gauss distribution, whose Allan deviation can be expressed as 3
E
-
12
/
τ, the detection performance of the method for phase and frequency jumping is shown in Tables 1 and 2. They separately give the PD, PFA, and detection delay in 1
/
1
0
8 period phase jumping and 4
/
10
11 frequency jumping. Detection delay is defined as Δ
t
=
m
·
T, where T is the sampling interval and m is the number of sampling points that lasted from the moment anomaly occurred to the moment they are detected by the algorithm. So it is in fact determined by the output frequency of VCO in Figure 2.Table 1
Detection performance for phase jumping.
Loop parameter
Threshold
PD
PFA
Delay (m)
ξ
=
2
/
2
ω
n
=
4
8
E
-
9
1.0000
6.131
e
-
4
1
9
E
-
9
1.0000
1.193
e
-
4
1
10
E
-
9
0.9999
2.224
e
-
5
1
11
E
-
9
0.9978
2.412
e
-
6
1
12
E
-
9
0.9833
8.300
e
-
8
1Table 2
Detection performance for frequency jumping.
Loop parameter
Threshold
PD
PFA
Delay (m)
ξ
=
2
/
2
ω
n
=
4
8
E
-
9
1.0000
6.131
e
-
4
234
9
E
-
9
1.0000
1.193
e
-
4
296
10
E
-
9
1.0000
2.224
e
-
5
337
11
E
-
9
0.9933
2.412
e
-
6
380
12
E
-
9
0.9873
8.300
e
-
8
429Simulation 2. During the second simulation, detection threshold and resolution change with the stability of atomic clock. Tables 3 and 4 show it.Table 3
Detection performance of phase jumping for different clock stability.
Clock stability
Resolution
Threshold
PD
PFA
Delay (m)
3
E
-
12
/
τ
1
E
-
8
11
E
-
9
0.9978
2.412
e
-
6
1
3
E
-
11
/
τ
1
E
-
7
11
E
-
8
0.9972
6.290
e
-
7
1
3
E
-
10
/
τ
1
E
-
6
11
E
-
7
0.9955
4.370
e
-
7
1Table 4
Detection performance of frequency jumping for different clock stability.
Clock stability
Resolution
Threshold
PD
PFA
Delay (m)
3
E
-
12
/
τ
4
E
-
11
11
E
-
9
0.9933
2.412
e
-
6
380
3
E
-
11
/
τ
4
E
-
10
11
E
-
8
0.9947
6.290
e
-
7
384
3
E
-
10
/
τ
4
E
-
9
11
E
-
7
0.9978
4.370
e
-
7
381Analysis. When 0
<
ξ
<
1, PLL is called underdamped system, in which phase and frequency jumping will result in drastic oscillation. If ξ
>
1, PLL is overdamped and usually more stable and slow to anomaly. In practice, we often set ξ
=
0.707, which is an acceptable compromise between stability and response speed. From expression (11) and Figure 4, we notice that detection for frequency jumping will become difficult when ω
n is too large, and suppression for noise will also become weaker. Conversely, if ω
n is too small, on one hand, the locking process will become difficult, and detection delay becomes longer; on the other hand, the loop will be too sensitive, which leads the Decision Module to regard bottom noise as jumping by mistake frequently, which results in rising in PFA. During simulation, ω
n
=
4, which is also a compromise between PD and detection delay and can be adjusted as required.The noise level of atomic clock directly determines the detection resolution, which we can see from Tables3 and 4. The relationship between resolution and stability can be described as Re
⁡
s
(
p
)
≈
(
1
E
4
/
3
)
·
σ
(
1
) and Re
⁡
s
(
f
)
≈
(
4
E
1
/
3
)
·
σ
(
1
), while threshold can be set as T
h
r
≈
(
11
E
3
/
3
)
·
σ
(
1
). From Tables 1 and 2, we can see that PD will be more than 99
% and PFA less than 0.001
% with appropriate threshold for both phase and frequency jumping. Moreover, it should be noted that PD and PFA listed in the tables are for least phase and frequency jumping that the method can distinguish; the detection performance improves with jumping size increasing. In fact, before the method was used in the satellite, we have tested it in real circuit board for a long time and it works well. Detection delay depends on Phase Detecting frequency, which is 1000 Hz. Delay for phase jumping is 1 ms, and it is less than 0.5 s for frequency jumping.The method based on PLL can realize Self-Monitoring for phase jumping, frequency jumping, and signal loss. The computation complexity is low and costs little time to detect anomaly. But if we want to enhance its weak anomaly detection performance, we need to lower the working frequency of Phase Detector, which will lead to longer detection delay. In practice, we pay more attention to large frequency jumping in satellite atomic clock, which will obviously affect the positioning accuracy and our PLL method is designed for it.


## 2.2. Self-Monitoring Method Based on Statistics

### 2.2.1. Allan Variance

We usually use Allan variance [15, 16] to evaluate the stability of atomic clock; it can be expressed as follows:(12)σy2τ=12M-1∑i=1M-1y¯i+1m-y¯im2,yi=xi+1-xiτ=fo-frfr,y¯im=1m∑j=1myj,where τ
=
m
τ
0 is the averaging time and M is the amount of y
¯
i
(
m
). What needs to be pointed out is that y
i is the relative frequency deviation. f
o is the instantaneous frequency, f
r is the nominal frequency, and x
i is the clock time error at the ith measuring instant.Power-law spectrum is used to analyze noise property in frequency domain:(13)Syf=h-2f-2+h-1f-1+h0+h1f+h2f2=12π2∑α=-2α=2hαfα,where S
y
(
f
) is spectrum density for relative frequency deviation and h
α is the amplitude corresponding to different noise type. The power-law spectrum model contains five kinds of noise (α
=
-
2
,
-
1,0
,
1,2); they are RW FM, Flicker FM, White FM, Flicker PM, and White PM. σ
y
2
(
τ
) is determined by these five kinds of noise:(14)σy2τ=3fhh22πτ2+1.038+3ln⁡2πfhτh12πτ2+h02τ+2ln⁡2h-1+2π2τh-23.The slope of Allan variance gives us a knowledge of noise distribution in different averaging time.After obtaining a sufficient number of measurement data, Allan variance can be used to calculate the stability of different averaging time. But when anomaly occurs, the results given by Allan variance may lose practical significance. As shown in Figure6(a), the frequency of clock signal jumped and then returned some time later. We cannot get correct judgment for noise distribution according to the computed result by (12) shown in Figure 6(b). Besides, we do not know the anomaly type and detection delay is also too long.Figure 6
Relative frequency deviation and its Allan deviation when anomaly occurs.
(a)
Measured value of relative frequency deviation
(b)
Allan deviation


### 2.2.2. Dynamic Allan Variance (DAVAR)

As can be seen from Figure6(a), the main noise type does not change, but Figure 6(b) gives wrong judgment. Hence the conclusion is not consistent with the actual situation, and we cannot find the sign of frequency jumping either. Therefore, the traditional Allan variance cannot give believable information on such anomaly. In view of this, Galleani and Tavella put forward DAVAR, which can be expressed as (15) and can be used to evaluate the performance of clock in real time:(15)σy2n,k=12k2τ021N/k-1∑i=0N/k-2y¯i+1,k,n-y¯i,k,n2,y¯i,k,n=1k∑m=n-N+1+ikm=n-N+1+i+1kym.When we calculate DAVAR, a sliding window is used to cut the data. The window length isN, and σ
y
2
(
n
,
k
) will be updated when new measurement data y
n comes, so it can tell us health condition of clock in real time. τ
0 is the least measurement period, and k
τ
0 is the averaging time.


### 2.2.3. Modified DAVAR

From expression (15), we know that DAVAR can be updated in real time, but in order to guarantee the reliability of long-term stability, N must be large enough, which will greatly reduce the detection probability of instantaneous anomaly. Because when frequency jumping occurs, …
,
y
i
-
1
-
y
i
-
2
≈
0
,
y
i
-
y
i
-
1
=
δ
,
y
i
+
1
-
y
i
≈
0
,
…, only one factor is not 0, DAVAR is not sensitive enough to weak frequency jumping. In this paper, we modify Dynamic Allan Variance to improve detection sensitivity for small frequency jumping:(16)σy2n=12m∑i=1my¯i+n-m-y¯i+n-2m2,(17)σy2n=σy2n-1+y¯n-y¯n-m2-y¯n-m-y¯n-2m22m.Expression (16) is the Modified DAVAR, and (17) is its iterative calculation method.


### 2.2.4. Detection Performance of Modified DAVAR

In this section, we will firstly analyze and compare the detection performances of DAVAR and Modified DAVAR when facing phase and frequency jumping and then show that Modified DAVAR is also effective in detecting instantaneous stability deterioration.The monitoring method for phase and frequency jumping is based on statistics.Δ
t
12, Δ
t
13, and Δ
t
23 from Figure 1 will be used here. Assume that only one of the three clocks breaks down. If phase or frequency jumping occurs in clock 1, Δ
t
12, Δ
t
13 will be abnormal, while Δ
t
23 is still normal. Because the anomaly in Δ
t
12 is the same as that in Δ
t
13, we only need to analyze Δ
t
12.In Figure7(a), the amplitude of phase jumping is 12 times the standard deviation of the relative frequency deviation data y
i. In Figure 7(b), the frequency jumping is 4 times the standard deviation of y
i.Figure 7
Phase and frequency jumping.
(a)
Phase jumping
(b)
Frequency jumpingWe simulated 1000 sampling points and phase and frequency jumping occurred at the 500th point. We simulated one realization and saved the response data of DAVAR and Modified DAVAR to jumping at every time instant; then we repeated 10000 realizations in the same way. Of course the 1000 sampling points’ data is different in every realization. Then we got the average response at every time instant that is shown in Figures8 and 9.Figure 8
Comparison of two variances for phase jumping.
(a)
Average response of DAVAR to phase jumping
(b)
Average response of Modified DAVAR to phase jumpingFigure 9
Comparison of two variances for frequency jumping.
(a)
Average response of DAVAR to frequency jumping
(b)
Average response of Modified DAVAR to frequency jumpingBecause the peak value of response of DAVAR and Modified DAVAR to jumping determines if the jumping could be detected, we focus on the peak value in each realization. Assuming that, in one realization, the maximum value of response of DAVAR to frequency jumping isσ
i
,
m
2, the maximum value of response of Modified DAVAR to frequency jumping is υ
i
,
m
2, and then we saved σ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
). We studied the saved data to give the minimum value, maximum value, mean value, and standard deviation of σ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
) in Figures 8 and 9.To make the statistical result more clear, we list the statistical characteristic ofσ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
) in frequency jumping in Table 5. From Table 5 and Figure 9 we know that σ
i
,
m
2 is not big enough to be distinguished from the base noise, and Modified DAVAR is more sensitive to weak frequency jumping.Table 5
Statistical characteristic ofσ
i
,
m
2 and υ
i
,
m
2.
Statistical item
DAVAR
Modified DAVAR
Minimum value
1.60
e
-
23
3.09
e
-
23
Maximum value
5.59
e
-
23
1.49
e
-
22
Mean value
2.89
e
-
23
8.14
e
-
23
Standard deviation
5.01
e
-
24
1.70
e
-
23In Tables6 and 7, the detection performance of two kinds of variance for phase and frequency jumping is given. Tables 8 and 9 tell us the detection performance for different clock stability.Table 6
Detection performance for phase jumping.
Detection item
DAVAR
Modified DAVAR
Resolution
12
12
Threshold
4.5
E
-
23
4.5
E
-
23
PD
1.0000
0.9978
PFA
1.61
e
-
5
3.00
e
-
7
Detection delay (points)
1
1Table 7
Detection performance for phase jumping.
Detection item
DAVAR
Modified DAVAR
Resolution
12
4
Threshold
4.5
E
-
23
4.5
E
-
23
PD
0.9874
0.9927
PFA
1.61
e
-
5
3.00
e
-
7
Detection delay (points)
1
5Table 8
Detection performance of phase jumping for different clock stability.
σ
(
τ
0
)
Resolution
Threshold
PD
PFA
Delay (points)
3
E
-
12
12
4.5
E
-
23
0.9978
3.00
e
-
7
1
3
E
-
11
12
4.5
E
-
21
0.9961
5.00
e
-
7
1
3
E
-
10
12
4.5
E
-
19
0.9977
6.00
e
-
7
1Table 9
Detection performance of frequency jumping for different clock stability.
σ
(
τ
0
)
Resolution
Threshold
PD
PFA
Delay (points)
3
E
-
12
4
4.5
E
-
23
0.9927
3.00
e
-
7
5
3
E
-
11
4
4.5
E
-
21
0.9923
5.00
e
-
7
5
3
E
-
10
4
4.5
E
-
19
0.9934
6.00
e
-
7
5What should be pointed out is that the unit of resolution in phase and frequency jumping in Tables6–9 is the standard deviation of the relative frequency deviation data y
i, namely, σ
y. Detection delay is the number of sampling points from anomaly occurring to be detected.During the simulation, we simulate 10000 realizations and use same threshold for both DAVAR and Modified DAVAR. The window lengthN
=
10, and k
τ
0
=
τ
0
=
1 s. It should be noted that because k
τ
0
=
τ
0
=
1 s, it is precise enough for us to do the approximation to only consider WFM in the simulation. In fact, the Modified DAVAR is still effective in detecting phase and frequency jumping in the presence of other noise types.Table6 tells us that PD of Modified DAVAR is almost as good as DAVAR and PFA is lower. From Table 7, we know that Modified DAVAR is more sensitive to weak frequency jumping but detection delay is longer.From Tables8 and 9, we notice that the resolution of Modified DAVAR for phase and frequency jumping is the same for different clocks. What we need to do is only to reset the threshold according to expression (18):(18)Thr=5·σ21.In addition, the window lengthN is an important parameter for Modified DAVAR. The longer the window, the weaker the detection performance. However, if the window length is too short, PFA will rise and detection resolution will also deteriorate.Figure10 shows that the Modified DAVAR can also monitor instantaneous deterioration effectively.Figure 10
Instantaneous deterioration of clock.
(a)
Instantaneous deterioration
(b)
Modified DAVARModified DAVAR can be considered as a statistical tool; it is effective to detect phase and frequency jumping. Compared with PLL method, we need to measure the time error data firstly and then calculate the statistical characteristics of clock. Modified DAVAR can monitor weaker frequency jumping compared to PLL method, but PLL method is independent of a second standard reference and time-comparison device, which will give us more flexibility. Taking their respective characteristics into account, cooperation between them may be a good choice to improve the Self-Monitoring reliability.


### 2.2.5. LS Method and Kalman Method on Detection for Frequency Jumping

In this section, we will introduce two existing methods, which are called LS method and Kalman filter method.LS Method. We may make use of LS algorithm to calculate the averaging frequency deviation y
¯ with M newest saved sampling points and then predict the time error Δ
t
^
k
+
1
=
Δ
t
k
+
y
¯
·
τ. After we can compare it with the real measurement Δ
t
k
+
1, if the difference ε
=
Δ
t
^
k
+
1
-
Δ
t
k
+
1 is beyond the threshold γ, we think that frequency jumping occurs.Kalman Filter Method. We can make use of Kalman filter to predict the next state y
^
n of clock and then compare it with the real measurement y
n. If the difference ε
=
y
^
n
-
y
n is larger than the configurable threshold, we think that anomaly occurs.Simulations. Because the sampling interval τ
=
1 s, we only consider WFM noise, whose standard deviation is σ
0
=
3
E
-
12. During the simulations for LS method, we choose M
=
20, while, for Kalman filter method, the state transition matrix Φ
=
1, observation matrix H
=
1, system error covariance matrix Q
=
(
1.0
×
10
-
13
)
2, and observation error covariance matrix R
=
(
3.0
×
10
-
12
)
2.Figures11 and 12 show the frequency jumping detection of the two methods.Figure 11
LS method in detecting frequency jumping.
(a)
ε
=
Δ
t
^
k
+
1
-
Δ
t
k
+
1
(b)
y
- and yFigure 12
Kalman filter method in detecting frequency jumping.
(a)
ε
=
y
^
n
-
y
n
(b)
y
^
n and y
nAfter we have done numerical simulations, we give Table10 to show the detection performance of the methods, in which PLL method, DAVAR method, and Modified DAVAR method are included.Table 10
Detection performance of frequency jumping for different clock stability.
Detection method
Resolution
PD
PFA
Delay
PLL method
4.0
E
-
11
0.9933
2.41
E
-
6
0.4 s
DAVAR method
3.6
E
-
11
0.9874
1.61
E
-
5
1 s
Modified DAVAR method
1.2
E
-
11
0.9927
3.00
E
-
7
5 s
LS method
1.5
E
-
11
0.9944
2.53
E
-
5
1 s
Kalman filter method
1.8
E
-
11
0.9926
2.60
E
-
6
1 sDiscussion. It should be noted firstly that the detection performance will be different with different parameters. We use the same noise level to test different methods to give Table 10, which can be a reference to show different characteristics of different methods.Different observation quantity is needed for different methods. In our opinion, DAVAR, Modified DAVAR, LS, and Kalman filter are all effective in detecting weak frequency jumping. They need a second standard reference and time-comparison device to get the time error measurementΔ
t
k to be used as observation quantity to run the algorithm. PLL method can realize Self-Monitoring for frequency jumping without standard source, which can give us more flexibility and has been used on BeiDou satellite. The computation complexity is also different among them. LS method, Kalman filter method, and PLL method should be of less computation quantity. The resolution of Modified DAVAR is the best but at the cost of longer delay, while the PLL method has the shortest delay but at the cost of worst resolution. Sometimes we may have to compromise between resolution and delay.


### 2.2.6. Stability Evaluation of Satellite Clock

In this section, we use the well-known “three-cornered hat” approach [17–19] to evaluate the stability of atomic clocks.According to the references, we can getσ
1
2
(
τ
), σ
2
2
(
τ
), and σ
3
2
(
τ
) from Δ
t
12, Δ
t
13, and Δ
t
23 in Figure 1. For example, σ
1
2
(
τ
) can be expressed as(19)σ12τ=12σ122τ+σ132τ-σ232τ.Because Allan variance is convergent for the five kinds of noises at different averaging time, it is often used to evaluate the stability of clock. However Allan variance cannot rule out frequency drift. Especially when the drift is almost equal to Allan variance for certain averaging time, if we use Allan variance to calculateσ
12
2
(
τ
), σ
13
2
(
τ
), and σ
23
2
(
τ
) and then calculate σ
1
2
(
τ
), σ
2
2
(
τ
), and σ
3
2
(
τ
), the result cannot reflect the real situation of atomic clock. In order to avoid the influence of frequency drift, two-order difference for frequency data or three-order difference for phase data is needed, which is just the definition of Hadamard variance.In order to make full use of the measurement data and also track slow change of satellite clock timely, we use windowed overlapping Hadamard variance, as shown in (20); even though the additional overlapping differences are not all statistically independent, they nevertheless increase the number of degrees of freedom and thus improve the confidence in the estimation. Moreover, by using the latest data, the variance can evaluate the health condition of atomic clock in real time:(20)Hσy2n=16N-3kτ2·∑i=n-N+1n-3kxi+3k-3xi+2k+3xi+k-xi2,where N is the length of the window, namely, the amount of data used for each update, and b
τ
0 is the time interval between x
i and x
i
+
1, which is defined as sampling period. τ
0 is the measurement period, averaging time τ
=
m
τ
0, and k
=
m
/
b. We know that, from the characteristics of Hadamard variance, the longer the averaging time, the larger the amount of data needed.Expression (21) can be used as recursive algorithm to reduce computation complexity in the updating of H
σ
y
2
(
n
). Moreover, 6
(
N
-
3
m
)
τ
2 is a constant for each averaging time; division operation can be done only when necessary:(21)Hσy2n=Hσy2n-1+Δ1n-Δ2n6N-3mτ2,where Δ
1
n
=
x
n
-
3
x
n
-
m
+
3
x
n
-
2
m
-
x
n
-
3
m
2 and Δ
2
n
=
x
n
-
N
+
3
m
-
3
x
n
-
N
+
2
m
+
3
x
n
-
N
+
m
-
x
n
-
N
2.As opposed to the Allan variance, which makes use of a second difference, the Hadamard variance employs a third difference that leads to reduction in the degrees of freedom by one. The Hadamard variance requires more data to produce a single stability calculation, as compared to the Allan variance, given equal averaging timeτ. So it will be a better choice to use different statistical tool for different averaging time. When the averaging time is short, it is much more convenient to use an ADEV three-cornered hat method, while HDEV will be a better choice when the linear frequency drift is dominant.


### 2.2.7. Detection for Frequency Drift-Rate Anomaly of Clock

According to [20–24], several drift-rate estimators are discussed and compared. We will firstly compare six different estimators by simulations in the following:(a)
Two points:z
^
1
=
(
y
(
n
)
-
y
(
1
)
)
/
n
-
1
τ
0
.
(b)
Two groups of points:z
^
2
=
(
2
/
n
τ
)
[
(
2
/
n
)
∑
i
=
n
/
2
+
1
n
y
(
i
)
-
(
2
/
n
)
∑
i
=
1
n
/
2
y
(
i
)
]
.
(c)
LS:z
^
3
=
(
6
/
n
n
2
-
1
τ
0
)
∑
i
=
1
n
(
2
i
-
n
-
1
)
y
(
i
)
.
(d)
Three points:z
^
4
=
(
x
2
n
+
1
-
2
x
1
+
n
+
x
1
)
/
(
n
τ
0
)
2
.
(e)
ω
4: z
^
5
=
(
6
/
N
3
τ
0
2
r
1
1
-
r
1
)
w
N
-
w
0
-
(
w
N
-
n
1
-
w
n
1
)
/
(
1
-
2
r
1
)
.
(f)
Kalman filter:X
k
=
Φ
X
k
-
1
+
ε
k, Z
k
=
H
X
k
+
ν
k,where X
k
=
x
k
,
y
k
,
z
k is the phase data, frequency deviation and drift rate, and observation matrix H
=
1,0
,
0. The state transition matrix is (22)ϕ=1ττ2201τ001and state error covariance matrix Q can be expressed as [25].To compare these six estimators, we generate simulation data by well-known Stable 32 software. The parameters of phase data are shown in Table11, and we consider three kinds of noise type which are WFM, FFM, and RWFM.Table 11
Parameters of generated phase data.
σ
0
σ
-
1
σ
-
2
Drift rate/day
Frequency deviation
1
E
-
12
1
E
-
14
1
E
-
16
1
E
-
12
1
E
-
12From Table12 and Figures 13 and 14, we know that ω
4 method and Kalman filter should be good choice for drift-rate estimation. To make sure that the simulation result is reliable, we make use of different simulation data to test the performance of the six estimators, and the result is just similar.Table 12
Statistical property of the estimation result of six estimators.
Estimator
a
b
c
d
e
f
E
(
z
)
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
9.97
E
-
13
σ
(
z
)
3.48
E
-
13
2.40
E
-
13
2.51
E
-
13
1.74
E
-
13
1.58
E
-
13
3.71
E
-
14Figure 13
Comparison of different drift-rate estimators.
(a)
Drift-rate estimation of different estimators
(b)
Partial enlarged view of (a)Figure 14
The drift-rate estimation of Kalman filter.In fact, no matter which method we choose to evaluate the drift rate, we must know the time error datax
k, which is equivalent to Δ
t in this section.When the satellite is within the sight of ground station, we can calculate the drift rate with certain estimator by comparing the satellite time with the timescale on the ground. But when the satellite cannot contact with the station, the only available data is the intercomparison dataΔ
t
12, Δ
t
13, and Δ
t
23.In order to evaluate the drift rate, there may be two methods. (1) Three-state Kalman filter can be used to evaluate the drift rate directly with the observation quantity Δ
t
12, Δ
t
13, and Δ
t
23, which is similar to the Kalman timescale algorithm on the ground. (2) Two steps are needed. The first step is to predict Δ
t
1, Δ
t
2, and Δ
t
3 from Δ
t
12, Δ
t
13, and Δ
t
23. And the second step is to evaluate drift rate with Δ
t
^
1, Δ
t
^
2, and Δ
t
^
3 (Δ
t
^
1, Δ
t
^
2, and Δ
t
^
3 are the prediction values of Δ
t
1, Δ
t
2, and Δ
t
3).In the following, we will compare these two methods. The three-state Kalman filter method is called method 1, and the combination of two-state Kalman filter withω
4 is called method 2.The basic Kalman fiter equations are as follows.System equation is(23)Xk=ΦXk-1+εk.Observation equation is(24)Zk=HXk+νk.For the two-state Kalman filter,(25)Xk=x1,k,y1,k,x2,k,y2,k,x3,k,y3,kT,Φ=diag⁡ϕ,ϕ,ϕ,ϕ=1τ01,Zk=Δt12,k,Δt13,kT,H=10-10001000-10.The system error covariance matrix Q can be expressed as [26–28].For the three-state Kalman filter,(26)Xk=x1,k,y1,k,z1,k,x2,k,y2,k,z2,k,x3,k,y3,k,z3,kT,Φ=diag⁡ϕ,ϕ,ϕ,ϕ=1ττ2201τ001,Zk=Δt12,k,Δt13,kT,H=100-100000100000-100,and we can get the system error covariance matrix Q according to [29].To compare these two methods, we make use of Stable 32 to generate simulation data. The parameters of the phase data are shown in Table13.Table 13
Parameter setting.
Parameters
σ
0
σ
-
1
σ
-
2
Drift rate/day
Frequency deviation
Clock 1
1
e
-
12
2
e
-
14
3
e
-
16
1
e
-
13
1
e
-
12
Clock 2
2
e
-
12
3
e
-
14
1
e
-
16
1
e
-
12
1
e
-
12
Clock 3
3
e
-
12
1
e
-
14
2
e
-
16
-
2
e
-
11
1
e
-
12Firstly, we will show the prediction performance of both two-state and three-state Kalman filter forΔ
t
1, Δ
t
2, and Δ
t
3; they are similar, which can be shown in Figure 15.Figure 15
Prediction performance of Kalman filter forΔ
t
1 and Δ
t
12.
(a)
Δ
t
12 and Δ
t
^
12
(b)
Δ
t
1 and Δ
t
^
1For method 2, after we have got the prediction ofΔ
t
1, Δ
t
2, and Δ
t
3, we will use ω
4 to evaluate the drift rate.Figure15(a) is the comparison of the clock difference Δ
t
12
(
Δ
t
12
=
Δ
t
1
-
Δ
t
2
) and its prediction Δ
t
^
12. Figure 15(b) compares the real time errors Δ
t
1 and Δ
t
^
1. From Figure 15, we know that the prediction error ε  (
ε
=
Δ
t
1
-
Δ
t
^
1
) will increase gradually, while Δ
t
^
12 is unbiased. In fact, this phenomenon is inevitable owing to the lack of absolute standard reference. No matter which method we choose, the prediction error will increase with time going. That also indicates that we cannot get precise drift rate, but it does not mean we can do nothing. Next we will prove that method 1 can give an alarm when the drift rate of one clock is much larger than the best one. It should be noted that when we make use of ω
4 method to evaluate drift rate for method 2, a sliding window is used to evaluate the drift rate in real time instead of batch processing.Figures16 and 17 show the drift-rate estimation of these two methods for three clocks. From the simulation result, we know that if the drift rate of one clock is much larger than the best clock, its drift-rate estimation is very close to its real value, and its computed estimation is much larger than the best clock, which enables method 1 to give an alarm. Besides, method 1 gives more precise and stable computation result compared to method 2.Figure 16
Drift-rate estimation of method 1.
(a)
Drift rate of clock 1
(b)
Drift rate of clock 2
(c)
Drift rate of clock 3Figure 17
Drift-rate estimation of method 2.
(a)
Drift rate of clock 1
(b)
Drift rate of clock 2
(c)
Drift rate of clock 3


## 2.2.1. Allan Variance

We usually use Allan variance [15, 16] to evaluate the stability of atomic clock; it can be expressed as follows:(12)σy2τ=12M-1∑i=1M-1y¯i+1m-y¯im2,yi=xi+1-xiτ=fo-frfr,y¯im=1m∑j=1myj,where τ
=
m
τ
0 is the averaging time and M is the amount of y
¯
i
(
m
). What needs to be pointed out is that y
i is the relative frequency deviation. f
o is the instantaneous frequency, f
r is the nominal frequency, and x
i is the clock time error at the ith measuring instant.Power-law spectrum is used to analyze noise property in frequency domain:(13)Syf=h-2f-2+h-1f-1+h0+h1f+h2f2=12π2∑α=-2α=2hαfα,where S
y
(
f
) is spectrum density for relative frequency deviation and h
α is the amplitude corresponding to different noise type. The power-law spectrum model contains five kinds of noise (α
=
-
2
,
-
1,0
,
1,2); they are RW FM, Flicker FM, White FM, Flicker PM, and White PM. σ
y
2
(
τ
) is determined by these five kinds of noise:(14)σy2τ=3fhh22πτ2+1.038+3ln⁡2πfhτh12πτ2+h02τ+2ln⁡2h-1+2π2τh-23.The slope of Allan variance gives us a knowledge of noise distribution in different averaging time.After obtaining a sufficient number of measurement data, Allan variance can be used to calculate the stability of different averaging time. But when anomaly occurs, the results given by Allan variance may lose practical significance. As shown in Figure6(a), the frequency of clock signal jumped and then returned some time later. We cannot get correct judgment for noise distribution according to the computed result by (12) shown in Figure 6(b). Besides, we do not know the anomaly type and detection delay is also too long.Figure 6
Relative frequency deviation and its Allan deviation when anomaly occurs.
(a)
Measured value of relative frequency deviation
(b)
Allan deviation


## 2.2.2. Dynamic Allan Variance (DAVAR)

As can be seen from Figure6(a), the main noise type does not change, but Figure 6(b) gives wrong judgment. Hence the conclusion is not consistent with the actual situation, and we cannot find the sign of frequency jumping either. Therefore, the traditional Allan variance cannot give believable information on such anomaly. In view of this, Galleani and Tavella put forward DAVAR, which can be expressed as (15) and can be used to evaluate the performance of clock in real time:(15)σy2n,k=12k2τ021N/k-1∑i=0N/k-2y¯i+1,k,n-y¯i,k,n2,y¯i,k,n=1k∑m=n-N+1+ikm=n-N+1+i+1kym.When we calculate DAVAR, a sliding window is used to cut the data. The window length isN, and σ
y
2
(
n
,
k
) will be updated when new measurement data y
n comes, so it can tell us health condition of clock in real time. τ
0 is the least measurement period, and k
τ
0 is the averaging time.


## 2.2.3. Modified DAVAR

From expression (15), we know that DAVAR can be updated in real time, but in order to guarantee the reliability of long-term stability, N must be large enough, which will greatly reduce the detection probability of instantaneous anomaly. Because when frequency jumping occurs, …
,
y
i
-
1
-
y
i
-
2
≈
0
,
y
i
-
y
i
-
1
=
δ
,
y
i
+
1
-
y
i
≈
0
,
…, only one factor is not 0, DAVAR is not sensitive enough to weak frequency jumping. In this paper, we modify Dynamic Allan Variance to improve detection sensitivity for small frequency jumping:(16)σy2n=12m∑i=1my¯i+n-m-y¯i+n-2m2,(17)σy2n=σy2n-1+y¯n-y¯n-m2-y¯n-m-y¯n-2m22m.Expression (16) is the Modified DAVAR, and (17) is its iterative calculation method.


## 2.2.4. Detection Performance of Modified DAVAR

In this section, we will firstly analyze and compare the detection performances of DAVAR and Modified DAVAR when facing phase and frequency jumping and then show that Modified DAVAR is also effective in detecting instantaneous stability deterioration.The monitoring method for phase and frequency jumping is based on statistics.Δ
t
12, Δ
t
13, and Δ
t
23 from Figure 1 will be used here. Assume that only one of the three clocks breaks down. If phase or frequency jumping occurs in clock 1, Δ
t
12, Δ
t
13 will be abnormal, while Δ
t
23 is still normal. Because the anomaly in Δ
t
12 is the same as that in Δ
t
13, we only need to analyze Δ
t
12.In Figure7(a), the amplitude of phase jumping is 12 times the standard deviation of the relative frequency deviation data y
i. In Figure 7(b), the frequency jumping is 4 times the standard deviation of y
i.Figure 7
Phase and frequency jumping.
(a)
Phase jumping
(b)
Frequency jumpingWe simulated 1000 sampling points and phase and frequency jumping occurred at the 500th point. We simulated one realization and saved the response data of DAVAR and Modified DAVAR to jumping at every time instant; then we repeated 10000 realizations in the same way. Of course the 1000 sampling points’ data is different in every realization. Then we got the average response at every time instant that is shown in Figures8 and 9.Figure 8
Comparison of two variances for phase jumping.
(a)
Average response of DAVAR to phase jumping
(b)
Average response of Modified DAVAR to phase jumpingFigure 9
Comparison of two variances for frequency jumping.
(a)
Average response of DAVAR to frequency jumping
(b)
Average response of Modified DAVAR to frequency jumpingBecause the peak value of response of DAVAR and Modified DAVAR to jumping determines if the jumping could be detected, we focus on the peak value in each realization. Assuming that, in one realization, the maximum value of response of DAVAR to frequency jumping isσ
i
,
m
2, the maximum value of response of Modified DAVAR to frequency jumping is υ
i
,
m
2, and then we saved σ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
). We studied the saved data to give the minimum value, maximum value, mean value, and standard deviation of σ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
) in Figures 8 and 9.To make the statistical result more clear, we list the statistical characteristic ofσ
i
,
m
2 and υ
i
,
m
2
(
i
=
1,2
,
3
,
…
,
10000
) in frequency jumping in Table 5. From Table 5 and Figure 9 we know that σ
i
,
m
2 is not big enough to be distinguished from the base noise, and Modified DAVAR is more sensitive to weak frequency jumping.Table 5
Statistical characteristic ofσ
i
,
m
2 and υ
i
,
m
2.
Statistical item
DAVAR
Modified DAVAR
Minimum value
1.60
e
-
23
3.09
e
-
23
Maximum value
5.59
e
-
23
1.49
e
-
22
Mean value
2.89
e
-
23
8.14
e
-
23
Standard deviation
5.01
e
-
24
1.70
e
-
23In Tables6 and 7, the detection performance of two kinds of variance for phase and frequency jumping is given. Tables 8 and 9 tell us the detection performance for different clock stability.Table 6
Detection performance for phase jumping.
Detection item
DAVAR
Modified DAVAR
Resolution
12
12
Threshold
4.5
E
-
23
4.5
E
-
23
PD
1.0000
0.9978
PFA
1.61
e
-
5
3.00
e
-
7
Detection delay (points)
1
1Table 7
Detection performance for phase jumping.
Detection item
DAVAR
Modified DAVAR
Resolution
12
4
Threshold
4.5
E
-
23
4.5
E
-
23
PD
0.9874
0.9927
PFA
1.61
e
-
5
3.00
e
-
7
Detection delay (points)
1
5Table 8
Detection performance of phase jumping for different clock stability.
σ
(
τ
0
)
Resolution
Threshold
PD
PFA
Delay (points)
3
E
-
12
12
4.5
E
-
23
0.9978
3.00
e
-
7
1
3
E
-
11
12
4.5
E
-
21
0.9961
5.00
e
-
7
1
3
E
-
10
12
4.5
E
-
19
0.9977
6.00
e
-
7
1Table 9
Detection performance of frequency jumping for different clock stability.
σ
(
τ
0
)
Resolution
Threshold
PD
PFA
Delay (points)
3
E
-
12
4
4.5
E
-
23
0.9927
3.00
e
-
7
5
3
E
-
11
4
4.5
E
-
21
0.9923
5.00
e
-
7
5
3
E
-
10
4
4.5
E
-
19
0.9934
6.00
e
-
7
5What should be pointed out is that the unit of resolution in phase and frequency jumping in Tables6–9 is the standard deviation of the relative frequency deviation data y
i, namely, σ
y. Detection delay is the number of sampling points from anomaly occurring to be detected.During the simulation, we simulate 10000 realizations and use same threshold for both DAVAR and Modified DAVAR. The window lengthN
=
10, and k
τ
0
=
τ
0
=
1 s. It should be noted that because k
τ
0
=
τ
0
=
1 s, it is precise enough for us to do the approximation to only consider WFM in the simulation. In fact, the Modified DAVAR is still effective in detecting phase and frequency jumping in the presence of other noise types.Table6 tells us that PD of Modified DAVAR is almost as good as DAVAR and PFA is lower. From Table 7, we know that Modified DAVAR is more sensitive to weak frequency jumping but detection delay is longer.From Tables8 and 9, we notice that the resolution of Modified DAVAR for phase and frequency jumping is the same for different clocks. What we need to do is only to reset the threshold according to expression (18):(18)Thr=5·σ21.In addition, the window lengthN is an important parameter for Modified DAVAR. The longer the window, the weaker the detection performance. However, if the window length is too short, PFA will rise and detection resolution will also deteriorate.Figure10 shows that the Modified DAVAR can also monitor instantaneous deterioration effectively.Figure 10
Instantaneous deterioration of clock.
(a)
Instantaneous deterioration
(b)
Modified DAVARModified DAVAR can be considered as a statistical tool; it is effective to detect phase and frequency jumping. Compared with PLL method, we need to measure the time error data firstly and then calculate the statistical characteristics of clock. Modified DAVAR can monitor weaker frequency jumping compared to PLL method, but PLL method is independent of a second standard reference and time-comparison device, which will give us more flexibility. Taking their respective characteristics into account, cooperation between them may be a good choice to improve the Self-Monitoring reliability.


## 2.2.5. LS Method and Kalman Method on Detection for Frequency Jumping

In this section, we will introduce two existing methods, which are called LS method and Kalman filter method.LS Method. We may make use of LS algorithm to calculate the averaging frequency deviation y
¯ with M newest saved sampling points and then predict the time error Δ
t
^
k
+
1
=
Δ
t
k
+
y
¯
·
τ. After we can compare it with the real measurement Δ
t
k
+
1, if the difference ε
=
Δ
t
^
k
+
1
-
Δ
t
k
+
1 is beyond the threshold γ, we think that frequency jumping occurs.Kalman Filter Method. We can make use of Kalman filter to predict the next state y
^
n of clock and then compare it with the real measurement y
n. If the difference ε
=
y
^
n
-
y
n is larger than the configurable threshold, we think that anomaly occurs.Simulations. Because the sampling interval τ
=
1 s, we only consider WFM noise, whose standard deviation is σ
0
=
3
E
-
12. During the simulations for LS method, we choose M
=
20, while, for Kalman filter method, the state transition matrix Φ
=
1, observation matrix H
=
1, system error covariance matrix Q
=
(
1.0
×
10
-
13
)
2, and observation error covariance matrix R
=
(
3.0
×
10
-
12
)
2.Figures11 and 12 show the frequency jumping detection of the two methods.Figure 11
LS method in detecting frequency jumping.
(a)
ε
=
Δ
t
^
k
+
1
-
Δ
t
k
+
1
(b)
y
- and yFigure 12
Kalman filter method in detecting frequency jumping.
(a)
ε
=
y
^
n
-
y
n
(b)
y
^
n and y
nAfter we have done numerical simulations, we give Table10 to show the detection performance of the methods, in which PLL method, DAVAR method, and Modified DAVAR method are included.Table 10
Detection performance of frequency jumping for different clock stability.
Detection method
Resolution
PD
PFA
Delay
PLL method
4.0
E
-
11
0.9933
2.41
E
-
6
0.4 s
DAVAR method
3.6
E
-
11
0.9874
1.61
E
-
5
1 s
Modified DAVAR method
1.2
E
-
11
0.9927
3.00
E
-
7
5 s
LS method
1.5
E
-
11
0.9944
2.53
E
-
5
1 s
Kalman filter method
1.8
E
-
11
0.9926
2.60
E
-
6
1 sDiscussion. It should be noted firstly that the detection performance will be different with different parameters. We use the same noise level to test different methods to give Table 10, which can be a reference to show different characteristics of different methods.Different observation quantity is needed for different methods. In our opinion, DAVAR, Modified DAVAR, LS, and Kalman filter are all effective in detecting weak frequency jumping. They need a second standard reference and time-comparison device to get the time error measurementΔ
t
k to be used as observation quantity to run the algorithm. PLL method can realize Self-Monitoring for frequency jumping without standard source, which can give us more flexibility and has been used on BeiDou satellite. The computation complexity is also different among them. LS method, Kalman filter method, and PLL method should be of less computation quantity. The resolution of Modified DAVAR is the best but at the cost of longer delay, while the PLL method has the shortest delay but at the cost of worst resolution. Sometimes we may have to compromise between resolution and delay.


## 2.2.6. Stability Evaluation of Satellite Clock

In this section, we use the well-known “three-cornered hat” approach [17–19] to evaluate the stability of atomic clocks.According to the references, we can getσ
1
2
(
τ
), σ
2
2
(
τ
), and σ
3
2
(
τ
) from Δ
t
12, Δ
t
13, and Δ
t
23 in Figure 1. For example, σ
1
2
(
τ
) can be expressed as(19)σ12τ=12σ122τ+σ132τ-σ232τ.Because Allan variance is convergent for the five kinds of noises at different averaging time, it is often used to evaluate the stability of clock. However Allan variance cannot rule out frequency drift. Especially when the drift is almost equal to Allan variance for certain averaging time, if we use Allan variance to calculateσ
12
2
(
τ
), σ
13
2
(
τ
), and σ
23
2
(
τ
) and then calculate σ
1
2
(
τ
), σ
2
2
(
τ
), and σ
3
2
(
τ
), the result cannot reflect the real situation of atomic clock. In order to avoid the influence of frequency drift, two-order difference for frequency data or three-order difference for phase data is needed, which is just the definition of Hadamard variance.In order to make full use of the measurement data and also track slow change of satellite clock timely, we use windowed overlapping Hadamard variance, as shown in (20); even though the additional overlapping differences are not all statistically independent, they nevertheless increase the number of degrees of freedom and thus improve the confidence in the estimation. Moreover, by using the latest data, the variance can evaluate the health condition of atomic clock in real time:(20)Hσy2n=16N-3kτ2·∑i=n-N+1n-3kxi+3k-3xi+2k+3xi+k-xi2,where N is the length of the window, namely, the amount of data used for each update, and b
τ
0 is the time interval between x
i and x
i
+
1, which is defined as sampling period. τ
0 is the measurement period, averaging time τ
=
m
τ
0, and k
=
m
/
b. We know that, from the characteristics of Hadamard variance, the longer the averaging time, the larger the amount of data needed.Expression (21) can be used as recursive algorithm to reduce computation complexity in the updating of H
σ
y
2
(
n
). Moreover, 6
(
N
-
3
m
)
τ
2 is a constant for each averaging time; division operation can be done only when necessary:(21)Hσy2n=Hσy2n-1+Δ1n-Δ2n6N-3mτ2,where Δ
1
n
=
x
n
-
3
x
n
-
m
+
3
x
n
-
2
m
-
x
n
-
3
m
2 and Δ
2
n
=
x
n
-
N
+
3
m
-
3
x
n
-
N
+
2
m
+
3
x
n
-
N
+
m
-
x
n
-
N
2.As opposed to the Allan variance, which makes use of a second difference, the Hadamard variance employs a third difference that leads to reduction in the degrees of freedom by one. The Hadamard variance requires more data to produce a single stability calculation, as compared to the Allan variance, given equal averaging timeτ. So it will be a better choice to use different statistical tool for different averaging time. When the averaging time is short, it is much more convenient to use an ADEV three-cornered hat method, while HDEV will be a better choice when the linear frequency drift is dominant.


## 2.2.7. Detection for Frequency Drift-Rate Anomaly of Clock

According to [20–24], several drift-rate estimators are discussed and compared. We will firstly compare six different estimators by simulations in the following:(a)
Two points:z
^
1
=
(
y
(
n
)
-
y
(
1
)
)
/
n
-
1
τ
0
.
(b)
Two groups of points:z
^
2
=
(
2
/
n
τ
)
[
(
2
/
n
)
∑
i
=
n
/
2
+
1
n
y
(
i
)
-
(
2
/
n
)
∑
i
=
1
n
/
2
y
(
i
)
]
.
(c)
LS:z
^
3
=
(
6
/
n
n
2
-
1
τ
0
)
∑
i
=
1
n
(
2
i
-
n
-
1
)
y
(
i
)
.
(d)
Three points:z
^
4
=
(
x
2
n
+
1
-
2
x
1
+
n
+
x
1
)
/
(
n
τ
0
)
2
.
(e)
ω
4: z
^
5
=
(
6
/
N
3
τ
0
2
r
1
1
-
r
1
)
w
N
-
w
0
-
(
w
N
-
n
1
-
w
n
1
)
/
(
1
-
2
r
1
)
.
(f)
Kalman filter:X
k
=
Φ
X
k
-
1
+
ε
k, Z
k
=
H
X
k
+
ν
k,where X
k
=
x
k
,
y
k
,
z
k is the phase data, frequency deviation and drift rate, and observation matrix H
=
1,0
,
0. The state transition matrix is (22)ϕ=1ττ2201τ001and state error covariance matrix Q can be expressed as [25].To compare these six estimators, we generate simulation data by well-known Stable 32 software. The parameters of phase data are shown in Table11, and we consider three kinds of noise type which are WFM, FFM, and RWFM.Table 11
Parameters of generated phase data.
σ
0
σ
-
1
σ
-
2
Drift rate/day
Frequency deviation
1
E
-
12
1
E
-
14
1
E
-
16
1
E
-
12
1
E
-
12From Table12 and Figures 13 and 14, we know that ω
4 method and Kalman filter should be good choice for drift-rate estimation. To make sure that the simulation result is reliable, we make use of different simulation data to test the performance of the six estimators, and the result is just similar.Table 12
Statistical property of the estimation result of six estimators.
Estimator
a
b
c
d
e
f
E
(
z
)
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
1.00
E
-
12
9.97
E
-
13
σ
(
z
)
3.48
E
-
13
2.40
E
-
13
2.51
E
-
13
1.74
E
-
13
1.58
E
-
13
3.71
E
-
14Figure 13
Comparison of different drift-rate estimators.
(a)
Drift-rate estimation of different estimators
(b)
Partial enlarged view of (a)Figure 14
The drift-rate estimation of Kalman filter.In fact, no matter which method we choose to evaluate the drift rate, we must know the time error datax
k, which is equivalent to Δ
t in this section.When the satellite is within the sight of ground station, we can calculate the drift rate with certain estimator by comparing the satellite time with the timescale on the ground. But when the satellite cannot contact with the station, the only available data is the intercomparison dataΔ
t
12, Δ
t
13, and Δ
t
23.In order to evaluate the drift rate, there may be two methods. (1) Three-state Kalman filter can be used to evaluate the drift rate directly with the observation quantity Δ
t
12, Δ
t
13, and Δ
t
23, which is similar to the Kalman timescale algorithm on the ground. (2) Two steps are needed. The first step is to predict Δ
t
1, Δ
t
2, and Δ
t
3 from Δ
t
12, Δ
t
13, and Δ
t
23. And the second step is to evaluate drift rate with Δ
t
^
1, Δ
t
^
2, and Δ
t
^
3 (Δ
t
^
1, Δ
t
^
2, and Δ
t
^
3 are the prediction values of Δ
t
1, Δ
t
2, and Δ
t
3).In the following, we will compare these two methods. The three-state Kalman filter method is called method 1, and the combination of two-state Kalman filter withω
4 is called method 2.The basic Kalman fiter equations are as follows.System equation is(23)Xk=ΦXk-1+εk.Observation equation is(24)Zk=HXk+νk.For the two-state Kalman filter,(25)Xk=x1,k,y1,k,x2,k,y2,k,x3,k,y3,kT,Φ=diag⁡ϕ,ϕ,ϕ,ϕ=1τ01,Zk=Δt12,k,Δt13,kT,H=10-10001000-10.The system error covariance matrix Q can be expressed as [26–28].For the three-state Kalman filter,(26)Xk=x1,k,y1,k,z1,k,x2,k,y2,k,z2,k,x3,k,y3,k,z3,kT,Φ=diag⁡ϕ,ϕ,ϕ,ϕ=1ττ2201τ001,Zk=Δt12,k,Δt13,kT,H=100-100000100000-100,and we can get the system error covariance matrix Q according to [29].To compare these two methods, we make use of Stable 32 to generate simulation data. The parameters of the phase data are shown in Table13.Table 13
Parameter setting.
Parameters
σ
0
σ
-
1
σ
-
2
Drift rate/day
Frequency deviation
Clock 1
1
e
-
12
2
e
-
14
3
e
-
16
1
e
-
13
1
e
-
12
Clock 2
2
e
-
12
3
e
-
14
1
e
-
16
1
e
-
12
1
e
-
12
Clock 3
3
e
-
12
1
e
-
14
2
e
-
16
-
2
e
-
11
1
e
-
12Firstly, we will show the prediction performance of both two-state and three-state Kalman filter forΔ
t
1, Δ
t
2, and Δ
t
3; they are similar, which can be shown in Figure 15.Figure 15
Prediction performance of Kalman filter forΔ
t
1 and Δ
t
12.
(a)
Δ
t
12 and Δ
t
^
12
(b)
Δ
t
1 and Δ
t
^
1For method 2, after we have got the prediction ofΔ
t
1, Δ
t
2, and Δ
t
3, we will use ω
4 to evaluate the drift rate.Figure15(a) is the comparison of the clock difference Δ
t
12
(
Δ
t
12
=
Δ
t
1
-
Δ
t
2
) and its prediction Δ
t
^
12. Figure 15(b) compares the real time errors Δ
t
1 and Δ
t
^
1. From Figure 15, we know that the prediction error ε  (
ε
=
Δ
t
1
-
Δ
t
^
1
) will increase gradually, while Δ
t
^
12 is unbiased. In fact, this phenomenon is inevitable owing to the lack of absolute standard reference. No matter which method we choose, the prediction error will increase with time going. That also indicates that we cannot get precise drift rate, but it does not mean we can do nothing. Next we will prove that method 1 can give an alarm when the drift rate of one clock is much larger than the best one. It should be noted that when we make use of ω
4 method to evaluate drift rate for method 2, a sliding window is used to evaluate the drift rate in real time instead of batch processing.Figures16 and 17 show the drift-rate estimation of these two methods for three clocks. From the simulation result, we know that if the drift rate of one clock is much larger than the best clock, its drift-rate estimation is very close to its real value, and its computed estimation is much larger than the best clock, which enables method 1 to give an alarm. Besides, method 1 gives more precise and stable computation result compared to method 2.Figure 16
Drift-rate estimation of method 1.
(a)
Drift rate of clock 1
(b)
Drift rate of clock 2
(c)
Drift rate of clock 3Figure 17
Drift-rate estimation of method 2.
(a)
Drift rate of clock 1
(b)
Drift rate of clock 2
(c)
Drift rate of clock 3


## 3. Conclusion

This paper puts forward a set of Self-Monitoring Methods for common anomalies. We use PLL to realize Self-Monitoring for signal loss and phase and frequency jumping. Based on the measurement data from intercomparison among three clocks, Modified DAVAR is used to detect phase and frequency jumping and instantaneous deterioration; we use windowed overlapping Hadamard variance to evaluate clock stability in real time and the three-state Kalman filter for large drift-rate anomaly.The method based on PLL has been proved effective and used in newest BeiDou satellite. And the other research on Self-Monitoring Method in this paper can be used in next generation navigation satellites after year 2019.



---

*Source: 1759512-2016-07-28.xml*