CN116847452B - Satellite signal timing synchronization system - Google Patents
Satellite signal timing synchronization system Download PDFInfo
- Publication number
- CN116847452B CN116847452B CN202310966519.4A CN202310966519A CN116847452B CN 116847452 B CN116847452 B CN 116847452B CN 202310966519 A CN202310966519 A CN 202310966519A CN 116847452 B CN116847452 B CN 116847452B
- Authority
- CN
- China
- Prior art keywords
- monitoring object
- clock
- clock difference
- time
- gene
- 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
- 238000012544 monitoring process Methods 0.000 claims abstract description 96
- 239000011159 matrix material Substances 0.000 claims abstract description 59
- 238000000034 method Methods 0.000 claims abstract description 34
- 239000013078 crystal Substances 0.000 claims abstract description 32
- 238000001914 filtration Methods 0.000 claims abstract description 29
- 238000004364 calculation method Methods 0.000 claims abstract description 24
- 238000005070 sampling Methods 0.000 claims abstract description 14
- 238000012937 correction Methods 0.000 claims abstract description 9
- 108090000623 proteins and genes Proteins 0.000 claims description 76
- 210000000349 chromosome Anatomy 0.000 claims description 24
- 230000002068 genetic effect Effects 0.000 claims description 11
- 230000035772 mutation Effects 0.000 claims description 11
- 238000005259 measurement Methods 0.000 claims description 8
- 238000009826 distribution Methods 0.000 claims description 7
- 238000005315 distribution function Methods 0.000 claims description 4
- 230000007704 transition Effects 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 4
- 230000001186 cumulative effect Effects 0.000 claims description 3
- 238000005516 engineering process Methods 0.000 abstract description 6
- 238000004891 communication Methods 0.000 abstract description 3
- 238000011161 development Methods 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 239000010453 quartz Substances 0.000 description 1
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N silicon dioxide Inorganic materials O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04W—WIRELESS COMMUNICATION NETWORKS
- H04W56/00—Synchronisation arrangements
- H04W56/001—Synchronization between nodes
-
- G—PHYSICS
- G04—HOROLOGY
- G04R—RADIO-CONTROLLED TIME-PIECES
- G04R20/00—Setting the time according to the time information carried or implied by the radio signal
- G04R20/02—Setting the time according to the time information carried or implied by the radio signal the radio signal being sent by a satellite, e.g. GPS
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/12—Computing arrangements based on biological models using genetic models
- G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04J—MULTIPLEX COMMUNICATION
- H04J3/00—Time-division multiplex systems
- H04J3/02—Details
- H04J3/06—Synchronising arrangements
- H04J3/0635—Clock or time synchronisation in a network
- H04J3/0638—Clock or time synchronisation among nodes; Internode synchronisation
- H04J3/0644—External master-clock
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04W—WIRELESS COMMUNICATION NETWORKS
- H04W84/00—Network topologies
- H04W84/02—Hierarchically pre-organised networks, e.g. paging networks, cellular networks, WLAN [Wireless Local Area Network] or WLL [Wireless Local Loop]
- H04W84/04—Large scale networks; Deep hierarchical networks
- H04W84/06—Airborne or Satellite Networks
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Computational Mathematics (AREA)
- Biophysics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Life Sciences & Earth Sciences (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Signal Processing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physiology (AREA)
- Genetics & Genomics (AREA)
- Artificial Intelligence (AREA)
- Biomedical Technology (AREA)
- Computational Linguistics (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention relates to the technical field of communication, in particular to a satellite signal timing synchronization system, which comprises: and a data sampling module: the method comprises the steps of sampling data of a first monitoring object, a second monitoring object and a third monitoring object; model building module: the observation model is used for respectively building the second monitoring object and the third monitoring object to observe the first monitoring object based on the crystal oscillator clock; and a clock difference calculation module: the method comprises the steps of calculating clock difference based on a built observation model, and calculating and outputting optimal estimated clock difference through a nonlinear Kalman filtering algorithm; clock correction module: for correcting the synchronization time of the first monitoring object according to the optimal estimated clock difference. According to the invention, time service is performed on the reference satellite and the base station through the crystal oscillator, then the clock difference between the satellites is calculated, so that a clock difference matrix updated in real time is built, the optimal estimated clock difference is obtained based on a nonlinear Kalman filtering algorithm to perform satellite time synchronization, the calculation time is reduced, the accuracy of measured data is improved, and various technologies can be developed in time conveniently.
Description
Technical Field
The invention relates to the technical field of communication, in particular to a satellite signal timing synchronization system.
Background
With the wide application of Beidou system time service products in the whole country, high-precision time service is increasingly wide and popular. While current satellite timing technology can meet most of the demands in work and life, some fine-working fields require microsecond or even nanosecond time accuracy. Because the inside of each satellite device is not provided with a real-time clock design, the self-time keeping function is not realized, the time between the satellite and the base station cannot be synchronized, and the measured data has larger error. In the prior art, time delay is reduced by performing time service on satellites and base stations through a crystal oscillator time service method, and the accuracy of measured data is greatly improved, however, along with the rapid development of aerospace technology, each item of data increasingly depends on a satellite system, time synchronization cannot be achieved among satellites, so that great inconvenience is brought to the development of each technology, and a great amount of time is consumed for performing time service calculation on each satellite respectively, so that the development of each technology is not facilitated.
Disclosure of Invention
The present invention aims to solve the above-mentioned drawbacks of the prior art by providing a satellite signal timing synchronization system.
The technical scheme adopted by the invention is as follows:
providing a satellite signal timing synchronization system, comprising:
and a data sampling module: the method comprises the steps of sampling data of a first monitoring object, a second monitoring object and a third monitoring object;
model building module: the observation model is used for respectively building the second monitoring object and the third monitoring object to observe the first monitoring object based on the crystal oscillator clock;
and a clock difference calculation module: the method comprises the steps of calculating clock difference based on a built observation model, and calculating and outputting optimal estimated clock difference through a nonlinear Kalman filtering algorithm;
clock correction module: for correcting the synchronization time of the first monitoring object according to the optimal estimated clock difference.
As a preferred technical scheme of the invention: the data sampling module samples the first monitoring object data, the second monitoring object data and the third monitoring object data every interval preset time.
As a preferred technical scheme of the invention: in the model building module, the model building module carries out time service on a first monitoring object and a second monitoring object based on a crystal oscillator clock, and the model building module specifically comprises the following steps:
setting a crystal oscillator clock to accurately time a first monitoring object, wherein the initial deviation of time synchronization of the crystal oscillator clock and the international standard time and the initial deviation of a second period are respectively a and b, and the fixed offset error in each second is tau; deviation deltat of x second first monitoring object from international standard time x The method comprises the following steps:
Δt x ＝b+τx
cumulative deviation delta between the x second of the first monitoring object and the international standard time 1 (x) The method comprises the following steps:
the time error of the second monitoring object obeys the normal distribution mu-N (0, sigma) 2 ) Wherein μ is a second monitoring object time random error, N is a normal distribution function, and σ is a standard deviation;
the x second clock deviation of the second monitor object output is:
δ 2 (x)＝x-μ x
wherein mu x A clock error of the second monitoring object representing the x second;
obtaining a deviation model delta (x):
performing deviation calculation based on the regression model to obtain a deviation value of the satellite clock and the crystal oscillator clock;
the second monitoring object and the third monitoring object observe the clock signal of the first monitoring object at the same time, and an observation model is built as follows:
h 1 ＝L 1 +ct-cT 1 +d 1 +d 2
h 2 ＝L 2 +ct-cT 2 +d 1 +d 2
wherein h is 1 And h 2 Observed amounts of the second monitoring object and the third monitoring object respectively, L 1 Pseudo-range from first monitoring object to second monitoring object, L 2 Pseudo range from the first monitoring object to the third monitoring object, c is signal propagation speed, T is clock difference of the first monitoring object, T 1 And T 2 Clock difference, d, of the second and third monitoring objects respectively 1 And d 2 Tropospheric delay and ionospheric delay, respectively.
As a preferred technical scheme of the invention: the clock difference calculation module performs difference based on an observation model of the model building module to obtain a difference observed quantity delta h of the second monitoring object and the third monitoring object as follows:
Δh＝ΔL-cΔT
wherein Δl is the differential pseudorange;
the method comprises the following steps:
wherein Δt is the differential clock difference.
As a preferred technical scheme of the invention: the clock difference calculation module is further used for building a clock difference matrix based on the acquired clock difference data, performing clock difference filtering based on a nonlinear Kalman filtering algorithm and outputting an optimal estimated clock difference.
As a preferred technical scheme of the invention: and the clock difference matrix is updated in real time according to the acquired and calculated clock difference data.
As a preferred technical scheme of the invention: in the nonlinear Kalman filtering algorithm, a clock difference matrix is used as a state variable X k+1 The method is characterized by comprising the following steps:
X k+1 ＝f(X k )+V k
Z k+1 ＝H k+1 X k+1 +W k+1
wherein f (X) k ) As a state function, Z k+1 To measure the vector, H k+1 To observe the matrix, V k And W is k+1 System noise and measurement noise, respectively;
the state one-step prediction is:
the state one-step prediction mean square error is:
the filtering gain is:
the state estimation is:
the state estimation mean square error is:
P k+1，k+1 ＝(1-K k+1 H k+1 )P k+1，k
wherein,a clock difference prediction value indicating the next time; />Representing an optimal clock difference estimated value at the current moment; p (P) k+1，k A predicted mean square error matrix for the next moment; p (P) k，k Representing an estimated mean square error matrix at the current time; phi (phi) k+1，k Is a state transition matrix; j (J) k+1，k Is a Jacobian matrix; />Is a transpose of the jacobian matrix; q (Q) k Representing a process noise covariance; k (K) k+1 Representing the Kalman filtering gain; />Is the transposition of the observation matrix; r is R k+1 Is observed noise covariance; p (P) k+1，k+1 Representing an estimation error matrix; />The state clock difference optimal estimated value at the next time is represented.
As a preferred technical scheme of the invention: in the clock difference calculation module, the layout position of the first monitoring object is optimized and obtained based on an improved genetic algorithm.
As a preferred technical scheme of the invention: the improved genetic algorithm is specifically as follows:
the population is formed by combining M chromosomes, the M chromosomes contain n gene sequences, the gene sequences are subjected to crossover and mutation for iteration, a clock error function is an fitness function value, and the probability of being selected is higher as the fitness value is higher;
probability of selection p for M chromosomes m The following are provided;
wherein f m Is the reciprocal of fitness value of the m-th chromosome.
As a preferred technical scheme of the invention: in the improved genetic algorithm, the crossover probability p of genes on M chromosomes a The method comprises the following steps:
wherein p is a Representing crossover probability of genes, p amax Represents the maximum crossover probability of a gene, p amin Representing the minimum crossover probability of a gene, f ai Indicating the fitness of the gene, f amin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the probability of variation of genes on the M chromosomes is:
wherein p is b Representing the probability of variation of the gene, p bmax Represents the maximum mutation probability of the gene, p bmin Representing the minimum probability of variation of the gene, f bi Indicating the fitness of the gene, f bmin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the fitness function value of the ith gene in the n gene sequences is f i The average fitness function value of the current gene isThe standard deviation ρ is defined as follows:
population standard deviation threshold value theta, maximum iteration times is R, times of standard deviation smaller than threshold value are C, C=C+1 when rho is smaller than or equal to theta, C is reset to 0 when rho is larger than theta, and algorithm is iterated toWhen the rho is less than or equal to theta in the continuous P iterations of the algorithm, gaussian disturbance is carried out on the current population data, and the variation probability P is increased b Let the current mutation probability p b Becomes 2p b 。
Compared with the prior art, the satellite signal timing synchronization system provided by the invention has the beneficial effects that:
according to the invention, the optimal base station position is obtained based on optimization of an improved genetic algorithm, time service is carried out on a reference satellite and the optimal base station through crystal oscillator, then clock difference among satellites is calculated, a clock difference matrix updated in real time is built, the optimal estimated clock difference is obtained based on a nonlinear Kalman filtering algorithm to carry out satellite time synchronization, the calculation time is reduced, the accuracy of measurement data is improved, and various technologies can be conveniently developed in time.
Drawings
Fig. 1 is a system block diagram of a preferred embodiment of the present invention.
The meaning of each label in the figure is: 100. a data sampling module; 200. a model building module; 300. a clock difference calculation module; 400. and a clock correction module.
Detailed Description
It should be noted that, under the condition of no conflict, the embodiments of the present embodiments and features in the embodiments may be combined with each other, and the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and obviously, the described embodiments are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
Referring to fig. 1, a preferred embodiment of the present invention provides a satellite signal timing synchronization system comprising:
the data sampling module 100: the method comprises the steps of sampling data of a first monitoring object, a second monitoring object and a third monitoring object;
model building module 200: the observation model is used for respectively building the second monitoring object and the third monitoring object to observe the first monitoring object based on the crystal oscillator clock;
clock difference calculation module 300: the method comprises the steps of calculating clock difference based on a built observation model, and calculating and outputting optimal estimated clock difference through a nonlinear Kalman filtering algorithm;
clock correction module 400: for correcting the synchronization time of the first monitoring object according to the optimal estimated clock difference.
The data sampling module 100 samples the first, second and third monitor object data every preset time interval.
Crystal oscillator time service is a process of transmitting an accurate clock signal to a device or system requiring synchronization by using a crystal oscillator (crystal oscillator) as a clock source. Crystal oscillator is an oscillator based on quartz crystal, and has stable oscillation frequency and high precision time reference. In crystal time service, a crystal is used as the primary clock source and the oscillating signal it generates is transferred to the target device or system for synchronizing its internal clocks. The frequency of the crystal is typically very stable and can provide a highly accurate time reference. Crystal oscillator time service can be used in various applications such as electronic equipment, communication systems, computer networks and the like.
The model building module 200 performs time service on the first monitoring object and the second monitoring object based on the crystal oscillator clock, and specifically includes the following steps:
setting a crystal oscillator clock to accurately time a first monitoring object, wherein the initial deviation of time synchronization of the crystal oscillator clock and the international standard time and the initial deviation of a second period are respectively a and b, and the fixed offset error in each second is tau; deviation deltat of x second first monitoring object from international standard time x The method comprises the following steps:
Δt x ＝b+τx
cumulative deviation delta between the x second of the first monitoring object and the international standard time 1 (x) The method comprises the following steps:
the time error of the second monitoring object obeys the normal distribution mu-N (0, sigma) 2 ) Wherein μ is a second monitoring object time random error, N is a normal distribution function, and σ is a standard deviation;
the x second clock deviation of the second monitor object output is:
δ 2 (x)＝x-μ x
wherein mu x A clock error of the second monitoring object representing the x second;
obtaining a deviation model delta (x):
performing deviation calculation based on the regression model to obtain a deviation value of the satellite clock and the crystal oscillator clock;
the second monitoring object and the third monitoring object observe the clock signal of the first monitoring object at the same time, and an observation model is built as follows:
h 1 ＝L 1 +ct-cT 1 +d 1 +d 2
h 2 ＝L 2 +ct-cT 2 +d 1 +d 2
wherein h is 1 And h 2 Observed amounts of the second monitoring object and the third monitoring object respectively, L 1 Pseudo-range from first monitoring object to second monitoring object, L 2 Pseudo range from the first monitoring object to the third monitoring object, c is signal propagation speed, T is clock difference of the first monitoring object, T 1 And T 2 Clock difference, d, of the second and third monitoring objects respectively 1 And d 2 Tropospheric delay and ionospheric delay, respectively.
The clock difference calculation module 300 performs difference based on the observation model of the model building module 200, and obtains a differential observed quantity Δh of the second monitored object and the third monitored object as follows:
Δh＝ΔL-cΔT
wherein Δl is the differential pseudorange;
the method comprises the following steps:
wherein Δt is the differential clock difference.
The clock difference calculation module 300 also builds a clock difference matrix based on the acquired clock difference data, performs clock difference filtering based on a nonlinear kalman filtering algorithm, and outputs an optimal estimated clock difference.
And the clock difference matrix is updated in real time according to the acquired and calculated clock difference data.
In the nonlinear Kalman filtering algorithm, a clock difference matrix is used as a state variable X k+1 The method is characterized by comprising the following steps:
X k+1 ＝f(X k )+V k
Z k+1 ＝H k+1 X k+1 +W k+1
wherein f (X) k ) As a state function, Z k+1 To measure the vector, H k+1 To observe the matrix, V k And W is k+1 System noise and measurement noise, respectively;
the state one-step prediction is:
the state one-step prediction mean square error is:
the filtering gain is:
the state estimation is:
the state estimation mean square error is:
P k+1，k+1 ＝(1-K k+1 H k+1 )P k+1，k
wherein,a clock difference prediction value indicating the next time; />Representing an optimal clock difference estimated value at the current moment; p (P) k+1，k A predicted mean square error matrix for the next moment; p (P) k，k Representing an estimated mean square error matrix at the current time; phi (phi) k+1，k Is a state transition matrix; j (J) k+1，k Is a Jacobian matrix; />Is a transpose of the jacobian matrix; q (Q) k Representing a process noise covariance; k (K) k+1 Representing the Kalman filtering gain; />Is the transposition of the observation matrix; r is R k+1 Is observed noise covariance; p (P) k+1，k+1 Representing an estimation error matrix; />The state clock difference optimal estimated value at the next time is represented.
In this embodiment, the inventors introduced adaptive techniques to dynamically adjust the parameters of the filtering algorithm due to the non-linear nature of the clock bias and the non-specific distribution of noise. Specifically, the noise covariance matrix of the Kalman filtering algorithm is adaptively updated based on real-time observation data to adapt to different noise distributions and changes.
For each instant k;
1. and (5) performing observation updating:
acquisition ofMeasuring vector Z k+1 。
Calculating an observed noise covariance matrix R k+1 。
Calculating Kalman filter gain K k+1 ：
Calculating intermediate variablesWherein H is k+1 Is an observation matrix.
Calculating Kalman filtering gain:
2. and (3) performing state estimation:
obtaining observed quantity Z k+1 。
Calculating a predicted value of the observed quantity
Updating the state clock difference estimated value:
updating a prediction error covariance matrix: p (P) k+1，k+1 ＝(I-K k+1 H k+1 )P k+1，k Where I is the identity matrix.
3. Updating the observed noise covariance matrix:
according to the real-time observation data and the current state clock difference estimated value, calculating an observation noise covariance matrix R k+1 。
Using least square estimation, estimating a new observed noise covariance matrix according to an observed residual and a measured noise model, wherein the method specifically comprises the following steps:
defining an observation residual:
in time step k+1, according to observed quantity Z k+1 And predicted valueCalculating an observation residual e k+1 ：
Establishing a measurement noise model:
assume that the observed residual follows a normal distribution of zero mean and has a measurement noise covariance matrix R k+1 ：
e k+1 ～N(0，R k+1 )
Least squares estimation:
for each observation residual e k+1 Calculate its covariance matrix P e，k+1 ：
P e，k+1 ＝cov(e k+1 )
Wherein cov represents covariance.
Updating the observed noise covariance matrix:
using the estimated observation residual covariance matrix P e，k+1 To update the observed noise covariance matrix R k+1 ：
R k+1 ＝P e，k+1
Repeating the steps 1-3 until all the observed data are processed.
In the clock difference calculation module 300, the layout position of the first monitoring object is also optimized based on the improved genetic algorithm.
The improved genetic algorithm is specifically as follows:
the population is formed by combining M chromosomes, the M chromosomes contain n gene sequences, the gene sequences are subjected to crossover and mutation for iteration, a clock error function is an fitness function value, and the probability of being selected is higher as the fitness value is higher;
probability of selection p for M chromosomes m The following are provided;
wherein f m Is the reciprocal of fitness value of the m-th chromosome.
In the improved genetic algorithm, the crossover probability p of genes on M chromosomes a The method comprises the following steps:
wherein p is a Representing crossover probability of genes, p amax Represents the maximum crossover probability of a gene, p amin Representing the minimum crossover probability of a gene, f ai Indicating the fitness of the gene, f amin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the probability of variation of genes on the M chromosomes is:
wherein p is b Representing the probability of variation of the gene, p bmax Represents the maximum mutation probability of the gene, p bmin Representing the minimum probability of variation of the gene, f bi Indicating the fitness of the gene, f bmin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the fitness function value of the ith gene in the n gene sequences is f i The average fitness function value of the current gene isThe standard deviation ρ is defined as follows:
population standard deviation threshold value theta, maximum iteration times is R, times of standard deviation smaller than threshold value are C, C=C+1 when rho is smaller than or equal to theta, C is reset to 0 when rho is larger than theta, and algorithm is iterated toWhen the rho is less than or equal to theta in the continuous P iterations of the algorithm, gaussian disturbance is carried out on the current population data, and the variation probability P is increased b Let the current mutation probability p b Becomes 2p b 。
In this embodiment, the data sampling module 100 samples the base station position data, and optimizes and obtains the layout position of the best base station based on the improved genetic algorithm: the population is formed by combining 6 chromosomes, the 6 chromosomes contain 150 gene sequences, the gene sequences are subjected to crossover and mutation for iteration, a clock error function is set as an fitness function value, and the probability of being selected is higher as the fitness value is higher;
selection probability p of 6 chromosomes m The following are provided;
wherein f m Reciprocal of fitness value for the m-th chromosome;
crossover probability p of genes on 6 chromosomes a The method comprises the following steps:
wherein p is a Representing crossover probability of genes, p amax Represents the maximum crossover probability of a gene, p amin Representing the minimum crossover probability of a gene, f ai Indicating the fitness of the gene, f amin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the probability of variation of the genes on the 6 chromosomes is:
wherein p is b Representing the probability of variation of the gene, p bmax Represents the maximum mutation probability of the gene, p bmin Representing the minimum probability of variation of the gene, f bi Indicating the fitness of the gene, f bmin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the fitness function value of the ith gene in the 150 gene sequences is f i The average fitness function value of the current gene isThe standard deviation ρ is defined as follows:
the maximum iteration number is set to be 1000 times, the number of times of standard deviation is smaller than the threshold value is C, C=C+1 when ρ is smaller than or equal to θ, C is reset to 0 when ρ is larger than or equal to θ, when the algorithm iterates to 500 times, if ρ is smaller than or equal to θ in 50 continuous iterations of the algorithm, the current algorithm is judged to fall into a local optimal solution, and Gaussian disturbance is carried out on current population data and variation probability p is increased b Probability of current variation p b Becomes 2p b The algorithm search can be made to jump out locally optimal.
The data sampling module 100 samples the base station signal data of the optimal layout position and the signal data of the satellite 1 and the satellite 2, and based on time service and correction of the base station and the satellite 1 based on the crystal oscillator clock, the crystal oscillator clock is set to accurately time service the base station, and the crystal oscillator clock and the international standard time are in time synchronizationThe initial deviation and the initial deviation of the second period are a and b respectively, and the fixed offset error in each second is tau; the deviation deltat of the x second base station from the international standard time is counted x The method comprises the following steps:
Δt x ＝b+τx
deviation delta of x second of base station from international standard time 1 (x) The method comprises the following steps:
let satellite 1 time error obey normal distribution mu-N (0, sigma) 2 ) Wherein mu is a satellite 1 time random error, N is a normal distribution function, and sigma is a standard deviation;
the x second clock bias output by satellite 1 is:
δ 2 (x)＝x-μ x
wherein mu x Representing the clock error of the reference satellite of the x second;
obtaining a deviation model delta (x):
and carrying out deviation calculation based on the regression model to obtain a deviation value of the satellite 1 clock and the crystal oscillator clock.
The satellite 1 and the satellite 2 simultaneously receive signals of the base station, and an observation model is built as follows:
h 1 ＝L 1 +ct-cT 1 +d 1 +d 2
h 2 ＝L 2 +ct-cT 2 +d 1 +d 2
wherein h is 1 And h 2 Observed quantity, L, of satellite 1 and satellite 2, respectively 1 For the pseudoranges of the base station to satellite 1 and satellite 2, L 2 The pseudo range from the base station to the satellite 2 is represented by c, the signal propagation speed is represented by T, the clock error of the base station is represented by T 1 And T 2 Clock differences, d, of satellite 1 and satellite 2, respectively 1 And d 2 Respectively tropospheric delayDelay and ionospheric delay;
the differential observations Δh for satellite 1 and satellite 2 were obtained as follows:
Δh＝ΔL-cΔT
wherein DeltaL is a differential pseudo range, deltaT is a differential clock difference;
the satellite closest to the base station can be used as a reference satellite, and the clock differences of the other satellites can be obtained according to the differential clock differences.
The clock difference calculation module 300 also builds a clock difference matrix based on the clock difference data in a plurality of time periods, and the clock difference matrix is updated in real time according to the acquired and calculated clock difference data, and performs clock difference filtering based on a nonlinear Kalman filtering algorithm, and takes the clock difference matrix as a state variable X k+1 Outputting the optimal estimated clock difference:
X k+1 ＝f(X k )+V k
Z k+1 ＝H k+1 X k+1 +W k+1
wherein f (X) k ) As a state function, Z k+1 To measure the vector, H k+1 To observe the matrix, V k And W is k+1 System noise and measurement noise, respectively;
the state one-step prediction is:
the state one-step prediction mean square error is:
the filtering gain is:
the state estimation is:
the state estimation mean square error is:
P k+1，k+1 ＝(1-K k+1 H k+1 )P k+1，k
wherein,a clock difference prediction value indicating the next time; />Representing an optimal clock difference estimated value at the current moment; p (P) k+1，k A predicted mean square error matrix for the next moment; p (P) k，k Representing an estimated mean square error matrix at the current time; phi (phi) k+1，k Is a state transition matrix; j (J) k+1，k Is a Jacobian matrix; />Is a transpose of the jacobian matrix; q (Q) k Representing a process noise covariance; k (K) k+1 Representing the Kalman filtering gain; />Is the transposition of the observation matrix; r is R k+1 Is observed noise covariance; p (P) k+1，k+1 Representing an estimation error matrix; />The state clock difference optimal estimated value at the next time is represented.
The nonlinear Kalman filtering algorithm is used for iteration, so that the measurement error and the prediction error can be balanced, and an optimal clock error estimated value is output.
The clock correction module 400 performs clock correction on the satellite 2 time based on the optimal correction value obtained by the calculation of the clock calculation module 300 and the time service deviation.
It will be evident to those skilled in the art that the invention is not limited to the details of the foregoing illustrative embodiments, and that the present invention may be embodied in other specific forms without departing from the spirit or essential characteristics thereof. The present embodiments are, therefore, to be considered in all respects as illustrative and not restrictive, the scope of the invention being indicated by the appended claims rather than by the foregoing description, and all changes which come within the meaning and range of equivalency of the claims are therefore intended to be embraced therein. Any reference sign in a claim should not be construed as limiting the claim concerned.
Furthermore, it should be understood that although the present disclosure describes embodiments, not every embodiment is provided with a separate embodiment, and that this description is provided for clarity only, and that the disclosure is not limited to the embodiments described in detail below, and that the embodiments described in the examples may be combined as appropriate to form other embodiments that will be apparent to those skilled in the art.
Claims (7)
1. Satellite signal timing synchronization system, its characterized in that: comprising the following steps:
a data sampling module (100): the method comprises the steps of sampling data of a first monitoring object, a second monitoring object and a third monitoring object;
model building module (200): the observation model is used for respectively building the second monitoring object and the third monitoring object to observe the first monitoring object based on the crystal oscillator clock;
clock difference calculation module (300): the method comprises the steps of calculating clock difference based on a built observation model, and calculating and outputting optimal estimated clock difference through a nonlinear Kalman filtering algorithm;
clock correction module (400): the synchronization time of the first monitoring object is corrected according to the optimal estimated clock difference;
the clock difference calculation module (300) is further used for building a clock difference matrix based on the acquired clock difference data, performing clock difference filtering based on a nonlinear Kalman filtering algorithm and outputting an optimal estimated clock difference;
the clock difference matrix is updated in real time according to the acquired and calculated clock difference data;
in the nonlinear Kalman filtering algorithm, a clock difference matrix is used as a state variable X k+1 The method is characterized by comprising the following steps:
X k+1 ＝f(X k )+V k
Z k+1 ＝H k+1 X k+1 +W k+1
wherein f (X) k ) As a state function, Z k+1 To measure the vector, H k+1 To observe the matrix, V k And W is k+1 System noise and measurement noise, respectively;
the state one-step prediction is:
the state one-step prediction mean square error is:
the filtering gain is:
the state estimation is:
the state estimation mean square error is:
P k+1，k+1 ＝(1-K k+1 H k+1 )P k+1，k
wherein,representing the optimal clock difference predicted value of the next moment; />Representing an optimal clock difference estimated value at the current moment; p (P) k+1，k A predicted mean square error matrix for the next moment; p (P) k，k Representing an estimated mean square error matrix at the current time; phi (phi) k+1，k Is a state transition matrix; j (J) k+1，k Is a Jacobian matrix; />Is a transpose of the jacobian matrix; q (Q) k Representing a process noise covariance; k (K) k+1 Representing the Kalman filtering gain; />Is the transposition of the observation matrix; r is R k+1 Is observed noise covariance; p (P) k+1，k+1 Representing an estimation error matrix; />The state clock difference optimal estimated value at the next time is represented.
2. The satellite signal timing synchronization system of claim 1, wherein: the data sampling module (100) samples the first, second and third monitoring object data every interval preset time.
3. The satellite signal timing synchronization system of claim 1, wherein: the model building module (200) is used for carrying out time service on the first monitoring object and the second monitoring object based on the crystal oscillator clock, and specifically comprises the following steps:
setting a crystal oscillator clock to accurately time a first monitoring object, wherein the initial deviation of time synchronization of the crystal oscillator clock and the international standard time and the initial deviation of a second period are respectively a and b, and the fixed offset error in each second is tau; first monitor of x secondDeviation deltat of measured object from international standard time x The method comprises the following steps:
Δt x ＝b+τx
cumulative deviation delta between the x second of the first monitoring object and the international standard time 1 (x) The method comprises the following steps:
the time error of the second monitoring object obeys the normal distribution mu-N (0, sigma) 2 ) Wherein μ is a second monitoring object time random error, N is a normal distribution function, and σ is a standard deviation;
the x second clock deviation of the second monitor object output is:
δ 2 (x)＝x-μ x
wherein mu x A clock error of the second monitoring object representing the x second;
obtaining a deviation model delta (x):
performing deviation calculation based on the regression model to obtain a deviation value of the satellite clock and the crystal oscillator clock;
the second monitoring object and the third monitoring object observe the clock signal of the first monitoring object at the same time, and an observation model is built as follows:
h 1 ＝L 1 +ct-cT 1 +d 1 +d 2
h 2 ＝L 2 +ct-cT 2 +d 1 +d 2
wherein h is 1 And h 2 Observed amounts of the second monitoring object and the third monitoring object respectively, L 1 Pseudo-range from first monitoring object to second monitoring object, L 2 Pseudo range from the first monitoring object to the third monitoring object, c is signal propagation speed, T is clock difference of the first monitoring object, T 1 And T 2 Respectively are secondClock difference between the monitoring object and the third monitoring object, d 1 And d 2 Tropospheric delay and ionospheric delay, respectively.
4. A satellite signal timing synchronization system in accordance with claim 3, wherein: the clock difference calculation module (300) performs difference based on an observation model of the model building module (200) to obtain a difference observed quantity delta h of the second monitoring object and the third monitoring object as follows:
Δh＝ΔL-cΔT
wherein Δl is the differential pseudorange;
the method comprises the following steps:
wherein Δt is the differential clock difference.
5. The satellite signal timing synchronization system of claim 1, wherein: in the clock difference calculation module (300), the layout position of the first monitoring object is optimized and obtained based on the improved genetic algorithm.
6. The satellite signal timing synchronization system of claim 5, wherein: the improved genetic algorithm is specifically as follows:
the population is formed by combining M chromosomes, the M chromosomes contain n gene sequences, the gene sequences are subjected to crossover and mutation for iteration, a clock error function is an fitness function value, and the probability of being selected is higher as the fitness value is higher;
probability of selection p for M chromosomes m The following are provided;
wherein f m Is the reciprocal of fitness value of the m-th chromosome.
7. The satellite signal timing synchronization system of claim 6, wherein: in the improved genetic algorithm, the crossover probability p of genes on M chromosomes a The method comprises the following steps:
wherein p is a Representing crossover probability of genes, p amax Represents the maximum crossover probability of a gene, p amin Representing the minimum crossover probability of a gene, f ai Indicating the fitness of the gene, f amin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the probability of variation of genes on the M chromosomes is:
wherein p is b Representing the probability of variation of the gene, p bmax Represents the maximum mutation probability of the gene, p bmin Representing the minimum probability of variation of the gene, f bi Indicating the fitness of the gene, f bmin Represents the minimum value of the fitness of the gene,respectively representing the fitness average value of genes, wherein e is a mathematical constant;
the fitness function value of the ith gene in the n gene sequences is f i The average fitness function value of the current gene isThe standard deviation ρ is defined as follows:
population standard deviation threshold value theta, maximum iteration times is R, times of standard deviation smaller than threshold value are C, C=C+1 when rho is smaller than or equal to theta, C is reset to 0 when rho is larger than theta, and algorithm is iterated toWhen the rho is less than or equal to theta in the continuous P iterations of the algorithm, gaussian disturbance is carried out on the current population data, and the variation probability P is increased b Let the current mutation probability p b Becomes 2p b 。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310966519.4A CN116847452B (en) | 2023-08-03 | 2023-08-03 | Satellite signal timing synchronization system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310966519.4A CN116847452B (en) | 2023-08-03 | 2023-08-03 | Satellite signal timing synchronization system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116847452A CN116847452A (en) | 2023-10-03 |
| CN116847452B true CN116847452B (en) | 2024-03-15 |
Family
ID=88167238
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310966519.4A Active CN116847452B (en) | 2023-08-03 | 2023-08-03 | Satellite signal timing synchronization system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116847452B (en) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103957095A (en) * | 2014-05-15 | 2014-07-30 | 北京航空航天大学 | Time synchronization method based on intersatellite bidirectional distance measurement |
| CN105071914A (en) * | 2015-07-27 | 2015-11-18 | 深圳思凯微电子有限公司 | Digital frequency modulation (FM) broadcast based clock synchronization method and FM broadcast receiver |
| CN105610534A (en) * | 2015-12-17 | 2016-05-25 | 北京无线电计量测试研究所 | Multi-station time synchronization method and equipment |
| CN109299496A (en) * | 2018-07-26 | 2019-02-01 | 国网湖南省电力有限公司 | A kind of high-precise synchronization clock generation method |
| CN111948686A (en) * | 2020-08-05 | 2020-11-17 | 航天恒星科技有限公司 | Time synchronization method and device |
| CN114565122A (en) * | 2021-12-29 | 2022-05-31 | 重庆驰知科技有限公司 | UWB-based satellite-signal-free area positioning system |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070218931A1 (en) * | 2006-03-20 | 2007-09-20 | Harris Corporation | Time/frequency recovery of a communication signal in a multi-beam configuration using a kinematic-based kalman filter and providing a pseudo-ranging feature |
-
2023
- 2023-08-03 CN CN202310966519.4A patent/CN116847452B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103957095A (en) * | 2014-05-15 | 2014-07-30 | 北京航空航天大学 | Time synchronization method based on intersatellite bidirectional distance measurement |
| CN105071914A (en) * | 2015-07-27 | 2015-11-18 | 深圳思凯微电子有限公司 | Digital frequency modulation (FM) broadcast based clock synchronization method and FM broadcast receiver |
| CN105610534A (en) * | 2015-12-17 | 2016-05-25 | 北京无线电计量测试研究所 | Multi-station time synchronization method and equipment |
| CN109299496A (en) * | 2018-07-26 | 2019-02-01 | 国网湖南省电力有限公司 | A kind of high-precise synchronization clock generation method |
| CN111948686A (en) * | 2020-08-05 | 2020-11-17 | 航天恒星科技有限公司 | Time synchronization method and device |
| CN114565122A (en) * | 2021-12-29 | 2022-05-31 | 重庆驰知科技有限公司 | UWB-based satellite-signal-free area positioning system |
Non-Patent Citations (4)
| Title |
|---|
| 刘亮等.基于Kalman滤波的持续卫星时间同步攻击防护方法.电力系统自动化.2020,(第10期),全文. * |
| 北斗广域高精度时间服务原型系统;施闯等;测绘学报;全文 * |
| 基于Kalman滤波的持续卫星时间同步攻击防护方法;刘亮等;电力系统自动化(第10期);全文 * |
| 施闯 ; 张东 ; 宋伟 ; 于佳亮 ; 郭文飞 ; .北斗广域高精度时间服务原型系统.测绘学报.2020,(第03期),全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116847452A (en) | 2023-10-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP3430429B1 (en) | Satellite navigation receiver with improved ambiguity resolution | |
| US7741996B2 (en) | Positioning system, information supply device, terminal device, control method of information supply device, control program of information supply device, and computer readable recording medium recording control program of information supply device | |
| CN104714244A (en) | Multi-system dynamic PPP resolving method based on robust self-adaption Kalman smoothing | |
| CN107678050A (en) | GLONASS phase inter-frequency deviation real-time tracings and precise estimation method based on particle filter | |
| CN111736185B (en) | Terminal positioning method and device, computer readable storage medium and terminal equipment | |
| CN107966722B (en) | GNSS clock error resolving method | |
| AU2008238113A1 (en) | Method of processing radionavigation signals | |
| WO2017160585A1 (en) | Satellite navigation receiver with improved ambiguity resolution | |
| Wang et al. | Atomic clock performance and combined clock error prediction for the new generation of BeiDou navigation satellites | |
| CN110940976B (en) | Error correction-based multi-station multi-external radiation source radar moving target positioning method | |
| CN108196267A (en) | A kind of uninterrupted time delivering method based on GNSS CP technologies | |
| CN106772510B (en) | A kind of frequency hopping distance measuring method based on carrier phase measurement | |
| CN112433235A (en) | Method, system, and medium for determining a time reference | |
| Ruan et al. | Between-satellite single-difference integer ambiguity resolution in GPS/GNSS network solutions | |
| CN116847452B (en) | Satellite signal timing synchronization system | |
| CN116125512A (en) | PPP self-adaptive clock difference model estimation method considering clock frequency time-varying characteristics | |
| CN116594046A (en) | Moving target positioning method based on low orbit satellite signal Doppler error compensation | |
| CN116224746B (en) | High-stability time reference establishment method for satellite-ground atomic clock fusion | |
| Wang et al. | Gaussian filtered rssi-based indoor localization in wlan using bootstrap filter | |
| CN112880659A (en) | Fusion positioning method based on information probability | |
| Minetto | A theoretical study on the benefits of integrating GNSS and collaborative relative ranges | |
| CN115267847B (en) | Single-frequency cycle slip detection and repair method, system, equipment and medium | |
| CN118277879A (en) | Satellite clock error forecasting method based on fusion prediction model | |
| CN116794688A (en) | Robust filtering method based on satellite quality evaluation | |
| Glushankov et al. | Methods to Improve the Accuracy of State Estimation in Nonlinear Kalman Filter in Case of a Priori Uncertainty |
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 |