Orbit maneuver compensation method based on constant acceleration and zero and velocity pulse constraints

Abstract

The application discloses an orbit maneuver compensation method based on constant acceleration and zero-sum speed pulse constraint, relates to the technical field of satellite orbit calculation, and comprises the following steps: acquiring GNSS observation data, satellite attitude data, ephemeris data and external source data, and performing orbit fitting to obtain a rough dynamic orbit; performing cycle slip detection and gross error elimination on the GNSS observation data to obtain clean observation values required for subsequent parameter estimation; using the clean observation values, satellite state, state transition matrix and parameter sensitivity matrix, adding a zero-sum speed pulse constraint equation to an original observation error equation to construct an augmented observation error equation; iteratively performing parameter estimation and orbit integration until convergence is achieved; and generating a high-precision dynamic orbit. The orbit maneuver compensation method based on constant acceleration and zero-sum speed pulse constraint improves the consistency between the dynamic orbit and the kinematic orbit and the KBR checking precision of the dynamic baseline.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of satellite orbit calculation, and particularly relates to an orbit maneuver compensation method based on constant acceleration and zero-sum velocity pulse constraints. BACKGROUND

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

[0003] At present, in the field of low-orbit satellite orbit determination, the parameter modeling methods of orbit maneuver mainly include the velocity pulse model and the constant acceleration model.

[0004] The traditional velocity pulse model absorbs the error of the dynamic model by adding an instantaneous velocity pulse parameter at a certain time or at several times during the orbit maneuver. This compensation method is only suitable for short-duration orbit maneuvers. If the prior constraints of the velocity pulse parameter are too strong, the condition number of the normal equation matrix will be too large, resulting in a pathological normal equation matrix. If the prior constraints are too weak, the parameters will not be de-correlated, which will also result in a pathological normal equation matrix. The pathological matrix is very sensitive to numerical errors, which will amplify the calculation errors of the inverse of the normal equation matrix, causing a serious deviation from the true value, and thus leading to the failure of orbit determination and the inability to generate a dynamic orbit. The constant acceleration model uses an average compensation strategy, which estimates the orbit maneuver acceleration as a constant during the satellite orbit maneuver. Compared with the velocity pulse model, this compensation model is not only suitable for short-duration orbit maneuvers, but also suitable for long-duration orbit maneuvers, and the compensation is more complete. However, the compensation effect is still not good for very long-duration orbit maneuvers.

[0005] Therefore, from the accuracy and robustness of orbit determination, it is necessary to design an orbit maneuver compensation method that has better compensation effect and does not lead to the failure of orbit determination. SUMMARY

[0006] To solve the above problems, the present application proposes an orbit maneuver compensation method based on constant acceleration and zero-sum velocity pulse constraints. In view of the defects and improvement needs of the existing orbit maneuver compensation methods, and considering different orbit maneuver durations and the rank deficiency problem of the normal equation matrix, the present application is based on least squares batch processing and uses the measured data under the orbit maneuver conditions of GRACE-FO to verify the new compensation method, aiming to solve the technical problem that the existing orbit maneuver compensation methods do not have robustness and accuracy.

[0007] The orbit maneuver compensation method based on constant acceleration and zero-sum velocity pulse constraint comprises the following steps:

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

[0009] S2, detecting cycle slip and eliminating gross errors of the GNSS observation data to obtain clean observation values required for subsequent parameter estimation;

[0010] S3, using the clean observation values, satellite state, state transition matrix and parameter sensitivity matrix, adding a zero-sum velocity pulse constraint equation to the original observation error equation to construct an augmented observation error equation;

[0011] S4, iteratively performing parameter estimation and orbit integration until convergence to generate a high-precision dynamic orbit.

[0012] Preferably, the specific content of S1 of acquiring GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain a rough dynamic orbit is as follows:

[0013] S101, acquiring GNSS observation data, satellite attitude data, ephemeris data and external data and performing preprocessing 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 first-order differential equations of satellite motion state, state transition matrix and parameter sensitivity matrix according to the satellite initial state parameters, dynamic parameters and perturbation force model, merging the first-order differential equations of the state transition matrix and the parameter sensitivity matrix to obtain variational equations, and jointly solving the first-order differential equations of the satellite motion state and the variational equations by a numerical integration method to obtain the satellite state, the state transition matrix and the parameter sensitivity matrix at any time, wherein the dynamic parameters include constant acceleration parameters and velocity pulse parameters;

[0016] S104, taking a satellite single point positioning orbit as a reference orbit for orbit fitting, constructing an original observation error equation using the satellite state, the state transition matrix and the parameter sensitivity matrix, and then solving the equation by a least square method to obtain satellite initial state parameters and dynamic parameters;

[0017] S105, iteratively performing steps S103 and S104 until a rough dynamic orbit is fitted.

[0018] Preferably, the specific content of S3 of adding a zero-sum velocity pulse constraint equation to the original observation error equation to construct an augmented observation error equation is as follows:

[0019] Increase the velocity pulse parameter ΔV to absorb the residual perturbation model error and impose a zero-sum constraint 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 expressions of the original observation error equation and the normal equation are respectively:

[0022] v l =Adx-l, 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 velocity pulse parameters is:

[0027] It is defined that s sets 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 velocity pulse parameter is: the zero-sum velocity pulse constraint equation is:

[0031] ε=Ddx,W ΔV ;

[0032] Where ε is the zero-sum constraint residual vector, W ΔV The zero-sum constraint weight matrix is ​​defined as follows: if a total of m orbital maneuvers occur in the orbit-fixing arc, 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 elements 1 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 to obtain a high-precision dynamical orbit, is as follows:

[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 follows:

[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 the 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 advantages over traditional technologies:

[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 parameter and the velocity pulse parameter during the parameter estimation process, and estimate these two parameters at the same time, thus avoiding 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 adopting the compensation method based on constant acceleration and zero and velocity pulse constraints in parameter estimation, the consistency between the dynamic orbit and the kinematic orbit has basically been improved to varying degrees. The improvement effect is particularly obvious for arc segments with very long orbital maneuvers, and 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 accompanying drawings and examples. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0055] Figure 2 This is a comparison chart of the RMS difference between the dynamical 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 chart of the RMS difference between the dynamical 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 comparison diagram of the dynamic baseline KBR check is obtained by using the present invention and the constant acceleration compensation method to determine the orbit. DETAILED DESCRIPTION

[0058] The technical method of the present invention is further described below through the accompanying drawings and embodiments. It should be noted that unless otherwise specifically stated, the relative arrangement of components and steps, numerical expressions and values ​​described in these embodiments do not limit the scope of this 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 disclosure, its application, or uses.

[0060] Technologies, systems, and devices known to those of ordinary skill in the relevant art may not be discussed in detail, but where appropriate, they 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 same meaning as commonly understood by one of ordinary skill in the art 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 to the error equation to more accurately compensate for the dynamic model error caused by orbital maneuvers.

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

[0065] S1. Obtain GNSS observation data, satellite attitude data, ephemeris data and external data and perform 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 the 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 files used are L1A-level data of the GRACE-FO satellite, and the attitude files are L1B-level data. The GRACE-FO orbital maneuver information for 2018-2021 is read from the official TN-01a_SCE.txt file, which contains a total of 8 orbital maneuver days. Each orbit determination arc is from 0:00 to 24:00 on the orbital maneuver day. For ease of description, the GRACE-FO C satellite is referred to as GFC, and the GRACE-FO D satellite is referred to as GFD. The orbital maneuver times for GRACE-FO from 2018 to 2021 are shown in the following table:

[0069] Table 1 GFC orbital 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 track maneuver 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, generating a satellite discrete orbit according to the data required for orbiting, obtaining satellite initial state parameters according to the satellite discrete orbit.

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

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

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

[0077]

[0078] In the formula, t is time, is the satellite acceleration vector, which can be calculated by the satellite motion equation, that is:

[0079]

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

[0081] Specifically, the state transition matrix is The first-order differential equation thereof is:

[0082]

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

[0084] The parameter sensitive matrix is The first-order differential equation thereof is:

[0085]

[0086] In the formula, n p is the number of dynamic parameters, and since the satellite initial state is irrelevant to the dynamic parameters, the initial condition is:

[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] Introduction:

[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, constructing an original observation error equation using the satellite state, state transition matrix, and parameter sensitivity matrix, and then solving the equation using the least squares method to obtain the satellite initial state parameters and dynamic parameters;

[0097] S105 , iterate 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), state transfer matrix Φ(t,t0), and parameter sensitivity matrix S(t) can all be obtained through numerical integration methods.

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

[0102] S3. Using clean observations, satellite states, state transition matrices, and parameter sensitivity matrices, the zero-sum velocity pulse constraint equation is added to the original observation error equation to construct an augmented observation error equation.

[0103] S4. Iteratively perform parameter estimation and orbit integration until convergence, and generate high-precision dynamic orbit.

[0104] Specifically, the estimation strategy of constant acceleration and zero-sum constrained velocity pulse parameters in S4 is as follows:

[0105] (1) For constant acceleration parameters, a set of (tangential, radial, normal) is estimated each time of orbit maneuver.

[0106] (2) For zero-sum constrained velocity pulse parameters, if the duration of orbit maneuver t is greater than 300 seconds, divide the orbit maneuver time by 300 and take the integer part to get k, set the interval of velocity pulse parameters τ=t / k, and estimate a set of velocity pulse parameters (tangential, radial, normal) at the end of each interval. If the duration of orbit maneuver is not greater than 300 seconds, estimate a set of velocity pulse parameters at the middle and end time respectively.

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

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

[0109] In the formula, Y represents the observation value, h(T,Q,B) is the model calculation value, T is the receiver clock bias parameter, B is the ambiguity parameter, Q includes the dynamic parameters p and the satellite initial state parameters y(t0), the constant acceleration parameters a man and the zero-sum constrained velocity pulse parameters ΔV all belong to the dynamic parameters p, and v l is the observation residual vector.

[0110] After Taylor expansion at the parameter initial value (T*,Q*,B*) and ignoring terms higher than the second order, we get:

[0111]

[0112] Let dT=(T-T*), dQ=(Q-Q*), and dB=(B-B*).

[0113] Let h i (T,Q,B) be the model calculation value at time t i , and Y i be the observation vector, then the observation error equation can be obtained as:

[0114]

[0115] In the formula, (dT,dQ,dB) is the correction of the initial value of the estimated parameters.

[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- l, W;

[0120]

[0121] where A is the design matrix, dx is the parameter correction vector, l is the observation minus the calculated value, W is the observation weight matrix, is the prior information weight matrix of dx, v l is the observation residual, is the parameter correction residual.

[0122] Specifically, considering that the constant acceleration parameter a man in the to-be-estimated dynamic parameters p has a strong correlation with the velocity pulse parameter ΔV:

[0123]

[0124] The shorter the orbit maneuver time, the more linearly correlated the constant acceleration parameter and the velocity pulse parameter are, and in order to decorrelate the two parameters, a zero-sum constraint is imposed on the velocity pulse parameter.

[0125] Defining that s groups of velocity pulses are estimated in the jth orbit maneuver, the zero-sum constraint imposed on the velocity pulse parameter is:

[0126]

[0127] where A, C and R are tangential, normal and radial respectively.

[0128] The zero-sum velocity pulse constraint equation is:

[0129] ε = Ddx, W ΔV ;

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

[0131]

[0132] Where s is the number of the velocity impulse set estimated by the orbit maneuver, which is equal to the number of 1 in any row. The three rows correspond to the zero and constraint of A, C and R respectively.

[0133] Considering the zero-sum velocity impulse constraint equation, the augmented error equation can be written as:

[0134]

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

[0136]

[0137] Solving the augmented equation gives:

[0138]

[0139] The initial value of the next parameter estimation:

[0140] x = x0 + dx

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

[0142]

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

[0144] The quadratic form of the post-fit residual is:

[0145]

[0146] The post-fit unit weight error is:

[0147]

[0148] Where n l , n x , n ΔV are the number of observations, the total number of parameters, and the number of added zero-sum constraint velocity impulses respectively. Generally, it is considered that Close to 1 indicates that the weight is reasonable. After several parameter estimations, until dx converges.

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

[0150] After the ambiguity is fixed, a parameter estimation and orbit integration are performed again to obtain a single-star dynamic orbit, and a dynamic baseline is obtained by differencing the double-star dynamic orbits, and a kinematic orbit is obtained according to the residual file generated by parameter estimation.

[0151] The dynamic orbit and kinematic orbit difference RMS of a plurality of arc segments is counted, the generated dynamic baseline is subjected to KBR checking, and the improvement effect of the application compared with the constant acceleration compensation method is compared.

[0152] Specifically, as shown in Figure 2 and Figure 3 The light purple bar represents the dynamic and kinematic orbit difference RMS generated by using the constant acceleration compensation method, and the red bar represents the dynamic and kinematic orbit difference RMS generated by using the constant acceleration and zero-sum velocity pulse constraint orbit maneuver compensation method. It can be seen that for all arc segments, the constant acceleration and zero-sum velocity pulse constraint orbit maneuver compensation method is better than the constant acceleration compensation method, especially for the 2021 year accumulation day 307, the GFC and the GFD respectively occur once for nearly 1800s. The application improves the consistency of the dynamic and kinematic orbits by 73.7%. As shown in Figure 4 The KBR checking comparison chart of the dynamic baseline generated by the two different orbit maneuver compensation methods, also for the 2021 year accumulation day 307, the application improves the KBR checking accuracy of the dynamic baseline by 67.1%.

[0153] The application considers that during the actual orbit maneuver, the orbit maneuver force will not be constant, so the 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, the zero-sum constraint is applied to the velocity pulse parameter in the parameter estimation part, which is de-correlated with the constant acceleration parameter, ensuring the success of orbit determination, and has robustness and accuracy. Finally, through the GRACE-FO measured data verification, compared with the dynamic orbit and baseline generated by using the constant acceleration compensation model, the consistency between the dynamic orbit and the kinematic orbit generated by using the application is significantly improved, and the KBR checking accuracy of the dynamic baseline is also improved, which proves that compared with the traditional orbit maneuver compensation method, the effect of the application is better.

[0154] Those skilled in the art will readily understand that the above is only a preferred embodiment of the application and does not limit the application, and any modifications, equivalent replacements and improvements made within the spirit and principles of the application shall be included in the protection scope of the application.

[0155] It should be pointed out finally that the above embodiments are only used to illustrate the technical method of the present application but not to limit it, and although the present application is described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical method of the present application can be modified or replaced equivalently, and these modifications or equivalent replacements cannot make the modified technical method deviate from the spirit and scope of the technical method of the present application.

Claims

1. Orbital maneuver compensation method based on constant acceleration and zero-sum velocity pulse constraints, characterized by: The following steps are involved: S1. Obtain GNSS observation data, satellite attitude data, ephemeris data and external data and perform orbit fitting to obtain a rough dynamic orbit; S2. Perform cycle slip detection and gross error elimination on the GNSS observation data to obtain clean observation values ​​required for subsequent parameter estimation; S3. Using clean observations, satellite states, state transition matrices, and parameter sensitivity matrices, the zero-sum velocity pulse constraint equation is added 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 acquiring GNSS observation data, satellite attitude data, ephemeris data and external data and performing orbit fitting to obtain the rough dynamic orbit are as follows: 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. Generate satellite discrete orbits based on the data required for orbit determination, and obtain satellite initial state parameters based on the satellite discrete orbits; S103. Constructing first-order differential equations of the satellite's motion state, state transfer matrix, and parameter sensitivity matrix based on the satellite's initial state parameters, dynamic parameters, and perturbation force model; combining the first-order differential equations of the state transfer matrix and the parameter sensitivity matrix to obtain a variational equation; jointly solving the first-order differential equation and the variational equation of the satellite's motion state using a numerical integration method to obtain the satellite state, state transfer matrix, and parameter sensitivity matrix at any time; wherein the dynamic parameters include a constant acceleration parameter and a velocity pulse parameter; S104, using the satellite single point positioning orbit as the reference orbit for orbit fitting, constructing an original observation error equation using the satellite state, state transition matrix, and parameter sensitivity matrix, and then solving the equation using the least squares method to obtain the satellite initial state parameters and dynamic parameters; S105 , iterate 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: Increase speed pulse parameter This is used to absorb the residual perturbation force model error and impose zero-sum constraints on the velocity pulse parameters; 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: ； ； ； in, is the design matrix, is the transposed matrix of the design matrix, is the parameter correction vector to be estimated, is the observed value minus the calculated value, is the observation weight matrix, for The prior information weight matrix, 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 4, characterized in that: The specific contents of imposing zero-sum constraints on velocity pulse parameters are as follows: Definition Sub-orbital maneuvers estimated If there is a group of velocity pulses, 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 5, characterized in that: The specific content of applying zero-sum constraint to velocity pulse parameters is as follows: The zero-sum velocity pulse constraint equation is: ； Where, is the zero-sum constrained residual vector, is the zero-sum constraint weight matrix, the matrix The definition is: if a total of Sub-orbital maneuvers, then the matrix That is OK A matrix of columns, is the total number of parameters to be estimated, and a certain orbital maneuver corresponds to the matrix The elements in are: ； Among them, the number of velocity pulse groups estimated for this orbital maneuver is It 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 6, 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. Iterate parameter estimation and orbital integration until the parameters converge to obtain high-precision dynamical orbits. 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: ； in, is the parameter initial value vector. According to the covariance propagation law, the parameter error covariance matrix after update is: ； in, , , 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, is the covariance matrix of the parameter correction; Perform multiple parameter estimations until the parameters converge.

Images