CN110554597B - Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering - Google Patents

Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering Download PDF

Info

Publication number
CN110554597B
CN110554597B CN201910787030.4A CN201910787030A CN110554597B CN 110554597 B CN110554597 B CN 110554597B CN 201910787030 A CN201910787030 A CN 201910787030A CN 110554597 B CN110554597 B CN 110554597B
Authority
CN
China
Prior art keywords
clock
time scale
cesium
hydrogen
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201910787030.4A
Other languages
Chinese (zh)
Other versions
CN110554597A (en
Inventor
姜萌
董绍武
武文俊
宋会杰
高喆
安卫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Chinese Academy of Sciences
National Time Service Center of CAS
Original Assignee
University of Chinese Academy of Sciences
National Time Service Center of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Chinese Academy of Sciences, National Time Service Center of CAS filed Critical University of Chinese Academy of Sciences
Priority to CN201910787030.4A priority Critical patent/CN110554597B/en
Publication of CN110554597A publication Critical patent/CN110554597A/en
Application granted granted Critical
Publication of CN110554597B publication Critical patent/CN110554597B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G04—HOROLOGY
    • G04F—TIME-INTERVAL MEASURING
    • G04F5/00—Apparatus for producing preselected time intervals for use as timing standards
    • G04F5/14—Apparatus for producing preselected time intervals for use as timing standards using atomic clocks

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Stabilization Of Oscillater, Synchronisation, Frequency Synthesizers (AREA)
  • Complex Calculations (AREA)

Abstract

The invention provides a hydrogen cesium time scale fusion method based on Vondark-Cepek filtering, which comprises the steps of firstly, respectively generating a clock group time scale for a hydrogen clock group and a cesium clock group by using an exponential filtering method, then selecting key parameters of Vondark-Cepek combined filtering according to the principle of least square, and further performing performance enhancement on the cesium clock group time scale through differential information of a hydrogen clock group time scale time sequence so as to obtain the hydrogen cesium fusion time scale. The invention can give consideration to the characteristics of two types of hydrogen and cesium atomic clocks while effectively inhibiting noise, and can be fused on the scale of a clock group to generate a time scale with good long-term and short-term stability.

Description

Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering
Technical Field
The invention belongs to the field of time frequency and signal processing, and relates to a time scale fusion method.
Background
At present, high-precision time becomes a vital parameter in national science and technology, economy, military and social life, is widely applied to the fields of navigation, electric power, communication, aviation, national defense and the like, and has a fundamental effect on improving the national scientific research level as the most basic physical quantity. For this reason, each time frequency laboratory in the world has its own set of atomic clocks, each of which can produce a time scale. But timekeeping cannot rely on only a single atomic clock because each physical device may have an exception. The time scale algorithm is to calculate a standard time according to the time of each clock of the time keeping clock group, and aims to find a mode to minimize the uncertainty and maximize the stability of the algorithm.
The hydrogen cesium atomic clock is two atomic frequency standards with different types, and the hydrogen cesium atomic clock has good short-term stability and long-term stability respectively, and is two punctual atomic clocks which are most common in laboratories at various times around the world currently. The fusion method of hydrogen cesium time scale mainly researches effective combination of two different atomic clocks to generate more stable and accurate time scale. A time laboratory emphasizes the stability of the time scale.
At present, many time-keeping laboratories adopt a mode of hydrogen and cesium for mutual reference to perform fusion of time scales of two different types of atomic clocks: respectively taking the time of a hydrogen clock, a cesium clock or a cesium clock group as a reference, performing clock error measurement on another atomic clock to form a time sequence, correcting the short-term fluctuation of the cesium clock in the cesium clock group or subtracting the speed and frequency drift of the hydrogen clock, and then obtaining a final time scale by weighted averaging. In the above research method, a single clock of another atomic clock is evaluated and corrected mainly based on hydrogen or cesium, and then all clocks in the time keeping clock group are weighted and averaged together to calculate the final time scale. The method has the following disadvantages: when the hydrogen atomic clock is used as a reference to measure the cesium atomic clock, the main noise is phase white noise, and after the phase white noise is filtered by a mathematical method, the short-term stability of the time scale of the phase white noise is still influenced by the noise of the cesium atomic clock. And the hydrogen clock is corrected by the time scale of the cesium clock group, and when the hydrogen clock fails, the reliability of the generated time scale cannot be guaranteed.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a hydrogen cesium time scale fusion method based on Vondark-Cepek filtering, which can give consideration to the characteristics of two types of atomic clocks of hydrogen cesium while effectively inhibiting noise, and can be used for fusion in a clock group scale level to generate a time scale with good long-term and short-term stability.
The technical scheme adopted by the invention for solving the technical problem comprises the following steps:
step 1, grouping cesium atomic clocks and hydrogen atomic clocks, and respectively passing through clock differencesPrediction, frequency estimation and weight estimation steps produce an exponentially filtered time scale TA of the cesium clock setCsAnd hydrogen time scale TAH；
Step 2: carrying out primary difference on the hydrogen clock group time scale sequence to obtain corresponding frequency data;
step 3, Fourier transform is carried out on the time scale sequence of the cesium clock set to convert the time scale sequence into a frequency domain, and the frequency f or the period P of a signal needing to be suppressed is determined through spectrum analysis; calculating a time scale smoothing factor of the cesium clock group according to the least square theory
Figure GDA0002895421270000021
And hydrogen clock group time scale smoothing factor
Figure GDA0002895421270000022
Wherein T and T' are each TACsAnd TAHA corresponding frequency response;
step 4, taking the cesium clock group time scale and the hydrogen clock group frequency sequence as input, using a smoothing factor in a Vondark-Cepek combined filtering method, generating a time scale fusion result, and realizing the fusion of the cesium clock group time scale and the hydrogen clock group time scale;
establishing a hydrogen cesium time scale fusion model, namely defining a cesium clock group time scale (TA)Cs) Data is m (t)i) Is a time tiI-1, 2,3, … at tiThe derivative of the time instant is defined as m' (t) ═ m (t)i+1)-m(ti)]/(ti+1-ti)；M(tj) Is the hydrogen clock time scale (TA)H) The derivative of which is defined as M' (t) ═ M (t)j+1)-M(tj)]/(tj+1-tj) (ii) a Wherein j is 1,2,3, …; the physical meanings of M '(t) and M' (t) are consistent, both being the rates of the time scale sequence; therefore, when t isi＝tjWhen M '(t) is M' (t);
m(MJD2)＝m(MJD1)+M′(tj)△tj (3)
wherein MJD is the simplified julian day;
because the time intervals of the time scales of the selected cesium clock group and the hydrogen clock group are the same, the clock difference estimation of the next moment of the cesium clock group with the formula (3) can be represented by the sum of the clock difference of the cesium clock group at the moment and the product of the time scale rate and the time interval of the hydrogen clock group;
suppose input observation sequence 1 is represented as y'jCorresponding time scale is xjWeight p ofj(ii) a First derivative of observation sequence 2 is expressed as
Figure GDA0002895421270000023
Corresponding time scale is
Figure GDA0002895421270000024
The weight is
Figure GDA0002895421270000025
The output sequence is then the smooth curve to be obtained, denoted yiCorresponding time scale is xi(ii) a Three quantities are defined:
1) smoothness of the curve:
Figure GDA0002895421270000026
in the formula (I), the compound is shown in the specification,
Figure GDA0002895421270000027
is unknown, based on the smoothed data, on its third derivative
Figure GDA0002895421270000028
Carrying out estimation; at two points [ x ]i+1,yi+1]And [ x ]i+2,yi+2]The smooth curve between is defined as a lagrange polynomial L of third order derived from four adjacent points i, i +1, i +2, i +3i(x)：
Figure GDA0002895421270000031
The third derivative is:
Figure GDA0002895421270000032
the third derivative between each pair of data points is set to a constant, and equation (1) can be expressed as:
Figure GDA0002895421270000033
wherein
Figure GDA0002895421270000034
Figure GDA0002895421270000035
Figure GDA0002895421270000036
Figure GDA0002895421270000037
This smooth definition means that the ideal smoothing function is a quadratic function in time, the first derivative of which is a linear function;
2) fidelity of the smoothed curve to the observed values:
Figure GDA0002895421270000038
3) fidelity of the smooth curve to the first derivative of the observed value:
Figure GDA0002895421270000039
first derivative of
Figure GDA00028954212700000310
Can be based on the smoothing function value yiAnd (4) showing.
The method comprises the following steps that 1, atomic clock difference data measured in a laboratory are preprocessed, atomic clocks with the best long-term and short-term stability in a cesium atomic clock group and a hydrogen atomic clock group are respectively selected as reference clocks, and clock difference data of other atomic clocks in the same clock group and the reference clocks are measured by using a time interval counter.
Removing abnormal points from the clock error data obtained by preprocessing by adopting a 3 sigma rule, detecting the continuity of the clock error data, and repairing a clock error sequence with continuous missing data less than 5 measurement periods by adopting an interpolation method; and if the clock error data continuously lack for more than 5 measurement periods, eliminating the atomic clock in the time keeping clock group.
In the step 1, if the weight of each atomic clock is greater than the set maximum weight upper limit, the weight of the atomic clock is set as the maximum weight; h 'is'i(t) is the atomic clock H at time tiTime correction quantity is weight, atomic clock data, h'i(t) by the formula
Figure GDA0002895421270000041
Calculation of where xij(t) is an atomic clock HjAnd atomic clock HiThe clock difference at time t, N being the number of atomic clocks involved, xi(t) is the difference between the mean time scales of the laboratory atomic clock and the laboratory atomic clock set, wiIs the weight of each atomic clock.
In the step 1, the maximum weight is set to be 2.5/N.
In the step 4, a time scale TA of the cesium clock group is definedCsData is m (t)i) Is a time tiI-1, 2,3, … at tiThe derivative of time m' (t) ═ m (t)i+1)-m(ti)]/(ti+1-ti) (ii) a Definition of M (t)j) Is hydrogen clock group time scale TAHThe derivative M' (t) thereof is [ M (t) ]j+1)-M(tj)]/(tj+1-tj) Wherein, in the step (A),j ═ 1,2,3, …; when t isi＝tjWhen M ' (t) is M ' (t), M (MJD2) is M (MJD1) + M ' (t)j)△tjMJD denotes reduced julian days; suppose input observation sequence 1 is represented as y'jCorresponding time scale is xjWeight p ofj(ii) a First derivative of observation sequence 2 is expressed as
Figure GDA0002895421270000042
Corresponding time scale is
Figure GDA0002895421270000043
The weight is
Figure GDA0002895421270000044
The output sequence is then the smooth curve y to be obtainediCorresponding time scale is xiDefining the smoothness of the curve
Figure GDA0002895421270000045
In the formula (I), the compound is shown in the specification,
Figure GDA0002895421270000046
is unknown, based on the smoothed data, on its third derivative
Figure GDA0002895421270000047
Carrying out estimation; at two points [ x ]i+1,yi+1]And [ x ]i+2,yi+2]The smooth curve between is defined as a lagrange polynomial of third order derived from four adjacent points i, i +1, i +2, i +3
Figure GDA0002895421270000048
Third derivative thereof
Figure GDA0002895421270000049
Defining the fidelity of a smooth curve to an observed value
Figure GDA00028954212700000410
Defining the fidelity of the smooth curve to the first derivative of the observed value
Figure GDA00028954212700000411
The invention has the beneficial effects that: fusion of hydrogen and cesium time scale layers is realized, and optimal weighting under three conditions of absolute smoothness and absolute fidelity of the time scale of the cesium clock group and absolute fitting of the first-order derivative of the time scale of the hydrogen clock group is discussed. On the basis of an exponential filtering algorithm, filtering on a clock group time scale is further realized, and the influence of short-term fluctuation of the cesium clock on a fusion result is effectively improved. Meanwhile, the reliability of the weighted average time scale only once is further improved by fusing the two time scales. The method can effectively inhibit noise, and fully utilize the characteristics of the hydrogen cesium atomic clock, so that the long-term and short-term stability of the fusion time scale is good.
Drawings
Fig. 1 is a schematic flow chart of the present invention for generating a fusion time scale (taking three cesium clocks and three hydrogen clocks as an example).
Fig. 2 is a schematic time scale diagram, in which (a) is a time scale of a cesium clock set and (b) is a time scale of a hydrogen clock set.
FIG. 3 is a diagram illustrating the fusion results of the time scale of the present invention, wherein (a) is the fusion result and (b) is the Allan bias.
Detailed Description
The present invention will be further described with reference to the following drawings and examples, which include, but are not limited to, the following examples.
The invention provides a hydrogen cesium time scale fusion method based on Vondark-Cepek filtering, which comprises the steps of firstly, respectively generating a clock group time scale for a hydrogen clock group and a cesium clock group by using an exponential filtering method, then selecting key parameters of Vondark-Cepek combined filtering according to the principle of least square, and further performing performance enhancement on the cesium clock group time scale through difference information of a hydrogen clock group time scale time sequence so as to obtain the hydrogen cesium fusion time scale.
The invention specifically comprises the following steps:
step 1: and preprocessing the atomic clock error data measured by the laboratory. Firstly, selecting an atomic clock with the best long-term and short-term stability in a clock group as a reference clock, and measuring clock difference data of other atomic clocks and the reference clock by using a time interval counter. The clock error is measured once per hour. And removing abnormal points by adopting a 3 sigma rule, and simultaneously detecting the continuity of the clock error data. Repairing a clock error sequence with continuous missing data for less than 5 hours by adopting an interpolation method; and if the clock error data are continuously missing for more than 5 hours, removing the clock from the time keeping clock group.
Step 2: grouping cesium atomic clocks and hydrogen atomic clocks, and generating exponential filtering cesium clock time scale TA through clock error prediction, frequency estimation and weight estimation steps (once per hour)CsAnd hydrogen time scale TAH. (calculation of a time scale of 3 months is exemplified below)
And if the weight of each atomic clock is greater than the set maximum weight upper limit, setting the weight of the atomic clock as the maximum weight. According to a preferred embodiment of the method, the maximum weighting is set to 2.5/N. Wherein N is the number of atomic clocks involved in the calculation. H 'is'i(t) is the atomic clock H at time tiTime correction quantity is weight, atomic clock data, h'i(t) is calculated by the following formula:
Figure GDA0002895421270000051
wherein x isij(t) is an atomic clock HjAnd atomic clock HiThe clock difference at time t; x is the number ofi(t) is the difference between the mean time scales of the laboratory atomic clock and the laboratory atomic clock set, wiIs the weight of each atomic clock.
And step 3: and (3) carrying out primary difference on the hydrogen clock group time scale sequence obtained by calculation in the step (2) to obtain corresponding frequency data.
And 4, step 4: and transforming the time scale sequence of the cesium clock set into a frequency domain through Fourier transform, and determining a signal f or a period P to be suppressed through spectrum analysis. Determining TA based on the frequency f or period P, experience and actual algorithmCsAnd TAHCorresponding frequency responses T and T', calculating cesium clock set according to least square theoryTime scale smoothing factor epsilon and hydrogen clock time scale smoothing factor
Figure GDA0002895421270000061
Figure GDA0002895421270000062
And 5: and (3) taking the time scale of the cesium clock set generated by calculation in the step (2) and the frequency sequence of the hydrogen clock set generated in the step (3) as input, using the smoothing factor obtained in the step (4) in a Vondark-Cepek combined filtering method, generating a time scale fusion result, and realizing the fusion of the time scales of the cesium clock set and the hydrogen clock set.
Establishing a hydrogen cesium time scale fusion model, namely defining a cesium clock group time scale (TA)Cs) Data is m (t)i) Is a time tiI-1, 2,3, … at tiThe derivative of the time instant is defined as m' (t) ═ m (t)i+1)-m(ti)]/(ti+1-ti)。M(tj) Is the hydrogen clock time scale (TA)H) The derivative of which is defined as M' (t) ═ M (t)j+1)-M(tj)]/(tj+1-tj). Wherein j is 1,2,3, …. The physical meanings of M '(t) and M' (t) are consistent and are both the rates of the time scale sequence. Therefore, when t isi＝tjWhen M '(t) is M' (t).
m(MJD2)＝m(MJD1)+M′(tj)△tj (3)
Wherein MJD is a simplified julian day.
Since the time intervals of the time scales of the selected cesium clock set and the hydrogen clock set are the same, the clock difference estimation at the next moment of the cesium clock set can be expressed by the sum of the clock difference of the cesium clock set at the moment and the product of the time scale rate and the time interval of the hydrogen clock set according to the formula (3).
Suppose input observation sequence 1 is represented as y'jCorresponding time scale is xjWeight p ofj(ii) a First derivative of observation sequence 2 is expressed as
Figure GDA0002895421270000063
Corresponding time scale is
Figure GDA0002895421270000064
The weight is
Figure GDA0002895421270000065
The output sequence is then the smooth curve to be obtained, denoted yiCorresponding time scale is xi. Three quantities are defined:
1) smoothness of the curve:
Figure GDA0002895421270000066
in the formula (I), the compound is shown in the specification,
Figure GDA0002895421270000067
is unknown, based on the smoothed data, on its third derivative
Figure GDA0002895421270000068
And (6) estimating. At two points [ x ]i+1,yi+1]And [ x ]i+2,yi+2]The smooth curve between is defined as a lagrange polynomial L of third order derived from four adjacent points i, i +1, i +2, i +3i(x)：
Figure GDA0002895421270000071
The third derivative is:
Figure GDA0002895421270000072
the third derivative between each pair of data points is set to a constant, and equation (1) can be expressed as:
Figure GDA0002895421270000073
wherein
Figure GDA0002895421270000074
Figure GDA0002895421270000075
Figure GDA0002895421270000076
Figure GDA0002895421270000077
This smooth definition means that the ideal smoothing function is a quadratic function of time (the first derivative of which is a linear function).
2) Fidelity of the smoothed curve to the observed values:
Figure GDA0002895421270000078
3) fidelity of the smooth curve to the first derivative of the observed value:
Figure GDA0002895421270000079
first derivative of
Figure GDA00028954212700000710
Can be based on the smoothing function value yiAnd (4) showing.
Using lagrange polynomials L to define Si(x) Is denoted as L'i(x)：
L′i(x)＝Ai(x)yi+Bi(x)yi+1+Ci(x)yi+2+Di(x)yi+3 (10)
Wherein the content of the first and second substances,
Figure GDA0002895421270000081
Figure GDA0002895421270000082
Figure GDA0002895421270000083
Figure GDA0002895421270000084
to use the smoothing function value yiSmoothed values representing first derivatives
Figure GDA0002895421270000085
For each time instant xiThere is a certain freedom in the choice of the surrounding four points. The constraint of the following formula ensures
Figure GDA0002895421270000086
The value of (b) is on a smooth curve.
1) To calculate
Figure GDA0002895421270000087
Using the time x1，x2，x3，x4That is, the formula (11) can be obtained
Figure GDA0002895421270000088
Wherein the content of the first and second substances,
Figure GDA0002895421270000089
2) for the first half of the input data
Figure GDA00028954212700000810
I is 2,3, … N/2, using time xi-1，xi，xi+1，xi+2I.e. by
Figure GDA00028954212700000811
Wherein the content of the first and second substances,
Figure GDA00028954212700000812
3) for the second half of the input data
Figure GDA00028954212700000813
I ═ N/2+1, N/2+2, … N-1, using the time xi-2，xi-1，xi，xi+1I.e. by
Figure GDA00028954212700000814
Wherein the content of the first and second substances,
Figure GDA00028954212700000815
4) for the
Figure GDA00028954212700000816
Using the time xN-3，xN-2，xN-1，xNI.e. by
Figure GDA00028954212700000817
Wherein the content of the first and second substances,
Figure GDA00028954212700000818
then
Figure GDA00028954212700000819
Can be expressed as:
Figure GDA0002895421270000091
the V-C combined filtering method aims at determining a smoothing value yiEqualization under three different conditions, each different setting will lead to different results:
(a) the curve needs to be smoothed (minimize S)
(b) Smooth values require observations that are close to a function (minimize F)
(c) Smoothing the first derivative of the curve requires an observation that is close to the first derivative (minimizing
Figure GDA0002895421270000092
)
The least square method is adopted to minimize the above constraint conditions to adjust each parameter, and the method is expressed as follows:
Figure GDA0002895421270000093
wherein epsilon is more than or equal to 0,
Figure GDA0002895421270000094
the degree of equalization for the three conditions can be determined by selecting the values of these two parameters: smoothing coefficients ε and ∈
Figure GDA0002895421270000095
The larger the smoothing coefficient value, the larger the weight of the fidelity with respect to the observation function value or its first derivative, and the closer the smoothing value is to the observation value.
The embodiment of the invention provides a method for realizing hydrogen cesium time scale fusion by using Vondark-Cepek combined filtering. The fusion of the hydrogen cesium clock group time scale layer is put forward for the first time, and the problem that white noise generated by hydrogen cesium mutual reference influences a time scale result is avoided. The hydrogen clock group time scale and the cesium clock group time scale respectively have good short-term stability and long-term stability, so that the hydrogen clock group time scale and the cesium clock group time scale can be subjected to performance enhancement by utilizing the excellent short-term stability of the hydrogen clock group time scale under the condition of not influencing the long-term stability of the cesium clock group time scale, and a fusion time scale with good long stability and good short stability can be generated.
As shown in fig. 1, the method of generating a fusion timescale comprises the steps of:
step 1, preprocessing atomic clock difference data measured in a laboratory, and removing abnormal points by adopting a 3 sigma rule.
The gross error elimination is realized by using a 3 sigma criterion, and the detected atomic clock frequency difference data is { x1,x2,x3,…,xN-1,xNMean of samples is
Figure GDA0002895421270000101
Standard deviation of
Figure GDA0002895421270000102
The 3 sigma criterion is the measured value
Figure GDA0002895421270000103
The value is singular and should be eliminated.
Step 2, grouping the cesium atomic clocks and the hydrogen atomic clocks, and calculating time scales of the respective clock groups by using the clock difference data processed in the step 1 and adopting an exponential filtering method to obtain TACsAnd TAHAs shown in fig. 2.
Clock error estimation: clock error data X of each watch clock and reference clockij(t), the measurement interval is τ. The clock difference estimate at time (t + τ) is calculated from equation (1) based on the clock difference and clock speed at time t. The optimal estimate of the reference clock relative to the time scale clock difference is calculated from equation (2) based on a weighted average principle.
Figure GDA0002895421270000104
Figure GDA0002895421270000105
Wherein, Xi(t),Yi(t) estimates of the time difference and frequency difference, respectively, for the clock i at time t with respect to a reference time scale;
Figure GDA0002895421270000106
is the clock error forecast for clock i at time t.
Frequency estimation: the average frequency difference estimation of the clock i at the time of (t + tau) is obtained by adopting the first difference of clock difference data and then calculating through an exponential filter, such as the formulas (3) and (4).
Figure GDA0002895421270000107
Figure GDA0002895421270000108
Wherein
Figure GDA0002895421270000109
Is the frequency difference forecast of the clock i over the time interval from time t to t + τ; m isiIs an exponential frequency-averaged time constant, is the error introduced in the minimized time difference estimation equation (1), and is estimated according to equation (5).
Figure GDA00028954212700001010
Wherein tau isminiTaking the Allen variance curve as the most stable time interval of the clock i
Figure GDA00028954212700001011
The value of τ at the minimum.
And (3) weight estimation: the weight is determined by the difference between the clock difference of each clock relative to the page clock and the calculated value of the current time. As in equation (6):
Figure GDA00028954212700001012
wherein KiThe deviation for error estimation is calculated according to equation (10). The mean square time error per clock is determined by the exponential filter of equation (7).
Figure GDA00028954212700001013
Wherein epsiloni(τ) is the accumulated error of the moveout estimate of clock i over the time interval τ;<ε2 x(τ)>is the mean square error of time t on the paper surface with time interval tau, and the initial value is generally taken
Figure GDA0002895421270000111
NτIs the time constant of an exponential filter used to estimate the current mean square error, NτGenerally, it is taken for 20 days.
Figure GDA0002895421270000112
And n is the number of atomic clocks. The weights are applied in equation (2), which is calculated as shown in equation (9).
Figure GDA0002895421270000113
Figure GDA0002895421270000114
First reading the clock error file, for Xi(t)、Yi(t) taking an initial value. For the selection of the initial value, the clock of the atomic clock relative to UTC at the starting time is selected for simplicityThe difference value. The BIPM's T publication discloses the clock difference between UTC and UTC (ntsc), and it is known that the clock difference between three hydrogen clocks and UTC (ntsc) clock difference at that time are subtracted from each other to obtain the clock difference between three hydrogen clocks and UTC. Clock speed YiAnd (t) and the initial value of the drift parameter d are selected according to the clock error data of the previous ten days by secondary fitting. Secondly, the next time forecast clock error is calculated according to the formula (1)
Figure GDA0002895421270000115
Then calculate
Figure GDA0002895421270000116
Is generally taken as the initial value of
Figure GDA0002895421270000117
Then the allen variances for three clocks should be calculated first. Because the measured clock difference between each species and the main clock needs to be solved by using a triangular cap method to obtain the respective Allen variance of each clock, the known clock difference data needs to be subtracted from each other, the main clock is removed, and the clock difference between the three hydrogen atomic clocks is obtained. Then calculating the Allan variance between each two atomic clocks, solving the Allan variance of each clock according to the triangular cap calculation formula (11), and solving the Allan variance of each clock according to the formula
Figure GDA0002895421270000118
Computing
Figure GDA0002895421270000119
And (5) initial value. And calculating the initial weight value according to the formula (8-9). Then, according to the formula (2) and the data of each atomic clock, calculating the clock error of the main clock relative to the paper clock at the next moment, and carrying out weighted average to obtain the final clock error estimation Xj(t + τ) is the difference between the calculated average time at time (t + τ) and the reference clock, TACsOr TAH. Then, based on the estimated clock difference, the clock speed at the next time (t + τ) is predicted. According to Xi＝Xj+XijAnd calculating the clock difference between each clock and the paper clock, and substituting the clock difference into the formula (3) to calculate the frequency difference estimation. Then calculating the exponential frequency constant according to the formula (5), substituting the exponential frequency constant into the formula (4) to calculate Yi(t+τ)，For the calculation of the clock error prediction at the next time instant (t +2 τ). Finally, the clock difference, weight, and frequency estimation at the next time (t +2 τ) are performed. In this case, the time difference estimation is the same as the previous calculation method, the weight calculation is calculated according to equations (6) - (10), and the parameter K is determined according to equation (10)iThen, the current weighting is calculated according to the formulas (6) to (9). And circulating in sequence, and performing prediction calculation at the next moment.
Figure GDA0002895421270000121
And 3, carrying out primary difference on the hydrogen clock group time scale sequence obtained by calculation in the step 2 to obtain corresponding frequency data.
Figure GDA0002895421270000129
Step 4, determining the frequency response T/T' according to the frequency f or the period P of the noise to be weakened, and calculating the smoothing factors epsilon and epsilon
Figure GDA0002895421270000122
Time scale sequence TA of cesium clock groupCsTransforming the frequency domain to a frequency domain through Fourier transform, determining a noise period P or a frequency f through spectral analysis, simultaneously determining a frequency response T/T', and finally calculating a smoothing factor according to a formula.
Figure GDA0002895421270000123
And 5, taking the time scale of the cesium clock set generated by calculation in the step 2 and the frequency sequence of the hydrogen clock set generated in the step 3 as input, and using the smoothing factor obtained in the step 4 in a Vondark-Cepek combined filtering method to generate a time scale fusion result.
First, system equation
Figure GDA0002895421270000124
In a matrix mannerCan be expressed as
Figure GDA0002895421270000125
Wherein the content of the first and second substances,
Figure GDA0002895421270000126
Figure GDA0002895421270000127
wherein the content of the first and second substances,
Figure GDA0002895421270000128
since A is a positive definite matrix, the Cholesky decomposition can be used to solve system equation (14). Matrix a is decomposed by Cholesky to obtain:
A＝UTDU (15)
let UTThe vector w can be found By taking w as By'. Then according to Uy ═ D-1w, y can be obtained. Namely:
y＝A-1By' (16)
therefore, the time scale result of hydrogen-cesium fusion, i.e. the sequence y, can be obtained, and the result is shown in FIG. 3, wherein (a) is the time scale fusion result, and (b) is the Allan deviation.
The above embodiments are merely exemplary embodiments of the present invention, and are not intended to limit the present invention. The scope of the invention is defined by the claims. Various modifications and equivalents may be made by those skilled in the art within the spirit and scope of the invention, and such modifications and equivalents should also be considered as falling within the scope of the invention.

Claims (6)

1. A hydrogen cesium time scale fusion method based on Vondark-Cepek filtering is characterized by comprising the following steps:
step 1: grouping cesium atomic clocks and hydrogen atomic clocks, and generating exponential filtering cesium clock time scale TA through clock error prediction, frequency estimation and weight estimation stepsCsAnd hydrogen time scale TAH；
Step 2: carrying out primary difference on the hydrogen clock group time scale sequence to obtain corresponding frequency data;
and step 3: transforming the time scale sequence of the cesium clock set into a frequency domain through Fourier transform, and determining the frequency f or the period P of a signal to be suppressed through spectrum analysis; calculating a time scale smoothing factor of the cesium clock group according to the least square theory
Figure FDA0002895421260000011
And hydrogen clock group time scale smoothing factor
Figure FDA0002895421260000012
Wherein T and T' are each TACsAnd TAHA corresponding frequency response;
and 4, step 4: taking the cesium clock set time scale and the hydrogen clock set frequency sequence as input, using a smoothing factor in a Vondark-Cepek combined filtering method to generate a time scale fusion result, and realizing the fusion of the cesium clock set time scale and the hydrogen clock set time scale:
establishing a hydrogen cesium time scale fusion model, namely defining a cesium clock group time scale (TA)Cs) Data is m (t)i) Is a time tiI-1, 2,3, … at tiThe derivative of the time instant is defined as m' (t) ═ m (t)i+1)-m(ti)]/(ti+1-ti)；M(tj) Is the hydrogen clock time scale (TA)H) The derivative of which is defined as M' (t) ═ M (t)j+1)-M(tj)]/(tj+1-tj) (ii) a Wherein j is 1,2,3, …; the physical meanings of M '(t) and M' (t) are consistent, both being the rates of the time scale sequence; therefore, when t isi＝tjWhen M '(t) is M' (t);
m(MJD2)＝m(MJD1)+M′(tj)△tj (3)
wherein MJD is the simplified julian day;
because the time intervals of the time scales of the selected cesium clock group and the hydrogen clock group are the same, the clock difference estimation of the next moment of the cesium clock group with the formula (3) can be represented by the sum of the clock difference of the cesium clock group at the moment and the product of the time scale rate and the time interval of the hydrogen clock group;
suppose input observation sequence 1 is represented as y'jCorresponding time scale is xjWeight p ofj(ii) a First derivative of observation sequence 2 is expressed as
Figure FDA0002895421260000014
Corresponding time scale is
Figure FDA0002895421260000015
The weight is
Figure FDA0002895421260000016
The output sequence is then the smooth curve to be obtained, denoted yiCorresponding time scale is xi(ii) a Three quantities are defined:
1) smoothness of the curve:
Figure FDA0002895421260000013
in the formula (I), the compound is shown in the specification,
Figure FDA00028954212600000210
is unknown, based on the smoothed data, on its third derivative
Figure FDA00028954212600000211
Carrying out estimation; at two points [ x ]i+1,yi+1]And [ x ]i+2,yi+2]The smooth curve between is defined as a lagrange polynomial L of third order derived from four adjacent points i, i +1, i +2, i +3i(x)：
Figure FDA0002895421260000021
The third derivative is:
Figure FDA0002895421260000022
the third derivative between each pair of data points is set to a constant, and equation (1) can be expressed as:
Figure FDA0002895421260000023
wherein
Figure FDA0002895421260000024
Figure FDA0002895421260000025
Figure FDA0002895421260000026
Figure FDA0002895421260000027
This smooth definition means that the ideal smoothing function is a quadratic function in time, the first derivative of which is a linear function;
2) fidelity of the smoothed curve to the observed values:
Figure FDA0002895421260000028
3) fidelity of the smooth curve to the first derivative of the observed value:
Figure FDA0002895421260000029
first derivative of
Figure FDA00028954212600000212
Can be based on the smoothing function value yiAnd (4) showing.
2. The Vondark-Cepek filtering-based hydrogen cesium time scale fusion method according to claim 1, characterized in that: the method comprises the following steps that 1, atomic clock difference data measured in a laboratory are preprocessed, atomic clocks with the best long-term and short-term stability in a cesium atomic clock group and a hydrogen atomic clock group are respectively selected as reference clocks, and clock difference data of other atomic clocks in the same clock group and the reference clocks are measured by using a time interval counter.
3. The Vondark-Cepek filtering-based hydrogen cesium time scale fusion method according to claim 2, characterized in that: removing abnormal points from the clock error data obtained by preprocessing by adopting a 3 sigma rule, detecting the continuity of the clock error data, and repairing a clock error sequence with continuous missing data less than 5 measurement periods by adopting an interpolation method; and if the clock error data continuously lack for more than 5 measurement periods, eliminating the atomic clock in the time keeping clock group.
4. The Vondark-Cepek filtering-based hydrogen cesium time scale fusion method according to claim 1, characterized in that: in the step 1, if the weight of each atomic clock is greater than the set maximum weight upper limit, the weight of the atomic clock is set as the maximum weight; h 'is'i(t) is the atomic clock H at time tiTime correction quantity is weight, atomic clock data, h'i(t) by the formula
Figure FDA0002895421260000031
Calculation of where xij(t) is an atomic clock HjAnd atomic clock HiThe clock difference at time t, N being the number of atomic clocks involved, xi(t) is a laboratory single atomDifference between average time scales of clock and laboratory atomic clock set, wiIs the weight of each atomic clock.
5. The Vondark-Cepek filter-based hydrogen cesium time scale fusion method according to claim 4, characterized in that: in the step 1, the maximum weight is set to be 2.5/N.
6. The Vondark-Cepek filtering-based hydrogen cesium time scale fusion method according to claim 1, characterized in that: in the step 4, a time scale TA of the cesium clock group is definedCsData is m (t)i) Is a time tiI-1, 2,3, … at tiThe derivative of time m' (t) ═ m (t)i+1)-m(ti)]/(ti+1-ti) (ii) a Definition of M (t)j) Is hydrogen clock group time scale TAHThe derivative M' (t) thereof is [ M (t) ]j+1)-M(tj)]/(tj+1-tj) Wherein j is 1,2,3, …; when t isi＝tjWhen M ' (t) is M ' (t), M (MJD2) is M (MJD1) + M ' (t)j)△tjMJD denotes reduced julian days; suppose input observation sequence 1 is represented as y'jCorresponding time scale is xjWeight p ofj(ii) a First derivative of observation sequence 2 is expressed as
Figure FDA0002895421260000032
Corresponding time scale is
Figure FDA0002895421260000034
The weight is
Figure FDA0002895421260000033
The output sequence is then the smooth curve y to be obtainediCorresponding time scale is xiDefining the smoothness of the curve
Figure FDA0002895421260000035
In the formula (I), the compound is shown in the specification,
Figure FDA0002895421260000036
is unknown, based on the smoothed data, on its third derivative
Figure FDA0002895421260000037
Carrying out estimation; at two points [ x ]i+1,yi+1]And [ x ]i+2,yi+2]The smooth curve between is defined as a lagrange polynomial of third order derived from four adjacent points i, i +1, i +2, i +3
Figure FDA0002895421260000041
Third derivative thereof
Figure FDA0002895421260000042
Defining the fidelity of a smooth curve to an observed value
Figure FDA0002895421260000043
Defining the fidelity of the smooth curve to the first derivative of the observed value
Figure FDA0002895421260000044
CN201910787030.4A 2019-08-25 2019-08-25 Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering Expired - Fee Related CN110554597B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910787030.4A CN110554597B (en) 2019-08-25 2019-08-25 Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910787030.4A CN110554597B (en) 2019-08-25 2019-08-25 Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering

Publications (2)

Publication Number Publication Date
CN110554597A CN110554597A (en) 2019-12-10
CN110554597B true CN110554597B (en) 2021-03-26

Family

ID=68737950

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910787030.4A Expired - Fee Related CN110554597B (en) 2019-08-25 2019-08-25 Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering

Country Status (1)

Country Link
CN (1) CN110554597B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110989326B (en) * 2019-12-26 2021-03-30 中国计量科学研究院 Local high-precision time frequency real-time comprehensive device
CN112269311B (en) * 2020-09-22 2022-04-08 中国计量科学研究院 Method and device for realizing real-time atomic time scale of remote distributed union
CN112329197A (en) * 2020-09-23 2021-02-05 北京无线电计量测试研究所 Comprehensive atomic time establishing method based on gray model
CN114047684B (en) * 2021-10-21 2023-05-30 中国人民解放军61081部队 Atomic clock joint time keeping method and device

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9252795B2 (en) * 2010-12-15 2016-02-02 Raytheon Company Distribution system for optical reference
CN105676627A (en) * 2015-12-25 2016-06-15 中国科学院国家授时中心 Time keeping system primary and standby main clock seamless switching system and method
CN106506136A (en) * 2016-11-24 2017-03-15 上海市计量测试技术研究院 A kind of network time transmission method and device based on atomic clock group
CN106773610A (en) * 2016-11-28 2017-05-31 北京工业大学 A kind of cesium-beam atomic clock and hydrogen clock frequency difference predictor method

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3443534B2 (en) * 1998-12-17 2003-09-02 日本電信電話株式会社 Atomic frequency standard laser pulse oscillator
CN2672698Y (en) * 2003-12-19 2005-01-19 中国科学院国家授时中心 Novel time interval measurer
CN104579340B (en) * 2015-02-04 2018-01-30 上海航天测控通信研究所 A kind of passive hydrogen clock digital servosystem based on FPGA
CN105718642B (en) * 2016-01-18 2018-10-23 中国科学院国家授时中心 A kind of reference time scale production method based on threshold autoregressive model
CN105974777B (en) * 2016-07-19 2018-07-31 北京工业大学 It is a kind of to generate atomic time calibration method using Algos and Kalman combinations

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9252795B2 (en) * 2010-12-15 2016-02-02 Raytheon Company Distribution system for optical reference
CN105676627A (en) * 2015-12-25 2016-06-15 中国科学院国家授时中心 Time keeping system primary and standby main clock seamless switching system and method
CN106506136A (en) * 2016-11-24 2017-03-15 上海市计量测试技术研究院 A kind of network time transmission method and device based on atomic clock group
CN106773610A (en) * 2016-11-28 2017-05-31 北京工业大学 A kind of cesium-beam atomic clock and hydrogen clock frequency difference predictor method

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
一种基于小波变换的氢铯联合守时算法研究;宋会杰等;《天文学报》;20180131;第59卷(第1期);46-55 *
基于Vondark模型平滑精密单点定位（PPP）授时钟差解;管健安等;《测绘工程》;20120430;第21卷(第2期);25-28 *

Also Published As

Publication number Publication date
CN110554597A (en) 2019-12-10

Similar Documents

Publication Publication Date Title
CN110554597B (en) Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering
CN111650617B (en) Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering
CN106026919B (en) The punctual compensation method of crystal oscillator
Heo et al. Improving prediction accuracy of GPS satellite clocks with periodic variation behaviour
CN109581856B (en) Time synchronization and time keeping method based on high-performance crystal oscillator frequency calibration
CN105718642B (en) A kind of reference time scale production method based on threshold autoregressive model
CN109508510B (en) Improved Kalman filtering-based rubidium atomic clock parameter estimation algorithm
KR20110133436A (en) Sensor system
CN112433235B (en) Method, system and medium for determining time reference
Weiss et al. A study of the NBS time scale algorithm
CN107086901A (en) A kind of BDT method for building up and UTC (NTSC) method for building up
US11035902B2 (en) Advanced fuel gauge
CN111711446B (en) Method, system and medium for taming crystal oscillator frequency by using GPS signal
CN110865530A (en) Atomic time calculation method
CN115904000A (en) Real-time clock crystal oscillator compensation method and system based on orthogonal least square method curve fitting
Krasheninnikov et al. IRI-2001 model efficiency in ionospheric radiowave propagation forecasting
CN110858309B (en) Multi-reference time clock weighting synthesis method
Lepek Clock prediction and characterization
Griffin et al. Drift correction for active hydrogen MASERs
CN112329197A (en) Comprehensive atomic time establishing method based on gray model
Wang et al. Enhanced stability for local atomic clock ensemble time scale using weighted moving average filter
CN111865267A (en) Temperature measurement data prediction method and device
Xue et al. An enhanced prediction model for BDS ultra-rapid clock offset that combines singular spectrum analysis, robust estimation and gray model
Tseng et al. Introduction of GNSS Up-Sampled All-in-View Time Transfers
CN109873640A (en) Oscillator frequency control method and device

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210326