CN111711446B - Method, system and medium for taming crystal oscillator frequency by using GPS signal - Google Patents
Method, system and medium for taming crystal oscillator frequency by using GPS signal Download PDFInfo
- Publication number
- CN111711446B CN111711446B CN202010523171.8A CN202010523171A CN111711446B CN 111711446 B CN111711446 B CN 111711446B CN 202010523171 A CN202010523171 A CN 202010523171A CN 111711446 B CN111711446 B CN 111711446B
- Authority
- CN
- China
- Prior art keywords
- time
- frequency difference
- formula
- vector
- crystal oscillator
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 159
- 239000013078 crystal Substances 0.000 title claims abstract description 104
- 238000001914 filtration Methods 0.000 claims abstract description 115
- 230000008569 process Effects 0.000 claims abstract description 114
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 36
- 238000013178 mathematical model Methods 0.000 claims abstract description 25
- 238000012545 processing Methods 0.000 claims abstract description 25
- 238000012937 correction Methods 0.000 claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims description 62
- 230000006870 function Effects 0.000 claims description 58
- 238000005259 measurement Methods 0.000 claims description 51
- 108010076504 Protein Sorting Signals Proteins 0.000 claims description 28
- 238000004590 computer program Methods 0.000 claims description 14
- 230000007704 transition Effects 0.000 claims description 12
- 125000003275 alpha amino acid group Chemical group 0.000 claims description 8
- 150000001875 compounds Chemical class 0.000 claims description 8
- 230000035772 mutation Effects 0.000 claims description 7
- 238000003860 storage Methods 0.000 claims description 7
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 238000009826 distribution Methods 0.000 claims description 4
- 210000000349 chromosome Anatomy 0.000 claims description 3
- 230000006835 compression Effects 0.000 claims description 3
- 238000007906 compression Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims description 3
- 108090000623 proteins and genes Proteins 0.000 claims description 3
- 238000012549 training Methods 0.000 claims 1
- 238000005516 engineering process Methods 0.000 description 10
- 238000010586 diagram Methods 0.000 description 5
- 230000007774 longterm Effects 0.000 description 3
- 230000032683 aging Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 229910052792 caesium Inorganic materials 0.000 description 1
- TVFDJXOCXUVLDH-UHFFFAOYSA-N caesium atom Chemical compound [Cs] TVFDJXOCXUVLDH-UHFFFAOYSA-N 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 229910052701 rubidium Inorganic materials 0.000 description 1
- IGLNJRXAVVLDKE-UHFFFAOYSA-N rubidium atom Chemical compound [Rb] IGLNJRXAVVLDKE-UHFFFAOYSA-N 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
Classifications
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03L—AUTOMATIC CONTROL, STARTING, SYNCHRONISATION OR STABILISATION OF GENERATORS OF ELECTRONIC OSCILLATIONS OR PULSES
- H03L7/00—Automatic control of frequency or phase; Synchronisation
- H03L7/06—Automatic control of frequency or phase; Synchronisation using a reference signal applied to a frequency- or phase-locked loop
- H03L7/16—Indirect frequency synthesis, i.e. generating a desired one of a number of predetermined frequencies using a frequency- or phase-locked loop
- H03L7/18—Indirect frequency synthesis, i.e. generating a desired one of a number of predetermined frequencies using a frequency- or phase-locked loop using a frequency divider or counter in the loop
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03L—AUTOMATIC CONTROL, STARTING, SYNCHRONISATION OR STABILISATION OF GENERATORS OF ELECTRONIC OSCILLATIONS OR PULSES
- H03L7/00—Automatic control of frequency or phase; Synchronisation
- H03L7/06—Automatic control of frequency or phase; Synchronisation using a reference signal applied to a frequency- or phase-locked loop
- H03L7/08—Details of the phase-locked loop
- H03L7/085—Details of the phase-locked loop concerning mainly the frequency- or phase-detection arrangement including the filtering or amplification of its output signal
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D30/00—Reducing energy consumption in communication networks
- Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks
Landscapes
- Oscillators With Electromechanical Resonators (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
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
Technical Field
The invention relates to a frequency source error filtering technology of a GPS (Global positioning System) tame crystal oscillator signal, in particular to a method, a system and a medium for taming the crystal oscillator frequency by utilizing the GPS signal, which are used for realizing high-precision tame of the GPS signal to the crystal oscillator frequency.
Background
With the rapid development of network technology, the requirements of the system on time references are also increasing. The time reference in the power system has important significance for accident analysis, electric energy quality improvement, power grid dispatching operation and electric energy time-sharing charging.
The current time (frequency) reference can be divided into a primary frequency standard (cesium atomic frequency standard), a secondary frequency standard (including rubidium atomic frequency standard and high-stability crystal oscillator) and other frequency standards (other crystal oscillators unexpected to high-stability crystal oscillator) according to application fields and performance indexes. Although the primary frequency standard has higher long-term stability and accuracy, the primary frequency standard is quite expensive and is difficult to apply to civil occasions with strict requirements on cost; the price of the secondary frequency scale is relatively low, and the long-term stability and accuracy are not as good as the former, so that the secondary frequency scale is difficult to apply to high-precision occasions.
The principle of the prior art for locking a secondary frequency scale (high stable crystal oscillator) by using a GPS is shown in the accompanying figure 1: firstly, measuring the relative frequency difference between a second pulse signal PPS output by a GPS receiver and a pulse signal output after frequency division of a crystal oscillator, then carrying out filtering processing on the frequency difference signal, and finally generating a control signal, and continuously adjusting the frequency of the crystal oscillator so as to output frequency output with high accuracy and high stability. The technique solves the above-mentioned problems to a certain extent, namely, achieving a frequency standard of high long-term stability and average accuracy with lower cost. However, in practical application, frequency drift caused by random error, jump wild value, aging, temperature and other characteristics of the 1PPS signal of the GPS causes certain difficulty in frequency calibration.
In order to solve the above problems, some scholars propose to process the measured frequency difference data by adopting a method of sampling and averaging for a plurality of times, and the calculation is simple but the frequency correction time is long. The error of the frequency difference mathematical model between the crystal oscillator signal and the GPS signal is calculated and compensated by a scholars through a least square method, but the method is complex in calculation and does not consider the influence of the wild value of the GPS signal on the frequency correction precision. Still other scholars propose to eliminate noise in the frequency difference signal by using a kalman filtering algorithm, but the existing mode of eliminating noise in the frequency difference signal by using the kalman filtering algorithm does not consider the influence of crystal oscillator frequency drift, so that the frequency correction accuracy is limited.
Disclosure of Invention
The invention aims to solve the technical problems: aiming at the problems in the prior art, the invention provides a method, a system and a medium for taming the frequency of a crystal oscillator by utilizing a GPS signal.
In order to solve the technical problems, the invention adopts the following technical scheme:
a method for taming crystal oscillator frequency using GPS signals, the method comprising the steps of:
1) Establishing a frequency difference mathematical model of the GPS second signal and the crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal;
2) Performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model, so as to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency and obtain a filtering model greatly influenced by a jump wild value and system noise;
3) 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;
4) And integrating the new sequence weighting 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 of the frequency difference signal.
Optionally, the detailed steps of step 1) include:
1.1 Acquiring a GPS second signal sequence Y, wherein a function expression of universal time UTC of GPS second signals at any k moment in the GPS second signal sequence Y is shown in a formula (1);
Y k ＝k-ε k (1)
in the above, Y k Universal time UTC, epsilon for GPS second signal representing time k in GPS second signal sequence Y k Is the error of the GPS second signal at the k moment and satisfies epsilon-N (0, sigma) 2 )；
Acquiring a crystal oscillator signal sequence Y ', wherein a function expression of universal time UTC of crystal oscillator signals at any k moment in the crystal oscillator signal sequence Y' is shown in a formula (2);
Y′ k ＝k+a+bk+ck 2 (2)
in the above, Y' k The universal time UTC of the crystal oscillator signal at the moment k in the crystal oscillator signal sequence Y', a is an initial phase error, b is an error coefficient considering frequency deviation, and c is an error coefficient considering frequency linear drift;
1.2 According to the GPS second signal sequence Y and the crystal oscillator signal sequence Y', obtaining a frequency difference sequence, wherein the function 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 2 +ε k (3)
in the above, X k And represents the frequency difference signal at time k in the frequency difference sequence.
Optionally, the step of performing conventional extended kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model in step 2) includes:
2.1 A function expression of a state equation and an observation equation of the frequency difference is established as shown in a formula (4);
in the above formula, x (k+1) represents a frequency difference state vector at k+1, x (k) represents a frequency difference state vector at k, f and h represent a state function and a measurement function, respectively, and z (k+1) represents a measurement value at k+1; omega (k) represents a process noise vector at the moment k, and the covariance matrix of the process noise vector is Q; gamma (k) represents a measurement noise vector at the moment k, and the covariance matrix of the measurement noise vector is R;
2.2 Let the function expression of the state transition matrix F (k+1, k) from the k moment to the k+1 moment and the observation matrix H (k) from the k moment be shown as the formula (5);
in the above formula, f represents a state function,the estimation of the frequency difference state vector at the moment k is represented, h represents a measurement function, and x (k) represents the frequency difference state vector at the moment k;
2.3 Establishing a predictive equation of a conventional extended Kalman filtering model as shown in formulas (6) to (8), and an update equation as shown in formulas (9) to (10), wherein the process noise is shown in formula (11), and the process noise variance is shown in formula (12);
in the formula (6), the amino acid sequence of the compound,a k+1 time frequency difference state vector predicted value based on the k time predicted value, f represents a state function,/->Representing the optimal estimated value of the frequency difference state vector at the k moment;
p(k+1,k)＝F(k+1,k)p(k,k)F T (k+1,k)+Q(k) (7)
in the formula (7), p (k+1, k) represents a covariance vector predictor at time k+1, F (k+1, k) represents a state transition matrix from time k to time k+1, p (k, k) represents a covariance vector predictor at time k, F T (k+1, k) represents the transpose of the state transition matrix from time k to time k+1, and 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 the formula (8), K (k+1) represents a Kalman filtering gain matrix at the time of k+1, p (k+1, K) represents a covariance vector predicted value at the time of k+1, and H T (k+1) represents a transpose of the observation matrix H (k+1) at time k+1, p (k+1, k) represents a covariance vector predictor at time k+1, and R (k+1) represents an observation noise variance vector at time k+1;
In the formula (9), the amino acid sequence of the compound,optimal estimate representing the state vector of the frequency difference at time k+1,/and>a predicted value of a frequency difference state vector at time k+1, K (k+1) represents a Kalman filtering gain matrix at time k+1, z (k+1) represents a measured value at time k+1, +.>Representing the state vector predictor +.>The obtained k+1 moment measurement vector predicted value;
p(k+1)＝[I-K(k+1)H(k+1)]p(k+1,k) (10)
in the formula (10), p (k+1) represents an updated covariance vector estimated value at time k+1, I represents an identity matrix with the same dimension as K (k+1) H (k+1), K (k+1) represents a kalman filter gain matrix at time k+1, H (k+1) represents an observation matrix at time k+1, and p (k+1, K) represents a covariance vector predicted value at time k+1;
in the formula (11), ω (k+1) represents process noise at time k+1,the optimal estimated value of the frequency difference state vector at the moment k+1 is represented, and z (k+1) represents the measured value at the moment k+1;
Q k+1 ＝var(ω(k+1))(12)
in the formula (12), Q k+1 Represents the process noise variance at time k+1, ω (k+1) represents the process noise at time k+1, and var represents the operation of solving the process noise variance.
Optionally, when the differential evolution algorithm is introduced in the step 3) to dynamically adjust the process noise in the conventional extended kalman filtering process, the step of generating the process noise variance at any k moment through the differential evolution algorithm includes:
3.1 Population initialization is performed, and the functional expression is shown as the following formula (13):
x i,k (0)＝l k +rand(0,1)*(u k -l k ) (13)
in the formula (13), x i,k (0) A kth gene representing the ith "chromosome" of the 0 th generation in the population; the rand (0, 1) is a random number uniformly distributed between 0 and 1; u (u) k And l k Is the upper and lower bounds of the search, and u k And l k Taking the maximum value Q of the process noise variance respectively max And minimum value Q min ；
3.2 Performing mutation operation on the population, wherein the function expression of the mutation operation is shown as the following formula (14):
X i (g+1)＝x r1 (g)+F*[x r2 (g)-x r3 (g)] (14)
in the formula (14), X i (g+1) is the resulting g+1 generation variant; f is a compression proportion factor, and the value range is 0-1; x is x r1 (g)、x r2 (g) And x r3 (g) 3 parents;
3.3 Cross-operating the population to retain better variables, and the cross-operating adopts a two-term cross-mode, wherein the function expression is shown as the following formula (15):
in the formula (15), y i,j (g+1) represents the result of the cross operation corresponding to the current individual, X i,j (g+1) represents the component corresponding to the variant, x i,j (g) Representing the components corresponding to the current individuals, wherein r is a random number which is uniformly distributed between 0 and 1 and is generated by each variable; cr is the crossover probability of the variable; rnd is an integer uniformly distributed between 1 and d, and if r is less than cr, the component X corresponding to the variant is accepted i,j (g+1), otherwise, the component x corresponding to the current individual is reserved i,j (g)；
3.4 Aiming at the result obtained by the cross operation, adopting a greedy selection mode shown in the following formula (16) to select an optimal solution and taking the optimal solution as a process noise variance at the moment k;
in formula (16), x i (g+1) is the result of the selection, y i (g+1) represents an individual corresponding to the result of the cross operation, x i (g) Representing the current individual, fy i (g+1)]Representing individual y i (g+1) corresponding population fitness value, fx i (g)]Representing individual x i (g) Corresponding population fitness values.
Optionally, step 4) integrating the weight of the innovation sequence into the extended Kalman filtering process after the dynamic adjustment of the process noise is completed, and the obtained function expression of the update equation of the extended Kalman filtering process is shown in the following formulas (17) - (21);
D k+1 ＝P(k+1,k)+R k+1 (21)
in the formulas (17) to (21),representing the optimal estimate of the adjusted k+1 moment frequency difference state vector,/for>Represents the predicted value of the frequency difference state vector at the time of k+1, and K (k+1) represents the Kalman filtering gain matrix at the time of k+1, phi k+1 Representing a correction function τ k+1 A variable e for measuring the measurement error at time k+1 k+1 Indicating the measurement error at time k+1, z (k+1) indicating the measurement value of the frequency difference vector at time k+1,/>Representing the predicted value of the measurement vector at time k+1, τ k+1 Obeying χ 2 Distribution, τ when outliers occur k+1 Increase phi k+1 (. Cndot.) decrease; c (C) k+1 Lambda is the selected threshold or constant sequence k+1 Is thatMaximum eigenvalue of matrix, K k+1 Kalman filter gain, D, representing time k+1 k+1 The weight matrix representing the measurement error at time k+1, P (k+1, k) representing the covariance predicted value at time k+1, and R (k+1) representing the measurement noise variance at time k+1.
In addition, the invention also provides a system for taming the crystal oscillator frequency by using the GPS signal, which comprises:
the frequency difference model generation program unit is used for establishing a frequency difference mathematical model of the GPS second signal and the crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal;
the conventional filtering program unit is used for performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model, so that random errors and errors caused by nonlinear drift of crystal oscillator frequency are eliminated, and a filtering model greatly influenced by a jump wild value and system noise is obtained;
the noise adjusting program unit is used for introducing a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process so as to reduce the influence of the process noise on the filtering result of the frequency difference signal;
and the filtering integration program unit is used for integrating the weight of the innovation sequence into the extended Kalman filtering process after the dynamic adjustment of the process noise is completed.
In addition, the invention also provides a system for taming the crystal oscillator frequency by using the GPS signal, which comprises a computer device, wherein the computer device is programmed or configured to execute the steps of the method for taming the crystal oscillator frequency by using the GPS signal.
In addition, the invention also provides a system for taming the crystal oscillator frequency by using the GPS signal, which comprises a computer device, wherein a memory of the computer device is stored with a computer program which is programmed or configured to execute the method for taming the crystal oscillator frequency by using the GPS signal.
In addition, the invention also provides an intelligent terminal with GPS, which comprises a GPS module, a microprocessor and a memory, wherein the microprocessor is programmed or configured to execute the steps of the method for using the GPS signal to tame the crystal oscillator frequency, or the memory is stored with a computer program programmed or configured to execute the method for using the GPS signal to tame the crystal oscillator frequency.
Furthermore, the invention provides a computer readable storage medium having stored thereon a computer program programmed or configured to perform the method of using GPS signals to tame crystal oscillator frequencies.
Compared with the prior art, the invention has the following advantages: the method for taming the crystal oscillator frequency by using the GPS signals mainly comprises four parts: firstly, establishing a frequency difference mathematical model of the 1PPS signal and the crystal oscillator signal according to the error characteristics of the GPS. Providing a foundation for establishing a state equation and a measurement equation of the frequency difference signal in the next step through mathematical modeling of the frequency difference signal; performing conventional extended Kalman filtering processing on the frequency difference signal according to the mathematical model of the frequency difference signal, so as to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency and obtain a filtering model greatly influenced by jump wild values and system noise; thirdly, introducing a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering, so as to reduce the influence of the process noise on a frequency difference signal filtering result; and fourthly, integrating an innovation sequence weighting technology into the extended Kalman filtering algorithm, and reducing the influence of a jump wild value contained in the GPS output PPS signal on a filtering result by introducing a correction function into an optimal estimation equation of a frequency difference vector. Therefore, the innovation sequence weighting extended Kalman filtering algorithm based on the differential evolution algorithm can carry out filtering processing on jump wild values and random errors in GPS signals, accumulated errors in crystal oscillator signals and nonlinear frequency drift, and can improve the final frequency correction precision of the GPS locking secondary frequency standard.
Drawings
Fig. 1 shows the principle of locking a secondary frequency scale (high stable crystal oscillator) by using a GPS in the prior art.
Fig. 2 is a schematic diagram of a basic flow of a method according to an embodiment of the present invention.
Fig. 3 is a schematic diagram of a basic flow of introducing a differential evolution algorithm in an embodiment of the present invention.
Detailed Description
The invention will now be described in more detail by means of examples, which are given for illustrative purposes only and are not intended to limit the scope of the invention.
As shown in fig. 2, the method for taming the crystal oscillator frequency by using the GPS signal in this embodiment includes the following steps:
1) Establishing a frequency difference mathematical model of the GPS second signal (1 PPS signal) and the crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal;
2) Performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model, so as to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency and obtain a filtering model greatly influenced by a jump wild value and system noise;
3) 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;
4) And integrating the new sequence weighting 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 of the frequency difference signal.
The method for taming the crystal oscillator frequency by using the GPS signal aims at filtering noise of a frequency difference signal in the GPS locked crystal oscillator frequency technology, and the method for weighting and expanding Kalman filtering by the innovation sequence based on the differential evolution algorithm can process jump wild values and random errors in the GPS signal and random errors and nonlinear frequency drift in the crystal oscillator signal, so that final frequency correction precision is improved. The method mainly comprises four parts: firstly, establishing a frequency difference mathematical model of the 1PPS signal and the crystal oscillator signal according to the error characteristics of the GPS. Providing a foundation for establishing a state equation and a measurement equation of the frequency difference signal in the next step through mathematical modeling of the frequency difference signal; performing conventional extended Kalman filtering processing on the frequency difference signal according to the mathematical model of the frequency difference signal, so as to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency and obtain a filtering model greatly influenced by jump wild values and system noise; thirdly, introducing a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering, so as to reduce the influence of the process noise on a frequency difference signal filtering result; and fourthly, integrating an innovation sequence weighting technology into the extended Kalman filtering algorithm. And reducing the influence of the jump wild value contained in the GPS output PPS signal on the filtering result by introducing a correction function into the optimal estimation equation of the frequency difference vector.
In this embodiment, when setting the frequency difference signal and the filter correlation parameter, the process noise vector ω (k) is set with a covariance matrix of q=0.0005. Setting a measurement noise vector gamma (k), and setting a covariance matrix of the measurement noise vector gamma (k) as R=0.01; in addition, the initial value of the frequency difference is set to x 0 ＝5×10 -9 s。
When acquiring the GPS second signal sequence Y, the universal time UTC of the GPS second signal in the GPS second signal sequence Y (with the size of n) can be expressed as:
1-ε 1 ，2-ε 2 ，…，k-ε k ，…，n-ε n
when the crystal oscillator signal sequence Y 'is acquired, the universal time UTC of the crystal oscillator signal in the crystal oscillator signal sequence Y' (the size is n) can be expressed as:
1+a+b+c，1+a+2b+4c，…，k+a+bk+ck 2 ，…，n+a+bn+cn 2
wherein a is an initial phase error, b is an error coefficient considering frequency deviation, c is an error coefficient considering frequency linear drift; thus, in this embodiment, the detailed steps of step 1) include:
1.1 Acquiring a GPS second signal sequence Y, wherein a function expression of universal time UTC of GPS second signals at any k moment in the GPS second signal sequence Y is shown in a formula (1);
Y k ＝k-ε k (1)
in the above, Y k Universal time UTC, epsilon for GPS second signal representing time k in GPS second signal sequence Y k Is the error of the GPS second signal at the k moment and satisfies epsilon-N (0, sigma) 2 )；
Acquiring a crystal oscillator signal sequence Y ', wherein a function expression of universal time UTC of crystal oscillator signals at any k moment in the crystal oscillator signal sequence Y' is shown in a formula (2);
Y′ k ＝k+a+bk+ck 2 (2)
In the above, Y' k The universal time UTC of the crystal oscillator signal at the moment k in the crystal oscillator signal sequence Y', a is an initial phase error, b is an error coefficient considering frequency deviation, and c is an error coefficient considering frequency linear drift;
1.2 According to the GPS second signal sequence Y and the crystal oscillator signal sequence Y', obtaining a frequency difference sequence, wherein the function 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 2 +ε k (3)
in the above, X k Frequency difference signal representing k moment in frequency difference sequence, the frequency difference sequence is X 1 ，X 2 ，…，X k ，…，X n 。
In this embodiment, step 2) is configured to perform conventional extended kalman filtering processing on the frequency difference signal according to the mathematical model of the frequency difference signal, so as to eliminate random error and error caused by nonlinear drift of the crystal oscillator frequency, and obtain a filtering model that is greatly affected by the jump wild value and system noise. The step of performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model in the step 2) comprises the following steps:
2.1 A function expression of a state equation and an observation equation of the frequency difference is established as shown in a formula (4);
in the above formula, x (k+1) represents a frequency difference state vector at k+1, x (k) represents a frequency difference state vector at k, f and h represent a state function and a measurement function, respectively, and z (k+1) represents a measurement value at k+1; omega (k) represents a process noise vector at the moment k, and the covariance matrix of the process noise vector is Q; gamma (k) represents a measurement noise vector at the moment k, and the covariance matrix of the measurement noise vector is R;
2.2 Let the function expression of the state transition matrix F (k+1, k) from the k moment to the k+1 moment and the observation matrix H (k) from the k moment be shown as the formula (5);
in the above formula, f represents a state function,estimation of state vector representing frequency difference at k timeH (k) represents an observation matrix at the time of k, H represents a measurement function, and x (k) represents a frequency difference state vector at the time of k;
2.3 Establishing a predictive equation of a conventional extended Kalman filtering model as shown in formulas (6) to (8), and an update equation as shown in formulas (9) to (10), wherein the process noise is shown in formula (11), and the process noise variance is shown in formula (12);
in the formula (6), the amino acid sequence of the compound,a k+1 time frequency difference state vector predicted value based on the k time predicted value, f represents a state function,/->Representing the optimal estimated value of the frequency difference state vector at the k moment;
p(k+1,k)＝F(k+1,k)p(k,k)F T (k+1,k)+Q(k) (7)
in the formula (7), p (k+1, k) represents a covariance vector predictor at time k+1, F (k+1, k) represents a state transition matrix from time k to time k+1, p (k, k) represents a covariance vector predictor at time k, F T (k+1, k) represents the transpose of the state transition matrix from time k to time k+1, and 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 the formula (8), K (k+1) represents a Kalman filtering gain matrix at the time of k+1, p (k+1, K) represents a covariance vector predicted value at the time of k+1, and H T (k+1) represents a transpose of the observation matrix H (k+1) at time k+1, p (k+1, k) represents a covariance vector predictor at time k+1, and R (k+1) represents an observation noise variance vector at time k+1;
in the formula (9), the amino acid sequence of the compound,optimal estimate representing the state vector of the frequency difference at time k+1,/and>a predicted value of a frequency difference state vector at time k+1, K (k+1) represents a Kalman filtering gain matrix at time k+1, z (k+1) represents a measured value at time k+1, +.>Representing the state vector predictor +.>The obtained k+1 moment measurement vector predicted value;
p(k+1)＝[I-K(k+1)H(k+1)]p(k+1,k) (10)
in the formula (10), p (k+1) represents an updated covariance vector estimated value at time k+1, I represents an identity matrix with the same dimension as K (k+1) H (k+1), K (k+1) represents a kalman filter gain matrix at time k+1, H (k+1) represents an observation matrix at time k+1, and p (k+1, K) represents a covariance vector predicted value at time k+1;
in the formula (11), ω (k+1) represents process noise at time k+1,the optimal estimated value of the frequency difference state vector at the moment k+1 is represented, and z (k+1) represents the measured value at the moment k+1;
Q k+1 ＝var(ω(k+1)) (12)
in the formula (12), Q k+1 Represents the process noise variance at time k+1, ω (k+1) represents the process noise at time k+1, and var represents the operation of solving the process noise variance.
In this embodiment, step 3) introduces a differential evolution algorithm to dynamically adjust the process noise in the conventional extended kalman filtering algorithm, so as to reduce the influence of the process noise on the filtering result of the frequency difference signal. Conventional extended kalman generally assumes a process noise mean of zero and a variance of Q during filtering. Because of uncertain factors and abnormal disturbance of the system model, the process noise is also changed; fixed process noise will increase the error and reduce the filtering accuracy. Therefore, a differential evolution algorithm is introduced to solve the optimal variance of the process noise, and then the optimal variance is applied to an extended Kalman filtering algorithm to improve the filtering precision, and finally the frequency correction precision is improved. As shown in fig. 3, when the differential evolution algorithm is introduced in step 3) to dynamically adjust the process noise in the conventional extended kalman filter process, the step of generating the process noise variance at any k time by the differential evolution algorithm includes:
3.1 Population initialization is performed, and the functional expression is shown as the following formula (13):
x i,k (0)＝l k +rand(0,1)*(u k -l k ) (13)
in the formula (13), x i,k (0) A kth gene representing the ith "chromosome" of the 0 th generation in the population; the rand (0, 1) is a random number uniformly distributed between 0 and 1; u (u) k And l k Is the upper and lower bounds of the search, and u k And l k Taking the maximum value Q of the process noise variance respectively max And minimum value Q min The method comprises the steps of carrying out a first treatment on the surface of the When setting the initial value of the differential evolution algorithm in this embodiment, the mutation rate f=0.5, the crossover probability cr=0.9, the population size np=10, the iteration number g=100, and the search boundary u k ＝0.005,l k ＝-0.005；
3.2 Performing mutation operation on the population, wherein the function expression of the mutation operation is shown as the following formula (14):
X i (g+1)＝x r1 (g)+F*[x r2 (g)-x r3 (g)] (14)
in the formula (14), X i (g+1) is the resulting g+1 generation variant; f is a compression proportion factor, and the value range is 0-1; x is x r1 (g)、x r2 (g) And x r3 (g) For 3 parents；
3.3 Cross-operating the population to retain better variables, and the cross-operating adopts a two-term cross-mode, wherein the function expression is shown as the following formula (15):
in the formula (15), y i,j (g+1) represents the result of the cross operation corresponding to the current individual, X i,j (g+1) represents the component corresponding to the variant, x i,j (g) Representing the components corresponding to the current individuals, wherein r is a random number which is uniformly distributed between 0 and 1 and is generated by each variable; cr is the crossover probability of the variable; rnd is an integer uniformly distributed between 1 and d, and if r is less than cr, the component X corresponding to the variant is accepted i,j (g+1), otherwise, the component x corresponding to the current individual is reserved i,j (g)；
3.4 Aiming at the result obtained by the cross operation, adopting a greedy selection mode shown in the following formula (16) to select an optimal solution and taking the optimal solution as a process noise variance at the moment k;
in formula (16), x i (g+1) is the result of the selection, y i (g+1) represents an individual corresponding to the result of the cross operation, x i (g) Representing the current individual, fy i (g+1)]Representing individual y i (g+1) corresponding population fitness value, fx i (g)]Representing individual x i (g) Corresponding population fitness values. Assigning an optimal solution to Q i Substituting the covariance prediction equation can improve the filtering precision of the extended Kalman.
In this embodiment, step 4) integrates the innovation sequence weighting technique into the extended kalman filtering algorithm to reduce the influence of the hop wild value contained in the GPS output PPS signal on the filtering result. In practical application, because the GPS signal is interfered by noise and generates a larger jump wild value, the frequency difference signal with the wild value can offset or even diverge the filtering result. In order to solve the problem of outlier resistance of the extended Kalman algorithm, a scholars introduces a Kalman filtering method for weighting an innovation sequence to eliminate the influence of outliers. When the observed data does not contain the wild value, the filtering algorithm can normally operate, and the observed noise is filtered; when the observed data contains the wild value, the influence of the wild value on the filtering result can be overcome or reduced, so that the filtering result is more accurate. In this embodiment, step 4) integrates the weighting of the innovation sequence to the extended kalman filter process after the dynamic adjustment of the process noise is completed, and the obtained function expression of the update equation of the extended kalman filter process is shown in the following formulas (17) to (21);
D k+1 ＝P(k+1,k)+R k+1 (21)
In the formulas (17) to (21),representing the optimal estimate of the adjusted k+1 moment frequency difference state vector,/for>Represents the predicted value of the frequency difference state vector at the time of k+1, and K (k+1) represents the Kalman filtering gain matrix at the time of k+1, phi k+1 Representing a correction function τ k+1 A variable representing the measurement error magnitude at time k +1,e k+1 indicating the measurement error at time k+1, z (k+1) indicating the measurement value of the frequency difference vector at time k+1,/>Representing the predicted value of the measurement vector at time k+1, τ k+1 Obeying χ 2 Distribution, τ when outliers occur k+1 Increase phi k+1 (. Cndot.) decrease; c (C) k+1 Lambda is the selected threshold or constant sequence k+1 Is thatMaximum eigenvalue of matrix, K k+1 Kalman filter gain, D, representing time k+1 k+1 The weight matrix representing the measurement error at time k+1, P (k+1, k) representing the covariance predicted value at time k+1, and R (k+1) representing the measurement noise variance at time k+1.
In summary, the main technology conventionally applied to the low-cost high-precision frequency standard is the GPS locking secondary frequency standard (crystal oscillator) technology, so as to reduce the influence of random errors and hopping wild values of GPS signals and nonlinear frequency offset of crystal oscillator signals on frequency difference signals in the GPS locking secondary frequency standard technology. The embodiment provides a method for processing a frequency difference signal in a GPS locked crystal oscillator frequency technology by using a differential evolution algorithm-based innovation sequence weighted extension Kalman filtering method, which mainly comprises the following four parts: firstly, establishing a frequency difference mathematical model of the 1PPS signal and the crystal oscillator signal according to the error characteristics of the GPS. Providing a foundation for establishing a state equation and a measurement equation of the frequency difference signal in the next step through mathematical modeling of the frequency difference signal; performing conventional extended Kalman filtering processing on the frequency difference signal according to the mathematical model of the frequency difference signal, so as to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency and obtain a filtering model greatly influenced by jump wild values and system noise; thirdly, introducing a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering, so as to reduce the influence of the process noise on a frequency difference signal filtering result; and fourthly, integrating an innovation sequence weighting technology into the extended Kalman filtering algorithm. And reducing the influence of the jump wild value contained in the GPS output PPS signal on the filtering result by introducing a correction function into the optimal estimation equation of the frequency difference vector. In addition, the method of the embodiment adopts an extended Kalman filtering algorithm to process the frequency difference signal in the process of locking the crystal oscillator frequency by the GPS, and considers nonlinear errors of the crystal oscillator due to aging and temperature influence while processing the random errors of the GPS signal. The method introduces an innovation sequence weighting and differential evolution algorithm into an extended Kalman filtering process, and can reduce the influence of larger jump wild values in GPS signals and continuous change of system noise on filtering results.
In addition, the embodiment also provides a system for taming the crystal oscillator frequency by using the GPS signal, which comprises:
the frequency difference model generation program unit is used for establishing a frequency difference mathematical model of the GPS second signal and the crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal;
the conventional filtering program unit is used for performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model, so that random errors and errors caused by nonlinear drift of crystal oscillator frequency are eliminated, and a filtering model greatly influenced by a jump wild value and system noise is obtained;
the noise adjusting program unit is used for introducing a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process so as to reduce the influence of the process noise on the filtering result of the frequency difference signal;
and the filtering integration program unit is used for integrating the weight of the innovation sequence into the extended Kalman filtering process after the dynamic adjustment of the process noise is completed.
In addition, the embodiment also provides a system for taming the crystal oscillator frequency by using the GPS signal, which comprises a computer device, wherein the computer device is programmed or configured to execute the steps of the method for taming the crystal oscillator frequency by using the GPS signal.
In addition, the embodiment also provides a system for taming the crystal oscillator frequency by using the GPS signal, which comprises a computer device, wherein a memory of the computer device is stored with a computer program programmed or configured to execute the method for taming the crystal oscillator frequency by using the GPS signal.
In addition, the embodiment also provides an intelligent terminal with GPS, which comprises a GPS module, a microprocessor and a memory, wherein the microprocessor is programmed or configured to execute the steps of the method for using the GPS signal to tame the crystal oscillator frequency, or the memory is stored with a computer program programmed or configured to execute the method for using the GPS signal to tame the crystal oscillator frequency. The intelligent terminal can be a smart phone, a tablet computer or other types of terminal equipment.
In addition, the embodiment also provides a computer readable storage medium, and the computer readable storage medium stores a computer program programmed or configured to execute the method for taming the crystal oscillator frequency by using the GPS signal.
It will be appreciated by those skilled in the art that embodiments of the present application may be provided as a method, system, or computer program product. 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, and the like) having computer-usable program code embodied therein. The present application is directed to methods, apparatus (systems), and computer program products in accordance with embodiments of the present application that produce means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks by reference to the instructions that execute in the flowchart and/or processor. These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks. These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
The above description is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above examples, and all technical solutions belonging to the concept of the present invention belong to the protection scope of the present invention. It should be noted that modifications and adaptations to the present invention may occur to one skilled in the art without departing from the principles of the present invention and are intended to be within the scope of the present invention.
Claims (8)
1. A method for taming crystal oscillator frequency using GPS signals, the method comprising the steps of:
1) Establishing a frequency difference mathematical model of the GPS second signal and the crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal;
2) Performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model, so as to eliminate random errors and errors caused by nonlinear drift of crystal oscillator frequency; the step of performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model comprises the following steps:
2.1 A function expression of a state equation and an observation equation of the frequency difference is established as shown in a formula (4);
in the above formula, x (k+1) represents a frequency difference state vector at k+1, x (k) represents a frequency difference state vector at k, f and h represent a state function and a measurement function, respectively, and z (k+1) represents a measurement value at k+1; omega (k) represents a process noise vector at the moment k, and the covariance matrix of the process noise vector is Q; gamma (k) represents a measurement noise vector at the moment k, and the covariance matrix of the measurement noise vector is R;
2.2 Let the function expression of the state transition matrix F (k+1, k) from the k moment to the k+1 moment and the observation matrix H (k) from the k moment be shown as the formula (5);
in the above formula, f represents a state function,the estimation of the frequency difference state vector at the moment k is represented, h represents a measurement function, and x (k) represents the frequency difference state vector at the moment k;
2.3 Establishing a predictive equation of a conventional extended Kalman filtering model as shown in formulas (6) to (8), and an update equation as shown in formulas (9) to (10), wherein the process noise is shown in formula (11), and the process noise variance is shown in formula (12);
in the formula (6), the amino acid sequence of the compound,a k+1 time frequency difference state vector predicted value based on the k time predicted value, f represents a state function,/->Representing the optimal estimated value of the frequency difference state vector at the k moment;
p(k+1,k)＝F(k+1,k)p(k,k)F T (k+1,k)+Q(k) (7)
in the formula (7), p (k+1, k) represents a covariance vector predictor at time k+1, F (k+1, k) represents a state transition matrix from time k to time k+1, p (k, k) represents a covariance vector predictor at time k, F T (k+1, k) represents the transpose of the state transition matrix from time k to time k+1, and 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 the formula (8), K (k+1) is as followsThe Kalman filter gain matrix at time k+1 is shown, p (k+1, k) represents the covariance vector predictor at time k+1, H T (k+1) represents a transpose of the observation matrix H (k+1) at time k+1, p (k+1, k) represents a covariance vector predictor at time k+1, and R (k+1) represents an observation noise variance vector at time k+1;
In the formula (9), the amino acid sequence of the compound,optimal estimate representing the state vector of the frequency difference at time k+1,/and>a predicted value of a frequency difference state vector at time k+1, K (k+1) represents a Kalman filtering gain matrix at time k+1, z (k+1) represents a measured value at time k+1, +.>Representing the state vector predictor +.>The obtained k+1 moment measurement vector predicted value;
p(k+1)＝[I-K(k+1)H(k+1)]p(k+1,k) (10)
in the formula (10), p (k+1) represents an updated covariance vector estimated value at time k+1, I represents an identity matrix with the same dimension as K (k+1) H (k+1), K (k+1) represents a kalman filter gain matrix at time k+1, H (k+1) represents an observation matrix at time k+1, and p (k+1, K) represents a covariance vector predicted value at time k+1;
in the formula (11), ω (k+1) represents time k+1The noise of the process is generated by the process,the optimal estimated value of the frequency difference state vector at the moment k+1 is represented, and z (k+1) represents the measured value at the moment k+1;
Q k+1 ＝var(ω(k+1)) (12)
in the formula (12), Q k+1 Representing the process noise variance at time k+1, ω (k+1) representing the process noise at time k+1, var representing the operation of solving the process noise variance;
3) 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;
4) 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 of the frequency difference signal; and integrating the weight of the innovation sequence into the extended Kalman filtering process after the dynamic adjustment of the process noise is completed, and the obtained function expression of the update equation of the extended Kalman filtering process is shown in the following formulas (17) to (21);
D k+1 ＝P(k+1,k)+R k+1 (21)
In the formulas (17) to (21),represents the optimal estimated value of the adjusted k+1 time frequency difference state vector,represents the predicted value of the frequency difference state vector at the time of k+1, and K (k+1) represents the Kalman filtering gain matrix at the time of k+1, phi k+1 Representing a correction function τ k+1 A variable e for measuring the measurement error at time k+1 k+1 Indicating the measurement error at time k+1, z (k+1) indicating the measurement value of the frequency difference vector at time k+1,/>Representing the predicted value of the measurement vector at time k+1, τ k+1 Obeying χ 2 Distribution, τ when outliers occur k+1 Increase phi k+1 (. Cndot.) decrease; c (C) k+1 Lambda is the selected threshold or constant sequence k+1 Is thatMaximum eigenvalue of matrix, K k+1 Kalman filter gain, D, representing time k+1 k+1 Weight matrix representing measurement error magnitude at time k+1, P (k+1, k) represents covariance predicted value at time k+1, R k+1 The measured noise variance at time k+1 is shown.
2. The method for taming a crystal oscillator frequency using GPS signals according to claim 1, wherein the detailed steps of step 1) include:
1.1 Acquiring a GPS second signal sequence Y, wherein a function expression of universal time UTC of GPS second signals at any k moment in the GPS second signal sequence Y is shown in a formula (1);
Y k ＝k-ε k (1)
in the above, Y k GPS seconds representing time k in GPS second signal sequence Y Universal time UTC, epsilon of signal k Is the error of the GPS second signal at the k moment and satisfies epsilon-N (0, sigma) 2 )；
Acquiring a crystal oscillator signal sequence Y ', wherein a function expression of universal time UTC of crystal oscillator signals at any k moment in the crystal oscillator signal sequence Y' is shown in a formula (2);
Y′ k ＝k+a+bk+ck 2 (2)
in the above, Y' k The universal time UTC of the crystal oscillator signal at the moment k in the crystal oscillator signal sequence Y', a is an initial phase error, b is an error coefficient considering frequency deviation, and c is an error coefficient considering frequency linear drift;
1.2 According to the GPS second signal sequence Y and the crystal oscillator signal sequence Y', obtaining a frequency difference sequence, wherein the function 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 2 +ε k (3)
in the above, X k And represents the frequency difference signal at time k in the frequency difference sequence.
3. The method for training the crystal oscillator frequency by using the GPS signal according to claim 1, wherein when the differential evolution algorithm is introduced in the step 3) to dynamically adjust the process noise in the conventional extended Kalman filtering process, the step of generating the process noise variance at any k time by the differential evolution algorithm comprises the following steps:
3.1 Population initialization is performed, and the functional expression is shown as the following formula (13):
x i,k (0)＝l k +rand(0,1)*(u k -l k ) (13)
In the formula (13), x i,k (0) A kth gene representing the ith "chromosome" of the 0 th generation in the population; the rand (0, 1) is a random number uniformly distributed between 0 and 1; u (u) k And l k Is the upper and lower bounds of the search, and u k And l k Taking the maximum value Q of the process noise variance respectively max And minimum value Q min ；
3.2 Performing mutation operation on the population, wherein the function expression of the mutation operation is shown as the following formula (14):
X i (g+1)＝x r1 (g)+F*[x r2 (g)-x r3 (g)] (14)
in the formula (14), X i (g+1) is the resulting g+1 generation variant; f is a compression proportion factor, and the value range is 0-1; x is x r1 (g)、x r2 (g) And x r3 (g) 3 parents;
3.3 Cross-operating the population to retain better variables, and the cross-operating adopts a two-term cross-mode, wherein the function expression is shown as the following formula (15):
in the formula (15), y i,j (g+1) represents the result of the cross operation corresponding to the current individual, X i,j (g+1) represents the component corresponding to the variant, x i,j (g) Representing the components corresponding to the current individuals, wherein r is a random number which is uniformly distributed between 0 and 1 and is 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 receives the component X corresponding to the variant i,j (g+1), otherwise, the component x corresponding to the current individual is reserved i,j (g)；
3.4 Aiming at the result obtained by the cross operation, adopting a greedy selection mode shown in the following formula (16) to select an optimal solution and taking the optimal solution as a process noise variance at the moment k;
In formula (16), x i (g+1) is the result of the selection, y i (g+1) represents an individual corresponding to the result of the cross operation, x i (g) Representing the current individual, fy i (g+1)]Representing individual y i (g+1) corresponding population fitness value, fx i (g)]Representing individual x i (g) Corresponding toPopulation fitness value.
4. A system for taming crystal oscillator frequency using GPS signals, comprising:
the frequency difference model generation program unit is used for establishing a frequency difference mathematical model of the GPS second signal and the crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal;
the conventional filtering program unit is used for performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model so as to eliminate random errors and errors caused by nonlinear drift of the crystal oscillator frequency; the step of performing conventional extended Kalman filtering processing on the frequency difference signal according to the frequency difference mathematical model comprises the following steps:
2.1 A function expression of a state equation and an observation equation of the frequency difference is established as shown in a formula (4);
in the above formula, x (k+1) represents a frequency difference state vector at k+1, x (k) represents a frequency difference state vector at k, f and h represent a state function and a measurement function, respectively, and z (k+1) represents a measurement value at k+1; omega (k) represents a process noise vector at the moment k, and the covariance matrix of the process noise vector is Q; gamma (k) represents a measurement noise vector at the moment k, and the covariance matrix of the measurement noise vector is R;
2.2 Let the function expression of the state transition matrix F (k+1, k) from the k moment to the k+1 moment and the observation matrix H (k) from the k moment be shown as the formula (5);
in the above formula, f represents a state function,the estimation of the frequency difference state vector at the moment k is represented, h represents a measurement function, and x (k) represents the frequency difference state vector at the moment k;
2.3 Establishing a predictive equation of a conventional extended Kalman filtering model as shown in formulas (6) to (8), and an update equation as shown in formulas (9) to (10), wherein the process noise is shown in formula (11), and the process noise variance is shown in formula (12);
in the formula (6), the amino acid sequence of the compound,a k+1 time frequency difference state vector predicted value based on the k time predicted value, f represents a state function,/->Representing the optimal estimated value of the frequency difference state vector at the k moment;
p(k+1,k)＝F(k+1,k)p(k,k)F T (k+1,k)+Q(k) (7)
in the formula (7), p (k+1, k) represents a covariance vector predictor at time k+1, F (k+1, k) represents a state transition matrix from time k to time k+1, p (k, k) represents a covariance vector predictor at time k, F T (k+1, k) represents the transpose of the state transition matrix from time k to time k+1, and 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 the formula (8), K (k+1) represents a Kalman filtering gain matrix at the time of k+1, p (k+1, K) represents a covariance vector predicted value at the time of k+1, and H T (k+1) represents a transpose of the observation matrix H (k+1) at time k+1, p (k+1, k) represents a covariance vector predictor at time k+1, and R (k+1) represents an observation noise variance vector at time k+1;
In the formula (9), the amino acid sequence of the compound,optimal estimate representing the state vector of the frequency difference at time k+1,/and>a predicted value of a frequency difference state vector at time k+1, K (k+1) represents a Kalman filtering gain matrix at time k+1, z (k+1) represents a measured value at time k+1, +.>Representing the state vector predictor +.>The obtained k+1 moment measurement vector predicted value;
p(k+1)＝[I-K(k+1)H(k+1)]p(k+1,k) (10)
in the formula (10), p (k+1) represents an updated covariance vector estimated value at time k+1, I represents an identity matrix with the same dimension as K (k+1) H (k+1), K (k+1) represents a kalman filter gain matrix at time k+1, H (k+1) represents an observation matrix at time k+1, and p (k+1, K) represents a covariance vector predicted value at time k+1;
in the formula (11), ω (k+1) represents process noise at time k+1,the optimal estimated value of the frequency difference state vector at the moment k+1 is represented, and z (k+1) represents the measured value at the moment k+1;
Q k+1 ＝var(ω(k+1)) (12)
in the formula (12), Q k+1 Representing the process noise variance at time k+1, ω (k+1) representing the process noise at time k+1, var representing the operation of solving the process noise variance;
the noise adjusting program unit is used for introducing a differential evolution algorithm to dynamically adjust the process noise in the conventional extended Kalman filtering process so as to reduce the influence of the process noise on the filtering result of the frequency difference signal;
The filtering integration program unit is used for integrating 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 the jump wild value on the filtering result of the frequency difference signal; and integrating the weight of the innovation sequence into the extended Kalman filtering process after the dynamic adjustment of the process noise is completed, and the obtained function expression of the update equation of the extended Kalman filtering process is shown in the following formulas (17) to (21);
D k+1 ＝P(k+1,k)+R k+1 (21)
in the formulas (17) to (21),represents the optimal estimated value of the adjusted k+1 time frequency difference state vector,represents the predicted value of the frequency difference state vector at the time of k+1, and K (k+1) represents the Kalman filtering gain at the time of k+1Matrix phi k+1 Representing a correction function τ k+1 A variable e for measuring the measurement error at time k+1 k+1 Indicating the measurement error at time k+1, z (k+1) indicating the measurement value of the frequency difference vector at time k+1,/>Representing the predicted value of the measurement vector at time k+1, τ k+1 Obeying χ 2 Distribution, τ when outliers occur k+1 Increase phi k+1 (. Cndot.) decrease; c (C) k+1 Lambda is the selected threshold or constant sequence k+1 Is thatMaximum eigenvalue of matrix, K k+1 Kalman filter gain, D, representing time k+1 k+1 Weight matrix representing measurement error magnitude at time k+1, P (k+1, k) represents covariance predicted value at time k+1, R k+1 The measured noise variance at time k+1 is shown.
5. A system for frequency-tuning a crystal oscillator using GPS signals, comprising computer equipment programmed or configured to perform the steps of the method of frequency-tuning a crystal oscillator using GPS signals as claimed in any of claims 1 to 3.
6. A system for frequency-domesticating a crystal oscillator using GPS signals, comprising a computer device, wherein a memory of the computer device has stored thereon a computer program programmed or configured to perform the method of frequency-domesticating a crystal oscillator using GPS signals as claimed in any one of claims 1 to 3.
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 the steps of the method for frequency-tuning a crystal oscillator using GPS signals according to any one of claims 1 to 3, or the memory has stored thereon a computer program programmed or configured to perform the method for frequency-tuning a crystal oscillator using GPS signals according to any one of claims 1 to 3.
8. A computer readable storage medium having stored thereon a computer program programmed or configured to perform the method of using GPS signals to tame crystal oscillator frequencies of any of claims 1 to 3.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010523171.8A CN111711446B (en) | 2020-06-10 | 2020-06-10 | Method, system and medium for taming crystal oscillator frequency by using GPS signal |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010523171.8A CN111711446B (en) | 2020-06-10 | 2020-06-10 | Method, system and medium for taming crystal oscillator frequency by using GPS signal |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111711446A CN111711446A (en) | 2020-09-25 |
| CN111711446B true CN111711446B (en) | 2023-08-15 |
Family
ID=72539968
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010523171.8A Active CN111711446B (en) | 2020-06-10 | 2020-06-10 | Method, system and medium for taming crystal oscillator frequency by using GPS signal |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111711446B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114362140B (en) * | 2021-12-01 | 2023-11-21 | 国网内蒙古东部电力有限公司通辽供电公司 | High-precision time keeping method and device suitable for power distribution network measuring device |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106026919A (en) * | 2016-05-16 | 2016-10-12 | 南京理工大学 | Time-keeping compensation method for high-precision crystal oscillator |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7317361B2 (en) * | 2003-07-23 | 2008-01-08 | The Johns Hopkins University | Ensemble oscillator and related methods |
-
2020
- 2020-06-10 CN CN202010523171.8A patent/CN111711446B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106026919A (en) * | 2016-05-16 | 2016-10-12 | 南京理工大学 | Time-keeping compensation method for high-precision crystal oscillator |
Non-Patent Citations (1)
| Title |
|---|
| 宁玉磊.GPS信号校准晶振信号频率源误差在线修正方法.《中国测试》.2016,第15-19页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111711446A (en) | 2020-09-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111650617B (en) | Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering | |
| Hou et al. | Gray-box parsimonious subspace identification of Hammerstein-type systems | |
| Weiss et al. | AT2, a new time scale algorithm: AT1 plus frequency variance | |
| CN109508510B (en) | Improved Kalman filtering-based rubidium atomic clock parameter estimation algorithm | |
| EP0461557B1 (en) | Time scale computation system including complete and weighted ensemble definition | |
| CN107465393B (en) | System and method for frequency compensation of real time clock system | |
| CN109299496B (en) | High-precision synchronous clock generation method | |
| CN110554597B (en) | Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering | |
| CN114966766B (en) | Method, device and system for constructing navigation constellation time reference | |
| CN109902882A (en) | Atomic clock clock deviation prediction model training method and device | |
| CN114966767B (en) | Method, device and system for constructing autonomous time-frequency reference of navigation satellite | |
| CN111711446B (en) | Method, system and medium for taming crystal oscillator frequency by using GPS signal | |
| CN112433235B (en) | Method, system and medium for determining time reference | |
| CN114647178A (en) | Automatic atomic clock calibration method and system based on Beidou and ground reference transmission | |
| CN116980065A (en) | Clock calibration method, clock calibration device, terminal equipment and storage medium | |
| CN110309593B (en) | Device and method for predicting aging rate of constant-temperature crystal oscillator | |
| CN109474276B (en) | CPT atomic clock frequency synchronization control method and system | |
| CN112505386B (en) | Method and system for detecting current value of direct current charging pile | |
| Cortés et al. | Adaptive techniques in scalar tracking loops with direct-state Kalman-filter | |
| CN109655081B (en) | On-orbit adaptive correction method and system for star sensor optical system parameters | |
| CN114567288B (en) | Distribution collaborative nonlinear system state estimation method based on variable decibels | |
| Allan et al. | An accuracy algorithm for an atomic time scale | |
| CN115406426A (en) | MEMS gyroscope interface circuit and modulation method | |
| CN115904000A (en) | Real-time clock crystal oscillator compensation method and system based on orthogonal least square method curve fitting | |
| CN111950123B (en) | Gyroscope error coefficient curve fitting prediction method and system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |