Time scale algorithm based on Kalman filtering and fusion weight

Abstract

The invention discloses a time scale algorithm based on Kalman filtering and fusion weight. The time scale algorithm comprises the following steps: acquiring time deviation data and frequency deviation data of an atomic clock; calculating a time deviation optimal estimation value and a frequency deviation optimal estimation value of the atomic clocks in the clock group by adopting a Kalman filter, and calculating a time deviation optimal estimation value residual error and a frequency deviation optimal estimation value residual error of the atomic clocks in the clock group; calculating a fusion weight according to the time optimal estimation value residual error and the frequency optimal estimation value residual error of the atomic clocks in the clock group; and calculating a time scale through a time scale basic equation based on the fusion weight. According to the method, prediction values of time deviation and frequency deviation are improved through introduction of a Kalman filtering process, the noise influence of an atomic clock is reduced, and the real-time performance of a time scale is improved; and the time deviation residual error and the frequency deviation residual error are fused to calculate the weight of the atomic clock, so that the time scale with higher stability is established, and the medium and short term stability of the time scale is effectively improved.

Description

Technical Field

[0001] The invention belongs to the technical field of atomic clock data processing, and in particular relates to a time scale algorithm based on Kalman filtering and fusion weights. Background Art

[0002] With the current development of digitalization, informatization, and intelligence, time standards are becoming increasingly important as the basis for the normal operation of power, transportation, communications, networks, and other fields. Especially in the military field, high-precision time sources are an important support for most information equipment. Maintaining a stable and accurate time standard is very important for improving the time synchronization capability of the equipment system. Generally speaking, the key to building a high-performance time standard is how to establish and maintain time. In addition to being affected by the performance of hardware such as atomic clocks in the system, the time scale will also have different performance levels due to the different time scale algorithms selected. Therefore, studying the time scale algorithm is of great significance to improving the time scale performance.

[0003] The purpose of the time scale algorithm is to organically combine the data of each atomic clock in the clock-keeping group based on statistical methods, improve the algorithm stability, reduce uncertainty, and thus establish a time scale with better performance. The main time scale algorithms currently include the ALGOS algorithm, the AT1 algorithm, and the Kalman algorithm. Among them, the ALGOS algorithm measures the data of more than 400 atomic clocks in the world and calculates the weighted average to calculate the free atomic time (Echelle Atomique Libre, referred to as EAL), and then calibrates the EAL frequency to obtain the International Atomic Time (International Atomic Time, referred to as TAI). The AT1 algorithm uses the clock difference data of high-performance cesium clocks and hydrogen clocks with standard clocks to calculate the time scale. The AT1 algorithm only considers the change in frequency and has good real-time performance. The Kalman algorithm optimizes the time deviation between the reference clock and the ideal clock, superimposes all the single clocks in the clock group into the Kalman filter of the clock group, inputs the observed values ​​of the clock difference between all the clocks in the clock group and the reference clock, and outputs the optimal estimate of the phase and frequency of all the clocks.

[0004] The core idea of ​​the ALGOS algorithm and the AT1 algorithm is the weighted average method. By setting a weight for each atomic clock in the clock group, a more stable time scale is obtained after weighted averaging. However, the atomic clock noise cannot be eliminated or suppressed by weighted averaging, and the accuracy and reliability of the time scale cannot be further improved. In addition, the ALGOS algorithm does not have real-time performance and is difficult to meet the real-time requirements of the timekeeping system. The Kalman algorithm essentially uses a series of filter algorithms to mathematically model the various noises of the atomic clock to achieve a partial suppression effect. However, the existing Kalman algorithm will cause Kalman filtering divergence due to the non-objective nature of the system, and the improvement of short-term stability in the time scale still needs to be improved. Summary of the invention

[0005] The purpose of the present invention is to provide a time scale algorithm based on Kalman filtering and fusion weights with higher short- to medium-term stability.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A time-scaling algorithm based on Kalman filtering and fusion weights includes the following steps:

[0008] S1. Obtaining time deviation data and frequency deviation data of an atomic clock;

[0009] S2, using Kalman filter to calculate the optimal estimate of the time deviation and the optimal estimate of the frequency deviation of the atomic clocks in the clock group, and calculate the residual of the optimal estimate of the time deviation and the residual of the optimal estimate of the frequency deviation of the atomic clocks in the clock group;

[0010] S3, calculating the fusion weight according to the residual of the optimal estimated value of time and the residual of the optimal estimated value of frequency of the atomic clocks in the clock group;

[0011] S4. Based on the fusion weight, the time scale is calculated through the time scale basic equation.

[0012] Furthermore, the basic equation of the time scale is:

[0013] h in the formula i (t) is the clock time of the ith atomic clock, h i '(t) is the predicted clock time of the i-th atomic clock, p i is the fusion weight of the i-th atomic clock in the clock group;

[0014] Fusion weight In the formula is the weight calculated based on the residual of the optimal estimate of the time deviation of the i-th atomic clock, is the weight calculated based on the residual of the optimal estimate of the frequency deviation of the i-th atomic clock, m and n are weight coefficients, and m+n=1.

[0015] Furthermore,

[0016]

[0017] ε in the formula i is the residual of the optimal estimate of the time deviation of the i-th atomic clock, ε j is the residual of the optimal estimate of the time bias of the jth atomic clock, λ i is the residual of the optimal estimate of the frequency deviation of the i-th atomic clock, λ j is the residual of the optimal estimate of the frequency deviation of the jth atomic clock, i,j=1,2,3,…,N, N is the number of atomic clocks in the clock group.

[0018] Furthermore, the values ​​of m and n are determined according to the following formula:

[0019]

[0020] In the formula is the weight calculated based on the residual of the optimal estimate of the time bias of the jth atomic clock, is the weight calculated based on the residual of the optimal estimate of the frequency deviation of the jth atomic clock, is the optimal estimate of the i+1th frequency deviation of the jth atomic clock, is the optimal estimate of the i-th frequency deviation of the j-th atomic clock;

[0021] The fusion weight p is further determined according to the following conditions: i :

[0022]

[0023] Furthermore, the residual of the optimal estimate of the time deviation of the i-th atomic clock is

[0024] The residual of the optimal estimate of the frequency bias of the i-th atomic clock

[0025] are the time deviation and frequency deviation of the i-th atomic clock at the j-th moment, are the expected value of time deviation and frequency deviation of the i-th atomic clock at the j-th moment respectively.

[0026] Furthermore, the Kalman filter is:

[0027]

[0028] Where X(k / k-1) represents the predicted value at time k, X(k-1 / k-1) is the optimal estimate at time k-1, P(k / k-1) is the predicted value of the error variance matrix at time k, P(k-1 / k-1) is the optimal estimate of the error variance matrix at time k-1, X(k / k) is the optimal estimate at time k, P(k / k) is the optimal estimate of the error variance matrix at time k, R is the observation noise matrix, K is the Kalman gain, I is the unit matrix, Q i is the system noise matrix of the i-th atomic clock, H is the measurement matrix, Z(k) is the observed value of the k-th measurement of the clock group, T is the measurement time interval.

[0029] Furthermore, Z(k) = HX(k) + V(k), where X(k) is the dynamic model of the clock group and V(k) is the observation noise;

[0030] The dynamic model of the clock group is X(k)=AX(k-1)+W(k-1), where X(k-1)=[x1(k-1)y1(k-1)x2(k-1)y2(k-1)...x i (k-1)y i (k-1)] T , W(k-1)=[ξ1(k-1)η1(k-1)ξ2(k-1)η2(k-1)...ξ i (k-1)η i (k-1)] T , where x1(k-1), x2(k-1), x i (k-1) are the time deviations of the 1st, 2nd, and i-th atomic clocks during the k-1th measurement, y1(k-1), y2(k-1), y i where (k-1) are the frequency deviations of the 1st, 2nd, and i-th atomic clocks during the k-1th measurement, ξ1(k-1), ξ2(k-1), ξ i (k-1) are the frequency white noises of the 1st, 2nd, and i-th atomic clocks during the k-1th measurement, η1(k-1), η2(k-1), η i (k-1) are the frequency random walk noises of the 1st, 2nd and i-th atomic clocks during the k-1th measurement.

[0031] Furthermore, an atomic clock output model and a noise model are established, and the atomic clock output model and the noise model are used to simulate and generate the time deviation data and the frequency deviation data of the atomic clock.

[0032] Furthermore, the atomic clock output model is: V(t)=(V0+ε(t))sin(2πv0t+φ(t)), where V0 represents the nominal voltage amplitude, ε(t) represents the random amplitude deviation at time t, v0 represents the nominal frequency, and φ(t) represents the phase deviation at time t relative to 2πv0t.

[0033] Furthermore, the noise model includes frequency random walk noise, frequency flicker noise, frequency white noise and phase white noise;

[0034] The frequency white noise is Gaussian random white noise; the frequency white noise is integrated to obtain the frequency random walk noise; the frequency white noise is differentiated to obtain the phase white noise; the frequency flicker noise is simulated by using a uniformly distributed random variable or a Markov random process or a noise modeling filter transfer function.

[0035] It can be seen from the above technical solutions that the present invention proposes a new fusion weight calculation method to address the fact that the traditional weighted average scaling algorithm does not pay enough attention to the noise of the clock group and the timekeeping system's demand for real-time performance. The Kalman filter is used to calculate the optimal estimate of the time deviation and the optimal estimate of the frequency deviation of all the clocks in the clock group. The Kalman filtering process is introduced to improve the predicted values ​​of the time deviation and the frequency deviation, reduce the noise impact of the atomic clock, and improve the real-time performance of the time scale. Moreover, the time deviation residual and the frequency deviation residual are fused to calculate the atomic clock weight, thereby establishing a more stable time scale, which effectively improves the medium- and short-term stability of the time scale. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] In order to more clearly illustrate the embodiments of the present invention, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0037] Figure 1 is a flow chart of the method of the present invention;

[0038] Figure 2 It is the calculation relationship diagram of WFM, WPM and RWFM;

[0039] Figure 3a to Figure 3d They are WFM, WPM, RWFM and FFM simulated frequency data;

[0040] Figures 4a to 4d Allen deviation diagrams for simulated WFM, WPM, RWFM, and FFM, respectively;

[0041] Figure 5The clock error data of 5 simulated clocks and the comparison chart of 2 time scales;

[0042] Figure 6 The clock error diagrams are of four time scales;

[0043] Figure 7 Allan variance diagrams for four time scales;

[0044] Figure 8 Allan variance diagrams for five simulated clocks and two time scales.

[0045] The specific implementation modes of the present invention are further described in detail below with reference to the accompanying drawings. DETAILED DESCRIPTION

[0046] The present invention is described in detail below in conjunction with the accompanying drawings. When describing the embodiments of the present invention in detail, for the convenience of explanation, the drawings representing the device structure will not be partially enlarged according to the general proportion, and the schematic diagram is only an example, which should not limit the scope of protection of the present invention. It should be noted that the drawings are simplified in form and all use non-precise proportions, which are only used to facilitate and clearly assist in explaining the purpose of the embodiments of the present invention. At the same time, in the description of this application, the terms "first", "second", etc. are only used to distinguish the description, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated; the terms "positive", "negative", "bottom", "upper", "lower", etc. indicate the orientation or position relationship based on the orientation or position relationship shown in the drawings, which is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation on the present invention.

[0047] The time scale algorithm can comprehensively process several atomic clocks with different performances to form a "virtual clock" with better performance than a single atomic clock in the original clock group. This "virtual clock" has better stability and reliability. In view of the problem that the traditional ALGOS algorithm and AT1 algorithm based on the weighted average method lack attention to noise, and the timekeeping system's own demand for real-time performance, the present invention is based on Kalman filtering and the weight calculation of the fusion time deviation and frequency deviation residuals to construct a new time scale algorithm. After the method of the present invention calculates the optimal estimate of the time deviation and frequency deviation of the atomic clock based on the Kalman filter, the time deviation residual and the frequency deviation residual are fused to calculate the atomic clock weight, and finally the time scale is obtained by weighted averaging using the basic equation of the time scale.

[0048] Combine the following Figure 1 , the method of the present invention is described. For a clock group having N atomic clocks, the steps of the method of the present invention are as follows:

[0049] S1. Obtaining time deviation data and frequency deviation data of an atomic clock;

[0050] S2, using a Kalman filter to calculate the optimal estimate of the time deviation and the optimal estimate of the frequency deviation of the atomic clock, and calculate the residual of the optimal estimate of the time deviation and the residual of the optimal estimate of the frequency deviation;

[0051] S3, calculating the fusion weight according to the residual of the optimal time estimate and the residual of the optimal frequency estimate;

[0052] S4. Based on the fusion weight, the time scale is calculated through the time scale basic equation.

[0053] In step S1, for a real clock group, the time deviation data and frequency deviation data of the atomic clock can be obtained by measurement, or by simulation. In this embodiment, the time deviation data and frequency deviation data of the atomic clock are obtained by simulation. When obtaining data by simulation, an atomic clock output model and a noise model are established, and the atomic clock output model and the noise model are used to simulate and generate the time deviation data and frequency deviation data of the atomic clock.

[0054] The atomic clock is a stable frequency source. Its output instantaneous voltage signal can be expressed using a sine function. The output model of the atomic clock is:

[0055] V(t)=(V0+ε(t))sin(2πv0t+φ(t)),

[0056] Where V0 represents the nominal voltage amplitude, ε(t) represents the random amplitude deviation at time t, v0 represents the nominal frequency, and φ(t) represents the phase deviation at time t relative to 2πv0t.

[0057] Oscillation period fluctuations and phase fluctuations will affect the stability of the atomic clock's time measurement. The phase deviation x(t), frequency deviation y(t), and frequency drift deviation ω(t) between the reference clock and the ideal clock are:

[0058]

[0059] Where τ is the time interval, ξ is the random noise error in the phase deviation, η is the random noise error in the frequency deviation, and ζ is the random noise error in the frequency drift deviation.

[0060] The deviation (phase deviation, frequency deviation, frequency drift deviation) between the reference clock and the ideal clock consists of two parts: one is the deviation caused by the initial deviation (initial phase deviation, initial frequency deviation, initial frequency drift deviation), and the other is the deviation caused by noise. Therefore, when obtaining the time deviation data and frequency deviation data of the atomic clock, a noise model must also be established.

[0061] Usually, a superposition model based on five kinds of noises, namely, frequency random walk noise (hereinafter referred to as RWFM), frequency flicker noise (hereinafter referred to as FFM), frequency white noise (hereinafter referred to as WFM), phase flicker noise (hereinafter referred to as FPM) and phase white noise (hereinafter referred to as WPM) is used as a noise model, and a power spectrum density function is used to describe the statistical relationship between the five kinds of noises. It is found through research that FPM has little effect on the time deviation data and frequency deviation data of the simulated atomic clock, so the noise model established by the present invention does not include FPM.

[0062] WFM is a Gaussian random noise with zero mean and unit standard deviation. Its power spectral density S(f) = f 0 , f represents the Fourier frequency. Figure 2 As shown, integrating and differentiating WFM can obtain RWFM and WPM respectively, so the noise models of WFM, RWFM and WPM are as follows:

[0063]

[0064] f in the formula G represents Gaussian random white noise, represents colored noise with power law i, iη2,0,2.

[0065] For FFM, simulation can be performed using methods such as uniformly distributed random variables, Markov random processes, and noise modeling filter transfer functions. This embodiment uses the linear superposition of four independent Markov frequency processes to simulate FFM, and the FFM noise model is:

[0066] y i (t) represents an independent Markov frequency process, i = 1, 2, 3, 4;

[0067] The Markov frequency process y(t) is:

[0068] Where x(t) represents the time deviation (phase deviation) corresponding to the Markov frequency process, and w(t) is the double-sided spectral density function The white noise process, R m is a constant.

[0069] right Solving in discrete time with a time interval of τ0, we obtain the difference equation:

[0070] In the formula, ω x (t),ω y (t) is a white noise process with a mean of 0, and (ωx (t)ω y (t)) T The covariance matrix satisfies Q M ;

[0071]

[0072] According to the above noise model simulation, four types of noise are generated: frequency white noise WFM, phase white noise WPM, frequency random walk noise RWFM and frequency flicker noise FFM.

[0073] In step S2, when using the Kalman filter to calculate the optimal estimate of the time deviation and the optimal estimate of the frequency deviation of the atomic clock, a clock group observation equation is constructed. The clock group observation equation of this embodiment is:

[0074] Z(k)=HX(k)+V(k),

[0075] Where Z(k) is the observed value of the kth measurement of the clock group, H is the measurement matrix, X(k) is the dynamic model of the clock group, and V(k) is the observation noise.

[0076] Measurement Matrix N is the number of atomic clocks in the clock group.

[0077] The state space model of the i-th atomic clock in the clock group is:

[0078] The x in the formula i (k) is the time deviation (phase deviation) of the i-th atomic clock during the k-th measurement, y i (k) is the frequency deviation of the k-th measurement of the i-th atomic clock. The time deviation of the present invention refers to the difference between the time of the atomic clock and the time of the ideal time scale. The time deviation is actually the difference between the phase of the atomic clock and the phase of the ideal time scale. The frequency deviation refers to the difference between the frequency of the atomic clock and the frequency of the ideal time scale. By analogy, x i (k-1) is the time deviation of the i-th atomic clock at the k-1th measurement, y i (k-1) is the frequency deviation of the i-th atomic clock measured for the k-1th time, T is the measurement time interval, ξ i (k-1) is the variance of the i-th atomic clock at the k-1th measurement, which is E i The frequency of white noise, η i (k-1) is the variance of the i-th atomic clock at the k-1th measurement, which is H i The frequency random walk noise of the i-th atomic clock is

[0079] Based on the state space model of each atomic clock, the dynamic model of the clock group is obtained:

[0080] X(k)=AX(k-1)+W(k-1),

[0081] Where X(k-1)=[x1(k-1)y1(k-1)x2(k-1)y2(k-1)...x i (k-1)y i (k-1)] T , W(k-1)=[ξ1(k-1)η1(k-1)ξ2(k-1)η2(k-1)...ξ i (k-1)η i (k-1)] T ,

[0082]

[0083] The Kalman filter is used to calculate the optimal estimated value of the time deviation and the optimal estimated value of the frequency deviation of each atomic clock i in the clock group. The Kalman filter of the present invention is:

[0084]

[0085] In the formula, X(k / k-1) represents the predicted value at time k (including the predicted value of time deviation and the predicted value of frequency deviation). The predicted value is used for iterative update in the Kalman filter to calculate the optimal estimate. X(k-1 / k-1) is the optimal estimate at time k-1, X(k / k) is the optimal estimate at time k, P(k / k-1) is the predicted value of the error variance matrix at time k, P(k-1 / k-1) is the optimal estimate of the error variance matrix at time k-1, P(k / k) is the optimal estimate of the error variance matrix at time k, R is the observation noise matrix, K is the Kalman gain, and I is the unit matrix. The essence of the error variance matrix is ​​the covariance matrix used for iterative calculation.

[0086] After obtaining the optimal estimate of the time deviation and frequency deviation of each atomic clock through the Kalman filter, the residual ε of the optimal estimate of the time deviation of each atomic clock is calculated. i and the optimal estimate of the frequency deviation residual λ i .

[0087] The residual of the optimal estimate of the time bias of the i-th atomic clock

[0088] The residual of the optimal estimate of the frequency bias of the i-th atomic clock

[0089] are the time deviation and frequency deviation of the i-th atomic clock at the j-th moment, The specific values ​​​​substituted into are the optimal estimates of the time deviation and frequency deviation of the atomic clock after passing through the Kalman filter. are the expected time deviation and frequency deviation of the i-th atomic clock at the j-th moment, respectively. Since the ideal clock has both 0 time deviation and frequency deviation, The value of is 0.

[0090] In step S3, the fusion weight p is calculated according to the following formula: i :

[0091]

[0092] p in the formula i is the fusion weight of the i-th atomic clock, is the weight calculated based on the residual of the optimal estimate of the time deviation of the i-th atomic clock, is the weight calculated based on the residual of the optimal estimate of the frequency deviation of the i-th atomic clock, m and n are weight coefficients, and m and n satisfy m+n=1;

[0093] ε in the formula i is the residual of the optimal estimate of the time deviation of the i-th atomic clock, ε j is the residual of the optimal estimate of the time bias of the jth atomic clock, i,j=1,2,3,…,N;

[0094] In the formula is the residual square value of the optimal estimate of the frequency deviation of the i-th atomic clock, is the residual sum of squares of the optimal estimate of the clock group frequency deviation, i,j=1,2,3,…,N.

[0095] The present invention uses the ratio of the square value of the residual of the optimal estimate of the time deviation of each atomic clock to the sum of the square values ​​of the residual of the optimal estimate of the time deviation of the clock group as the weight of the clock. At the same time, the residual of the optimal estimate of the frequency deviation is introduced when calculating the weight of the atomic clock, which can optimize the stability of the clock group time.

[0096] The present invention uses Allan variance To express the stability of the time scale. For the clock group, in order to improve the stability of the clock group time scale and minimize the Allan variance of the clock group time scale frequency deviation data, the values ​​of m and n are determined according to the following formula:

[0097]

[0098] In the formula is the weight calculated based on the residual of the optimal estimate of the time bias of the jth atomic clock, is the weight calculated based on the residual of the optimal estimate of the frequency deviation of the jth atomic clock, is the optimal estimate of the frequency deviation of the jth atomic clock (i+1) calculated by the Kalman filter, represents the optimal estimate of the frequency deviation of the jth atomic clock calculated by the Kalman filter.

[0099] Allen Variance t is the sampling time, J is the number of sampling points, and They respectively represent the i-th and i+1-th frequency deviation data of the time scale obtained after fusion weight calculation, i=1,2,3,…,J-1. Where p j is the fusion weight of the jth atomic clock, is the ith original frequency deviation data of the jth atomic clock, It is the original frequency deviation data that has not been processed by the Kalman filter. It can generally be measured by measuring instruments such as frequency counters, or simulated data can be obtained through simulation.

[0100] After calculating the values ​​of m and n, the fusion weight is determined according to the following conditions to adapt to clock groups with different performances:

[0101]

[0102] Right now:

[0103] (1) When the calculated value of m is in the interval [0,1], the calculated value of n must also be in the interval [0,1]. In this case, the calculated values ​​of m and n are directly substituted into Calculate the fusion weight;

[0104] (2) When the calculation result of m is in the interval (-∞,0), it is considered that the calculation of the fusion weight introduces The fusion effect will be reduced. Select m=0, n=1 and substitute Calculate the fusion weight;

[0105] (3) When the calculation result of m is in the interval (1, +∞), it is considered that the calculation of the fusion weight introduces The fusion effect will be reduced, so m=1,n=0 is selected and substituted into Calculate the fusion weight.

[0106] In step S4, the basic time scale equation of the present invention is:

[0107] h in the formula i (t) is the clock time of the ith atomic clock, hi '(t) is the predicted clock time of the i-th atomic clock.

[0108] The present invention uses the clock face time values ​​of two time scales (two atomic clocks) to calculate the time scale TA(t). For the i-th atomic clock and the j-th atomic clock in the clock group, the time deviation of the two atomic clocks at time t is x ij (t) = h i (t)-h j (t), i,j=1,2,3,4,...,N,h i (t), h j (t) are the clock times of clock i and clock j at time t. The time deviation x between the i-th atomic clock and the j-th atomic clock ij (t) can be measured or calculated to obtain the predicted clock time h of the i-th atomic clock. i '(t) can be predicted by Kalman filter.

[0109] In order to verify the effect of the method of the present invention, the method of this embodiment is compared with the traditional method of calculating weights based on the residual of the optimal estimate of the time deviation. The difference between the traditional method of calculating weights based on the residual of the optimal estimate of the time deviation and the method of the present invention is that the traditional method only calculates weights by time residuals, that is, weights The method of the present invention also introduces frequency deviation residual to achieve a fusion effect.

[0110] During simulation, the time deviation data and frequency deviation data of the atomic clock are obtained by simulation. The time deviation data and frequency deviation data of the atomic clock are simulated by the atomic clock output model and noise model. According to the principle of power-law spectrum noise, the simulated frequency data of 5 hydrogen clocks with different performances are generated, which contain 4 types of noise: WFM, WPM, RWFM, and FFM. Each atomic clock has 10,000 data points, and the time interval τ between adjacent data points is set to 1s. Figure 3a to Figure 3d The simulation data of clock 1 is shown. The Allan deviation diagram of the corresponding noise is calculated based on the simulation data, as shown in FIG. Figures 4a to 4d The slope of the linear fitting curve in the Allan deviation diagram of the simulated atomic clock is shown in Table 1.

[0111] Table 1 Allen deviation slope of simulated random noise

[0112]

[0113] The simulated frequency data is integrated to obtain the time deviation data and frequency deviation data of the simulated clock. According to the Kalman filter, the time deviation data and frequency deviation data of the five simulated clocks are substituted to calculate the optimal estimated values ​​of the time deviation and the optimal estimated values ​​of the frequency deviation of the five simulated clocks.

[0114] The traditional method is to calculate the residual based on the optimal estimate of the time deviation of each atomic clock, calculate the weight of each atomic clock, and then calculate the time scale TA1.

[0115] The method of the present invention is to calculate the time deviation optimal estimate value residual and the frequency deviation optimal estimate value residual respectively according to the time deviation optimal estimate value and the frequency deviation optimal estimate value residual of each atomic clock, and then calculate the fusion weight according to the time deviation optimal estimate value residual and the frequency deviation optimal estimate value residual, when calculating the fusion weight, select appropriate m and n according to the performance of the simulated clock in the clock group, determine the proportion of the time deviation and the frequency deviation in the fusion weight by the values ​​of m and n, and finally calculate the time scale TA2 according to the fusion weight.

[0116] The clock errors of the five simulated clocks and the clock error data of the time scales TA1 and TA2 simulated by the method of the present invention are compared and analyzed, and the root mean square value (RMS) is calculated as follows: Figure 5 And as shown in Table 2.

[0117] Table 2 Accuracy of five simulated clocks and two time scales

[0118]

[0119] from Figure 5 It can be seen from the results in Table 3 that the absolute value of the clock error of TA2 is significantly different from that of TA1, among which the maximum clock error of TA2 is 32.5746ns, while the maximum clock error of TA1 is 0.6464ns. It is found from the clock error data of each atomic clock in the clock group that the clock error of clock 5 changes more with time than the other four atomic clocks. However, by calculating the Allen variance of each atomic clock in the clock group, it is found that clock 5 has better short-term stability. If clock 5 is simply eliminated when calculating the time scale, the short-term stability of the time scale calculated at this time will be lower than that of TA2. This shows that although the short-term accuracy can be improved by eliminating atomic clocks with large clock errors, the error accumulation caused by high stability will be smaller as time accumulates. If an atomic clock with a large clock error deviation but good stability appears in the clock group, the method of the present invention uses fusion weights to give play to its stability advantages, essentially improving the time stability of the clock group more, which is beneficial to the long-term accuracy improvement. Figure 6 , Figure 7 As shown, after clock 5 is removed from the clock group, TA1' and TA2' are calculated using the same method as that for calculating time scales TA1 and TA2, and compared with TA1 and TA2.

[0120] from Figure 6 and Figure 7It can be seen that when calculating the time scale clock error, by comparing the clock error data of TA2 and TA2', it is found that removing clock 5 with a large error can effectively reduce the clock error of TA2', while there is no obvious difference between TA1' and TA1. This is because only the time deviation is used to determine the weight when calculating TA1' and TA1, and the weight set for clock 5 with a large time deviation is small. Therefore, there is almost no difference in the calculation of the clock error of time TA1' and TA1 with or without clock 5. However, the fusion weight used in the TA2 and TA2' time scales contains frequency deviation factors. At the same time, when calculating the fusion weight, considering that clock 5 has better short-term frequency stability, the weight value set for clock 5 through the fusion weight is greater than the weight value calculated by time deviation. Therefore, the accuracy of the clock error data of TA2 and TA2' is lower than that of TA1' and TA1, and after removing clock 5, the clock error performance of TA2' is significantly higher than that of TA2. The root mean square errors of the clock errors of TA1' and TA2' are 0.41617ns and 2.5935ns respectively. The above results show that if there is an atomic clock with a large clock error but good stability in the clock group, the fusion weight method of the present invention can be used to take advantage of its good stability to improve long-term accuracy, rather than simply removing it from the clock group.

[0121] When calculating the Allan variance, if clock 5 is removed, it is found that the frequency stability of TA1' is lower than that of TA1, and the frequency stability of TA2' is also lower than that of TA2. In order to achieve better short-term stability of the time scale, clock 5 does not need to be removed from the clock group. This also verifies the rationality of the method of determining the fusion weight value by calculating the minimum Allan variance of the frequency deviation of the time scale.

[0122] Next, the stability of the time scale is analyzed by Allan variance. The corresponding Allan variance is calculated based on the frequency data of five simulated clocks and two time scales with weights determined by different methods, as shown in Figure 8 As shown. Figure 8It can be seen that the TA2 time scale is more stable in the short term. Relatively speaking, the TA1 time scale is less stable than TA2 because the clock weight is determined only by the time deviation, and even lower than the stability of a single atomic clock in the clock group. This is because the TA1 calculated by determining the atomic clock weight only by the time deviation ignores the factor of frequency stability, and the weight is not set reasonably. The Allan variance is mainly calculated using frequency deviation data. If the weight is determined only by the time deviation, the factor of frequency deviation will be ignored, resulting in inaccurate calculated weights, and even making the time scale stability obtained by weighted average lower than that of a single atomic clock in the original clock group, which sometimes causes the Allan variance of the time scale TA1 to be greater than the Allan variance of the atomic clock in the original clock group. Time scale TA2 also takes frequency deviation as a factor in the calculation of atomic clock weights. The calculated Allan variance is smaller than that of TA1. Therefore, the overall stability of TA2 is better than that of TA1. Similarly, by comparing the Allan variance with that of a single atomic clock in the original clock group, it can be found that within the same time interval, the Allan variance of TA2 is smaller than the corresponding Allan variance of a single atomic clock in the original clock group. Therefore, the stability of TA2 is better than the stability of a single atomic clock in the original clock group.

[0123] Through the above simulation verification, the method of the present invention adopts fusion weights, and the time scale of calculation can effectively improve the short-term stability of the time scale and has real-time performance.

[0124] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not used to limit the present invention. Any technician familiar with this profession can make some changes or modify the technical contents disclosed above into equivalent embodiments without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A time scale algorithm based on Kalman filtering and fusion weights, characterized in that: The following steps are involved: S1. Obtaining time deviation data and frequency deviation data of an atomic clock; S2, using Kalman filter to calculate the optimal estimate of the time deviation and the optimal estimate of the frequency deviation of the atomic clocks in the clock group, and calculate the residual of the optimal estimate of the time deviation and the residual of the optimal estimate of the frequency deviation of the atomic clocks in the clock group; S3, calculating the fusion weight according to the residual of the optimal estimated value of time and the residual of the optimal estimated value of frequency of the atomic clocks in the clock group; S4. Based on the fusion weight, the time scale is calculated through the time scale basic equation.

2. The time scale algorithm of claim 1, wherein: The basic equation for the time scale is: h in the formula i (t) is the clock time of the ith atomic clock, h i '(t) is the predicted clock time of the i-th atomic clock, p i is the fusion weight of the i-th atomic clock in the clock group; Fusion weight In the formula is the weight calculated based on the residual of the optimal estimate of the time deviation of the i-th atomic clock, is the weight calculated based on the residual of the optimal estimate of the frequency deviation of the i-th atomic clock, m and n are weight coefficients, and m+n=1.

3. The time scale algorithm as claimed in claim 2, characterized in that: ε in the formula i is the residual of the optimal estimate of the time deviation of the i-th atomic clock, ε j is the residual of the optimal estimate of the time bias of the jth atomic clock, λ i is the residual of the optimal estimate of the frequency deviation of the i-th atomic clock, λ j is the residual of the optimal estimate of the frequency deviation of the jth atomic clock, i,j=1,2,3,…,N, N is the number of atomic clocks in the clock group.

4. The time scale algorithm of claim 2, wherein: The values ​​of m and n are determined according to the following formula: In the formula is the weight calculated based on the residual of the optimal estimate of the time bias of the jth atomic clock, is the weight calculated based on the residual of the optimal estimate of the frequency deviation of the jth atomic clock, is the optimal estimate of the i+1th frequency deviation of the jth atomic clock, is the optimal estimate of the ith frequency deviation of the jth atomic clock, and J is the number of sampling points; The fusion weight p is further determined according to the following conditions: i :

5. The time scale algorithm of claim 1, wherein: The residual of the optimal estimate of the time bias of the i-th atomic clock The residual of the optimal estimate of the frequency bias of the i-th atomic clock are the time deviation and frequency deviation of the i-th atomic clock at the j-th moment, are the expected value of time deviation and frequency deviation of the i-th atomic clock at the j-th moment respectively.

6. The time scale algorithm of claim 1, wherein: The Kalman filter is: Where X(k / k-1) represents the predicted value at time k, X(k-1 / k-1) is the optimal estimate at time k-1, P(k / k-1) is the predicted value of the error variance matrix at time k, P(k-1 / k-1) is the optimal estimate of the error variance matrix at time k-1, X(k / k) is the optimal estimate at time k, P(k / k) is the optimal estimate of the error variance matrix at time k, R is the observation noise matrix, K is the Kalman gain, I is the unit matrix, Q i is the system noise matrix of the i-th atomic clock, H is the measurement matrix, Z(k) is the observed value of the k-th measurement of the clock group, T is the measurement time interval.

7. The time scale algorithm of claim 6, wherein: Z(k)=HX(k)+V(k), X(k) is the dynamic model of the clock group, V(k) is the observation noise; The dynamic model of the clock group is X(k)=AX(k-1)+W(k-1), where X(k-1)=[x1(k-1)y1(k-1)x2(k-1)y2(k-1)...x i (k-1)y i (k-1)] T , W(k-1)=[ξ1(k-1)η1(k-1)ξ2(k-1)η2(k-1)...ξ i (k-1)th i (k-1)] T , among them, x1(k-1), x2(k-1), x i (k-1) are the time deviations of the 1st, 2nd, and i-th atomic clocks during the k-1th measurement, y1(k-1), y2(k-1), y i (k-1) are the frequency deviations of the 1st, 2nd, and i-th atomic clocks during the k-1th measurement, ξ1(k-1), ξ2(k-1), ξ i (k-1) are the frequency white noises of the 1st, 2nd, and i-th atomic clocks during the k-1th measurement, η1(k-1), η2(k-1), η i (k-1) are the frequency random walk noises of the 1st, 2nd and i-th atomic clocks during the k-1th measurement.

8. The time scale algorithm of claim 1, wherein: An atomic clock output model and a noise model are established, and the time deviation data and frequency deviation data of the atomic clock are simulated by the atomic clock output model and the noise model.

9. The time scale algorithm of claim 8, wherein: The atomic clock output model is: V(t)=(V0+ε(t))sin(2πv0t+φ(t)), where V0 represents the nominal voltage amplitude, ε(t) represents the random amplitude deviation at time t, v0 represents the nominal frequency, and φ(t) represents the phase deviation at time t relative to 2πv0t.

10. The time scale algorithm of claim 8, wherein: The noise model includes frequency random walk noise, frequency flicker noise, frequency white noise and phase white noise; The frequency white noise is Gaussian random white noise; the frequency white noise is integrated to obtain the frequency random walk noise; the frequency white noise is differentiated to obtain the phase white noise; the frequency flicker noise is simulated by using a uniformly distributed random variable or a Markov random process or a noise modeling filter transfer function.

Images