A GNSS autonomous orbit determination method for a spacecraft

Abstract

The application provides a GNSS autonomous orbit determination method for a spacecraft, which obtains observation data in real time after GNSS signal acquisition and tracking, combines an orbit dynamics model and an observation model, realizes real-time autonomous orbit determination of the spacecraft in orbit, and improves the positioning and orbit determination precision of the spacecraft. Meanwhile, the application proposes a GNSS observation data preprocessing method based on polynomial fitting and a convergence judgment condition of the autonomous orbit determination filter, so as to ensure the stability of the autonomous orbit determination filter of the spacecraft. The GNSS autonomous orbit determination method for the spacecraft realized by the application can be widely applied to high-precision navigation, positioning and autonomous orbit determination service scenes of high, medium and low orbit spacecrafts, and has a wide application prospect.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of satellite application, and relates to a GNSS autonomous orbit determination method for a spacecraft, which can be applied to the field of autonomous navigation and high-precision orbit determination of high, medium and low spacecrafts in different orbits. BACKGROUND

[0002] Spacecraft autonomous orbit determination refers to real-time determination of the position of the spacecraft by using the equipment carried by the spacecraft without relying on ground support. The spacecraft autonomous orbit determination technology can realize autonomous position keeping and orbit control of the spacecraft, can reduce the dependence on manpower, material resources and equipment of the ground, and more importantly, can guarantee strong survivability of the spacecraft when information transmission between the spacecraft and the ground measurement and control system is interrupted, and has very important significance for improving the autonomous management capability of the task.

[0003] Traditional spacecraft generally adopts a geometric positioning method to determine the satellite orbit in real time, and the solving precision is not high. Especially in the scene of low measurement precision under poor geometric distribution of the GNSS constellation and low signal-to-noise ratio, only relying on the geometric positioning method cannot meet the task requirements of the spacecraft. In recent years, in order to improve the autonomous orbit determination precision of the spacecraft, some scholars have researched a method of using GNSS observations and an orbit dynamics model to perform real-time filtering orbit determination. In addition to being capable of improving the position and velocity measurement precision of the spacecraft, the method can also realize orbit recursion through the dynamics model when the number of GNSS satellites is insufficient, so as to ensure uninterrupted navigation and positioning data.

[0004] During the flight of the spacecraft, in addition to the effects of the earth gravity, the earth non-spherical gravity, the sun and moon two-body gravity, the solar pressure and the solid tide, the spacecraft will also be affected by non-conservative forces such as the solar pressure and the atmospheric resistance. Therefore, in view of the high-precision autonomous orbit determination requirement, how to establish an accurate and reasonable system model and an observation usage strategy, and how to design a spaceborne GNSS autonomous orbit determination filtering estimation algorithm suitable for the flight environment of the spacecraft, have important research significance. SUMMARY

[0005] The GNSS autonomous orbit determination method for a spacecraft provided by the application gives a GNSS real-time observation + dynamics model autonomous orbit determination filtering method in view of the high-precision autonomous orbit determination requirement of the spacecraft, proposes a GNSS observation preprocessing based on polynomial fitting and a convergence judgment condition of the autonomous orbit determination filter, and is helpful to improve the filtering precision and stability of the autonomous orbit determination of the spacecraft.

[0006] A GNSS autonomous orbit determination method for a spacecraft, comprising the following steps:

[0007] S1, a satellite receiver completes signal acquisition, tracking and demodulation, and real-time calculation and update of GNSS observations and ephemeris data;

[0008] S2, according to the GNSS observation and ephemeris data, the geometric positioning solution based on the least square method is completed, and the position, velocity and clock information of the receiver is obtained. If the autonomous orbit determination filter is not initialized, the position, velocity and clock information obtained by the geometric positioning solution is used as the filter initial value, and the autonomous orbit determination filter initialization is completed;

[0009] S3, the orbit dynamics model and the observation model are established, the state transition matrix calculation and time update of the autonomous orbit determination filter are completed, and the orbit prediction result at the new time is obtained;

[0010] S4, the current GNSS observation is preprocessed, and the measurement update of the autonomous orbit determination filter is completed;

[0011] S5, whether the updated autonomous orbit determination filter converges or not is judged;

[0012] S6, if the filter has converged, the autonomous orbit determination filter result is output; otherwise, the orbit prediction result after the time update of the filter in step S3 is output;

[0013] S7, steps S1-S7 are repeated.

[0014] The preprocessing process of the GNSS observation mentioned in step S1 is to judge the quality of the observation data. By polynomial fitting on the continuous observation and calculating the square error of each observation, when the error value is greater than the decision threshold, it is determined that the observation contains gross error.

[0015] The observation data quality judgment has the technical features that it includes the following specific steps:

[0016] S4.1, selecting the current observation time t of a certain navigation star k And the corresponding M observation values (y k y k-1 …y k-M-1 ) of the continuous M-1 observation time (t k t k-1 …t k-M-1 ) before it;

[0017] S4.2, according to the M observation values (y k y k-1 …y k-M-1 ) in step S4.1 and the corresponding observation time (t k t k-1 …t k-M-1 ), polynomial fitting is carried out, the difference between each observation value and the calculated value based on the fitted polynomial is calculated, and then the square of the difference value is calculated, that is, the square error;

[0018] S4.3, judging whether the square error value of the current observation time t k is greater than a decision threshold. If greater than the threshold, it is determined that the observation contains a gross error, and the observation is excluded. Otherwise, the observation is introduced into the autonomous orbit determination filter.

[0019] S4.4, after the observation update in the receiver, repeating steps S4.1-S4.4.

[0020] The polynomial fitting in step S4.2 has the technical features that the polynomial fitting is performed by the least square method, and the polynomial order is preferably 2.

[0021] In the execution process of the autonomous orbit determination filter, it is necessary to judge whether the filter converges. When the autonomous orbit determination filter convergence judgment condition satisfies one or more of the following, it is considered that the autonomous orbit determination filter has diverged, and needs to be reinitialized. The specific judgment conditions are as follows:

[0022] Judgment condition 1: the position information in the filter results after time update and measurement update is substituted into the GNSS observation equation to calculate the observation residual values of all navigation stars participating in the orbit determination filter, and the sum of the residual values is taken as the judgment quantity. When the judgment quantity continuously exceeds the set decision threshold, it is determined that the orbit determination filter diverges.

[0023] Judgment condition 2: the covariance values about the position vector in the state covariance matrices before and after the autonomous filter update are summed, and are denoted as the pre-test variance sum σ apr and the post-test variance sum σ pos . First, it is judged whether the post-test variance sum σ pos continuously exceeds the set threshold for multiple times, and then the difference between σ apr and σ pos is taken as the judgment quantity. When the difference continuously exceeds the set threshold for multiple times, it is also determined that the orbit determination filter diverges.

[0024] Judgment condition 3: the difference between the position obtained by geometric positioning at the current time in step S2 and the position obtained by autonomous orbit determination filter calculation. When the difference between the two continuously exceeds the set threshold for multiple times, it is determined that the orbit determination filter diverges.

[0025] The present application has the following beneficial effects:

[0026] (1) The present application calculates the square error of the observation by the polynomial fitting method, realizes the GNSS observation preprocessing method, can effectively detect and identify the gross error in the GNSS observation, and thus improves the autonomous orbit determination accuracy and reliability.

[0027] (2) The application monitors convergence state of the autonomous orbit determination filter in execution process in real time, and proposes three filter orbit determination filter convergence judgment conditions, which is helpful to improve reliability of the autonomous orbit determination filter. BRIEF DESCRIPTION OF DRAWINGS

[0028] Figure 1 A flow chart of the method of the application is shown in the figure.

[0029] Figure 2 A flow chart of the work of the observation gross error judgment of the application is shown in the figure. DETAILED DESCRIPTION

[0030] The application will be further described in detail below in combination with the annexes.

[0031] The spacecraft GNSS autonomous positioning and orbit determination method mainly includes the following two kinds: geometric single point positioning method, and autonomous orbit determination filter orbit determination method based on GNSS observation and dynamic model. The geometric single point positioning method is used in the spacecraft, which is seriously affected by GNSS observation accuracy and constellation geometric distribution. Therefore, the orbit determination method using GNSS observation and fusing orbit dynamic model information can further improve the navigation solution accuracy of high-orbit spacecraft.

[0032] The single point positioning solution based on the geometric method belongs to single-epoch solution of single observation moment, and the solution results of the front and rear epochs are independent and are not affected by spacecraft orbit maneuver. The geometric single point positioning method has low accuracy, but the geometric positioning result can be used to set the initial value of the filter when initializing the orbit determination filter or when reinitializing after detecting filter divergence. In the autonomous orbit determination filter process, the orbit determination filter result and the geometric positioning result can be compared for consistency, which is one of the bases for judging the effectiveness of the orbit determination filter result.

[0033] The application proposes a GNSS autonomous orbit determination method for spacecraft, and the processing flow is as shown in Figure 1 The specific steps include the following:

[0034] S1, the satellite receiver completes signal acquisition, tracking and demodulation, and real-time calculation and update of GNSS observation and ephemeris data;

[0035] S2, geometric positioning solution based on least square method is completed according to GNSS observation and ephemeris data, and position, velocity and clock information of the receiver is obtained. If the autonomous orbit determination filter is not initialized, the position, velocity and clock information obtained by the geometric positioning solution is used as the initial value of the filter, and the autonomous orbit determination filter initialization is completed;

[0036] S3, the orbit dynamic model and the observation model are established, the state transition matrix calculation and time update of the autonomous orbit determination filter are completed, and the orbit prediction result at the new time is obtained;

[0037] S4, preprocessing the current GNSS observation, completing the measurement update of the autonomous orbit determination filter;

[0038] S5, judging whether the updated autonomous orbit determination filter converges;

[0039] S6, if the filter has converged, outputting the autonomous orbit determination filter result; otherwise, outputting the orbit prediction result after the filter time update in step S3;

[0040] S7, repeating steps S1-S7.

[0041] The autonomous orbit determination algorithm for a spacecraft is closely related to the characteristics of GNSS observations and orbit motion. The orbit motion characteristics of a spacecraft are related to the dynamic model, and the motion equation is established for the spacecraft as follows:

[0042]

[0043] In the formula, (r, v, t) represent the position, velocity and time of the spacecraft in the inertial system respectively. m For the perturbation force model, modeling can be performed, including the earth's gravity, the earth's non-spherical gravity, the sun and moon two-body gravity, atmospheric resistance, solar pressure and solid tide, etc. For the influence of residual perturbation force acceleration model which is difficult to model, a compensation acceleration w(t) is introduced, i.e. dynamic model compensation (DMC), and its model is assumed to be a first-order Gauss-Markov process.

[0044] The spacecraft GNSS autonomous orbit determination filter is mainly based on the extended Kalman filter (EKF) algorithm. In the specific implementation process of the filter, taking GPS, GLONASS and BDS three systems as an example, first, the state vector is established, which is as follows:

[0045]

[0046] In the formula, (x, y, z) are three components of the receiver position in the inertial system, are three components of the receiver velocity in the inertial system, G , b R , b B represent the clock errors of GPS, GLONASS and BDS respectively, is the receiver clock drift, R is the light pressure parameter, R , w T , w N is the compensation acceleration of the empirical force model in the radial (Radial), normal (Normal) and tangential (Transverse), i.e. RTN.

[0047] The state vector X can be considered in two parts, one part is the state quantity related to the orbit dynamics model denoted as X1; the other part is the state quantity related to the receiver clock characteristics denoted as X2.

[0048] For the state vector related to the dynamics model, the satellite motion equation given by equation (2) is a first-order differential equation. If the initial state vector X1(t0) at time t0 is known, the state vector X1(t) at next time t can be calculated as follows:

[0049]

[0050] In the formula, T 3×3 is the conversion matrix of the RTN coordinate system to the space rectangular coordinate system, and Φ1(t, t0) is the state transition matrix, which can be expressed as:

[0051]

[0052] In the formula, Φ rr , Φ rv , Φ vr , Φ vv are 3x3 order spacecraft position and velocity vector transition matrices about themselves, respectively, are 3x1 order spacecraft position and velocity partial derivatives about the solar pressure coefficient, which can be obtained by the dynamics orbit integral method, and Φ rw , Φ vw , Φ ww are 3x3 order spacecraft position, velocity, and compensation acceleration transition matrices about the compensation acceleration.

[0053] The first term on the right side of equation (3) is the deterministic component of the spacecraft state, and the second term on the right side is the random term, which is the integral result of the compensation acceleration random term in the filter update time interval. Because u(t) is a zero-mean Gaussian white noise, although its result does not affect the estimated value of the state quantity, it affects the system noise covariance matrix.

[0054] Considering all the estimated parameters in the state vectors X1 and X2, the state equation Φ1(t, t0) is extended to obtain the complete state transition matrix of the orbit determination filter as follows:

[0055]

[0056] In the formula, the matrix B is the transition matrix of the GNSS navigation receiver clock error and clock drift itself.

[0057] The observation equation of the GNSS autonomous orbit determination filter based on EKF in the application is represented as:

[0058]

[0059] In the formula, Y k is an observation vector, V k is a measurement noise vector. At time t k , the element k in the vector Y is calculated from each navigation star pseudorange observation, and according to a pseudorange observation equation, the pseudorange of the i th navigation star at the observation epoch k is The calculation expression is:

[0060]

[0061] In the formula, δt s and δt u are navigation star clock errors and receiver clock errors respectively, the navigation star clock error can be corrected by a parameter in navigation data, c is the speed of light, δρ trop and δρ iono are troposphere and ionosphere errors respectively, ε ρ is pseudorange and carrier observation noise, r is the geometric distance between the navigation star and the receiver, and ρ is the pseudorange observation.

[0062] The observation matrix for the i th navigation star is obtained after linearization of the pseudorange observation equation, and the linearization process is:

[0063]

[0064] In the formula, is the position of the i th navigation star, r i is the geometric distance between the i th navigation star and the receiver, is the state vector predicted by the EKF , and T 1×3 is the result in the state vector, for GPS, [1, 0, 0], for GLONASS, [0, 1, 0], and for BDS, [0, 0, 1].

[0065] The GNSS autonomous orbit determination filter based on the EKF in the application needs to perform two processes of time update and measurement update in the execution process, and the specific processes are as follows:

[0066] (1) Time update

[0067] The time update process of the autonomous orbit determination filter includes state update and state error covariance matrix update, that is, according to the state quantity X k-1 and the covariance matrix P k-1 at the last time t k-1 , the predicted value of the state quantity at the current time t k is calculated Covariance Matrix The calculation process includes:

[0068] (2) Status update

[0069] Based on the spacecraft's equations of motion and dynamic model, with t k-1 The state variable X at time t is k-1 As initial values, a fourth-order Runge-Kutta integral (RK4) is used for numerical integration to obtain the dynamic integral prediction value. Then, the prediction value is corrected based on the influence of the deterministic component in the compensated acceleration, and t is calculated. k State quantity prediction value and t k-1 to t k The state transition matrix Φ at any given time k,k-1 .

[0070] (3) Covariance matrix update

[0071] Based on the calculation methods of the dynamic compensation model and the receiver clock bias model, the dynamic noise covariance matrix Q is calculated. k-1 Then, the covariance matrix is ​​updated over time as follows:

[0072]

[0073] In the formula, For t k The state error covariance matrix updated after time step 1. For t k-1 The state error covariance matrix after measurement and update at each time step.

[0074] (4) Measurement update

[0075] To reduce errors caused by linearization and word length limitations in the spacecraft CPU, and to avoid singular matrix problems that may be encountered in the inversion of the state error covariance matrix, the UD decomposition method is used to update the observations for EKF orbit determination filtering. During the measurement update process, measurement update calculations are performed sequentially for each GNSS pseudorange observation data point. Coefficient weighting is applied for different navigation systems. The calculation process is as follows, where the superscript i represents the i-th navigation satellite.

[0076] Step 1: Based on t k Time-based measurement of updated state quantities The observation matrix is ​​calculated by linearizing the observations using formula (8).

[0077] Step 2: Based on t k Covariance matrix updated at time measurement And formula (10), the filter gain matrix is calculated as follows. Wherein, R is the measurement error, the error weight corresponding to different navigation system is R g , r , b .

[0078]

[0079] Step 3: state vector update is carried out, and the state error covariance matrix update is calculated as follows:

[0080]

[0081]

[0082] In the navigation star observation update process, steps 1-3 are repeatedly executed until all GNSS observations are updated, the autonomous orbit determination filtering process is completed, and the state vector k at time t and the state error covariance matrix are output.

[0083] In the autonomous orbit determination filtering measurement update process, the state estimation value k at time t (after i-1 GNSS star measurement update state) can be calculated according to formula (7), and the value of , i.e., the pseudo-range residual value, is calculated, and the filter error covariance is also calculated.

[0084] In the autonomous orbit determination filtering execution time and measurement update process, it is also necessary to perform real-time inspection on whether the filtering converges, and the filtering convergence judgment and inspection method includes:

[0085] 1) According to the filtering results after time update and measurement update, the residual values of all participating orbit determination observations are calculated by substituting the GNSS observation equation, the residual values are summed and the variance of the residual values is calculated, and the sum is taken as the judgment quantity, when the judgment quantity continuously exceeds the set decision threshold, it is determined that the orbit determination filter diverges.

[0086] 2) The covariance values of the position vector in the state covariance matrix before and after the autonomous filter update are summed, denoted as the pre-test variance sum σ apr and the post-test variance sum σ pos . First, it is judged that when the post-test variance sum σ pos continuously exceeds the set threshold for multiple times, it is determined that the orbit determination filter diverges. Then, the difference between σ apr and σ pos is taken as the judgment quantity, when the difference continuously exceeds the set threshold for multiple times, it is also determined that the orbit determination filter diverges.

[0087] 3) The difference between the position result obtained by the current time geometric positioning and the position result obtained by the autonomous orbit determination filter solution, when the difference between the two is continuously more than the set threshold, the orbit determination filter is determined to diverge.

[0088] In the orbit determination filter working process, if the GNSS observation is abnormal, it will directly affect the orbit determination accuracy, and even cause the orbit determination filter to diverge. Therefore, it is necessary to develop the abnormal observation data judgment method, and the observation rough error judgment method in the application is shown as follows: Figure 2 As shown in the figure, by polynomial fitting on the continuous observation, and calculating the square error of each observation, when the error value is greater than the judgment threshold, it is determined that the observation contains rough error. The specific execution steps are as follows:

[0089] S4.1, selecting the current observation time t k of a certain navigation star and the continuous M-1 observation times before it, denoted as (t k t k-1 …t k-M-1 ), the M observation values corresponding to the M observation times are (y k y k-1 …y k-M-1 );

[0090] S4.2, according to the M observation values (y k y k-1 …y k-M-1 ) in step S4.1 and the observation times (t k t k-1 …t k-M-1 ), polynomial fitting is carried out, and the square error of the difference between each observation value and the calculated value based on the fitted polynomial is calculated;

[0091] S4.3, judging whether the square error value of the current observation time t k is greater than the judgment threshold, if yes, it is determined that the observation contains rough error, and it is removed;

[0092] S4.4, after the observation in the receiver is updated, steps S4.1-S4 are repeated.

[0093] In the application, polynomial fitting is carried out by least square method, and the square error of each observation is calculated. Taking a second-order polynomial as an example, the specific processing method is as follows:

[0094] The specific expression of the second-order polynomial is:

[0095] f(t)=b0+b1t+b2t 2

[0096] The observation values of M points continuously stored by the receiver form the following equation group:

[0097]

[0098] The above formula is simplified as follows by using matrix expression:

[0099]

[0100] Let

[0101]

[0102] The above formula is simplified as follows:

[0103] Y = HB

[0104] The solution is obtained as follows:

[0105] B = (H T H) -1 H T Y

[0106] After the coefficient matrix B is solved according to the observation values, the expression of the polynomial is determined. The square of the difference between each observation value and the calculated value based on the fitted polynomial, i.e., the square error, is obtained:

[0107]

[0108] If the square error value of the current observation time t k is greater than the decision threshold, it is determined that the observation contains a pseudorange gross error, and it is excluded. The value of the decision threshold is related to the square errors of the other M-1 times except the current time t k , and the mean value of the other M-1 square errors is preferably N times, and N is preferably 3.

[0109] Although the present application has been disclosed with the above preferred embodiments, it is not intended to limit the present application, and any person skilled in the art can make possible changes and modifications to the technical solutions of the present application by using the disclosed methods and technical contents without departing from the spirit and scope of the present application. Therefore, any simple modification, equivalent change and modification made according to the technical essence of the present application to the above embodiments, which does not depart from the technical solutions of the present application, belongs to the protection scope of the technical solutions of the present application.

Claims

1. A GNSS autonomous orbit determination method for spacecraft, characterized in that, Includes the following steps: S1. The satellite receiver completes signal acquisition, tracking, and demodulation, and calculates and updates GNSS observation and ephemeris data in real time; S2. Based on GNSS observations and ephemeris data, perform geometric positioning calculations based on the least squares method to obtain the receiver's position, velocity, and clock information. If the autonomous orbit determination filter is not initialized, use the position, velocity, and clock information obtained from the geometric positioning calculations as the initial values ​​of the filter to complete the initialization of the autonomous orbit determination filter. S3. Establish the orbit dynamics model and the observation model, complete the calculation of the state transition matrix and time update of the autonomous orbit determination filter, and obtain the orbit prediction results at the new time. S4. Preprocess the current GNSS observations and complete the measurement update of the autonomous orbit determination filter; The aforementioned preprocessing of the current GNSS observations involves judging the quality of the observation data by performing polynomial fitting on continuous observations and calculating the squared error of each observation. If the error value is greater than the decision threshold, the observation is considered to contain gross errors. The specific method is as follows: S4.1 Select the current observation time t of a certain navigation satellite k and the preceding M-1 consecutive observation times (t) k t k-1 …t k-M-1 The M observations (y) corresponding to ) k y k-1 … y k-M-1 ); S4.2, Based on the M observations (y) in step S4.1 k y k-1 … y k-M-1 ) and their corresponding observation time (t) k t k-1 …t k-M-1 ), perform polynomial fitting, calculate the difference between each observed value and the calculated value based on the fitted polynomial, and then square the difference, i.e., the squared error; S4.3 Determine the current observation time t k The squared error value is checked against the decision threshold. If it is greater than the threshold, the observation is considered to contain gross errors and is discarded. Otherwise, the observation is introduced into the autonomous orbit determination filter; S4.4 After the receiver's observations are updated, repeat steps S4.1 to S4.

4. S5. Determine whether the updated autonomous orbit determination filter has converged. S6. If the filter has converged, output the autonomous orbit determination filtering result; otherwise, output the orbit prediction result after the filter time update in step S3. S7. Repeat steps S1 to S7.

2. The GNSS autonomous orbit determination method for spacecraft according to claim 1, characterized in that, The polynomial fitting described in step S4.2 is performed using the least squares method, with the polynomial order being 2.

3. The GNSS autonomous orbit determination method for spacecraft according to claim 1, characterized in that, The decision threshold mentioned in step S4.2 is determined by removing the current time t. k It is related to the other M-1 squared errors, which are N times the mean of the other M-1 squared errors, where N is 3.

4. The GNSS autonomous orbit determination method for spacecraft according to claim 1, characterized in that, The specific method for S5 to determine whether the updated autonomous orbit determination filter has converged is as follows: if one or more of the following convergence criteria are met, the autonomous orbit determination filter is considered to have diverged and needs to be re-initialized. The specific criteria are as follows: Judgment Condition 1: Substitute the position information in the filtered results after time update and measurement update into the GNSS observation equation to calculate the residual values ​​of all navigation star observations participating in orbit determination filtering. Sum the residual values ​​as the judgment quantity. When the judgment quantity continuously exceeds the set judgment threshold, it is considered that the orbit determination filter is diverging. Decision condition 2: Sum the covariance values ​​of the position vectors in the state covariance matrix before and after the autonomous filter update, denoted as the empirical covariance sum σ. apr and post-test variance and σ pos First, determine the posterior variance and σ. pos If the threshold is exceeded multiple times consecutively, the orbital filter is considered to be diverging; then, σ... apr and σ pos The difference is used as a judgment quantity. When the difference exceeds the set threshold multiple times in a row, it is also considered that the orbit determination filter is diverging. Judgment condition 3: If the difference between the position obtained by geometric positioning at the current time in step S2 and the position obtained by autonomous orbit determination filtering exceeds the set threshold multiple times consecutively, it is considered that the orbit determination filter is diverging.

Images