Method and system for acquiring two-row elements by numerical orbit recursion device

Abstract

The invention provides a method and a system for acquiring two rows of elements by a numerical orbit recursion device. The method comprises the following steps of: inputting a two-row element template; determining orbit time reference points corresponding to the two rows of elements; according to a set numerical value recursion device model, a numerical value recursion device obtains a position value sequence of the track in the set space coordinate system within a period of time; converting the position value sequence obtained by the numerical recursion device into a true equatorial mean equinox coordinate system; obtaining two rows of elements by using a single-point method through the two rows of element templates; taking the two rows of elements obtained by the single-point method as an initial orbit, and solving an orbit parameter estimation value in the two rows of elements under an optimal criterion by utilizing a least square method; calculating an orbit residual error sequence corresponding to an orbit parameter estimation value in the two rows of elements; and fitting the residual error by using a polynomial model to obtain a compensation item. The method is easy to implement, the numerical orbit recursion device can be effectively applied to be converted into two rows of elements, and the method can be compatible with time periods with orbit maneuvering.

Description

Technical Field

[0001] The present invention relates to the technical field of orbit prediction, and specifically, to a method and system for obtaining two-line elements by a numerical orbit propagator. Background Art

[0002] Two-line elements (TLE) is a commonly used international form for publishing orbit elements. It is a data format formulated by the US Department of Defense, and an orbit prediction model SGP4 (Simplified general perturbations) based on two-line elements has been published, which is widely used in the orbit expression and prediction of low Earth orbit satellites.

[0003] With the development of space technology, the number of on-orbit satellites has increased rapidly, including many satellites with orbital maneuvers. Taking low Earth orbit satellites as an example, at an altitude of 500 to 2000 kilometers above the ground, there are nearly ten thousand satellites. For the purpose of orbit maintenance and other goals, many satellites are subject to non-natural forces (i.e., orbital maneuver forces). Most orbital perturbation factors are considered in the traditional orbit prediction model SGP4, but the maneuver forces are not defined, resulting in certain accuracy limitations of the two-line elements and the SGP4 model during the maneuver section of the satellite.

[0004] For the solution of the satellite motion equation, except for the most basic two-body problem, the analytical solutions considering perturbation terms are all approximate solutions under certain conditions. With the development of computer technology, the numerical method for orbit determination and orbit propagation is more and more widely used. Under a given differential equation, a discrete solution with a certain accuracy is used to approximate the exact solution, and it has the advantages of being able to define any form of perturbation force, etc. The numerical propagator is the basic recurrence model for studying orbital motion by the numerical method. In order to be compatible with the commonly used two-line element format, in the field of orbit prediction technology, it is often necessary to convert the numerical propagator into two-line elements. Document 1 (Analysis of TLE Two-Line Ephemeris Data Parameters, Digital Communication World, No. 8, 2020) gives the format and parameter definitions of TLE data. Document 2 (Research on the Orbit Error of the Two-Line Element Cataloging System, Acta Astronomica Sinica, No. 3, 2018) analyzes the influence of the zonal harmonic perturbation term of the Earth's non-spherical gravitational potential in the dynamic model on the accuracy of the cataloged orbit, and finds that there are systematic errors that vary periodically with time in the TLE cataloged orbit.

[0005] In terms of obtaining the two-line element set, Patent Document 1 (Method for generating two-line element set based on model error compensation, CN103514362A) discloses a method for generating two-line element set based on model error compensation. This method first converts the observed position and velocity vectors from the geocentric equatorial inertial coordinate system to the true equatorial mean equinox coordinate system, then calculates the corresponding instantaneous orbital elements and uses them as the initial estimates for the iterative calculation of the two-line element set; uses the initial estimates to predict the position and velocity vectors at the corresponding sampling times respectively, calculates the differences between them and the position and velocity vectors at the actual observed sampling times, and determines whether the termination condition is satisfied; calculates the partial derivative matrix of the simplified universal perturbation orbital prediction function with respect to the two-line element set at the sampling time by the numerical differentiation method; obtains the correction amount of the two-line element set by solving the compensated least squares problem; calculates the new two-line element set, and repeats the iteration until the accuracy requirement is met. Document 3 (Obtaining two-line element set by numerical method, Journal of Wuhan University (Information Science Edition), No. 9, 2009) replaces the differentiation with the difference quotient, avoiding the cumbersome solution of partial derivatives, and gives the method for obtaining two-line element set by numerical method. Document 4 (Research on the conversion method from instantaneous orbital elements to two-line element set, Modern Radar, No. 12, 2012) gives the method for obtaining two-line element set from instantaneous orbital elements by numerical iteration method based on SGP4 and SDP4 models. Document 5 (TLE orbit determination based on simplex method, Chinese Journal of Space Science, No. 6, 2020) proposes a local search algorithm - simplex optimization method to realize TLE orbit determination. To avoid the optimal solution obtained by the constructed initial simplex belonging to the local optimum, the Monte Carlo method is introduced to sample the initial simplex, obtain the statistical distribution of a series of solutions, and obtain the final result by calculating the expectation and variance of this distribution. Document 6 (Spacecraft two-line element set conversion method based on UKF, Spacecraft Engineering, No. 2, 2021) takes the spacecraft two-line element set as the state variable and the flight position and velocity as the measurement variables, and proposes a spacecraft two-line element set conversion method based on the Unscented Kalman Filter (UKF).

[0006] The existing methods for obtaining two-line element set do not consider the satellite maneuvering section and the accuracy limitations of the two-line element set description. To solve this problem, the present invention provides a method for obtaining two-line element set by a numerical orbit propagator. Summary of the Invention

[0007] Aiming at the defects in the prior art, the purpose of the present invention is to provide a method and system for obtaining two-line element set by a numerical orbit propagator.

[0008] According to a method for obtaining two-line element set by a numerical orbit propagator provided by the present invention, it includes:

[0009] Step S1, input the two-line element set template;

[0010] Step S2, determining the orbital time reference points corresponding to the two rows of elements;

[0011] Step S3, according to the set numerical recursor model, the numerical recursor obtains the position value sequence of the track in the set space coordinate system within a period of time;

[0012] Step S4, converting the position value sequence obtained by the numerical recursor into the true equatorial mean equinox coordinate system;

[0013] Step S5, using the two-row root number template, using the single-point method to obtain two-row root numbers;

[0014] Step S6, using the two rows of elements obtained by the single point method as the initial orbit, and using the least squares method to solve the orbit parameter estimation values in the two rows of elements under the optimal criterion;

[0015] Step S7, calculating the orbit residual sequence corresponding to the orbit parameter estimation values in the two rows of elements;

[0016] Step S8, fitting the residual using a polynomial model to obtain a compensation term.

[0017] Further, in step S3, the orbital time reference point t corresponding to the two rows of elements epoch Within the time range corresponding to the position sequence obtained by the numerical recursive device, if the time series in S3 are t0, t1, t2, ..., t N-1 , then t0≤t epoch ≤t N-1 , N is the number of time series points.

[0018] Further, step S5 includes:

[0019] Step S5.1, according to the set numerical recursor model, the numerical recursor obtains the orbit time reference point t corresponding to the two row roots. epoch The orbital position at time p num , speed v num ;

[0020] Step S5.2, set the track position p num , speed v num Transformed to the position p of the TEMEE coordinate system state and speed v state ;

[0021] Step S5.3, from position p state and speed v state Calculate the Kepler orbit elements and equatorial orbit elements. The Kepler orbit elements include the semi-major axis a, eccentricity e, and orbit inclination i nc, the right ascension of the ascending node Ω, the argument of perigee ω, the mean anomaly M, and the equatorial orbital elements include the semi-major axis a, the first element e of the eccentricity vector x , the second element e of the eccentricity vector y , the first element h of the inclination vector x , the second element h of the inclination vector y , the true longitude angular separation l v ;

[0022] Step S5.4, using the equatorial orbital elements corresponding to the position p state and the velocity v state as the starting values of the equatorial orbital elements to be updated, and the equatorial orbital elements to be updated include the semi-major axis a n , the first element e of the eccentricity vector xn , the second element e of the eccentricity vector yn , the first element h of the inclination vector xn , the second element h of the inclination vector yn , the true longitude angular separation l vn ;

[0023]

[0024] Step S5.5, generating the current two-line elements from the two-line element template and the Keplerian orbital elements;

[0025] Step S5.6, calculating the corresponding equatorial orbital elements based on the SGP4 model from the current two-line elements, and the corresponding equatorial orbital elements include the semi-major axis a c , the first element e of the eccentricity vector xc , the second element e of the eccentricity vector yc , the first element h of the inclination vector xc , the second element h of the inclination vector yc , the true longitude angular separation l vc ;

[0026] Step S5.7, calculating the equatorial orbital element residuals, including the semi-major axis residual a d , the first element residual e of the eccentricity vector xd , the second element residual e of the eccentricity vector yd , the first element residual h of the inclination vector xd , the second element residual h of the inclination vector yd , the true longitude angular separation residual l vd ;

[0027]

[0028] Step S5.8, if the equatorial orbital element residuals do not reach the threshold, then update the equatorial orbital elements:

[0029]

[0030] Step S5.8: Convert the updated equatorial orbital elements into Keplerian elements and update them to the current two-line elements.

[0031] Among them, repeat steps S5.6 to S5.8 until the equatorial orbital element residuals reach the threshold, and output the current two-line elements.

[0032] Furthermore, in step S5.5, the method for generating the current two-line elements from the two-line element template and the Keplerian orbital elements is as follows:

[0033] If the Keplerian elements are the semi-major axis a, eccentricity e, orbital inclination i nc , the right ascension of the ascending node Ω, argument of perigee ω, mean anomaly M, convert the semi-major axis a into the number of orbits around the Earth per day q, and then fill the corresponding values into the six mean element positions of the two-line elements.

[0034]

[0035] Among them, μ is the geocentric gravitational constant, with the unit of m 3 s -2 , and the semi-major axis a has the unit of m.

[0036] Furthermore, in step S6, the optimal criterion of the least squares method is that the finally output two-line elements are also recursively derived to the same moment in step S4, and the root mean square difference between the obtained orbital position sequence and the numerical integrator position sequence in step S4 is the smallest.

[0037] If the corresponding time series t in step S4 i , the position sequences are p i , and the position sequence obtained by recursively deriving the two-line elements is k i , i = 0, 1, 2,..., N - 1, the root mean square difference σ is defined as:

[0038]

[0039] Among them, |||| represents calculating the 2-norm of the vector.

[0040] Furthermore, in step S6, the iterative solution method for the least squares problem is solved using the LM algorithm.

[0041] 7. The method for obtaining two-line elements by the numerical orbital integrator according to claim 1, wherein in step S7, the calculation method of the orbital residual sequence Δ i is as follows:

[0042] Δ i = p i - k i , i = 0, 1, 2,..., N - 1.

[0043] Further, in step S8, using a polynomial model to fit the residuals includes:

[0044] Step S8.1, converting the time series t in S3 i into the relative time corresponding to the orbital time reference point t of the relative two-line elements, i = 0, 1, 2, ……, N-1; epoch of the relative two-line elements, i = 0, 1, 2, ……, N-1;

[0045] Step S8.2, respectively using a polynomial model to fit the relative time in the three-dimensional spatial directions of the orbital residual sequence Δ i to obtain the polynomial coefficients as compensation terms.

[0046] Further, in step S1, the two-line element template includes: satellite catalog number, international designator, initial value of the first derivative of the mean motion, initial value of the second derivative of the mean motion, initial value of the BSTAR coefficient of the drag term, number of orbits flown since launch.

[0047] According to a two-line element system obtained by a numerical orbit propagator provided by the present invention, it includes:

[0048] Module M1, inputting the two-line element template;

[0049] Module M2, determining the orbital time reference point corresponding to the two-line elements;

[0050] Module M3, according to the set numerical propagator model, obtaining the position value sequence of the orbit in the set space coordinate system within a period of time by the numerical propagator;

[0051] Module M4, converting the position value sequence obtained by the numerical propagator into the true equator mean equinox coordinate system;

[0052] Module M5, obtaining the two-line elements by using the single-point method from the two-line element template;

[0053] Module M6, using the two-line elements obtained by the single-point method as the initial orbit, and using the least squares method to solve the estimated value of the orbit parameters in the two-line elements under the optimal criterion;

[0054] Module M7, calculating the orbital residual sequence corresponding to the estimated value of the orbit parameters in the two-line elements;

[0055] Module M8, using a polynomial model to fit the residuals to obtain the compensation terms.

[0056] Compared with the prior art, the present invention has the following beneficial effects:

[0057] In view of the problem that the conversion accuracy of the numerical orbit propagator is limited during the process of converting to two-line elements, the present invention proposes a new method for obtaining two-line elements. First, a set of two-line elements is obtained as an initial value through the single-point method, then the least squares method is used to solve for the two-line elements with the minimum root mean square error over a period of time, and finally a compensation term is generated by polynomial fitting of the residuals. The method of the present invention is reasonable, simple in calculation, and easy to implement, and can be effectively applied to obtaining two-line elements from a numerical orbit propagator. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Other features, objects, and advantages of the present invention will become more apparent by reading the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0059] Figure 1 It is a flowchart of the present invention.

[0060] Figure 2 It is the residual between the two-line elements and the numerical propagator without a compensation term.

[0061] Figure 3 It is the residual between the two-line elements and the numerical propagator with a compensation term. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0062] The present invention will be described in detail below with reference to specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that those of ordinary skill in the art can make several changes and improvements without departing from the concept of the present invention. These all fall within the protection scope of the present invention.

[0063] As Figure 1 shown, a method for obtaining two-line elements from a numerical orbit propagator provided by the present invention includes:

[0064] Step S1, input a two-line element template.

[0065] Step S2, determine the orbital time reference point corresponding to the two-line elements.

[0066] Step S3, according to the set numerical propagator model, obtain a sequence of position values of the orbit in a set space coordinate system over a period of time from the numerical propagator.

[0067] Step S4, transform the sequence of position values obtained by the numerical propagator into the true equator mean equinox coordinate system.

[0068] Step S5, obtain two-line elements using the single-point method from the two-line element template.

[0069] Step S6: Use the number of two lines obtained by the single-point method as the initial orbit, and use the least squares method to solve the estimated values of the orbital parameters in the two-line elements under the optimal criterion.

[0070] Step S7: Calculate the orbital residual sequence corresponding to the estimated values of the orbital parameters in the two-line elements.

[0071] Step S8: Use the polynomial model to fit the residuals to obtain the compensation term.

[0072] The data format of the two-line elements TLE is fixed, also known as 3LE (with an additional satellite name in the 0th line). The first line includes the satellite catalog number, orbital reference time, international designator, initial value of the first derivative of the mean motion, initial value of the second derivative of the mean motion, initial value of the BSTAR coefficient of the drag term, and checksum. The second line includes the satellite catalog number, orbital inclination, eccentricity, argument of perigee, right ascension of the ascending node, number of revolutions around the Earth per day, number of revolutions since launch, and checksum. When using the SGP4 model to perform orbital propagation based on TLE, it involves the orbital element information in the second line of TLE, while information such as the satellite target number, international designator, and number of revolutions since launch is managed as part of the target catalog (Catalogue). When the numerical orbit propagator obtains the two-line elements, the goal is to obtain the orbital element information in the second line. Therefore, other information in the two-line element format needs to be provided in the two-line element template in advance.

[0073] The two-line elements are defined in the true equator mean equinox coordinate system, while the numerical orbit propagator may not be in the true equator mean equinox coordinate system according to specific settings. For example, the J2000 geocentric inertial coordinate system is a commonly used coordinate system. When obtaining the two-line elements based on the numerical orbit propagator, it is necessary to ensure that when performing propagation based on the SGP4 model with the finally obtained two-line elements, the difference between the propagation results and the propagation results of the numerical orbit propagator within a certain period of time should be within the allowable range of the conversion accuracy.

[0074] To compare the differences in the conversion results, it is first necessary to unify the time and space of the propagation results. When the numerical propagator does not define the true equator mean equinox coordinate system, it is necessary to first convert the position vector obtained by its propagation to the mean equinox coordinate system.

[0075] If the numerical propagator obtains the positions in the true equator mean equinox coordinate system at a certain period of time (the time series is t0, t1, t2,..., t N-1 ) as p i (i = 0, 1, 2,..., N - 1), where N is the number of time series points. The position sequence obtained by propagating the two-line elements is k i (i = 0, 1, 2,..., N - 1). Use the root mean square difference between the two as the evaluation criterion. The root mean square difference σ is defined as

[0076]

[0077] Among them, |||| represents the 2-norm of the computational vector.

[0078] To achieve this least-squares optimal criterion, it is necessary to estimate the two-line elements. If a reliable initial value can be obtained and the least-squares problem is transformed into a "small residual" problem, the LM algorithm can be used to perform iterative solution based on the initial value.

[0079] To obtain a reliable initial value, the single-point method can be used. The difference between the single-point method and the least-squares method is that the optimal criterion of the least-squares method is established on the position sequence over a period of time, while the single-point method only considers the orbital state at a certain moment. The single-point method only needs to input the orbital state at a certain moment, and its result only guarantees that the orbits are close at that moment. If the state at that moment is used for recursion, there may be a large orbital difference.

[0080] The single-point method is also based on iterative calculation. To avoid the singularity problem existing in the Keplerian orbital elements, the non-singular equatorial orbital elements are used for orbital estimation calculation.

[0081] The single-point method iteration also requires an initial value. The orbital time reference point t corresponding to the two-line elements epoch The orbital position p corresponding to the time numerical integrator num , velocity v num can be transformed into the TEME coordinate system to obtain the position p state and velocity v state .

[0082] From the orbital position p state and velocity v state , the Keplerian orbital elements (semi-major axis a, eccentricity e, orbital inclination i nc , right ascension of the ascending node Ω, argument of perigee ω, mean anomaly M) and the equatorial orbital elements (semi-major axis a, first element e of the eccentricity vector x , second element e of the eccentricity vector y , first element h of the inclination vector x , second element h of the inclination vector y , true longitude angular distance l v ) can be calculated.

[0083] The equatorial orbital elements at this time can be used as the initial value of the iteration. The parameters to be updated include: semi-major axis a n , first element e of the eccentricity vector xn , second element e of the eccentricity vector yn , first element h of the inclination vector xn , second element h of the inclination vector yn , true longitude angular distance lvn Its initial value is set as:

[0084]

[0085] Based on the initial value and the parameters in the two-line element set template, the current two-line element set can be generated. At this time, the orbital elements in the two-line element set are only a rough estimate, and after simple conversion from Keplerian elements, they are filled into the corresponding positions, involving the conversion of the semi-major axis a to the number of orbits q around the Earth per day.

[0086]

[0087] where, μ is the geocentric gravitational constant, with the unit of m 3 s -2 , and the semi-major axis a has the unit of m.

[0088] Based on the current two-line element set, the corresponding actual equatorial orbital elements (semi-major axis a c , the first element e of the eccentricity vector xc , the second element e of the eccentricity vector yc , the first element h of the inclination vector xc , the second element h of the inclination vector yc , the true longitude angular distance l vc ) can be calculated based on the SGP4 model.

[0089] Calculate the equatorial orbital element residuals (semi-major axis residual a d , the first element residual e of the eccentricity vector xd , the second element residual e of the eccentricity vector yd , the first element residual h of the inclination vector xd , the second element residual h of the inclination vector yd , the true longitude angular distance residual l vd ).

[0090]

[0091] If the residuals meet certain threshold conditions, it means that the orbital states at single points are the same. If the equatorial orbital element residuals do not reach the threshold, then update the equatorial orbital elements.

[0092]

[0093] Then convert the updated equatorial orbital elements to Keplerian elements and update them to the current two-line element set. Repeat the iteration until the equatorial orbital element residuals reach the threshold, and the current two-line element set can be output as the initial value for the least squares solution.

[0094] Within a not very long period of time (usually not exceeding one orbital period) where the time span covers the reference time of the two rows of roots, the difference between the two rows of roots obtained by the single-point method and the orbital parameters to be estimated is not large, satisfying the "small remainder" condition. The iterative solution method of this least squares problem adopts the LM (Levenberg-Marquardt) algorithm to solve it.

[0095] After solving the final two rows of elements, although the estimated value meets the optimal criterion within a given time period, there may still be a certain gap with the numerical orbit recurser, especially when the orbit is maneuvering. At this time, we further obtain the orbit residual sequence corresponding to the orbit parameter estimates in the two rows of elements, fit the residual sequence, and use the fitted coefficients as compensation terms. This can ultimately ensure that the difference between the numerical orbit recurser and the two rows of elements after adding the compensation term is always within the range allowed by the conversion accuracy.

[0096] Orbital residual sequence Δ i The calculation method for (i=0, 1, 2, ..., N-1) is

[0097] Δ i =p i -k i (Formula 6)

[0098] In order to simplify the number of compensation parameters, the absolute time series t i (i = 0, 1, 2, ..., N-1) is converted into the orbital time reference point t corresponding to the two row elements epoch After the relative time series, the orbital residual sequence Δ i In the three-dimensional directions of space, polynomial models are used for fitting.

[0099] The following test data is used to verify the effectiveness of the method of the present invention. For a numerical orbit recursion device defined in the J2000 geocentric inertial coordinate system, the perturbation force definition includes: Earth's gravity field, atmospheric drag, solar light pressure, third body gravity, relativistic effect, and maneuvering force. Maneuvering force is defined as the X direction of the satellite orbit, 1mm / s 2 , action time 2000s. Using the method of the present invention, firstly, the number of roots of the two rows obtained by the single point method is

[0100]

[0101] Using this result as the initial value, the final two rows of roots are

[0102]

[0103] At this time, the residuals of the two rows of roots and the numerical recursion are as shown in the attached Figure 2As shown in the figure, after fitting with a cubic polynomial, the residuals of the two rows of roots and the numerical orbit recursion considering the compensation term are shown in the attached figure. Figure 3 It is shown that the difference between the two rows of roots obtained by the method of the present invention and the numerical recursion is at the meter level, which meets the accuracy requirement of the conversion.

[0104] The present invention also provides a numerical orbit recurser for obtaining a two-row element system, which can be implemented by executing the process steps of the method for obtaining two-row elements using the numerical orbit recurser, that is, those skilled in the art can understand the method for obtaining two-row elements using the numerical orbit recurser as a preferred implementation of the method for obtaining two-row elements using the numerical orbit recurser. The system includes:

[0105] Module M1, input two lines of root number template.

[0106] Module M2 determines the orbital time reference points corresponding to the two rows of elements.

[0107] Module M3, according to the set numerical recursor model, obtains the position value sequence of the track in the set space coordinate system within a period of time by the numerical recursor.

[0108] Module M4 converts the position value sequence obtained by the numerical recursive device into the true equatorial vernal equinox coordinate system.

[0109] Module M5, using two rows of root number templates, uses the single-point method to obtain two rows of root numbers.

[0110] Module M6 uses the two rows of roots obtained by the single-point method as the initial orbit and uses the least squares method to solve the orbital parameter estimates in the two rows of roots under the optimal criterion.

[0111] Module M7 calculates the orbital residual sequence corresponding to the orbital parameter estimates in the two rows of elements.

[0112] Module M8 uses a polynomial model to fit the residual and obtain a compensation term.

[0113] Those skilled in the art know that, in addition to realizing the system and its various devices, modules, and units provided by the present invention in a purely computer-readable program code, it is entirely possible to realize the same functions in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps. Therefore, the system and its various devices, modules, and units provided by the present invention can be considered as a hardware component, and the devices, modules, and units included therein for realizing various functions can also be regarded as structures within the hardware component; the devices, modules, and units for realizing various functions can also be regarded as both software modules for realizing the method and structures within the hardware component.

[0114] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. Without conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. A method for obtaining two-line elements by a numerical orbit propagator, characterized in that, Including: Step S1, input the two-line element set template; Step S2, determine the orbital time reference points corresponding to the two-line elements; Step S3, according to the set numerical integrator model, obtain the position value sequence of the orbit in the set space coordinate system within a period of time by the numerical integrator; Step S4, transform the position value sequence obtained by the numerical integrator into the true equator mean equinox coordinate system; Step S5, obtain the two-line elements by using the single-point method from the two-line element set template; Step S6, take the two-line elements obtained by the single-point method as the initial orbit, and use the least squares method to solve the estimated values of the orbital parameters in the two-line elements under the optimal criterion; Step S7, calculate the orbital residual sequence corresponding to the estimated values of the orbital parameters in the two-line elements; Step S8, use the polynomial model to fit the residuals to obtain the compensation term.

2. The method for obtaining two-line elements of a numerical orbit propagator according to claim 1, characterized in that, In step S3, the orbital time reference point t corresponding to the number of roots in two rows epoch Within the time range corresponding to the position sequence obtained by the numerical recursive device, if the time series in S3 are t0, t1, t2,..., t N-1 , then t0 ≤ t epoch ≤ t N-1 , where N is the number of time series points.

3. The method for obtaining two-line elements by the numerical orbit propagator according to claim 1, characterized in that, Step S5 includes: Step S5.1: According to the set numerical recurrence model, obtain the orbital time reference point \(t\) corresponding to the orbital two-line elements by the numerical recurrence, epoch the orbital position \(p\) at the num moment, and the velocity \(v\) num ; Step S5.2, convert the orbital position p num , velocity v num to the position p state and velocity v state in the TEME coordinate system; Step S5.3, from position p state and velocity v state calculate the Keplerian orbital elements and equatorial orbital elements. The Keplerian orbital elements include the semi-major axis a, eccentricity e, orbital inclination i nc , right ascension of the ascending node Ω, argument of perigee ω, mean anomaly M. The equatorial orbital elements include the semi-major axis a, the first element e of the eccentricity vector x , the second element e of the eccentricity vector y , the first element h of the inclination vector x , the second element h of the inclination vector y , true longitude angular distance l v ; Step S5.4, from position p state and velocity v state the corresponding equatorial orbit elements are used as the starting values of the equatorial orbit elements to be updated. The equatorial orbit elements to be updated include semi-major axis a n , the first element e of the eccentricity vector xn , the second element e of the eccentricity vector yn , the first element h of the inclination vector xn , the second element h of the inclination vector yn , the true longitude angular distance l vn ; Step S5.5, generate the current two-line elements from the two-line element set template and the Keplerian orbital elements; Step S5.6, calculate the corresponding equatorial orbital elements based on the SGP4 model from the current two-line elements, and the corresponding equatorial orbital elements include the semi-major axis a c , the first element e of the eccentricity vector xc , the second element e of the eccentricity vector yc , the first element h of the inclination vector xc , the second element h of the inclination vector yc , the true longitude angular distance l vc ; Step S5.7, calculate the equatorial orbital element residuals, including the semi-major axis residual a d , the first element residual e of the eccentricity vector xd , the second element residual e of the eccentricity vector yd , the first element residual h of the inclination vector xd , the second element residual h of the inclination vector yd , the true longitude angular distance residual l vd ; Step S5.8, if the equatorial orbital element residuals do not reach the threshold, update the equatorial orbital elements: Step S5.8, transform the updated equatorial orbital elements into Keplerian elements and update them to the current two-line elements; Among them, repeat steps S5.6 - 5.8 until the equatorial orbital element residuals reach the threshold, and output the current two-line elements.

4. The method for obtaining two-row elements by a numerical orbit recursion device according to claim 3, characterized in that: In step S5.5, the method for generating the current two-line elements from the two-line element set template and the Keplerian orbital elements is: If the Keplerian elements are the semi-major axis a, eccentricity e, inclination i nc , the right ascension of the ascending node Ω, argument of perigee ω, mean anomaly M, the semi-major axis a will be converted to the number of orbits q around the Earth per day, and then the corresponding values will be filled into the six mean element positions of the two rows of elements; where μ is the geocentric gravitational constant, with the unit of m 3 s -2 , and the semi-major axis a has the unit of m.

5. The method for obtaining two-line elements by the numerical orbit propagator according to claim 1, characterized in that In step S6, the optimal criterion of the least squares method is that the two-line elements finally output are also propagated to the same moment in step S4, and the root mean square difference between the obtained orbital position sequence and the numerical integrator position sequence in step S4 is the smallest; If the corresponding position sequences for the corresponding time series t in step S4 i are successively p i and the position sequence obtained by recursive derivation of the number of rows in two rows is k i where i = 0, 1, 2, ……, N - 1, the root mean square deviation σ is defined as: Among them, || || represents calculating the 2-norm of the vector.

6. The method for obtaining two-line elements of a numerical orbit propagator according to claim 1, characterized in that, In step S6, the iterative solution method for the least squares problem uses the Levenberg-Marquardt (LM) algorithm to solve.

7. The method for obtaining two-line elements by the numerical orbit propagator according to claim 1, characterized in that In step S7, the calculation method of the orbit residual sequence Δ i is as follows: Δ i = p i - k i , i = 0, 1, 2, ……, N - 1.

8. The method for obtaining two-line elements of a numerical orbit propagator according to claim 1, characterized in that, In step S8, using the polynomial model to fit the residuals includes: Step S8.1, convert the time series t in S3 i into the relative time corresponding to the orbital time reference point t of the number of roots in two adjacent rows epoch , where i = 0, 1, 2, ……, N - 1; Step S8.

2. In the three-dimensional spatial directions of the orbital residual sequence Δ i , use a polynomial model to fit the relative time respectively to obtain the polynomial coefficients as compensation terms.

9. The method for obtaining two-line elements of a numerical orbit propagator according to claim 1, characterized in that, In step S1, the two-line element set template includes: satellite catalog number, international designator, initial value of the first derivative of the mean motion, initial value of the second derivative of the mean motion, initial value of the BSTAR coefficient of the drag term, number of orbits flown since launch.

10. A two-line element set acquisition system for a numerical orbit propagator, characterized in that, Including: Module M1, input the two-line element set template; Module M2, determine the orbital time reference points corresponding to the two-line elements; Module M3, according to the set numerical integrator model, obtain the position value sequence of the orbit in the set space coordinate system within a period of time by the numerical integrator; Module M4, transform the position value sequence obtained by the numerical integrator into the true equator mean equinox coordinate system; Module M5, obtain the two-line elements by using the single-point method from the two-line element set template; Module M6, take the two-line elements obtained by the single-point method as the initial orbit, and use the least squares method to solve the estimated values of the orbital parameters in the two-line elements under the optimal criterion; Module M7, calculate the orbital residual sequence corresponding to the estimated values of the orbital parameters in the two-line elements; Module M8, use the polynomial model to fit the residuals to obtain the compensation term.

Images