Method for accurately and efficiently calculating dense ephemeris of high-eccentricity orbit

Abstract

A method for accurately and efficiently calculating a dense ephemeris of a high-eccentricity orbit is provided. With respect to the ephemeris calculation of the high-eccentricity orbit, the method constructs uneven interpolation nodes through time transformation and interpolates by an interpolation polynomial based on uneven interpolation nodes to obtain a dense ephemeris, which significantly improves the calculation efficiency and accuracy. Based on a large-scale numerical experiment, the method derives an optimal universal value (that is, 0.3) of a transformation parameter for all orbital eccentricities and various interpolation polynomials. In the case of using the optimal universal value of the transformation parameter δ, the method further verifies the Hermite interpolation polynomial as the preferable one among various interpolation polynomials.

Claims

1. A method for accurately and efficiently calculating a dense ephemeris of a high-eccentricity orbit, wherein the method implements ephemeris calculations of a space object by solving a dynamic equation with initial values and uses a construction algorithm for uneven interpolation nodes, and the method comprises the following steps: step 1: calculating the initial values, comprising a semi-major axis a.sub.0, an orbital period P.sub.t.sup.0 and an eccentricity e.sub.0 of an osculating orbit of the space object at a time t.sub.0; step 2: introducing a time transformation parameter δ to construct an orbital period transformation factor β; step 3: determining a number N of the uneven interpolation nodes in one orbital period by an interpolation method; step 4: constructing a set of uneven interpolation nodes {T.sub.l} from an initial node, and generating ephemeris values on the set of interpolation nodes by numerical integration; and step 5: performing interpolations on dense ephemeris times {t.sub.i} according to the set of uneven interpolation nodes {T.sub.l} and the ephemeris values of the set of interpolation nodes, so as to complete all ephemeris calculations.

2. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 1, wherein the ephemeris calculation of the space object corresponds to the solving of the dynamic equation as follows:
{right arrow over ({umlaut over (r)})}={right arrow over (F)}({right arrow over (r)},{right arrow over ({dot over (r)})},t)={right arrow over (F)}.sub.0({right arrow over (r)})+{right arrow over (F)}.sub.ε({right arrow over (r)},{right arrow over ({dot over (r)})},t)  (1) wherein, {right arrow over (F)} represents a force acting on the space object; {right arrow over (F)}.sub.0 represents a central gravity of the earth; {right arrow over (F)}.sub.ε represents perturbations; {right arrow over (r)} represents a position vector; {right arrow over ({dot over (r)})} represents a velocity vector; t represents a time variable; a position vector {right arrow over (r)}(t.sub.0)={right arrow over (r)}.sub.0 and a velocity vector {right arrow over ({dot over (r)})}(t.sub.0)={right arrow over ({dot over (r)})}.sub.0 of the space object at t.sub.0 are known; t.sub.i, i=1, 2, 3, . . . m represents a set of designated densely distributed time points; the calculation of the dense ephemeris comes down to calculating a position vector {right arrow over (r)}(t.sub.i)={right arrow over (r)}.sub.i and a velocity vector {right arrow over ({dot over (r)})}(t.sub.i)={right arrow over ({dot over (r)})}.sub.i of the space object at t.sub.i from the known {right arrow over (r)}.sub.0 and {right arrow over ({dot over (r)})}.sub.0; considering that t.sub.0 may be in an intermediate position of each range of t.sub.i, {t.sub.i} is divided into {t.sub.j.sup.b} and {t.sub.k.sup.c}; {t.sub.j.sup.b} comprises m.sub.j elements increasing monotonously under a subscript j, wherein m.sub.j≥0 and t.sub.0<t.sub.j.sup.b; {t.sub.k.sup.c} comprises m.sub.k elements decreasing monotonically under a subscript k, wherein m.sub.k≥0, t.sub.0>t.sub.k.sup.c, and m.sub.j+m.sub.k=m.

3. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 2, wherein in step 1, the initial values comprising the semi-major axis a.sub.0, the period P.sub.t.sup.0 and the eccentricity e.sub.0 of the osculating orbit of the space object at t.sub.0 are calculated as follows: $\begin{array}{cc}{a}_0=\frac{\mu r_0}{-{\stackrel{\underset{.fwdarw.}{.}}{r}}_0^2{r}_0+2\mu }& \left(2\right)\\ P_t^0=2\pi a_0\sqrt{a_0}/\sqrt{\mu }& \left(3\right)\\ e_0\sqrt{1-{.Math.{\stackrel{.fwdarw.}{r}}_0\times {\stackrel{\underset{.fwdarw.}{.}}{r}}_0.Math.}^2/\mu a_0}& \left(4\right)\end{array}$ wherein, μ represents a geocentric constant of gravitation; r.sub.0 is a modulus of {right arrow over (r)}.sub.0; {right arrow over (r)}.sub.0 and {right arrow over ({dot over (r)})}.sub.0 represent a position vector and a velocity vector of the space object at the time t.sub.0, respectively.

4. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 2, wherein in step 2, the time transformation parameter δ is introduced to express the orbital period transformation factor β, and β is calculated by performing a definite integration as follows: $\begin{array}{cc}\beta =\frac{1}{2\pi }{\int }_0^{2\pi }\frac{1}{{\left(1-e\cos E\right)}^{\delta }}\mathrm{dE}& \left(5\right)\end{array}$ wherein, β is a constant related to an orbital eccentricity e and δ; e is substituted by e.sub.0, and δ is 0.3; E represents an orbital eccentric anomaly.

5. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 4, wherein in step 3, by taking e and δ as discrete independent variables, a numerical experiment is carried out based on the selected interpolation method to determine a minimum value of N capable of meeting an orbital calculation accuracy requirement; a table is made, and in use, a value of N corresponding to known e and δ is determined by consulting the table.

6. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 2, wherein in step 4, the set of uneven interpolation nodes are constructed by taking t.sub.0 as the initial node; T.sub.0=t.sub.0, wherein T.sub.0 is a time corresponding to the initial node; ephemeris values {right arrow over (R)}.sub.0={right arrow over (r)}.sub.0 and {right arrow over ({dot over (R)})}.sub.0={right arrow over ({dot over (r)})}.sub.0 on the initial node are known; the set of uneven interpolation nodes {T.sub.l} are constructed according to a requirement for the interpolation based on the dense ephemeris times {t.sub.i}; the ephemeris values on the set of uneven interpolation nodes are generated along with the construction process.

7. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 6, wherein in step 4, the set of uneven interpolation nodes {T.sub.l} are specifically constructed as follows: considering that an even number q of interpolation nodes is required by an odd-order interpolation polynomial for one interpolation, the set of uneven interpolation nodes are constructed in three cases according to different distribution characteristics of the time points {t.sub.i}relative to t.sub.0: 1) if {t.sub.j.sup.b} and {t.sub.k.sup.c} are both non-empty, then m.sub.j>0, m.sub.k>0, and there are time points on both sides of t.sub.0; in this case, forward interpolation node construction is performed first; in a forward interpolation node construction process, if T.sub.l>t.sub.m.sub.j.sup.b appears for the first time, wherein t.sub.m.sub.j.sup.b is an element in {t.sub.j.sup.b}, then the construction continues q/2−1 steps along the forward direction, and then switches to backward interpolation node construction; in a backward interpolation node construction process, if T.sub.l<t.sub.m.sub.k.sup.c appears for the first time, wherein t.sub.m.sub.k.sup.c is an element in {t.sub.k.sup.c}, then the construction continues q/2−1 steps along the backward direction, and the entire interpolation node construction process ends; 2) if {t.sub.k.sup.c} is empty, then m.sub.k=0, and the time points {t.sub.i} are all located in the forward direction of t.sub.0; in this case, forward interpolation node construction is performed first; in the forward interpolation node construction process, if T.sub.l>t.sub.m.sub.j.sup.b appears for the first time, then the construction continues q/2−1 steps along the forward direction; then a number s, s≥0 of interpolation nodes in the backward direction of an interpolation node t.sub.1.sup.b is checked; if s≥q/2, then the entire interpolation node construction process ends; otherwise, the construction switches to backward interpolation node construction; in the backward node construction process, the construction goes q/2−s steps along the backward direction, and the entire interpolation node construction process ends; and 3) if {t.sub.j.sup.b} is empty, then m.sub.j=0, and the interpolation nodes {t.sub.i} are all located in the backward direction of t.sub.0; in this case, backward node construction is performed first; in the backward node construction process, if T.sub.l<t.sub.m.sub.k.sup.c, appears for the first time, then the construction continues q/2−1 steps in the backward direction; then a number s, s≥0 of nodes in the forward direction of an interpolation node t.sub.1.sup.c is checked; if s≥q/2, then the entire node construction process ends; otherwise, the construction switches to forward node construction; in the forward node construction process, the construction goes q/2−s steps in the forward direction, and the entire interpolation node construction process ends; these three interpolation node construction cases cover all the distribution characteristics of {t.sub.i} relative to t.sub.0, and uneven forward or backward interpolation nodes {T.sub.l} involved are derived by: assuming that an interpolation node T.sub.l and corresponding ephemeris values {right arrow over (R)}.sub.l and {right arrow over ({dot over (R)})}.sub.l are known, then calculating the semi-major axis and the orbital period of the osculating orbit respectively according to Eqs. (6) and (7): $\begin{array}{cc}{a}_l=\frac{\mu R_l}{-{\stackrel{\underset{.fwdarw.}{.}}{R}}_l^2{R}_l+2\mu }& \left(6\right)\\ P_t^l=2\pi a_l\sqrt{a_l}/\sqrt{\mu }& \left(7\right)\end{array}$ wherein, μ is a geocentric constant of gravitation; R.sub.l is a modulus of {right arrow over (R)}.sub.l; P.sub.t.sup.l is an orbital period measured by t and calculated from the ephemeris values of a l-th node; values of a.sub.0 and P.sub.t.sup.0, when l=0, are calculated in step 1; calculating an interval between T.sub.l and an adjacent interpolation node according to Eq. (8): $\begin{array}{cc}\Delta T_l={\beta \left(\frac{R_l}{a_l}\right)}^{1+\delta }\frac{P_t^l}{N}& \left(8\right)\end{array}$ constructing a forward node T.sub.l+1 as:
T.sub.l+1=T.sub.l+ΔT.sub.l  (9) constructing a backward node T.sub.l−1 as:
T.sub.l−1=T.sub.l−ΔT.sub.l  (10).

8. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 7, wherein in step 4, the ephemeris values on the interpolation node are specifically generated by a numerical integration as follows: the ephemeris values on the interpolation node are generated along with the construction process, and ephemeris values {right arrow over (R)}.sub.l+1, {right arrow over ({dot over (R)})}.sub.l+1 or {right arrow over (R)}.sub.l−1, {right arrow over ({dot over (R)})}.sub.l−1 on the node T.sub.l+1 or T.sub.l−1 are derived from T.sub.l to T.sub.l+1 or T.sub.l−1 by a single-step method for achieving variable-step integration.

9. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 7, wherein in step 5, for any ephemeris time t.sub.i, q/2 interpolation nodes are respectively taken from both sides of t.sub.i; ephemeris values on these nodes are used to interpolate by a selected interpolation method to obtain ephemeris values {right arrow over (r)}.sub.i and {right arrow over ({dot over (r)})}.sub.i at the time t.sub.i; then {t.sub.i} is traversed to complete all ephemeris calculations.

10. The method for accurately and efficiently calculating the dense ephemeris of the high-eccentricity orbit according to claim 9, wherein in step 5, the selected interpolation method is a Hermite interpolation method.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows a flowchart of the method for dense ephemeris computation.

(2) FIG. 2 shows interpolation errors of object positions within an orbital period based on even nodes.

(3) FIG. 3 shows interpolation errors of object velocities within an orbital period based on even nodes.

(4) FIG. 4 shows interpolation errors of object positions within an orbital period based on uneven nodes.

(5) FIG. 5 shows interpolation errors of object velocities within an orbital period based on uneven nodes.

(6) FIG. 6 shows a position-error comparison of Lagrange and Hermite methods in Case 2.

(7) FIG. 7 shows a velocity-error comparison of Lagrange and Hermite methods in Case 2.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(8) The present disclosure is described in more detail below with reference to the drawings.

(9) A dense ephemeris calculation method, as shown in FIG. 1, implements the ephemeris calculation of a space object by solving a dynamic equation as follows:
{right arrow over ({umlaut over (r)})}={right arrow over (F)}({right arrow over (r)},{right arrow over ({dot over (r)})},t)={right arrow over (F)}.sub.0({right arrow over (r)})+{right arrow over (F)}.sub.ε({right arrow over (r)},{right arrow over ({dot over (r)})},t)

(10) where, {right arrow over (F)} represents a synthetic force acting on the object; {right arrow over (F)}.sub.0 represents the central gravity of the earth; {right arrow over (F)}.sub.ε represents the perturbations. A position vector {right arrow over (r)}(t.sub.0)={right arrow over (r)}.sub.0 and a velocity vector {right arrow over ({dot over (r)})}(t.sub.0)={right arrow over ({dot over (r)})}.sub.0 of the object at an epoch t.sub.0 are known. t.sub.i (i=1, 2, 3, . . . m) represents a set of designated densely distributed time points. The calculation of the dense ephemeris comes down to calculating a position vector {right arrow over (r)}(t.sub.i)={right arrow over (r)}.sub.i and a velocity vector {right arrow over ({dot over (r)})}(t.sub.i)={right arrow over ({dot over (r)})}.sub.i of the object at t.sub.i from the known {right arrow over (r)}.sub.0 and {right arrow over ({dot over (r)})}.sub.0. Generally, since t.sub.0 may be between the minimum value of {t.sub.i} and the maximum value of {t.sub.i}, for convenience of description, {t.sub.i} is divided into {t.sub.j.sup.b} and {t.sub.k.sup.c}. {t.sub.j.sup.b} includes m.sub.j elements increasing monotonously under a subscript j, wherein m.sub.j≥0 and t.sub.0<t.sub.j.sup.b. {t.sub.k.sup.c} includes m.sub.k elements decreasing monotonically under a subscript k, wherein m.sub.k≥0, t.sub.0>t.sub.k.sup.c, and m.sub.j+m.sub.k=m.

(11) It is difficult to balance the calculation accuracy and calculation efficiency of the above-mentioned dense ephemeris of the high-eccentricity orbit through the existing interpolation methods. When the number of interpolation nodes is fixed, the interpolation nodes evenly spaced by time intervals may cause a reduction in the interpolation accuracy near the perigee due to the sparse interpolation nodes, and an increase in the interpolation accuracy near the apogee due to the dense interpolation nodes, which shows an imbalance problem. Increasing the number of interpolation nodes can solve this imbalance problem, but the numerical integration time will be prolonged, which is likely to cause low calculation efficiency. In order to balance the calculation accuracy and calculation efficiency of the ephemeris of the high-eccentricity orbit, the present disclosure proposes an uneven interpolation node construction method to coordinate with the Hermite interpolation polynomial to carry out dense ephemeris interpolation. This technical solution can be used in coordination with a general single-step integration method in the ephemeris calculation process, and can meet the high-accuracy and high-efficiency requirements of the dense ephemeris calculation of the high-eccentricity orbit. Based on an uneven interpolation node construction theory, the present disclosure proposes an interpolation node construction equation as:

(12) $\Delta T_l={\beta \left(\frac{R_l}{a}\right)}^{1+\delta }\frac{P_t}{N}$

(13) Correspondingly, an even interpolation node construction equation is;

(14) $\Delta T_l=T_{l+1}-T_l=\frac{P_t}{N}$

(15) The meaning of each symbol is described in the derivation process, and will not be repeated here. It needs to be pointed out that δ and N are constants predetermined before the interpolation nodes are constructed. δ and N affect the construction method of the interpolation nodes, and further affect the calculation accuracy and efficiency of the dense ephemeris points, so they need to be described in detail, but for convenience of description, the present disclosure only provides the numerical experimental results of δ and the determination method of N.

(16) 1. The value of δ determines the difference in the interpolation accuracy near the perigee and near the apogee. By reasonably selecting the value of 8 in a suitable value range of [−1,1], the ephemeris calculation error can be more uniformly distributed on the entire orbit. A large number of numerical experimental results show that the optimal universal value for all orbital eccentricities and various interpolation polynomials is 0.3 (see Table 1).

(17) 2. Under a certain interpolation accuracy requirement, corresponding to the known e and δ, there must be a minimum N capable of maximizing the calculation efficiency. Therefore, a numerical experiment is carried out to find the minimum N and make a table. The value of N can be determined by consulting the table.

(18) TABLE-US-00001 TABLE 1 Optimal δ and minimum N for position interpolations of different eccentricities by different interpolation methods (Runge-Kutta-Felhberg (RKF8 (7)) integrator, requiring position interpolation accuracy greater than 10.sup.−7 a.sub.e) Lagrange (9-degree) Newton (9-degree) Neville (9-degree) Hermite (7 degree) e δ N δ N δ N δ N 0.00 −1.0-1.0   16 −1.0-1.0   16 −1.0-1.0   16 −1.0-1.0   12 0.05 −0.2-1.0   19 −0.2-1.0   19 −0.2-1.0   19 −0.4-1.0   13 0.10 0.7-1.0 22 0.7-1.0 22 0.7-1.0 22 0.1-1.0 14 0.15 0.5-1.0 24 0.5-1.0 24 0.5-1.0 24 −0.1-1.0   16 0.20 0.3-1.0 26 0.3-1.0 26 0.3-1.0 26 −0.1-1.0   17 0.25 0.5-0.9 28 0.5-0.9 28 0.5-0.9 28 0.2-0.9 19 0.30 0.2-0.4 31 0.2-0.4 31 0.2-0.4 31 0.1-0.8 20 0.35 0.4-0.5 40 0.4-0.5 40 0.4-0.5 40 0.3-0.5 21 0.40 0.2-0.6 33 0.2-0.6 33 0.2-0.6 33 0.2-0.4 23 0.45 0.3-0.4 38 0.3-0.4 38 0.3-0.4 38 0.2-0.4 24 0.50 0.3-0.5 41 0.3-0.5 41 0.3-0.5 41 0.3 26 0.55 0.3 43 0.3 43 0.3 43 0.3 27 0.60 0.3 46 0.3 46 0.3 46 0.3 29 0.65 0.3 50 0.3 50 0.3 50 0.3 32 0.70 0.3 54 0.3 54 0.3 54 0.3 35 0.75 0.3 59 0.3 59 0.3 59 0.3 40 0.80 0.3 65 0.3 65 0.3 65 0.3 43 0.85 0.3 77 0.3 77 0.3 77 0.3 53 0.90 0.3 115 0.3 115 0.3 115 0.3 54 Note: a.sub.e = 6378136 m. The optimal δ for the smaller-eccentricity orbit is given in a range, and any δ in this range can lead to the minimum N capable of meeting the accuracy requirement.

(19) The technical solution of dense ephemeris calculation is described below with reference to an embodiment.

(20) Embodiment: The orbital elements of a space object at an epoch t.sub.0=58362.0 (modified Julian date (MJD)) are known: the semi-major axis a=5.25a.sub.e (a.sub.e is an equatorial radius of a reference ellipsoid, which is 6378136 m), the eccentricity e=0.8 (the values of a and e in Case 3 are different from these values, see the description of Case 3 for details), the orbital inclination angle: i=45°; the right ascension of ascending node: Ω=0°; the argument of perigee: ω=0°; the mean anomaly: M=0°. The orbital elements can be converted to obtain a position vector {right arrow over (r)}.sub.0 and a velocity vector {right arrow over ({dot over (r)})}.sub.0 at the epoch t.sub.0, and the conversion method will not be repeated here. The time points to generate the dense ephemeris are {t.sub.i}(i=1, 2, 3, . . . m), where t.sub.1=t.sub.0+P.sub.0/2, t.sub.m=t.sub.0+3P.sub.0/2, lasting one orbital period P.sub.0 and the interval between intermediate points is 1 s.

(21) Case 1: The even time interval method (commonly used in navigation satellite ephemerides) and the uneven interval method proposed by the present disclosure are used to construct interpolation nodes, wherein the number N of interpolation nodes is 80, and the optimal universal value of δ is 0.3. After that, the 9.sup.th-degree Lagrange polynomial is used for interpolation, so as to compare the effects of the two interpolation node construction methods on the ephemeris calculation accuracy.

(22) Case 2: δ is set as the optimal universal value of 0.3. The uneven interpolation node construction method proposed by the present disclosure in coordination with the 9.sup.th-degree Lagrange polynomial (N=114) and the 7.sup.th-degree Hermite polynomial (N=60), respectively, is used for interpolation, so as to compare the ephemeris calculation accuracy of the two interpolation methods.

(23) Case 3: δ is set as the optimal universal value of 0.3. A point-to-point integration method and the uneven node interpolation method (using the Hermite interpolation polynomial) proposed by the present disclosure are used to calculate the ephemeris of orbits with different eccentricities, wherein the eccentricity e increases from 0 to 0.9 by an increment of 0.05, a is obtained by the constraint a(1−e)=1.05a.sub.e; the other orbital parameters remain unchanged, and a total of 19 sets of orbital elements are generated. The two methods are compared in the calculation times of the integral right function and the calculation time by a central processing unit (CPU).

(24) The dense ephemeris calculation in the embodiment comprehensively integrates dynamic equation solving and ephemeris interpolating two technologies. The position vector {dot over (r)}.sub.0 and velocity vector {right arrow over ({dot over (r)})}.sub.0 at t.sub.0 are known, and dynamic models used for the orbital calculation include: the earth's gravitational field being GRIM4_S4 model (8×8); the atmospheric density model being DTM94 (disclosed by C. Berger, R. Biancale, M. Ill, et al. Improvement of the empirical thermospheric model DTM: DTM94—a comparative review of various temporal variations and prospects in space geodesy applications [J]. Journal of Geodesy, 72(3):161-178); and the position calculation of the sun and moon being Jean Meeus' analysis formula (disclosed by Meeus J. Astronomical formulae for calculators [M]. Tweede Druk, 1979). During interpolation, the ephemeris values on the interpolation node are generated by RKF8 (7) integrator capable of achieving variable-step integration. The truncation error of each step of integration is controlled to be 10.sup.−10a.sub.e for position and 10.sup.−10 a.sub.e/day for velocity. In addition, in order to evaluate the interpolation error, the standard ephemeris values {right arrow over (r)}.sub.i.sup.s and {right arrow over ({dot over (r)})}.sub.i.sup.s on {t.sub.i} are obtained by variable-step point-to-point integration by RKF8 (7). In Case 1 and Case 2, the truncation error of the point-to-point integration is controlled to be 10.sup.−10a.sub.e for position and 10.sup.−10 a.sub.e/day for velocity. In Case 3, the truncation error of point-to-point integration is controlled to be 10.sup.−8a.sub.e for position and 10.sup.−7 a.sub.e/day for velocity. The interpolation errors of the position vector {right arrow over (r)}.sub.i and the velocity vector {right arrow over ({dot over (r)})}.sub.i at each ephemeris point are:

(25) $\{\begin{array}{c}\Delta {\stackrel{.fwdarw.}{r}}_i=.Math.{\stackrel{.fwdarw.}{r}}_i-{\stackrel{.fwdarw.}{r}}_i^s.Math.\\ \Delta {\stackrel{\underset{.fwdarw.}{.}}{r}}_i=.Math.{\stackrel{\underset{.fwdarw.}{.}}{r}}_i-{\stackrel{\underset{.fwdarw.}{.}}{r}}_i^s.Math.\end{array}\hspace{1em}$

(26) The even interpolation node construction method involved in Case 1 has nothing to do with the technical solution of the present disclosure, and its implementation process and steps will not be described in detail here. The following only elaborates on the technical solution proposed by the present disclosure:

(27) Step 1: Calculating initial values, including a semi-major axis a.sub.0, an orbital period P.sub.t.sup.0 and an eccentricity e.sub.0 of an osculating orbit of a space object at a time t.sub.0 as follows.

(28) $a_0=\frac{\mu r_0}{-{\stackrel{\underset{.fwdarw.}{.}}{r}}_0^2{r}_0+2\mu }$$P_t^0=2\pi a_0\sqrt{a_0}/\sqrt{\mu }$$e_0\sqrt{1-{.Math.{\stackrel{.fwdarw.}{r}}_0\times {\stackrel{\underset{.fwdarw.}{.}}{r}}_0.Math.}^2/\mu a_0}$

(29) where, μ is a geocentric constant of gravitation, which is 0.39860043770442×10.sup.15 m.sup.3/s.sup.2 In this embodiment, the orbital eccentricity of the space object is 0.8, which indicates that the orbit is a high-eccentricity orbit.

(30) Step 2: Introducing a transformation parameter δ to express an orbital period transformation factor β, and obtaining β by performing a definite integration as follows:

(31) $\beta =\frac{1}{2\pi }{\int }_0^{2\pi }\frac{1}{{\left(1-e\cos E\right)}^{\delta }}\mathrm{dE}$

(32) In this embodiment, the optimal universal value of δ is 0.3, and the value of eccentricity e is substituted by e.sub.0. β can be calculated by using any definite integration method, and the Simpson method is used in this embodiment.

(33) Step 3: Determining the number N of uneven interpolation nodes in one orbital period.

(34) The number N of interpolation nodes has been given in Case 1 and Case 2, so no calculation is needed therein. In Case 3, δ is set as the optimal universal value of 0.3, and the Hermite 7.sup.th-degree polynomial is used for interpolation. For the 19 sets of eccentricities e, the minimum N capable of meeting the ephemeris calculation accuracy requirement is determined through a numerical experiment (see Table 2)

(35) TABLE-US-00002 TABLE 2 Numerical experimental results of minimum N for orbits of different eccentricities in Case 3 (Requiring position interpolation accuracy greater than 10.sup.−7 a.sub.e and velocity interpolation accuracy greater than 10.sup.−6 a.sub.e/day) e 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 N 18 19 22 24 26 28 30 33 35 37 42 43 46 50 55 59 65 73 75

(36) Step 4: Constructing a set of uneven interpolation nodes {T.sub.l} from an initial node t.sub.0, and generating ephemeris values on the set of interpolation nodes by numerical integration.

(37) In this embodiment, all the ephemeris points {t.sub.i} are located on the right of t.sub.0, so forward interpolation node construction is performed first. t.sub.0 is taken as an initial interpolation node, and ephemeris values at the initial interpolation node are known, T.sub.0=t.sub.0, {right arrow over (R)}.sub.0={right arrow over (r)}.sub.0, {right arrow over ({dot over (R)})}.sub.0={right arrow over ({dot over (r)})}.sub.0. The rest of the interpolation nodes are constructed recursively as follows.

(38) Assuming that the l-th node T.sub.l and corresponding ephemeris values {right arrow over (R)}.sub.l and {right arrow over ({dot over (R)})}.sub.l thereon are known, then the semi-major axis and orbital period of the osculating orbit are calculated according to the following equations:

(39) 0$a_l=\frac{\mu R_l}{-{\stackrel{\underset{.fwdarw.}{.}}{R}}_l^2{R}_l+2\mu }$$P_t^l=2\pi a_l\sqrt{a_l}/\sqrt{\mu }$

(40) where, R.sub.l is the modulus of {right arrow over (R)}.sub.l. When l=0, a.sub.0 and p.sub.t.sup.0 are calculated in Step 1. According to the known N, an interval between the l-th node T.sub.l and an (l+1)-th node T.sub.l+1 is calculated as follows:

(41) $\Delta T_l={\beta \left(\frac{R_l}{a_l}\right)}^{1+\delta }\frac{P_t^l}{N}$

(42) The (l+1)-th node can be derived according to the following equation:
T.sub.l+1=T.sub.l+ΔT.sub.l

(43) By taking the ephemeris values {right arrow over (R)}.sub.l and {right arrow over ({dot over (R)})}.sub.l on the l-th node as initial values, a variable-step integration is performed through RKF8(7) from T.sub.l to T.sub.l+1 to obtain the ephemeris values {right arrow over (R)}.sub.l+1 and {right arrow over ({dot over (R)})}.sub.l+1 on the (l+1)-th node. At the same time, {right arrow over ({umlaut over (R)})}.sub.l={right arrow over (F)} ({right arrow over (R)}.sub.l, {right arrow over ({dot over (R)})}.sub.l, T.sub.l) is calculated. In the forward recursion process, every time a new interpolation node T.sub.l is generated, it is necessary to determine whether T.sub.l>t.sub.m is satisfied. If not, the forward construction process continues. If yes, (q/2−1) interpolation nodes are constructed forward based on the current interpolation node, and the forward construction process ends. Then the number s of nodes on the left of t.sub.1 is checked. Because t.sub.1 and t.sub.0 are separated by half an orbital period, and s>>q/2, there is no need to carry out backward node construction, and the entire interpolation node construction process ends.

(44) When the 9.sup.th-degree Lagrange polynomial is used for interpolation, the number q of nodes required for one interpolation is 10. When the 7.sup.th-degree Hermite interpolation polynomial is used, the number q of nodes required for one interpolation is only 4.

(45) Step 5: performing interpolations on dense time points {t.sub.i} according to the set of uneven interpolation nodes {T.sub.l} and the ephemeris values thereof, so as to complete all ephemeris calculations.

(46) For any time point t.sub.i, q/2 interpolation nodes on the left and right of t.sub.i, respectively, are taken, and the ephemeris values {right arrow over (R)}.sub.l and {right arrow over ({dot over (R)})}.sub.l on these nodes are used for interpolation to obtain ephemeris values {right arrow over (r)}.sub.i and {right arrow over ({dot over (r)})}.sub.i at t.sub.i ({right arrow over ({umlaut over (R)})}.sub.l is further used in case of the Hermite interpolation polynomial). Similarly, {t.sub.i} is traversed to complete all ephemeris calculations by the above interpolation process.

(47) The calculation results of Cases 1 to 3 are shown in FIGS. 2 to 7 and Table 3.

(48) FIGS. 2 to 5 show the calculation results of the two interpolation node construction methods in one orbital period. According to the interpolation results based on evenly spaced nodes in FIGS. 2 and 3, the position and velocity errors near the perigee of the orbit have great fluctuations. The maximum position error reaches an order of 10.sup.−2, that is, the position error can reach at least tens of kilometers, and the maximum velocity error can reach an order of 10°, that is, the velocity error is greater than 100 m/s. FIGS. 4 and 5 show the interpolation results based on unevenly spaced nodes constructed according to the technical solution of the present disclosure. The position and velocity errors are respectively above an order of 10.sup.−7 and 10.sup.−5 at any time within one orbital period, that is, the position error does not exceed 1 m, and the velocity error does not exceed 1 mm/s. The error comparison results show that the ephemeris calculation method based on unevenly spaced interpolation nodes is superior to the navigation satellite ephemeris calculation method based on evenly spaced interpolation nodes and suitable for the calculation of high-eccentricity orbits, and can significantly improve the calculation accuracy without reducing the calculation efficiency.

(49) FIGS. 6 and 7 show changes in the position error and velocity error of the dense ephemeris obtained by using the 9.sup.th-degree Lagrange (N=114) and 7.sup.th-degree Hermite (N=60) interpolation methods within one orbital period. It can be seen from the comparison that the calculation errors of the two interpolation methods are roughly equivalent, which indirectly shows the superiority of the Hermite method, that is, the Hermite method can still guarantee the calculation accuracy with less interpolation nodes (meaning high calculation efficiency).

(50) Table 3 demonstrates the operation results of the point-to-point integration method and the technical solution of the present disclosure on the same computer. It can be seen from the table that compared to the point-to-point integration method, although the local truncation error control accuracy of the present technical solution is higher, because the interpolation nodes are sparsely distributed on the orbit, the calculation times of the integral right function and the integration time are both greatly reduced. The time taken by this technical solution of the present disclosure is concentrated on the operation of interpolations, which is a simple algebraic operation with extremely high calculation efficiency. Therefore, the overall calculation time of the present disclosure is greatly reduced compared to the point-to-point integration method. In Case 3, the calculation efficiency of the two methods differs by about 30 to 100 times, indicating that the dense ephemeris calculation method proposed by the present disclosure can greatly improve the calculation efficiency while ensuring a certain calculation accuracy.

(51) This embodiment illustrates the superiority of the technical solution of the present disclosure from multiple aspects. The uneven interpolation node construction method based on the time transformation theory can significantly improve the overall interpolation accuracy of the high-eccentricity orbit, and it can be used in coordination with the Hermite method to further improve the calculation efficiency. Finally, the present disclosure illustrates the efficiency of the technical solution by comparing it with the point-to-point integration method.

(52) TABLE-US-00003 TABLE 3 Comparison of calculation efficiency between point-to-point integration and the technical solution of the present disclosure Point-to-point integration Integration + interpolation Calculation CPU Calculation CPU times of integral calculation times of integral calculation e right function time (s) right function time (s) 0.00 71187 1.125000 1283 0.031250 0.05 76843 1.218750 1219 0.031250 0.10 83317 1.265625 1233 0.031250 0.15 90675 1.421875 1287 0.046875 0.20 99384 1.578125 1391 0.062500 0.25 109447 1.750000 1495 0.031250 0.30 121379 1.953125 1599 0.046875 0.35 135590 2.234375 1599 0.046875 0.40 152891 2.187500 1597 0.046875 0.45 174135 2.781250 1697 0.062500 0.50 200875 2.968750 1675 0.046875 0.55 235234 3.156250 1777 0.062500 0.60 280605 4.203125 1766 0.062500 0.65 342797 5.515625 1831 0.078125 0.70 431886 6.531250 1844 0.109375 0.75 567658 8.531250 2063 0.109375 0.80 793182 11.812500 2208 0.109375 0.85 1220986 17.093750 2325 0.187500 0.90 2242747 31.671875 3190 0.343750

(53) The above described are only preferred implementations of the present disclosure, and the scope of the present disclosure is not limited thereto. All technical solutions based on the idea of the present disclosure should fall within the claimed scope of the present disclosure. It should be pointed out that for a person of ordinary skill in the art, several improvements and modifications made without departing from the principle of the present disclosure should be deemed as falling within the protection scope of the present disclosure.