A semi-analytical and semi-numerical orbit calculation method

Abstract

The present invention relates to a semi-analytical and semi-numerical orbit calculation method. The method is mainly used to calculate the position vector at any time from the number of oscillating roots at a given time. First, the idea of ​​the analytical method is used to iteratively calculate the average number of roots at a given time, and then two critical points are constructed near the initial time, and the average number of roots and the rate of change of the average number of roots at the critical points are calculated; combined with the Chebyshev polynomials, the numerical integration of the average number of roots is completed, and the average number of roots at any given time is obtained, and finally the number of oscillating roots at the time is calculated, and then the target position vector is obtained. The present invention completes the orbital characteristic analysis and precise calculation method of the resonance region, improves the adaptability of the orbital target determination in the resonance region; takes into account the advantages of the analytical method and the numerical method, and has a better balance between calculation accuracy and speed in orbit cataloging.

Description

Technical Field

[0001] The invention belongs to the field of aerospace measurement and control data processing, relates to spacecraft orbit calculation, and can be used for initial orbit calculation based on remote external measurement data of aerospace measurement ships and for prediction using orbits. Background Art

[0002] When launching a spacecraft, the measurement ship is responsible for the spacecraft's orbital entry mission. One of its important tasks is to calculate the initial orbital roots after the spacecraft separates from the rocket. Its main function is to quickly determine whether the rocket launch is successful or not, and to grasp the initial operating status of the spacecraft after entering the space orbit, laying the foundation for subsequent measurement and control stations to conduct measurement and control of the spacecraft.

[0003] When survey vessels are operating at sea, orbit calculations are based on data sources such as external measurements from the vessel, rocket telemetry trajectories, and GNSS trajectories. Orbit calculation methods are generally categorized into analytical and numerical methods. The analytical method treats the target's motion around the Earth as a two-body problem, using algebraic methods to calculate the target's orbital elements and approximate the perturbations in the form of harmonics. The numerical method calculates the forces acting on the target as it orbits the Earth, combines these with initial values, and then calculates the target's orbit through numerical integration.

[0004] The above two orbit calculation methods have their own advantages and disadvantages when implemented. The computational complexity is small when using the analytical method for orbit calculation, but the accuracy is limited in theory because the various influencing factors of perturbation are not fully considered. When using the numerical method for calculation, different orders of gravity models and resistance models can be used according to different accuracy requirements, which can effectively improve the calculation accuracy. However, its disadvantage is that under general accuracy requirements, the computational complexity is large, which affects the stability and timeliness of the calculation. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a semi-analytical and semi-numerical orbit calculation method for the above-mentioned existing technologies, which combines the respective advantages of analytical and numerical methods and better achieves a balance between calculation accuracy and computational complexity.

[0006] The technical solution adopted by the present invention to solve the above problems is: a semi-analytical and semi-numerical orbit calculation method, which includes the following steps:

[0007] Step A: Calculate the square root number at the starting time.

[0008] Since the orbital elements used during input are generally instantaneous elements (kissing elements), the spacecraft motion they describe includes not only gravity but also various perturbations. When performing orbital calculations, to simplify the analytical calculations, the long-period and short-period perturbation terms are typically separated, and the orbital elements from the two-body problem are used for analytical calculations. This means that subsequent calculations are performed using flat elements.

[0009] Step B: Determine the critical point of the interval and calculate the approximate square root number at the critical point.

[0010] To reduce the computational effort during numerical calculations and improve the speed of the numerical integration process, this method uses the analytical results as the initial conditions for the numerical integration. Therefore, during each numerical integration, the initial value of the integral at the starting point must be independently calculated. The starting point of the integration is the critical point of each interval, and calculating the initial value of the integral requires calculating the approximate square root at the critical point.

[0011] Step C: Calculate the rate of change of the mean root number at the critical point.

[0012] The square root of the number can basically reflect the instantaneous position of the aircraft's motion, but because it discards the long and short period terms, its speed at the instantaneous moment needs to take into account the comprehensive forces caused by the earth's gravity and periodic perturbations. When performing numerical integration, it is necessary to calculate a more accurate aircraft motion speed, and it is necessary to consider the change of the square root of the number, that is, the square root of the number variation rate.

[0013] Step D: Calculate the Chebyshev polynomial coefficients.

[0014] Since the interpolation polynomials selected as the interpolation nodes when performing Lagrange interpolation are optimal in terms of the upper bound of the error, this solution uses Chebyshev polynomials as numerical integration nodes.

[0015] Step E: Calculate the number of square roots at the end point by numerical integration.

[0016] By performing numerical integration in each integral interval, the square root number at the end point of the integral is calculated, and then the final square root number at the end point is obtained through Chebyshev polynomial fitting, which is used to calculate the aircraft position and velocity vector at the end point.

[0017] Step F: Calculate the field harmonic terms.

[0018] When calculating the actual position and velocity of the aircraft by square root number, the influence of the field harmonic term caused by the perturbation also needs to be considered.

[0019] Step G: Calculate the position and velocity vector of the aircraft at the end point.

[0020] In the actual calculation process, users generally use the position and velocity vector of the aircraft as the final product. Therefore, the target square root, the first-stage periodic term, and the field harmonic term are combined to calculate the target position and velocity vector as the final result.

[0021] Compared with the prior art, the advantages of the present invention are:

[0022] The present invention is primarily used to calculate the position vector at any given moment from the number of oscillating roots at a given moment. First, the average number of roots at a given moment is iteratively calculated using analytical methods. Then, two critical points are constructed near the initial moment, and the average number of roots and the rate of change of the average number of roots at the critical points are calculated. Combined with Chebyshev polynomials, the numerical integration of the average number of roots is performed to obtain the average number of roots at any given moment. Finally, the number of oscillating roots at that moment is calculated, thereby obtaining the target position vector.

[0023] The present invention completes the analysis and precise calculation method of the resonance region orbit characteristics, improves the adaptability of the resonance region orbit area target orbit determination; takes into account the advantages of the analytical method and the numerical method, and has a better balance between the calculation accuracy and speed in the orbit cataloging. DETAILED DESCRIPTION

[0024] The present invention is described in further detail below with reference to the examples.

[0025] Step A: Calculate the square root number at time t0

[0026] This calculation process involves the use of square roots many times. To prevent singularities in the calculation of the perturbation motion equation, the square roots used in this invention use first-class singularity-free variables to represent orbital parameters, that is, the parameters ξ, η, and λ are introduced:

[0027] ξ=esinω

[0028] η=ecosω

[0029] λ=M+ω

[0030] The first type of singularity-free orbital parameter representation method is equivalent to the traditional Kepler six-element representation method, but can ensure that no singularity will appear in the perturbation motion equation under arbitrary eccentricity (0≤e<1).

[0031] The input parameters of the algorithm are the osculating root number σ(a, i, e, Ω, ω, M) at time t0, let a0 = a, i0 = i, Ω0 = Ω, ξ0 = esinω, η0 = ecosω, λ0 = M + ω, where

[0032] a is the semi-major axis of the orbit;

[0033] i is the orbital inclination;

[0034] e is the orbital eccentricity, which refers to the dihedral angle between the elliptical orbit and the reference plane;

[0035] Ω is the right ascension of the ascending node, which is the angle from the x-axis on the reference plane to the center of the Earth pointing to the ascending node;

[0036] ω is the argument of perigee, which is the angle from the center of the Earth pointing to the ascending node to the center of the Earth pointing to the perigee on the reference plane;

[0037] M is the mean anomaly, and each element is substituted into the following square root equation:

[0038]

[0039]

[0040]

[0041] r, u and f are intermediate parameters, where f is the true anomaly, and its relationship with other orbital parameters is:

[0042]

[0043] u=f+M

[0044]

[0045] M=E-esinE;

[0046]

[0047]

[0048]

[0049] make Repeat the cycle in sequence.

[0050] If |a i -a i-1 |<10 -8 , then the number of roots at the end of the cycle can be taken as the number of average roots at time t0

[0051] Step B: Determine the critical point t i , and calculate the approximate number of square roots at the critical point:

[0052] Select 2 points in [t0-h,t0+h]:

[0053]

[0054] Here, h is the step size for calculating the perturbations, which is generally 1-2 days for cataloging and orbit determination.

[0055] Calculate t i Approximate square root of

[0056] For Ω, ω,

[0057] in, is the first-order long-term variability of the root, that is:

[0058]

[0059] Among them, J2 refers to the first-order coefficient of the earth belt harmonic term, which is a constant, J2≈1.08263e-3, n is the mean motion speed, p=a(1-e 2 ), all calculated using square roots.

[0060] For a, i, e, we have:

[0061] at last

[0062] Using t i The square root of Calculate the coefficient of the m-daily term spare.

[0063] Step C: Calculate t i The average root rate of change at :

[0064] Constant term, long period term,

[0065] in, To the first order osculating root, is a first-order short-period term. Its expression is as follows:

[0066]

[0067]

[0068]

[0069]

[0070]

[0071]

[0072] In the formula, the roots are all flat roots. Where n is the flat motion speed, p=a(1-e 2), f is the true anomaly, which is calculated using square roots, and p means the semi-diameter of the orbit.

[0073] In the expression of the short-period term, as well as in the subsequent perturbation motion equations and calculation formulas such as (S, T, W), there is e k sin(nu±kf), It must be expanded and calculated, such as:

[0074] e sin f=ηsinu-ξcosu,

[0075] And so on.

[0076] f in the formula (0) (σ (1) ),f (1) (σ (1) ),f (2) (σ (1) ) are all vectors:

[0077]

[0078] is the square root rate of change caused by J2; is the square root rate of change caused by other perturbations, f (1) ,f (2) The expression form is the same, which is the result of substituting the perturbation acceleration into the Gaussian equation. The expression is as follows:

[0079]

[0080]

[0081]

[0082]

[0083]

[0084]

[0085] Calculate f (1) When , (S, T, W) is substituted with the perturbation acceleration caused by J2, that is:

[0086]

[0087]

[0088]

[0089] Calculate f (2)When , (S, T, W) is substituted with the sum of accelerations caused by other perturbations, that is:

[0090] S2=∑S 2i ;

[0091] T2=∑T 2i ;

[0092] W2=∑W 2i

[0093] (S,T,W), (Ω,Ω′,W), (Ω,Ω1,k), and (i,j,k) are all vector systems used in satellite orbit calculations. The origin of the reference system is the center of the earth. In terms of physical meaning, they can be simply described as follows:

[0094] k points to the North Pole, i points to the direction of the vernal equinox, and (i, j, k) forms a right-handed system;

[0095] With k as the rotation axis, i and j rotate around k by an angle Ω (the right ascension of the ascending node) so that Ω points to the direction of the ascending node of the orbit, and (Ω, Ω1, k) forms a right-handed system;

[0096] With Ω as the rotation axis, Ω1,k rotates around Ω by an angle i so that W points to the direction of the orbital normal, and (Ω, Ω′, W) forms a right-handed system;

[0097] With W as the rotation axis, Ω, Ω′ rotate around W by an angle ω+f (true latitude angle) so that S points to the radial direction of the satellite, and (S, T, W) forms a right-handed system;

[0098] Perturbations are considered here, including high-order harmonic perturbations, atmospheric drag perturbations, solar pressure perturbations, solar and lunar perturbations, and tidal perturbations.

[0099] Average root rate of change It is obtained by numerical averaging, the specific method is as follows: Select [NF is optional. For near-circular orbits, it is generally taken as 18. It changes according to the size of the eccentricity. NF=18(e≤0.1), NF=18+int[(e-0.1)*20](e>0.1)], according to the given The first-order short-period term can be calculated, which can be calculated So the algorithm for the square root rate of change is:

[0100]

[0101] Step D: Integrate to find the Chebyshev polynomial coefficient a k (i)

[0102] Use Chebyshev iteration to integrate the perturbed equations of motion with square roots:

[0103]

[0104] Step E: Find the square root at any time t

[0105]

[0106] Among them, T k (τ) is the k-th order Chebyshev polynomial,

[0107] Step F: Calculate the field harmonic perturbation

[0108] There are three types of harmonic perturbation terms that need to be considered in this method:

[0109] (a) m-daily terms, i.e., terms with j = 0

[0110] (b) Possible primary and secondary resonance terms, i.e., terms with j = 1, 2

[0111] (c) For the larger short-period term, only the largest field harmonic coefficient J needs to be considered. 22 =2.81×10 -6 The specific calculation process is as follows:

[0112]

[0113]

[0114]

[0115]

[0116]

[0117]

[0118] Step G: Find the osculating elements and satellite coordinates at any time t

[0119]

[0120] in, is the first-order short-period term, is the harmonic perturbation term,

[0121] Then the coordinates of the satellite in the orbital coordinate system

[0122]

[0123] In addition to the above embodiments, the present invention also includes other implementation methods. Any technical solutions formed by equivalent transformation or equivalent replacement should fall within the scope of protection of the claims of the present invention.

Claims

1. A semi-analytical and semi-numerical orbit calculation method, characterized by: The method comprises the following steps: Step A: Calculate the number of square roots at the starting time When performing orbit calculations, the long-period and short-period perturbation terms are separated, and the orbital roots under the two-body problem are used for analysis, that is, the flat roots are used for subsequent calculations. Step B: Determine the critical point of the interval and calculate the approximate square root number at the critical point The results of the analytical method are used as the initial conditions for numerical integration. In each numerical integration process, the initial value of the integral at the starting point of the integral is independently calculated. The starting point of the integral is the critical point of each interval. To calculate the initial value of the integral, it is necessary to calculate the approximate square root number at the critical point. Step C: Calculate the rate of change of the average root number at the critical point When performing numerical integration, if you need to calculate the aircraft's motion speed more accurately, you need to consider the change of the square root number, that is, the square root number variation rate; Step D: Calculate Chebyshev polynomial coefficients Since the interpolation polynomials selected as the interpolation nodes of the Chebyshev polynomial roots are optimal in terms of the upper bound of the error, the Chebyshev polynomials are used as the numerical integration nodes. Step E: Calculate the number of square roots at the end point by numerical integration By performing numerical integration in each integral interval, the square root number at the end point of the integral is calculated, and then the final square root number at the end point is obtained by fitting the Chebyshev polynomial, which is used to calculate the position and velocity vector of the aircraft at the end point. Step F: Calculate the harmonic terms When calculating the actual position and velocity of the aircraft from the square root number, the influence of the field harmonic term caused by the perturbation needs to be considered; Step G: Calculate the aircraft position and velocity vector at the end point In the actual calculation process, the user uses the position and velocity vector of the aircraft as the final result. Therefore, the target square root, the first-stage periodic term, and the field harmonic term are combined to calculate the target position and velocity vector as the final result.

2. The semi-analytical and semi-numerical orbit calculation method according to claim 1, characterized in that: The input parameters in step A are the osculating element σ(a, i, e, Ω, ω, M) at time t0, and the output target is the position vector of the aircraft at time t The specific meanings of each variable are as follows: a: orbital semi-major axis; i: orbital inclination; e: orbital eccentricity, which refers to the dihedral angle between the elliptical orbit and the reference plane, 0≤e<1; Ω: right ascension of the ascending node, the angle from the x-axis on the reference plane to the center of the Earth pointing to the ascending node; ω: Argument of perigee, the angle from the center of the Earth pointing to the ascending node to the center of the Earth pointing to the perigee on the reference plane; M: mean anomaly, and substitute each element into the following square root equation: Among them, r, u and f are intermediate parameters, f is the true anomaly, J2 refers to the first-order coefficient of the earth belt harmonic term, which is a constant, J2≈1.08263e -3 , and the relationship between it and other orbital parameters is: u=f+M, M=Ee sin E.

3. The semi-analytical and semi-numerical orbit calculation method according to claim 2, characterized in that: In order to prevent singularities from occurring in the calculation of the perturbation equation of motion, the square roots used are first-order non-singularity variables to represent the orbital parameters, that is, the parameters ξ, η, and λ are introduced, where ξ=esinω η=ecosω λ=M+ω The first type of singularity-free orbital parameter representation method is equivalent to the traditional Kepler six-element representation method, which can ensure that no singularity will appear in the perturbation motion equation under any eccentricity. At the same time, the square root formula is used to calculate make a0=a,i0=i,Ω0=Ω,ξ0=esinω,η0=ecosω,λ0=M+ω,repeat in sequence. If |a i -a i-1 |<10 -8 , then the number of roots at the end of the cycle can be taken as the number of average roots at time t0 4. The semi-analytical and semi-numerical orbit calculation method according to claim 3, characterized in that: The specific implementation of step B is as follows: Select 2 points in [t0-h,t0+h]: Here, h is the step size for calculating perturbations, which is 1-2 days for cataloging orbit determination. Calculate t i Approximate square root of For Ω, ω, in, is the first-order long-term variability of the root, that is: Where n is the horizontal motion speed, p=a(1-e 2 ), are calculated using square roots, For a, i, e, we have: at last Using t i The square root of Calculate the coefficient of the m-daily term spare.

5. The semi-analytical and semi-numerical orbit calculation method according to claim 4, characterized in that: In step C, the calculation of the average root number variation at the critical point is specifically calculated using the following formula: i The average root rate of change at : Constant term, long period term, in, To the first order osculating root, is a first-order short-period term.

6. The semi-analytical and semi-numerical orbit calculation method according to claim 5, characterized in that: Average root rate of change It is obtained by numerical averaging, the specific method is as follows: Select According to the given The first-order short-period term can be calculated, which can be calculated So the algorithm for the square root rate of change is:

7. The semi-analytical and semi-numerical orbit calculation method according to claim 5, characterized in that: The perturbed equation of motion using Chebyshev iteration to integrate the square roots is as follows:

8. The semi-analytical and semi-numerical orbit calculation method according to claim 7, characterized in that: In step E, the square root number at any time t is calculated using the following formula: Among them, T k (τ) is the k-th order Chebyshev polynomial, 9. The semi-analytical and semi-numerical orbit calculation method according to claim 8, characterized in that: In step F, the harmonic perturbation is obtained by calculation:

10. The semi-analytical and semi-numerical orbit calculation method according to claim 9, characterized in that: In step G, the following formula is used to calculate the number of osculating elements and satellite coordinates at any time t: in, is the first-order short-period term, is the field harmonic perturbation term calculated in step F, and the coordinates of the satellite in the orbital coordinate system

Images