CN111565084B - Satellite time service time keeping system and method based on frequency estimation - Google Patents

Abstract

The satellite time service and time keeping system based on frequency estimation comprises a main control computer, a satellite receiver, a constant temperature crystal oscillator and an FPGA module. The FPGA module consists of a frequency multiplier, a continuous time counter, a 1PPS signal time interval measurer and a local 1PPS signal generator, and the modules are all realized through a programmable logic unit in the FPGA. A working method of the satellite time service time keeping system based on frequency estimation is also provided. The system and the method of the invention receive the 1PPS second pulse signal from the satellite receiver, filter and eliminate non-points, fill up missing points and reduce random fluctuation, and then output the stable 1PPS second pulse signal; the frequency of a clock signal output by a constant-temperature crystal oscillator is measured and estimated, and when the 1PPS signal output by a satellite receiver is in abnormal states such as interruption, translation, out-of-tolerance and the like, the system utilizes the estimated frequency of the clock signal to keep time so as to ensure the stability of the local 1PPS signal, thereby realizing the method. The invention has the advantages of strong universality, simple circuit structure, high locking speed of 1PPS signals and the like.

Description

Satellite time service time keeping system and method based on frequency estimation

Technical Field

The invention relates to the technical field of communication, in particular to a satellite time service time keeping system and method based on frequency estimation.

Background

The essence of using a satellite for time service and keeping time is the process of taming a clock. By disciplining the clock, the time service and time keeping system can realize the following of the satellite time and output UTC time (coordinated universal time). The methods generally used are: the method comprises the following steps of (1) carrying out a crystal oscillator frequency calibration method, namely comparing a fixed frequency signal provided by a satellite time service receiver with an oscillation signal generated by a local crystal oscillator to obtain a frequency difference; and then, the oscillation frequency is basically consistent with the satellite oscillation frequency by adjusting the local crystal oscillator (the satellite time service principle and application, 1 st edition, Yangjun, single celebration, national defense industry press, 2013, 8 th page 129), so that a pulse per second signal (1PPS signal) is generated, and time synchronization is realized.

The satellite timing system using the method generally comprises a satellite receiver, a main control chip (usually an FPGA), a crystal oscillator, a D/a converter, a clock interface circuit and the like (satellite timing principle and application, 1 st edition yangjun, single dawn national defense industry press 2013, page 143).

Disclosure of Invention

Aiming at the problems in the prior art, the invention provides a satellite time service and keeping-in-time system based on frequency estimation, which comprises a main control computer, a satellite receiver, a constant-temperature crystal oscillator and an FPGA module; wherein

A satellite receiver for processing and converting a satellite signal received by the antenna; the acquisition, tracking, carrier recovery and demodulation of satellite signals are realized; outputting '1 PPS satellite signals' and time service telegraph text information to the FPGA module;

the constant-temperature crystal oscillator provides a clock reference for the FPGA module;

the FPGA module is used for realizing integration of a frequency multiplier and a continuous time counter to form a local continuous time standard; measuring the time interval of the 1PPS signal sent by the satellite receiver; outputting the filtered local 1PPS signal, i.e., "1 PPS filtering";

the main control computer is used for providing calculation and display control services for the whole system and realizing a time service and time keeping system filtering algorithm; the main control computer acquires a '1 PPS time scale' output by a 1PPS signal time interval measurer in the FPGA module, wherein the '1 PPS time scale' is correspondingly counted by a local time scale when a rising edge of a '1 PPS satellite signal' arrives; according to the characteristic that the mean error of the ' 1PPS ' pulse is zero in a long steady state, the main control computer carries out statistical filtering processing on time scale information of each second, estimates the error size of the central frequency and the nominal frequency of the constant-temperature crystal oscillator, predicts the output position of the next ' 1PPS ' second pulse, forms new local 1PPS predicted second data ' (n +1) second data), and outputs the new local 1PPS predicted second data to a local 1PPS signal generator in the FPGA module; when the satellite signal is abnormal, the main control computer solidifies the parameters of each register, enters a time-keeping state and keeps stable output of the '1 PPS filtering' signal.

In one embodiment of the invention, the FPGA module consists of a frequency multiplier, a continuous time counter, a 1PPS signal time interval measurer and a local 1PPS signal generator, which are all realized by a programmable logic unit in the FPGA, wherein

The frequency multiplier is a digital phase-locked loop which multiplies the 10MHz frequency output by the constant-temperature crystal oscillator into a high-frequency clock signal CLK _HR And outputting;

continuous time counter for CLK from frequency multiplier _HR Calculating the continuous time of the frequency signal to form a local continuous time reference; to reduce the timing step, the multiplier output CLK _HR 0 and CLK _HR Two-way _piclock signals with 180 deg phase difference and two independent continuous time counters for counting to form local time mark 1 and local time mark 2 with minimum step length of 1/CLK _HR The joint step size is 1/2CLK _HR (ii) a The time scale 1 is divided into two paths, one path is output to a 1PPS signal time interval measurer, and the other path is output to a local 1PPS signal generator; the time scale 2 is output to a 1PPS signal time interval measurer;

the 1PPS signal time interval measurer measures the leading edge of a 1PPS satellite signal from the satellite receiver to form 1PPS time mark information and sends the information to the main control computer;

a local 1PPS signal generator that receives the "time stamp 1" information from the continuous-time counter and the "(n +1) seconds of data" information from the host computer, and forms a "1 PPS filtered" signal when the "time stamp 1" is equal to the "(n +1) seconds of data".

In an embodiment of the present invention, the local time scale 1 and the local time scale 2 have a minimum step size of 5ns for each time scale and a joint step size of 2.5 ns.

The working method of the satellite time service time keeping system based on the frequency estimation is also provided, and comprises the following steps:

step 1: connection system

Connecting parts of a system according to any of claims 1 to 3;

and 2, step: initialization

After the system is powered on, the system initialization is carried out, and the step is subdivided into 5 steps:

(a2) through output

After power-on and before convergence of a filtering algorithm, in order to quickly give UTC time, a direct connection mode is adopted, and a 1PPS (pulse per second) pulse signal sent by a satellite receiver is directly used as a '1 PPS filtering' signal to be output;

(b2) parameter initialization

The time service filter algorithm model is as follows:

x(n|n)＝x(n|n-1)+α[y(n)-x(n|n-1)]

wherein x (n | n) is the filtered value of the nth beat, x (n | n-1) is the predicted value of the nth beat relative to the nth beat,

the filtered value of the nth beat frequency, y (n) is the measured value of the nth beat, and x (n +1| n) is the predicted value of the nth period to the n +1 th period; alpha and beta are respectively the weight coefficients of the time difference and the frequency difference of the signal of the 1PPS filtering;

T ₀ the time filtering value is as follows: x (0|0) ═ y (0);

T ₀ time frequency filtering value:

T ₁ moment predicted value:

(c)T ₀ validation of the measured value y (0) at a time

The method specifically comprises the following steps:

1) the satellite receiver sends out a 'normal time service' message;

2) after the constant temperature crystal oscillator is stabilized, continuously taking 3 beats of measured values, and defining as: y (-3), y (-2), y (-1);

3) the filtering algorithm is started in advance, and values of x (-2| -3), x (-1| -2) and x (0| -1) are extrapolated according to the algorithm model in the step (b) in the step 2;

4) y (0) valid decision, e.g. T _-3 、T _-2 、T _-1 Taking the 4 th observation value as y (0) when no non-point appears at any moment;

note that: firstly, a non-point judgment method is shown in step 4; during the 'y (0) effective judgment', the 1PPS pulse keeps direct output until the filtering is normally started, and the 1PPS predicts the value according to the filtering algorithm; considering that the crystal oscillator is not sufficiently stable during initialization, the frequency filtering value of the crystal oscillator may not be accurate enough, so that the non-point judgment condition is determined in y (0): "| y (n) | ≧ 3 σ/τ" is relaxed to "y (n) | (n) -x (n | n-1) | ≧ 5 σ/τ"), σ is the random error of the satellite receiver 1PPS pulse, τ is the quantization precision of y (n);

and step 3: reading in the measured values y (n)

Reading a local time scale corresponding to the rising edge of a 1PPS pulse signal output by a satellite receiver; this step is subdivided into 2 steps:

(a3) y (n) calculating weights

The system measures 1PPS pulse by two clocks CLK with 180 DEG phase difference _HR Counting, in order to prevent the continuous time counter from overflowing, the word width of the measured value information is more than or equal to 48 bits, and two groups of measured values are respectively output by two time scales: y is ₁ (n)，y ₂ (n)；

When y ₁ (n)-y ₂ When (n) | is less than or equal to 1, y (n) ═ y ₁ (n)+y ₂ (n)]/2；

When y ₁ (n)-y ₂ When the (n) | > 1 is used,

(b3) non-point culling

According to the pulse precision output by the PPS of the north fighter 1, when | y (n) -x (n | n-1) | ≧ 3 σ/tau, y (n) is abandoned, x (n | n) | x (n | n-1),

and 4, step 4: updating the n-th beat filter value x (n | n)

Obtaining a local time scale x (n | n) of the smoothed second pulse signal of the system 1PPS through a measurement value y (n) of the satellite 1PPS pulse signal according to the algorithm model in the step (b) in the step 2;

and 5: updating the predicted value x (n +1| n) of the nth period to the (n +1) th period

According to the algorithm model in step2 (b), the n-th beat filter value x (n | n) and the n-th beat filter value of the crystal oscillator frequency are obtained

Calculating a predicted value x (n +1| n) of the nth period to the (n +1) th period;

step 6: updating estimated value f (n) of crystal oscillator frequency

According to the measured value y (n) of the satellite 1PPS pulse signal, a unary linear regression model is adopted to solve f (n):

wherein:

f (n) is the estimated value of the nth beat crystal oscillator frequency;

x _i for the ith valid 1PPS pulse number, i.e. x _i I, i is 1, 2, 3, 4, … …, n, if there is a k-th beat non-dot, the serial number x _k Discard, i.e. x _i The sequence becomes: … … x _k-2 ，x _k-1 ，x _k+1 ，x _k+2 ……；

size (x) isThe total number of sequence numbers, if there are j non-dots, the total number of sequence numbers size (x) n-j,

the serial numbers corresponding to j non-points are also removed;

y _i an observed value y (i) of 1PPS pulse of the ith beat; if the k-th beat is not dot, y _k Abandon, at the same time, the corresponding sequence number x _k Are also discarded, i.e. if y ₅ When it is not a point, x _i The sequence becomes: 1, 2, 3, 4, 6, 7 … …; size (x) n-1;

size (y) is the total number of observed points, size (y) size (x);

note that:

(ii) involving in calculations

Are all valid y (n) values;

secondly, in the timekeeping state f (n), the calculation is suspended, and after the timekeeping is quitted, the historical point can still be normally used and the calculation is continued;

if limited by the calculation scale, when y (n) needs to be sampled, the method can be carried out according to the following steps:

when the calculation scale is k points and the diffusion times are q times, the diffusion completion duration r is:

the method comprises the following specific steps:

step 1: defining the first valid observation as y ₁ Y (1) with the number x ₁ ＝x(1)＝1，i＝1；

Step 2: when i is less than or equal to k, calculating y (n) in real time according to the actual points;

step 3: when k is more than i and less than or equal to 2k, the sliding window calculates y (n);

step 4: when i is more than 2k and less than or equal to 3k, starting the 1 st diffusion, wherein the diffusion rule is as follows:

when i is 2k +1, taking k +1, k +3, k +4, k +5, k +6, ·, 2k and 2k + 1;

when i is 2k +2, k +1, k +3, k +5, k +6, · 2k, 2k +1, 2k +2 are taken;

thirdly, the rest is done by analogy;

when i is 3k-1 and i is 3k, finishing the 1 st diffusion of k +1, k +3, 3k-3 and 3 k-1;

step 5: i is more than 3k and less than or equal to 4k, starting the 2 nd diffusion, wherein the diffusion rule is as follows:

when i is 3k +1 and i is 3k +2, k +1, k +4, k +7, k +9, k +11, k +13, ·, 3k-1 and 3k +1 are taken;

when i is 3k +3 and i is 3k +4, k +1, k +4, k +7, k +10, k +13,. loganj, 3k-1, 3k +1 and 3k +3 are taken;

thirdly, the rest is done by analogy;

when i is 4K-1 and i is 4K, finishing the 2 nd diffusion of K +1, K +4, K +7, 4K-5 and 4K-2;

step 6: according to the method described in the steps 4-5, when i is more than 4k and less than or equal to 5k, beginning to perform diffusion for the 3 rd time;

step 7: according to the method described in the steps 4-5, when i is more than 5k and less than or equal to 6k, starting to perform the 4 th diffusion;

step 8: by analogy, when (q +1) k is less than i and less than or equal to (q +2) k, performing the q-th diffusion;

step 9: when i is equal to (q +2) k, q times of diffusion are completed, equal-interval sampling with the dot pitch of q +1 is formed, and when i is equal to or more than (q +2) k, y (n) is calculated by a sliding window according to the sampling state after q times of diffusion, namely, the fixed dot pitch is q + 1;

the following points are explained for the above sampling algorithm:

firstly, if the algorithm enters the time-keeping state in the process, the algorithm is suspended, and after the algorithm exits the time-keeping state, the effective y can be continuously used _n Participate in y (n) calculation, note: number x _i Always corresponding to absolute time;

secondly, the 1 st diffusion is carried out after 1h, and the stability of the crystal oscillator reaches a nominal value at the moment, so that the unstable period of the crystal oscillator is passed;

after the 1 st diffusion starts, all data participate in the calculation of y (n) until the 14 th diffusion is completed, so that the sampling logic of y (n) is strictly checked, and the influence of non-points and the like on the calculation precision of y (n) is avoided;

and 7: time keeping state

When the 1PPS pulse signal output by the satellite receiver is unstable or interrupted, the system enters a timekeeping state; the timekeeping state can be subdivided into 3 steps as follows:

(a7) time-keeping state parameter assignment

x(n|n)＝x(n|n-1)，

(b7) Entry condition of time keeping state

The state of keeping watch is entered when one of the following conditions is satisfied:

firstly, the satellite receiver sends out a message of abnormal time service;

secondly, when non-points appear in 3 continuous beats, entering a time keeping mode;

(c7) principle of use of f (n)

1) F (n) updating values by using factory binding frequency or f (n) measured for 24h continuously after the system is started and before the temperature of the constant-temperature crystal oscillator is stable and meets the time-keeping condition;

2) after the system is started, if the valid y (n) is greater than 3600 points and the time keeping condition is met, the current frequency estimation value f (n) is used for time keeping, otherwise, the value of Y (n) of factory binding or the value of Y (n) measured for 24h continuously is still used for updating.

In a specific embodiment of the invention, α, β are 0.4 and 0.1, respectively.

In another embodiment of the present invention, k is 1800, q is 14, and r is 8 h.

The system and the method estimate the constant-temperature crystal oscillator frequency through the 1PPS pulse output by the satellite receiver, and then directly form a local second pulse signal according to the estimated frequency value, thereby realizing UCT time tracking and synchronization and carrying out time service and time-keeping filtering on satellite time. The system and the method do not need to use a D/A converter to acclimate the constant-temperature crystal oscillator, the crystal oscillator can quickly enter a stable period, the locking of 1PPS signals can be completed within about 30s, and the locking speed of the system and the method is more than one order of magnitude higher than that of a crystal oscillator frequency calibration method.

After receiving the 1PPS second pulse signal from the satellite receiver, the system and the method filter the pulse signal, eliminate non-points, fill up missing points and reduce random fluctuation, and output a stable 1PPS second pulse signal; and measuring and estimating the frequency of a clock signal output by the constant-temperature crystal oscillator, and when the 1PPS signal output by the satellite receiver is in abnormal states such as interruption, translation, out-of-tolerance and the like, the system utilizes the estimated frequency of the clock signal to keep time so as to ensure the stability of the local 1PPS signal. The invention has the advantages of strong universality, simple circuit structure, high locking speed of 1PPS signals and the like.

Drawings

FIG. 1 is a schematic diagram of a satellite time service time keeping system for implementing frequency estimation;

FIG. 2 shows a satellite time service time keeping system method work flow based on frequency estimation.

Detailed description of the preferred embodiment

The satellite time service and keeping-in-time system based on frequency estimation according to the present invention will be described in detail with reference to the accompanying drawings.

As shown in fig. 1, the invention provides a satellite time service and keeping-in-time system based on frequency estimation, which comprises a main control computer, a satellite receiver, a constant temperature crystal oscillator and an FPGA module, and is described in detail as follows.

The satellite receiver processes and transforms the satellite signals received by the antenna; the acquisition, tracking, carrier recovery and demodulation of satellite signals are realized; and outputting information such as a 1PPS second pulse signal (hereinafter, referred to as a "1 PPS satellite signal") and a time service message to the FPGA module.

The constant-temperature crystal oscillator provides a clock reference for the FPGA module, the common output frequency is 10MHz, and the frequency stability is 0.1-10 ppb.

The FPGA module is used for realizing integration of a frequency multiplier and a continuous time counter to form a local continuous time standard; measuring the time interval of the 1PPS signal sent by the satellite receiver; a filtered local 1PPS signal, referred to as "1 PPS filtering" for short, is output.

In the system of the invention, the FPGA module consists of a frequency multiplier, a continuous time counter, a 1PPS signal time interval measurer and a local 1PPS signal generator, which are all realized by a programmable logic unit in the FPGA, wherein

The frequency multiplier is a digital phase-locked loop which multiplies the 10MHz frequency output by the constant-temperature crystal oscillator into a high-frequency clock signal CLK _HR And outputs, in this example, CLK _HR The oscillation frequency of (2) is 200MHz, and other suitable clock frequencies can be set according to the design requirements of the system.

Continuous time counter for CLK from frequency multiplier _HR And carrying out continuous time calculation on the frequency signal to form a local continuous time reference. To reduce the timing step, the multiplier output CLK _HR 0 and CLK _HR Two-way _piclock signals with 180 deg phase difference and two independent continuous time counters for counting to form local time mark 1 and local time mark 2 with minimum step length of 1/CLK _HR (e.g., 5ns) with a joint step size of 1/2CLK _HR (e.g., 2.5 ns). The local time scale is generated by a continuous time counter (the counting clock of the continuous time counter is a constant temperature crystal oscillator), and the stability is high. The time scale 1 is divided into two paths, one path is output to a 1PPS signal time interval measurer, and the other path is output to a local 1PPS signal generator; the time scale 2 is output to the 1PPS signal time interval measurer.

The 1PPS signal time interval measurer measures the leading edge of a 1PPS satellite signal from the satellite receiver to form 1PPS time mark information and sends the information to the main control computer.

And the local 1PPS signal generator receives the time scale 1 information from the continuous time counter and the (n +1) second data information from the main control computer, and forms a 1PPS filtering signal when the time scale 1 is equal to the (n +1) second data.

And the main control computer is used for providing calculation and display control services for the whole system and realizing a time service and time keeping system filtering algorithm. In consideration of the requirements of system real-time performance, robustness and realizability of programming, the time service filter adopts an alpha-beta filtering algorithm (an improved alpha-beta filtering algorithm, New year modern electronic technology, No. 21 of 2012) and the time service filter adopts a unary linear regression algorithm (application regression analysis, 5 th edition of Hexuan swarm, Liuwenqing, China university Press, 2019, No. 15-24). And the main control computer acquires a 1PPS time scale output by the 1PPS signal time interval measurer in the FPGA module, wherein the 1PPS time scale is correspondingly counted by a local time scale when the rising edge of the 1PPS satellite signal arrives. According to the characteristic that the mean error of the '1 PPS' pulse is zero in a long steady state, the main control computer carries out statistical filtering processing on time scale information of each second, estimates the error between the central frequency and the nominal frequency of the constant-temperature crystal oscillator, predicts the output position of the next '1 PPS' second pulse, forms new local 1PPS predicted second data (n +1) second data), and outputs the new local 1PPS predicted second data to a local 1PPS signal generator in the FPGA module. When the satellite signal is abnormal, the main control computer solidifies the parameters of each register, enters a time-keeping state and keeps stable output of the '1 PPS filtering' signal. Due to the high stability of the frequency source, the time keeping precision of the local time can be effectively ensured.

The satellite time service timekeeping system carries out filtering optimization on a pulse per second signal (1PPS satellite) output by a satellite receiver through a high-stability frequency source (constant-temperature crystal oscillator), and outputs a local pulse per second signal (1PPS filtering) with stability meeting the system index requirement through a main control computer and an FPGA module after leak points are filled up, non-points are removed and phase noise is filtered. Meanwhile, the central frequency of the constant-temperature crystal oscillator is measured and estimated through long-term statistics of the '1 PPS satellite' signals, so that when the satellite signals of the system are abnormal, the time-keeping module continues to normally output the pulse-per-second signals, the output stability of the '1 PPS filtering' signals is kept, and the time-keeping precision requirement is met.

As shown in fig. 2, the working method of the satellite time service and keeping in time system based on frequency estimation of the present invention comprises the following steps:

step 1: connection system

Connecting the parts of the system of the invention according to the relationship shown in figure 1;

step 2: initialization

After the system is powered on, the system initialization is performed, and the step can be subdivided into 5 steps:

(a) through output

After power-on and before convergence of a filtering algorithm, in order to quickly give UTC time, a direct connection mode is adopted, and a 1PPS (pulse per second) pulse signal sent by a satellite receiver is directly used as a '1 PPS filtering' signal to be output;

(b) parameter initialization

The time service filter algorithm model is as follows:

x(n|n)＝x(n|n-1)+α[y(n)-x(n|n-1)]

wherein x (n | n) is the filtered value of the nth beat, x (n | n-1) is the predicted value of the nth beat relative to the nth beat,

and (d) the filtered value of the nth beat frequency, y (n) is the measured value of the nth beat, and x (n +1| n) is the predicted value of the nth period to the n +1 th period. Alpha and beta are weight coefficients of the time difference and the frequency difference of the signal filtered by the 1PPS respectively, and the recommended values are 0.4 and 0.1.

T ₀ The time filtering value is as follows: x (0|0) ═ y (0);

T ₀ time frequency filtering value:

T ₁ moment predicted value:

(c)T ₀ validation of the measured value y (0) at a time

1) The satellite receiver sends a 'normal time service' message;

2) after the constant temperature crystal oscillator is stabilized (generally 180-600 s), continuously taking 3 beats of measured values, and defining as: y (-3), y (-2), y (-1);

3) the filtering algorithm is started in advance, and values of x (-2| -3), x (-1| -2) and x (0| -1) are extrapolated according to the algorithm model in the step (b) in the step 2;

4) y (0) valid decision, e.g. T _-3 、T _-2 、T _-1 No non-point appeared at any time, and the observed value of the 4 th beat was taken as y (0).

(note: a non-point judgment method is shown in a step 4; during the period of ' y (0) effective judgment ', 1PPS pulse keeps direct output all the time, and the value of 1PPS is predicted according to a filter algorithm until the filtering is normally started; and thirdly, considering that the crystal oscillator is not sufficiently stable and the filtering value of the crystal oscillator frequency is not accurate enough during initialization, in the judgment of y (0), the non-point judgment condition is that | y (n) < n > -x (n | n-1) | > 3 σ/τ ' can be widened to ' y (n) < x (n | n-1) | > 5 σ/τ '), wherein σ is the random error of the satellite receiver 1PPS pulse, τ is the quantization precision of y (n), and in the example, the joint step length is 2.5 ns.

And step 3: reading in the measured values y (n)

And reading a local time scale corresponding to the rising edge of the 1PPS pulse signal output by the satellite receiver. This step can be subdivided into 2 steps:

(a) y (n) calculating weights

The system measures 1PPS pulse by two clocks CLK with 180 DEG phase difference _HR Counted as CLK _HR When the clock frequency is 200MHz, theoretically, the measurement accuracy is ± 2.5ns, in order to prevent the overflow of the continuous time counter, the width of the measured value information word should be as wide as possible, preferably equal to or greater than 48 bits (if 48 bits are selected, namely 390h unidirectional cycle), and two groups of measured values are respectively output by two time scales: y is ₁ (n)，y ₂ (n)；

When y ₁ (n)-y ₂ When (n) | is less than or equal to 1, y (n) ═ y ₁ (n)+y ₂ (n)]/2；

When y ₁ (n)-y ₂ When the (n) | > 1,

(b) non-point culling

According to the pulse precision output by the PPS of the north fighter 1, when | y (n) -x (n | n-1) | ≧ 3 σ/tau, y (n) is abandoned, x (n | n) | x (n | n-1),

and 4, step 4: updating the n-th beat filter value x (n | n)

And (c) acquiring a local time scale x (n | n) of the smoothed second pulse signal of the system 1PPS through the measurement value y (n) of the satellite 1PPS pulse signal according to the algorithm model in the step (b) in the step 2.

And 5: updating the predicted value x (n +1| n) of the nth period to the (n +1) th period

According to the algorithm model in step2 (b), the n-th beat filter value x (n | n) and the n-th beat filter value of the crystal oscillator frequency are obtained

The nth cycle to nth +1 cycle prediction value x (n +1| n) is calculated.

Step 6: updating estimated value f (n) of crystal oscillator frequency

According to the measured value y (n) of the satellite 1PPS pulse signal, a unary linear regression model is adopted to solve f (n):

wherein:

f (n) is the estimated value of the nth beat crystal oscillator frequency;

x _i for the ith valid 1PPS pulse number, i.e. x _i I, i is 1, 2, 3, 4, … …, n, if there is a k-th beat non-dot, the serial number x _k Discard, i.e. x _i The sequence becomes: … … x _k-2 ，x _k-1 ，x _k+1 ，x _k+2 ……；

size (x) is the total number of sequence numbers, if there are j non-dots, the sequenceNumber size (x) n-j,

the serial numbers corresponding to j non-points are also removed;

y _i observed value y (i) of 1PPS pulse at the ith beat; if the k-th beat is not dot, y _k Abandon, at the same time, the corresponding sequence number x _k Are also discarded, i.e. if y ₅ When it is not a point, x _i The sequence becomes: 1, 2, 3, 4, 6, 7 … …; size (x) n-1;

size (y) is the total number of observed points, size (y) size (x);

note that:

y (n) participating in calculation is all effective y (n) values;

secondly, in the timekeeping state f (n), the calculation is suspended, and after the timekeeping is quitted, the historical point can still be normally used and the calculation is continued;

if limited by the calculation scale, when y (n) needs to be sampled, the method can be carried out according to the following steps:

in the following, the calculation scale k is 1800 points, the number of diffusion times is 14, and the diffusion completion duration is 8 h.

1. Defining the first valid observation as y ₁ Y (1) with the number x ₁ ＝x(1)＝1，i＝1；

2. When i is less than or equal to k, calculating y (n) in real time according to the actual points;

3. when k is more than i and less than or equal to 2k, the sliding window calculates y (n);

4. when i is more than 2k and less than or equal to 3k, starting the 1 st diffusion, wherein the diffusion rule is as follows:

when i is 2k +1, taking k +1, k +3, k +4, k +5, k +6, ·, 2k and 2k + 1;

when i is 2k +2, k +1, k +3, k +5, k +6, · 2k, 2k +1, 2k +2 are taken;

thirdly, the rest is done by analogy;

when i is 3k-1 and i is 3k, finishing the 1 st diffusion of k +1, k +3, 3k-3 and 3 k-1;

5. i is more than 3k and less than or equal to 4k, starting the 2 nd diffusion, wherein the diffusion rule is as follows:

when i is 3k +1 and i is 3k +2, k +1, k +4, k +7, k +9, k +11, k +13, ·, 3k-1 and 3k +1 are taken;

when i is 3k +3 and i is 3k +4, k +1, k +4, k +7, k +10, k +13,. loganj, 3k-1, 3k +1 and 3k +3 are taken;

thirdly, the rest is done by analogy;

when i is 4K-1 and i is 4K, finishing the 2 nd diffusion of K +1, K +4, K +7, 4K-5 and 4K-2;

6. according to the method described in the steps 4-5, when i is more than 4k and less than or equal to 5k, beginning to perform diffusion for the 3 rd time;

7. according to the method described in the steps 4-5, when i is more than 5k and less than or equal to 6k, starting to perform the 4 th diffusion;

8. by analogy, when i is more than 15k and less than or equal to 16k, performing the 14 th diffusion;

9. when i is equal to 16k, the 14 th diffusion is completed, equal-interval sampling with the dot pitch of 15 is formed, and when i is equal to or more than 16k, y (n) is calculated according to a sliding window of the 14 th diffused sampling state (namely, the fixed dot pitch is 15);

the following points are explained for the above sampling algorithm:

firstly, if the algorithm enters a time-keeping state in the process, the algorithm is suspended, and after the algorithm exits the time-keeping state, the effective y can be continuously used _n Participate in y (n) calculation, note: number x _i Always corresponding to absolute time;

secondly, the 1 st diffusion is carried out after 1h, and the stability of the crystal oscillator reaches a nominal value at the moment, so that the unstable period of the crystal oscillator is passed;

and thirdly, after the 1 st diffusion is started, all data participate in the calculation of y (n) until the 14 th diffusion is completed, so that the sampling logic of y (n) is strictly checked, and the influence of non-points and the like on the calculation precision of y (n) is avoided.

And 7: time keeping state

When the 1PPS pulse signal output by the satellite receiver is unstable or interrupted, the system enters a timekeeping state. The timekeeping state can be subdivided into 3 steps as follows:

(a) time-keeping state parameter assignment

x(n|n)＝x(n|n-1)，

(b) Entry condition of time keeping state

The state of keeping watch is entered when one of the following conditions is satisfied:

firstly, the satellite receiver sends out a message of abnormal time service;

and secondly, when non-points appear in 3 continuous beats, entering a time keeping mode.

(c) Usage principle of f (n)

1) F (n) updating values by using factory binding frequency or f (n) measured for 24h continuously when the temperature of the constant-temperature crystal oscillator is stable after the system is started and meets the time-keeping condition.

2) After the system is started, if the valid y (n) is greater than 3600 points and the time keeping condition is met, the current frequency estimation value f (n) is used for time keeping, otherwise, the value of Y (n) of factory binding or the value of Y (n) measured for 24h continuously is still used for updating.

The method of the invention has the following advantages:

the system and the method of the invention receive the 1PPS second pulse signal from the satellite receiver, filter and eliminate non-points, fill up missing points and reduce random fluctuation, and then output a stable 1PPS second pulse signal; the frequency of a clock signal output by a constant-temperature crystal oscillator is measured and estimated, and when the 1PPS signal output by a satellite receiver is in abnormal states such as interruption, translation, out-of-tolerance and the like, the system utilizes the estimated frequency of the clock signal to keep time so as to ensure the stability of the local 1PPS signal, thereby realizing the method. The method has the advantages of strong universality, no need of using a D/A converter to taminate the constant-temperature crystal oscillator, high locking speed of the 1PPS signal and the like.

1. High versatility

The system can be used for time filtering of a satellite time service system, can also be used for time filtering and smoothing of other time service systems, and has strong universality.

2. Simple circuit structure

Because the time filtering and the prediction are carried out by adopting the frequency estimation method, the crystal oscillator frequency does not need to be domesticated by using a D/A converter, and the circuit structure is simpler than that of the traditional time service system.

3. 1PPS signal locking speed is fast

When the traditional time service system tamines the crystal oscillator frequency, the output voltage of the D/A converter needs to be continuously adjusted, the constant-temperature crystal oscillator has long stabilization time, and the locking speed of 1PPS signals is slow. By adopting the system, after the satellite receiver 1PPS signal is output, the system algorithm can quickly complete the locking of the 1PPS signal (within 30 s) by utilizing the characteristic of high short-term stability of the crystal oscillator; by utilizing the characteristic of high long-term stability of the satellite 1PPS signal, the frequency estimation is carried out by the system algorithm, and the time-keeping precision can be ensured to meet the requirement.

Claims (5)

1. A satellite time service time keeping system based on frequency estimation is characterized by comprising a main control computer, a satellite receiver, a constant temperature crystal oscillator and an FPGA module; wherein

A satellite receiver for processing and converting a satellite signal received by the antenna; the acquisition, tracking, carrier recovery and demodulation of satellite signals are realized; outputting '1 PPS satellite signals' and time service telegraph text information to the FPGA module;

the constant-temperature crystal oscillator provides a clock reference for the FPGA module;

the FPGA module is used for realizing integration of a frequency multiplier and a continuous time counter to form a local continuous time standard; measuring the time interval of the 1PPS signal sent by the satellite receiver; outputting the filtered local 1PPS signal, i.e., "1 PPS filtering"; and wherein

The FPGA module consists of a frequency multiplier, a continuous time counter, a 1PPS signal time interval measurer and a local 1PPS signal generator, and the modules are realized by a programmable logic unit in the FPGA, wherein

The frequency multiplier is a digital phase-locked loop which multiplies the 10MHz frequency output by the constant-temperature crystal oscillator into a high-frequency clock signal CLK _HR And outputting;

continuous time counter for CLK from frequency multiplier _HR Calculating the continuous time of the frequency signal to form a local continuous time reference; to reduce the time scale step, the frequency multiplier output CLK _{HR_0} And CLK _{HR_pi} Two paths of clock signals, the phase difference of the two signals is 180 degrees, two independent continuous time counters are respectively controlled to count to form a local time scale 1 and a local time scale 2, and the minimum step length of each path of time scale is 1/CLK _HR The joint step size is 1/2CLK _HR (ii) a The time scale 1 is divided into two paths, one path is output to a 1PPS signal time interval measurer, and the other path is output to a local 1PPS signal generator; the time scale 2 is output to a 1PPS signal time interval measurer;

the 1PPS signal time interval measurer measures the leading edge of a 1PPS satellite signal from the satellite receiver to form 1PPS time mark information and sends the information to the main control computer;

a local 1PPS Signal Generator that receives "time Scale 1" information from the continuous time counter and "(n +1) seconds data" information from the host computer, and that forms a "1 PPS filtered" signal when the "time Scale 1" is equal to the "(n +1) seconds data";

the main control computer is used for providing calculation and display control services for the whole system and realizing a time service and time keeping system filtering algorithm; the main control computer acquires a '1 PPS time scale' output by a 1PPS signal time interval measurer in the FPGA module, wherein the '1 PPS time scale' is correspondingly counted by a local time scale when a rising edge of a '1 PPS satellite signal' arrives; according to the characteristic that the mean error of the ' 1PPS ' pulse is zero in a long steady state, the main control computer carries out statistical filtering processing on time scale information of each second, estimates the error size of the central frequency and the nominal frequency of the constant-temperature crystal oscillator, predicts the output position of the next ' 1PPS ' second pulse, forms new local 1PPS predicted second data ' (n +1) second data), and outputs the new local 1PPS predicted second data to a local 1PPS signal generator in the FPGA module; when the satellite signal is abnormal, the main control computer solidifies the parameters of each register, enters a time-keeping state and keeps stable output of the '1 PPS filtering' signal.

2. The frequency estimation based satellite time service time keeping system of claim 1, wherein the local time scale 1 and the local time scale 2 have a minimum step size of 5ns per path and a joint step size of 2.5 ns.

3. A working method of a satellite time service time keeping system based on frequency estimation is characterized by comprising the following steps:

step 1: connection system

Connecting parts of a system according to any of claims 1 to 2;

step 2: initialization

After the system is powered on, the system initialization is carried out, and the step is subdivided into 5 steps:

(a2) through output

After power-on and before convergence of a filtering algorithm, in order to quickly give UTC time, a direct connection mode is adopted, and a 1PPS (pulse per second) pulse signal sent by a satellite receiver is directly used as a '1 PPS filtering' signal to be output;

(b2) parameter initialization

The time service filter algorithm model is as follows:

x(n|n)＝x(n|n-1)+α[y(n)-x(n|n-1)]

wherein x (n | n) is the filtered value of the nth beat, x (n | n-1) is the predicted value of the nth beat relative to the nth beat,

the filtered value of the nth beat frequency is obtained, y (n) is a measured value of the nth beat, and x (n +1| n) is a predicted value of the nth period to the n +1 th period; alpha and beta are respectively the weight coefficients of the time difference and the frequency difference of the signal of the 1PPS filtering;

T ₀ the time filtering value is as follows: x (0|0) ═ y (0);

T ₀ time frequency filtering value:

T ₁ moment predicted value:

(c)T ₀ validation of the measured value y (0) at a time

The method specifically comprises the following steps:

1) the satellite receiver sends a 'normal time service' message;

2) after the constant temperature crystal oscillator is stabilized, continuously taking 3 beats of measured values, and defining as: y (-3), y (-2), y (-1);

3) the filtering algorithm is started in advance, and values of x (-2| -3), x (-1| -2) and x (0| -1) are extrapolated according to the algorithm model in the step (b) in the step 2;

4) y (0) valid decision, e.g. T _-3 、T _-2 、T _-1 Taking the 4 th observation value as y (0) when no non-point appears at any moment;

note that: firstly, a non-point judgment method is shown in step 4; during the 'y (0) effective judgment', the 1PPS pulse keeps direct output until the filtering is normally started, and the 1PPS predicts the value according to the filtering algorithm; considering that the crystal oscillator is not sufficiently stable during initialization, the frequency filtering value of the crystal oscillator may not be accurate enough, so that the non-point judgment condition is determined in y (0): "| y (n) | ≧ 3 σ/τ" is relaxed to "y (n) | (n) -x (n | n-1) | ≧ 5 σ/τ"), σ is the random error of the satellite receiver 1PPS pulse, τ is the quantization precision of y (n);

and step 3: reading in measured values y (n)

Reading a local time scale corresponding to the rising edge of a 1PPS pulse signal output by a satellite receiver; this step is subdivided into 2 steps:

(a3) y (n) calculating weights

The system measures 1PPS pulse by two clocks CLK with 180 DEG phase difference _HR Counting, in order to prevent the continuous time counter from overflowing, the word width of the measured value information is more than or equal to 48 bits, and two groups of measured values are respectively output by two time scales: y is ₁ (n)，y ₂ (n)；

When y ₁ (n)-y ₂ (n) when ≦ 1, y (n)＝[y ₁ (n)+y ₂ (n)]/2；

When y ₁ (n)-y ₂ When the (n) | > 1 is used,

(b3) non-point culling

According to the pulse precision output by the PPS of the north fighter 1, when | y (n) -x (n | n-1) | ≧ 3 σ/tau, y (n) is abandoned, x (n | n) | x (n | n-1),

and 4, step 4: updating the n-th beat filter value x (n | n)

Obtaining a local time scale x (n | n) of the smoothed second pulse signal of the system 1PPS through a measurement value y (n) of the satellite 1PPS pulse signal according to the algorithm model in the step (b) in the step 2;

and 5: updating the predicted value x (n +1| n) of the nth period to the (n +1) th period

According to the algorithm model in step2 (b), the n-th beat filter value x (n | n) and the n-th beat filter value of the crystal oscillator frequency are obtained

Calculating a predicted value x (n +1| n) of the nth period to the (n +1) th period;

step 6: updating estimated value f (n) of crystal oscillator frequency

According to the measured value y (n) of the satellite 1PPS pulse signal, a unary linear regression model is adopted to solve f (n):

wherein:

f (n) is the estimated value of the nth beat crystal oscillator frequency;

x _i for the ith valid 1PPS pulse number, i.e. x _i I, i-1, 2, 3, 4, … …, n, if presentWhen the beat is not dot, the serial number x _k Discard, i.e. x _i The sequence becomes: … … x _k-2 ，x _k-1 ，x _k+1 ，x _k+2 ……；

size (x) is the total number of indices, and if there are j non-dots, the total number of indices size (x) is n-j,

the serial numbers corresponding to j non-points are also removed;

y _i an observed value y (i) of 1PPS pulse of the ith beat; if the k-th beat is not dot, y _k Abandoning and simultaneously corresponding serial number x _k Are also discarded, i.e. if y ₅ When it is not a point, x _i The sequence becomes: 1, 2, 3, 4, 6, 7 … …; size (x) n-1;

size (y) is the total number of observed points, size (y) size (x);

note that:

y (n) participating in calculation is effective y (n) value;

secondly, in the timekeeping state f (n), the calculation is suspended, and after the timekeeping is quitted, the historical point can still be normally used and the calculation is continued;

if limited by the calculation scale, when y (n) needs to be sampled, the method can be carried out according to the following steps:

when the calculation scale is k points and the diffusion times are q times, the diffusion completion duration r is:

the method comprises the following specific steps:

step 1: defining the first valid observation as y ₁ Y (1) with the number x ₁ ＝x(1)＝1，i＝1；

Step 2: when i is less than or equal to k, calculating y (n) in real time according to the actual points;

step 3: when k is more than i and less than or equal to 2k, the sliding window calculates y (n);

step 4: when i is more than 2k and less than or equal to 3k, starting the 1 st diffusion, wherein the diffusion rule is as follows:

when i is 2k +1, taking k +1, k +3, k +4, k +5, k +6, ·, 2k and 2k + 1;

when i is 2k +2, k +1, k +3, k +5, k +6, · 2k, 2k +1, 2k +2 are taken;

thirdly, the rest is done by analogy;

when i is 3k-1 and i is 3k, finishing the 1 st diffusion of k +1, k +3, 3k-3 and 3 k-1;

step 5: i is more than 3k and less than or equal to 4k, starting the 2 nd diffusion, wherein the diffusion rule is as follows:

when i is 3k +1 and i is 3k +2, k +1, k +4, k +7, k +9, k +11, k +13, ·, 3k-1 and 3k +1 are taken;

when i is 3k +3 and i is 3k +4, k +1, k +4, k +7, k +10, k +13,. loganj, 3k-1, 3k +1 and 3k +3 are taken;

thirdly, the rest is done by analogy;

when i is 4K-1 and i is 4K, finishing the 2 nd diffusion of K +1, K +4, K +7, 4K-5 and 4K-2;

step 6: according to the method described in the steps 4-5, when i is more than 4k and less than or equal to 5k, beginning to perform diffusion for the 3 rd time;

step 7: according to the method described in the steps 4-5, when i is more than 5k and less than or equal to 6k, starting to perform the 4 th diffusion;

step 8: by analogy, when (q +1) k is less than i and less than or equal to (q +2) k, performing the q-th diffusion;

step 9: when i is equal to (q +2) k, the q diffusion is completed, equal-interval sampling with the dot pitch of q +1 is formed, and when i is equal to or more than (q +2) k, according to the sampling state after the q diffusion, namely using the fixed dot pitch of q +1, y (n) is calculated by a sliding window;

the following points are explained for the above sampling algorithm:

firstly, if the algorithm enters the time-keeping state in the process, the algorithm is suspended, and after the algorithm exits the time-keeping state, the effective y can be continuously used _n Participate in y (n) calculation, note: number x _i Always corresponding to absolute time;

secondly, the 1 st diffusion is carried out after 1h, and the stability of the crystal oscillator reaches a nominal value at the moment, so that the unstable period of the crystal oscillator is passed;

after the 1 st diffusion starts, all data participate in the calculation of y (n) until the 14 th diffusion is completed, so that the sampling logic of y (n) is strictly checked, and the influence of non-points and the like on the calculation precision of y (n) is avoided;

and 7: time keeping state

When the 1PPS pulse signal output by the satellite receiver is unstable or interrupted, the system enters a timekeeping state;

the timekeeping state can be subdivided into 3 steps as follows:

(a7) time-keeping state parameter assignment

x(n|n)＝x(n|n-1)，

(b7) Entry condition of time keeping state

The state of keeping watch is entered when one of the following conditions is satisfied:

firstly, the satellite receiver sends out a message of abnormal time service;

secondly, when non-points appear in 3 continuous beats, entering a time keeping mode;

(c7) principle of use of f (n)

1) F (n) updating values measured by using factory binding frequency or continuous 24h before the temperature of the constant-temperature crystal oscillator is stable after the system is started and meets the time keeping condition;

2) after the system is started, if the valid y (n) is greater than 3600 points and the time keeping condition is met, the current frequency estimation value f (n) is used for time keeping, otherwise, the value of Y (n) of factory binding or the value of Y (n) measured for 24h continuously is still used for updating.

4. The method as claimed in claim 3, wherein α and β are 0.4 and 0.1, respectively.

5. The method as claimed in claim 3, wherein k is 1800, q is 14, and r is 8 h.

Images