CN111565084A - A satellite timing 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 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 invention has the advantages of strong universality, simple circuit structure, high locking speed of 1PPS signals and the like.

Description

A satellite timing system and method based on frequency estimation

technical field

The present invention relates to the technology in the field of communication technologies, in particular to a satellite timing system and method based on frequency estimation.

Background technique

The use of satellites for timing and punctuality is essentially the process of taming the clock. By taming the clock, the timing system can follow the satellite time and output UTC time (Coordinated Universal Time). The commonly used method is: the crystal oscillator frequency calibration method, that is, the fixed frequency signal provided by the satellite timing receiver is used to compare with the oscillating signal generated by the local crystal oscillator to obtain the frequency difference; and then through the adjustment of the local crystal oscillator, Make the oscillation frequency basically the same as the oscillation frequency of the satellite (the first edition of "The Principle and Application of Satellite Timing" by Yang Jun and Shan Qingxiao, National Defense Industry Press, August 2013, p. 129), and then generate a second pulse signal (1PPS signal) , to achieve time synchronization.

The satellite timing system using the above method is usually composed of satellite receiver, main control chip (usually FPGA), crystal oscillator, D/A converter, clock interface circuit, etc. Jun, Shan Qingxiao, National Defense Industry Press, August 2013, p. 143).

SUMMARY OF THE INVENTION

In view of the problems existing in the prior art, the present invention provides a satellite timing system based on frequency estimation, including a main control computer, a satellite receiver, a constant temperature crystal oscillator and an FPGA module; wherein

The satellite receiver, which processes and transforms the satellite signal received by the antenna; realizes the capture, tracking, carrier recovery and demodulation of the satellite signal; outputs "1PPS satellite signal" and timing message information to the FPGA module;

An oven-controlled crystal oscillator, which provides a clock reference to the FPGA module;

The FPGA module is used to realize the integrated frequency multiplier and continuous time counter to form a local continuous time standard; to measure the time interval of the 1PPS signal sent by the satellite receiver; to output the filtered local 1PPS signal, that is, "1PPS filter" ";

The main control computer is used to provide calculation and display control services for the whole system, and realize the filtering algorithm of the timing system; the main control computer obtains the "1PPS time scale" and "1PPS time scale" output by the 1PPS signal time interval measurer in the FPGA module. Corresponds to the local time scale count when the rising edge of the "1PPS satellite signal" arrives; according to the characteristic that the mean error of the "1PPS" pulse is zero in the long steady state, the main control computer performs statistical filtering on the "time scale" information of each second. , estimate the error between the center frequency of the constant temperature crystal oscillator and the nominal frequency, and predict the output position of the next "1PPS" second pulse to form a new local 1PPS predicted second data "(n+1) second data", output To the 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 the punctual state, and maintains the stable output of the "1PPS filter" signal.

In one embodiment of the present invention, the FPGA module is composed of a frequency multiplier, a continuous time counter, a 1PPS signal time interval measurer, and a local 1PPS signal generator, and these modules are all implemented by programmable logic units 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 it;

Continuous time counter, which performs continuous time calculation on the CLK _HR frequency signal from the frequency multiplier to form a local continuous time reference; in order to reduce the time scale step, the frequency multiplier outputs two clock signals CLK _HR _0 and CLK _HR _pi, which are The two signals differ by 180° in phase, respectively control two independent continuous time counters to count, forming local time stamp 1 and local time stamp 2, the minimum step size of each time stamp is 1/CLK _HR , and the joint step size is 1/2CLK _HR ; Time stamp 1 is divided into two channels, one is output to the 1PPS signal time interval measurer, and the other is output to the local 1PPS signal generator; Time stamp 2 is output to the 1PPS signal time interval measurer;

1PPS signal time interval measurer, which measures the leading edge of the "1PPS satellite" signal from the satellite receiver, forms "1PPS time scale" information and sends it to the main control computer;

The local 1PPS signal generator, which receives the "time stamp 1" information from the continuous time counter and the "(n+1) second data" information from the main control computer, when the "time stamp 1" and "(n+1) second data" information When equal, a "1PPS filtered" signal is formed.

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

A working method of a satellite timing system based on frequency estimation is also provided, and the steps are as follows:

Step 1: Connect the System

connecting parts of a system as claimed in any one of claims 1 to 3;

Step 2: Initialize

After the system is powered on, perform system initialization, which is subdivided into 5 steps:

(a2) Thru output

After power-on and before the convergence of the filtering algorithm, in order to quickly give the UTC time, the 1PPS second pulse signal sent by the satellite receiver is directly output as a "1PPS filtering" signal by using a pass-through method;

(b2) Parameter initialization

The timing filter algorithm model is as follows:

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

为第n拍频率滤波值，y(n)为第n拍测量值，x(n+1|n)为第n周期对第n+1周期预测值；α、β分别是“1PPS滤波”信号时差和频差的权系数；Among them, x(n|n) is the filter value of the nth beat, x(n|n-1) is the prediction value of the n-1th beat to the nth beat,

is the frequency filter value of the nth beat, y(n) is the measured value of the nth beat, x(n+1|n) is the predicted value of the nth cycle for the n+1th cycle; α and β are the “1PPS filter” signals respectively Weights of time difference and frequency difference;

Filter value at time T ₀ : x(0|0)=y(0);

Frequency filter value at time T ₀ :

Predicted value at time T1 _:

(c) Effective decision of the measured value y(0) at time T ₀

Specifically include the following steps:

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

2) After the constant temperature crystal oscillator is stable, take 3 consecutive beats of the measured values, which are defined as: y(-3), y(-2), y(-1);

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

4) y(0) is a valid judgment, if no non-point occurs at T _-3 , T _-2 , and T _-1 , the observation value of the fourth beat is taken as y(0);

Note: ① For the non-point judgment method, see step 4; ② During the period of "y(0) valid judgment", the 1PPS pulse is always output straight through, and the 1PPS does not predict the value according to the filtering algorithm until the filtering starts normally; ③ Considering the initialization, The crystal oscillator is not sufficiently stable, and the crystal oscillator frequency filter value may not be accurate enough. Therefore, in the judgment of y(0), the non-point judgment condition: "|y(n)-x(n|n-1)|≥3σ/τ" is relaxed is "y(n)-x(n|n-1)|≥5σ/τ"), σ is the random error of the 1PPS pulse of the satellite receiver, and τ is the quantization accuracy of y(n);

Step 3: Read in the measured value y(n)

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

(a3)y(n) calculation weight

The measurement of the 1PPS pulse by the system is obtained by counting two clocks CLK _HR with a phase difference of 180°. In order to prevent the continuous time counter from overflowing, the width of the measured value information word is ≥48 bits, and the two time stamps output two sets of measured values: y ₁ (n), y ₂ (n);

When |y ₁ (n)-y ₂ (n)|≤1, y(n)=[y ₁ (n)+y ₂ (n)]/2;

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

(b3) Non-point culling

According to the pulse precision of the Beidou board 1PPS output, when |y(n)-x(n|n-1)|≥3σ/τ, discard y(n), let x(n|n)=x(n|n -1),

Step 4: Update the filter value x(n|n) for the nth beat

According to the algorithm model in step (b) in step 2, through the measured value y(n) of the satellite 1PPS pulse signal, obtain the smoothed local time scale x(n|n) of the system 1PPS second pulse signal;

Step 5: Update the predicted value x(n+1|n) of the nth cycle for the n+1th cycle

计算第n周期对第n+1周期预测值x(n+1|n)；According to the algorithm model in step (b) in step 2, the filter value x(n|n) of the nth beat and the filter value of the nth beat of the crystal frequency are

Calculate the predicted value x(n+1|n) of the nth cycle to the n+1th cycle;

Step 6: Update the crystal frequency estimate f(n)

According to the measured value y(n) of the satellite 1PPS pulse signal, f(n) is calculated by using a linear regression model:

in:

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

x _i is the i-th effective 1PPS pulse sequence number, that is, x _i =i, i=1, 2, 3, 4,..., n, if there is a non-point in the k-th beat, the sequence number x _k is discarded, that is, x _{The i} sequence becomes: ... x _k-2 , x _k-1 , x _k+1 , x _k+2 ......;

size(x)为序号总数，如果存在j个非点，序号总数size(x)＝n-j，

也应对应剔除j个非点所对应的序号；

size(x) is the total number of serial numbers. If there are j non-points, the total number of serial numbers size(x)=nj,

The serial numbers corresponding to j non-points should also be eliminated correspondingly;

y _i is the observed value y(i) of the 1PPS pulse in the i-th beat; if the k-th beat is a non-point, y _k is discarded, and at the same time, the corresponding serial number x _k is also discarded, that is, if y ₅ is a non-point, the x _i sequence Becomes: 1, 2, 3, 4, 6, 7...; size(x)=n-1;

size(y)为观测值总点数，size(y)＝size(x)；

size(y) is the total number of observations, size(y)=size(x);

Notice:

均为有效y(n)值；①Participate in the calculation

All are valid y(n) values;

② In the punctuality state f(n), the calculation is suspended, and after exiting the punctuality, the historical point can still be used normally, and the calculation continues;

③ If it is limited by the calculation scale, when y(n) needs to be sampled, the following steps can be performed:

When the calculation scale is k points and the number of diffusion is q times, the diffusion completion duration r is:

Specific steps are as follows:

Step1: Define the first valid observation as y ₁ =y(1), the serial number is x ₁ =x(1)=1, i=1;

Step2: When i≤k, calculate y(n) in real time according to the actual number of points;

Step3: When k＜i≤2k, the sliding window calculates y(n);

Step4: When 2k<i≤3k, start the first diffusion, and the diffusion rule is:

①When i=2k+1, take k+1, k+3, k+4, k+5, k+6,..., 2k, 2k+1;

②When i=2k+2, take k+1, k+3, k+5, k+6,..., 2k, 2k+1, 2k+2;

③ and so on;

④ When i=3k-1, i=3k, take k+1, k+3, ..., 3k-3, 3k-1 to complete the first diffusion;

Step5: 3k＜i≤4k, start the second diffusion, the diffusion rule is:

①When i=3k+1, i=3k+2, take k+1, k+4, k+7, k+9, k+11, k+13,..., 3k-1, 3k+1;

②When i=3k+3, i=3k+4, take k+1, k+4, k+7, k+10, k+13,..., 3k-1, 3k+1, 3k+3;

③ and so on;

④ When i=4k-1, i=4k, take k+1, K+4, k+7..., 4k-5, 4k-2 and the second diffusion is completed;

Step6: According to the method described in steps 4 to 5, when 4k<i≤5k, start the third diffusion;

Step7: According to the method described in steps 4 to 5, when 5k<i≤6k, start the fourth diffusion;

Step8: By analogy, when (q+1)k<i≤(q+2)k, start the qth diffusion;

Step9: When i=(q+2)k, the qth diffusion is completed, and the sampling at equal intervals with a point distance of q+1 has been formed. When i≥(q+2)k, according to the qth diffusion Sampling state, that is, using a fixed point distance of q+1, the sliding window calculates y(n);

A few notes on the sampling algorithm above:

①If you enter the punctual state during the above process, the algorithm suspends operation. After exiting punctuality, you can continue to use the valid y _n to participate in the y(n) calculation. Note: the serial number x _i always corresponds to the absolute time;

②The first diffusion is carried out after 1h, at this time, the stability of the crystal oscillator should reach the nominal value, and the unstable period of the crystal oscillator will pass;

③ After the first diffusion starts, until the 14th diffusion is completed, all data participate in the calculation of y(n), so the sampling logic of y(n) should be strictly checked to avoid non-points and other influences on the calculation of y(n) precision;

Step 7: Punctuality Status

When the 1PPS pulse signal output by the satellite receiver is unstable or interrupted, the system enters the punctual state; the punctual state can be subdivided into the following three steps:

(a7) Parameter assignment in punctual state

x(n|n)=x(n|n-1),

(b7) Entry conditions for the punctual state

One of the following conditions is met to enter the punctual state:

①The satellite receiver sends out the message "Time service is abnormal";

②When there is a non-point in 3 consecutive beats, enter the time-keeping mode;

(c7) Principles of use of f(n)

1) After the system is started and before the temperature of the constant temperature crystal oscillator is stable, when the punctuality conditions are met, f(n) is updated using the factory binding frequency or f(n) measured continuously for 24 hours;

2) After the system is started, when the effective y(n) is greater than 3600 points and the punctuality condition is met, use the current frequency estimate f(n) for punctuality, otherwise, still use the factory-bound Y(n) value or continue for 24h Measured Y(n) update value.

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

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

The system and method of the present invention estimate the frequency of the constant temperature crystal oscillator 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 the tracking and synchronization of the UCT time, and the timing and synchronization of the satellite time. Punctuality filtering. The system and method do not need to use a D/A converter to tame the constant temperature crystal oscillator, the crystal oscillator can quickly enter a stable period, and the 1PPS signal can be locked in about 30s, and the locking speed is more than an order of magnitude better than the crystal oscillator frequency calibration method.

After receiving the 1PPS second pulse signal from the satellite receiver, the system and method of the invention outputs a stable 1PPS second pulse signal after filtering and eliminating the non-points, filling the missing points and reducing random fluctuations; The frequency of the output clock signal is measured and estimated. When the 1PPS signal output by the satellite receiver has abnormal states such as interruption, translation, and out-of-tolerance, the system uses the estimated frequency of the clock signal to keep time to ensure the stability of the local 1PPS signal. The invention has the advantages of strong versatility, simple circuit structure, fast locking speed of 1PPS signal and the like.

Description of drawings

Fig. 1 is the working principle diagram of realizing the satellite time service punctuality system based on frequency estimation;

FIG. 2 shows the workflow of the satellite timing system method based on frequency estimation.

Specific implementation method

The satellite timing system based on frequency estimation of the present invention will be described in detail below with reference to the accompanying drawings.

As shown in FIG. 1 , the present invention provides a satellite timing system based on frequency estimation, including a main control computer, a satellite receiver, a constant temperature crystal oscillator and an FPGA module, as described below.

The satellite receiver processes and transforms the satellite signal received by the antenna; realizes the capture, tracking, carrier recovery and demodulation of the satellite signal; outputs the 1PPS second pulse signal (hereinafter referred to as "1PPS satellite signal") and the timing message to the FPGA module and other information.

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

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

In the system of the present invention, the FPGA module is composed of a frequency multiplier, a continuous time counter, a 1PPS signal time interval measurer and a local 1PPS signal generator, and these modules are all implemented by programmable logic units 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 it. In this example, the oscillation frequency of CLK _HR is 200MHz. Others can also be set according to the system design requirements. appropriate clock frequency.

A continuous time counter that performs continuous time calculations on the CLK _HR frequency signal from the frequency multiplier to form a local continuous time reference. In order to reduce the time scale step size, the frequency multiplier outputs two clock signals CLK _HR _0 and CLK _HR _pi, which are 180° out of phase in phase, respectively control two independent continuous time counters to count, forming a local time scale 1 and local time stamp 2, the minimum step size of each time stamp is 1/CLK _HR (for example, 5ns), and the combined step size is 1/2CLK _HR (for example, 2.5ns). The local time scale is generated by a continuous time counter (its counting clock is a constant temperature crystal oscillator), and the stability is very high. Time stamp 1 is divided into two channels, one is output to the 1PPS signal time interval measurer, and the other is output to the local 1PPS signal generator; the time stamp 2 is output to the 1PPS signal time interval measurer.

The 1PPS signal time interval measurer measures the leading edge of the "1PPS satellite" signal from the satellite receiver to form the "1PPS time scale" information and send it to the main control computer.

The local 1PPS signal generator receives the "time stamp 1" information from the continuous time counter and the "(n+1) second data" information from the main control computer, when the "time stamp 1" and "(n+1) second data" information When equal, a "1PPS filtered" signal is formed.

The main control computer is used to provide calculation and display control services for the whole system, and realize the filtering algorithm of the timing system. Considering the system's real-time performance, robustness and program writing achievability requirements, the timing filter adopts α-β filtering algorithm ("An Improved α-β Filtering Algorithm" Holliman Modern Electronic Technology, 2012, No. 21), The punctual filter adopts the univariate linear regression algorithm (the fifth edition of "Applied Regression Analysis" edited by He Xiaoqun and Liu Wenqing, Renmin University of China Press, July 2019, pp. 15-24). The main control computer obtains the 1PPS time stamp output by the 1PPS signal time interval measurer in the FPGA module, and "1PPS time stamp" corresponds to the local time stamp count when the rising edge of the 1PPS satellite signal arrives. According to the characteristic that the mean error of "1PPS" pulse is zero in the long steady state, the main control computer performs statistical filtering processing on the "time scale" information of each second, and estimates the error between the center frequency of the constant temperature crystal oscillator and the nominal frequency, and Predict the output position of the next "1PPS" second pulse, form new local 1PPS predicted second data "(n+1) second data", and output it to the 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 the punctual state, and maintains the stable output of the "1PPS filter" signal. Due to the high stability of the frequency source, the punctuality accuracy of the local time can be effectively guaranteed.

The satellite timing system of the present invention filters and optimizes the second pulse signal "1PPS satellite" output by the satellite receiver through a high-stable frequency source (constant temperature crystal oscillator). The main control computer and the FPGA module output the local pulse-per-second signal "1PPS filter" whose stability meets the requirements of the system index. At the same time, through the long-term statistics of the "1PPS satellite" signal, the center frequency of the constant temperature crystal oscillator is measured and estimated to ensure that when the system satellite signal is abnormal, the timing module continues to output the second pulse signal normally, and maintains the output of the "1PPS filtered" signal. Stable and meet the requirements of punctual accuracy.

As shown in Figure 2, the working method steps of the satellite timing system based on frequency estimation of the present invention are as follows:

Step 1: Connect the System

Each part of the system of the present invention is connected according to the relationship shown in Figure 1;

Step 2: Initialize

After the system is powered on, perform system initialization, which can be subdivided into 5 steps:

(a) Thru output

After power-on and before the convergence of the filtering algorithm, in order to quickly give the UTC time, the 1PPS second pulse signal sent by the satellite receiver is directly output as a "1PPS filtering" signal by using a pass-through method;

(b) parameter initialization

The timing filter algorithm model is as follows:

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

为第n拍频率滤波值，y(n)为第n拍测量值，x(n+1|n)为第n周期对第n+1周期预测值。α、β分别是“1PPS滤波”信号时差和频差的权系数，推荐值为0.4和0.1。Among them, x(n|n) is the filter value of the nth beat, x(n|n-1) is the prediction value of the n-1th beat to the nth beat,

is the frequency filtering value of the nth beat, y(n) is the measured value of the nth beat, and x(n+1|n) is the predicted value of the nth cycle to the n+1th cycle. α and β are the weight coefficients of the time difference and frequency difference of the "1PPS filtering" signal, respectively, and the recommended values are 0.4 and 0.1.

Filter value at time T ₀ : x(0|0)=y(0);

Frequency filter value at time T ₀ :

Predicted value at time T1 _:

(c) Effective decision of the measured value y(0) at time T ₀

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

2) After the constant temperature crystal oscillator is stable (usually 180-600s), take 3 consecutive beats of the measured value, which are defined as: y(-3), y(-2), y(-1);

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

4) y(0) is a valid judgment. For example, there is no non-point at time T _-3 , T _-2 , and T _-1 , and the observation value of the fourth beat is taken as y(0).

(Note: ① For the non-point judgment method, see step 4; ② During the period of "y(0) valid judgment", the 1PPS pulse is always output straight through, and the 1PPS will not predict the value according to the filtering algorithm until the filtering starts normally; ③ Considering the initialization , the crystal oscillator is not sufficiently stable, and the crystal oscillator frequency filter value may not be accurate enough, so in the judgment of y(0), the non-point judgment condition: "|y(n)-x(n|n-1)|≥3σ/τ" It can be relaxed to "y(n)-x(n|n-1)|≥5σ/τ"), σ is the random error of the satellite receiver 1PPS pulse, and τ is the quantization accuracy of y(n). In this example, for example is the joint step size of 2.5ns.

Step 3: Read in the measured value y(n)

Read in the 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) calculation weight

The measurement of the 1PPS pulse by the system is obtained by counting the two clocks CLK _HR with a phase difference of 180°. When the CLK _HR clock frequency is 200MHz, the theoretical measurement accuracy is ±2.5ns. In order to prevent the continuous time counter from overflowing, the measurement value information The word width should be as wide as possible, preferably ≥48bit (if 48bit is selected, that is, a one-way cycle of 390h), the two time scales output two sets of measurement values: y ₁ (n), y ₂ (n);

When |y ₁ (n)-y ₂ (n)|≤1, y(n)=[y ₁ (n)+y ₂ (n)]/2;

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

(b) Non-point culling

According to the pulse precision of the Beidou board 1PPS output, when |y(n)-x(n|n-1)|≥3σ/τ, discard y(n), let x(n|n)=x(n|n -1),

Step 4: Update the filter value x(n|n) for the nth beat

According to the algorithm model in step (b) of step 2, through the measured value y(n) of the satellite 1PPS pulse signal, obtain the smoothed local time scale x(n|n) of the system 1PPS second pulse signal.

Step 5: Update the predicted value x(n+1|n) of the nth cycle for the n+1th cycle

计算第n周期对第n+1周期预测值x(n+1|n)。According to the algorithm model in step (b) in step 2, the filter value x(n|n) of the nth beat and the filter value of the nth beat of the crystal frequency are

Calculate the predicted value x(n+1|n) for the nth cycle versus the n+1th cycle.

Step 6: Update the crystal frequency estimate f(n)

According to the measured value y(n) of the satellite 1PPS pulse signal, f(n) is calculated by using a linear regression model:

in:

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

x _i is the i-th effective 1PPS pulse sequence number, that is, x _i =i, i=1, 2, 3, 4,..., n, if there is a non-point in the k-th beat, the sequence number x _k is discarded, that is, x _{The i} sequence becomes: ... x _k-2 , x _k-1 , x _k+1 , x _k+2 ......;

size(x)为序号总数，如果存在j个非点，序号总数size(x)＝n-j，

也应对应剔除j个非点所对应的序号；

size(x) is the total number of serial numbers. If there are j non-points, the total number of serial numbers size(x)=nj,

The serial numbers corresponding to j non-points should also be eliminated correspondingly;

y _i is the observed value y(i) of the 1PPS pulse in the i-th beat; if the k-th beat is a non-point, y _k is discarded, and at the same time, the corresponding serial number x _k is also discarded, that is, if y ₅ is a non-point, the x _i sequence Becomes: 1, 2, 3, 4, 6, 7...; size(x)=n-1;

size(y)为观测值总点数，size(y)＝size(x)；

size(y) is the total number of observations, size(y)=size(x);

Notice:

①The y(n) involved in the calculation are all valid y(n) values;

② In the punctuality state f(n), the calculation is suspended, and after exiting the punctuality, the historical point can still be used normally, and the calculation continues;

③ If it is limited by the calculation scale, when y(n) needs to be sampled, the following steps can be performed:

The following takes the calculation scale k=1800 points, the diffusion times 14 times as an example, and the diffusion completion duration is 8h.

1. Define the first valid observation as y ₁ =y(1), the serial number is x ₁ =x(1)=1, i=1;

2. When i≤k, calculate y(n) in real time according to the actual number of points;

3. When k<i≤2k, the sliding window calculates y(n);

4. When 2k<i≤3k, the first diffusion starts, and the diffusion rule is:

①When i=2k+1, take k+1, k+3, k+4, k+5, k+6,..., 2k, 2k+1;

②When i=2k+2, take k+1, k+3, k+5, k+6,..., 2k, 2k+1, 2k+2;

③ and so on;

④ When i=3k-1, i=3k, take k+1, k+3, ..., 3k-3, 3k-1 to complete the first diffusion;

5. 3k<i≤4k, start the second diffusion, the diffusion rule is:

①When i=3k+1, i=3k+2, take k+1, k+4, k+7, k+9, k+11, k+13,..., 3k-1, 3k+1;

②When i=3k+3, i=3k+4, take k+1, k+4, k+7, k+10, k+13,..., 3k-1, 3k+1, 3k+3;

③ and so on;

④ When i=4k-1, i=4k, take k+1, K+4, k+7..., 4k-5, 4k-2 and the second diffusion is completed;

6. According to the method described in steps 4 to 5, when 4k<i≤5k, start the third diffusion;

7. According to the method described in steps 4 to 5, when 5k<i≤6k, start the fourth diffusion;

8. By analogy, when 15k<i≤16k, the 14th diffusion begins;

9. When i=16k, the 14th diffusion is completed, and the sampling at equal intervals with a point distance of 15 has been formed. When i≥16k, the sliding window is performed according to the sampling state after the 14th diffusion (that is, a fixed point distance of 15 is used). compute y(n);

A few notes on the sampling algorithm above:

①If you enter the punctual state during the above process, the algorithm suspends operation. After exiting punctuality, you can continue to use the valid y _n to participate in the y(n) calculation. Note: the serial number x _i always corresponds to the absolute time;

②The first diffusion is carried out after 1h. At this time, the stability of the crystal oscillator should reach the nominal value, and the unstable period of the crystal oscillator will pass;

③ After the first diffusion starts, until the 14th diffusion is completed, all data participate in the calculation of y(n), so the sampling logic of y(n) should be strictly checked to avoid non-points and other influences on the calculation of y(n) precision.

Step 7: Punctuality Status

When the 1PPS pulse signal output by the satellite receiver is unstable or interrupted, the system enters the punctual state. The punctuality status can be subdivided into the following 3 steps:

(a) Parameter assignment in punctual state

x(n|n)=x(n|n-1),

(b) Entry conditions for the punctual state

One of the following conditions is met to enter the punctual state:

①The satellite receiver sends out the message "Time service is abnormal";

②When there is a non-point in 3 consecutive beats, enter the time-keeping mode.

(c) Principles for the use of f(n)

1) After the system is started and before the temperature of the constant temperature crystal oscillator is stable, when the punctuality condition is met, f(n) is updated using the factory binding frequency or f(n) measured continuously for 24 hours.

2) After the system is started, when the effective y(n) is greater than 3600 points and the punctuality condition is met, use the current frequency estimate f(n) for punctuality, otherwise, still use the factory-bound Y(n) value or continue for 24h Measured Y(n) update value.

The method of the present invention has the following advantages:

The system and method of the present invention receive the 1PPS second pulse signal from the satellite receiver, filter and eliminate the non-points, fill in the missing points and reduce random fluctuations, and then output a stable 1PPS second pulse signal; The frequency is measured and estimated. When the 1PPS signal output by the satellite receiver has abnormal states such as interruption, translation, and out of tolerance, the system uses the estimated frequency of the clock signal to be punctual to ensure the stability of the local 1PPS signal, thereby realizing the present invention. . The method has the advantages of strong versatility, no need to use a D/A converter to tame the constant temperature crystal oscillator, and fast locking speed of the 1PPS signal.

1. Strong versatility

The system can not only be used for time filtering of satellite timing systems, but also for time filtering and smoothing of other timing systems, with strong versatility.

2, the circuit structure is simple

Since the frequency estimation method is used for time filtering and prediction, there is no need to use a D/A converter to tame the frequency of the crystal oscillator, and the circuit structure is simpler than the traditional timing system.

3. 1PPS signal locking speed is fast

When the traditional timing system tames the frequency of the crystal oscillator, it is necessary to continuously adjust the output voltage of the D/A converter, the constant temperature crystal oscillator has a long stabilization time, and the 1PPS signal locking speed is slow. Using the system of the invention, after the satellite receiver 1PPS signal is output, the system algorithm can quickly (within 30s) complete the 1PPS signal locking by utilizing the high short-term stability of the crystal oscillator; using the high long-term stability of the satellite 1PPS signal, the system algorithm performs frequency It can also ensure that the punctuality accuracy meets the requirements.

Claims (6)

1. a satellite timing system based on frequency estimation, is characterized in that, comprises main control computer, satellite receiver, constant temperature crystal oscillator and FPGA module; Wherein The satellite receiver, which processes and transforms the satellite signal received by the antenna; realizes the capture, tracking, carrier recovery and demodulation of the satellite signal; outputs "1PPS satellite signal" and timing message information to the FPGA module; An oven-controlled crystal oscillator, which provides a clock reference to the FPGA module; The FPGA module is used to realize the integrated frequency multiplier and continuous time counter to form a local continuous time standard; measure the time interval of the 1PPS signal sent by the satellite receiver; output the filtered local 1PPS signal, that is, "1PPS filter" "; The main control computer is used to provide calculation and display control services for the whole system, and realize the filtering algorithm of the timing system; the main control computer obtains the "1PPS time scale" and "1PPS time scale" output by the 1PPS signal time interval measurer in the FPGA module. Corresponds to the local time scale count when the rising edge of the "1PPS satellite signal" arrives; according to the characteristic that the mean error of the "1PPS" pulse is zero in the long steady state, the main control computer performs statistical filtering on the "time scale" information of each second. , estimate the error between the center frequency of the constant temperature crystal oscillator and the nominal frequency, and predict the output position of the next "1PPS" second pulse to form a new local 1PPS predicted second data "(n+1) second data", output To the 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 the punctual state, and maintains the stable output of the "1PPS filter" signal. 2. the satellite timing system based on frequency estimation as claimed in claim 1, is characterized in that, FPGA module is made up of frequency multiplier, continuous time counter, 1PPS signal time interval measurer and local 1PPS signal generator, these modules All are implemented by programmable logic units in FPGA, where 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 it; Continuous time counter, which performs continuous time calculation on the CLK _HR frequency signal from the frequency multiplier to form a local continuous time reference; in order to reduce the time scale step, the frequency multiplier outputs two clock signals CLK _HR _0 and CLK _HR _pi, which are The two signals differ by 180° in phase, respectively control two independent continuous time counters to count, forming local time stamp 1 and local time stamp 2, the minimum step size of each time stamp is 1/CLK _HR , and the joint step size is 1/2CLK _HR ; Time stamp 1 is divided into two channels, one is output to the 1PPS signal time interval measurer, and the other is output to the local 1PPS signal generator; Time stamp 2 is output to the 1PPS signal time interval measurer; 1PPS signal time interval measurer, which measures the leading edge of the "1PPS satellite" signal from the satellite receiver, forms "1PPS time scale" information and sends it to the main control computer; The local 1PPS signal generator, which receives the "time stamp 1" information from the continuous time counter and the "(n+1) second data" information from the main control computer, when the "time stamp 1" and "(n+1) second data" information When equal, a "1PPS filtered" signal is formed. 3 . The satellite timing system based on frequency estimation as claimed in claim 2 , wherein the local time scale 1 and the local time scale 2 have a minimum step size of 5ns for each channel, and a combined step size of 2.5ns. 4 . 4. a working method of the satellite timing system based on frequency estimation, is characterized in that, step is as follows: Step 1: Connect the System connecting parts of a system as claimed in any one of claims 1 to 3; Step 2: Initialize After the system is powered on, perform system initialization, which is subdivided into 5 steps: (a2) Thru output After power-on and before the convergence of the filtering algorithm, in order to quickly give the UTC time, the 1PPS second pulse signal sent by the satellite receiver is directly output as a "1PPS filtering" signal by using a pass-through method; (b2) Parameter initialization The timing filter algorithm model is as follows: x(n|n)=x(n|n-1)+α[y(n)-x(n|n-1)]

Among them, x(n|n) is the filter value of the nth beat, x(n|n-1) is the prediction value of the n-1th beat to the nth beat,

is the frequency filter value of the nth beat, y(n) is the measured value of the nth beat, x(n+1|n) is the predicted value of the nth cycle for the n+1th cycle; α and β are the “1PPS filter” signals respectively Weights of time difference and frequency difference; Filter value at time T ₀ : x(0|0)=y(0); Frequency filter value at time T ₀ :

Predicted value at time T1 _:

(c) Effective decision of the measured value y(0) at time T ₀ Specifically include the following steps: 1) The satellite receiver sends a "time service normal" message; 2) After the constant temperature crystal oscillator is stable, take 3 consecutive beats of the measured values, which are defined as: y(-3), y(-2), y(-1); 3) The filtering algorithm is pre-started, and the values of x(-2|-3), x(-1|-2), and x(0|-1) are extrapolated according to the algorithm model in step (b) in step 2; 4) y(0) is a valid judgment, if no non-point occurs at T _-3 , T _-2 , and T _-1 , the observation value of the fourth beat is taken as y(0); Note: ① For the non-point judgment method, see step 4; ② During the period of "y(0) valid judgment", the 1PPS pulse is always output straight through, and the 1PPS will not predict the value according to the filtering algorithm until the filtering starts normally; ③ Considering the initialization, The crystal oscillator is not sufficiently stable, and the crystal oscillator frequency filter value may not be accurate enough. Therefore, in the judgment of y(0), the non-point judgment condition: "|y(n)-x(n|n-1)|≥3σ/τ" is relaxed is "y(n)-x(n|n-1)|≥5σ/τ"), σ is the random error of the 1PPS pulse of the satellite receiver, and τ is the quantization accuracy of y(n); Step 3: Read in the measured value y(n) Read in the local time scale corresponding to the rising edge of the 1PPS pulse signal output by the satellite receiver; this step is subdivided into 2 steps: (a3)y(n) calculation weight The measurement of the 1PPS pulse by the system is obtained by counting two clocks CLK _HR with a phase difference of 180°. In order to prevent the continuous time counter from overflowing, the width of the measured value information word is ≥48 bits, and the two time stamps output two sets of measured values: y ₁ (n), y ₂ (n); When |y ₁ (n)-y ₂ (n)|≤1, y(n)=[y ₁ (n)+y ₂ (n)]/2; When |y ₁ (n)-y ₂ (n)|>1,

(b3) Non-point culling According to the pulse precision of the Beidou board 1PPS output, when |y(n)-x(n|n-1)|≥3σ/τ, discard y(n), let x(n|n)=x(n|n -1),

Step 4: Update the filter value x(n|n) for the nth beat According to the algorithm model in step (b) in step 2, through the measured value y(n) of the satellite 1PPS pulse signal, obtain the smoothed local time scale x(n|n) of the system 1PPS second pulse signal; Step 5: Update the predicted value x(n+1|n) of the nth cycle for the n+1th cycle According to the algorithm model in step (b) in step 2, the filter value x(n|n) of the nth beat and the filter value of the nth beat of the crystal frequency are

Calculate the predicted value x(n+1|n) of the nth cycle to the n+1th cycle; Step 6: Update the crystal frequency estimate f(n) According to the measured value y(n) of the satellite 1PPS pulse signal, f(n) is calculated by using a linear regression model:

in: f(n) is the estimated value of the crystal oscillator frequency of the nth beat; x _i is the i-th effective 1PPS pulse sequence number, that is, x _i =i, i=1, 2, 3, 4,..., n, if there is a non-point in the k-th beat, the sequence number x _k is discarded, that is, x _{The i} sequence becomes: ... x _k-2 , x _k-1 , x _k+1 , x _k+2 ......;

size(x) is the total number of serial numbers. If there are j non-points, the total number of serial numbers size(x)=nj,

The serial numbers corresponding to j non-points should also be eliminated correspondingly; y _i is the observed value y(i) of the 1PPS pulse in the i-th beat; if the k-th beat is a non-point, y _k is discarded, and at the same time, the corresponding serial number x _k is also discarded, that is, if y ₅ is a non-point, the x _i sequence Becomes: 1, 2, 3, 4, 6, 7...; size(x)=n-1;

size(y) is the total number of observations, size(y)=size(x); Notice: ①The y(n) involved in the calculation are all valid y(n) values; ② In the punctuality state f(n), the calculation is suspended, and after exiting the punctuality, the historical point can still be used normally, and the calculation continues; ③ If it is limited by the calculation scale, when y(n) needs to be sampled, the following steps can be performed: When the calculation scale is k points and the number of diffusion is q times, the diffusion completion duration r is:

Specific steps are as follows: Step1: Define the first valid observation as y ₁ =y(1), the serial number is x ₁ =x(1)=1, i=1; Step2: When i≤k, calculate y(n) in real time according to the actual number of points; Step3: When k＜i≤2k, the sliding window calculates y(n); Step4: When 2k<i≤3k, start the first diffusion, and the diffusion rule is: ⑤ When i=2k+1, take k+1, k+3, k+4, k+5, k+6,..., 2k, 2k+1; ⑥ When i=2k+2, take k+1, k+3, k+5, k+6, ..., 2k, 2k+1, 2k+2; ⑦ And so on; ⑧ When i=3k-1, i=3k, take k+1, k+3, ..., 3k-3, 3k-1 to complete the first diffusion; Step5: 3k＜i≤4k, start the second diffusion, the diffusion rule is: ⑤ When i=3k+1, i=3k+2, take k+1, k+4, k+7, k+9, k+11, k+13,..., 3k-1, 3k+1; ⑥ When i=3k+3, i=3k+4, take k+1, k+4, k+7, k+10, k+13,..., 3k-1, 3k+1, 3k+3; ⑦ And so on; ⑧ i=4k-1, when i=4k, take k+1, K+4, k+7..., 4k-5, 4k-2 and the second diffusion is completed; Step6: According to the method described in steps 4 to 5, when 4k<i≤5k, start the third diffusion; Step7: According to the method described in steps 4 to 5, when 5k<i≤6k, start the fourth diffusion; Step8: By analogy, when (q+1)k<i≤(q+2)k, start the qth diffusion; Step9: When i=(q+2)k, the qth diffusion is completed, and the sampling at equal intervals with a point distance of q+1 has been formed. When i≥(q+2)k, the sampling after the qth diffusion is performed. state, that is, using a fixed point distance of q+1, the sliding window calculates y(n); A few notes on the sampling algorithm above: ④If in the above process, enter the punctual state, the algorithm suspends operation, and after exiting punctuality, can continue to use the valid y _n to participate in the y(n) calculation, note: the serial number x _i always corresponds to the absolute time; ⑤ The first diffusion is carried out after 1h, at this time, the stability of the crystal oscillator should reach the nominal value, and the unstable period of the crystal oscillator will pass; ⑥After the first diffusion starts, until the 14th diffusion is completed, all data participate in the calculation of y(n), so the sampling logic of y(n) should be strictly checked to avoid non-points and other influences on the calculation of y(n) precision; Step 7: Punctuality Status When the 1PPS pulse signal output by the satellite receiver is unstable or interrupted, the system enters the punctual state; the punctual state can be subdivided into the following three steps: (a7) Parameter assignment in punctual state x(n|n)=x(n|n-1),

(b7) Entry conditions for the punctual state One of the following conditions is met to enter the punctual state: ①The satellite receiver sends out the message of "not normal timing"; ②When there is a non-point in 3 consecutive beats, enter the time-keeping mode; (c7) Principles of use of f(n) 1) After the system is started and before the temperature of the constant temperature crystal oscillator is stable, when the punctuality condition is met, f(n) uses the factory binding frequency or the updated value of f(n) measured continuously for 24 hours; 2) After the system is started, when the effective y(n) is greater than 3600 points and the punctuality condition is met, use the current frequency estimate f(n) for punctuality, otherwise, still use the factory-bound Y(n) value or continue for 24h Measured Y(n) update value. 5 . The working method of the satellite timing system based on frequency estimation according to claim 4 , wherein α and β are respectively 0.4 and 0.1. 6 . 6 . The working method of the satellite timing system based on frequency estimation according to claim 4 , wherein k=1800, q=14, and r=8h. 7 .

Images