# Research on self-monitoring method for anomalies of satellite atomic clock.

1. IntroductionThe 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 j udge 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 Figure 1.

Firstly, we define that [DELTA][t.sub.1] is the time error of clock 1, which is the difference between clock time and the standard time. [DELTA][t.sub.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 [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.23] among them. The Signal Processing Module uses the time difference data [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.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, [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.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 Figure 1 can be described in Figure 2.

2.1. Self-Monitoring Method Based on PLL. Figure 2 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[pi][f.sub.r]t + [[phi].sub.1](t) and that of Direct Digital Synthesizer (DDS) output is 2[pi][f.sub.r]t + [[phi].sub.2](t), then we get

K ([[phi].sub.1] (s) - [[phi].sub.2] (s)) F(s) 1/s = [[phi].sub.2] (s), (1)

where F(s) = (1 + s[[tau].sub.2])/s[[tau].sub.1] is the transfer function of the twoorder ideal loop filter, 1/s is the normalized transfer function of the DDS, and K is the loop gain. From expression (1), we get

[[phi].sub.e](s) = [s/s + KF(s)] [[phi].sub.1](s), (2)

where [[phi].sub.e](s) = [[phi].sub.1](s) - [[phi].sub.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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

where [[omega].sub.n] = [square root of (K/[[tau].sub.1])] is the undamped oscillation frequency and [xi] = ([[tau].sub.2]/2) [square root of (K/[[tau].sub.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 [[phi].sub.1](t) = [DELTA][phi] x [epsilon](t), whose Laplace transform can be expressed as [[phi].sub.1](s) = [DELTA][phi]/s, then the error response is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

Through factorization, (4) is equivalent to

[[phi].sub.e](s) = A/s - [s.sub.1] + B/s - [s.sub.2], (5)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

Considering (6) and (7), the inverse Laplace transform of (5) can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

From (8) we notice that the phase difference at t = 0 reaches its peakvalue 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 = 500 s which leads to obvious jumping in the output of Phase Detector. In Figure 3, the loop parameters of three PLLs are [[omega].sub.n] = 1,4, 10 and [xi] = [square root of 2]/2 and the amplitude of phase jumping is equal to 1/[10.sup.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.

2.1.2. Frequency Jumping. Assuming that frequency jumping is [[phi].sub.2](t) = [DELTA][phi]t x [epsilon](t), whose Laplace transform is [[phi].sub.2](s) = [DELTA][phi]/[s.sup.2], then the error response can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

After factorization,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

[s.sub.1] and [s.sub.2] can also be described by expression (6). According to (6) and (10), the inverse Laplace transform of (5) can be expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

As can be seen from (11), the maximum amplitude of the phase difference in tracking is inversely proportional to [[omega].sub.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.sup.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.

2.1.3. Signal Loss. As shown in Figure 5, 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.

2.1.4. Summary. It can be seen that, from Figures 3, 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, [xi] = [square root of 2]/2, [[omega].sub.n] = 4, and we use two-order ideal loop filter. Assuming that the relative frequency deviation [y.sub.i] (12) of clock signal follows Gauss distribution, whose Allan deviation can be expressed as 3E - 12/[square root of [tau]], 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/[10.sup.8] period phase jumping and 4/[10.sup.11] frequency jumping. Detection delay is defined as [DELTA]t = m x 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.

Simulation 2. During the second simulation, detection threshold and resolution change with the stability of atomic clock. Tables 3 and 4 show it.

Analysis. When 0 < [xi] < 1, PLL is called underdamped system, in which phase and frequency jumping will result in drastic oscillation. If [xi] > 1, PLL is overdamped and usually more stable and slow to anomaly. In practice, we often set [xi] = 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 [[omega].sub.n] is too large, and suppression for noise will also become weaker. Conversely, if [[omega].sub.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, [[omega].sub.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 Tables 3 and 4. The relationship between resolution and stability can be described as Re s(p) [approximately equal to] (1E4/3) x [sigma](1) and Re s(f) [approximately equal to] (4E1/3) x [sigma](1), while threshold can be set as Thr [approximately equal to] (11E3/3) x [sigma](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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

where [tau] = m[[tau].sub.0] is the averaging time and M is the amount of [[bar.y].sub.i](m). What needs to be pointed out is that [y.sub.i] is the relative frequency deviation. [f.sub.o] is the instantaneous frequency, [f.sub.r] is the nominal frequency, and [x.sub.i] is the clock time error at the ith measuring instant.

Power-law spectrum is used to analyze noise property in frequency domain:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where [S.sub.y](f) is spectrum density for relative frequency deviation and [h.sub.[alpha]] is the amplitude corresponding to different noise type. The power-law spectrum model contains five kinds of noise ([alpha] = -2, -1,0, 1, 2); they are RW FM, Flicker FM, White FM, Flicker PM, and White PM. [[sigma].sup.2.sub.y]([tau]) is determined by these five kinds of noise:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

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 Figure 6(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.

2.2.2. Dynamic Allan Variance (DAVAR). As can be seen from Figure 6(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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

When we calculate DAVAR, a sliding window is used to cut the data. The window length is N, and [[sigma].sup.2.sub.y](n, k) will be updated when new measurement data [y.sub.n] comes, so it can tell us health condition of clock in real time. [[tau].sub.0] is the least measurement period, and k[[tau].sub.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.sub.i-1] - [y.sub.i-2] [approximately equal to] 0, [y.sub.i] - [y.sub.i-1] = [delta], [y.sub.i+1] - [y.sub.i] [approximately equal to] 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

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. [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.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, [DELTA][t.sub.12], [DELTA][t.sub.13] will be abnormal, while [DELTA][t.sub.23] is still normal. Because the anomaly in [DELTA][t.sub.12] is the same as that in [DELTA][t.sub.13], we only need to analyze [DELTA][t.sub.12].

In Figure 7(a), the amplitude of phase jumping is 12 times the standard deviation of the relative frequency deviation data [y.sub.i]. In Figure 7(b), the frequency jumping is 4 times the standard deviation of [y.sub.i].

We 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 Figures 8 and 9.

Because 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 [[sigma].sup.2.sub.i,m], the maximum value of response of Modified DAVAR to frequency jumping is [v.sup.2.sub.i,m], and then we saved [[sigma].sup.2.sub.i,m] and [v.sup.2.sub.i,m] (i = 1, 2, 3, ..., 10000). We studied the saved data to give the minimum value, maximum value, mean value, and standard deviation of [[sigma].sup.2.sub.i,m] and [v.sup.2.sub.i,m] (i = 1, 2,3, ..., 10000) in Figures 8 and 9.

To make the statistical result more clear, we list the statistical characteristic of [[sigma].sup.2.sub.i,m] and [v.sup.2.sub.i,m] (i = 1, 2, 3, ..., 10000) in frequency jumping in Table 5. From Table 5 and Figure 9 we know that [[sigma].sup.2.sub.i,m] is not big enough to be distinguished from the base noise, and Modified DAVAR is more sensitive to weak frequency jumping.

In Tables 6 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.

What should be pointed out is that the unit of resolution in phase and frequency jumping in Tables 6-9 is the standard deviation of the relative frequency deviation data [y.sub.i], namely, [[sigma].sub.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 length N = 10, and k[[tau].sub.0] = [[tau].sub.0] = 1 s. It should be noted that because k[[tau].sub.0] = [[tau].sub.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.

Table 6 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 Tables 8 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):

Thr = 5 x [[sigma].sup.2] (1). (18)

In addition, the window length N 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.

Figure 10 shows that the Modified DAVAR can also monitor instantaneous deterioration effectively.

Modified 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 [bar.y] with M newest saved sampling points and then predict the time error [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. After we can compare it with the real measurement [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], if the difference [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is beyond the threshold y, we think that frequency jumping occurs.

Kalman Filter Method. We can make use of Kalman filter to predict the next state [[??].sub.n] of clock and then compare it with the real measurement [y.sub.n]. If the difference [epsilon] = [[??].sub.n] - [y.sub.n] is larger than the configurable threshold, we think that anomaly occurs.

Simulations. Because the sampling interval [tau] = 1 s, we only consider WFM noise, whose standard deviation is [[sigma].sub.0] = 3E - 12. During the simulations for LS method, we choose M = 20, while, for Kalman filter method, the state transition matrix [PHI] = 1, observation matrix H = 1, system error covariance matrix Q = [(1.0 x [10.sup.-13]).sup.2], and observation error covariance matrix R = [(3.0 x [10.sup.-12]).sup.2].

Figures 11 and 12 show the frequency jumping detection of the two methods.

After we have done numerical simulations, we give Table 10 to show the detection performance of the methods, in which PLL method, DAVAR method, and Modified DAVAR method are included.

Discussion. 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 [DELTA][t.sub.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 [[sigma].sup.2.sub.1]([tau]), [[sigma].sup.2.sub.2]([tau]), and [[sigma].sup.2.sub.3]([tau]) from [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.23] in Figure 1. For example, [[sigma].sup.2.sub.1]([tau]) can be expressed as

[[sigma].sup.2.sub.1]([tau]) = 1/2 [[[sigma].sup.2.sub.12]([tau]) + [[sigma].sup.2.sub.13]([tau]) - [[sigma].sup.2.sub.23]([tau])]. (19)

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 [[sigma].sup.2.sub.12]([tau]), [[sigma].sup.2.sub.13]([tau]), and [[sigma].sup.2.sub.23]([tau]) and then calculate [[sigma].sup.2.sub.1]([tau]), [[sigma].sup.2.sub.2]([tau]), and [[sigma].sup.2.sub.3]([tau]), 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

where N is the length of the window, namely, the amount of data used for each update, and b[[tau].sub.0] is the time interval between [x.sub.i] and [x.sub.i+1], which is defined as sampling period. [[tau].sub.0] is the measurement period, averaging time [tau] = m[[tau].sub.0], and k = m/fc. 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[[sigma].sup.2.sub.y](m). Moreover, 6(N - 3m)[[tau].sup.2] is a constant for each averaging time; division operation can be done only when necessary:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

where [[DELTA].sub.1](n) = [([x.sub.n] - 3[x.sub.n-m] + 3[x.sub.n-2m] - [x.sub.n-3m]).sup.2] and [[DELTA].sub.2](n) = [([x.sub.n-N+3m] - 3[x.sub.n-N+2m] + 3[x.sub.n-N+m] - [x.sub.n-N]).sup.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 [tau]. 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: [[??].sub.1] = (y(n) - y(l))/(n - 1)[[tau].sub.0].

(b) Two groups of points: [[??].sub.2] = (2/n[tau])[(2/n) [[summation].sup.n.sub.i=n/2+1] y(i)z - (2/n) [[summation].sup.n/2.sub.i=1] y(i)].

(c) LS: [[??].sub.3] = (6/n([n.sup.2] - 1)[[tau].sub.0]) [[summation].sup.n.sub.i=1](2i - n - 1)y(i).

(d) Three points: [z.sub.4] = (x(2n+1)-2x(1+n)+x(1))/[(n[[tau].sub.0]).sup.2].

(e) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

(f) Kalman filter: [X.sub.k] = [PHI][X.sub.k-1] + [[epsilon].sub.k], [Z.sub.k] = H[X.sub.k] + [v.sub.k],

where [X.sub.k] = [[x.sub.k], [y.sub.k], [z.sub.k]] is the phase data, frequency deviation and drift rate, and observation matrix H = [l, 0, 0]. The state transition matrix is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

and 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 Table 11, and we consider three kinds of noise type which are WFM, FFM, and RWFM.

From Table 12 and Figures 13 and 14, we know that [[omega].sub.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.

In fact, no matter which method we choose to evaluate the drift rate, we must know the time error data [x.sub.k], which is equivalent to [DELTA]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 [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.23].

In order to evaluate the drift rate, there may be two methods. (1) Three-state Kalman filter canbe used to evaluate the drift rate directly with the observation quantity [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.23], which is similar to the Kalman timescale algorithm on the ground. (2) Two steps are needed. The first step is to predict [DELTA][t.sub.1], [DELTA][t.sub.2], and [DELTA][t.sub.3] from [DELTA][t.sub.12], [DELTA][t.sub.13], and [DELTA][t.sub.23]. And the second step is to evaluate drift rate with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are the prediction values of [DELTA][t.sub.1], [DELTA][t.sub.2], and [DELTA][t.sub.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 [[omega].sub.4] is called method 2.

The basic Kalman fiter equations are as follows.

System equation is

[X.sub.k] = [PHI][X.sub.k-1] + [[epsilon].sub.k]. (23)

Observation equation is

[Z.sub.k] = H[X.sub.k] + [v.sub.k]. (24)

For the two-state Kalman filter,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

The system error covariance matrix Q can be expressed as [26-28].

For the three-state Kalman filter,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (26)

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 Table 13.

Firstly, we will show the prediction performance of both two-state and three-state Kalman filter for [DELTA][t.sub.1], [DELTA][t.sub.2], and [DELTA][t.sub.3]; they are similar, which can be shown in Figure 15.

For method 2, after we have got the prediction of [DELTA][t.sub.1], [DELTA][t.sub.2], and [DELTA][t.sub.3], we will use [[omega].sub.4] to evaluate the drift rate.

Figure 15(a) is the comparison of the clock difference [DELTA][t.sub.12] ([DELTA][t.sub.12] = [DELTA][t.sub.1] - [DELTA][t.sub.2]) and its prediction [??][t.sub.12]. Figure 15(b) compares the real time errors [DELTA][t.sub.1] and [DELTA][[??].sub.1]. From Figure 15, we know that the prediction error [epsilon] ([epsilon] = [DELTA][t.sub.1] - [DELTA][[??].sub.1]) will increase gradually, while [DELTA][[??].sub.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 [[omega].sub.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.

Figures 16 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.

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.

http://dx.doi.org/10.1155/2016/1759512

Competing Interests

The authors declare that they have no competing interests.

References

[1] Y. C. Chan, W. A. Johnson, S. K. Karuza, A. M. Young, and J. C. Camparo, "Self-monitoring and self-assessing atomic clocks," IEEE Transactions on Instrumentation and Measurement, vol. 59, no. 2, pp. 330-334, 2010.

[2] Y. Liu and S. Tang, "Studying method of monitoring composite clocks performance," in Proceedings of the 3th China Satellite Navigation Conference (CSNC '12), Guangzhou, China, June 2012.

[3] S. Tang, Y. Liu, and X.-H. Li, "A study on onboard satellite atomic clock autonomous integrity monitoring," Journal of Astronautics, vol. 34, no. 1, pp. 39-45, 2013.

[4] E. Nunzi, P. Carbone, and P. Tavella, "Fault detection in atomic clock frequency standards affected by mean and variance changes and by an additive periodic component: the GLRT approach," in Proceedings of the 2008 IEEE International Instrumentation and Measurement Technology Conference (IMTC '08), vol. 34, no 1, pp. 1594-1597, IEEE, Victoria, Canada, May 2008.

[5] E. Nunzi, L. Galleani, P Tavella, and P Carbone, "Detection of anomalies in the behavior of atomic clocks," IEEE Transactions on Instrumentation and Measurement, vol. 56, no. 2, pp. 523-528, 2007.

[6] U. Bartoccini, G. Barchi, and E. Nunzi, "Methods and tools for frequency jump detection," in Proceedings of the IEEE International Workshop on Advanced Methods for Uncertainty Estimation in Measurement (AMUEM '09), pp. 109-112, Bucharest, Romania, July 2009.

[7] L. Galleani and P Tavella, "Detection of atomic clock frequency jumps with the Kalman filter," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 59, no. 3, pp. 504-509, 2012.

[8] X. Huang, H. Gong, and G. Ou, "Detection of weak frequency jumps for GNSS onboard clocks," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, no. 5, pp. 747-755, 2014.

[9] S.-W. Lee, J. Kim, and Y. J. Lee, "Protecting signal integrity against atomic clock anomalies on board GNSS satellites," IEEE Transactions on Instrumentation and Measurement, vol. 60, no. 7, pp. 2738-2745, 2011.

[10] G. Signorile, "Analysis of a Kalman filter detector for atomic clock anomalies," in Proceedings of the 31st URSI General Assembly and Scientific Symposium (URSI GASS '14), pp. 1-4, Beijing, China, August 2014.

[11] L. Galleani and P Tavella, "The dynamic allan variance," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 56, no. 3, pp. 450-464, 2009.

[12] I. Sesia, L. Galleani, and P. Tavella, "Implementation of the dynamic allan variance for the galileo system test bed V2," in Proceedings of the IEEE International Frequency Control Symposium Joint with the 21st European Frequency and Time Forum, pp. 946-949, Geneva, Switzerland, June 2007

[13] A. Cenigliaro and S. Valloreia, "Analysis on GNSS space clocks performances," in Proceedings of the Joint European Frequency and Time Forum & International Frequency Control Symposium (EFTF/IFC '13), pp. 835-837, Prague, Czech Republic, July 2013.

[14] Z. Ji, Synchronization Technology and Application in Communication, Tsinghua University Press, Beijing, China, 2008.

[15] D. W. Allan, "Statistics of Atomic Frequency Standards," Proceedings of the IEEE, vol. 54, no. 2, pp. 221-230, 1966.

[16] D. Allan, "Time and frequency (time-domain) characterization, estimation, and prediction of precision clocks and oscillators," IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, vol. 34, no. 6, pp. 647-654, 1987

[17] J. E. Gray and D. W. Allan, "A method for estimating the frequency stability of an individual oscillator," in Proceedings of the 28th Annual Symposium on Frequency Conrol, pp. 243-246, May 1974.

[18] J. Groslambert, D. Fest, M. Olivier, and J. J. Gagnepain, "Characterization of frequency fluctuations by crosscorrelations and by using three or more oscillators," in Proceedings of the 35th Annual Frequency Control Symposium, Philadelphia, Pa, USA, May 1981.

[19] A. Premoli and P Tavella, "A revisited three-cornered hat method for estimating frequency standard instability," IEEE Transactions on Instrumentation and Measurement, vol. 42, no. 1, pp. 7-13, 1993.

[20] M. A. Weiss and C. Hackman, "Confidence on the three-point estimator of frequency drift," in Proceedings of the 24th Annual PTTI Meeing, pp. 451-460, December 1992.

[21] V. A. Logachev and G. P Pashev, "Estimation of linear frequency drift coefficient of frequency standards," in Proceedings of the 50th IEEE International Frequency Control Symposium, pp. 960-963, Honolulu, Hawaii, USA, June 1996.

[22] Q. Wei, "Estimations of frequency and its drift rate," IEEE Transactions on Instrumentation and Measurement, vol. 46, no. 1, pp. 79-82, 1997

[23] C. A. Greenhall, "Frequency-drift estimator and its removal from modified Allan variance," in Proceedings of the IEEE International Frequency Control Symposium, pp. 428-432, May 1997.

[24] W. J. Riley, Handbook of Frequency Stability Analysis, U.S. Department of Commerce, 2008.

[25] G. Hairong, "Determination of covariance matrix of kalman filter used for time prediction of atomic clocks of navigation satellites," Acta Geodaeicaet Cartographica Sinica, vol. 39, no. 2, pp. 146-150, 2010.

[26] L. Galleani, L. Sacerdote, P Tavella, and C. Zucca, "A mathematical model for the atomic clock error," Metrologia, vol. 40, no. 3, p. S257, 2003.

[27] L. Galleani and P Tavella, "Time and the Kalman filter," IEEE Control Systems Magazine, vol. 30, no. 2, pp. 44-65, 2010.

[28] C. A. Greenhall, "Forming stable timescales from the Jones-Tryon Kalman filter," Metrologia, vol. 40, no. 3, pp. S335-S341, 2003.

[29] C. Zucca and P Tavella, "The clock model and its relationship with the allan and related variances," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 52, no. 2, pp. 289-295, 2005.

Lei Feng (1,2) and Guotong Li (2,3)

(1) Shanghai Institute of Microsystem and Information Technology, Chinese Academy of Sciences, Shanghai 200050, China

(2) Shanghai Engineering Center for Microsatellite, Shanghai 201203, China

(3) Shanghai Tech University, Shanghai 201203, China

Correspondence should be addressed to Lei Feng; fenglei6036@126.com

Received 25 February 2016; Revised 14 June 2016; Accepted 29 June 2016

Academic Editor: Paolo Tortora

Caption: Figure 1: Global schematic diagram of Self-Monitoring Method.

Caption: Figure 2: Schematic diagram of Self-Monitoring Method based on PLL.

Caption: Figure 3: Response of Phase Detector to phase jumping during simulation.

Caption: Figure 4: Response of Phase Detector to frequency jumping during simulation.

Caption: Figure 5: Response of Phase Detector to signal loss during simulation.

Caption: Figure 6: Relative frequency deviation and its Allan deviation when anomaly occurs.

Caption: Figure 7: Phase and frequency jumping.

Caption: Figure 8: Comparison of two variances for phase jumping.

Caption: Figure 9: Comparison of two variances for frequency jumping.

Caption: Figure 10: Instantaneous deterioration of clock.

Caption: Figure 11: LS method in detecting frequency jumping.

Caption: Figure 12: Kalman filter method in detecting frequency jumping.

Caption: Figure 13: Comparison of different drift-rate estimators.

Caption: Figure 14: The drift-rate estimation of Kalman filter.

Caption: Figure 15: Prediction performance of Kalman filter for [DELTA][t.sub.1] and [DELTA][t.sub.12].

Caption: Figure 16: Drift-rate estimation of method 1.

Caption: Figure 17: Drift-rate estimation of method 2.

Table 1: Detection performance for phase jumping. Loop parameter Threshold PD PFA [epsilon] = [square root or 2]/2 8E - 9 1.0000 6.131e - 4 [[omega].sub.n] = 4 9E - 9 1.0000 1.193e - 4 10E - 9 0.9999 2.224e - 5 11E - 9 0.9978 2.412e - 6 12E - 9 0.9833 8.300e - 8 Loop parameter Delay (m) [epsilon] = [square root or 2]/2 1 [[omega].sub.n] = 4 1 1 1 1 Table 2: Detection performance for frequency jumping. Loop parameter Threshold PD PFA [xi] = [square root or 2]/2 8E -9 1.0000 6.131e - 4 [[omega].sub.n] = 4 9E -9 1.0000 1.193e - 4 10E -9 1.0000 2.224e - 5 11E -9 0.9933 2.412e - 6 12E -9 0.9873 8.300e - 8 Loop parameter Delay (m) [xi] = [square root or 2]/2 234 [[omega].sub.n] = 4 296 337 380 429 Table 3: Detection performance of phase jumping for different clock stability. Clock stability Resolution Threshold PD 3E - 12/[square root or [tau]] 1E - 8 11E - 9 0.9978 3E - 11/[square root or [tau]] 1E - 7 11E - 8 0.9972 3E - 10/[square root or [tau]] 1E - 6 11E - 7 0.9955 Clock stability PFA Delay (m) 3E - 12/[square root or [tau]] 2.412e - 6 1 3E - 11/[square root or [tau]] 6.290e - 7 1 3E - 10/[square root or [tau]] 4.370e - 7 1 Table 4: Detection performance of frequency jumping for different clock stability. Clock stability Resolution Threshold PD 3E - 12/[square root or [tau]] 4E - 11 11E - 9 0.9933 3E - 11/[square root or [tau]] 4E - 10 11E - 8 0.9947 3E - 10/[square root or [tau]] 4E - 9 11E - 7 0.9978 Clock stability PFA Delay (m) 3E - 12/[square root or [tau]] 2.412e - 6 380 3E - 11/[square root or [tau]] 6.290e - 7 384 3E - 10/[square root or [tau]] 4.370e - 7 381 Table 5: Statistical characteristic of of and [[sigma].sup.2.sub.i,m] and [v.sup.2.sub.i,m]. Statistical item DAVAR Modified DAVAR Minimum value 1.60e - 23 3.09e - 23 Maximum value 5.59e - 23 1.49e - 22 Mean value 2.89e - 23 8.14e - 23 Standard deviation 5.01e - 24 1.70e - 23 Table 6: Detection performance for phase jumping. Detection item DAVAR Modified DAVAR Resolution 12 12 Threshold 4.5E - 23 4.5E - 23 PD 1.0000 0.9978 PFA 1.61e - 5 3.00e - 7 Detection delay (points) 1 1 Table 7: Detection performance for phase jumping. Detection item DAVAR Modified DAVAR Resolution 12 4 Threshold 4.5E - 23 4.5E - 23 PD 0.9874 0.9927 PFA 1.61e - 5 3.00e - 7 Detection delay (points) 1 5 Table 8: Detection performance of phase jumping for different clock stability. [sigma]([[tau].sub.0]) Resolution Threshold PD 3E - 12 12 4.5E - 23 0.9978 3E - 11 12 4.5E - 21 0.9961 3E - 10 12 4.5E - 19 0.9977 [sigma]([[tau].sub.0]) PFA Delay (points) 3E - 12 3.00e - 7 1 3E - 11 5.00e - 7 1 3E - 10 6.00e - 7 1 Table 9: Detection performance of frequency jumping for different clock stability. [sigma]([[tau].sub.0]) Resolution Threshold PD 3E - 12 4 4.5E - 23 0.9927 3E - 11 4 4.5E - 21 0.9923 3E - 10 4 4.5E - 19 0.9934 [sigma]([[tau].sub.0]) PFA Delay (points) 3E - 12 3.00e - 7 5 3E - 11 5.00e - 7 5 3E - 10 6.00e - 7 5 Table 10: Detection performance of frequency jumping for different clock stability. Detection method Resolution PD PFA Delay PLL method 4.0E - 11 0.9933 2.41E - 6 0.4s DAVAR method 3.6E - 11 0.9874 1.61E - 5 1s Modified DAVAR method 1.2E - 11 0.9927 3.00E - 7 5s LS method 1.5E - 11 0.9944 2.53E - 5 1s Kalman filter method 1.8E - 11 0.9926 2.60E - 6 1s Table 11: Parameters of generated phase data. [[sigma].sub.0] [[sigma].sub.-1] [[sigma].sub.-2] 1E - 12 1E - 14 1E - 16 Drift rate/day Frequency deviation 1E - 12 1E - 12 Table 12: Statistical property of the estimation result of six estimators. Estimator a b c d E(z) 1.00E - 12 1.00E - 12 1.00E - 12 1.00E - 12 [sigma](z) 3.48E - 13 2.40E - 13 2.51E - 13 1.74E - 13 Estimator e f E(z) 1.00E - 12 9.97E - 13 [sigma](z) 1.58E - 13 3.71E - 14 Table 13: Parameter setting. Parameters [[sigma].sub.0] [[sigma].sub.-1] [[sigma].sub.-2] Clock 1 1e - 12 2e - 14 3e - 16 Clock 2 2e - 12 3e - 14 1e - 16 Clock 3 3e - 12 1e - 14 2e - 16 Drift Frequency Parameters rate/day deviation Clock 1 1e - 13 1e - 12 Clock 2 1e - 12 1e - 12 Clock 3 -2e - 11 1e - 12

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Feng, Lei; Li, Guotong |

Publication: | International Journal of Aerospace Engineering |

Date: | Jan 1, 2016 |

Words: | 7858 |

Previous Article: | A Bayesian classifier for X-ray pulsars recognition. |

Next Article: | Investigation of an autofocusing method for visible aerial cameras based on image processing techniques. |

Topics: |