CN102998681B - A kind of high-frequency clock error estimation method of satellite navigation system - Google Patents

A kind of high-frequency clock error estimation method of satellite navigation system Download PDF

Info

Publication number
CN102998681B
CN102998681B CN201210541485.6A CN201210541485A CN102998681B CN 102998681 B CN102998681 B CN 102998681B CN 201210541485 A CN201210541485 A CN 201210541485A CN 102998681 B CN102998681 B CN 102998681B
Authority
CN
China
Prior art keywords
mrow
msub
clock
msubsup
satellite
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
Application number
CN201210541485.6A
Other languages
Chinese (zh)
Other versions
CN102998681A (en
Inventor
陈智
张小强
张雪辉
宋杨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CHINA AEROSPACE SCIENCE & INDUSTRY ACADEMY OF INFORMATION TECHNOLOGY
Original Assignee
CHINA AEROSPACE SCIENCE & INDUSTRY ACADEMY OF INFORMATION TECHNOLOGY
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by CHINA AEROSPACE SCIENCE & INDUSTRY ACADEMY OF INFORMATION TECHNOLOGY filed Critical CHINA AEROSPACE SCIENCE & INDUSTRY ACADEMY OF INFORMATION TECHNOLOGY
Priority to CN201210541485.6A priority Critical patent/CN102998681B/en
Publication of CN102998681A publication Critical patent/CN102998681A/en
Application granted granted Critical
Publication of CN102998681B publication Critical patent/CN102998681B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention discloses a kind of high-frequency clock error estimation method of satellite navigation system, comprising: the observation data obtaining tracking station; The reliability of check observation data; Observation data through certificate authenticity is utilized to the linear combination of observed quantity, reject outlier, determine initial blur level, detect and repair cycle slip, after smoothing processing, obtain code and carrier phase observable; Obtain satellite orbit parameter, polar motion of globe parameter, troposphere parameter, survey station coordinate parameters; According to orbit parameter, polar motion of globe parameter, troposphere parameter, the survey station coordinate parameters of original observed data and satellite; Obtain the clock bias estimation value of satellite and receiver; Utilize code and carrier phase observable, interpolation is carried out to the clock bias estimation obtained, obtain the high frequency clock bias estimation of satellite navigation system.The present invention can improve the estimated accuracy of satellite navigation system high frequency clock correction, has that reliability is high, precision advantages of higher.

Description

High-frequency clock error estimation method of satellite navigation system
Technical Field
The invention relates to the technical field of satellite navigation positioning, in particular to a high-frequency clock error estimation method of a satellite navigation system.
Background
The satellite clock error refers to the difference between the satellite clock time and the standard time of the navigation system caused by the frequency drift of the satellite clock. The observation quantity of the navigation system is obtained by taking clock frequency signals of a satellite and a receiver as reference, the positioning precision of a user is closely related to the precision of the clock error and the orbit of the satellite, which can be obtained by the user, and the estimation precision of the clock error is closely related to the precision of the orbit estimation (also called a coupling relation).
Taking a GPS (Global Positioning System) navigation satellite as an example, The accuracy of a GPS satellite clock difference currently used as a broadcast ephemeris for transmitting a navigation message is only 7ns, which cannot meet The requirements of a decimeter-level Positioning user, and The fast and ultra-fast GPS satellite clock differences currently provided by an IGS (The International GNSS Service, International GPS Service organization) center are only 15min sampling intervals, while The 30s GPS satellite clock difference provided by a CODE center is delayed by 14-17 hours, which sometimes hardly meets The engineering requirements. On the other hand, if a high-frequency and high-precision satellite clock error product is to be obtained, theoretically, the more observation data, the higher the precision of the resolved satellite clock error, but if the ground observation data are too much, the more corresponding parameters are resolved, and further, a higher requirement is put forward on the running speed of the computer. In the most important aspect, if the IGS center is unable to offer some GPS satellite products for some natural or political reasons, it will certainly affect precise orbit determination for satellite-borne GPS low-orbit satellites.
Disclosure of Invention
The invention aims to solve the technical problem that a high-frequency clock error estimation method of a satellite navigation system can effectively solve the problem that the satellite clock error with high precision of 30s or higher frequency is difficult to obtain.
In order to solve the above technical problem, the present invention provides a high frequency clock error estimation method for a satellite navigation system, comprising:
s101, acquiring observation data of a tracking station;
s102, checking the reliability of the observation data;
s103, eliminating outliers from observation data subjected to reliability test by utilizing linear combination of observation quantities, determining initial ambiguity, detecting and repairing cycle slip, and obtaining a code and phase observation value after smoothing;
s104, acquiring satellite orbit parameters, earth polar motion parameters, troposphere parameters and survey station coordinate parameters;
s105, the original observation data obtained in the step S103 and the orbit parameters, the earth polar motion parameters, the troposphere parameters and the survey station coordinate parameters of the satellite obtained in the step S104 are obtained; bringing the data into a non-differential pseudo range and phase observation value error equation with ionosphere influence eliminated, selecting a reference clock to apply constraint on a clock to be estimated, and performing unified adjustment processing to obtain clock error estimation values of a satellite and a receiver;
and S106, interpolating the clock error estimation obtained in the step S105 by using the code and phase observed value in the step S103 to obtain a high-frequency clock error estimation of the satellite navigation system.
Further, non-differential phase and pseudo-range observed values influenced by an ionized layer are eliminated, and an observed value error equation is as follows:
<math><mrow> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>&Phi;</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&Delta;t</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&Delta;t</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&delta;</mi> <msubsup> <mi>&rho;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&lambda;</mi> <mo>&CenterDot;</mo> <msubsup> <mi>N</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>/</mo> <mi>C</mi> <mo>+</mo> <msubsup> <mi>&epsiv;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>&Phi;</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>&lambda;</mi> <mo>&CenterDot;</mo> <msubsup> <mi>&Phi;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow></math>
<math><mrow> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>p</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&Delta;t</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&Delta;t</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&delta;</mi> <msubsup> <mi>&rho;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <msubsup> <mi>&epsiv;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>p</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>P</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow></math>
in formulas (4) and (5): k is the station number, j is the satellite number, i is the corresponding observation epoch, and C is the speed of light in vacuum; Δ tk(i) For receiver clock error, Δ tj(i) In order to be the clock error of the satellite,is tropospheric delay effects;for unmodeled error effects;is an ambiguity parameter;the pseudo range and phase combination observed value without the influence of an ionized layer is obtained;lambda is the wavelength of the ionosphere-free combined observed value as the observation error;the geometric distance between the satellite position at the time of signal transmission to the receiver position at the time of signal reception.
Further, the pseudorange to clock error and phase error equations may be expressed as:
<math><mrow> <mi>V</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>V</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mi>&Phi;</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>A</mi> <mi>p</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>A</mi> <mi>&Phi;</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>L</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>L</mi> <mi>&Phi;</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow></math>
in formula (6): <math><mrow> <msub> <mi>A</mi> <mi>p</mi> </msub> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <munder> <mi>A</mi> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> <mtd> <munder> <mn>0</mn> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow></math> <math><mrow> <msub> <mi>A</mi> <mi>&Phi;</mi> </msub> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <munder> <mi>A</mi> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> <mtd> <munder> <mrow> <mi>&lambda;</mi> <mo>/</mo> <mi>C</mi> <mo>&CenterDot;</mo> <mi>I</mi> </mrow> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow></math> let I be an n-order identity matrix, A be a design matrix for determining clock offset of a receiver and a satellite, be an (m.n) (m + n) matrix, V (I) be an observation error of an ith observation position, L (I) be a free term of the ith observation position, phi and p respectively represent phase and pseudo-range observation, X (I) be an ith parameter to be solved, V (I), Lp(i)、LΦ(i) Is a known amount; the number of m reference stations and n is the number of synchronously tracking satellites.
Furthermore, one receiver clock in the tracking network is used as a reference clock, and the clock difference precision of the reference clock is better than 10-6s; the relative clock error of all satellites in the same navigation system is zero mean.
Further, the reference clock is a hydrogen atomic clock.
Further, step S106 specifically includes:
s1061, clock estimation based on the code observation;
s1062, clock estimation based on the phase observation;
and S1063, performing clock combination to obtain high-frequency clock difference estimation of the navigation satellite and the receiver.
The invention has the following beneficial effects:
the method can improve the estimation precision of the high-frequency clock error of the satellite navigation system, and has the advantages of high reliability, high precision and the like.
Drawings
Fig. 1 is a flowchart of a high frequency clock error estimation method of a satellite navigation system according to an embodiment of the present invention.
Detailed Description
The present invention will be described in further detail below with reference to the drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not limit the invention.
The estimation of the satellite clock error can have two ideas according to different parameter processing modes: firstly, putting clock parameters and other parameters such as satellite orbit parameters, survey station coordinates, troposphere parameters and the like into a block, and estimating simultaneously; the other method is to firstly solve the non-epoch parameters such as the satellite orbit parameter, the coordinate parameter of the measuring station and the like, solve and fix the non-epoch parameters, then solve the epoch parameters, reduce the number of the parameters by the epoch, and finally solve the clock error parameter. For the two methods, if the epoch parameter is less, the first parameter can be adopted, and the second method is used for solving the parameter step by step and in stages, the method can be used for solving the multi-epoch parameter, and the method is more effective in efficiency, so that the second method is adopted in the embodiment of the invention when the clock difference parameter is solved.
Since the GPS observations are the relative time delays between the survey station and the satellites. Thus, all satellite and receiver clock differences cannot be determined simultaneously, and must be determined by fixing the clock difference of a reference clock (either the receiver clock or the satellite clock), then determining the relative clock differences of other receivers and satellites,but the clock difference precision of the reference clock must be ensured to be better than 10-6And s. The relative clock difference and the absolute clock difference are equivalent to the positioning result of the user, that is, the systematic deviation of the relative clock difference can be completely absorbed by the clock difference of the user receiver in the positioning model of the user without influencing the positioning accuracy of the user.
As shown in fig. 1, an embodiment of the present invention relates to a high-precision, high-reliability and high-frequency satellite clock error estimation method for a satellite navigation system, including:
step S101, acquiring observation data of a tracking station;
and downloading and acquiring raw observation data provided by the ground tracking station, wherein the raw observation data comprises RINEX observation data and broadcast ephemeris of the satellite.
Step S102, checking the reliability of the observation data;
the reliability indicator being adapted to verify both the code observations and the phase observations, i.e.
<math><mrow> <mi>&gamma;</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>v</mi> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>v</mi> </munderover> <msub> <mrow> <mo>(</mo> <mi>PR</mi> <mo>)</mo> </mrow> <mi>ij</mi> </msub> <mo>/</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>v</mi> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>v</mi> </munderover> <msub> <mi>P</mi> <mi>ij</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow></math>
In the formula (1), gamma is a reliability measure of the system at the ith observation position; for the code observed value, v is 1, the larger gamma is, the better the reliability is, and the smaller gamma is, the worse the reliability is; p is an observation weight matrix, and R is an orthogonal projection matrix. If the component of the redundant observation is approximately equal to zero, the reliability of the ith observation position is poor, and the reliability of the system at the ith observation position is poor when the redundant observation component is smaller; j is the satellite number. The internal reliability index MDB and the external reliability MDE can be expressed as
<math><mrow> <msub> <mrow> <mo>|</mo> <mi>MDB</mi> <mo>|</mo> </mrow> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <msubsup> <mi>&lambda;</mi> <mn>0</mn> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </msubsup> <msub> <mi>&sigma;</mi> <mn>0</mn> </msub> <msup> <mrow> <mo>[</mo> <mi>&gamma;</mi> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>v</mi> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>v</mi> </munderover> <msub> <mi>P</mi> <mi>ij</mi> </msub> <mo>]</mo> </mrow> <mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow></math>
MDEc(p)0-1-1) (3)
In the formulas (2) and (3), c (p) represents a code or a phase observed value; lambda [ alpha ]0It is a non-central parameter under a certain probability condition, and is a function of the false rejection probability, the test efficacy and the number of model parameters (gross errors) to be detected. Once these three parameters are selected, λ0Is a constant; usually, these three parameters are taken to be 0.001, 0.80 and 1, then λ0Was 17.07. When comparing the reliability of different observation positions of the same system, only the relative ratio of the two is concerned, and the constant lambda is0And does not work. On the other hand, the false rejection probability and the test efficacy parameter are generally selected according to empirical information, and have certain subjectivity and arbitrariness. Sigma0For mid-error of observations, v is the number of consecutive observations for a quantity. P, R is as defined above.
Step S103, preprocessing data to obtain a code and phase observed value;
and (4) if the reliability of the data subjected to reliability inspection in the step (S102) is very poor, removing, and selecting a threshold value according to actual observation experience. And (4) further utilizing the linear combination of the observed quantities to eliminate outliers aiming at the data obtained in the step S102, determining initial ambiguity, detecting and repairing cycle slip, and obtaining a code and phase observed value after smoothing. And estimating the receiver clock error of each epoch by using least squares according to the pseudo-range observed value and the satellite clock error in the broadcast ephemeris.
Step S104, calculating parameters such as orbit, earth rotation and the like of the satellite;
and forming differential information by using the codes and the phase observation values preprocessed in the step S103, and solving satellite orbit parameters, earth polar motion parameters, troposphere parameters, survey station coordinate parameters and the like. It is of course also possible to obtain the forecast part of the ultrafast track product from the IGS center, when the track parameters may not have to be solved, directly to the next step.
Step S105, estimating the satellite clock error with low sampling rate;
and taking the original observation data obtained by the processing of the step S103 and the orbit parameters, the earth polar motion parameters, the troposphere parameters and the coordinate parameters of the survey station of the satellite obtained by the step S104 as known values, bringing the known values into a non-differential pseudo range and a phase observation value error equation which eliminate the influence of an ionosphere, selecting a reference clock so as to apply constraint on the clock to be estimated, and carrying out unified adjustment processing to obtain the clock error estimation values of the satellite and the receiver with high precision.
In the precision satellite clock error estimation, non-differential phase and pseudo-range observed values with ionosphere influence eliminated are generally adopted, and an observed value error equation is as follows:
<math><mrow> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>&Phi;</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&Delta;t</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&Delta;t</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&delta;</mi> <msubsup> <mi>&rho;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&lambda;</mi> <mo>&CenterDot;</mo> <msubsup> <mi>N</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>/</mo> <mi>C</mi> <mo>+</mo> <msubsup> <mi>&epsiv;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>&Phi;</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>&lambda;</mi> <mo>&CenterDot;</mo> <msubsup> <mi>&Phi;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow></math>
<math><mrow> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>p</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&Delta;t</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&Delta;t</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&delta;</mi> <msubsup> <mi>&rho;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <msubsup> <mi>&epsiv;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>p</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>P</mi> <mi>k</mi> <mi>i</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow></math>
in formulas (4) and (5): k is the station number, j is the satellite number, i is the corresponding observation epoch, and C is the speed of light in vacuum;for receiver clock error, Δ tj(i) In order to be the clock error of the satellite,is tropospheric delay effects;for unmodeled error effects;is an ambiguity parameter;the pseudo range and phase combination observed value without the influence of an ionized layer is obtained;lambda is the wavelength of the ionosphere-free combined observed value as the observation error;the geometric distance between the satellite position at the time of signal transmission to the receiver position at the time of signal reception.
Setting m reference stations, synchronously tracking n satellites, and expressing the pseudo range and phase error equation of clock error as follows:
<math><mrow> <mi>V</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>V</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mi>&Phi;</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>A</mi> <mi>p</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>A</mi> <mi>&Phi;</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>L</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>L</mi> <mi>&Phi;</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow></math>
in formula (6):
<math><mrow> <msub> <mi>A</mi> <mi>p</mi> </msub> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <munder> <mi>A</mi> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> <mtd> <munder> <mn>0</mn> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow></math> <math><mrow> <msub> <mi>A</mi> <mi>&Phi;</mi> </msub> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <munder> <mi>A</mi> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> <mtd> <munder> <mrow> <mi>&lambda;</mi> <mo>/</mo> <mi>C</mi> <mo>&CenterDot;</mo> <mi>I</mi> </mrow> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </munder> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow></math> let I be an n-order identity matrix, A be a design matrix for determining clock offset of a receiver and a satellite, be an (m.n) (m + n) matrix, V (I) be an observation error of an ith observation position, L (I) be a free term of the ith observation position, phi and p respectively represent phase and pseudo-range observation, X (I) be an ith parameter to be solved, V (I), Lp(i)、LΦ(i) In known amounts. If the above equation is singular for solving the satellite clock error parameters for the observation equation, the normal equation must introduce a reference clock to solve the clock error parameters, and the clock errors of the other receiver clocks and the satellite clock with respect to the reference clock must be solved. When a reference clock is introduced, the invention adopts two measures aiming at a ground receiver clock and a satellite clock: (1) one receiver clock in the trace network is used as a reference clock, but the premise is to ensure that the reference clock is stable and reliable (such as a hydrogen atomic clock). (2) Assume that the relative clock bias of all satellites in the same navigation system is a zero-mean condition. The equations then solve for the relative clock difference of the satellites.
Step S106, acquiring a high-precision and high-frequency satellite clock error;
the low frequency satellite clock offset is interpolated by reusing the code and phase observations from step S103. The conventional interpolation method is relatively poor in interpolation accuracy due to the influence of the data time series. In order to obtain higher interpolation precision, the step relates to a special satellite clock difference interpolation method, and the algorithm is a high-precision and high-efficiency clock difference interpolation method; has the following characteristics: 1. the method is used for satellite clock error interpolation, so that the clock error result of each satellite is more stable, the jump of very large even a few ns can not occur, and the interpolation precision can be improved; 2. the accuracy of a few satellites after clock difference interpolation is slightly reduced, but the magnitude of the influence is small, and the accuracy of the whole satellite interpolation method is higher than that of the general interpolation method; 3. with the increase of the number of stations, the trend of the clock interpolation result tends to be smooth, the corresponding precision is improved to a certain extent, and unstable and jumping with larger magnitude cannot occur.
Taking a GPS satellite as an example, in order to ensure the effectiveness of the algorithm, it is necessary to use a GPS satellite orbit, an earth polar motion parameter, a survey station coordinate, a troposphere parameter, and a low-sampling-rate satellite clock error, which are already obtained in the double-difference processing of the precise orbit determination of the GPS satellite in step S104 and in the non-difference data processing of the low-frequency clock error in step S105, and according to these information, the interpolation method of this step is used to correct the clock error at a high frequency, and the algorithm can perform high-frequency clock error correction interpolation at sampling intervals of 30S and 5S.
For describing the interpolation method, the present invention describes the algorithm by taking the interpolation of 30s clock difference in 300s clock difference as an example. The parameter estimation is typically based on a least squares method. The interpolation method comprises the following steps: 1. performing clock estimation based on the code observations; 2. performing clock error estimation based on the phase difference observation value; 3. the reference clocks are selected and combined. It should be emphasized here that if a clock difference with a high sampling rate is to be obtained, the raw observation data must have the same sampling frequency as the high-frequency clock difference that was finally obtained. In step S105, in order to avoid excessive solution parameters, we perform thinning processing on the original observed data, so that the obtained clock difference frequency is the same as the frequency of the thinned original data (for example, the calculated original observed sampling rate is 300S, and the calculated clock difference sampling rate is also 300S). Therefore, when solving for the high frequency clock error, the sampling frequency of the original data must be increased to obtain the satellite clock error with a sampling interval of 30s or more.
Step S1061, clock estimation based on the code observation value;
the receiver clock error must be synchronized with the GPS system time in order to ensure that the geometric distance error between the satellite and the receiver due to the receiver clock error is reduced to less than 1mm, and the receiver clock error correction and the satellite clock error correction are separately estimated for each epoch when synchronization is performed using code observations.
The observation equation of the GPS code is constructed by using the observation data (data sampling rate 30S) processed in step S103, and expressed as follows (unit: m, simplified equation, and no error term):
<math><mrow> <msubsup> <mi>p</mi> <mi>fk</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msup> <mi>&delta;</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msub> <mi>&delta;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mrow> <mi>fk</mi> <mo>,</mo> <mi>iono</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mrow> <mi>fk</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msup> <mi>b</mi> <mi>j</mi> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow></math>
wherein,is the observed code pseudorange;is the pitch;jis satellite clock error correction;kreceiver clock error correction;is ionospheric refraction;is the tropospheric reflection; bjIs a Differential Code Bias (DCB); c is the speed of light; the index k denotes the receiver, f denotes the frequency, the index j denotes the satellite, tiRefers to the observation epoch.
Using the satellite orbit parameters, the earth polar motion parameters, the survey station coordinate parameters, the troposphere parameters and the DCB obtained in step S104 as known numbers, the observation equation is:
<math><mrow> <msubsup> <mi>p</mi> <mi>fk</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msup> <mi>&delta;</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msub> <mi>&delta;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mrow> <mi>fk</mi> <mo>,</mo> <mi>iono</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow></math>
if constructedAnddeionizing ionosphere LC combinationThen only satellite and receiver clock error corrections remain:
<math><mrow> <msubsup> <mi>Cp</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msup> <mi>&delta;</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msub> <mi>&delta;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow></math>
based on such observation equations, a Normal Equation (NEQ) is established for each epoch for all relevant satellites and receivers, respectively. Meanwhile, the step of screening is included in the clock error estimation, and because a parameter is estimated by using a large amount of data, the receiver and the satellite clock error correction estimation algorithm are very stable, and unqualified data can be marked through iteration according to the size of a residual error and eliminated in subsequent iteration.
Step S1062, clock estimation based on the phase observation value;
for the generation of high frequency clock error correction, the phase observation is the main observation, and in contrast to the code observation process, the GPS phase observation equation is constructed using the observation data (code and phase observation, data sampling rate 30S) processed in step S103 as expressed below (unit: m, simplified formula, and no error term):
<math><mrow> <msubsup> <mi>&phi;</mi> <mi>fk</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msup> <mi>&delta;</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mo>&CenterDot;</mo> <msub> <mi>&delta;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>&rho;</mi> <mrow> <mi>fk</mi> <mo>,</mo> <mi>iono</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mrow> <mi>fk</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&lambda;</mi> <mi>f</mi> </msub> <mo>&CenterDot;</mo> <msubsup> <mi>N</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow></math>
wherein,is a phase observation; lambda [ alpha ]fIs the wavelength of frequency f;is the degree of ambiguity; compared with the code observation equation, the phase observation equation has the following differences: term of degree of ambiguityAnd code bias information (DCB) and ionospheric refraction is of opposite sign. To eliminate the ambiguity term in the phase observation equation, a phase observation operation is performed on subsequent epochs:
<math><mrow> <msubsup> <mi>&Delta;&phi;</mi> <mi>fk</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&phi;</mi> <mi>fk</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>&phi;</mi> <mi>fk</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow></math>
further, the same procedure as in the code observation processing was used to obtain a deionizing layer LC combinationThe phase observation equation of (1) is as follows:
<math><mrow> <mi>&Delta;</mi> <msubsup> <mi>C&phi;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <mi>C</mi> <mo>&CenterDot;</mo> <mi>&Delta;</mi> <msup> <mi>&delta;</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <msub> <mrow> <mo>&CenterDot;</mo> <mi>&Delta;&delta;</mi> </mrow> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow></math>
the difference result of the satellite and receiver clock difference correction values estimated with the single epoch algorithm is used as the prior value of the phase difference process. Thus, only the correction values of the a priori clock difference are estimated. If the code and the phase observed value are inconsistent, the corresponding phase observed value is removed, which is equivalent to the step S1061 of constraining the result calculated in the step S1062; while the equation is adjusted using the zero mean condition of the satellite clock correction values involved in the calculation. From the variance-covariance matrix of each epoch difference, the variance of each satellite j is obtainedAnd for each receiver kThese values are used in subsequent calculations, subject to high frequency clock error corrections of the satellite and the receiver.
Step S1063, performing clock combination to obtain high-frequency clock difference estimation of the navigation satellite and the receiver;
in order to obtain a correction value of the clock difference, clock difference combination is performed for each satellite and each receiver according to the code observation value obtained in step S1061 or the estimated values of the clock difference of the satellite and the receiver obtained from the high-precision clock difference time series (such as 5min clock difference solution) in step S105; in order to obtain the clock correction value of the "absolute" phase clock difference, interpolation for fixed clock correction is performed using this clock difference information.
Using the difference in phase clock as an observed value in step S1062i =1, L, n-1 (n is the number of epochs) establishes a new observation equation. Correcting value by known clock differenceAs a pseudo observed value, it is included in the clock difference set of the 5min solution to form a new observation equation set, as follows:
wherein (t)i) For the phase clock difference of the i-th observation,andthe clock difference term of the previous term is subtracted from the phase clock difference term,for fixed clock error correction, i is more than or equal to 1 and less than or equal to n. Neglecting the correlation between the epoch difference decompositions and having dimensions d = (n-1) + nfix×(n-1)+nfixIs a diagonal matrix (n)fixIs a clock errorThe number of epochs where the solution interval is relatively low). Obtaining a pseudo-observed value from the posterior RMS value of the last epoch Difference decompositionThe right of (1).
<math><mrow> <mi>P&Delta;&delta;</mi> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>P</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mfrac> <msubsup> <mi>&sigma;</mi> <mn>0</mn> <mn>2</mn> </msubsup> <msubsup> <mi>&sigma;</mi> <mrow> <mi>&Delta;&delta;</mi> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> <mn>2</mn> </msubsup> </mfrac> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>.</mo> <mo>.</mo> <mo>.</mo> <mo>,</mo> <mi>n</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow></math>
Wherein σ0Pseudorange differencing for ith observationCorresponding error in unit weight, P is weight matrix, has been fixed forCorrecting the known clock error value by applying a pseudo-observed value to the known clock error value according to the weight matrix informationAnd (6) carrying out constraint.
Finally, based on equation (13) and constraints (14), an effective solution to the system of equations, i.e., the satellite and receiver clock differences for the high reliability, high accuracy 30s or higher sampling intervals to be found, is found.
Table 1 shows the estimated GPS satellite clock difference between three days and 30s sampling intervals calculated by using 40 IGS reference stations from 197 th to 199 th of 2012, and compares the calculated result with the CODE center 30s clock difference, and as can be seen from table 1, the average value of the satellite clock difference estimated by three days of 31 GPS satellites is about 0.55ns, and the RMS is within 0.3ns, so the test result shows that the accuracy can reach about 0.2ns when the satellite clock difference processed by using the method is compared with the CODE center 30s clock difference.
TABLE 1 statistical table of GPS satellite clock error estimation results
The embodiment shows that the method can improve the estimation precision of the high-frequency clock error of the satellite navigation system, and has the advantages of high reliability, high precision and the like.
Although the preferred embodiments of the present invention have been disclosed for illustrative purposes, those skilled in the art will appreciate that various modifications, additions and substitutions are possible, and the scope of the invention should not be limited to the embodiments described above.

Claims (5)

1. A high frequency clock error estimation method of a satellite navigation system is characterized by comprising the following steps:
s101, acquiring observation data of a tracking station;
s102, checking the reliability of the observation data;
s103, eliminating outliers from observation data subjected to reliability test by utilizing linear combination of observation quantities, determining initial ambiguity, detecting and repairing cycle slip, and obtaining a code and phase observation value after smoothing;
s104, acquiring satellite orbit parameters, earth polar motion parameters, troposphere parameters and survey station coordinate parameters;
s105, the original observation data obtained in the step S103 and the orbit parameters, the earth polar motion parameters, the troposphere parameters and the survey station coordinate parameters of the satellite obtained in the step S104 are obtained; bringing the data into a non-differential pseudo range and phase observation value error equation with ionosphere influence eliminated, selecting a reference clock to apply constraint on a clock to be estimated, and performing unified adjustment processing to obtain clock error estimation values of a satellite and a receiver;
s106, interpolating the clock error estimation obtained in the step S105 by using the code and phase observation value in the step S103 to obtain a high-frequency clock error estimation of the satellite navigation system;
step S106 specifically includes:
s1061, clock estimation based on the code observation;
s1062, clock estimation based on the phase observation;
and S1063, performing clock combination to obtain high-frequency clock difference estimation of the navigation satellite and the receiver.
Performing clock difference combination on each satellite and the receiver according to the code observation value obtained in the step S1061 or the clock difference estimation values of the satellite and the receiver obtained from the high-precision clock difference time sequence in the step S105; in order to obtain the clock correction value of the "absolute" phase clock difference, interpolation for fixed clock correction is performed using this clock difference information.
2. The method of high frequency clock bias estimation for a satellite navigation system of claim 1, wherein the ionospheric effects of the non-differential phase and pseudorange observations are removed, and wherein the observation error equation is as follows:
<math> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>&Phi;</mi> </mrow> <mi>j</mi> </msubsup> <mi></mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&Delta;t</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&Delta;t</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&delta;</mi> <msubsup> <mi>&rho;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&lambda;</mi> <mo>&CenterDot;</mo> <msubsup> <mi>N</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mo>/</mo> <mi>C</mi> <mo>+</mo> <msubsup> <mi>&epsiv;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>&Phi;</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>&lambda;</mi> <mo>&CenterDot;</mo> <msubsup> <mi>&Phi;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>v</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>p</mi> </mrow> <mi>j</mi> </msubsup> <mi></mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&Delta;t</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>&Delta;t</mi> <mi>j</mi> </msup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&rho;</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <mi>&delta;</mi> <msubsup> <mi>&rho;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>trop</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>+</mo> <msubsup> <mi>&epsiv;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mi>p</mi> </mrow> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>P</mi> <mi>k</mi> <mi>j</mi> </msubsup> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>C</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>
in formulas (4) and (5): k is the station number, j is the satellite number, i is the corresponding observation epoch, and C is the speed of light in vacuum; Δ tk(i) For receiver clock error, Δ tj(i) In order to be the clock error of the satellite,is tropospheric delay effects;for unmodeled error effects;is an ambiguity parameter;the pseudo range and phase combination observed value without the influence of an ionized layer is obtained;lambda is the wavelength of the ionosphere-free combined observed value as the observation error;the geometric distance between the satellite position at the time of signal transmission to the receiver position at the time of signal reception.
3. The method for estimating the high frequency clock offset of a satellite navigation system according to claim 2, wherein the pseudo range of the clock offset and the phase error equation are expressed as:
<math> <mrow> <mi>V</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>V</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>V</mi> <mi>&Phi;</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>A</mi> <mi>p</mi> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>A</mi> <mi>&Phi;</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfenced open='(' close=')'> <mtable> <mtr> <mtd> <msub> <mi>L</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>L</mi> <mi>&Phi;</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>
in formula (6): <math> <mrow> <msub> <mi>A</mi> <mi>p</mi> </msub> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mover> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> <mi>A</mi> </mover> </mtd> <mtd> <mover> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> <mn>0</mn> </mover> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> <msub> <mi>A</mi> <mi>&Phi;</mi> </msub> <mo>=</mo> <mfenced open='[' close=']'> <mtable> <mtr> <mtd> <mover> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mn>2</mn> <mi>m</mi> <mo>+</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> <mi>A</mi> </mover> </mtd> <mtd> <mover> <mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&times;</mo> <mrow> <mo>(</mo> <mi>m</mi> <mo>&CenterDot;</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>&lambda;</mi> <mo>/</mo> <mi>C</mi> <mo>&CenterDot;</mo> <mi>I</mi> </mrow> </mover> </mtd> </mtr> </mtable> </mfenced> <mo>;</mo> </mrow> </math> let I be an n-order identity matrix, A be a design matrix for determining the clock offset of the receiver and the satellite, and be an (m.n) (m + n) matrix, V (I) be the observation error of the ith observation position, and L (I) be the ithThe free terms of observed position, phi and p respectively represent phase and pseudo range observation, X (i) is the ith parameter to be solved, V (i), Lp(i)、LΦ(i) Is a known amount; the number of m reference stations and n is the number of synchronously tracking satellites.
4. A high frequency clock error estimation method for a satellite navigation system according to claim 3, characterized in that a receiver clock in the tracking network is used as a reference clock, and the clock error accuracy of the reference clock is better than 10-6s; the relative clock error of all satellites in the same navigation system is zero mean.
5. The high frequency clock error estimation method for a satellite navigation system according to claim 4, wherein the reference clock is a hydrogen atomic clock.
CN201210541485.6A 2012-12-13 2012-12-13 A kind of high-frequency clock error estimation method of satellite navigation system Active CN102998681B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210541485.6A CN102998681B (en) 2012-12-13 2012-12-13 A kind of high-frequency clock error estimation method of satellite navigation system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210541485.6A CN102998681B (en) 2012-12-13 2012-12-13 A kind of high-frequency clock error estimation method of satellite navigation system

Publications (2)

Publication Number Publication Date
CN102998681A CN102998681A (en) 2013-03-27
CN102998681B true CN102998681B (en) 2015-09-09

Family

ID=47927489

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210541485.6A Active CN102998681B (en) 2012-12-13 2012-12-13 A kind of high-frequency clock error estimation method of satellite navigation system

Country Status (1)

Country Link
CN (1) CN102998681B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11841443B2 (en) * 2021-03-08 2023-12-12 Microchip Technology Incorporated Scalable common view time transfer and related apparatuses and methods

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104749587B (en) * 2013-12-31 2017-03-29 清华大学 Receiver pseudorange fault monitoring method and receiver
CN104503223B (en) * 2014-12-17 2017-01-25 同济大学 GNSS (Global Navigation Satellite System) three-frequency high-precision satellite clock correction estimating and service method
CN108107455A (en) * 2017-10-30 2018-06-01 千寻位置网络(浙江)有限公司 A kind of satellite clock correction Real-time Forecasting Method based on phase hit
CN108196279B (en) * 2017-12-23 2021-10-15 航天恒星科技有限公司 Satellite clock error calculating and forecasting method based on real-time data flow
CN107966722B (en) * 2018-01-20 2021-07-13 中国人民解放军61540部队 GNSS clock error resolving method
CN109001771B (en) * 2018-06-04 2020-10-23 北京未来导航科技有限公司 Navigation satellite and low-orbit satellite real-time clock error determining and forecasting method and system
CN109239750B (en) * 2018-09-10 2021-07-02 贵州省水利水电勘测设计研究院有限公司 Method and device for monitoring track parameters and electronic terminal
CN109061696B (en) * 2018-09-28 2022-12-09 中国人民解放军61540部队 Method for determining orbit and clock error of navigation satellite
CN110488328B (en) * 2019-07-18 2022-03-08 北京未来导航科技有限公司 Message receiving and sending method and system for low-orbit satellite navigation enhancement platform
CN110398758A (en) * 2019-07-24 2019-11-01 广州中海达卫星导航技术股份有限公司 Detection of Gross Errors method, apparatus, equipment and storage medium in real-time clock bias estimation
CN111045034B (en) * 2019-12-13 2020-09-29 北京航空航天大学 GNSS multi-system real-time precise time transfer method and system based on broadcast ephemeris
CN111045066A (en) * 2019-12-30 2020-04-21 威海欧瑞亚信息科技有限公司 Method for determining GNSS position change based on parameter equivalence reduction principle
CN111323796B (en) * 2020-03-18 2021-11-09 中国科学院国家空间科学中心 GNSS receiver high-sampling clock error resolving method
CN111413719B (en) * 2020-03-21 2022-07-15 哈尔滨工程大学 Beidou real-time precise clock prediction method based on neural network
CN111478725B (en) * 2020-05-08 2021-11-23 中国人民解放军63921部队 Satellite clock difference adjustment correction method based on inter-satellite link closed residual error detection
CN112540387A (en) * 2020-11-11 2021-03-23 湖南星南信息技术有限公司 Carrier phase cycle slip detection and restoration method based on clock hopping receiver
CN115598676B (en) * 2022-10-17 2023-05-05 北京航天飞行控制中心 Satellite-borne multimode GNSS fusion precise orbit determination method and device

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6778136B2 (en) * 2001-12-13 2004-08-17 Sirf Technology, Inc. Fast acquisition of GPS signal

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
GPS非差相位精密单点定位理论与实现;叶世榕;《中国学位论文全文数据库》;20021231;第I~III、14~22、50~57、71、78页 *
季善标,朱文耀,熊永清.星载GPS定轨中精密卫星钟差的改正和应用.《导航》.1999,(第3期),第100-105页. *
张益泽,陈俊平,王解先,.基于历元间差分的GNSS精密卫星钟差加密.《第三届中国卫星导航学术年会电子文集》.2012,第1-2页. *
李星星,徐运,王磊.非差导航卫星实时/事后精密钟差估计.《武汉大学学报》.2010,第35卷(第6期),第660-664页. *
汪志明,花向红,刘炎炎.精密单点定位钟卫星钟差影响分析.《武汉大学学报》.2010,第35卷(第11期),第1339-1341页. *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11841443B2 (en) * 2021-03-08 2023-12-12 Microchip Technology Incorporated Scalable common view time transfer and related apparatuses and methods

Also Published As

Publication number Publication date
CN102998681A (en) 2013-03-27

Similar Documents

Publication Publication Date Title
CN102998681B (en) A kind of high-frequency clock error estimation method of satellite navigation system
CN109581452B (en) GNSS reference station carrier phase integer ambiguity resolution method
CN111045034B (en) GNSS multi-system real-time precise time transfer method and system based on broadcast ephemeris
CN104570011A (en) Relative positioning device for satellite navigation and carrier phase cycle-slip repairing method of device
CN107966722B (en) GNSS clock error resolving method
CN105699999A (en) Method for fixing narrow lane ambiguity of Beidou ground based augmentation system base station
CN114355418A (en) Beidou foundation enhancement system-based posterior data quality assessment method and system
CN111044972B (en) GNSS precision time synchronization-based aircraft time difference positioning method and system
CN113325446B (en) Multimode common-frequency GNSS carrier phase time transfer method and system
CN111766616A (en) Beidou second-order time transfer satellite-side multipath error correction method
CN112235033A (en) GNSS common-view system based on RDSS communication link
Defraigne et al. GPS Time and Frequency Transfer: PPP and Phase-Only Analysis.
CN115993617B (en) GNSS system time deviation monitoring method
CN112146557A (en) GNSS-based real-time bridge deformation monitoring system and method
CN114879239B (en) Regional three-frequency integer clock error estimation method for enhancing instantaneous PPP fixed solution
CN115933356B (en) High-precision time synchronization system and method for virtual atomic clock
CN110068848B (en) High-performance RTK processing technical method
Yin et al. A novel cycle slips detection model for the high precision positioning
CN112731487A (en) GNSS (Global navigation satellite System) co-seismic displacement determination method based on high-stability atomic clock
CN116577815A (en) Multi-frequency multi-GNSS precise single-point positioning method, device and equipment
Li et al. Assessment and analysis of the four-satellite QZSS precise point positioning and the integrated data processing with GPS
CN115308781A (en) BDGIM assistance-based phase smoothing pseudorange high-precision time transfer method
Esteban et al. A GPS calibration trip experience between ROA and PTB
CN106707311A (en) GPS based enhanced GLONASS RTK (Real-time Kinematic) positioning method
CN118091718B (en) Method for improving UT1 calculation accuracy through low orbit satellite downlink navigation signal

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant