CN111947667B - Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination - Google Patents
Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination Download PDFInfo
- Publication number
- CN111947667B CN111947667B CN202010591804.9A CN202010591804A CN111947667B CN 111947667 B CN111947667 B CN 111947667B CN 202010591804 A CN202010591804 A CN 202010591804A CN 111947667 B CN111947667 B CN 111947667B
- Authority
- CN
- China
- Prior art keywords
- orbit
- satellite
- precision
- gnss
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/02—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by astronomical means
-
- 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/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
- G01S19/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
-
- 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/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
- G01S19/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
- G01S19/072—Ionosphere corrections
-
- 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
-
- 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/393—Trajectory determination or predictive tracking, e.g. Kalman filtering
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Astronomy & Astrophysics (AREA)
- Automation & Control Theory (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The application relates to a low-orbit satellite real-time high-precision orbit determination method based on a combination of kinematics and dynamics, which comprises the steps of firstly, receiving and processing high-precision GNSS orbits, high-precision clock errors, earth rotation parameter forecasting data and solar radiation flux parameter data by a low-orbit satellite through a feed uplink link or an inter-satellite link, and then carrying out orbit recursion from the moment k by utilizing a high-precision satellite orbit dynamics model and combining a Longkutta method to obtain the position speed at the moment k + 1; performing cycle slip detection and repair on GNSS observation data at the k +1 moment; then correcting the orbit clock error at the time of k +1 according to the noted high-precision orbit clock error data product, and correcting errors such as GNSS antenna phase center deviation, low orbit satellite antenna phase center deviation, GNSS antenna phase center change correction, ionosphere and the like; and finally, after an observation matrix at the moment k +1 and a pseudo code carrier residual error matrix are obtained, Kalman filtering measurement updating is carried out, the position, the speed and the integer ambiguity after filtering are obtained, a state vector and a variance matrix are updated, and the calculation of the next epoch is carried out.
Description
Technical Field
The invention relates to the technical field of low-orbit satellite orbit determination, in particular to a low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination.
Background
At present, with the development of aerospace technology, the application field of low-orbit satellites is wider and wider, systems such as low-orbit surveying and mapping satellites, low-orbit navigation enhanced satellites, low-orbit gravity measurement satellites and low-orbit magnetic measurement satellites are gradually developed, and the premise key of the systems for realizing load tasks is to realize high-precision orbit determination.
Currently, for example, a resource three-dimensional surveying and mapping satellite carries out post-event high-precision orbit determination through a post-event high-precision orbit determination technology, then inversion of load data is carried out, load application and popularization are greatly affected by load processing duration, and a real-time high-precision orbit determination technology is urgently needed.
In addition, researchers also propose that the GNSS observation data is utilized, GNSS broadcast messages are used for calculating the positions and clock errors of the navigation satellites, and real-time orbit determination is carried out by combining a dynamic model.
The students also propose to download GNSS observation data on the low-orbit navigation enhancement satellite to the ground, upload the orbit determination result to the low-orbit satellite after the ground performs real-time orbit determination, and package and issue the navigation enhancement message after the low-orbit satellite receives the orbit determination result.
Disclosure of Invention
The invention aims to overcome one of the defects of the prior art, and provides a low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination, which can combine high-precision satellite orbit dynamics orbit determination with non-differential precision single-point positioning.
The invention discloses a low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination, which comprises the following steps:
the method comprises the following steps: the low earth orbit satellite receives high-precision GNSS orbit, high-precision clock error and earth rotation parameter forecasting data through a feed uplink or an inter-satellite link, and then lagrangian interpolation is carried out to obtain a corresponding value at the current moment;
step two: calculating the position speed of a recursion k +1 epoch satellite by combining the earth rotation parameter forecast data, the solar radiation flux forecast data and the like of the upper notes through a Longge Kutta method according to a high-precision satellite orbit dynamics model;
step three: detecting and repairing cycle slip gross error of the k +1 epoch GNSS original observation data;
step four: and (3) performing k +1 moment orbit and clock errors according to high-precision GNSS orbits and clock error products, and correcting various errors such as antenna phase center deviation, solid tide, ionized layer and the like according to other error models.
Step five: and updating Kalman filtering measurement to obtain parameters such as position, speed, integer ambiguity and the like after filtering, updating the state vector and the variance matrix, and calculating the next epoch.
Compared with the prior art, the orbit determination method has the following advantages:
the short-term recursion precision high advantage of a high-precision satellite orbit dynamics model and the error correction advantage of non-differential precision single-point positioning kinematics are combined, so that real-time high precision of orbit determination is realized. On the basis, the key technical support is realized for the low-orbit navigation enhancement satellite to develop the navigation enhancement service, and compared with the existing after-the-fact and ground orbit determination technology, the satellite-ground transmission link pressure can be reduced, and the low-orbit navigation enhancement service efficiency is improved.
Drawings
FIG. 1 is a flow chart of a low-orbit satellite real-time high-precision orbit determination method based on a combination of kinematics and dynamics;
FIG. 2 is a flowchart of a satellite orbit dynamics recursion calculation;
FIG. 3 is a schematic diagram of the transition between time systems;
FIG. 4 is a flow chart of data preprocessing;
FIG. 5 is a flow chart of GNSS error correction;
FIG. 6 is a flow chart of Kalman filtering parameter estimation;
Detailed Description
The invention is described in detail below with reference to the figures and specific embodiments. It is to be understood that the embodiments described are only a few embodiments of the present invention, and not all embodiments. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, shall fall within the scope of protection of the present invention.
A real-time high-precision orbit determination method for low-orbit satellites can combine high-precision dynamic orbit determination with non-differential precision single-point positioning, as shown in figure 1, and comprises the following steps:
the method comprises the following steps: and receiving and processing high-precision information of the ground upper notes.
The low-orbit satellite receives a high-precision GNSS orbit, a high-precision clock error, forecast earth rotation parameter data and solar radiation flux through a feed uplink or an inter-satellite link, and then lagrangian interpolation or extrapolation is respectively carried out to obtain a corresponding value at the current moment; the corresponding data provides input for satellite orbit dynamics recursion calculation and non-differential precise single point positioning calculation.
Step two: and (5) satellite orbit dynamics recursion calculation.
The method comprises the step of calculating the position speed of a recursion K +1 epoch satellite by a Longge Kutta method according to the position speed of the K epoch satellite and a high-precision satellite orbit dynamics model and by combining the earth rotation parameter forecast data, the solar radiation flux forecast data and the like. And the recursion calculation result is used as the input of the subsequent Kalman filtering measurement updating.
Step three: and preprocessing the k +1 epoch GNSS raw observation data.
The method mainly comprises cycle slip detection and restoration of the carrier phase observed value to obtain corrected original observation data, and provides guarantee for obtaining a high-precision result by subsequently performing error correction.
Step four: and correcting the GNSS error.
And eliminating the k +1 time orbit and the clock error according to the high-precision GNSS orbit and the clock error product. Furthermore, various errors such as antenna phase center deviation, solid tide, ionized layer and the like can be corrected according to other error models. And the data after error elimination or correction is used as the input of parameter estimation Kalman filtering measurement updating.
Step five: and performing Kalman filtering parameter estimation to obtain a filtered orbit determination result.
And updating Kalman filtering measurement to obtain the filtered position, speed and integer ambiguity, updating the state vector and the variance matrix, and calculating the next epoch.
The first step comprises the following steps:
step 101: receiving high-precision GNSS orbit, high-precision clock error, earth rotation parameter forecasting and solar radiation flux data, preferably receiving the data through a feed uplink or an inter-satellite link, and updating at regular time, wherein the high-precision GNSS orbit data product and the clock error data product are acquired from a ground data processing center in real time, the orbit data product updating time interval is less than 1min, the clock error data product updating time interval is less than 10s, the earth rotation parameter data is updated once a day, and the solar radiation flux data is updated once a day.
Step 102: adopting a Lagrange interpolation method or a linear extrapolation method to carry out interpolation to obtain GNSS orbit data and precision clock error data, earth rotation parameters and solar radiation flux at the current moment, wherein all four items of data can adopt Lagrange interpolation, and the Lagrange interpolation formula is as follows:
in the formula: f (x) k ) Is the function value at the interpolation node; l k (x) The basis function is interpolated n times.
The GNSS orbit data and the clock error data can also adopt a linear extrapolation method, and the formula for acquiring the precise GNSS orbit at the current moment by the GNSS orbit data linear extrapolation method is as follows:
δX＝[e radial e along e cross ]δO
e radial ＝e along ×e cross
X orbit ＝X broadcast -δX
in the formula, X orbit For corrected GNSS satellite geostationary position, X broadcast Calculating the position of the earth-fixed system for the GNSS satellite broadcast message, δ X being the correction number under the earth-fixed system, r being the broadcast messageThe position of the satellite being calculated is calculated,
the satellite velocity calculated for the broadcast message, δ O is the radial tangential normal deviation for the corresponding time instant, and t0 is the reference time instant in the correction information.
The calculation formula of the GNSS clock error linear extrapolation is as follows:
δC＝C 0 +C 1 (t-t 0 )+C 2 (t-t 0 ) 2
wherein t is satellite For corrected GNSS satellite time, t broadcast The on-satellite time calculated by the GNSS satellite broadcast message, t is the current time, t0 is the reference time in the correction information, and C0, C1 and C2 are the reference time clock error correction coefficients.
Through the first step, the GNSS precise orbit, the precise clock error, the earth rotation information and the solar radiation flux corresponding to the current moment are obtained, and input is provided for satellite orbit dynamics recursion calculation and non-error precise single-point positioning calculation.
As shown in fig. 2, the second step includes:
step 201: and initializing a state vector and a variance matrix by using the pseudo-range positioning result for the first epoch, wherein the non-first epoch takes the value after the last epoch is filtered.
Wherein the state vector comprises: the position, the speed, the clock error of a low-orbit receiver, an atmospheric resistance coefficient, a sunlight pressure coefficient and the whole-cycle ambiguity of each GNSS satellite;
the variance matrix of the state vector is a diagonal matrix, which is generally set to be larger, so that the filtering can be converged faster, in the example, the position corresponding variance is set to be 100, the speed corresponding variance is set to be 1000, the receiver clock error corresponding variance is set to be 100, the atmospheric resistance coefficient and the sunlight pressure coefficient corresponding variance are set to be 5, and the whole-cycle ambiguity corresponding variance is set to be 10.
Step 202: and (4) time system conversion, namely converting the current navigation system time into julian days, earth dynamics time and world time, and calculating a polar shift matrix, a time difference matrix, a nutation matrix and an earth rotation matrix to obtain a current time earth-fixed system and inertia system conversion matrix.
The time system conversion is shown in fig. 3, wherein leap seconds are periodically updated, and the current time UT1-UTC is obtained by interpolation of the daily terrestrial rotation parameter data.
The calculation formula of the conversion matrix from the earth fixed system to the inertia system is as follows:
E＝Prec*Nut*GAST*Pole
where Prec is a time matrix, Nut is a nutation matrix, GAST is an earth rotation matrix, and Pole is a polar shift matrix calculated from earth polar shift data annotated daily.
Step 203: the current satellite earth fixed system position is converted into an inertial system position, the aspheric perturbation acceleration under the earth fixed system is calculated according to the earth fixed system position and the earth gravity field model, the aspheric perturbation acceleration is converted into the aspheric perturbation acceleration of the inertial system, and the earth central gravity acceleration under the inertial system is calculated according to the inertial system position.
The earth gravity field model is an JGM3 model of 70 th order, and the gravity field potential function formula is as follows:
wherein Re and mu are respectively the reference radius of the earth and the gravitational constant, the value is 6378136.3 and 3.986004415 × 1014 according to JGM3,
λ is the satellite geocentric latitude and geocentric longitude.
Step 204: and calculating the position of the sun and moon at the current moment, calculating the gravitational acceleration of the sun and moon according to the position of the inertial system of the satellite, and calculating the perturbation acceleration of the earth solid tide according to the solid tide model.
Step 205: and calculating atmospheric resistance perturbation acceleration at the current moment according to the input solar radiation flux, the atmospheric density model and the current satellite position, and calculating sunlight pressure acceleration according to the sunlight pressure model.
Step 206: and adding all current perturbation accelerations to obtain the acceleration of an inertial system borne by the satellite at the current moment, and converting the acceleration into a ground fixed system.
The conversion formula of the acceleration of the inertial system to the acceleration of the ground fixed system is as follows:
wherein E is a transformation matrix from an inertial system to a ground-based system, and omega is the rotational angular velocity of the earth.
Step 207: and (3) repeating the step 202 to the step 206 according to a fourth-order five-level Runge Kutta method, calculating the satellite position velocity acceleration at the moment k +1, and obtaining a state transition matrix from the moment k to the moment k + 1.
The method comprises the step of calculating the position speed of the next epoch satellite by combining the earth rotation parameter forecast data, the solar radiation flux forecast data and the like according to the current satellite position speed and the high-precision satellite orbit dynamics model. And the recursion calculation result is used as the input of the subsequent Kalman filtering measurement updating.
Under a dynamic observation environment, due to the complex dynamic change of a carrier GNSS observation environment and multipath influence, cycle slip is easy to occur on an observation value of a signal by a GNSS receiver, wherein observation data mainly comprise observation data of carrier phase of the GNSS signal. Therefore, as shown in fig. 4, the following step three includes:
step 301: and carrying out data preprocessing on the original observation data to realize cycle slip detection on the carrier phase observation value.
The low-orbit GNSS receiver is a high-dynamic receiver, the cycle slip probability of observation data is high, and therefore, the MW combination and the ionosphere TEC change rate (TECR) are preferably adopted for double cycle slip detection.
Constructing an M-W combination using the pseudoranges and the phases such that the widelane ambiguity is calculated as:
in the formula phi 1 And phi 2 Respectively phase observations (unit: cycle), P at different frequencies 1 And P 2 Respectively pseudorange observations (in meters), L at different frequencies WL Is a wide lane combination (unit: meter), lambda WL Wide-lane wavelengths.
When N is present WL (k) If the following expression is satisfied, it is assumed that cycle slip has occurred in epoch k.
In the formula (I), the compound is shown in the specification,
representing the mean value of the widelane ambiguity, σ 2 Representing the variance of the widelane ambiguity.
The ionospheric change rate is used to detect the cycle slip of the observed value, and in a short time, the ionospheric (TECR) value can be regarded as a constant, and the ionospheric formula of epoch k-1 is as follows:
wherein γ ═ f 1 2 /f 2 2 ，f 1 And f 2 For the corresponding carrier frequency, b i And B p The signal inter-frequency offset, respectively at the receiver side and at the satellite side, can be considered constant over a period of time. Thus, the calculation of the ionospheric TEC rate of change TECR for epoch k is
In addition, because the change value of the ionosphere TECR is smooth in a short time, the prior epoch can be used for predicting the TECR of the k epoch, and the formula is as follows:
according to the characteristic that the change of the ionosphere TECR is gentle in a short time, the difference value between the TECR calculated value and the TECR predicted value of the current epoch is small under the condition that the cycle slip does not occur theoretically. Therefore, when the difference exceeds the threshold (0.15TECU/s), the epoch is considered to have a cycle slip.
Step 302: and restoring the original observation data by cycle slip calculation of the carrier phase observation value.
When a week jump in epoch k is detected by step 301, the following calculation equation may be listed according to the MW combination and TECR combination detection formula:
wherein a is an integer and b is a real number. The real value Δ N obtained by the above formula 1 、ΔN 2 Rounding to obtain L 1 Frequency sum L 2 And restoring the observed value by the whole cycle slip value on the frequency.
The method mainly comprises cycle slip detection and restoration of the carrier phase observed value to obtain corrected original observation data, and provides guarantee for obtaining a high-precision result by subsequently performing error correction.
As shown in fig. 5, the fourth step includes:
step 401: according to the position and clock error of the high-precision GNSS satellite at the moment of k +1, correcting the orbit error and clock error of the broadcast ephemeris;
step 402: correcting the phase deviation PCO of the GNSS antenna of the GNSS satellite observed at the k +1 moment;
step 403: correcting earth solid tide errors at the moment k +1, and correcting phase center deviation of a GNSS antenna of an on-satellite receiver according to the satellite position speed;
step 404: and calculating an observation elevation according to the low-orbit satellite position and the GNSS satellite position at the moment of k +1, and eliminating the low-elevation satellite to ensure the quality of observation data. Wherein the elevation threshold value preferably takes 15 °.
Step 405: and carrying out observation value double-frequency combination to eliminate ionosphere errors, and carrying out GNSS satellite antenna phase center change correction according to the current state.
According to the method, k +1 moment orbit and clock errors are eliminated according to high-precision GNSS orbits and clock error products, errors such as antenna phase center deviation, solid tide and ionosphere are corrected according to other error models, and data after errors are eliminated or corrected serve as input of parameter estimation Kalman filtering measurement updating.
As shown in fig. 6, the fifth step includes:
step 501: obtaining an observation matrix, a pseudo-range carrier residual error matrix and the like by eliminating satellite positions and pseudo-range carrier observation data of various errors;
the observation matrix is:
wherein nsat is the number of GNSS stars observed by the epoch;
the pseudo code carrier residual matrix is:
the means1 pseudo code is a pseudo code combined observed value of the low-orbit satellite receiver and the first GNSS satellite after all errors are eliminated, the means1 carrier is a carrier combined phase observed value of the low-orbit satellite receiver and the first GNSS satellite after all errors are eliminated, and N1 is the current whole-cycle ambiguity of the first satellite; r1 is the distance between the first GNSS satellite and the low-earth satellite after the errors are eliminated.
Step 502: bringing a Kalman filtering measurement updating equation into the Kalman filtering measurement updating equation to obtain the position, the speed, the integer ambiguity and the like after filtering;
the measurement update equation is:
wherein
Predicted value, Φ, from kinetic extrapolation in step 207 k/k-1 For the state transition matrix obtained in step 207, Rk is the observation noise matrix and Qk-1 is the state noise covariance matrix.
Step 503: and (4) repeating the steps 501 and 502 for 10 times, and after 10 iterations of Kalman filtering, performing equivalent position and velocity updating on the state vector and the variance matrix to calculate the next epoch.
In the step, the position, the speed and the integer ambiguity after filtering are obtained through Kalman filtering measurement updating, the state vector and the variance matrix are updated, and a real-time high-precision orbit determination result is obtained.
Therefore, according to the scheme, the advantages of high short-term recursion precision of the high-precision satellite orbit dynamics model are utilized, the recursion of the state variables among the epochs is carried out, then the high-precision information injected on the ground is combined to correct all errors of GNSS observation, and the parameter estimation is carried out through Kalman filtering on the basis of data preprocessing. The advantages of dynamic orbit determination and kinematic orbit determination are integrated, so that the low-orbit satellite can be subjected to real-time high-precision orbit determination.
While the invention has been described with reference to specific embodiments, the invention is not limited thereto, and those skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope of the invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
Claims (5)
1. A low-orbit satellite real-time high-precision orbit determination method based on a combination of kinematics and dynamics is characterized by comprising the following steps of:
the method comprises the following steps: the low-orbit satellite receives high-precision GNSS orbits, high-precision clock errors and earth rotation parameter forecasting data through a feed uplink or an inter-satellite link, and then lagrangian interpolation or extrapolation is carried out to obtain a corresponding value at the current moment;
step two: calculating the position and speed of a recursion k +1 epoch satellite by combining the earth rotation parameter forecast data and the solar radiation flux forecast data which are injected upwards through a four-order five-stage Runge Kutta method according to a high-precision satellite orbit dynamics model, wherein the high-precision satellite orbit dynamics model comprises an earth gravitational field model adopting 70 orders;
step three: detecting and repairing cycle slip of the k +1 epoch GNSS original observation data;
step four: correcting the orbit and clock error at the moment k +1 according to the GNSS orbit and the clock error product, and correcting the phase center deviation of the antenna, the solid tide and the error of an ionized layer according to a preset error model;
step five: performing Kalman filtering measurement updating to obtain filtered position, speed and integer ambiguity parameters, updating state vector and variance matrix, and calculating next epoch;
wherein, step four includes:
step 401: according to the position and clock error of the high-precision GNSS satellite at the moment of k +1, correcting the orbit error and clock error of the broadcast ephemeris;
step 402: correcting the phase deviation PCO of the GNSS antenna of the GNSS satellite observed at the k +1 moment;
step 403: correcting earth solid tide errors at the moment of k +1, and correcting phase center deviation of a GNSS antenna of an on-satellite receiver according to the satellite position speed;
step 404: calculating an observation elevation according to the low-orbit satellite position and the GNSS satellite position at the moment of k +1, and rejecting a low-elevation satellite;
step 405: and carrying out observation value double-frequency combination to eliminate ionosphere errors, and carrying out GNSS satellite antenna phase center change correction according to the current state.
2. The method for real-time high-precision orbit determination of low-orbit satellites based on the combination of kinematics and dynamics as claimed in claim 1, wherein the first step is specifically as follows:
step 101: acquiring a high-precision GNSS orbit data product and a clock error data product from a ground data processing center in real time;
step 102: interpolation or linear extrapolation is carried out by adopting a Lagrange interpolation method of more than nine orders to obtain GNSS orbit data and precision clock error data, earth rotation parameters and solar radiation flux at the current moment.
3. The method for real-time high-precision orbit determination of the low-orbit satellite based on the combination of kinematics and dynamics as claimed in claim 1, wherein the second step is specifically as follows:
step 201: initializing a state vector and a variance matrix by using a pseudo-range positioning result for a first epoch, and taking a value after filtering a previous epoch for a non-first epoch;
step 202: calculating a coordinate system conversion matrix, calculating current julian day time, earth dynamics time and world time by using rotation parameters, then calculating a rotation parameter matrix, a time difference matrix, a nutation matrix and an earth rotation matrix, and multiplying the four matrixes to obtain a current time earth-fixed system and inertia system conversion matrix;
step 203: converting the current satellite earth fixed system position into an inertial system position, and calculating the satellite earth center gravity inertial system acceleration and the non-spherical perturbation acceleration;
step 204: calculating the position of the sun and moon at the current moment, calculating the three-body gravitational acceleration of the sun and moon according to the position of the inertial system of the satellite, and calculating the perturbation acceleration of the earth solid tide according to the solid tide model;
step 205: calculating atmospheric resistance perturbation acceleration at the current moment according to the input solar radiation flux, the atmospheric density model and the current satellite position, and calculating sunlight pressure acceleration according to the sunlight pressure model;
step 206: adding all current perturbation accelerations to obtain the acceleration of an inertial system borne by the satellite at the current moment, and converting the acceleration into a ground fixed system;
step 207: and (4) repeating the step 202 to the step 206 according to a fourth-order five-level Runge Kutta method, and calculating the satellite position velocity acceleration at the moment of k + 1.
4. The method for real-time high-precision orbit determination of the low-orbit satellite based on the combination of kinematics and dynamics as claimed in claim 1, wherein the third step is specifically:
step 301: carrying out data preprocessing on the original observation data to realize cycle slip detection on the carrier phase observation value;
step 302: and the original observation data is repaired through cycle slip calculation of the carrier phase observation value.
5. The method for real-time high-precision orbit determination of low-orbit satellites based on the combination of kinematics and dynamics as claimed in claim 1, wherein the fifth step is specifically as follows:
step 501: obtaining an observation matrix and a pseudo-range carrier residual error matrix by eliminating satellite position and pseudo-range carrier observation data of various errors
Step 502: bringing a Kalman filtering measurement updating equation into the Kalman filtering measurement updating equation to obtain the position, the speed and the integer ambiguity after filtering;
step 503: and repeating the steps 501 and 502 for 10 times, updating the state vector and the variance matrix according to the position velocity value after 10 iterations of Kalman filtering, and calculating the next epoch.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010591804.9A CN111947667B (en) | 2020-06-24 | 2020-06-24 | Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010591804.9A CN111947667B (en) | 2020-06-24 | 2020-06-24 | Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111947667A CN111947667A (en) | 2020-11-17 |
| CN111947667B true CN111947667B (en) | 2022-08-12 |
Family
ID=73337290
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010591804.9A Active CN111947667B (en) | 2020-06-24 | 2020-06-24 | Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111947667B (en) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113009537B (en) * | 2021-02-18 | 2023-10-31 | 中国人民解放军国防科技大学 | Inertial navigation assisted defensive navigation relative positioning single epoch part ambiguity solving method |
| CN112558125B (en) * | 2021-02-22 | 2021-05-25 | 腾讯科技（深圳）有限公司 | Vehicle positioning method, related device, equipment and storage medium |
| CN113514856B (en) * | 2021-04-12 | 2024-01-09 | 国网上海市电力公司 | Cycle slip detection method, device, electronic equipment and readable storage medium |
| CN115308779B (en) * | 2021-05-07 | 2023-11-03 | 华为技术有限公司 | Ephemeris forecasting method and ephemeris forecasting device |
| CN113608247A (en) * | 2021-08-11 | 2021-11-05 | 中国科学院上海天文台 | Orbit determination method for satellite |
| CN114136325B (en) * | 2021-11-29 | 2024-02-13 | 中国科学院国家授时中心 | Beidou orbit maneuver satellite orbit determination method and system |
| CN114383619B (en) * | 2021-12-07 | 2023-09-05 | 上海航天控制技术研究所 | High-precision track calculation method |
| CN114488227B (en) * | 2022-01-26 | 2023-10-20 | 西南交通大学 | Multipath error correction method based on spatial correlation |
| CN115356749B (en) * | 2022-08-16 | 2023-09-08 | 广州爱浦路网络技术有限公司 | Satellite position locating method for low orbit satellite, computer device and storage medium |
| CN117194869B (en) * | 2023-11-07 | 2024-03-19 | 中国科学院国家授时中心 | Attitude-considered low-orbit satellite antenna phase center forecasting and fitting method |
| CN118129767B (en) * | 2024-05-07 | 2024-07-19 | 武汉大学 | Dynamic constraint low-orbit satellite GNSS real-time precise orbit determination method and system |
Family Cites Families (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7869948B2 (en) * | 2007-04-27 | 2011-01-11 | Sirf Technology, Inc. | Method and apparatus in positioning without broadcast ephemeris |
| US8260551B2 (en) * | 2008-01-10 | 2012-09-04 | Trimble Navigation Limited | System and method for refining a position estimate of a low earth orbiting satellite |
| CN101858981A (en) * | 2009-04-10 | 2010-10-13 | 马维尔国际贸易有限公司 | Method for realizing high sensitivity and quick first positioning of satellite navigation receiver |
| US20110238308A1 (en) * | 2010-03-26 | 2011-09-29 | Isaac Thomas Miller | Pedal navigation using leo signals and body-mounted sensors |
| CN103542854B (en) * | 2013-11-02 | 2016-03-23 | 中国人民解放军国防科学技术大学 | Based on the autonomous orbit determination method of satellite-borne processor |
| CN107356947B (en) * | 2017-05-31 | 2019-06-18 | 中国科学院测量与地球物理研究所 | The method for determining satellite difference pseudorange biases based on single-frequency navigation satellite data |
| CN107421550B (en) * | 2017-07-25 | 2020-08-28 | 北京航空航天大学 | Earth-Lagrange combined constellation autonomous orbit determination method based on inter-satellite ranging |
| CN108196272A (en) * | 2017-12-29 | 2018-06-22 | 中国电子科技集团公司第二十研究所 | A kind of satellite navigation positioning device and method based on real-time accurate One-Point Location |
| CN109001786B (en) * | 2018-06-04 | 2020-06-16 | 北京未来导航科技有限公司 | Positioning method and system based on navigation satellite and low-orbit augmentation satellite |
| CN110058282B (en) * | 2019-04-03 | 2023-06-09 | 南京航空航天大学 | PPP high-precision positioning method based on dual-frequency GNSS smart phone |
| CN110187364B (en) * | 2019-06-14 | 2023-06-09 | 火眼位置数智科技服务有限公司 | Low-rail navigation enhanced precision correction data generation and uploading system and method |
| CN111272336B (en) * | 2020-03-23 | 2021-02-19 | 中国科学院空间应用工程与技术中心 | Method for realizing mass center displacement estimation of large-scale low-orbit spacecraft based on GNSS observation |
-
2020
- 2020-06-24 CN CN202010591804.9A patent/CN111947667B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN111947667A (en) | 2020-11-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111947667B (en) | Low-orbit satellite real-time high-precision orbit determination method based on kinematics and dynamics combination | |
| CN109061696B (en) | Method for determining orbit and clock error of navigation satellite | |
| CN100501331C (en) | Navigation satellite autonomous navigation system and method based on X-ray pulsar | |
| Allahvirdi-Zadeh et al. | POD of small LEO satellites based on precise real-time MADOCA and SBAS-aided PPP corrections | |
| CN102844679B (en) | Gnss signal processing with regional augmentation positioning | |
| CN102498414B (en) | Gnss signal processing to estimate orbits | |
| CN111596321B (en) | Multi-GNSS multi-path error star day filtering method and system using non-difference correction | |
| CN101403790A (en) | Accurate one-point positioning method for single-frequency GPS receiver | |
| Bock et al. | GPS single-frequency orbit determination for low Earth orbiting satellites | |
| Hwang et al. | GPS‐Based Orbit Determination for KOMPSAT‐5 Satellite | |
| CN112817023A (en) | Satellite-based enhanced service-based unsupported low-orbit navigation enhancement system and method | |
| CN113624243B (en) | On-satellite real-time orbit forecasting method for near-earth orbit satellite | |
| Li et al. | Precise orbit determination for the FY-3C satellite using onboard BDS and GPS observations from 2013, 2015, and 2017 | |
| CN108196284A (en) | One kind is into the poor fixed GNSS network datas processing method of fuzziness single between planet | |
| Svehla | Geometrical theory of satellite orbits and gravity field | |
| CN107966722A (en) | A kind of GNSS satellite clock solutions method | |
| CN116338742A (en) | Real-time high-precision PNT service method based on combined space and earth observation resources | |
| CN110554443A (en) | Method for determining earth gravity field based on carrier phase observed value and point acceleration method | |
| Fernández et al. | Copernicus POD service operations—Orbital accuracy of Sentinel-1A and Sentinel-2A | |
| Defeng et al. | Reduced dynamic orbit determination using differenced phase in adjacent epochs for spaceborne dual-frequency GPS | |
| Li et al. | Precise orbit determination for LEO satellites: single-receiver ambiguity resolution using GREAT products | |
| CN112731504B (en) | Method and device for automatically determining orbit of lunar probe | |
| Ge et al. | Zero-reconvergence PPP for real-time low-earth satellite orbit determination in case of data interruption | |
| Gaur et al. | One-second GPS orbits: A comparison between numerical integration and interpolation | |
| Li et al. | Performance assessment of real-time orbit determination for the Haiyang-2D using onboard BDS-3/GPS observations |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |