Orbit maneuver compensation method based on constant acceleration and zero sum speed pulse constraint

Abstract

The invention discloses an orbital maneuver compensation method based on constant acceleration and zero sum velocity pulse constraints, and relates to the technical field of satellite orbit calculation, comprising the steps of acquiring GNSS observation data, satellite attitude data, ephemeris data and exogenous data, performing orbit fitting to obtain a rough dynamic orbit, performing cycle slip detection and gross error elimination on the GNSS observation data, and performing orbit maneuver compensation on the rough dynamic orbit. Obtaining a clean observation value required by subsequent parameter estimation, adding a zero sum speed pulse constraint equation on the basis of an original observation error equation by using the clean observation value, a satellite state, a state transition matrix and a parameter sensitive matrix, constructing an augmented observation error equation, and iteratively performing parameter estimation and orbit integration until convergence, thereby obtaining a final observation value. According to the method, the orbit maneuvering compensation method based on the constant acceleration and the zero sum velocity pulse constraint is adopted, and the consistency between the dynamic orbit and the kinematic orbit and the KBR checking precision of a dynamic baseline are improved.

Description

Technical Field

[0001] The invention relates to the technical field of satellite orbit calculation, in particular to an orbit maneuver compensation method based on constant acceleration and zero-sum velocity pulse constraints. Background Art

[0002] Satellite orbital maneuvers refer to the process in which a spacecraft uses its own propulsion device to exert thrust (chemical propulsion, electric propulsion, cold gas propulsion, etc.) to change orbital parameters. This process can achieve orbit adjustment, orbit maintenance, and other specific mission objectives of the spacecraft. The purpose of orbital maneuvers is to overcome the impact of orbital perturbations (such as the non-spherical gravity of the earth, atmospheric drag, solar radiation pressure, etc.) on the orbit, or to adapt to mission requirements (such as orbit transfer or rendezvous and docking).

[0003] Currently, in the field of low-orbit satellite orbit determination, the parameter modeling methods for orbital maneuvers mainly include velocity pulse model and constant acceleration model.

[0004] The traditional velocity pulse model absorbs the error of the dynamic model by adding instantaneous velocity pulse parameters at one or several moments during the orbital maneuver. This compensation method is only applicable to orbital maneuvers of short duration. If the prior constraints of the velocity pulse parameters are too strong, the condition number of the normal matrix will be too large, resulting in the ill-conditioned normal matrix. If the prior constraints are too weak, the parameters cannot be decorrelated, which will also cause the normal matrix to be ill-conditioned. The ill-conditioned matrix is ​​very sensitive to numerical errors, which will cause the calculation error of the inversion of the normal matrix to be amplified and seriously deviate from the true value, thus causing orbit determination failure and failure to generate dynamic orbits. The constant acceleration model adopts an average compensation strategy, estimating the orbital maneuver acceleration as a constant during the satellite orbital maneuver. Compared with the velocity pulse model, this compensation model is not only applicable to short-term orbital maneuvers, but also to orbital maneuvers of longer duration, and the compensation is more adequate. However, the compensation effect is still poor for orbital maneuvers of particularly long duration.

[0005] Therefore, from the perspective of orbit determination accuracy and robustness, it is very necessary to design an orbit maneuver compensation method that has a better orbit maneuver compensation effect and will not cause orbit determination failure. Summary of the invention

[0006] To solve the above problems, this application proposes an orbital maneuvering compensation method based on constant acceleration and zero-sum velocity pulse constraints. In response to the defects and improvement needs of the existing orbital maneuvering compensation methods, taking into account the different durations of orbital maneuvers and the rank deficiency problem of the normal equation matrix, based on least squares batch processing and using the measured data under GRACE-FO orbital maneuvering conditions, the new compensation method is verified, aiming to solve the technical problem that the existing orbital maneuvering compensation methods do not have both robustness and accuracy.

[0007] The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints includes the following steps:

[0008] S1, obtaining GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain a rough dynamic orbit;

[0009] S2, perform cycle slip detection and gross error elimination on GNSS observation data to obtain clean observation values ​​required for subsequent parameter estimation;

[0010] S3, using clean observations, satellite states, state transfer matrices, and parameter sensitivity matrices, add zero and velocity pulse constraint equations to the original observation error equation to construct an augmented observation error equation;

[0011] S4. Iterate parameter estimation and orbital integration until convergence to generate high-precision dynamical orbits.

[0012] Preferably, the specific contents of acquiring GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain a rough dynamic orbit in S1 are:

[0013] S101, obtaining GNSS observation data, satellite attitude data, ephemeris data and external source data and preprocessing them to obtain data required for orbit determination;

[0014] S102, generating a satellite discrete orbit according to the data required for orbit determination, and obtaining satellite initial state parameters according to the satellite discrete orbit;

[0015] S103, constructing the first-order differential equations of the satellite motion state, state transfer matrix and parameter sensitivity matrix according to the initial state parameters, dynamic parameters and perturbation force model of the satellite, combining the first-order differential equations of the state transfer matrix and the parameter sensitivity matrix to obtain the variational equation, jointly solving the first-order differential equation and the variational equation of the satellite motion state by a numerical integration method, and obtaining the satellite state, state transfer matrix and parameter sensitivity matrix at any time, wherein the dynamic parameters include constant acceleration parameters and velocity pulse parameters;

[0016] S104, using the satellite single point positioning orbit as the reference orbit for orbit fitting, using the satellite state, state transfer matrix, and parameter sensitivity matrix to construct an original observation error equation, and then using the least squares method to solve the equation to obtain the satellite initial state parameters and dynamic parameters;

[0017] S105, iterating steps S103 and S104 until a rough dynamic trajectory is fitted.

[0018] Preferably, in S3, the zero and velocity pulse constraint equations are added to the original observation error equation, and the specific content of constructing the augmented observation error equation is:

[0019] The velocity pulse parameter ΔV is increased to absorb the residual perturbation force model error and a zero-sum constraint is imposed on the velocity pulse parameter;

[0020] The zero and velocity pulse constraint equations are added to the original observation error equation to obtain the augmented observation error equation.

[0021] Preferably, the expression of the original observation error equation and the expression of the normal equation are respectively:

[0022] v l =Adx-1, W;

[0023] v x0 =dx,W x0 ;

[0024] (A T WA+W x0 )dx=A T Wl;

[0025] Where A is the design matrix, dx is the vector of parameter corrections to be estimated, l is the observed value minus the calculated value, W is the observed value weight matrix, and W x0 is the prior information weight matrix of dx;

[0026] Preferably, the specific content of applying zero-sum constraint to the speed pulse parameter is:

[0027] It is defined that s groups of velocity pulses are estimated in the jth orbital maneuver, and the zero-sum constraint imposed on the velocity pulse parameters is:

[0028]

[0029] Among them, A, C, and R are tangential, normal, and radial directions, respectively;

[0030] Preferably, the specific content of applying zero-sum constraint to the speed pulse parameter is: the zero-sum speed pulse constraint equation is:

[0031] ε=Ddx,W ΔV ;

[0032] Where ε is the zero-constrained residual vector, W ΔV is the zero-sum constraint weight matrix, and the matrix D is defined as follows: if a total of m orbital maneuvers occur in the orbit determination arc segment, then the matrix D is 3m rows and n x A matrix of columns, n x is the total number of parameters to be estimated, and the element in the matrix D corresponding to a certain orbital maneuver is:

[0033]

[0034] Among them, the number of velocity pulse groups s estimated for this orbital maneuver is equal to the number of element 1s in any row, and the three rows correspond to the zero-sum constraints in the three directions of A, C, and R, respectively.

[0035] Preferably, the augmented observation error equation is:

[0036]

[0037] According to the weighted least squares principle The augmented method equation is:

[0038]

[0039] Preferably, S4, iteratively performing parameter estimation and orbital integration until the parameters converge, and obtaining a high-precision dynamical orbit, the specific content is:

[0040] Perform parameter estimation and solve the augmented method equation to obtain:

[0041]

[0042] Preferably, the parameter estimation process is performed several times, and the initial value of the next parameter estimation is defined as:

[0043] x=x0+dx;

[0044] Among them, x0 is the initial value vector of the parameter. According to the covariance propagation law, the covariance matrix of the parameter error after update is:

[0045]

[0046] Where N1 = A T WA, N2=D T W ΔV D.

[0047] Perform multiple parameter estimations until the parameters converge.

[0048] In summary, the orbital maneuvering compensation method based on constant acceleration and zero-sum velocity pulse constraints provided by the present invention has the following beneficial effects compared with the traditional technology:

[0049] (1) Compared with the traditional constant acceleration compensation method, the orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints more realistically reflects the physical characteristics of orbital maneuvers;

[0050] (2) The velocity pulse model is suitable for orbital maneuvers with a very short duration. The constant acceleration compensation method can more accurately compensate for orbital maneuvers with a longer duration, but it is not effective for orbital maneuvers with a particularly long duration. The compensation method based on constant acceleration and zero-sum velocity pulse constraints is more universal.

[0051] (3) The compensation method based on constant acceleration and zero and velocity pulse constraints can decorrelate the constant acceleration parameters and velocity pulse parameters during the parameter estimation process, estimate these two parameters at the same time, and avoid the rank deficiency of the normal equation matrix;

[0052] (4) Based on the measured data of GRACE-FO, compared with the traditional constant acceleration compensation method, after the compensation method based on constant acceleration and zero and velocity pulse constraints was adopted in parameter estimation, the consistency between the dynamic orbit and the kinematic orbit was basically improved to varying degrees. The improvement effect is particularly obvious for arc segments with very long orbital maneuvers. The KBR verification results also show that the baseline accuracy has been improved.

[0053] The technical method of the present invention is further described in detail below through the drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 A flow chart of an embodiment of an orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints;

[0055] Figure 2 This is a comparison diagram of the RMS difference between the dynamic orbit and the kinematic orbit of the GFACE-FO C star obtained by orbit determination using the present invention and the constant acceleration compensation method;

[0056] Figure 3 This is a comparison diagram of the RMS difference between the dynamic orbit and the kinematic orbit of the GFACE-FO D star obtained by orbit determination using the present invention and the constant acceleration compensation method;

[0057] Figure 4 A dynamic baseline KBR check comparison diagram is obtained for orbit determination using the present invention and the constant acceleration compensation method. DETAILED DESCRIPTION

[0058] The technical method of the present invention is further described below by means of the accompanying drawings and embodiments. It should be noted that unless otherwise specifically stated, the relative arrangement of the components and steps, numerical expressions and numerical values ​​described in these embodiments do not limit the scope of the present application.

[0059] The following description of at least one exemplary embodiment is merely illustrative in nature and is in no way intended to limit the present application, its application, or uses.

[0060] Technologies, systems, and devices known to ordinary technicians in the relevant art may not be discussed in detail, but where appropriate, the technologies, systems, and devices should be considered part of the specification.

[0061] In all examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not limiting. Therefore, other examples of the exemplary embodiments may have different values.

[0062] Unless otherwise defined, technical or scientific terms used in the present invention shall have the common meanings understood by one having ordinary skills in the field to which the present invention belongs.

[0063] The present invention provides an orbital maneuver compensation method based on constant acceleration and zero-sum velocity pulse constraints. The method adds the zero-sum velocity pulse constraint equation into the error equation to more accurately compensate for the dynamic model error caused by the orbital maneuver.

[0064] This embodiment provides an orbital maneuver compensation method based on constant acceleration and zero and velocity pulse constraints, comprising the following steps:

[0065] S1, obtaining GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain a rough dynamic orbit;

[0066] Furthermore, the specific contents of obtaining GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain a rough dynamic orbit in S1 are as follows:

[0067] S101, obtaining GNSS observation data, satellite attitude data, ephemeris data and external source data and preprocessing them to obtain data required for orbit determination;

[0068] Specifically, the navigation satellite orbit, clock error and bias products used are CODE post-processing products No. 3, the GNSS observation value file used is the L1A data of the GRACE-FO satellite, and the attitude file is the L1B data. The GRACE-FO orbital maneuver information for 2018-2021 is read from the official TN-01a_SCE.txt file, for a total of 8 orbital maneuver days, and each orbit determination arc is from 0:00 to 24:00 on the orbital maneuver day. For the convenience of description, the GRACE-FO C star is referred to as GFC, the GRACE-FO D star is referred to as GFD, and the orbital maneuver time of GRACE-FO from 2018 to 2021 is shown in the following table:

[0069] Table 1 GFC orbit maneuvering time

[0070] Serial number Start time Duration (seconds) 1 2018-05-24 06:19:44.60 127.00 2 2018-05-25 23:23:01.60 416.00 3 2018-05-26 02:32:58.60 423.00 4 2018-06-22 19:57:34.60 49.00 5 2018-10-22 05:03:40.59 735.00 6 2021-11-03 04:08:18.60 1790.00 7 2021-11-03 16:19:03.60 31.00

[0071] Table 2 GFD orbit maneuvering time

[0072] Serial number Start time Duration (seconds) 1 2018-05-24 10:56:15.60 126.00 2 2018-05-26 00:58:06.60 406.00 3 2018-05-26 04:08:02.60 413.00 4 2018-06-22 21:37:42.58 84.00 5 2019-03-19 13:30:22.58 17.00 6 2020-09-09 08:15:15.60 5.00 7 2021-11-03 04:02:48.60 1790.00

[0073] S102: Generate satellite discrete orbits according to data required for orbit determination, and obtain satellite initial state parameters according to the satellite discrete orbits.

[0074] S103, constructing the first-order differential equations of the satellite motion state, state transfer matrix and parameter sensitivity matrix according to the initial state parameters, dynamic parameters and perturbation force model of the satellite, combining the first-order differential equations of the state transfer matrix and the parameter sensitivity matrix to obtain the variational equation, jointly solving the first-order differential equation and the variational equation of the satellite motion state by a numerical integration method, and obtaining the satellite state, state transfer matrix and parameter sensitivity matrix at any time, wherein the dynamic parameters include constant acceleration parameters and velocity pulse parameters;

[0075] Specifically, the satellite state y(t) includes the position vector r(t) and the velocity vector

[0076] The first-order differential equation of the satellite state y(t) can be expressed as:

[0077]

[0078] Where t is time, is the satellite acceleration vector, which can be calculated by the satellite motion equation, namely:

[0079]

[0080] Where p is the dynamic parameter, constant acceleration parameter a man The velocity pulse parameter ΔV with zero and constraint is part of the dynamic parameter p.

[0081] Specifically, the state transfer matrix is Its first-order differential equation is:

[0082]

[0083] Where I is the unit matrix and the initial condition is Φ(t0,t0)=I 6×6 , t0 is the initial time.

[0084] The parameter sensitivity matrix is Its first-order differential equation is:

[0085]

[0086] Where n p is the number of dynamical parameters. Since the initial state of the satellite is independent of the dynamical parameters, the initial conditions are:

[0087] Specifically, the variational equation can be obtained by combining the first-order differential equation of the state transfer matrix and the parameter sensitivity matrix of claim 3:

[0088]

[0089] Import:

[0090]

[0091] After sorting:

[0092]

[0093] The above formula is a variational equation, and its initial conditions are:

[0094]

[0095] By jointly solving the equations of motion and variational equations using numerical integration methods, the satellite state, state transfer matrix and parameter sensitivity matrix at any time can be obtained.

[0096] S104, using the satellite single point positioning orbit as the reference orbit for orbit fitting, using the satellite state, state transfer matrix, and parameter sensitivity matrix to construct an original observation error equation, and then using the least squares method to solve the equation to obtain the satellite initial state parameters and dynamic parameters;

[0097] S105, iterating steps S103 and S104 until a rough dynamic trajectory is fitted.

[0098] Specifically, when establishing the observation error equation, the partial derivative of the satellite observation model h(T,Q,B) with respect to the satellite initial state y(t0) and the dynamic parameter p can be obtained by the following calculation formula:

[0099]

[0100] Among them, the satellite state y(t), the state transfer matrix Φ(t,t0), and the parameter sensitivity matrix S(t) can all be obtained by numerical integration methods.

[0101] S2, perform cycle slip detection and gross error elimination on GNSS observation data to obtain clean observation values ​​for subsequent parameter estimation process;

[0102] S3, using clean observations, satellite states, state transfer matrices, and parameter sensitivity matrices, add zero and velocity pulse constraint equations to the original observation error equation to construct an augmented observation error equation;

[0103] S4. Iterate parameter estimation and orbital integration until convergence to generate high-precision dynamical orbits.

[0104] Specifically, the estimation strategy of the constant acceleration and zero-sum-constrained velocity pulse parameters used for track maneuver force compensation in S4 is:

[0105] (1) For constant acceleration parameters, one set (tangential, radial, and normal) is estimated for each orbital maneuver.

[0106] (2) For the velocity pulse parameters with zero-sum constraints, if the orbital maneuver duration t is greater than 300 seconds, the orbital maneuver duration is divided by 300 and rounded to get k, and the interval of the velocity pulse parameters is set to τ = t / k. A ​​set of velocity pulse parameters (tangential, radial, and normal) is estimated at the end of each interval. If the orbital maneuver duration is not greater than 300 seconds, a set of velocity pulse parameters is estimated at the middle and end times, respectively.

[0107] Specifically, the general expression of the observation equation is:

[0108] Y=h(T,Q,B)-v l ;

[0109] Where Y represents the observed value, h(T,Q,B) is the model calculated value, T is the receiver clock error parameter, B is the ambiguity parameter, Q includes the dynamic parameter p and the satellite initial state parameter y(t0), and the constant acceleration parameter a man and the zero-constrained velocity pulse parameter ΔV are all part of the dynamic parameters p, v l is the observation residual vector.

[0110] Taylor expansion is performed at the initial values ​​of the parameters (T*, Q*, B*), and the terms above the second order are ignored to obtain:

[0111]

[0112] Let dT=(TT*), dQ=(QQ*), dB=(BB*).

[0113] Note i The observation model calculation value at time h i (T,Q,B), the observation vector is Y i , then the observation error equation is:

[0114]

[0115] Where (dT, dQ, dB) is the correction of the initial value of the parameter to be estimated.

[0116] Let the observation weight matrix be W, then the normal equation is:

[0117]

[0118] The original error equation and normal equation can be written as:

[0119] v l =Adx-1, W;

[0120]

[0121] Among them, A is the design matrix, dx is the vector of parameter corrections to be estimated, l is the observed value minus the calculated value, and W is the observed value weight matrix. is the prior information weight matrix of dx, v l is the observed residual, is the residual of the correction value of the parameter to be estimated.

[0122] Specifically, considering the constant acceleration parameter a in the dynamic parameter p to be estimated man There is a strong correlation with the speed pulse parameter ΔV:

[0123]

[0124] The shorter the orbital maneuvering time, the closer the constant acceleration parameter is to the velocity pulse parameter in linear correlation. In order to decorrelate these two parameters, we need to impose a zero-sum constraint on the velocity pulse parameter.

[0125] It is defined that s groups of velocity pulses are estimated in the jth orbital maneuver, and the zero-sum constraint imposed on the velocity pulse parameters is:

[0126]

[0127] Among them, A, C, and R are tangential, normal, and radial directions, respectively.

[0128] The zero and velocity pulse constraint equations are:

[0129] ε=Ddx,W ΔV ;

[0130] Where ε is the zero-constrained residual vector, W ΔV is the zero-sum constraint weight matrix, and the matrix D is defined as follows: if a total of m orbital maneuvers occur in the orbit determination arc segment, then the matrix D is 3m rows and n x A matrix of columns, n x is the total number of parameters to be estimated, and the element in the matrix D corresponding to a certain orbital maneuver is:

[0131]

[0132] Among them, the number of velocity pulse groups s estimated for this orbital maneuver is equal to the number of element 1s in any row, and the three rows correspond to the zero-sum constraints in the three directions of A, C, and R, respectively.

[0133] After considering the zero and velocity pulse constraint equations, the augmented error equation can be written as:

[0134]

[0135] According to the weighted least squares principle The augmented method equation is:

[0136]

[0137] Solving the augmented method equation yields:

[0138]

[0139] Initial value for the next parameter estimate:

[0140] x=x0+dx;

[0141] Among them, x0 is the initial value vector of the parameter. According to the covariance propagation law, the covariance matrix of the parameter error after update is:

[0142]

[0143] Where N1 = A T WA, N2=D T W ΔV D;

[0144] The quadratic form of the posterior residual is:

[0145]

[0146] The post-test unit weighted mean error is:

[0147]

[0148] Among them, n l ,n x ,n ΔV are the number of observations, the total number of parameters, the number of added zeros and the number of constrained velocity pulses. If it is close to 1, it means that the weighting is reasonable. After multiple parameter estimations, dx converges.

[0149] After the parameters converge, the ambiguity is fixed to further improve the orbit determination accuracy.

[0150] After the ambiguity is fixed, parameter estimation and orbit integration are performed again to obtain the single-star dynamic orbit. The dynamic baseline is obtained by subtracting the binary star dynamic orbit, and the kinematic orbit is obtained based on the residual file generated by the parameter estimation.

[0151] The RMS difference between the dynamic trajectory and the kinematic trajectory of multiple arc segments is counted, and the generated dynamic baseline is checked by KBR to compare the improvement effect of the present invention with that of the constant acceleration compensation method.

[0152] Specifically, Figure 2 and Figure 3 As shown, the light purple bar represents the RMS difference between the dynamic and kinematic orbits generated by the constant acceleration compensation method, and the red bar represents the RMS difference between the dynamic and kinematic orbits generated by the orbital maneuver compensation method using constant acceleration and zero-sum velocity pulse constraints. It can be seen that for all arc segments, the orbital maneuver compensation method with constant acceleration and zero-sum velocity pulse constraints is superior to the constant acceleration compensation method, especially for the 2021 annual accumulation day 307, GFC and GFD each had an orbital maneuver lasting nearly 1800s. This invention improved the consistency of the dynamic and kinematic orbits by 73.7%. Figure 4 Shown is a comparison chart of the dynamic baseline KBR verification generated by two different orbital maneuver compensation methods. Similarly, for the annual accumulation day 307 in 2021, the invention improves the dynamic baseline KBR verification accuracy by 67.1%.

[0153] The present invention takes into account that the orbital maneuvering force will not remain constant during the actual orbital maneuver, so a velocity pulse parameter is added on the basis of the constant acceleration model, which is closer to the real physical process. At the same time, a zero-sum constraint is imposed on the velocity pulse parameter in the parameter estimation part, and the velocity pulse parameter is decorrelated with the constant acceleration parameter, which ensures the success of orbit determination and has both robustness and accuracy. Finally, it is verified by GRACE-FO measured data that compared with the dynamic orbit and baseline generated by the constant acceleration compensation model, the consistency between the dynamic orbit and the kinematic orbit generated by the present invention has been significantly improved, and the KBR verification accuracy of the dynamic baseline has also been improved, proving that compared with the traditional orbital maneuver compensation method, the compensation effect of the present invention is better.

[0154] It will be easily understood by those skilled in the art that the above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

[0155] Finally, it should be noted that the above embodiments are only used to illustrate the technical method of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical method of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical method to deviate from the spirit and scope of the technical method of the present invention.

Claims

1. Orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints, characterized in that: The following steps are involved: S1, obtaining GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain a rough dynamic orbit; S2, perform cycle slip detection and gross error elimination on GNSS observation data to obtain clean observation values ​​required for subsequent parameter estimation; S3, using clean observations, satellite states, state transfer matrices, and parameter sensitivity matrices, add zero and velocity pulse constraint equations to the original observation error equation to construct an augmented observation error equation; S4. Iterate parameter estimation and orbital integration until convergence to generate high-precision dynamical orbits.

2. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 1, characterized in that: The specific contents of S1 in which GNSS observation data, satellite attitude data, ephemeris data and external data are obtained and orbit fitting is performed to obtain the rough dynamic orbit are: S101, obtaining GNSS observation data, satellite attitude data, ephemeris data and external source data and preprocessing them to obtain data required for orbit determination; S102, generating a satellite discrete orbit according to the data required for orbit determination, and obtaining satellite initial state parameters according to the satellite discrete orbit; S103, constructing the first-order differential equations of the satellite motion state, state transfer matrix and parameter sensitivity matrix according to the initial state parameters, dynamic parameters and perturbation force model of the satellite, combining the first-order differential equations of the state transfer matrix and the parameter sensitivity matrix to obtain the variational equation, jointly solving the first-order differential equation and the variational equation of the satellite motion state by a numerical integration method, and obtaining the satellite state, state transfer matrix and parameter sensitivity matrix at any time, wherein the dynamic parameters include constant acceleration parameters and velocity pulse parameters; S104, using the satellite single point positioning orbit as the reference orbit for orbit fitting, using the satellite state, state transfer matrix, and parameter sensitivity matrix to construct an original observation error equation, and then using the least squares method to solve the equation to obtain the satellite initial state parameters and dynamic parameters; S105, iterating steps S103 and S104 until a rough dynamic trajectory is fitted.

3. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 1, characterized in that: In S3, the zero and velocity pulse constraint equations are added to the original observation error equation, and the specific content of the augmented observation error equation is constructed as follows: The velocity pulse parameter ΔV is increased to absorb the residual perturbation force model error and a zero-sum constraint is imposed on the velocity pulse parameter; The zero and velocity pulse constraint equations are added to the original observation error equation to obtain the augmented observation error equation.

4. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 3, characterized in that: The expressions of the original observation error equation and the normal equation are respectively: v l =Adx-l,W; Among them, A is the design matrix, A T is the transposed matrix of the design matrix, dx is the vector of parameter corrections to be estimated, l is the observed value minus the calculated value, W is the observed value weight matrix, is the prior information weight matrix of dx, v l is the observed residual, is the residual of the correction value of the parameter to be estimated.

5. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 3, characterized in that: The specific contents of imposing zero-sum constraints on velocity pulse parameters are as follows: It is defined that s groups of velocity pulses are estimated in the jth orbital maneuver, and the zero-sum constraint imposed on the velocity pulse parameters is: Among them, A, C, and R are tangential, normal, and radial directions, respectively.

6. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 3, characterized in that: The specific content of applying zero-sum constraint to speed pulse parameters is: The zero-sum speed pulse constraint equation is: ε=Ddx,W ΔV ; Where ε is the zero-constrained residual vector, W ΔV is the zero-sum constraint weight matrix, and the matrix D is defined as follows: if a total of m orbital maneuvers occur in the orbit determination arc segment, then the matrix D is 3m rows and n x A matrix of columns, n x is the total number of parameters to be estimated, and the element in the matrix D corresponding to a certain orbital maneuver is: Among them, the number of velocity pulse groups s estimated for this orbital maneuver is equal to the number of element 1s in any row, and the three rows correspond to the zero-sum constraints in the three directions of A, C, and R, respectively.

7. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 3, characterized in that: The augmented observation error equation is: According to the weighted least squares principle The augmented method equation is:

8. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 7, characterized in that: S4, iteratively perform parameter estimation and orbital integration until the parameters converge, and obtain the high-precision dynamical orbit. The specific content is: Perform parameter estimation and solve the augmented method equation to obtain:

9. The orbital maneuvering compensation method based on constant acceleration and zero and velocity pulse constraints according to claim 8, characterized in that: The process of parameter estimation is repeated several times, and the initial value of the next parameter estimation is defined: x=x0+dx; Among them, x0 is the initial value vector of the parameter. According to the covariance propagation law, the covariance matrix of the parameter error after update is: Where N1 = A T WA, N2=D T W ΔV D, C x is the covariance matrix of the initial value of the next parameter estimation, is the covariance matrix of the initial value of the current parameter estimate, C dx is the covariance matrix of parameter correction; Perform multiple parameter estimations until the parameters converge.

Images