Real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit constraints

Abstract

The application discloses to belong to the technical field of navigation satellite real-time precise orbit determination, and particularly discloses a GNSS satellite real-time precise orbit determination method using ultrafast orbit constraints, which comprises the following steps: acquiring an orbit state at a current time, constructing an observation equation and a state equation containing a state transition matrix at the current time, and performing time updating by using the state equation; acquiring state parameters of an ultrafast orbit at the current time and a standard deviation of a constraint equation, calculating correction numbers of the ultrafast orbit relative to a filtering orbit, and constructing the constraint equation; performing measurement updating according to the observation equation and the constraint equation to obtain an evaluation of the orbit state correction numbers at the current time, and calculating updated orbit state parameters; and the method has the following advantages: the relationship between the prediction part precision of the ultrafast orbit and the filtering orbit precision is determined, and a mathematical model of the relationship is established; the constraint time is set and the constraint equation is added, so that the orbit convergence time can be further shortened and the speed can be improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of real-time precise orbit determination of navigation satellites, in particular to a real-time precise orbit determination method of GNSS satellites using ultra-fast orbit constraints. BACKGROUND

[0002] Global Navigation Satellite System (GNSS) is a large infrastructure that uses artificial satellite broadcast radio wave signals to navigate and position ground or space users, with the advantages of globality, all-weather, high precision, etc. Countries around the world are developing their own GNSS systems, such as the GPS system of the United States, the Galileo system of the European Union, the GLONASS system of Russia, and the BDS system of China. Satellite orbit is the spatial reference of navigation system, and real-time navigation and positioning (PNT) service based on broadcast ephemeris is the most basic form of GNSS real-time PNT service. At present, most real-time positioning and navigation users use broadcast ephemeris as real-time orbit, and the orbit accuracy of GPS broadcast ephemeris is better than 1m, and the orbit accuracy of BDS broadcast ephemeris is about 10m for GEO satellites and about 2m for IGSO / MEO satellites. As can be seen, the accuracy of broadcast ephemeris as real-time orbit is low and cannot meet the high-precision requirements of users. In order to improve the accuracy of real-time PNT service, precise point positioning (PPP) technology can be used to correct the observations of the flow station in the state domain. The appearance of PPP technology can make the data processing of the user end free from the dependence on reference stations, but since the influence of orbit and satellite clock error cannot be eliminated through the difference between station observations, therefore, real-time high-precision satellite orbit and clock error products are the key to realizing high-precision real-time PPP, and therefore, the research on GNSS real-time precise orbit determination has important military and commercial significance.

[0003] At present, real-time precise orbit determination mainly has batch processing prediction and extended filtering two modes. The defect of batch processing prediction mode is that if the state parameters and dynamic parameters of the satellite at the initial time are not accurate enough, the orbit accuracy will decrease sharply after a long time of integration. And batch processing prediction mode needs to store a large amount of observation data for post-processing, which has low calculation efficiency and long time consumption, and cannot meet the real-time needs of users. The extended filtering mode does not estimate the state of the orbit at the reference time, but updates the state and dynamic parameters of the orbit every epoch, which can reduce the data storage amount, improve the calculation efficiency, and also can reduce the error caused by the linearization of the equation due to the inaccuracy of the initial parameters of the orbit in the post-processing mode. And when the satellite encounters a maneuver or a failure, the user can be informed in time to take appropriate measures. The commonly used filtering models of this mode include extended Kalman filter, square root information filter (SRIF), adaptive robust filter, etc. The problem of extended filtering mode is that it takes a long convergence time to converge to centimeter level, and the satellite accuracy is poor during the convergence period.

[0004] At present, there are studies that use ultra-quick orbit to improve the convergence of real-time filtering orbit, mainly using the ultra-quick orbit to constrain the real-time orbit parameters, and adding the constraint in the real-time solution process. However, the size of the standard deviation of the constraint equation is determined by experience, which lacks universality, and the use of ultra-quick orbit to improve the convergence of real-time filtering orbit has been proposed, which mainly uses the ultra-quick orbit to constrain the real-time orbit parameters, and does not play a constraining role in the orbit parameter solution process. SUMMARY

[0005] The application aims to provide a GNSS satellite real-time precise orbit determination method using ultra-quick orbit constraint, determines the size of the standard deviation of the constraint equation of BDS-3 MEO satellite, GPS satellite and Galileo satellite, and uses the prediction part of the ultra-quick orbit to constrain the orbit parameters every epoch to enhance the strength of the orbit results and greatly reduce the convergence time of the initial real-time filtering orbit.

[0006] Therefore, the first aspect of the application provides a GNSS satellite real-time precise orbit determination method using ultra-quick orbit constraint.

[0007] The first aspect of the application provides a GNSS satellite real-time precise orbit determination method using ultra-quick orbit constraint, comprising the following steps: S1, performing dynamic fitting on precise ephemeris to obtain the orbit state and state transition matrix at the start time of filtering, and setting the constraint time length and the total processing time length; S2, obtaining the orbit state at the current time, constructing the observation equation and the state equation containing the state transition matrix at the current time; S3, judging whether the difference between the current time and the last time is less than the constraint time length, if yes, executing S4, and if not, executing S5; S4, obtaining the ultra-quick orbit state parameter at the current time and the standard deviation of the constraint equation, calculating the correction number of the ultra-quick orbit relative to the reference orbit and the correction number of the filtering orbit relative to the reference orbit, to construct the constraint equation and modify the estimation of the orbit state correction number at the current time in the measurement update of S5; S5, obtaining the correction number of the orbit state parameter at the current time according to the state equation, and obtaining the estimation of the orbit state correction number at the current time according to the observation equation, and calculating the updated orbit state parameter and outputting; S6, judging whether the difference between the current time and the last time is less than the total processing time length, if yes, returning to S2 and setting the next time as the current time, and if not, ending.

[0008] Further, the selection time of the precise ephemeris is 48 hours, and the total processing time length is 72 hours; the constraint time length is 14 hours for GPS satellite, 16 hours for Galileo satellite, or 18 hours for Beidou satellite.

[0009] Further, the constraint equation standard deviation is obtained by using the following formula: wherein, represents the standard deviation of the difference between the state correction number of the ultra-fast orbit relative to the reference orbit and the state correction number of the filter orbit relative to the reference orbit in the ECI coordinate system, G represents the rotation matrix, represents the standard deviation of the difference between the state correction number of the ultra-fast orbit relative to the reference orbit and the state correction number of the filter orbit relative to the reference orbit in the RTN coordinate system, G T represents the transpose of the rotation matrix.

[0010] Further, the rotation matrix G is obtained by using the following formula: G=(g1, g2, g3) T ; wherein, T represents the matrix transpose operator, g1 represents the first column vector of the rotation matrix, and g2 represents the second column vector of the rotation matrix, and g3 represents the third column vector of the rotation matrix, and g2=g1*g3, r represents the position vector of the satellite at the current time in the ECI coordinate system, represents the velocity vector of the satellite at the current time in the ECI coordinate system.

[0011] Further, the constraint equation standard deviation is the standard deviation of the difference between the state correction number of the ultra-fast orbit relative to the reference orbit and the state correction number of the filter orbit relative to the reference orbit in the satellite orbit coordinate system.

[0012] Further, the constraint equation is specifically the following formula:

[0013]

[0014] wherein, represents the state correction number of the ultra-fast product relative to the reference orbit at the current time, and x(t i ) represents the orbit state correction number at the current time, represents the ultra-fast orbit state parameter at the current time, X * (t i ) represents the reference orbit state at the current time.

[0015] Further, the constraint equation is used to make the filter orbit approach the ultra-fast orbit, and the degree of approach depends on so as to accelerate the convergence speed of the filter orbit after starting.

[0016] ​Further, the step S2 specifically comprises: S201, performing dynamic integration on the orbit state of the previous moment to obtain the orbit state of the current moment and a state transition matrix, and constructing a state equation; S202, reading the pseudorange and phase observation value of each observation station at the current moment, and constructing an observation equation.

[0017] Further, the step of correcting in the S4 is specifically: in the filter, by setting the standard deviation size of , the change range of is limited to reduce the orbit state correction number estimation deviation of the current moment.

[0018] Compared with the prior art, the application has the beneficial effects of:

[0019] Firstly, the square root information filter (SRIF) is used as a filter estimator, which has higher numerical accuracy and more stable filter solution than the Kalman filter; secondly, the relationship between the prediction accuracy of the ultra-fast orbit and the accuracy of the converged filter orbit is determined, and the relationship is mathematically modeled, which determines the universality of the conclusion for selecting the standard deviation of the constraint equation of different satellite systems; finally, the standard deviation of the constraint equation is added in the RTN coordinate system, and then converted into the Earth-Centered Inertial (ECI) coordinate system, and the constraint is added every epoch in the filter orbit calculation process.

[0020] Compared with the prior art, the orbit convergence time can be further shortened, the real-time orbit convergence time of GPS, Galileo and BDS-3 MEO satellites is reduced to less than one hour, the orbit accuracy during the convergence period is improved, and a scientific basis is provided for the selection of the standard deviation of the constraint equation of GPS, Galileo and BDS-3 MEO satellites, the size of the standard deviation of the constraint equation of different satellite systems is accurately determined, and a universal conclusion is formed.

[0021] Additional aspects and advantages of embodiments of the application will become apparent from the following description, or will be learned by practice of embodiments of the application. BRIEF DESCRIPTION OF DRAWINGS

[0022] The accompanying drawings are included to provide a further understanding of the application and are incorporated in and constitute a part of this specification, illustrate embodiments of the application and together with the description serve to explain the principles of the application.

[0023] Figure 1 The flowchart of the real-time precise orbit determination of GNSS satellites using the ultra-fast orbit constraint of the application;

[0024] Figure 2 The observation station distribution map used in the real-time filter orbit determination of the application;

[0025] Figure 3The timing diagram for comparing the real-time filtering orbit accuracy of the application using the traditional unconstrained real-time filtering orbit and the real-time filtering orbit using the ultrafast orbit as a constraint;

[0026] Figure 4 The accuracy improvement of the real-time filtering orbit of the application using the ultrafast orbit as a constraint in one dimension and three dimensions during the constraint period relative to the traditional unconstrained real-time filtering orbit;

[0027] Figure 5 The convergence time of the traditional unconstrained real-time filtering orbit of the application in the tangential, normal, and radial directions;

[0028] Figure 6 The convergence time of the real-time filtering orbit of the application using the ultrafast orbit as a constraint in the tangential, normal, and radial directions. DETAILED DESCRIPTION

[0029] In order to more clearly understand the above-mentioned purposes, features and advantages of the application, the application will be further described in detail below in combination with the drawings and specific embodiments. It should be noted that the embodiments of the application and the features in the embodiments can be combined with each other without conflict.

[0030] In the following description, many specific details are set forth in order to provide a thorough understanding of the application, but the application can also be practiced without other different ways from those described herein, and therefore, the scope of protection of the application is not limited by the specific embodiments disclosed below.

[0031] Please refer to Figures 1-6 The first aspect of the application provides a GNSS satellite real-time precise orbit determination method using an ultrafast orbit constraint, comprising the following steps:

[0032] Step 1: Perform dynamic fitting on the precise ephemeris (WUM) of the previous 48 hours, and extrapolate to obtain the orbit state X0 at the start time of the filter and the state transition matrix Set the constraint duration as len and the total processing duration as T;

[0033] The constraint duration is determined according to the convergence time of the traditional filtering orbit, and the constraint duration of the GPS satellite is 14 hours, the constraint duration of the Galileo satellite is 16 hours, and the constraint duration of the Beidou satellite is 18 hours. The total processing duration T is 72 hours.

[0034] Step 2: Perform dynamic integration on the orbit state i-1 at time t to obtain the reference orbit state X i at time t * (t i ) and the state transition matrix The state equation is constructed and time is updated.

[0035] When GNSS satellites are in motion, they are subject to different forces on the trajectory of motion, which is called the perturbed motion of the satellite. If the satellite is simplified as a particle, its motion equation is as follows:

[0036]

[0037] wherein, represents the acceleration vector of the satellite, r represents the three-dimensional position vector of the satellite, v represents the three-dimensional velocity vector of the satellite, p represents the dynamic model parameters of the satellite, and a represents the acceleration of the satellite under the influence of various forces. In order to solve the motion equation, the above formula is converted into a first-order differential equation, as shown in the following formula:

[0038]

[0039] wherein, represents the first derivative of the three-dimensional position vector of the satellite, represents the first derivative of the velocity vector of the satellite, represents the first derivative of the dynamic model parameters of the satellite;

[0040] which is further converted into:

[0041]

[0042] wherein, X represents the state vector of the satellite, which is composed of the three-dimensional position vector r, the three-dimensional velocity vector v, and the dynamic parameter vector p, X0 represents the initial condition at the reference time, represents the first derivative of the state vector of the satellite, t represents the current time, represents the state vector of the satellite at t0.

[0043] Since the medium-high orbit satellite is relatively stable under the constraint of orbit mechanics, the orbit decay is small in a short time, so the real-time orbit is obtained by integrating the satellite state parameters.

[0044] Under the condition of given initial state X0, the reference state X of the satellite can be obtained by numerical integration of the orbit * The motion equation of the satellite is linearized at the reference state to obtain the linear expression of the motion equation:

[0045]

[0046]

[0047] wherein, x represents the correction amount of the state of the satellite relative to the reference state, A represents the partial derivative of the right function relative to the reference state of the satellite, and I represents the unit matrix, denotes the partial derivative of the acceleration vector of the satellite with respect to the state position vector of the satellite, denotes the partial derivative of the acceleration vector of the satellite with respect to the state velocity vector of the satellite, denotes the partial derivative of the acceleration vector of the satellite with respect to the state dynamics model parameter of the satellite.

[0048] The state equation is derived according to the above formula, and is specifically as follows:

[0049]

[0050] wherein x(t i ) denotes the satellite state correction number at the current time, denotes the state transition matrix at the current time, x(t i-1 ) denotes the satellite state correction number at the previous time, and The following formula is used for solving:

[0051]

[0052] wherein, denotes the partial derivative of the satellite position vector at the current time with respect to the position vector at the previous time, denotes the partial derivative of the satellite position vector at the current time with respect to the velocity vector at the previous time, denotes the partial derivative of the satellite position vector at the current time with respect to the dynamics model parameter at the previous time, denotes the partial derivative of the satellite velocity vector at the current time with respect to the position vector at the previous time, denotes the partial derivative of the satellite velocity vector at the current time with respect to the velocity vector at the previous time, denotes the partial derivative of the satellite velocity vector at the current time with respect to the dynamics model parameter at the previous time, the state parameter is estimated by associating each satellite position information with the state parameter at the previous time and eliminating individual coarse error orbit position data.

[0053] Step 3, reading the GNSS satellite observation data of all observation stations at the current time, obtaining the observation equation by performing ionosphere-free combination, and performing error correction.

[0054] After the state equation of the GNSS satellite is established, the observation value of the ground observation station needs to be used to estimate the orbit parameter, and therefore the observation equation of the GNSS satellite needs to be established. The pseudorange and phase observation values of each station are obtained by using the pseudorange and phase observation values received by the ground tracking station to establish the observation equation, and there are two kinds of combination models and non-combination models, and the ionosphere-free combination model is adopted in the present application.

[0055] The observation equation is represented as follows by selecting the observation values at two frequencies:

[0056]

[0057]

[0058] where omc represents the linearized pre-observation residual corrected by the error model, represents the ionosphere-free combined pseudorange observation, u r represents the partial derivative of the observation with respect to the station coordinates in the Earth-fixed coordinate system, x r represents the station coordinate correction vector in the Earth-fixed coordinate system, represents the station-to-satellite direction cosine in the inertial coordinate system, x j,s represents the orbit parameter correction vector in the inertial coordinate system with respect to the reference orbit, u erp represents the partial derivative of the observation with respect to the ERP parameters, x erp represents the ERP parameter correction vector with respect to the initial value, represents the receiver clock bias with respect to GPST at the signal reception time, represents the inter-system bias with respect to GPST, represents the satellite clock bias at the signal transmission time, represents the tropospheric projection function, T w,r represents the zenith tropospheric wet component, represents the ionosphere-free combined pseudorange residual error, represents the ionosphere-free combined phase observation, represents the carrier phase ambiguity, represents the ionosphere-free combined phase residual error.

[0059] Step 4, determine whether the difference between the current time and the previous time exceeds the constraint time length, if not, execute step 5, if yes, execute step 6.

[0060] Step 5, if the constraint time length is not exceeded, obtain the ultrafast orbit state parameter i at time t and the standard deviation of the constraint equation Calculate the correction number of the ultrafast orbit with respect to the reference orbit: Construct the constraint equation: Add constraints to the orbit parameters epoch by epoch. Jointly update the measurement with the constraint equation and the observation equation to obtain the estimated value of the orbit state correction.

[0061] After the filter is started, the orbit parameters need a long time to converge due to the limitations of the priori accuracy of parameters, geometric structure and other factors. The orbit parameters are constrained by the dynamic model and have strong regularity, and have high prediction accuracy when the dynamic model is relatively accurate. Therefore, the prediction part of the post-solution can be used as an ultra-fast orbit to constrain the filter orbit to enhance the strength of the equation solution. The constraint equation is expressed as:

[0062]

[0063] wherein, represents the correction number of the orbit state parameter of the ultra-fast product at the current time relative to the reference orbit state parameter, x(t i represents the correction number of the filter orbit state parameter at the current time relative to the reference orbit state parameter, is the difference between the state correction number of the filter orbit relative to the reference orbit in the satellite orbit coordinate system (RTN) and the state correction number of the ultra-fast orbit relative to the reference orbit, and the standard deviation is

[0064] Since:

[0065]

[0066]

[0067] wherein, represents the ultra-fast orbit state parameter at the current time, X * (t i represents the reference orbit state parameter at the current time, represents the orbit state parameter at the current time.

[0068] Therefore, the constraint equation can be rewritten as:

[0069]

[0070] As can be seen, is actually the difference between the filter orbit state parameter at the current time and the ultra-fast orbit state parameter.

[0071] By constructing the constraint equation, the filter orbit approaches the ultra-fast orbit, and the degree of approach is determined by , which ultimately improves the convergence performance of the filter orbit after starting. By changing the size of , the standard deviation of the constraint equation of different sizes is achieved.

[0072] The traditional real-time filtering orbit determination only uses the state equation to perform time updating and uses the observation equation to perform measurement updating to solve the orbit state correction number. In the application, the state equation is used to perform time updating, and then the observation equation and the constraint equation are combined to perform measurement updating to calculate the orbit state correction number, so that the real-time filtering orbit converges quickly.

[0073] The currently published ultra-fast product contains 24-hour measured part and 24-hour predicted part of BDS, GPS and Galileo satellites, and the sampling rate is 300s, i.e. i With the increase of i by 300s, the predicted part in the product is used to constrain the real-time precise orbit determination of GNSS satellites. In order to explore the relationship between the ultra-fast orbit prediction part and the filtering orbit accuracy, a quadratic function is used to model the accuracy difference between the predicted orbit and the filtering orbit within 24h:

[0074]

[0075] Among them, ACR_Dif represents the variance of the orbit correction number, ACR_Dif pre ACR_Dif represents the difference between the predicted part of the ultra-fast orbit and the WUM precise ephemeris in the tangential, normal and radial directions in each epoch, ACR_Dif SRIF ACR_Dif represents the difference between the filtering orbit and the WUM precise ephemeris in the tangential, normal and radial directions in each epoch, and f(*) represents the variance model function. For the results in the tangential, normal and radial directions, a quadratic function is used for modeling.

[0076] When adding the constraint, the standard deviation of the constraint equation is added in the tangential, normal and radial directions of the satellite orbit coordinate system (RTN). Since the satellite precise orbit determination usually uses the Earth-centered inertial coordinate system (ECI), it is necessary to convert the constraint vector in the RTN coordinate system to the ECI coordinate system through a rotation matrix, and then constrain the satellite orbit. The conversion process is as follows:

[0077]

[0078] Among them, ACR_Dif represents the standard deviation of the difference between the state correction number of the filtering orbit relative to the reference orbit and the state correction number of the ultra-fast orbit relative to the reference orbit in the ECI coordinate system, G represents the rotation matrix, and G T G represents the transpose of the rotation matrix, ACR_Dif represents the standard deviation of the difference between the state correction number of the filtering orbit relative to the reference orbit and the state correction number of the ultra-fast orbit relative to the reference orbit in the RTN coordinate system.

[0079] G=(g1,g2,g3) T

[0080] wherein T denotes a matrix transposition operator, g1 denotes a first column vector of a rotation matrix, g2 denotes a second column vector of a rotation matrix, g3 denotes a third column vector of a rotation matrix, and g1, g2 and g3 are calculated using the following formula:

[0081]

[0082] wherein r denotes a position vector of the satellite at a current time instant in an ECI coordinate system, denotes a velocity vector of the satellite at a current time instant in an ECI coordinate system.

[0083] Step 6, if the constraint time length is exceeded, measurement update is performed according to the observation equation, and an estimated value of the orbit state correction number x(t i ) is calculated The orbit state parameter is updated epoch by epoch i.e. the position information of the satellite in the ECI coordinate system, and real-time fast orbit determination of the satellite is completed.

[0084] Example 1

[0085] As Figures 3-6 , the data from October 1, 2020 to October 28, 2020 for 28 days is processed, and the filtered orbit is calculated continuously for three days. The ground observation stations are shown in Figure 3 .

[0086] First, the GNSS satellite single system real-time filtering orbit under no constraint is solved, and the accuracy variation time series of different system satellites in the tangential, normal and radial directions are obtained;

[0087] Then, the prediction part of the ultra-fast orbit is used as external constraint to perform GNSS satellite single system real-time filtering orbit determination, and the satellite accuracy variation time series under constraint is obtained, and the two are compared, as shown in Figure 3 .

[0088] From Figures 3-6 , it can be seen that using the ultra-fast orbit as external constraint can greatly shorten the convergence performance of the GNSS satellite single system filtering orbit. Compared with the filtering orbit under no constraint, the addition of external constraint makes the one-dimensional and three-dimensional accuracy of GPS satellites during the constraint period increase by 86.7%, 92.3% respectively; the one-dimensional and three-dimensional accuracy of Galileo satellites increase by 84.7%, 93.6% respectively; the one-dimensional and three-dimensional accuracy of BDS-3 MEO satellites increase by 96.9%, 96.9% respectively;

[0089] Without adding external constraints, the convergence time of GPS, Galileo and BDS-3 MEO satellites is: 5.00 / 10.25 / 13.75 hours, 6.25 / 15.00 / 15.25 hours, 17.50 / 15.50 / 17.75 hours. After adding constraints, there is no convergence process in the tangential and radial directions of GPS and Galileo satellites, and the radial convergence time is 0.75 hours; the convergence time of BDS-3 MEO satellites in the tangential, normal and radial directions is 0.75 / 0.5 / 0.75 hours.

[0090] After the orbit converges, the orbit accuracy of the orbit with constraints is basically the same as the orbit accuracy without constraints, which can show that the present application only improves the orbit accuracy during the convergence period and shortens the convergence time, and has no effect on the orbit after the convergence.

[0091] In the description of the present application, 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 the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present application, and do not indicate or imply that the device or element referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as a limitation on the present application.

[0092] The above-described embodiments are only preferred modes of the present application, and do not limit the scope of the present application, and various modifications and improvements to the technical solutions of the present application made by those skilled in the art without departing from the design spirit of the present application shall fall within the protection scope of the present application as defined by the claims.

Claims

1. A real-time precise orbit determination method using ultra-rapid orbit constraints of GNSS satellites, characterized in that, Comprise the following steps: S1, the precise ephemeris is carried out dynamic fitting, obtains the orbit state and state transition matrix of filter starting time, sets the constraint time length and total processing time length; S2, the orbit state of current time is obtained, the observation equation of current time is constructed and the state equation containing the state transition matrix of current time is constructed; S3, judge whether the difference between current time and last time is less than the constraint time length, if yes, execute S4, if not, execute S5; S4, the ultra-quick orbit state parameter of current time and the constraint equation standard deviation are obtained, the correction number of ultra-quick orbit relative to reference orbit and the correction number of filter orbit relative to reference orbit are calculated, the constraint equation is constructed and the estimation of orbit state correction number of current time is modified in the measurement update of S5; S5, according to the state equation, time update is carried out to obtain the correction number of orbit state parameter of current time, and according to the observation equation, measurement update is carried out to obtain the estimation of orbit state correction number of current time, the updated orbit state parameter is calculated and output; S6, judge whether the difference between current time and last time is less than the total processing time length, if yes, return to S2 and make the next time as current time, if not, end; The constraint equation standard deviation is obtained by the following formula: ； wherein denotes the difference between the state correction of the super-fast orbit with respect to the reference orbit and the state correction of the filter orbit with respect to the reference orbit in the ECI coordinate system the standard deviation of denotes the rotation matrix denotes the difference between the state correction of the super-fast orbit with respect to the reference orbit and the state correction of the filter orbit with respect to the reference orbit in the RTN coordinate system the standard deviation of denotes the transpose of the rotation matrix.

2. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 1, characterized in that, The selection time of the precise ephemeris is 48 hours, and the total processing time length is 72 hours; The constraint time length is 14 hours for GPS satellite, or 16 hours for Galileo satellite, or 18 hours for Beidou satellite.

3. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 1, characterized in that, the rotation matrix is obtained using the following equation: ； wherein denotes the matrix transpose operator, denotes a first column vector of the rotation matrix, and denotes a second column vector of the rotation matrix, and denotes a third column vector of the rotation matrix, and denotes a position vector of the satellite at the current time instant in the ECI coordinate system, denotes a velocity vector of the satellite at the current time instant in the ECI coordinate system.​​​ 4. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 1, characterized in that, the constraint equation standard deviation the standard deviation of the difference between the state correction of the superfast orbit relative to the reference orbit in the satellite orbital coordinate system and the state correction of the filtered orbit relative to the reference orbit the standard deviation of the difference between the state correction of the superfast orbit relative to the reference orbit in the satellite orbital coordinate system and the state correction of the filtered orbit relative to the reference orbit 5. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 4, characterized in that, The constraint equation is specifically the following formula: ； wherein, denotes the state correction of the ultra-rapid product at the current time instant with respect to the reference orbit, denotes the orbit state correction at the current time instant.

6. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 5, characterized in that, The construction constraint equation is used to approach the filter orbit to the super fast orbit, and the degree of approach depends on to make the convergence speed after the filter orbit starts faster.

7. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 1, characterized in that, The step of S2 specifically comprises: S201, the orbit state of previous time is carried out dynamic integral, the orbit state and state transition matrix of current time are obtained, and the state equation is constructed; S202, the pseudorange and phase observation value of each observation station of current time are read, and the observation equation is constructed.

8. The real-time precise orbit determination method of GNSS satellites using ultra-rapid orbit perturbations according to claim 5, characterized in that, The step of modifying the estimation of orbit state correction number of current time in S4 is specifically: In the filter, by setting the standard deviation size , the variation range of is limited to reduce the deviation of the orbit state correction estimate at the current time.

Images