CN111711446B - Method, system and medium for taming crystal oscillator frequency by using GPS signal - Google Patents

Abstract

The invention discloses a method, a system and a medium for taming crystal oscillator frequency by using GPS signals, which comprises the following steps: establishing a frequency difference mathematical model according to error characteristics of the GPS second signal and the crystal oscillator signal; performing conventional extended Kalman filtering processing on the frequency difference signal to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency, and obtaining a filtering model greatly influenced by a jump wild value and system noise; introducing a differential evolution algorithm to dynamically adjust process noise in conventional extended Kalman filtering processing so as to reduce the influence of the process noise on a frequency difference signal filtering result; and integrating the weight of the innovation sequence into the extended Kalman filtering process after the dynamic noise adjustment in the process is completed so as to reduce the influence of the jump wild value on the filtering result. The invention can carry out filtering treatment on the jump wild value and random error in the GPS signal, the accumulated error in the crystal oscillator signal and nonlinear frequency drift, and can improve the final frequency correction precision of the GPS locking secondary frequency standard.

Description

A method, system and medium for taming crystal oscillator frequency by using GPS signal

technical field

The invention relates to a frequency source error filtering technology for GPS taming crystal oscillator signals, in particular to a method, system and medium for taming crystal oscillator frequencies using GPS signals, which are used to realize high-precision taming of crystal oscillator frequencies by GPS signals.

Background technique

With the rapid development of network technology, the system has higher and higher requirements for time reference. Time base in power system is of great significance to accident analysis, improvement of power quality, power grid dispatching operation and time-of-use billing of power energy.

Current time (frequency) references can be divided into primary frequency standards (cesium atomic frequency standards), secondary frequency standards (including rubidium atomic frequency standards and high-stability crystal oscillators) and other frequency standards (high-stability crystal oscillators) according to application fields and performance indicators. other crystal oscillators except crystal oscillators). Although the long-term stability and accuracy of the first-level frequency standard are relatively high, its price is very expensive, and it is difficult to apply it in civil occasions that have strict requirements on cost; the price of the second-level frequency standard is relatively low, but its long-term stability and accuracy are greater. Not as good as the former, so it is difficult to apply to high-precision occasions.

The principle of using GPS to lock the secondary frequency standard (high stability crystal oscillator) in the prior art is shown in Figure 1: First, measure the relative frequency difference between the second pulse signal PPS output by the GPS receiver and the pulse signal output by the crystal oscillator after frequency division, Then filter the frequency difference signal, and finally generate a control signal to continuously adjust the frequency of the crystal oscillator to output a high-accuracy and high-stability frequency output. This technology solves the above-mentioned problems to a certain extent, that is, realizes frequency standards with high long-term stability and average accuracy at a lower cost. However, in practical applications, due to the random error and jump outliers of the 1PPS signal of GPS, and the frequency drift of the crystal oscillator due to aging and temperature characteristics, it has caused certain difficulties in frequency calibration.

In order to solve the above problems, some scholars proposed to use the method of taking the average of multiple samples to process the measured frequency difference data. The calculation is simple but the frequency calibration takes a long time. Some scholars have established a mathematical model of the frequency difference between the crystal oscillator signal and the GPS signal, and used the least square method to calculate and compensate its error. However, this method is more complicated to calculate and does not consider the influence of the outlier value of the GPS signal on the frequency calibration accuracy. Some scholars have proposed to use the Kalman filter algorithm to eliminate the noise in the frequency difference signal, but the existing method of using the Kalman filter algorithm to eliminate the noise in the frequency difference signal does not consider the influence of the crystal oscillator frequency drift, resulting in limited frequency calibration accuracy.

Contents of the invention

The technical problem to be solved by the present invention: Aiming at the above-mentioned problems of the prior art, a method, system and medium for taming the frequency of the crystal oscillator by using the GPS signal are provided. Filtering outliers and random errors in the GPS signal, cumulative errors and nonlinear frequency drift in the crystal oscillator signal can improve the final frequency calibration accuracy of the GPS locked secondary frequency standard.

In order to solve the problems of the technologies described above, the technical solution adopted in the present invention is:

A kind of method utilizing GPS signal to tame crystal oscillator frequency, the step of this method comprises:

1) According to the error characteristics of the GPS second signal and the crystal signal, the mathematical model of the frequency difference between the two is established;

2) Perform conventional extended Kalman filter processing on the frequency difference signal according to the frequency difference mathematical model, thereby eliminating random errors and errors caused by non-linear drift of the crystal oscillator frequency, and obtaining a filtering model that is greatly affected by jump outliers and system noise;

3) Introduce the differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process to reduce its influence on the filtering results of the frequency difference signal;

4) Integrate the weighting of the innovation sequence into the extended Kalman filter processing after the process noise dynamic adjustment is completed, so as to reduce the influence of jump outliers on the filtering results of the frequency difference signal.

Optionally, the detailed steps of step 1) include:

1.1) obtain the GPS second signal sequence Y, the functional expression of the universal time UTC of the GPS second signal at any k moment in the GPS second signal sequence Y as shown in formula (1);

Y _k = k-ε _k (1)

In the above formula, Y _k represents the universal time UTC of the GPS second signal at time k in the GPS second signal sequence Y, ε _k is the error of the GPS second signal at k time, and satisfies ε～N(0,σ ² );

Obtain the crystal oscillator signal sequence Y ', the functional expression of the universal time UTC of the crystal oscillator signal at any k moment in the crystal oscillator signal sequence Y ' as shown in formula (2);

Y′ _k ＝k+a+bk+ck ² (2)

In the above formula, Y' _k represents the universal time UTC of the crystal oscillator signal at time k in the crystal oscillator signal sequence Y', a is the initial phase error, b is the error coefficient considering the frequency deviation, and c is the error coefficient considering the frequency linear drift ;

1.2) According to the GPS second signal sequence Y, the crystal oscillator signal sequence Y', the difference is obtained to obtain the frequency difference sequence, and the functional expression of the frequency difference signal at any k moment in the frequency difference sequence is shown in the following formula;

X _k ＝Y′ _k -Y _k ＝a+bk+ck ² +ε _k (3)

In the above formula, X _k represents the frequency difference signal at time k in the frequency difference sequence.

Optionally, the step of performing conventional extended Kalman filter processing on the frequency difference signal according to the frequency difference mathematical model in step 2) includes:

2.1) The functional expressions of the state equation and the observation equation of the frequency difference are established as shown in formula (4);

In the above formula, x(k+1) represents the frequency difference state vector at time k+1, x(k) represents the frequency difference state vector at time k, f and h represent the state function and measurement function respectively, z(k+ 1) represents the measurement value at time k+1; ω(k) represents the process noise vector at time k, and its covariance matrix is Q; γ(k) represents the measurement noise vector at time k, and its covariance matrix is R ;

2.2) Let the functional expressions of the state transition matrix F(k+1,k) from k time to k+1 time and the observation matrix H(k) at k time be shown in formula (5);

In the above formula, f represents the state function, Represents the estimation of the frequency difference state vector at time k, h represents the measurement function, and x(k) represents the frequency difference state vector at time k;

2.3) The prediction equations for establishing the conventional extended Kalman filter model are shown in equations (6) to (8), the update equations are shown in equations (9) to (10), and the process noise is shown in equation (11). The noise variance is shown in formula (12);

In formula (6), Represents the predicted value of the frequency difference state vector at time k+1 based on the predicted value at time k, f represents the state function, /> Represents the optimal estimated value of the frequency difference state vector at time k;

p(k+1,k)=F(k+1,k)p(k,k)F ^T (k+1,k)+Q(k) (7)

In formula (7), p(k+1,k) represents the predicted value of the covariance vector at time k+1, F(k+1,k) represents the state transition matrix from time k to time k+1, p(k , k) represents the predicted value of the covariance vector at time k, F ^T (k+1,k) represents the transposition of the state transition matrix from time k to time k+1, Q(k) represents the process noise variance at time k;

K(k+1)=p(k+1,k)H ^T (k+1)[H(k+1)p(k+1,k)H ^T (k+1)+R(k+1 )] ^-1 (8)

In formula (8), K(k+1) represents the Kalman filter gain matrix at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, H ^T (k+1) Represents the transposition of the observation matrix H(k+1) at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, R(k+1) represents the predicted value of the covariance vector at time k+1 Observation noise variance vector;

In formula (9), Indicates the optimal estimated value of the frequency difference state vector at time k+1, /> Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, z(k+1) represents the measured value at time k+1, /> Indicates the predicted value of the frequency difference state vector at time k+1 /> The obtained k+1 time measurement vector prediction value;

p(k+1)=[I-K(k+1)H(k+1)]p(k+1,k) (10)

In formula (10), p(k+1) represents the covariance vector estimated value updated at time k+1, I represents the identity matrix with the same dimension as K(k+1)H(k+1), and K(k +1) represents the Kalman filter gain matrix at time k+1, H(k+1) represents the observation matrix at time k+1, and p(k+1,k) represents the predicted value of the covariance vector at time k+1;

In formula (11), ω(k+1) represents the process noise at time k+1, Represents the optimal estimated value of the frequency difference state vector at time k+1, and z(k+1) represents the measured value at time k+1;

Q _k+1 = var(ω(k+1))(12)

In formula (12), Q _k+1 represents the variance of process noise at time k+1, ω(k+1) represents the process noise at time k+1, and var represents the operation of calculating the variance of process noise.

Optionally, when the differential evolution algorithm is introduced in step 3) to dynamically adjust the process noise in the conventional extended Kalman filtering process, the step of generating the variance of the process noise at any time k by the differential evolution algorithm includes:

3.1) Initialize the population, and the function expression is shown in the following formula (13):

x _i,k (0)＝l _k +rand(0,1)*(u _k -l _k ) (13)

In formula (13), x _i,k (0) represents the kth gene of the i-th "chromosome" of the 0th generation in the population; rand(0,1) is a random number uniformly distributed between 0 and 1; u _k and l _k are the upper bound and lower bound of the search, and u _k and l _k respectively take the maximum value Q _max and the minimum value Q _min of the process noise variance;

3.2) The mutation operation is performed on the population, and the function expression of the mutation operation is shown in the following formula (14):

X _i (g+1)=x _r1 (g)+F*[x _r2 (g)-x _r3 (g)] (14)

In formula (14), _Xi (g+1) is the obtained g+1 generation mutant individual; F is the compression scale factor, the value range is 0~1; x _r1 (g), x _r2 (g) and x _r3 (g) being 3 parents;

3.3) Perform a crossover operation on the population to retain better variables, and the crossover operation adopts the binomial crossover method, and its function expression is shown in the following formula (15):

In formula (15), y _i,j (g+1) represents the crossover operation result corresponding to the current individual, Xi _,j (g+1) represents the component corresponding to the mutant individual, and xi _,j (g) represents the current individual The corresponding component, r is a random number uniformly distributed between 0 and 1 generated by each variable; cr is the crossover probability of the variable; rnd is an integer uniformly distributed between 1 and d, if r<cr, the mutant individual is accepted The corresponding component X _i,j (g+1), otherwise keep the component x _i,j (g) corresponding to the current individual;

3.4) For the results obtained by the crossover operation, the optimal solution is selected as the process noise variance at time k by using the greedy selection method shown in the following formula (16);

In formula (16), x _i (g+1) is the selection result, y _i (g+1) represents the individual corresponding to the crossover operation result, _xi (g) represents the current individual, f[y _i (g+1) ] represents the population fitness value corresponding to individual y _i (g+1), and f[ _xi (g)] represents the population fitness value corresponding to individual _xi (g).

Optionally, after step 4) integrating the weighting of the innovation sequence into the extended Kalman filter processing after the process noise is dynamically adjusted, the function expression of the update equation of the extended Kalman filter processing obtained is as follows: (17)～( 21) as shown;

D _k+1 ＝P(k+1,k)+R _k+1 (21)

In formula (17)～(21), Indicates the optimal estimated value of the frequency difference state vector at time k+1 after adjustment, /> Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, φ _k+1 represents the correction function, and τ _k+1 represents the measurement error at time k+1 The variable of size, e _k+1 represents the measurement error at k+1 time, z(k+1) represents the measured value of the frequency difference vector at k+1 time, /> Indicates the predicted value of the measurement vector at time k+1, τ _k+1 obeys the χ ² distribution, when an outlier occurs, τ _k+1 increases and φ _k+1 ( ) decreases; C _k+1 is the selected Threshold value or constant sequence, λ _k+1 is The largest eigenvalue of the matrix, K _k+1 represents the Kalman filter gain at time k+1, D _k+1 represents the weight matrix that measures the measurement error at time k+1, and P(k+1,k) represents k+ The covariance prediction value at time 1, R(k+1) represents the measurement noise variance at time k+1.

In addition, the present invention also provides a system for using GPS signals to tame the frequency of the crystal oscillator, including:

A frequency difference model generation program unit is used to set up both frequency difference mathematical models according to the error characteristics of the GPS second signal and the crystal oscillator signal;

The conventional filter program unit is used to perform conventional extended Kalman filter processing on the frequency difference signal according to the frequency difference mathematical model, thereby eliminating random errors and errors caused by non-linear drift of the crystal oscillator frequency, and obtaining a relatively low frequency value affected by jump outliers and system noise. Large filtering model;

The noise adjustment program unit is used to introduce a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process to reduce its influence on the frequency difference signal filtering result;

The filter integration program unit is used to integrate the weight of the innovation sequence into the extended Kalman filter processing after the process noise is dynamically adjusted.

In addition, the present invention also provides a system for taming the frequency of the crystal oscillator using GPS signals, including computer equipment, which is programmed or configured to execute the steps of the method for taming the frequency of the crystal oscillator using GPS signals.

In addition, the present invention also provides a system for taming the crystal oscillator frequency by using GPS signals, including a computer device, and a computer program programmed or configured to execute the method for taming the crystal oscillator frequency by using GPS signals is stored in the memory of the computer device.

In addition, the present invention also provides an intelligent terminal with GPS, including a GPS module, a microprocessor and a memory, the microprocessor is programmed or configured to perform the steps of the method for using GPS signals to tame the frequency of the crystal oscillator, or the A computer program programmed or configured to perform the method of taming a frequency of a crystal oscillator using a GPS signal is stored on the memory.

In addition, the present invention also provides a computer-readable storage medium, on which a computer program programmed or configured to execute the method for taming the crystal oscillator frequency by using GPS signals is stored.

Compared with the prior art, the present invention has following advantages: the method that the present invention utilizes GPS signal to tame crystal oscillator frequency mainly comprises four parts: the one, establish the frequency difference mathematics of both according to the 1PPS signal of GPS and the error characteristic of crystal oscillator signal Model. The mathematical modeling of the frequency difference signal provides the basis for establishing the state equation and measurement equation of the frequency difference signal in the next step; the second is to perform conventional extended Kalman filter processing on the frequency difference signal according to the mathematical model of the frequency difference signal, thereby eliminating The random error and the error caused by the non-linear drift of the crystal frequency obtain a filtering model that is greatly affected by the jump outliers and system noise; the third is to introduce a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filter, thereby reducing its The impact on the filtering results of the frequency difference signal; the fourth is to integrate the innovation sequence weighting technology into the above-mentioned extended Kalman filter algorithm, and to reduce the number contained in the GPS output PPS signal by introducing a correction function into the optimal estimation equation of the frequency difference vector. The impact of jumping outliers on filtering results. Therefore, the present invention is based on the innovation sequence weighted extended Kalman filter algorithm of the differential evolution algorithm, which can filter the jump outliers and random errors in the GPS signal, the cumulative error and the nonlinear frequency drift in the crystal oscillator signal, and can improve the accuracy of the GPS signal. Lock the final frequency calibration accuracy of the secondary frequency standard.

Description of drawings

Figure 1 shows the principle of using GPS to lock the secondary frequency standard (high stability crystal oscillator) in the prior art.

Fig. 2 is a schematic flow chart of the basic method of the embodiment of the present invention.

Fig. 3 is a schematic flow chart of introducing a differential evolution algorithm in an embodiment of the present invention.

Detailed ways

The present invention is described in more detail below by means of examples, but the following examples are only illustrative, and the protection scope of the present invention is not limited by these examples.

As shown in Figure 2, the steps of the method for using the GPS signal to tame the crystal oscillator frequency in this embodiment include:

1) According to the error characteristics of the GPS second signal (1PPS signal) and the crystal signal, the mathematical model of the frequency difference between the two is established;

2) Perform conventional extended Kalman filter processing on the frequency difference signal according to the frequency difference mathematical model, thereby eliminating random errors and errors caused by non-linear drift of the crystal oscillator frequency, and obtaining a filtering model that is greatly affected by jump outliers and system noise;

3) Introduce the differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process to reduce its influence on the filtering results of the frequency difference signal;

4) Integrate the weighting of the innovation sequence into the extended Kalman filtering process after the dynamic adjustment of the process noise is completed, so as to reduce the influence of jump outliers on the filtering results of the frequency difference signal.

This embodiment uses the method of GPS signal to tame the crystal oscillator frequency to filter out the noise of the frequency difference signal in the technology of GPS locked oscillator frequency as the goal, and proposes a weighted extended Kalman filter method based on the differential evolution algorithm, which can be used for GPS The jump outliers and random errors in the signal, as well as the random errors and nonlinear frequency drift in the crystal oscillator signal are processed, so as to improve the final frequency calibration accuracy. The method mainly includes four parts: first, according to the error characteristics of the GPS 1PPS signal and the crystal oscillator signal, the mathematical model of the frequency difference between the two is established. The mathematical modeling of the frequency difference signal provides the basis for establishing the state equation and measurement equation of the frequency difference signal in the next step; the second is to perform conventional extended Kalman filter processing on the frequency difference signal according to the mathematical model of the frequency difference signal, thereby eliminating The random error and the error caused by the non-linear drift of the crystal frequency obtain a filtering model that is greatly affected by the jump outliers and system noise; the third is to introduce a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filter, thereby reducing its The impact on the filtering results of the frequency difference signal; the fourth is to integrate the innovation sequence weighting technology into the above-mentioned extended Kalman filtering algorithm. By introducing a correction function into the optimal estimation equation of the frequency difference vector, the influence of the jump wild value included in the GPS output PPS signal on the filtering result is reduced.

In this embodiment, when setting the frequency difference signal and filtering related parameters, the process noise vector ω(k) is set, and its covariance matrix is Q=0.0005. Set the measurement noise vector γ(k), and set its covariance matrix as R=0.01; in addition, the initial value of the frequency difference is set to x ₀ =5×10 ^-9 s.

When the GPS second signal sequence Y is obtained, the universal time UTC of the GPS second signal in the GPS second signal sequence Y (size is n) can be expressed as:

1-ε ₁ , 2-ε ₂ , ..., k-ε _k , ..., n-ε _n

When the crystal oscillator signal sequence Y' is obtained, the universal time UTC of the crystal oscillator signal in the crystal oscillator signal sequence Y' (size n) can be expressed as:

1+a+b+c, 1+a+2b+4c, ..., k+a+bk+ck ² , ..., n+a+bn+cn ²

Wherein, a is the initial phase error, b is an error coefficient considering frequency deviation, and c is an error coefficient considering frequency linear drift; therefore, in this embodiment, the detailed steps of step 1) include:

1.1) obtain the GPS second signal sequence Y, the functional expression of the universal time UTC of the GPS second signal at any k moment in the GPS second signal sequence Y as shown in formula (1);

Y _k = k-ε _k (1)

In the above formula, Y _k represents the universal time UTC of the GPS second signal at time k in the GPS second signal sequence Y, ε _k is the error of the GPS second signal at k time, and satisfies ε～N(0,σ ² );

Obtain the crystal oscillator signal sequence Y ', the functional expression of the universal time UTC of the crystal oscillator signal at any k moment in the crystal oscillator signal sequence Y ' as shown in formula (2);

Y′ _k ＝k+a+bk+ck ² (2)

In the above formula, Y' _k represents the universal time UTC of the crystal oscillator signal at time k in the crystal oscillator signal sequence Y', a is the initial phase error, b is the error coefficient considering the frequency deviation, and c is the error coefficient considering the frequency linear drift ;

1.2) According to the GPS second signal sequence Y, the crystal oscillator signal sequence Y', the difference is obtained to obtain the frequency difference sequence, and the functional expression of the frequency difference signal at any k moment in the frequency difference sequence is shown in the following formula;

X _k ＝Y′ _k -Y _k ＝a+bk+ck ² +ε _k (3)

In the above formula, X _k represents the frequency difference signal at time k in the frequency difference sequence, and the frequency difference sequence is X ₁ , X ₂ , . . . , X _k , . . . , X _n .

In this embodiment, step 2) is used to perform conventional extended Kalman filter processing on the frequency difference signal according to the mathematical model of the frequency difference signal, thereby eliminating random errors and errors caused by non-linear drift of the crystal oscillator frequency, and obtaining the outlier value affected by the jump and filtering models that have a greater influence on system noise. In step 2), the steps of carrying out conventional extended Kalman filter processing to the frequency difference signal according to the frequency difference mathematical model include:

2.1) The functional expressions of the state equation and the observation equation of the frequency difference are established as shown in formula (4);

In the above formula, x(k+1) represents the frequency difference state vector at time k+1, x(k) represents the frequency difference state vector at time k, f and h represent the state function and measurement function respectively, z(k+ 1) represents the measurement value at time k+1; ω(k) represents the process noise vector at time k, and its covariance matrix is Q; γ(k) represents the measurement noise vector at time k, and its covariance matrix is R ;

2.2) Let the functional expressions of the state transition matrix F(k+1,k) from k time to k+1 time and the observation matrix H(k) at k time be shown in formula (5);

In the above formula, f represents the state function, Represents the estimation of the frequency difference state vector at time k, H(k) represents the observation matrix at time k, h represents the measurement function, x(k) represents the frequency difference state vector at time k;

2.3) The prediction equations for establishing the conventional extended Kalman filter model are shown in equations (6) to (8), the update equations are shown in equations (9) to (10), and the process noise is shown in equation (11). The noise variance is shown in formula (12);

In formula (6), Represents the predicted value of the frequency difference state vector at time k+1 based on the predicted value at time k, f represents the state function, /> Represents the optimal estimated value of the frequency difference state vector at time k;

p(k+1,k)=F(k+1,k)p(k,k)F ^T (k+1,k)+Q(k) (7)

In formula (7), p(k+1,k) represents the predicted value of the covariance vector at time k+1, F(k+1,k) represents the state transition matrix from time k to time k+1, p(k , k) represents the predicted value of the covariance vector at time k, F ^T (k+1,k) represents the transposition of the state transition matrix from time k to time k+1, Q(k) represents the process noise variance at time k;

K(k+1)=p(k+1,k)H ^T (k+1)[H(k+1)p(k+1,k)H ^T (k+1)+R(k+1 )] ^-1 (8)

In formula (8), K(k+1) represents the Kalman filter gain matrix at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, H ^T (k+1) Represents the transposition of the observation matrix H(k+1) at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, R(k+1) represents the predicted value of the covariance vector at time k+1 Observation noise variance vector;

In formula (9), Indicates the optimal estimated value of the frequency difference state vector at time k+1, /> Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, z(k+1) represents the measured value at time k+1, /> Indicates the predicted value of the frequency difference state vector at time k+1 /> The obtained k+1 time measurement vector prediction value;

p(k+1)=[I-K(k+1)H(k+1)]p(k+1,k) (10)

In formula (10), p(k+1) represents the covariance vector estimated value updated at time k+1, I represents the identity matrix with the same dimension as K(k+1)H(k+1), and K(k +1) represents the Kalman filter gain matrix at time k+1, H(k+1) represents the observation matrix at time k+1, and p(k+1,k) represents the predicted value of the covariance vector at time k+1;

In formula (11), ω(k+1) represents the process noise at time k+1, Represents the optimal estimated value of the frequency difference state vector at time k+1, and z(k+1) represents the measured value at time k+1;

Q _k+1 = var(ω(k+1)) (12)

In formula (12), Q _k+1 represents the variance of process noise at time k+1, ω(k+1) represents the process noise at time k+1, and var represents the operation of calculating the variance of process noise.

In this embodiment, step 3) introduces a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filter algorithm, thereby reducing its influence on the filtering result of the frequency difference signal. Conventional extended Kalman generally assumes that the mean value of the process noise is zero and the variance is Q in the filtering process. Due to the uncertain factors and abnormal disturbances in the system model, the process noise also changes; the fixed process noise will increase the error and reduce the filtering accuracy. To this end, the differential evolution algorithm is introduced to solve the optimal variance of the process noise, and then it is applied to the extended Kalman filter algorithm to improve the filtering accuracy and finally improve the frequency calibration accuracy. As shown in Figure 3, when the differential evolution algorithm is introduced in step 3) to dynamically adjust the process noise in the conventional extended Kalman filter processing, the steps of generating the variance of the process noise at any time k by the differential evolution algorithm include:

3.1) Initialize the population, and the function expression is shown in the following formula (13):

x _i,k (0)＝l _k +rand(0,1)*(u _k -l _k ) (13)

In formula (13), x _i,k (0) represents the kth gene of the i-th "chromosome" of the 0th generation in the population; rand(0,1) is a random number uniformly distributed between 0 and 1; u _k and l _k are the upper bound and lower bound of the search, and u _k and l _k respectively take the maximum value Q _max and the minimum value Q _min of the process noise variance; when setting the initial value of the differential evolution algorithm in this embodiment, the mutation rate F=0.5, crossover probability cr=0.9, population size NP=10, iteration number G=100, search boundary u _k =0.005, l _k =-0.005;

3.2) The mutation operation is performed on the population, and the function expression of the mutation operation is shown in the following formula (14):

X _i (g+1)=x _r1 (g)+F*[x _r2 (g)-x _r3 (g)] (14)

In formula (14), _Xi (g+1) is the obtained g+1 generation mutant individual; F is the compression scale factor, the value range is 0~1; x _r1 (g), x _r2 (g) and x _r3 (g) being 3 parents;

3.3) Perform a crossover operation on the population to retain better variables, and the crossover operation adopts the binomial crossover method, and its function expression is shown in the following formula (15):

In formula (15), y _i,j (g+1) represents the crossover operation result corresponding to the current individual, Xi _,j (g+1) represents the component corresponding to the mutant individual, and xi _,j (g) represents the current individual The corresponding component, r is a random number uniformly distributed between 0 and 1 generated by each variable; cr is the crossover probability of the variable; rnd is an integer uniformly distributed between 1 and d, if r<cr, the mutant individual is accepted The corresponding component X _i,j (g+1), otherwise keep the component x _i,j (g) corresponding to the current individual;

3.4) For the results obtained by the crossover operation, the optimal solution is selected as the process noise variance at time k by using the greedy selection method shown in the following formula (16);

In formula (16), x _i (g+1) is the selection result, y _i (g+1) represents the individual corresponding to the crossover operation result, _xi (g) represents the current individual, f[y _i (g+1) ] represents the population fitness value corresponding to individual y _i (g+1), and f[ _xi (g)] represents the population fitness value corresponding to individual _xi (g). By assigning the optimal solution to Q _i and substituting it into the covariance prediction equation, the filtering accuracy of Extended Kalman can be improved.

In this embodiment, step 4) integrates the innovation sequence weighting technology into the above-mentioned extended Kalman filter algorithm, so as to reduce the impact of the jump outliers contained in the GPS output PPS signal on the filter result. In practical applications, because the GPS signal is interfered by noise, it will also produce large jump outliers, and the frequency difference signal with outliers will cause the filtering results to shift or even diverge. In order to solve the anti-outlier problem of the extended Kalman algorithm mentioned above, some scholars introduce the method of innovation sequence weighting into the Kalman filter to eliminate the influence of outliers. When the observation data does not contain outliers, the filtering algorithm can run normally and filter out the observation noise; when the observation data contains outliers, it can overcome or reduce the influence of the outliers on the filtering results, making the filtering results more accurate. In this embodiment, after step 4) integrates the weighting of the innovation sequence into the extended Kalman filter processing after the process noise is dynamically adjusted, the functional expression of the update equation of the extended Kalman filter processing obtained is as follows (17)～ (21);

D _k+1 ＝P(k+1,k)+R _k+1 (21)

In formula (17)～(21), Indicates the optimal estimated value of the frequency difference state vector at time k+1 after adjustment, /> Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, φ _k+1 represents the correction function, and τ _k+1 represents the measurement error at time k+1 The variable of size, e _k+1 represents the measurement error at k+1 time, z(k+1) represents the measured value of the frequency difference vector at k+1 time, /> Indicates the predicted value of the measurement vector at time k+1, τ _k+1 obeys the χ ² distribution, when an outlier occurs, τ _k+1 increases and φ _k+1 ( ) decreases; C _k+1 is the selected Threshold value or constant sequence, λ _k+1 is The largest eigenvalue of the matrix, K _k+1 represents the Kalman filter gain at time k+1, D _k+1 represents the weight matrix that measures the measurement error at time k+1, and P(k+1,k) represents k+ The covariance prediction value at time 1, R(k+1) represents the measurement noise variance at time k+1.

To sum up, the main technology traditionally applied to low-cost and high-precision frequency standards is GPS locking secondary frequency standard (crystal oscillator) technology. The influence of the nonlinear frequency offset of the crystal oscillator signal on the frequency difference signal. This embodiment provides a method for processing the frequency difference signal in the GPS-locked crystal oscillator frequency technology based on the differential evolution algorithm based on the weighted extended Kalman filter of the innovation sequence. The method of this embodiment mainly includes four parts: one is based on GPS Based on the error characteristics of the 1PPS signal and the crystal oscillator signal, a mathematical model of the frequency difference between the two is established. The mathematical modeling of the frequency difference signal provides the basis for establishing the state equation and measurement equation of the frequency difference signal in the next step; the second is to perform conventional extended Kalman filter processing on the frequency difference signal according to the mathematical model of the frequency difference signal, thereby eliminating The random error and the error caused by the non-linear drift of the crystal frequency obtain a filtering model that is greatly affected by the jump outliers and system noise; the third is to introduce a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filter, thereby reducing its The impact on the filtering results of the frequency difference signal; the fourth is to integrate the innovation sequence weighting technology into the above-mentioned extended Kalman filtering algorithm. By introducing a correction function into the optimal estimation equation of the frequency difference vector, the influence of the jump wild value included in the GPS output PPS signal on the filtering result is reduced. Moreover, the method of this embodiment adopts the extended Kalman filter algorithm to process the frequency difference signal in the process of GPS locking the frequency of the crystal oscillator. While processing the random error of the GPS signal, it also takes into account the abnormality caused by the aging and temperature of the crystal oscillator. linearity error. The method of this embodiment introduces the innovation sequence weighting and the differential evolution algorithm into the extended Kalman filtering process, which can reduce the influence of large jump outliers in the GPS signal and continuous changes of system noise on the filtering results.

In addition, this embodiment also provides a system for using GPS signals to tame the frequency of the crystal oscillator, including:

A frequency difference model generation program unit is used to set up both frequency difference mathematical models according to the error characteristics of the GPS second signal and the crystal oscillator signal;

The conventional filter program unit is used to perform conventional extended Kalman filter processing on the frequency difference signal according to the frequency difference mathematical model, thereby eliminating random errors and errors caused by non-linear drift of the crystal oscillator frequency, and obtaining a relatively low frequency value affected by jump outliers and system noise. Large filtering model;

The noise adjustment program unit is used to introduce a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process to reduce its influence on the frequency difference signal filtering result;

The filter integration program unit is used to integrate the weight of the innovation sequence into the extended Kalman filter processing after the process noise is dynamically adjusted.

In addition, this embodiment also provides a system for taming the frequency of the crystal oscillator by using GPS signals, including a computer device programmed or configured to execute the steps of the method for taming the frequency of the crystal oscillator by using GPS signals.

In addition, this embodiment also provides a system for taming the frequency of the crystal oscillator using GPS signals, including a computer device, and a computer program programmed or configured to execute the aforementioned method for taming the frequency of the crystal oscillator using GPS signals is stored in the memory of the computer device.

In addition, this embodiment also provides an intelligent terminal with GPS, including a GPS module, a microprocessor and a memory, the microprocessor is programmed or configured to perform the steps of the aforementioned method of using GPS signals to tame the frequency of the crystal oscillator, or the A computer program programmed or configured to execute the aforementioned method of taming the frequency of the crystal oscillator using GPS signals is stored on the memory. The smart terminal may be a smart phone, a tablet computer or other types of terminal devices.

In addition, this embodiment also provides a computer-readable storage medium, on which a computer program programmed or configured to execute the aforementioned method for taming the frequency of the crystal oscillator by using a GPS signal is stored.

Those skilled in the art should understand that the embodiments of the present application may be provided as methods, systems, or computer program products. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) having computer-usable program code embodied therein. The present application refers to the flowcharts of the methods, devices (systems), and computer program products according to the embodiments of the present application and/or the instructions executed by the processor are used to implement one or more procedures in the flowchart and/or one of the block diagrams means for the function specified in the box or boxes. These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing apparatus to operate in a specific manner, such that the instructions stored in the computer-readable memory produce an article of manufacture comprising instruction means, the instructions The device realizes the function specified in one or more procedures of the flowchart and/or one or more blocks of the block diagram. These computer program instructions can also be loaded onto a computer or other programmable data processing device, causing a series of operational steps to be performed on the computer or other programmable device to produce a computer-implemented process, thereby The instructions provide steps for implementing the functions specified in the flow chart or blocks of the flowchart and/or the block or blocks of the block diagrams.

The above descriptions are only preferred implementations of the present invention, and the scope of protection of the present invention is not limited to the above examples, and all technical solutions that fall under the idea of the present invention belong to the scope of protection of the present invention. It should be pointed out that for those skilled in the art, some improvements and modifications without departing from the principles of the present invention should also be regarded as the protection scope of the present invention.

Claims (8)

1. a method utilizing GPS signal to tame crystal oscillator frequency, is characterized in that, the step of this method comprises: 1) According to the error characteristics of the GPS second signal and the crystal signal, the mathematical model of the frequency difference between the two is established; 2) Carry out conventional extended Kalman filter processing to the frequency difference signal according to the frequency difference mathematical model, thereby eliminating the error caused by the random error and the non-linear drift of the crystal oscillator frequency; The steps of Mann filter processing include: 2.1) The functional expressions of the state equation and the observation equation of the frequency difference are established as shown in formula (4); In the above formula, x(k+1) represents the frequency difference state vector at time k+1, x(k) represents the frequency difference state vector at time k, f and h represent the state function and measurement function respectively, z(k+ 1) represents the measurement value at time k+1; ω(k) represents the process noise vector at time k, and its covariance matrix is Q; γ(k) represents the measurement noise vector at time k, and its covariance matrix is R ; 2.2) Let the functional expressions of the state transition matrix F(k+1,k) from k time to k+1 time and the observation matrix H(k) at k time be shown in formula (5); In the above formula, f represents the state function, Represents the estimation of the frequency difference state vector at time k, h represents the measurement function, and x(k) represents the frequency difference state vector at time k; 2.3) The prediction equations for establishing the conventional extended Kalman filter model are shown in equations (6) to (8), the update equations are shown in equations (9) to (10), and the process noise is shown in equation (11). The noise variance is shown in formula (12); In formula (6), Represents the predicted value of the frequency difference state vector at time k+1 based on the predicted value at time k, f represents the state function, /> Represents the optimal estimated value of the frequency difference state vector at time k; p(k+1,k)=F(k+1,k)p(k,k)F ^T (k+1,k)+Q(k) (7) In formula (7), p(k+1,k) represents the predicted value of the covariance vector at time k+1, F(k+1,k) represents the state transition matrix from time k to time k+1, p(k , k) represents the predicted value of the covariance vector at time k, F ^T (k+1,k) represents the transposition of the state transition matrix from time k to time k+1, Q(k) represents the process noise variance at time k; K(k+1)=p(k+1,k)H ^T (k+1)[H(k+1)p(k+1,k)H ^T (k+1)+R(k+1 )] ^-1 (8) In formula (8), K(k+1) represents the Kalman filter gain matrix at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, H ^T (k+1) Represents the transposition of the observation matrix H(k+1) at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, R(k+1) represents the predicted value of the covariance vector at time k+1 Observation noise variance vector; In formula (9), Indicates the optimal estimated value of the frequency difference state vector at time k+1, /> Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, z(k+1) represents the measured value at time k+1, /> Indicates the predicted value of the frequency difference state vector at time k+1 /> The obtained k+1 time measurement vector prediction value; p(k+1)=[I-K(k+1)H(k+1)]p(k+1,k) (10) In formula (10), p(k+1) represents the covariance vector estimated value updated at time k+1, I represents the identity matrix with the same dimension as K(k+1)H(k+1), and K(k +1) represents the Kalman filter gain matrix at time k+1, H(k+1) represents the observation matrix at time k+1, and p(k+1,k) represents the predicted value of the covariance vector at time k+1; In formula (11), ω(k+1) represents the process noise at time k+1, Represents the optimal estimated value of the frequency difference state vector at time k+1, and z(k+1) represents the measured value at time k+1; Q _k+1 = var(ω(k+1)) (12) In formula (12), Q _k+1 represents the variance of process noise at time k+1, ω(k+1) represents the process noise at time k+1, and var represents the operation of calculating the variance of process noise; 3) Introduce the differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process to reduce its influence on the filtering results of the frequency difference signal; 4) Integrate the weighting of the innovation sequence into the extended Kalman filter processing after the dynamic adjustment of the process noise to reduce the influence of jump outliers on the filtering results of the frequency difference signal; and integrate the weighting of the innovation sequence into the dynamic adjustment of the process noise After the extended Kalman filter processing, the functional expression of the update equation of the extended Kalman filter processing obtained is as shown in formula (17)～(21); D _k+1 ＝P(k+1,k)+R _k+1 (21) In formula (17)～(21), Represents the optimal estimated value of the frequency difference state vector at time k+1 after adjustment, Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, φ _k+1 represents the correction function, and τ _k+1 represents the measurement error at time k+1 The variable of size, e _k+1 represents the measurement error at k+1 time, z(k+1) represents the measured value of the frequency difference vector at k+1 time, /> Indicates the predicted value of the measurement vector at time k+1, τ _k+1 obeys the χ ² distribution, when an outlier occurs, τ _k+1 increases and φ _k+1 ( ) decreases; C _k+1 is the selected Threshold value or constant sequence, λ _k+1 is The largest eigenvalue of the matrix, K _k+1 represents the Kalman filter gain at time k+1, D _k+1 represents the weight matrix that measures the measurement error at time k+1, and P(k+1,k) represents k+ The covariance prediction value at time 1, R _k+1 represents the measurement noise variance at time k+1. 2. the method utilizing GPS signal to tame crystal frequency according to claim 1, is characterized in that, the detailed steps of step 1) comprise: 1.1) obtain the GPS second signal sequence Y, the functional expression of the universal time UTC of the GPS second signal at any k moment in the GPS second signal sequence Y as shown in formula (1); Y _k = k-ε _k (1) In the above formula, Y _k represents the universal time UTC of the GPS second signal at time k in the GPS second signal sequence Y, ε _k is the error of the GPS second signal at k time, and satisfies ε～N(0,σ ² ); Obtain the crystal oscillator signal sequence Y ', the functional expression of the universal time UTC of the crystal oscillator signal at any k moment in the crystal oscillator signal sequence Y ' as shown in formula (2); Y′ _k ＝k+a+bk+ck ² (2) In the above formula, Y' _k represents the universal time UTC of the crystal oscillator signal at time k in the crystal oscillator signal sequence Y', a is the initial phase error, b is the error coefficient considering the frequency deviation, and c is the error coefficient considering the frequency linear drift ; 1.2) According to the GPS second signal sequence Y, the crystal oscillator signal sequence Y', the difference is obtained to obtain the frequency difference sequence, and the functional expression of the frequency difference signal at any k moment in the frequency difference sequence is shown in the following formula; X _k ＝Y′ _k -Y _k ＝a+bk+ck ² +ε _k (3) In the above formula, X _k represents the frequency difference signal at time k in the frequency difference sequence. 3. the method utilizing GPS signal to tame the crystal oscillator frequency according to claim 1, is characterized in that, in step 3) introduces differential evolution algorithm when the process noise in conventional Extended Kalman filter process is dynamically adjusted, by differential evolution The steps of the algorithm to generate the process noise variance at any time k include: 3.1) Initialize the population, and the function expression is shown in the following formula (13): x _i,k (0)＝l _k +rand(0,1)*(u _k -l _k ) (13) In formula (13), x _i,k (0) represents the kth gene of the i-th "chromosome" of the 0th generation in the population; rand(0,1) is a random number uniformly distributed between 0 and 1; u _k and l _k are the upper bound and lower bound of the search, and u _k and l _k respectively take the maximum value Q _max and the minimum value Q _min of the process noise variance; 3.2) The mutation operation is performed on the population, and the function expression of the mutation operation is shown in the following formula (14): X _i (g+1)=x _r1 (g)+F*[x _r2 (g)-x _r3 (g)] (14) In formula (14), _Xi (g+1) is the obtained g+1 generation mutant individual; F is the compression scale factor, the value range is 0~1; x _r1 (g), x _r2 (g) and x _r3 (g) being 3 parents; 3.3) Perform a crossover operation on the population to retain better variables, and the crossover operation adopts the binomial crossover method, and its function expression is shown in the following formula (15): In formula (15), y _i,j (g+1) represents the crossover operation result corresponding to the current individual, Xi _,j (g+1) represents the component corresponding to the mutant individual, and xi _,j (g) represents the current individual The corresponding component, r is a random number uniformly distributed between 0 and 1 generated by each variable; cr is the crossover probability of the variable; rnd is an integer uniformly distributed between 1 and d, if r<cr, the mutant individual is accepted The corresponding component X _i,j (g+1), otherwise keep the component x _i,j (g) corresponding to the current individual; 3.4) For the results obtained by the crossover operation, the optimal solution is selected as the process noise variance at time k by using the greedy selection method shown in the following formula (16); In formula (16), x _i (g+1) is the selection result, y _i (g+1) represents the individual corresponding to the crossover operation result, _xi (g) represents the current individual, f[y _i (g+1) ] represents the population fitness value corresponding to individual y _i (g+1), and f[ _xi (g)] represents the population fitness value corresponding to individual _xi (g). 4. A system utilizing GPS signals to tame the crystal oscillator frequency, characterized in that it comprises: A frequency difference model generation program unit is used to set up both frequency difference mathematical models according to the error characteristics of the GPS second signal and the crystal oscillator signal; The conventional filter program unit is used to perform conventional extended Kalman filter processing on the frequency difference signal according to the frequency difference mathematical model, thereby eliminating random errors and errors caused by the non-linear drift of the crystal oscillator frequency; the frequency difference signal is processed according to the frequency difference mathematical model The steps for conventional extended Kalman filter processing include: 2.1) The functional expressions of the state equation and the observation equation of the frequency difference are established as shown in formula (4); In the above formula, x(k+1) represents the frequency difference state vector at time k+1, x(k) represents the frequency difference state vector at time k, f and h represent the state function and measurement function respectively, z(k+ 1) represents the measurement value at time k+1; ω(k) represents the process noise vector at time k, and its covariance matrix is Q; γ(k) represents the measurement noise vector at time k, and its covariance matrix is R ; 2.2) Let the functional expressions of the state transition matrix F(k+1,k) from k time to k+1 time and the observation matrix H(k) at k time be shown in formula (5); In the above formula, f represents the state function, Represents the estimation of the frequency difference state vector at time k, h represents the measurement function, and x(k) represents the frequency difference state vector at time k; 2.3) The prediction equations for establishing the conventional extended Kalman filter model are shown in equations (6) to (8), the update equations are shown in equations (9) to (10), and the process noise is shown in equation (11). The noise variance is shown in formula (12); In formula (6), Represents the predicted value of the frequency difference state vector at time k+1 based on the predicted value at time k, f represents the state function, /> Represents the optimal estimated value of the frequency difference state vector at time k; p(k+1,k)=F(k+1,k)p(k,k)F ^T (k+1,k)+Q(k) (7) In formula (7), p(k+1,k) represents the predicted value of the covariance vector at time k+1, F(k+1,k) represents the state transition matrix from time k to time k+1, p(k , k) represents the predicted value of the covariance vector at time k, F ^T (k+1,k) represents the transposition of the state transition matrix from time k to time k+1, Q(k) represents the process noise variance at time k; K(k+1)=p(k+1,k)H ^T (k+1)[H(k+1)p(k+1,k)H ^T (k+1)+R(k+1 )] ^-1 (8) In formula (8), K(k+1) represents the Kalman filter gain matrix at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, H ^T (k+1) Represents the transposition of the observation matrix H(k+1) at time k+1, p(k+1,k) represents the predicted value of the covariance vector at time k+1, R(k+1) represents the predicted value of the covariance vector at time k+1 Observation noise variance vector; In formula (9), Indicates the optimal estimated value of the frequency difference state vector at time k+1, /> Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, z(k+1) represents the measured value at time k+1, /> Indicates the predicted value of the frequency difference state vector at time k+1 /> The obtained k+1 time measurement vector prediction value; p(k+1)=[I-K(k+1)H(k+1)]p(k+1,k) (10) In formula (10), p(k+1) represents the covariance vector estimated value updated at time k+1, I represents the identity matrix with the same dimension as K(k+1)H(k+1), and K(k +1) represents the Kalman filter gain matrix at time k+1, H(k+1) represents the observation matrix at time k+1, and p(k+1,k) represents the predicted value of the covariance vector at time k+1; In formula (11), ω(k+1) represents the process noise at time k+1, Represents the optimal estimated value of the frequency difference state vector at time k+1, and z(k+1) represents the measured value at time k+1; Q _k+1 = var(ω(k+1)) (12) In formula (12), Q _k+1 represents the variance of process noise at time k+1, ω(k+1) represents the process noise at time k+1, and var represents the operation of calculating the variance of process noise; The noise adjustment program unit is used to introduce a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process to reduce its influence on the frequency difference signal filtering result; The filter integration program unit is used to integrate the weighting of the innovation sequence into the extended Kalman filter processing after the dynamic adjustment of the process noise, so as to reduce the influence of the jump wild value on the filtering result of the frequency difference signal; and integrate the weighting of the innovation sequence into the completion After the process noise is dynamically adjusted in the extended Kalman filtering process, the functional expression of the update equation of the extended Kalman filtering process obtained is shown in the following formulas (17)-(21); D _k+1 ＝P(k+1,k)+R _k+1 (21) In formula (17)～(21), Represents the optimal estimated value of the frequency difference state vector at time k+1 after adjustment, Represents the predicted value of the frequency difference state vector at time k+1, K(k+1) represents the Kalman filter gain matrix at time k+1, φ _k+1 represents the correction function, and τ _k+1 represents the measurement error at time k+1 The variable of size, e _k+1 represents the measurement error at k+1 time, z(k+1) represents the measured value of the frequency difference vector at k+1 time, /> Indicates the predicted value of the measurement vector at time k+1, τ _k+1 obeys the χ ² distribution, when an outlier occurs, τ _k+1 increases and φ _k+1 ( ) decreases; C _k+1 is the selected Threshold value or constant sequence, λ _k+1 is The largest eigenvalue of the matrix, K _k+1 represents the Kalman filter gain at time k+1, D _k+1 represents the weight matrix that measures the measurement error at time k+1, and P(k+1,k) represents k+ The covariance prediction value at time 1, R _k+1 represents the measurement noise variance at time k+1. 5. A system utilizing GPS signals to tame the frequency of crystal oscillators, comprising computer equipment, characterized in that the computer equipment is programmed or configured to perform the method for taming crystal oscillator frequencies using GPS signals according to any one of claims 1 to 3 step. 6. A system for using GPS signals to tame the frequency of the crystal oscillator, comprising computer equipment, characterized in that the memory of the computer equipment is programmed or configured to perform taming using GPS signals according to any one of claims 1 to 3. A computer program for methods of crystal oscillator frequencies. 7. An intelligent terminal with GPS, comprising a GPS module, a microprocessor and a memory, characterized in that the microprocessor is programmed or configured to perform any one of claims 1 to 3 using GPS signals to tame The steps of the method for crystal oscillator frequency, or the computer program programmed or configured to execute the method for taming the crystal oscillator frequency by using GPS signals in any one of claims 1-3 is stored on the memory. 8. A computer-readable storage medium, characterized in that, the computer-readable storage medium is stored with a computer programmed or configured to perform the method for using GPS signals to tame the frequency of the crystal oscillator according to any one of claims 1 to 3 program.