CN116184464A - GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint - Google Patents

GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint Download PDF

Info

Publication number
CN116184464A
CN116184464A CN202310304121.4A CN202310304121A CN116184464A CN 116184464 A CN116184464 A CN 116184464A CN 202310304121 A CN202310304121 A CN 202310304121A CN 116184464 A CN116184464 A CN 116184464A
Authority
CN
China
Prior art keywords
track
orbit
time
constraint
state
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.)
Pending
Application number
CN202310304121.4A
Other languages
Chinese (zh)
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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202310304121.4A priority Critical patent/CN116184464A/en
Publication of CN116184464A publication Critical patent/CN116184464A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/25—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS
    • G01S19/258—Acquisition or tracking or demodulation of signals transmitted by the system involving aiding data received from a cooperating element, e.g. assisted GPS relating to the satellite constellation, e.g. almanac, ephemeris data, lists of satellites in view
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/40—Correcting position, velocity or attitude

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention belongs to the technical field of real-time precise orbit determination of navigation satellites, and particularly discloses a GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint, which comprises the following steps: acquiring the track state at the current moment, constructing an observation equation and a state equation containing a state transition matrix at the current moment, and updating time by using the state equation; acquiring state parameters of an ultra-fast track at the current moment and standard deviation of a constraint equation, calculating to obtain the correction of the ultra-fast track relative to a filtering track, and constructing the constraint equation; measuring and updating according to the observation equation and the constraint equation to obtain an estimated value of the track state correction at the current moment, and calculating to obtain updated track state parameters; has the following advantages: the relation between the precision of the forecasting part of the ultra-fast track and the precision of the filtering track is defined, and mathematical modeling is carried out on the relation; and the constraint time is set, and a constraint equation is added, so that the track convergence time can be further shortened, and the speed is improved.

Description

GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint
Technical Field
The invention relates to the technical field of real-time precise orbit determination of navigation satellites, in particular to a GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint.
Background
The Global Navigation Satellite System (GNSS) is a large-scale infrastructure for performing navigation positioning on the ground or space users by using radio wave signals broadcast by artificial earth satellites, and has the advantages of globality, all weather, high precision and the like. The world is developing its own GNSS systems, such as the GPS system in the united states, the Galileo system in the european union, the GLONASS system in russia, the BDS system in china, etc. Satellite orbit is a spatial reference of a navigation system, and a real-time navigation positioning time service (PNT) service based on broadcast ephemeris is the most basic form of GNSS real-time PNT service. At present, most real-time positioning navigation users adopt broadcast ephemeris as real-time orbits, the orbit precision of GPS broadcast ephemeris is better than 1m, BDS broadcast ephemeris orbit precision GEO satellites are about 10m, and IGSO/MEO satellites are about 2m. It follows that the accuracy of broadcast ephemeris as a real-time orbit is low and cannot meet the high accuracy requirement of users. To improve the accuracy of real-time PNT services, precise Point Positioning (PPP) techniques may be used to correct the observations of a flow station in the state domain. The PPP technology can make the data processing of the user end get rid of the dependence on the reference station, but because the influence of the orbit and the satellite clock error can not be eliminated by the difference of the observables between stations, the real-time high-precision satellite orbit and clock error products are the key for realizing the high-precision real-time PPP, so that the research on the GNSS real-time precise orbit determination has important military and commercial significance.
At present, the real-time precise orbit determination mainly comprises two modes of batch processing prediction and expansion filtering, and the defect of the batch processing prediction mode is that if the state parameters and the dynamics parameters of a satellite at the initial moment are not accurate enough, the orbit precision can be drastically reduced after long-time integration. In addition, the batch processing forecasting mode needs to store a large amount of observation data for post-batch processing, so that the calculation efficiency is low, the time consumption is long, and the real-time requirement of a user cannot be met. The extended filtering mode does not estimate the state of the track at the reference moment, but updates the state and the dynamic parameters of the track epoch by epoch, so that the data storage capacity can be reduced, the calculation efficiency is improved, and meanwhile, errors generated by equation linearization due to inaccurate initial parameters of the track in the post mode can be reduced. And when the satellite encounters maneuver or fault, the user can be informed in time to take a proper coping mode. The common filtering models of the mode are extended Kalman filtering, square Root Information Filtering (SRIF), adaptive robust filtering and the like. The problem with the extended filtering mode is that the orbit converges to the centimeter level with a long convergence time during which satellite accuracy is poor.
At present, studies have been proposed to improve the convergence problem of the real-time filtering orbit by using ultra-fast orbit constraint, mainly using the ultra-fast orbit to constrain real-time orbit parameters, and adding constraint epoch by epoch in the real-time resolving process. However, the standard deviation of the constraint equation is determined by means of an empirical value, universality is lacking, and the problem of convergence of a real-time filtering orbit is solved by utilizing an ultra-fast orbit, wherein the prior constraint is mainly performed on real-time orbit parameters by utilizing the ultra-fast orbit, and the constraint function is not realized in the epoch-by-epoch solution of the orbit parameters.
Disclosure of Invention
The invention aims to provide a real-time precise orbit determination method for GNSS satellites by utilizing ultra-fast orbit constraint, which is used for determining the standard deviation of constraint equations of BDS-3MEO satellites, GPS satellites and Galileo satellites, and restraining orbit parameters by utilizing a forecasting part of the ultra-fast orbit epoch by epoch so as to enhance the intensity of an orbit result and greatly reduce the convergence time of a real-time filtering orbit in the initial stage.
In view of this, a first aspect of the present invention is to provide a real-time precise orbit determination method for GNSS satellites using ultra-fast orbit constraints.
The first aspect of the present invention provides a real-time precise orbit determination method for a GNSS satellite using ultra-fast orbit constraint, comprising the steps of: s1, performing dynamics fitting on a precise ephemeris to obtain a track state and a state transition matrix at the initial moment of filtering, and setting constraint time and total processing time; s2, acquiring a track state at the current moment, and constructing an observation equation and a state equation containing a state transition matrix at the current moment; s3, judging whether the difference between the current moment and the previous moment is smaller than the constraint time length, if so, executing S4, and if not, executing S5; s4, obtaining the standard deviation of the ultra-fast track state parameter and the constraint equation at the current moment, calculating to obtain the correction of the ultra-fast track relative to the reference track and the correction of the filter track relative to the reference track, constructing the constraint equation, and adding the constraint equation into the measurement update of S5 to correct the estimation value of the track state correction at the current moment; s5, acquiring the correction of the track state parameter at the current moment according to the state equation, measuring and updating according to the observation equation to obtain an estimated value of the track state correction at the current moment, calculating to obtain an updated track state parameter, and outputting the updated track state parameter; and S6, judging whether the difference between the current time and the previous time is smaller than the total processing time, if so, returning to S2 and changing the next time into the current time, and if not, ending.
Further, the selection time of the precise ephemeris is 48 hours before, and the total processing time is 72 hours; the constraint duration is 14 hours for the GPS satellite, 16 hours for the Galileo satellite, or 18 hours for the Beidou satellite.
Further, the constraint equation standard deviation is obtained by adopting the following formula:
Figure BDA0004146111760000031
Figure BDA0004146111760000032
wherein ,/>
Figure BDA0004146111760000033
Standard deviation representing the difference between the state correction of the ultra-fast track relative to the reference track and the state correction of the filter track relative to the reference track in the ECI coordinate system, G representing the rotation matrix,/o>
Figure BDA0004146111760000041
Representing the difference between the state correction of the ultra-fast track relative to the reference track and the state correction of the filter track relative to the reference track in the RTN coordinate systemStandard deviation of difference G T Representing the transpose of the rotation matrix.
Further, the rotation matrix G is obtained by the following formula: g= (G) 1 ,g 2 ,g 3 ) T The method comprises the steps of carrying out a first treatment on the surface of the Wherein T represents a matrix transpose operator, g 1 Representing the first column vector of the rotation matrix, and
Figure BDA0004146111760000042
g 2 a second column vector representing the rotation matrix, and +.>
Figure BDA0004146111760000043
g 3 A third column vector representing a rotation matrix, and g 2 ＝g 1 *g 3 R represents the position vector of the satellite at the current moment in the ECI coordinate system, < >>
Figure BDA0004146111760000044
Representing the velocity vector of the satellite at the current moment in the ECI coordinate system.
Further, the constraint equation standard deviation
Figure BDA0004146111760000045
For the difference between the state correction of the ultra-fast orbit relative to the reference orbit and the state correction of the filtered orbit relative to the reference orbit in the satellite orbit coordinate system +.>
Figure BDA0004146111760000046
Standard deviation of (2).
Further, the constraint equation is specifically the following formula:
Figure BDA0004146111760000047
wherein ,
Figure BDA0004146111760000048
representing state correction of ultrafast products relative to a reference track at a current timeNumber and
Figure BDA0004146111760000049
x(t i ) Representing the current track state correction,/-at the present moment>
Figure BDA00041461117600000410
Representing the ultrafast track state parameter, X, at the current instant * (t i ) Representing the reference track state at the current instant.
Further, the construction constraint equation is used to make the filtered orbit approach to the ultra-fast orbit, and the approach degree depends on
Figure BDA00041461117600000411
So that the convergence speed after the start of the filter track is increased.
Further, the step S2 specifically includes: s201, carrying out kinetic integration on the track state at the previous moment, obtaining the track state and a state transition matrix at the current moment, and constructing a state equation; s202, the pseudo range and the phase observation value of each observation station at the current moment are read, and an observation equation is constructed.
Further, the step of correcting in S4 specifically includes: in the filter, by setting
Figure BDA0004146111760000051
Standard deviation size of (2), restriction->
Figure BDA0004146111760000052
To reduce the track state correction estimate bias at the current time.
Compared with the prior art, the invention has the following beneficial effects:
firstly, a Square Root Information Filter (SRIF) is used as a filtering estimator, and compared with a Kalman filter, the method has higher numerical precision and more stable filtering solution; secondly, the relation between the precision of the forecasting part of the ultra-fast orbit and the precision of the converged filtering orbit is clarified, mathematical modeling is carried out on the relation, and a conclusion with universality is selected and determined for the standard deviation of constraint equations of different satellite systems; and finally, adding the standard deviation of the constraint equation under an RTN coordinate system, converting the standard deviation into an earth center inertial coordinate system (ECI), and adding the constraint epoch by epoch in the process of solving the filtering orbit.
Compared with the existing research, the orbit convergence time can be further shortened, wherein the real-time orbit convergence time of a GPS, galileo, BDS-3MEO satellite single system is reduced to be within one hour, the orbit precision is improved during the convergence period, scientific basis is provided for the selection of GPS, galileo, BDS-3MEO satellite constraint equation standard deviation, the sizes of different satellite system constraint equation standard deviations are precisely determined, and a conclusion of universality is formed.
Additional aspects and advantages of embodiments according to the invention will be apparent from the description which follows, or may be learned by practice of embodiments according to the invention.
Drawings
The drawings are only for purposes of illustrating particular embodiments and are not to be construed as limiting the invention.
FIG. 1 is a flow chart of a real-time precise orbit determination implementation of a GNSS satellite using ultra-fast orbit constraints according to the present invention;
FIG. 2 is a diagram of an observation station used in the real-time filtered orbit determination of the present invention;
FIG. 3 is a timing diagram of the comparison of the accuracy of a real-time filtered track using a conventional unconstrained real-time filtered track with the accuracy of a real-time filtered track using an ultrafast track as a constraint in accordance with the present invention;
FIG. 4 shows the one-dimensional and three-dimensional precision improvement of the real-time filtering track using the ultra-fast track as the constraint in comparison with the conventional unconstrained real-time filtering track during the constraint period;
FIG. 5 illustrates the convergence time of the present invention in tangential, normal, and radial directions using a conventional unconstrained real-time filtering trajectory;
fig. 6 shows the convergence time of the real-time filtered trajectory in tangential, normal and radial directions using the ultrafast trajectory as a constraint according to the present invention.
Detailed Description
In order that the above-recited objects, features and advantages of the present invention will be more clearly understood, a more particular description of the invention will be rendered by reference to the appended drawings and appended detailed description. It should be noted that, in the case of no conflict, the embodiments of the present application and the features in the embodiments may be combined with each other.
In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention, however, the present invention may be practiced in other ways than those described herein, and therefore the scope of the present invention is not limited to the specific embodiments disclosed below.
Referring to fig. 1-6, a first aspect of the present invention provides a real-time precise orbit determination method for GNSS satellites using ultra-fast orbit constraint, comprising the steps of:
step 1, dynamics fitting is carried out on the precise ephemeris (WUM) for the previous 48 hours, and extrapolation is carried out to obtain a filtering initial moment orbit state X 0 And state transition matrix
Figure BDA0004146111760000071
Setting constraint time length to be len and total processing time length to be T;
the constraint time length is determined according to the traditional filtering orbit convergence time, the constraint time length of the GPS satellite is 14 hours, the constraint time length of the Galileo satellite is 16 hours, the constraint time length of the Beidou satellite is 18 hours, and the total processing time length is T, specifically 72 hours.
Step 2, for t i-1 Track state at time
Figure BDA0004146111760000072
Kinetic integration is carried out to obtain t i Reference track state X of time * (t i ) And state transition matrix->
Figure BDA0004146111760000073
And constructing a state equation and carrying out time updating.
GNSS satellites are subject to different forces on their trajectories during motion, known as the filmed motion of the satellites. Simplifying a satellite into one particle, its equation of motion is as follows:
Figure BDA0004146111760000074
wherein ,
Figure BDA0004146111760000075
the acceleration vector of the satellite is represented by r, v, p and a, respectively, the three-dimensional position vector of the satellite, the three-dimensional velocity vector of the satellite, the dynamic model parameter of the satellite and the acceleration of the satellite under the influence of various acting forces. To solve the equation of motion, the above equation is converted into a first order differential equation, as shown in the following equation:
Figure BDA0004146111760000076
wherein ,
Figure BDA0004146111760000077
first derivative representing the three-dimensional position vector of the satellite,/->
Figure BDA0004146111760000078
Representing the first derivative of the satellite velocity vector, +.>
Figure BDA0004146111760000079
Representing the first derivative of the satellite dynamics model parameters;
further convert into:
Figure BDA0004146111760000081
wherein X represents a satellite state vector, and consists of a three-dimensional position vector r, a three-dimensional speed vector v and a dynamics parameter vector p, X 0 Representing the initial conditions at the reference instant,
Figure BDA0004146111760000082
One representing satellite state vectorsThe derivative of order, t, represents the current moment,/-, and>
Figure BDA0004146111760000083
indicated at t 0 Satellite state vector of moment.
Because the medium-high orbit satellite is more stable under the constraint of orbit mechanics and the orbit attenuation is smaller in a short time, the real-time orbit is obtained through satellite state parameter integration.
In a given initial state X 0 In the case of (2), the reference state X of the satellite can be obtained by numerical integration of the orbit * Linearizing the motion equation at the reference state to obtain a linear expression of the motion equation of the satellite:
Figure BDA0004146111760000084
Figure BDA0004146111760000085
/>
wherein x represents the correction of the satellite state relative to the reference state, A represents the partial derivative of the right function relative to the satellite reference state, I represents the identity matrix,
Figure BDA0004146111760000086
Partial derivatives of the acceleration vector representing the satellite with respect to the satellite reference state position vector, and>
Figure BDA0004146111760000087
representing the partial derivative of the acceleration vector representing the satellite with respect to the satellite reference state velocity vector,/for the satellite>
Figure BDA0004146111760000088
Representing the partial derivatives of the acceleration vectors of the satellites with respect to the satellite reference state dynamics model parameters.
The state equation is obtained according to the above formula, specifically the following formula:
Figure BDA0004146111760000089
wherein ,x(ti ) Representing the satellite state correction at the current time,
Figure BDA0004146111760000091
Representing the current time state transition matrix, x (t) i-1 ) Representing the satellite state correction at the previous moment, and +.>
Figure BDA0004146111760000092
The following formula is adopted for solving:
Figure BDA0004146111760000093
wherein ,
Figure BDA0004146111760000094
representing the partial derivative of the satellite current time position vector relative to the previous time position vector,/for the satellite>
Figure BDA0004146111760000095
Representing the partial derivative of the satellite current time position vector relative to the previous time velocity vector,/for the satellite>
Figure BDA0004146111760000096
Representing the partial derivative of the satellite current time position vector relative to the dynamics model parameters of the previous time,/and the like>
Figure BDA0004146111760000097
Representing the partial derivative of the satellite current time velocity vector relative to the previous time position vector,/for the satellite>
Figure BDA0004146111760000098
Representing the partial derivative of the satellite current time velocity vector relative to the previous time velocity vector,/for the satellite>
Figure BDA0004146111760000099
And the partial derivative of the speed vector at the current moment of the satellite relative to the dynamic model parameter at the previous moment is expressed, each piece of satellite position information is connected with the state parameter at the previous moment, and meanwhile, the individual gross error orbit position data is removed, and the state parameter is estimated.
And step 3, reading GNSS satellite observation data of all observation stations at the current moment, carrying out ionosphere-free combination to obtain an observation equation, and carrying out error correction.
After the state equation of the GNSS satellite is established, the orbit parameter needs to be estimated by using the observation value of the ground observation station, so that the observation equation of the GNSS satellite needs to be established. Establishing an observation equation by using a pseudo-range and a phase observation value received by a ground tracking station, wherein the pseudo-range and the phase observation value are two types of combined models and non-combined models;
the observed quantity at two frequencies is selected, and the observation equation is expressed as:
Figure BDA0004146111760000101
/>
Figure BDA0004146111760000102
wherein omc represents a linear pre-assay residual error after error model correction,
Figure BDA0004146111760000103
Pseudo-range observations, u, representing ionosphere-free combinations r Representing partial derivative, x of observed value relative to station coordinates under ground fixed coordinate system r Coordinate correction vector representing station under ground fixed coordinate system>
Figure BDA0004146111760000104
Representing the direction cosine, x of the station to the satellite under the inertial coordinate system j,s Representing inertial coordinate systemCorrection vector of lower track parameter relative to reference track, u erp Representing the partial derivative, x, of the observed value with respect to ERP parameters erp Correction vector representing ERP parameters relative to initial values, +.>
Figure BDA0004146111760000105
Deviation of receiver clock relative to GPST, indicative of the moment of signal reception>
Figure BDA0004146111760000106
Representing the intersystem deviation of the different systems with respect to GPST +.>
Figure BDA0004146111760000107
Representing the satellite clock difference at the time of signal transmission,
Figure BDA0004146111760000108
Representing tropospheric projection function, T w,r Representing zenithal tropospheric moisture content, < > and +.>
Figure BDA0004146111760000109
Pseudo-range residual error indicative of ionosphere-free combination, +.>
Figure BDA00041461117600001010
Phase observations representing ionosphere-free combinations, +.>
Figure BDA00041461117600001011
Representing carrier phase ambiguity,/v>
Figure BDA00041461117600001012
Indicating the phase residual error of the ionosphere-free combination.
And step 4, judging whether the difference between the current moment and the previous moment exceeds the constraint time, if not, executing the step 5, and if so, executing the step 6.
Step 5, if the constraint time length is not exceeded, obtaining t i Ultrafast track state parameters for time of day
Figure BDA00041461117600001013
Standard deviation from constraint equation>
Figure BDA00041461117600001014
Calculating the correction of the ultrafast track relative to the reference track: />
Figure BDA00041461117600001015
Constructing a constraint equation: />
Figure BDA00041461117600001016
Constraints are added to the track parameters epoch by epoch. And (5) carrying out measurement updating on the parallel contract beam equation and the observation equation, and calculating to obtain an estimated value of the orbit state correction.
After the filtering is started, the orbit parameters can be converged in a longer time due to the limitation of parameters such as priori precision, geometric structures and the like. The orbit parameters are constrained by the dynamic model, so that the orbit parameters have strong regularity and have high prediction accuracy when the dynamic model is accurate. The prediction part of the post-solution can be used as an ultrafast orbit to constrain the filtering orbit so as to enhance the strength of the equation solution. The constraint equation is expressed as:
Figure BDA0004146111760000111
wherein ,
Figure BDA0004146111760000112
representing the correction of the track state parameter of the ultrafast product at the current moment relative to the reference track state parameter, x (t) i ) Representing the correction of the filtered track state parameter at the current time relative to the reference track state parameter,
Figure BDA0004146111760000113
Is the difference between the state correction of the filtering orbit relative to the reference orbit and the state correction of the ultra-fast orbit relative to the reference orbit under the satellite orbit coordinate system (RTN), and is markedThe standard deviation is +.>
Figure BDA0004146111760000114
As a result of:
Figure BDA0004146111760000115
Figure BDA0004146111760000116
wherein ,
Figure BDA0004146111760000117
representing the ultrafast track state parameter, X, at the current instant * (t i ) Reference track status parameter indicative of the current moment,/->
Figure BDA0004146111760000118
Representing the track state parameter at the current time.
The constraint equation can therefore be rewritten as:
Figure BDA0004146111760000119
it can be seen that in this way,
Figure BDA00041461117600001110
the actual difference between the current time filtering track state parameter and the ultra-fast track state parameter.
By constructing constraint equation, the filtered orbit approaches to ultra-fast orbit, the approach degree is defined by
Figure BDA00041461117600001111
To ultimately improve convergence performance after the start of the filter trajectory. By changing->
Figure BDA00041461117600001112
I.e. to implement constraint equation standard deviations of different sizes.
The traditional real-time filtering orbit determination only uses a state equation for time updating, and uses an observation equation for measurement updating, so as to calculate the orbit state correction. In the invention, firstly, a state equation is used for time updating, and then, an observation equation and a constraint equation are combined for measurement updating, and the track state correction is calculated, so that the real-time filtering track is quickly converged.
The currently released ultra-fast product comprises a 24-hour measured part and a 24-hour forecast part of a BDS, GPS, galileo satellite, and the sampling rate is 300s, namely t i With the increment of i being 300s, the forecasting part in the product is used for restraining the real-time precise orbit determination of the GNSS satellite. To explore the relationship between the ultrafast orbit prediction part and the accuracy of the filtered orbit, a quadratic function was used to model the difference in accuracy of the predicted orbit and the filtered orbit within 24 h:
Figure BDA0004146111760000121
wherein ,
Figure BDA0004146111760000122
variance representing track correction, acr_dif pre The prediction part representing the ultrafast orbit in each epoch is tangential, normal and radial differences with respect to WUM precise ephemeris, acr_dif SRIF Representing the difference in tangential, normal and radial of the filtered orbit relative to WUM refined ephemeris in each epoch, f (x) represents the variance model function. Modeling is performed by adopting a quadratic function aiming at tangential, normal and radial results.
In adding constraints, the standard deviation of constraint equations is selected to be added in tangential, normal and radial directions of a satellite orbit coordinate system (RTN), and because a geocentric inertial coordinate system (ECI) is generally used for precise orbit determination of a satellite, constraint vectors under the RTN coordinate system are required to be converted into the ECI coordinate system through a rotation matrix, so that the satellite orbit is constrained. The conversion process is as follows:
Figure BDA0004146111760000131
wherein ,
Figure BDA0004146111760000132
standard deviation representing difference between state correction of filter track relative to reference track and state correction of ultra-fast track relative to reference track under ECI coordinate system, G represents rotation matrix, G T Representing the transpose of the rotation matrix +.>
Figure BDA0004146111760000133
And the standard deviation of the difference between the state correction of the filtering track relative to the reference track and the state correction of the ultra-fast track relative to the reference track under the RTN coordinate system is shown.
G＝(g 1 ,g 2 ,g 3 ) T
Wherein T represents a matrix transpose operator, g 1 Representing the first column vector, g, of the rotation matrix 2 Representing the second column vector, g, of the rotation matrix 3 A third column vector representing a rotation matrix, and g 1 、g 2 and g3 The following formula is adopted for calculation:
Figure BDA0004146111760000134
wherein r represents a position vector of the satellite at the current moment under an ECI coordinate system,
Figure BDA0004146111760000137
Representing the velocity vector of the satellite at the current moment in the ECI coordinate system.
Step 6, if the constraint time is exceeded, measuring and updating are carried out according to an observation equation, and the track state correction x (t) is calculated i ) Estimate of (2)
Figure BDA0004146111760000135
Updating track state parameter from epoch to epoch>
Figure BDA0004146111760000136
The method is the position information of the satellite under the ECI coordinate system, and the quick orbit determination of the satellite is completed in real time.
Example 1
As shown in fig. 3-6, 28 days of data from 1 in 10 months 2020 to 28 days in 10 months 2020 are processed, three days of filtering tracks are continuously calculated, and the ground observation station is shown in fig. 3:
firstly, resolving a real-time filtering orbit of a GNSS satellite single system under unconstrained conditions to obtain precision change time sequence diagrams of satellites of different systems in tangential, normal and radial directions;
then, using a forecast part of the ultra-fast orbit as external constraint, carrying out GNSS satellite single-system real-time filtering orbit determination to obtain a satellite precision change time sequence diagram under the constraint, and comparing the two, as shown in figure 3;
from fig. 3-6, it is derived that the convergence performance of the single-system filtering orbit of the GNSS satellite can be greatly shortened by adopting the ultra-fast orbit as an external constraint. Compared with a filtering orbit without constraint, the addition of external constraint improves the one-dimensional and three-dimensional precision of the GPS satellite in the constraint period by 86.7 percent and 92.3 percent respectively; the one-dimensional and three-dimensional precision of the Galileo satellite is improved by 84.7 percent and 93.6 percent respectively; the one-dimensional and three-dimensional precision of the BDS-3MEO satellite is respectively improved by 96.9 percent and 96.9 percent;
when no external constraint is added, the GPS, galileo, BDS-3MEO satellite convergence time is: 5.00/10.25/13.75 hours, 6.25/15.00/15.25 hours, 17.50/15.50/17.75 hours. After constraint is added, the tangential and radial directions of the GPS and Galileo satellites have no convergence process, and the radial convergence time is 0.75 hour; the BDS-3MEO satellite has tangential, normal and radial convergence time of 0.75/0.5/0.75 hours.
After the track is converged, the track accuracy added with the constraint is basically consistent with the track accuracy without the constraint, so that the track accuracy is improved only during the convergence period, the convergence time is shortened, and no influence is caused on the track after the convergence.
In the description of the present invention, it should be understood that the terms "longitudinal," "transverse," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," "outer," and the like indicate or are based on the orientation or positional relationship shown in the drawings, merely to facilitate description of the present invention, and do not indicate or imply that the devices or elements referred to must have a particular orientation, be constructed and operated in a particular orientation, and thus should not be construed as limiting the present invention.
The above embodiments are only illustrative of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention, and various modifications and improvements made by those skilled in the art to the technical solutions of the present invention should fall within the protection scope defined by the claims of the present invention without departing from the design spirit of the present invention.

Claims (9)

1. The GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint is characterized by comprising the following steps:
s1, performing dynamics fitting on a precise ephemeris to obtain a track state and a state transition matrix at the initial moment of filtering, and setting constraint time and total processing time;
s2, acquiring a track state at the current moment, and constructing an observation equation and a state equation containing a state transition matrix at the current moment;
s3, judging whether the difference between the current moment and the previous moment is smaller than the constraint time length, if so, executing S4, and if not, executing S5;
s4, obtaining the standard deviation of the ultra-fast track state parameter and the constraint equation at the current moment, calculating to obtain the correction of the ultra-fast track relative to the reference track and the correction of the filter track relative to the reference track, constructing the constraint equation, and adding the constraint equation into the measurement update of S5 to correct the estimation value of the track state correction at the current moment;
s5, according to a state equation, time updating is carried out to obtain the correction of the track state parameter at the current moment, and according to an observation equation, measurement updating is carried out to obtain an estimated value of the track state correction at the current moment, and the updated track state parameter is calculated and output;
and S6, judging whether the difference between the current time and the previous time is smaller than the total processing time, if so, returning to S2 and changing the next time into the current time, and if not, ending.
2. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 1, wherein the precise ephemeris is selected for the first 48 hours and the total processing time is 72 hours;
the constraint duration is 14 hours for the GPS satellite, 16 hours for the Galileo satellite, or 18 hours for the Beidou satellite.
3. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 1, wherein the constraint equation standard deviation is obtained using the following formula:
Figure FDA0004146111750000021
wherein ,
Figure FDA0004146111750000022
representing the difference between the state correction of the ultra fast track relative to the reference track and the state correction of the filtered track relative to the reference track in ECI coordinate system>
Figure FDA0004146111750000023
Is represented by the standard deviation, G, of the rotation matrix, +.>
Figure FDA0004146111750000024
Representing the difference between the state correction of the ultra fast track relative to the reference track and the state correction of the filtered track relative to the reference track in the RTN coordinate system +.>
Figure FDA0004146111750000025
Standard deviation of G T Representing the transpose of the rotation matrix.
4. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 3, wherein the rotation matrix G is obtained by using the following formula:
G＝(g 1 ,g 2 ,g 3 ) T ；
wherein T represents a matrix transpose operator, g 1 Representing the first column vector of the rotation matrix, and
Figure FDA0004146111750000026
g 2 a second column vector representing the rotation matrix, and +.>
Figure FDA0004146111750000027
g 3 A third column vector representing a rotation matrix, and g 2 ＝g 1 *g 3 R represents the position vector of the satellite at the current moment in the ECI coordinate system, < >>
Figure FDA0004146111750000028
Representing the velocity vector of the satellite at the current moment in the ECI coordinate system.
5. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 3, wherein the constraint equation standard deviation
Figure FDA0004146111750000029
For the difference between the state correction of the ultra-fast orbit relative to the reference orbit and the state correction of the filtered orbit relative to the reference orbit in the satellite orbit coordinate system>
Figure FDA00041461117500000210
Standard deviation of (2). />
6. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 5, wherein the constraint equation is specifically the following formula:
Figure FDA0004146111750000031
wherein ,
Figure FDA0004146111750000032
representing the state correction of the ultrafast product relative to the reference track at the current moment, x (t) i ) Indicating the current track state correction.
7. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 6, wherein the constraint equation is constructed to approximate the filtered orbit to the ultra-fast orbit, and the degree of the approximation depends on
Figure FDA0004146111750000033
So that the convergence speed after the start of the filter track is increased.
8. The method for real-time precise orbit determination of GNSS satellites by using ultra-fast orbit constraint according to claim 1, wherein the step S2 specifically comprises:
s201, carrying out kinetic integration on the track state at the previous moment, obtaining the track state and a state transition matrix at the current moment, and constructing a state equation;
s202, the pseudo range and the phase observation value of each observation station at the current moment are read, and an observation equation is constructed.
9. The method for real-time precise orbit determination of a GNSS satellite using ultra-fast orbit constraint according to claim 6, wherein S4, specifically: in the filter, by setting
Figure FDA0004146111750000034
Standard deviation of size of (2) limit
Figure FDA0004146111750000035
To reduce the track state correction estimate bias at the current time. />
CN202310304121.4A 2023-03-27 2023-03-27 GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint Pending CN116184464A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310304121.4A CN116184464A (en) 2023-03-27 2023-03-27 GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310304121.4A CN116184464A (en) 2023-03-27 2023-03-27 GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint

Publications (1)

Publication Number Publication Date
CN116184464A true CN116184464A (en) 2023-05-30

Family

ID=86444496

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310304121.4A Pending CN116184464A (en) 2023-03-27 2023-03-27 GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint

Country Status (1)

Country Link
CN (1) CN116184464A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116840879A (en) * 2023-09-04 2023-10-03 中国科学院国家授时中心 Method and system for determining clock error of low-orbit satellite considering orbit constraint

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116840879A (en) * 2023-09-04 2023-10-03 中国科学院国家授时中心 Method and system for determining clock error of low-orbit satellite considering orbit constraint
CN116840879B (en) * 2023-09-04 2023-12-08 中国科学院国家授时中心 Method and system for determining clock error of low-orbit satellite considering orbit constraint

Similar Documents

Publication Publication Date Title
CN110275192B (en) High-precision single-point positioning method and device based on smart phone
CN111596321B (en) Multi-GNSS multi-path error star day filtering method and system using non-difference correction
Li et al. Review of PPP–RTK: Achievements, challenges, and opportunities
CN111947667A (en) Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination
CN103176188A (en) Single-epoch fixing method for enhancing PPP-RTK ambiguity of regional foundation
CN113687402B (en) Low-orbit navigation enhancement real-time positioning method considering satellite orbit errors
CN113253314B (en) Time synchronization method and system between low-orbit satellites
CN116840879B (en) Method and system for determining clock error of low-orbit satellite considering orbit constraint
CN111856534A (en) Dual-mode GNSS carrier precise single-point positioning method and system of intelligent terminal
CN113325446B (en) Multimode common-frequency GNSS carrier phase time transfer method and system
CN116184464A (en) GNSS satellite real-time precise orbit determination method utilizing ultra-fast orbit constraint
CN115902968A (en) PPP terminal positioning method based on Beidou third GEO broadcast enhancement information
CN114545458A (en) High-precision ionosphere real-time modeling method
CN114779301B (en) Satellite navigation real-time precise single-point positioning method based on broadcast ephemeris
CN112526564A (en) Precise single-point positioning re-convergence method
CN115616643A (en) City area modeling auxiliary positioning method
CN116203598A (en) Ionosphere modeling method, device and medium based on foundation and star-based enhanced fusion
CN115421172A (en) Beidou deformation monitoring method based on real-time and quasi-real-time combination
Bahadur et al. Real-time single-frequency multi-GNSS positioning with ultra-rapid products
CN113433576A (en) GNSS and V-SLAM fusion positioning method and system
CN116299586B (en) Precise single-point positioning method, receiver, equipment and medium based on broadcast ephemeris
CN109738912B (en) Method for realizing fixed point time service based on GNSS satellite signals
CN116147622A (en) Combined navigation system fusion positioning method based on graph optimization
Wang et al. GPS un-differenced ambiguity resolution and validation
Elsheikh Integration of GNSS precise point positioning and inertial technologies for land vehicle navigation

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