CN111650617B - Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering - Google Patents
Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering Download PDFInfo
- Publication number
- CN111650617B CN111650617B CN202010524159.9A CN202010524159A CN111650617B CN 111650617 B CN111650617 B CN 111650617B CN 202010524159 A CN202010524159 A CN 202010524159A CN 111650617 B CN111650617 B CN 111650617B
- Authority
- CN
- China
- Prior art keywords
- vector
- time
- frequency
- representing
- 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
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
- H03H17/00—Networks using digital techniques
- H03H17/02—Frequency selective networks
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
- H03H17/00—Networks using digital techniques
- H03H17/02—Frequency selective networks
- H03H17/0202—Two or more dimensional filters; Filters for complex signals
- H03H2017/0205—Kalman filters
Landscapes
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Mathematical Physics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention discloses a crystal oscillator frequency taming method, a crystal oscillator frequency taming system and a crystal oscillator frequency taming medium based on innovation weighted self-adaptive Kalman filtering, wherein the method comprises the steps of establishing a frequency difference mathematical model of a GPS (global positioning system) second signal and a crystal oscillator signal according to the error characteristics of the GPS second signal and the crystal oscillator signal; sigma sampling is carried out on the frequency difference mathematical model to obtain a sigma point set at the moment k; determining a state equation and an updating equation of the adaptive insensitive Kalman filtering algorithm aiming at a sigma point set at the moment k; and integrating innovation sequence weighting into a state equation and an update equation of the adaptive insensitive Kalman filtering algorithm to reduce the influence of a larger field of the GPS on the frequency correction precision. Based on the innovation weighting self-adaptive insensitive Kalman filtering technology, the invention carries out filtering processing on the jump wild value and the random error of the GPS signal, the accumulated error and the frequency drift of the crystal oscillator signal, and finally improves the final frequency correction precision of the GPS locking secondary frequency scale.
Description
Technical Field
The invention relates to a method for domesticating crystal oscillator signals by a GPS (global positioning system), in particular to a crystal oscillator frequency domesticating method, a crystal oscillator frequency domesticating system and a crystal oscillator frequency domesticating medium based on innovation weighted adaptive insensitive Kalman filtering, which are used for realizing high-precision domestication of the crystal oscillator frequency by the GPS signals.
Background
With the rapid development of network technology, the system has higher and higher requirements on the time reference. The time reference in the power system has important significance for accident analysis, improvement of power quality, power grid dispatching operation and power time-sharing charging.
The current time (frequency) standard can be divided into a first-level frequency standard (cesium atomic frequency standard), a second-level frequency standard (including rubidium atomic frequency standard and a high-stability crystal oscillator) and other frequency standards (other crystal oscillators unexpected by the high-stability crystal oscillator) according to the application field and performance indexes. Although the primary frequency scale has higher long-term stability and accuracy, the price is very expensive and is difficult to be applied 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 of the secondary frequency scale are much inferior to those of the former frequency scale, so that the secondary frequency scale is difficult to apply to high-precision occasions.
The principle of the prior art using a GPS locked secondary frequency scale (high stable crystal oscillator) is shown in fig. 1: firstly, the relative frequency difference between a pulse per second signal PPS output by a GPS receiver and a pulse signal output by a crystal oscillator after frequency division is measured, then the frequency difference signal is filtered, finally, a control signal is generated, and the crystal oscillator frequency is continuously adjusted so as to output frequency output with high accuracy and high stability. The technique solves the above problems to a certain extent, namely, the frequency standard of high long-term stability and average accuracy is realized with lower cost. However, in practical application, random errors and jump field values exist in the 1PPS signal of the GPS, and frequency drift caused by aging, temperature and other characteristics of the crystal oscillator causes certain difficulty in frequency calibration.
In order to solve the above problems, some scholars propose to process the measured frequency difference data by using a method of sampling and averaging for multiple times, and the calculation is simple but the frequency correction time is long. The learners calculate and compensate the error of the crystal oscillator signal and the GPS signal by establishing a frequency difference mathematical model and utilizing a least square method, but the method is more complex in calculation and does not consider the influence of the field value of the GPS signal on the frequency correction precision. Still, researchers propose to eliminate the noise in the frequency difference signal by using the kalman filtering algorithm, but the conventional way of eliminating the noise in the frequency difference signal by using the kalman filtering algorithm does not consider the influence of the crystal oscillator frequency drift, so that the frequency calibration precision is limited.
Disclosure of Invention
The technical problems to be solved by the invention are as follows: aiming at the problems in the prior art, the invention provides a crystal oscillator frequency domestication method, a crystal oscillator frequency domestication system and a crystal oscillator frequency domestication medium based on innovation weighted adaptive insensitive Kalman filtering.
In order to solve the technical problems, the invention adopts the technical scheme that:
a crystal oscillator frequency taming method based on innovation weighted adaptive insensitive Kalman filtering comprises the following steps:
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) Sigma sampling is carried out on the frequency difference mathematical model to obtain a sigma point set at the moment k;
3) Determining a state equation and an updating equation of a self-adaptive insensitive Kalman filtering algorithm according to a sigma point set at the moment k;
4) And integrating innovation sequence weighting into a state equation and an update equation of the self-adaptive insensitive Kalman filtering algorithm to reduce the influence of a larger field value of the GPS on frequency correction precision.
Optionally, the detailed steps of step 1) include:
1.1 A function expression of universal time UTC of the GPS second signal at any k time in the GPS second signal sequence Y is obtained as shown in formula (1);
Y k ＝k-ε k (1)
in the above formula, Y k Universal time UTC, epsilon of GPS-seconds signal representing time k in GPS-seconds signal sequence Y k Is the error of GPS second signal at the time of k, and satisfies epsilon-N (0, sigma) 2 )；
Acquiring a crystal oscillator signal sequence Y ', wherein the functional expression of the universal time UTC of the crystal oscillator signal at any k time in the crystal oscillator signal sequence Y' is shown as a formula (2);
Y′ k ＝k+a+bk+ck 2 (2)
in the above formula, Y' k The universal time UTC of the crystal oscillator signal at the time k in the crystal oscillator signal sequence Y' is represented, wherein 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', a frequency difference sequence is obtained, and a function expression of a frequency difference signal at any k time in the frequency difference sequence is shown as the following formula;
X k ＝Y′ k -Y k ＝a+bk+ck 2 +ε k (3)
in the above formula, X k Representing the frequency difference signal at time k in the sequence of frequency differences.
Optionally, the detailed steps of step 2) include:
2.1 Based on a mathematical model of frequency difference, a prediction equation for frequency difference is established:
X' k ＝X' k-1 +2ck+b-c+ω(k) (4)
in the above formula, X' k Denotes a predicted value, X ', of the frequency difference signal at time k' k-1 Representing a predicted value of the frequency difference signal at the time k-1, b being an error coefficient considering frequency deviation, c being an error coefficient considering frequency linear drift, and ω (k) representing process noise at the time k;
2.2 Let u (k) =2ck + b-c represent input vector, establish state equation and measurement equation of frequency difference:
in the above formula, x (k + 1) represents a frequency difference state vector at the time k +1, x (k) represents a frequency difference state vector at the time k, u (k) represents an input vector, and z (k + 1) represents a measurement value at the time k + 1; omega (k) represents a process noise vector at the time k, and a covariance matrix of the process noise vector is Q; gamma (k) represents the measurement noise vector at time k, and the covariance matrix is R;
2.3 X) performing dimension expansion on the frequency difference state vector a ＝[x T ω T γ T ]Wherein x is a Representing the expanded state vector, x T Representing transposes of the state vector of the frequency difference, ω T Representing the transpose of the process noise vector, gamma T Representing the transpose of the measurement noise vector; determining an initial state of a frequency difference state vector:
in the above formula, the first and second carbon atoms are,
representing the initial value of the frequency difference state vector, E (x) 0 ) Means, x, representing the initial value of the frequency difference vector 0 Representing the initial value of the frequency difference vector, P 0 An initial value of a covariance vector representing a frequency difference vector,
representing the covariance of the initial value of the frequency difference vector;
determining the initial state of the state vector after the dimension expansion of the frequency difference state vector:
in the above formula, the first and second carbon atoms are,
the initial value of the state vector after the dimension expansion is shown,
represents the mean value of the state vector after the dimension expansion,
representing dimension extensionsThe latter state vector is then used to determine the state,
represents the initial value of the state vector of the frequency difference,
representing the initial value of the covariance of the state vector after dimension expansion,
representing the covariance of the initial value of the state vector after dimension expansion, P 0 Representing an initial covariance vector value of the frequency difference vector, Q representing a covariance matrix of a process noise vector omega (k) at time k, and R representing a covariance matrix of a measurement noise vector gamma (k) at time k;
2.4 Symmetric sigma sampling is adopted for the frequency difference state vector to obtain a sigma point set of a k moment state
i =1, \ 8230, L, wherein
Indicating sigma points with the ith dimension a, L indicating the sampling number of the sigma points, L =2a +1, and the state dimension a = n + q + m, wherein n indicates the dimension of the frequency difference state vector, q indicates the dimension of the process noise vector, and m indicates the dimension of the measurement noise vector; and sigma point set
The function expression of (a) is as follows:
in the above formula, the first and second carbon atoms are,
means, P, representing state vectors xx Representing the covariance of the state vector, gamma is a coefficient,
a is the dimension of the state, λ = α 2 (a + k) -a, k is a proportionality coefficient, gamma P xx Is denoted as the ith row or column of the square root matrix
2.5 Determining the weight corresponding to each sigma point in the sigma point set of the state at the moment k:
in the above formula, W i m Denotes the mean weight, a is the state dimension, λ = α 2 (a + k) -a, k is a proportionality coefficient;
in the above formula, W i c Represents the covariance weight, W i m Denotes the mean weight, a is the state dimension, λ = α 2 (a + k) -a, k is a scaling factor and α is a mean value for determining the state vector
Coefficient of distribution of surrounding sigma points, 1e -4 ≤α<1; beta is a constant coefficient.
Optionally, when the state equation and the update equation of the adaptive Kalman filtering algorithm are determined for the sigma point set of the time k in the step 3), the method further comprises
Indicating particle
A column vector of the first n dimensions;
indicating particle
A column vector composed of n +1 dimensions to n + q dimensions;
indicating particle
The column vector formed by the dimensions from n + q +1 to n + q + m establishes the state equation of the prediction process at any k moment as shown in the formulas (12) to (14), and establishes the update equation at any k moment as shown in the formulas (15) to (21);
in the formulae (12) to (14),
represents the sigma point predicted value at the moment k +1,
represents an estimate of the sigma point of the state vector at time k, u (k) represents the input vector,
a sigma point estimation value representing a process noise vector at the k moment;
represents the predicted value of the state vector at the time k +1, W i m Means of expressionThe weight value L is the sampling number of sigma points; p (k +1, k) represents a covariance predicted value at the time of k + 1;
in formulae (15) to (21), z i (k +1, k) represents the estimated value of the measurement vector at the time k +1,
represents an estimated value of sigma point of the state vector at the moment k, u (k) is an input vector,
expressing the sigma point estimation value of the measurement noise vector at the k +1 moment;
the predicted value, W, of the measurement vector at time k +1 i m Expressing the average weight value, wherein L is the sampling number of sigma points; p is zz (k +1, k) represents the covariance of the estimated biases of the measurement vectors at time k +1, α k+1 For adaptive adjustment of the factor, W i c Representing a covariance weight; p xz (k +1, k) represents the cross-covariance of the state vector and the estimated bias of the metrology vector at time k +1,
represents the sigma point predicted value at the moment of k +1,
the state vector predicted value at the moment of k +1 is represented; k (K + 1) represents the filter gain at time K + 1;
represents the optimal estimated value of the state vector at the time k +1, and z (k + 1) represents the measured value at the time k + 1; p (k +1 ) represents the optimal estimation value of the covariance at the time k + 1.
Optionally, the factor α is adaptively adjusted k+1 Is expressed as follows:
in the above formula, tr (P) zz (k +1, k)) represents a matrix P zz Trace of (k +1, k), tr (. DELTA.Z) k+1 ΔZ k+1 T ) Representation matrix Δ Z k+1 ΔZ k+1 T Trace of (P) zz (k +1, k) represents the covariance of the estimated biases of the measurement vectors at time k +1, Δ Z k+1 The residual value is updated for the measurement time,
optionally, step 4) integrates the innovation sequence weighting into the state equation and the function expression in the update equation of the adaptive kalman filtering algorithm, as shown in equations (23) to (27):
D k+1 ＝P(k+1,k)+R k+1 (27)
in the formulae (23) to (27),
represents the best estimate of the state vector at time k +1,
represents the predicted value of the state vector at the time K +1, K (K + 1) represents the filter gain value at the time K +1, phi k+1 Representing a correction function, τ k+1 A variable representing the measurement error magnitude at the time k +1, e k+1 The measurement error at the k +1 moment is represented; z (k + 1) represents a measurement value at the time k + 1;
representing the predicted value of the measurement vector at the moment k + 1; tau is k+1 Obeying x degree of freedom of m 2 Distribution, when outliers are present, τ k+1 Increase of phi k+1 (. Decreased); c k+1 For a selected sequence of threshold values or constants, lambda k+1 Is composed of
Maximum eigenvalue of the matrix, K k+1 Represents the time k +1Kalman filter gain of, D k+1 A weight matrix for measuring the measurement error at the time k +1, P (k +1, k) represents the covariance predicted value at the time k +1, R k+1 Representing the measurement noise variance at time k + 1.
Furthermore, the present invention also provides an innovation-weighted adaptive kalman filter-based crystal frequency taming system, comprising a computer device programmed or configured to perform the steps of the innovation-weighted adaptive kalman filter-based crystal frequency taming method.
In addition, the invention also provides a crystal oscillator frequency taming system based on innovation weighted adaptive Kalman filtering, 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 crystal oscillator frequency taming method based on innovation weighted adaptive Kalman filtering.
In addition, the invention also provides an intelligent terminal with a GPS, which comprises a GPS module, a microprocessor and a memory, wherein the microprocessor is programmed or configured to execute the steps of the crystal oscillator frequency domestication method based on the innovation weighted adaptive Kalman filtering, or the memory is stored with a computer program which is programmed or configured to execute the crystal oscillator frequency domestication method based on the innovation weighted adaptive Kalman filtering.
Furthermore, the present invention also provides a computer-readable storage medium having stored thereon a computer program programmed or configured to perform the crystal frequency taming method based on innovation-weighted adaptive kalman filtering.
Compared with the prior art, the invention has the following advantages: the invention provides an adaptive unscented Kalman filtering method for weighting an innovation sequence for processing a frequency difference signal in a GPS (global positioning system) locked crystal oscillator frequency technology, which mainly comprises four parts: firstly, establishing a frequency difference mathematical model of a 1PPS signal and a crystal oscillator signal of a GPS according to the error characteristics of the two signals; secondly, determining a measurement and prediction equation of the system according to a frequency difference mathematical model, and carrying out sigma sampling; thirdly, determining a prediction equation and an updating equation of the self-adaptive insensitive Kalman algorithm; and fourthly, integrating innovation sequence weighting into the algorithm to reduce the influence of the larger field of the GPS on the frequency correction precision. Through the four steps, the adaptive insensitive Kalman filtering of innovation sequence weighting can be realized, and the method has the following technical effects: 1. processing a frequency difference signal in the process of locking the frequency of the GPS crystal oscillator by adopting an insensitive Kalman filtering algorithm, and considering the nonlinear error of the crystal oscillator caused by aging and temperature influence while processing the random error of the GPS signal; 2. the innovation sequence weighting and the self-adaptive filtering technology are introduced into the insensitive Kalman filtering process, so that the influence of large jump outliers in GPS signals and the continuous change of covariance matrixes of system noise and measurement noise on filtering results can be reduced.
Drawings
Fig. 1 is a schematic diagram of a GPS locked crystal oscillator frequency 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.
Detailed Description
As shown in fig. 2, the crystal oscillator frequency taming method based on innovation weighted adaptive kalman filtering in this embodiment includes:
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) Sigma sampling is carried out on the frequency difference mathematical model to obtain a sigma point set at the moment k;
3) Determining a state equation and an updating equation of the adaptive insensitive Kalman filtering algorithm aiming at a sigma point set at the moment k;
4) And integrating innovation sequence weighting into a state equation and an update equation of the adaptive insensitive Kalman filtering algorithm to reduce the influence of a larger field of the GPS on the frequency correction precision.
When the GPS second signal sequence Y is acquired, the universal time UTC of the GPS second signal in the GPS second signal sequence Y (with size n) can be respectively expressed as:
1-ε 1 ，2-ε 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' (with the size of n) can be respectively represented 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, and c is an error coefficient considering frequency linear drift; therefore, the detailed steps of step 1) in this embodiment include:
1.1 The function expression of the universal time UTC of the GPS second signal at any k time in the GPS second signal sequence Y is shown as the formula (1);
Y k ＝k-ε k (1)
in the above formula, Y k Universal time UTC, epsilon of GPS-seconds signal representing time k in GPS-seconds signal sequence Y k Is the error of GPS second signal at the time of k, and satisfies epsilon-N (0, sigma) 2 )；
Acquiring a crystal oscillator signal sequence Y ', wherein the functional expression of the universal time UTC of the crystal oscillator signal at any k time in the crystal oscillator signal sequence Y' is shown as a formula (2);
Y′ k ＝k+a+bk+ck 2 (2)
in the above formula, Y' k The universal time UTC of the crystal oscillator signal at the time k in the crystal oscillator signal sequence Y' is represented, 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', a frequency difference sequence is obtained, and a function expression of a frequency difference signal at any k time in the frequency difference sequence is shown as the following formula;
X k ＝Y′ k -Y k ＝a+bk+ck 2 +ε k (3)
in the above formula, X k Representing the frequency difference signal at time k in the sequence of frequency differences.
In this embodiment, the detailed steps of step 2) include:
2.1 Based on a mathematical model of frequency difference, a prediction equation for frequency difference is established:
X' k ＝X' k-1 +2ck+b-c+ω(k) (4)
in the above formula, X' k Denotes a predicted value, X ', of the frequency difference signal at time k' k-1 Representing a predicted value of the frequency difference signal at the time k-1, b being an error coefficient considering frequency deviation, c being an error coefficient considering frequency linear drift, and ω (k) representing process noise at the time k;
2.2 Let u (k) =2ck + b-c represent input vector, establish state equation and measurement equation of frequency difference:
in the above formula, x (k + 1) represents a frequency difference state vector at the time k +1, x (k) represents a frequency difference state vector at the time k, u (k) represents an input vector, and z (k + 1) represents a measurement value at the time k + 1; omega (k) represents a process noise vector at the time k, and a covariance matrix of the process noise vector is Q; gamma (k) represents the measurement noise vector at time k, and the covariance matrix is R;
2.3 Because the system has noise term, it needs to perform dimension expansion process on the state vector, and perform dimension expansion process on the frequency difference state vector x a ＝[x T ω T γ T ]Wherein x is a Representing the expanded state vector, x T Representing transposes of the state vector of the frequency difference, ω T Representing the transpose of the process noise vector, gamma T Representing the transpose of the measurement noise vector; determining an initial state of a frequency difference state vector:
in the above-mentioned formula, the compound has the following structure,
representing the initial value of the frequency difference state vector, E (x) 0 ) Means, x, representing the initial value of the frequency difference vector 0 Representing the initial value of the frequency difference vector, P 0 An initial value of a covariance vector representing a frequency difference vector,
representing the covariance of the initial value of the frequency difference vector;
determining the initial state of the state vector after the dimension expansion of the frequency difference state vector:
in the above formula, the first and second carbon atoms are,
the initial value of the state vector after the dimension expansion is shown,
represents the mean value of the state vector after the dimension expansion,
the representation represents the state vector after the dimension is expanded,
represents the initial value of the state vector of the frequency difference,
representing the initial value of the covariance of the state vector after dimension expansion,
representing the covariance of the initial value of the state vector after dimension expansion, P 0 Representing an initial value of a covariance vector of the frequency difference vector, Q representing a covariance matrix of a process noise vector omega (k) at the time k, and R representing a covariance matrix of a measurement noise vector gamma (k) at the time k;
2.4 Symmetric sigma sampling is adopted for the frequency difference state vector to obtain a sigma point set of a k moment state
i =1, \ 8230, L, wherein
Representing sigma points with the ith dimension a, L being the sampling number of the sigma points, L =2a +1, and the state dimension a = n + q + m at the moment, wherein n represents the dimension of the frequency difference state vector, q represents the dimension of the process noise vector, and m represents the dimension of the measurement noise vector; and sigma point set
The function expression of (a) is shown as follows:
in the above formula, the first and second carbon atoms are,
representing the mean value of the state vector, P xx Representing the covariance of the state vector, gamma is a coefficient,
a is the dimension of the state, λ = α 2 (a + k) -a, k is a scaling factor (the scaling factor k can be generally 0 or 3-a and can be used for adjusting the sigma point and the x mean value
Distance therebetween), γ P xx Is denoted as the ith row or column of the square root matrix
2.5 Determining the weight corresponding to each sigma point in the sigma point set of the state at the moment k:
in the above formula, W i m Denotes the mean weight, a is the state dimension, λ = α 2 (a + k) -a, k is a proportionality coefficient;
in the above formula, W i c Represents the covariance weight, W i m Denotes the mean weight, a is the state dimension, λ = α 2 (a + k) -a, k is a scaling factor and α is a mean value for determining the state vector
Coefficient of distribution of surrounding sigma points, 1e -4 ≤α<1; beta is a constant coefficient.
When the state equation and the update equation of the adaptive insensitive Kalman filtering algorithm are determined for the sigma point set at the time k in the step 3) of the embodiment, the method comprises the following steps
Indicating particle
A column vector of the first n dimensions;
indicating particle
A column vector composed of n +1 dimensions to n + q dimensions;
indicating particle
The column vector formed by the dimensions from n + q +1 to n + q + m establishes the state equation of the prediction process at any k moment as shown in the formulas (12) to (14), and establishes the update equation at any k moment as shown in the formulas (15) to (21);
in the formulae (12) to (14),
represents the sigma point predicted value at the moment k +1,
represents the estimated value of sigma point of the state vector at the moment k, u (k) represents the input vector,
a sigma point estimation value representing a process noise vector at the k moment;
represents the predicted value of the state vector at time k +1, W i m Expressing the average weight value, wherein L is the sampling number of sigma points; p (k +1, k) represents a covariance predicted value at the time of k + 1;
in formulae (15) to (21), z i (k +1, k) represents the estimated value of the measurement vector at the time k +1,
represents an estimated value of sigma point of the state vector at time k, u (k) is an input vector,
expressing the sigma point estimation value of the measurement noise vector at the k +1 moment;
the predicted value, W, of the measurement vector at time k +1 i m Expressing the average value weight, wherein L is the sampling number of sigma points; p zz (k +1, k) represents the covariance of the estimated biases of the measurement vectors at time k +1, α k+1 For adaptive adjustment of the factor, W i c Representing a covariance weight; p xz (k +1, k) represents the cross-covariance of the state vector and the estimated bias of the metrology vector at time k +1,
represents the sigma point predicted value at the moment of k +1,
the state vector predicted value at the moment of k +1 is represented; k (K + 1) represents the filter gain at time K + 1;
represents the optimal estimated value of the state vector at the time k +1, and z (k + 1) represents the measured value at the time k + 1; p (k +1 ) represents the optimal estimation value of the covariance at the time k + 1.
In this embodiment, the factor α is adaptively adjusted k+1 Is expressed as follows:
in the above formula, tr (P) zz (k +1, k)) represents a matrix P zz Trace of (k +1, k), tr (Δ Z) k+1 ΔZ k+1 T ) Representation matrix Δ Z k+1 ΔZ k+1 T Trace of (P), P zz (k +1, k) represents the covariance of the estimated bias of the measurement vector at time k +1, Δ Z k+1 The residual value is updated for the measurement time,
in practical application, because the GPS signal is interfered by noise and generates a large hop outlier, the frequency difference signal with the outlier will shift or even diverge the filtering result. In order to solve the above problem of outlier resistance of the insensitive kalman algorithm, some researchers introduce an innovation sequence weighting method into kalman filtering to eliminate the influence of outliers. When the observation data does not contain the outlier, the filtering algorithm can normally operate to filter the observation noise; when the observation data contains the outlier, the influence of the outlier on the filtering result can be overcome or reduced, and the filtering result is more accurate. In this embodiment, step 4) integrates the innovation sequence weighting into the state equation of the adaptive kalman filter algorithm and the function expression in the update equation as shown in formulas (23) to (27):
D k+1 ＝P(k+1,k)+R k+1 (27)
in the formulae (23) to (27),
represents the best estimate of the state vector at time k +1,
represents the predicted value of the state vector at the time K +1, K (K + 1) represents the filter gain value at the time K +1, phi k+1 Representing a correction function, τ k+1 A variable representing the measurement error magnitude at the time k +1, e k+1 The measurement error at the k +1 moment is represented; z (k + 1) represents a measurement value at the time k + 1;
representing the predicted value of the measurement vector at the moment k + 1; tau. k+1 Obeying x degree of freedom of m 2 Distribution, when outliers are present, τ k+1 Increase of phi k+1 (. O) decreases; c k+1 For a selected sequence of thresholds or constants, lambda k+1 Is composed of
Maximum eigenvalue of the matrix, K k+1 Kalman Filter gain, D, representing the k +1 time instant k+1 A weight matrix for measuring the measurement error at the time k +1, P (k +1, k) represents the covariance predicted value at the time k +1, R k+1 Representing the measurement noise variance at time k + 1.
To sum up, the main technology currently applied to the low-cost high-precision frequency standard is a GPS-locked secondary frequency standard (crystal oscillator) technology, and in order to reduce the random error and the hop outlier of a GPS signal in the GPS-locked secondary frequency standard technology and the influence of the nonlinear frequency offset of a crystal oscillator signal on a frequency offset signal, the method of this embodiment provides an adaptive insensitive kalman filtering method for weighting an innovation sequence for processing the frequency offset signal in the GPS-locked crystal oscillator frequency technology, and the method mainly includes four parts: firstly, establishing a frequency difference mathematical model of a 1PPS signal and a crystal oscillator signal of a GPS according to the error characteristics of the two signals; secondly, determining a measurement and prediction equation of the system according to the frequency difference mathematical model, and carrying out sigma sampling; thirdly, determining a prediction equation and an updating equation of the self-adaptive insensitive Kalman algorithm; and fourthly, integrating innovation sequence weighting into the algorithm to reduce the influence of the larger field of the GPS on the frequency correction precision. Through the four steps, the self-adaptive insensitive Kalman filtering of innovation sequence weighting can be realized, and the method has the following technical effects: 1. processing a frequency difference signal in the process of locking the frequency of the GPS crystal oscillator by adopting an insensitive Kalman filtering algorithm, and considering the nonlinear error of the crystal oscillator caused by aging and temperature influence while processing the random error of the GPS signal; 2. the innovation sequence weighting and the self-adaptive filtering technology are introduced into the insensitive Kalman filtering process, so that the influence of large jump outliers in GPS signals and the continuous change of covariance matrixes of system noise and measurement noise on filtering results can be reduced.
Furthermore, the present embodiment also provides a crystal oscillator frequency taming system based on innovation weighted adaptive kalman filtering, which includes a computer device, wherein the computer device is programmed or configured to execute the steps of the crystal oscillator frequency taming method based on innovation weighted adaptive kalman filtering as claimed in any one of claims 1 to 5.
In addition, the present embodiment also provides an innovation-weighted adaptive kalman filter-based crystal frequency taming system, which includes a computer device having a memory storing thereon a computer program programmed or configured to perform the aforementioned innovation-weighted adaptive kalman filter-based crystal frequency taming method.
In addition, the present embodiment also provides an intelligent terminal with GPS, which includes a GPS module, a microprocessor and a memory, wherein the microprocessor is programmed or configured to execute the steps of the aforementioned crystal oscillator frequency domestication method based on innovation weighted adaptive kalman filtering, or the memory stores thereon a computer program programmed or configured to execute the aforementioned crystal oscillator frequency domestication method based on innovation weighted adaptive kalman filtering.
Furthermore, the present embodiment also provides a computer-readable storage medium having stored thereon a computer program programmed or configured to execute the aforementioned method for crystal frequency taming based on innovation-weighted adaptive kalman filtering.
As will be appreciated by one skilled in the art, 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-readable 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 according to embodiments of the application wherein instructions, which execute via a flowchart and/or a processor of the computer program product, create means for implementing functions specified in the flowchart and/or block diagram block or blocks. 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 embodiments, and all technical solutions belonging to the idea of the present invention belong to the protection scope of the present invention. It should be noted that modifications and embellishments within the scope of the invention may occur to those skilled in the art without departing from the principle of the invention, and are considered to be within the scope of the invention.
Claims (9)
1. A crystal oscillator frequency taming method based on innovation weighted adaptive insensitive Kalman filtering is characterized by comprising the following steps:
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) Sigma sampling is carried out aiming at the frequency difference mathematical model to obtain a sigma point set at the k moment;
3) Determining a state equation and an updating equation of the adaptive insensitive Kalman filtering algorithm aiming at a sigma point set at the moment k;
4) Integrating innovation sequence weighting into a state equation and an update equation of the adaptive insensitive Kalman filtering algorithm to reduce the influence of a larger field of the GPS on frequency correction precision;
when the state equation and the update equation of the adaptive insensitive Kalman filtering algorithm are determined aiming at the sigma point set of the moment k in the step 3), so as to
Representing particles
A column vector of the first n dimensions;
indicating particle
A column vector consisting of n +1 dimension to n + q dimension;
indicating particle
The column vector formed by the dimensions from n + q +1 to n + q + m establishes the state equation of the prediction process at any k moment as shown in the formulas (12) to (14), and establishes the update equation at any k moment as shown in the formulas (15) to (21);
in the formulae (12) to (14),
represents the sigma point predicted value at the moment k +1,
represents the estimated value of sigma point of the state vector at the moment k, u (k) represents the input vector,
a sigma point estimation value representing a process noise vector at the k moment;
representing the state vector at time k +1Predicted value, W i m Expressing the average weight value, wherein L is the sampling number of sigma points; p (k +1, k) represents a covariance predicted value at the time of k + 1;
in formulae (15) to (21), z i (k +1, k) represents the estimated value of the measurement vector at the time k +1,
represents an estimated value of sigma point of the state vector at the moment k, u (k) is an input vector,
expressing the sigma point estimation value of the measurement noise vector at the k +1 moment;
represents the predicted value, W, of the measurement vector at time k +1 i m Expressing the average weight value, wherein L is the sampling number of sigma points; p zz (k +1, k) represents the covariance of the estimated biases of the measurement vectors at time k +1, α k+1 For adaptive adjustment of the factor, W i c Representing a covariance weight; p xz (k +1, k) represents the cross-covariance of the state vector and the estimated bias of the metrology vector at time k +1,
represents the sigma point predicted value at the moment k +1,
the state vector predicted value at the moment of k +1 is represented; k (K + 1) represents the filter gain at time K + 1;
represents the optimal estimated value of the state vector at the time k +1, and z (k + 1) represents the measured value at the time k + 1; p (k +1 ) represents the optimal estimation value of the covariance at the time k + 1.
2. The crystal oscillator frequency taming method based on innovation weighted adaptive Kalman filtering is characterized in that the detailed steps of the step 1) comprise:
1.1 The function expression of the universal time UTC of the GPS second signal at any k time in the GPS second signal sequence Y is shown as the formula (1);
Y k ＝k-ε k (1)
in the above formula, Y k Universal time UTC, epsilon of GPS-seconds signal representing time k in GPS-seconds signal sequence Y k Is the error of GPS second signal at the time of k, and satisfies epsilon-N (0, sigma) 2 )；
Acquiring a crystal oscillator signal sequence Y ', wherein the functional expression of the universal time UTC of the crystal oscillator signal at any k time in the crystal oscillator signal sequence Y' is shown as a formula (2);
Y′ k ＝k+a+bk+ck 2 (2)
in the above formula, Y' k The universal time UTC of the crystal oscillator signal at the time k in the crystal oscillator signal sequence Y' is represented, wherein 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', a frequency difference sequence is obtained, and a function expression of a frequency difference signal at any k time in the frequency difference sequence is shown as the following formula;
X k ＝Y′ k -Y k ＝a+bk+ck 2 +ε k (3)
in the above formula, X k Representing the frequency difference signal at time k in the sequence of frequency differences.
3. The crystal oscillator frequency taming method based on innovation weighted adaptive Kalman filtering is characterized in that the detailed steps of the step 2) comprise:
2.1 Based on a mathematical model of frequency difference, a prediction equation for frequency difference is established:
X' k ＝X' k-1 +2ck+b-c+ω(k) (4)
in the above formula, X' k Denotes a predicted value, X ', of the frequency difference signal at time k' k-1 Representing a predicted value of the frequency difference signal at the time k-1, b being an error coefficient considering frequency deviation, c being an error coefficient considering frequency linear drift, and ω (k) representing process noise at the time k;
2.2 Let u (k) =2ck + b-c represent input vector, establish state equation and measurement equation of frequency difference:
in the above equation, x (k + 1) represents a frequency offset state vector at time k +1, x (k) represents a frequency offset state vector at time k, u (k) represents an input vector, and z (k + 1) represents a measurement value at time k + 1; omega (k) represents a process noise vector at the time k, and a covariance matrix of the process noise vector is Q; gamma (k) represents the measurement noise vector at time k, and the covariance matrix is R;
2.3 X) performing dimension expansion on the frequency difference state vector a ＝[x T ω T γ T ]Wherein x is a Representing the expanded state vector, x T Representing transposes of the state vector of the frequency difference, ω T Representing the transpose of the process noise vector, gamma T Representing the transpose of the measurement noise vector; determining an initial state of the frequency difference state vector:
in the above-mentioned formula, the compound has the following structure,
representing the initial value of the frequency difference state vector, E (x) 0 ) Means, x, representing the initial value of the frequency difference vector 0 Representing the initial value of the frequency difference vector, P 0 An initial value of a covariance vector representing a frequency difference vector,
representing the covariance of the initial value of the frequency difference vector;
determining the initial state of the state vector after the dimension expansion of the frequency difference state vector:
in the above formula, the first and second carbon atoms are,
the initial value of the state vector after the dimension expansion is shown,
presentation expanderThe mean value of the state vector after the dimension is measured,
the representation represents the state vector after the dimension expansion,
represents the initial value of the state vector of the frequency difference,
representing the initial value of the covariance of the state vector after dimension expansion,
representing the covariance of the initial value of the state vector after dimension expansion, P 0 Representing an initial value of a covariance vector of the frequency difference vector, Q representing a covariance matrix of a process noise vector omega (k) at the time k, and R representing a covariance matrix of a measurement noise vector gamma (k) at the time k;
2.4 Symmetric sigma sampling is adopted for the frequency difference state vector to obtain a sigma point set of a k moment state
i =1, \ 8230, L, wherein
Indicating sigma points with the ith dimension a, L indicating the sampling number of the sigma points, L =2a +1, and the state dimension a = n + q + m, wherein n indicates the dimension of the frequency difference state vector, q indicates the dimension of the process noise vector, and m indicates the dimension of the measurement noise vector; and sigma point set
The function expression of (a) is shown as follows:
in the above formula, the first and second carbon atoms are,
representing the mean value of the state vector, P xx Representing the covariance of the state vector, gamma is a coefficient,
a is the dimension of the state, λ = α 2 (a + k) -a, k is a proportionality coefficient, gamma P xx Is denoted as the ith row or column of the square root matrix
2.5 Determining the weight corresponding to each sigma point in the sigma point set of the state at the moment k:
in the above formula, W i m Denotes the mean weight, a is the state dimension, λ = α 2 (a + k) -a, k is a proportionality coefficient;
in the above formula, W i c Represents the covariance weight, W i m Denotes the mean weight, a is the state dimension, λ = α 2 (a + k) -a, k is the scaling factor and a is the mean value used to determine the state vector
Coefficient of distribution of surrounding sigma points, 1e -4 ≤α<1; beta is a constant coefficient.
4. The crystal oscillator frequency taming method based on innovation-weighted adaptive Kalman filtering, as per claim 3, characterized in that the adaptationShould adjust the factor alpha k+1 Is expressed as follows:
in the above formula, tr (P) zz (k +1, k)) represents a matrix P zz Trace of (k +1, k), tr (Δ Z) k+1 ΔZ k+1 T ) Representation matrix Δ Z k+1 ΔZ k+1 T Trace of (P), P zz (k +1, k) represents the covariance of the estimated bias of the measurement vector at time k +1, Δ Z k+1 The residual value is updated for the measurement time,
5. the crystal oscillator frequency taming method based on innovation weighted adaptive kalman filtering according to claim 3, characterized in that step 4) integrates innovation sequence weighting into the functional expression in the state equation and the update equation of the adaptive kalman filtering algorithm as shown in formulas (23) to (27):
D k+1 ＝P(k+1,k)+R k+1 (27)
in the formulae (23) to (27),
represents the best estimate of the state vector at time k +1,
represents the predicted value of the state vector at the time K +1, K (K + 1) represents the filter gain value at the time K +1, phi k+1 Representing a correction function, τ k+1 A variable representing the measurement error magnitude at the time k +1, e k+1 The measurement error at the k +1 moment is represented; z (k + 1) represents a measurement value at the time k + 1;
representing the predicted value of the measurement vector at the k +1 moment; tau is k+1 Obeying x degree of freedom of m 2 Distribution, when outliers are present, [ tau ] k+1 Increase of phi k+1 (. O) decreases; c k+1 For a selected sequence of threshold values or constants, lambda k+1 Is composed of
Maximum eigenvalue of the matrix, K k+1 Kalman Filter gain, D, representing the k +1 time instant k+1 A weight matrix for measuring the measurement error at the time k +1, P (k +1, k) represents the covariance predicted value at the time k +1, R k+1 Representing the measurement noise variance at time k + 1.
6. A crystal frequency taming system based on innovation weighted adaptive Kalman filtering, comprising a computer device, characterized in that the computer device is programmed or configured to perform the steps of the crystal frequency taming method based on innovation weighted adaptive Kalman filtering as claimed in any one of claims 1 to 5.
7. A crystal frequency taming system based on innovation weighted adaptive Kalman filtering, comprising a computer device, characterized in that the memory of the computer device has stored thereon a computer program programmed or configured to perform the crystal frequency taming method based on innovation weighted adaptive Kalman filtering as claimed in any one of claims 1 to 5.
8. 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 crystal oscillator frequency discipline method based on innovation weighted adaptive kalman filtering according to any one of claims 1 to 5, or the memory has stored thereon a computer program programmed or configured to perform the crystal oscillator frequency discipline method based on innovation weighted adaptive kalman filtering according to any one of claims 1 to 5.
9. A computer-readable storage medium having stored thereon a computer program programmed or configured to perform the crystal frequency taming method based on innovation-weighted adaptive kalman filtering as claimed in any one of claims 1 to 5.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010524159.9A CN111650617B (en) | 2020-06-10 | 2020-06-10 | Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010524159.9A CN111650617B (en) | 2020-06-10 | 2020-06-10 | Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111650617A CN111650617A (en) | 2020-09-11 |
| CN111650617B true CN111650617B (en) | 2023-03-03 |
Family
ID=72343261
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010524159.9A Active CN111650617B (en) | 2020-06-10 | 2020-06-10 | Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111650617B (en) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112578415B (en) * | 2020-11-06 | 2023-10-13 | 中国科学院国家空间科学中心 | Digital frequency locking method and loop based on adaptive filter |
| CN112380774B (en) * | 2020-11-23 | 2022-04-15 | 青岛柯锐思德电子科技有限公司 | Dynamic modeling method and system based on residual echo state network |
| CN112713881B (en) * | 2020-12-10 | 2022-11-01 | 国网四川省电力公司电力科学研究院 | Synchronous clock maintaining system and method based on edge calculation |
| CN113419414B (en) * | 2021-07-13 | 2023-03-21 | 贵州省计量测试院 | Standard timing system with GNSS disciplining and interval time keeping capabilities |
| CN113959433B (en) * | 2021-09-16 | 2023-12-08 | 南方电网数字平台科技(广东)有限公司 | Combined navigation method and device |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6711230B1 (en) * | 2002-09-27 | 2004-03-23 | Nortel Networks Limited | Reference timing signal oscillator with frequency stability |
| US8643444B2 (en) * | 2012-06-04 | 2014-02-04 | Broadcom Corporation | Common reference crystal systems |
| KR20140070722A (en) * | 2012-11-26 | 2014-06-11 | 한국전자통신연구원 | Integration apparatus for multi-rate system and method thereof |
| JP2018165618A (en) * | 2017-03-28 | 2018-10-25 | セイコーエプソン株式会社 | Signal processing device, detection device, physical quantity measurement device, electronic apparatus and mobile body |
| CN110336279B (en) * | 2019-07-17 | 2020-11-20 | 国网湖南省电力有限公司 | Electric power system oscillation self-adaptive suppression method, system and medium based on grid-connected converter |
-
2020
- 2020-06-10 CN CN202010524159.9A patent/CN111650617B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN111650617A (en) | 2020-09-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111650617B (en) | Crystal oscillator frequency taming method, system and medium based on innovation weighted adaptive Kalman filtering | |
| CN106026919B (en) | The punctual compensation method of crystal oscillator | |
| US4899117A (en) | High accuracy frequency standard and clock system | |
| CN103033186B (en) | High-precision integrated navigation positioning method for underwater glider | |
| Weiss et al. | AT2, a new time scale algorithm: AT1 plus frequency variance | |
| CN110554597B (en) | Hydrogen cesium time scale fusion method based on Vondark-Cepek filtering | |
| CN109508510B (en) | Improved Kalman filtering-based rubidium atomic clock parameter estimation algorithm | |
| WO1995034850B1 (en) | Assured-integrity monitored-extrapolation navigation apparatus | |
| CN109581856B (en) | Time synchronization and time keeping method based on high-performance crystal oscillator frequency calibration | |
| CN101562451A (en) | Precise domestication conserving method of second-level frequency scale | |
| CN107465393B (en) | System and method for frequency compensation of real time clock system | |
| US10033390B2 (en) | Systems and methods for clock synchronization in a data acquisition system | |
| CN111694040B (en) | Positioning method and device for satellite/inertial integrated navigation system | |
| CN105718642B (en) | A kind of reference time scale production method based on threshold autoregressive model | |
| WO1985001807A1 (en) | Selective parametric self-calibrating control system | |
| US10028239B2 (en) | Method and system for precise temperature and timebase PPM error estimation using multiple timebases | |
| CN107607977B (en) | Self-adaptive UKF (unscented Kalman Filter) integrated navigation method based on minimum skewness simplex sampling | |
| Narasimhappa et al. | An innovation based random weighting estimation mechanism for denoising fiber optic gyro drift signal | |
| CN114584136A (en) | GNSS-based high-stability crystal oscillator taming and maintaining system and method | |
| Guinot | Some properties of algorithms for atomic time scales | |
| CN111711446B (en) | Method, system and medium for taming crystal oscillator frequency by using GPS signal | |
| CN112505386B (en) | Method and system for detecting current value of direct current charging pile | |
| US11035902B2 (en) | Advanced fuel gauge | |
| CN109655057B (en) | Filtering optimization method and system for accelerator measurement value of six-push unmanned aerial vehicle | |
| CN107421543B (en) | Implicit function measurement model filtering method based on state dimension expansion |
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 |