Kalman filtering time scale establishment method considering abnormal data

Abstract

The application discloses a Kalman filter time scale establishment method considering abnormal data, and comprises the following steps: acquiring a time difference measurement value sequence of an atomic clock group, processing the time difference measurement value sequence to obtain atomic clock state variable solution information of all epochs in the time difference measurement value sequence, and calibrating an arbitrary single time-keeping clock by using the atomic clock state variable solution information to obtain a calibrated time scale. Compared with an existing time scale algorithm, the algorithm can directly eliminate the influence of atomic clock time difference observation data interruption and abnormal conditions on the time scale, and the generated time scale is continuous and stable.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of time frequency technology, and in particular to a Kalman filter time scale establishment method considering abnormal data. BACKGROUND

[0002] Time plays a very important role in the economy, social life, military and other aspects of a country. The stability of time is related to the safety and stability of the country and society. As a major national infrastructure, the construction of time frequency system is of great significance. Satellite navigation, communication interconnection, joint operations, finance, power and other fields have higher and higher requirements for the accuracy and reliability of time frequency system. The measurement of time requires accurate time reference, and the appearance of atomic clock provides a high-accuracy and high-stability frequency source for the establishment of time reference. Each atomic clock can be theoretically regarded as an independent time reference, but a single atomic clock may fail, have sudden changes in time and rate, and other faults, and in fact cannot be used as a time reference. The time scale algorithm integrates multiple atomic clocks to reduce the noise of a single atomic clock and establish a more stable and reliable paper time.

[0003] The time scale algorithm usually includes weighted average algorithm (ALGOS algorithm and AT1 algorithm) and Kalman filter algorithm. The weighted average algorithm is essentially a weighted average according to the noise characteristics of the atomic clocks in the clock group at a certain smoothing time, without considering all noise characteristics at all smoothing times, so this method can only make the output time scale stability at the specified smoothing time point optimal. The Kalman filter time scale algorithm can smooth the atomic clock noise in real time by reasonably modeling the atomic clock noise model and accurately estimating the noise variance, and the established time scale has the characteristics of high short-term and long-term stability. The standard Kalman filter time scale algorithm needs to be linear optimal when the dynamic model noise and measurement variance noise conform to the Gaussian distribution. However, there are abnormal situations such as interruption, abnormal data, noise change, and inaccurate dynamic model of atomic clock data. At this time, the performance of the Kalman filter time scale algorithm is reduced. The existing research mainly focuses on the inaccurate model construction, and mainly increases the adaptive factor on the state estimation covariance matrix to limit the state estimation bias caused by model error. In the case of interruption and abnormality of atomic clock time difference observation data, there is less research on how to reduce the influence of abnormal time difference observation data on the performance of Kalman filter time scale. Therefore, the problem of the existing technology is that the existing Kalman filter time scale calculation method focuses on limiting the influence of model error on state estimation bias, and lacks the influence of reducing or even eliminating abnormal time difference observation data on the time scale calculation method. SUMMARY

[0004] The application discloses a Kalman filter time scale establishment method considering abnormal data.

[0005] To achieve the above object, the application discloses a Kalman filter time scale establishment method considering abnormal data, which comprises the following steps:

[0006] S1, two atomic clocks are selected from an atomic clock group as master clocks of a time scale algorithm, other atomic clocks are time-keeping clocks, a time difference measurement value sequence between the time-keeping clocks and the master clocks is obtained, and a sampling interval is set as τ; the atomic clock group comprises n atomic clocks, the n is a positive integer not less than 4; one of the two master clocks is a main master clock, and the other is a backup master clock;

[0007] S2, the time difference measurement value sequence is processed to obtain atomic clock state variable calculation information of all epochs in the time difference measurement value sequence; the state variable calculation information represents a single atomic clock and a virtual paper time difference;

[0008] S3, any single time-keeping clock is calibrated by using the atomic clock state variable calculation information to obtain a calibrated time scale.

[0009] As an optional implementation, the processing of the time difference measurement value sequence to obtain the atomic clock state variable calculation information of all epochs in the time difference measurement value sequence comprises the following steps:

[0010] S21, initial state information of the atomic clock group is constructed; the initial state information of the atomic clock group comprises initial state variables X k , an initial state variable covariance matrix P k and a covariance matrix Q0 of an initial driving noise η k of the n atomic clocks;

[0011] S22, an epoch corresponding to any time difference measurement value in the time difference measurement value sequence is selected as a calculation epoch;

[0012] S23, abnormal data identification processing is performed on the time difference measurement value at the calculation epoch moment in the time difference measurement value sequence to obtain abnormal atomic clock identification information at the calculation epoch moment;

[0013] S24, according to the abnormal atomic clock identification information of the calculation epoch, the time difference measurement value sequence and the atomic clock group initial state information are calculated and processed to obtain the normal atomic clock state information of the calculation epoch; the normal atomic clock state information of the calculation epoch includes state variable X k-ok , state variable covariance matrix P k-ok , covariance matrix Q of driving noise, state transition matrix Φ k,k-1 , observation matrix H nclkok-1,nclkok*3 and measurement value matrix Z k ;

[0014] S25, using the preset kalman filter model, the atomic clock state information of the calculation epoch, the observation matrix information and the measurement value matrix information are calculated and processed to obtain the atomic clock state variable calculation information, and the state variable estimation value and the covariance matrix estimation value corresponding to the next epoch of the calculation epoch; the atomic clock state variable calculation information of the calculation epoch includes state variable and state variable covariance matrix;

[0015] The expression of the preset kalman filter model is:

[0016]

[0017] In the formula, is the state variable estimation value of the kth epoch, and the state variable X k-ok value is used in iterative calculation; Φ k+1,k represents the state transition matrix; the kth epoch represents the calculation epoch; represents the predicted value obtained by predicting the state variable X of the kth epoch; P k+1 / k is the covariance matrix of ; P k is the covariance matrix of the state variable of the kth epoch, and the state variable covariance matrix P k-ok value is used in iterative calculation; H k represents the observation matrix of the kth epoch, and the observation matrix H nclkok-1,nclkok*3 value is used in iterative calculation; represents the transpose of H k matrix, R k represents the measurement noise matrix of the kth epoch, Z k represents the measurement value matrix; Q is the covariance matrix of driving noise; K k is the Kalman filter gain matrix; I is the unit matrix;

[0018] S26, updating the initial state information of the atomic clock group by using the preset assignment model, the state variable estimation value and the covariance matrix estimation value, to obtain updated initial state information of the atomic clock group;

[0019] The preset assignment model expression is:

[0020] P k ((ics-1)×3+1:(ics-1)×3+3,(jcs-1)×3+1:(jcs-1)×3+3)

[0021] =P k-ok ((i-1)×3+1:(i-1)×3+3,(j-1)×3+1:(j-1)×3+3)

[0022] X k (ics)=[X k-ok ((i-1)×3+1,1)X k-ok ((i-1)×3+2,1)X k-ok ((i-1)×3+3,1)] T

[0023] i=1,2…nclkok j=1,2…nclkok

[0024] In the formula, ics is the original atomic clock serial number corresponding to the i th normal atomic clock, jcs is the original satellite clock serial number corresponding to the j th normal atomic clock; X k represents the initial state variable; P k represents the initial state variable covariance matrix; X k-ok represents the clock state variable estimation value; P k-ok represents the covariance matrix estimation value;

[0025] S27, updating the next epoch of the resolving epoch in the time difference measurement value sequence to the resolving epoch; repeating steps S23-S27 to obtain the atomic clock state variable resolving information of all epochs in the time difference measurement value sequence.

[0026] As an optional implementation, the abnormal data identification processing on the time difference measurement value at the resolving epoch in the time difference measurement value sequence is performed to obtain abnormal atomic clock identification information at the resolving epoch, and the abnormal data identification processing includes:

[0027] The abnormal atomic clock identification clktag is assigned by analyzing the atomic clock data quality:

[0028]

[0029] In the formula, clktag iThe abnormal atomic clock identification information of the i th atomic clock is used to represent whether the i th atomic clock identification data is normal.

[0030] As an optional implementation, the calculation processing of the time difference measurement value sequence and the initial state information of the atomic clock group according to the abnormal atomic clock identification information of the calculation epoch time point includes:

[0031] S241, the abnormal atomic clock identification information of the calculation epoch time point is statistically processed to obtain the number nclkok of normal atomic clocks at the calculation epoch time point; the number nclkok of the normal atomic clocks is:

[0032]

[0033] S242, the initial state information of the atomic clock group is processed to obtain the state variable X k-ok and the state variable covariance matrix P k-ok .

[0034] S243, the covariance matrix of the driving noise of the normal atomic clock is assigned and processed to obtain the covariance matrix Q of the driving noise;

[0035] S244, the state transition matrix of the normal atomic clock is assigned and processed to obtain the state transition matrix Φ k,k-1 .

[0036] S245, the observation matrix of the normal atomic clock is assigned and processed to obtain the observation matrix H nclkok-1,nclkok*3 .

[0037] S246, the measurement value matrix of the normal atomic clock is assigned and processed according to the time difference measurement value sequence to obtain the measurement value matrix Z k .

[0038] As an optional implementation, the processing of the initial state information of the atomic clock group to obtain the state variable X k-ok and the state variable covariance matrix P k-ok includes:

[0039] S2421, the state variable X k-ok is calculated by using a state variable assignment model.

[0040] The state variable assignment model is expressed as:

[0041] X k-ok (i)=[X k ((ics-1)×3+1,1)X k ((ics-1)×3+2,1)X k((i-1) x 3 + 1:(i-1) x 3 + 3, (j-1) x 3 + 1:(j-1) x 3 + 3) T i = 1, 2,..., nclkok

[0042] wherein, X k-ok represents a state variable; X k represents an initial state variable; i represents a serial number of the i th normal atomic clock, and ics represents a serial number of an original satellite clock corresponding to the i th normal atomic clock;

[0043] S2422, calculating by using the state variable covariance matrix assignment model to obtain the state variable covariance matrix P k-ok ;

[0044] The state variable covariance matrix assignment model is represented as:

[0045] P k-ok ((i-1) x 3 + 1:(i-1) x 3 + 3, (j-1) x 3 + 1:(j-1) x 3 + 3)

[0046] = P k ((i-1) x 3 + 1:(i-1) x 3 + 3, (j-1) x 3 + 1:(j-1) x 3 + 3)

[0047] wherein, P k-ok represents a state variable covariance matrix; P k represents an initial state variable covariance matrix; i and j respectively represent serial numbers of the i th and j th normal atomic clocks, and ics and jcs respectively represent serial numbers of original satellite clocks corresponding to the i th and j th normal atomic clocks.

[0048] As an optional implementation, the assignment processing on the covariance matrix of the driving noise of the normal atomic clock is performed to obtain the covariance matrix Q of the driving noise, comprising:

[0049] The covariance matrix Q is represented as:

[0050]

[0051] wherein,

[0052]

[0053] wherein, ics represents a serial number of an original satellite clock corresponding to the i th normal atomic clock; represents a frequency white noise variance of the ics th original clock represents a frequency random walk noise variance of the i th original clock τ represents a sampling interval.

[0054] As an optional implementation, the state transition matrix of the normal atomic clock is valued to obtain a state transition matrix Φ k,k-1 , comprising:

[0055] The state transition matrix Φ k,k-1 is expressed as:

[0056]

[0057] Wherein,

[0058]

[0059] In the formula, τ represents a sampling interval.

[0060] As an optional implementation, the observation matrix of the normal atomic clock is valued to obtain an observation matrix H nclkok-1,nclkok*3 , comprising:

[0061] The observation matrix H of the normal atomic clock nclkok-1,nclkok*3 is expressed as:

[0062]

[0063] In the formula, the row number i represents the number of the normal satellite except the master clock,

[0064] As an optional implementation, the measurement matrix of the normal atomic clock is valued according to the time difference measurement value sequence to obtain a measurement matrix Z k , comprising:

[0065] When the master master clock is normal, the time difference between the other atomic clock and the master master clock is used as an observation, and the measurement matrix Z k is expressed as:

[0066]

[0067] In the formula, clkref1 i represents the time difference measurement value between the i+1th normal atomic clock and the master master clock;

[0068] When the master master clock is abnormal and the standby master clock is normal, the time difference between the atomic clock and the standby master clock is used as an observation, and the measurement matrix Z k is expressed as:

[0069]

[0070] In the formula, clkref2 i represents the time difference measurement value between the i+2th normal atomic clock and the standby master clock.

[0071] Compared with the prior art, the embodiments of the present invention have the following beneficial effects:

[0072] This invention discloses a method for establishing a Kalman filter timescale that takes into account anomalous data. This method can eliminate the impact of interruptions or anomalies in data from one or more atomic clocks on the performance of the Kalman filter timescale. Compared to the standard Kalman filter timescale calculation method, this method can eliminate the influence of interruptions and anomalous observations. Compared to methods that adjust the covariance matrix of prediction parameters, this method directly targets the observation domain. By processing anomalous observation information, it can eliminate the impact of interruptions and anomalies in the time difference observation data of the master clock or other atomic clocks on the timescale, thus obtaining a more stable timescale. The physical signal formed based on this timescale can serve as a time reference for fields requiring high time and frequency stability, such as satellite navigation systems, finance, and power. Attached Figure Description

[0073] Figure 1 This is a schematic diagram of a method for establishing the time scale of Kalman filtering that takes into account abnormal data, as disclosed in an embodiment of the present invention.

[0074] Figure 2 This is a simulation time difference sequence diagram of various clocks disclosed in an embodiment of the present invention;

[0075] Figure 3 This is a time difference diagram of the paper time calculated by the first master clock and the second master clock, as disclosed in an embodiment of the present invention.

[0076] Figure 4 This is an Allen variance plot of various clocks and paper time disclosed in an embodiment of the present invention. Detailed Implementation

[0077] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0078] Example 1

[0079] Please see Figure 1 . Figure 1 This is a schematic diagram illustrating a method for establishing the time scale of a Kalman filter that takes into account outlier data, as disclosed in an embodiment of the present invention. Figure 1The described Kalman filtering time scale establishment method considering abnormal data is applied to a management system, such as a local server or a cloud server for management, and embodiments of the present application are not limited. Figure 1 The Kalman filtering time scale establishment method considering abnormal data disclosed by the embodiments of the present application includes:

[0080] S1, two atomic clocks are selected from an atomic clock group as master clocks of a time scale algorithm, and other atomic clocks are time-keeping clocks, a time difference measurement value sequence between the time-keeping clocks and the master clocks is obtained, and a sampling interval is set as τ; the atomic clock group includes n atomic clocks, and the n is a positive integer not less than 4; one of the two master clocks is a primary master clock, and the other is a backup master clock;

[0081] S2, the time difference measurement value sequence is processed to obtain atomic clock state variable solution information of all epochs in the time difference measurement value sequence; the state variable solution information represents a single atomic clock and a virtual paper time difference;

[0082] S3, any single time-keeping clock is calibrated by using the atomic clock state variable solution information to obtain a calibrated time scale.

[0083] In another optional embodiment, the processing of the time difference measurement value sequence to obtain the atomic clock state variable solution information of all epochs in the time difference measurement value sequence includes:

[0084] S21, initial state information of the atomic clock group is constructed; the initial state information of the atomic clock group includes initial state variables X k , an initial state variable covariance matrix P k , and a covariance matrix Q0 of an initial driving noise η k ;

[0085] S22, an epoch corresponding to any time difference measurement value in the time difference measurement value sequence is selected as a solution epoch;

[0086] S23, the time difference measurement value at the solution epoch time of the time difference measurement value sequence is subjected to abnormal data identification processing to obtain abnormal atomic clock identification information at the solution epoch time;

[0087] S24, according to the abnormal atomic clock identification information at the solution epoch time, the time difference measurement value sequence and the initial state information of the atomic clock group are calculated and processed to obtain normal atomic clock state information at the solution epoch time; the normal atomic clock state information at the solution epoch time includes state variables X k-ok , a state variable covariance matrix P k-ok , a covariance matrix Q of a driving noise, a state transition matrix Φ k,k-1 , and an observation matrix Hnclkok-1,nclkok*3 and the measurement matrix Z k ;

[0088] S25, using a preset Kalman filter model, the atomic clock state information, observation matrix information and measurement value matrix information of the epoch are calculated and processed to obtain atomic clock state variable calculation information, and the state variable estimation value and covariance matrix estimation value corresponding to the next epoch of the calculation epoch; the atomic clock state variable calculation information of the calculation epoch includes state variable and covariance matrix of state variable;

[0089] The expression of the preset Kalman filter model is:

[0090]

[0091] In the formula, The state variable estimation value of the kth epoch is used in iterative calculation, and the state variable X k-ok Value; Φ k+1,k Indicates the state transition matrix; the kth epoch represents the calculation epoch; Indicates the predicted value obtained by predicting the state variable Of the kth epoch; P k+1 / k Is the covariance matrix of ; P k The covariance matrix of the state variable of the kth epoch is used in iterative calculation, and the state variable covariance matrix P k-ok Value; H k The observation matrix of the kth epoch is used in iterative calculation, and the observation matrix H nclkok-1,nclkok*3 Value; The transpose of H k Matrix is represented by R k The measurement noise matrix of the kth epoch is represented by Z k Indicates the measurement value matrix; Q is the covariance matrix of the driving noise; K k Is the Kalman filter gain matrix; I is the unit matrix;

[0092] S26, using a preset assignment model, the state variable estimation value and covariance matrix estimation value are used to update the initial state information of the atomic clock group to obtain the updated initial state information of the atomic clock group;

[0093] The expression of the preset assignment model is:

[0094] P k ((ics-1)×3+1:(ics-1)×3+3,(jcs-1)×3+1:(jcs-1)×3+3)

[0095] =Pk-ok ((i-1) x 3 + 1:(i-1) x 3 + 3, (j-1) x 3 + 1:(j-1) x 3 + 3)

[0096] X k (ics) = [X k-ok ((i-1) x 3 + 1, 1) X k-ok ((i-1) x 3 + 2, 1) X k-ok ((i-1) x 3 + 3, 1)] T

[0097] i = 1, 2…nclkok j = 1, 2…nclkok

[0098] In the formula, ics is the original atomic clock sequence number corresponding to the i th normal atomic clock, jcs is the original satellite clock sequence number corresponding to the j th normal atomic clock; X k represents an initial state variable; P k represents an initial state variable covariance matrix; X k-ok represents an atomic clock state variable estimation value; P k-ok represents a covariance matrix estimation value;

[0099] S27, the next epoch of the resolving epoch in the time difference measurement value sequence is updated to the resolving epoch; the steps S23-S27 are looped to obtain the atomic clock state variable resolving information of all epochs in the time difference measurement value sequence.

[0100] In yet another optional embodiment, the initial state information of the atomic clock group is constructed, comprising:

[0101] The initial state variable X k of the atomic clock group is:

[0102]

[0103] Wherein, x i and y i are the atomic clock i time difference and frequency difference respectively, d i is the frequency drift;

[0104] The covariance matrix Q0 of the initial driving noise η k is:

[0105]

[0106] Wherein,

[0107]

[0108] The initial state variable covariance matrix P k is:

[0109] P k = 1000 * Q0;

[0110] wherein, represents the frequency white noise variance of the i-th clock represents the frequency random walk noise variance of the i-th clock The calculation method can be found in the published doctoral thesis of Dr. Guo Hai Rong, entitled "Research on Theory and Method of Navigation Satellite Atomic Clock Time and Frequency Characteristics Analysis".

[0111] The abnormal data identification processing of the time difference measurement value at the calculation epoch time in the time difference measurement value sequence is performed to obtain abnormal atomic clock identification information at the calculation epoch time, and the abnormal data identification processing comprises the following steps:

[0112] The abnormal atomic clock identification clktag is valued through atomic clock data quality analysis:

[0113]

[0114] wherein, clktag i is the i-th atomic clock identification, which represents whether the i-th atomic clock identification data is normal.

[0115] In another optional embodiment, the abnormal data identification processing of the time difference measurement value at the calculation epoch time in the time difference measurement value sequence is performed to obtain abnormal atomic clock identification information at the calculation epoch time, and the abnormal data identification processing comprises the following steps:

[0116] Suppose the atomic clock time difference measurement value sequence is clk1, clk2, … clk m , and the corresponding observation time is t1, t2, … t m , the sampling interval is Δt, and the following steps are adopted for abnormal identification.

[0117] First-order difference is performed on the observation time to obtain a time difference sequence of the previous and subsequent observations:

[0118] dt i = t i+1 -t i (i = 1, 2 … m-1)

[0119] According to whether the i-th time difference is the sampling interval, it is judged whether the time difference data is missing at the i+1-th epoch:

[0120]

[0121] The frequency difference time sequence a1, a2, … a m-1 is calculated:

[0122]

[0123] A first-order polynomial is applied to the frequency difference time series, and the fitted residual series are calculated, denoted as resa1, resa2, ..., resa. m-1 ;

[0124] The median of the residual time series is calculated and denoted as md:

[0125] md = Median{resa i}

[0126] Calculate the critical parameter MAD:

[0127] MAD = Median{|resa i -md| / 0.6745}

[0128] Determine if atomic clock data is abnormal

[0129]

[0130] Assign values ​​to the anomaly flags of the atomic clock time difference data:

[0131]

[0132] In the formula, clktag i This is the identifier for the i-th atomic clock, indicating whether the data for that i-th atomic clock is normal.

[0133] In another optional embodiment, the step of calculating and processing the time difference measurement sequence and the initial state information of the atomic clock group based on the abnormal atomic clock identification information at the solved epoch time to obtain the normal atomic clock state information at the solved epoch time includes:

[0134] S241, perform statistical processing on the abnormal atomic clock identification information of the solved epoch time to obtain the number of normal atomic clocks nclkok at the solved epoch time; the value of the number of normal atomic clocks nclkok is:

[0135]

[0136] S242, Process the initial state information of the atomic clock group to obtain the state variable X. k-ok State variable covariance matrix P k-ok ;

[0137] S243, assign values ​​to the covariance matrix of the driving noise of a normal atomic clock to obtain the covariance matrix Q of the driving noise;

[0138] S244, assign values ​​to the state transition matrix of a normal atomic clock to obtain the state transition matrix Φ. k,k-1 ;

[0139] S245, the observation matrix of the normal atomic clock is valued to obtain the observation matrix H nclkok-1,nclkok*3 ;

[0140] S246, the measurement matrix of the normal atomic clock is valued according to the time difference measurement value sequence to obtain the measurement matrix Z k .

[0141] The initial state information of the atomic clock group is processed to obtain the state variable X k-ok , the state variable covariance matrix P k-ok , comprising:

[0142] S2421, the state variable X k-ok is calculated by using the state variable valuation model.

[0143] The state variable valuation model is represented as:

[0144] X k-ok (i)=[X k ((ics-1)×3+1,1)X k ((ics-1)×3+2,1)X k ((ics-1)×3+3,1)] T i=1,2…nclkok

[0145] In the formula, X k-ok represents the state variable; X k represents the initial state variable; i represents the serial number of the i-th normal atomic clock, and ics represents the original satellite clock serial number corresponding to the i-th normal atomic clock.

[0146] S2422, the state variable covariance matrix P k-ok is calculated by using the state variable covariance matrix valuation model.

[0147] The state variable covariance matrix valuation model is represented as:

[0148] P k-ok ((i-1)×3+1:(i-1)×3+3,(j-1)×3+1:(j-1)×3+3)

[0149] =P k ((ics-1)×3+1:(ics-1)×3+3,(jcs-1)×3+1:(jcs-1)×3+3)

[0150] In the formula, P k-ok represents the state variable covariance matrix; P krepresents an initial state variable covariance matrix; i and j represent the serial numbers of the i th and j th normal atomic clocks, respectively, and ics and jcs represent the original satellite clock serial numbers corresponding to the i th and j th normal atomic clocks, respectively.

[0151] In yet another optional embodiment, the covariance matrix of the driving noise of the normal atomic clocks is assigned to obtain a driving noise covariance matrix Q, including:

[0152] The covariance matrix Q is expressed as:

[0153]

[0154] wherein,

[0155]

[0156] In the formula, ics represents the original satellite clock serial number corresponding to the i th normal atomic clock; represents the frequency white noise variance of the ics th original clock represents the frequency random walk noise variance of the i th original clock τ represents a sampling interval.

[0157] In yet another optional embodiment, the state transition matrix of the normal atomic clocks is assigned to obtain a state transition matrix Φ k,k-1 , including:

[0158] The state transition matrix Φ k,k-1 is expressed as:

[0159]

[0160] wherein,

[0161]

[0162] In the formula, τ represents a sampling interval.

[0163] In yet another optional embodiment, the observation matrix of the normal atomic clocks is assigned to obtain an observation matrix H nclkok-1,nclkok*3 , including:

[0164] The observation matrix H nclkok-1,nclkok*3 of the normal atomic clocks is expressed as:

[0165]

[0166] In the formula, the row number i represents the number of normal satellites other than the master clock,

[0167] In yet another optional embodiment, the measurement value matrix Z of normal atomic clocks is assigned according to the time difference measurement value sequence k , comprising:

[0168] When the main master clock is normal, the time difference between the other atomic clocks and the main master clock is used as the observation value, and the measurement value matrix Z k is expressed as:

[0169]

[0170] In the formula, clkref1 i represents the time difference measurement value between the i+1th normal atomic clock and the main master clock;

[0171] When the main master clock is abnormal and the backup master clock is normal, the time difference between the atomic clocks and the backup master clock is used as the observation value, and the measurement value matrix Z k is expressed as:

[0172]

[0173] In the formula, clkref2 i represents the time difference measurement value between the i+2th normal atomic clock and the backup master clock.

[0174] It should be noted that when both master clocks are abnormal, step S1 is repeated, that is, a normal atomic clock is reselected as a master clock.

[0175] Embodiment two

[0176] The embodiment of the present application is used to perform time scale calculation by using the disclosed Kalman filter time scale establishment method considering abnormal data, so as to verify the effect of resisting abnormal time difference data.

[0177] Eleven hydrogen atomic clock time difference data are simulated, and the simulation parameters are shown in Table 1. Under normal circumstances, the first hydrogen atomic clock is set as the master clock by the time scale algorithm. In order to verify the continuity of switching the second hydrogen atomic clock as the master clock when the first hydrogen atomic clock is abnormal, the first hydrogen atomic clock is set to be interrupted after continuous operation for 120 days, and the second hydrogen atomic clock of the Beijing master station is automatically switched to the master clock.

[0178] Table 1 Hydrogen atomic clock time difference data simulation parameters

[0179]

[0180]

[0181] The simulation results are shown in Figure 2 , Figure 3 ,Figure 4 .

[0182] Figure 2 The simulation time difference sequence of each time-keeping clock is continuously operated for 180 days, and the time difference can reach the order of 1e-5s.

[0183] Figure 3 The time difference of paper time calculated by the first master clock and the second master clock is continuously operated for 120 days, the time difference of paper time calculated by the first master clock is in the order of 1e-8s; after the 120th day, the master clock is switched to the second master clock, the paper time calculated by the master clock is kept continuous, and is operated to the 180th day, the time difference of paper time calculated by the second master clock is in the order of 1e-8s.

[0184] Figure 4 The Allan variance of each time-keeping clock and paper time is shown, wherein the green circle is the Allan deviation of paper time, and the rest is the Allan deviation of each time-keeping atomic clock. It can be seen that the frequency stability of paper time is higher than that of a single hydrogen clock, the paper time ten-second stability of 11 hydrogen atomic clocks established by the Kalman filtering time scale calculation method considering abnormal data in the patent is in the order of 1e-15, the day stability is in the order of 1e-16~1e-15, and the 10-day stability is in the order of 1e-15.

[0185] The device embodiments described above are only schematic, wherein the modules illustrated as separate components can or can not be physically separated, and the components illustrated as modules can or can not be physical modules, that is, they can be located in one place, or can be distributed on multiple network modules. Part or all of the modules can be selected to achieve the purpose of the embodiment scheme according to actual needs. Those skilled in the art can understand and implement it without creative labor.

[0186] Those skilled in the art can clearly understand the implementation of the various embodiments by means of software and necessary general hardware platforms through the above specific description of the embodiments, and of course, the embodiments can also be implemented by hardware. Based on such understanding, the above technical solutions, essentially or in the sense of contribution to the prior art, can be embodied in the form of a software product. The computer software product can be stored in a computer readable storage medium, which includes a Read-Only Memory (ROM), a Random Access Memory (RAM), a Programmable Read-only Memory (PROM), an Erasable Programmable Read Only Memory (EPROM), a One-time Programmable Read-Only Memory (OTPROM), an Electrically-Erasable Programmable Read-Only Memory (EEPROM), a Compact Disc Read-Only Memory (CD-ROM), or other optical disk storage, a magnetic disk storage, a magnetic tape storage, or any other medium that can be used to carry or store data in a computer readable manner.

[0187] Finally, it should be noted that: the disclosed kalman filter time scale establishment method considering abnormal data disclosed by the embodiments of the present application is only the preferred embodiment of the present application, and is used to illustrate the technical solutions of the present application, but not to limit it; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that; the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced by equivalents; and these modifications or replacements do not make the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for establishing a time scale of Kalman filtering taking into account abnormal data, characterized in that, The method comprises: S1, selecting two atomic clocks from an atomic clock group as master clocks of a time scale algorithm, other atomic clocks as time keeping clocks, obtaining a time difference measurement value sequence between the time keeping clocks and the master clocks, and setting a sampling interval as τ; the atomic clock group comprises n atomic clocks, and n is a positive integer not less than 4; one of the two master clocks is used as a primary master clock, and the other is used as a backup master clock; S2, processing the time difference measurement value sequence to obtain atomic clock state variable calculation information of all epochs in the time difference measurement value sequence; the state variable calculation information represents a single atomic clock and a virtual paper time difference; specifically: S21, constructing an initial state information of an atomic clock group; the initial state information of the atomic clock group comprises initial state variables X of n atomic clocks k , an initial state variable covariance matrix P k , and a covariance matrix Q0 of initial driving noise η k ; S22, selecting an epoch corresponding to any time difference measurement value in the time difference measurement value sequence as a calculation epoch; S23, performing abnormal data identification processing on the time difference measurement value at the calculation epoch time in the time difference measurement value sequence to obtain abnormal atomic clock identification information at the calculation epoch time; S24, according to the abnormal atomic clock identification information of the calculation epoch time, the time difference measurement value sequence and the atomic clock group initial state information are calculated and processed to obtain the normal atomic clock state information of the calculation epoch time; the normal atomic clock state information of the calculation epoch time includes state variable X k-ok , state variable covariance matrix P k-ok , covariance matrix Q of driving noise, state transition matrix Φ k,k-1 , observation matrix H nclkok-1,nclkok*3 and measurement value matrix Z k ; S25, using a preset kalman filter model to perform calculation processing on the atomic clock state information, the observation matrix information and the measurement value matrix information at the calculation epoch time to obtain atomic clock state variable calculation information, and taking the atomic clock state variable calculation information as state variable estimation values and covariance matrix estimation values corresponding to the next epoch of the calculation epoch; The expression of the preset kalman filter model is: wherein, is the state variable estimate of the kth epoch, and the state variable X k-ok is updated in the iterative calculation; Φ k+1,k represents a state transition matrix; the kth epoch represents a calculation epoch; represents a predicted value obtained by performing a prediction on the state variable X of the kth epoch; P k+1 / k is a covariance matrix of X ; P k is the covariance matrix of the state variable of the kth epoch, and the state variable covariance matrix P k-ok is updated in the iterative calculation; H k represents an observation matrix of the kth epoch, and the observation matrix H nclkok-1,nclkok*3 is updated in the iterative calculation; represents a transpose of H k ; R k represents a measurement noise matrix of the kth epoch, and Z k represents a measurement value matrix; Q is a covariance matrix of driving noise; K k is a Kalman filter gain matrix; and I is a unit matrix; S26, using a preset assignment model to update the initial state information of the atomic clock group by using the state variable estimation values and the covariance matrix estimation values to obtain updated initial state information of the atomic clock group; The expression of the preset assignment model is: P k ((ics-1) x 3 + 1 : (ics-1) x 3 + 3, (jcs-1) x 3 + 1 : (jcs-1) x 3 + 3) = P k-ok ((i-1) x 3 + 1 : (i-1) x 3 + 3, (j-1) x 3 + 1 : (j-1) x 3 + 3) X k (ics) = [X k-ok ((i-1) x 3 + 1, 1) X k-ok ((i-1) x 3 + 2, 1) X k-ok ((i-1) x 3 + 3, 1)] T i = 1, 2…nclkok j = 1, 2…nclkok where ics is the original satellite clock sequence number corresponding to the ith normal atomic clock, and jcs is the original satellite clock sequence number corresponding to the jth normal atomic clock; X k characterizes the initial state variable; P k characterizes the initial state variable covariance matrix; X k-ok characterizes the clock state variable estimate; P k-ok characterizes the covariance matrix estimate; nclkok characterizes the number of normal atomic clocks; S27, updating the next epoch of the calculation epoch in the time difference measurement value sequence to the calculation epoch; repeating steps S23-S27 to obtain atomic clock state variable calculation information of all epochs in the time difference measurement value sequence; S3, calibrating any single time keeping clock by using the atomic clock state variable calculation information to obtain a calibrated time scale.

2. The method of claim 1, wherein, The abnormal data identification processing on the time difference measurement value at the calculation epoch time in the time difference measurement value sequence to obtain abnormal atomic clock identification information at the calculation epoch time comprises: The abnormal atomic clock identification clktag is assigned by analyzing the atomic clock data quality: In the formula, clktag(i) is the identification of the i th atomic clock, which represents whether the identification data of the i th atomic clock is normal.

3. The method of claim 2, wherein the method further comprises: The calculation processing on the time difference measurement value sequence and the initial state information of the atomic clock group according to the abnormal atomic clock identification information at the calculation epoch time to obtain normal atomic clock state information at the calculation epoch time comprises: S241, performing statistical processing on the abnormal atomic clock identification information at the calculation epoch time to obtain the number nclkok of normal atomic clocks at the calculation epoch time; the number nclkok of normal atomic clocks is: S242, processing the atomic clock group initial state information to obtain state variable X k-ok , state variable covariance matrix P k-ok ; S243, performing assignment processing on the covariance matrix of the drive noise of the normal atomic clock to obtain the covariance matrix Q of the drive noise; S244, the state transition matrix of normal atomic clock is valued, and the state transition matrix Φ is obtained k,k-1 ; S245, the observation matrix of normal atomic clock is valued, and the observation matrix H is obtained nclkok-1,nclkok*3 ; S246, according to the time difference measurement value sequence, the measurement value matrix of the normal atomic clock is assigned to obtain the measurement value matrix Z k .

4. The method for kalman filter time scale establishment considering abnormal data according to claim 3, characterized in that, The processing of the atomic clock group initial state information obtains a state variable X k-ok , a state variable covariance matrix P k-ok , comprising: S2421、Using the state variable assignment model to calculate, the state variable X is obtained k-ok ; The state variable assignment model is represented as: X k-ok (i) = [X k ((ics-1) x 3 + 1, 1) X k ((ics-1) x 3 + 2, 1) X k ((ics-1) x 3 + 3, 1)] T i = 1, 2...nclkok wherein X k-ok represents a state variable; X k represents an initial state variable; i represents a serial number of the i-th normal atomic clock, and ics represents a serial number of an original satellite clock corresponding to the i-th normal atomic clock. S2422. Calculate the state variable covariance matrix P using the state variable covariance matrix assignment model. k-ok ; The state variable covariance matrix assignment model is represented as: P k-ok ((i-1) x 3 + 1 : (i-1) x 3 + 3, (j-1) x 3 + 1 : (j-1) x 3 + 3) = P k ((ics-1) x 3 + 1 : (ics-1) x 3 + 3, (jcs-1) x 3 + 1 : (jcs-1) x 3 + 3) where P k-ok denotes the state variable covariance matrix; P k denotes the initial state variable covariance matrix; i, j represent the serial numbers of the i-th and j-th normal atomic clock, ics, jcs represent the original satellite clock serial numbers corresponding to the i-th and j-th normal atomic clock, respectively.

5. The method for kalman filter time scale establishment considering abnormal data according to claim 3, characterized in that, The covariance matrix of the driving noise of the normal atomic clock is valued to obtain a covariance matrix Q of the driving noise, comprising: The covariance matrix Q is expressed as: Wherein, In the formula, ics represents the original satellite clock sequence number corresponding to the i-th normal atomic clock; Frequency white noise variance of the i-th original clock Frequency random walk noise variance of the i-th original clock τ represents the sampling interval.

6. The method for kalman filter time scale establishment considering abnormal data according to claim 3, characterized in that, The state transition matrix of the normal atomic clock is valued to obtain a state transition matrix Φ k,k-1 , comprising: The state transition matrix Φ k,k-1 is represented as: Wherein, In the formula, τ represents a sampling interval.

7. The method for kalman filter time scale establishment considering abnormal data according to claim 3, characterized in that, The observation matrix of the normal atomic clock is valued to obtain an observation matrix H nclkok-1,nclkok*3 , comprising: The observation matrix H of the normal atomic clock nclkok-1,nclkok*3 is represented as: wherein the line number i represents the number of the normal satellite except for the master clock, 8. The method for kalman filter time scale establishment considering abnormal data according to claim 3, characterized in that, The measurement value matrix Z of the normal atomic clock is obtained by assigning the time difference measurement value sequence k , comprising: When the primary master clock is normal, the time difference between the other atomic clock and the primary master clock is used as an observation, and the measurement value matrix Z k is expressed as: In the formula, clkref1 i represents the measured time difference between the (i+1)th normal atomic clock and the master clock. When the primary master clock is abnormal and the backup master clock is normal, a time difference between the atomic clock and the backup master clock is used as an observation value, and the measurement value matrix Z k is represented as: where clkref2 i represents the measured time difference between the (i+2)th normal atomic clock and the backup master clock.

Images