Abstract
A scalable second-order analytical orbit propagator programme based on modern and classical perturbation methods is being developed. As a first step in the validation and verification of part of our orbit propagator programme, we only consider the perturbation produced by zonal harmonic coefficients in the Earth’s gravity potential, so that it is possible to analyze the behaviour of the mathematical expressions involved in Delaunay normalization and the Krylov-Bogoliubov-Mitropolsky method in depth and determine their limits.
1. Introduction
The collision between the Iridium and Cosmos satellites in 2009 has demonstrated, among other things, the necessity of improving orbit prediction methods. With this purpose in mind, a scalable second-order analytical orbit propagator programme (AOPP) is being developed. This AOPP combines modern perturbation methods based on Lie transforms and classical averaging techniques, depending on the orbit types or the requirements needed for a space mission, such as catalogue maintenance operations or long-period evolution, for example.
Most of the analytical theories in the artificial satellite problem start by removing the short-period terms [1–4]. In fact, the most frequently used AOPPs, SGP4 (simplified general perturbations 4) [5, 6] and PPT2 (position and partials as functions of time 2) [7, 8], are developed from the Brouwer-Lyddane theory. However, other analytical theories, like AOP [9], simplify the problem by means of a Lie transform known as the elimination of the parallax [10]. In this work, we analyze other alternatives to the Brouwer approach, which may arise from the elimination of the parallax. These are based on removing the long-period terms by means of the elimination of the perigee [11]. The transformed Hamiltonian has one degree of freedom and only depends on the radial distance and . Classically, this Hamiltonian is known as the radial intermediary [10]. Finally, the elimination of the short-period terms can be made using the classical Delaunay normalization (DN) [12] or by means of the Krylov-Bogoliubov-Mitropolsky method (KBM) [13, 14].
As a first step in the validation and verification of part of our AOPP, which we will call , we only consider the perturbation produced by zonal harmonic coefficients in the Earth’s gravity potential, such that it is possible to analyze the accuracy of the mathematical expressions involved in the corresponding analytical theories in depth, and thus determine its limits when classical problems, like critical inclination or small eccentricities and inclinations, appear.
The purpose of this study is to perform a comprehensive evaluation from the point of view of accuracy and the operative utility of the Delaunay normalization and Krylov-Bogoliubov-Mitropolsky method once the long-period terms have been removed. A methodology based on an exploratory data analysis (EDA) [15] is proposed in order to validate both integration methods and determine the values of the initial conditions where these methods are valid, in function of the duration of propagation. A space catalogue with two-line elements (TLE) is used as simulated data in order to validate our AOPP. After that, we identify the outliers and make a general study of the rest of the TLEs considering the zonal harmonic coefficients . Then we proceed to study the TLEs near critical inclination ( or ) and analyze the modification implemented in PPT2, HANDE (Hoots analytic dynamic ephemeris) [16], and AOP [9] for and . Finally, we summarize this study.
2. Analytical Theories
In polar-nodal variables , the Hamiltonian for an Earth satellite perturbed by zonal harmonic potential terms is given by This Hamiltonian defines a two-degree-of-freedom dynamical system. The coordinate is the radial distance from the Earth’s center of mass to the satellite, represents the argument of latitude ( is the true anomaly and is the argument of the perigee), and is the argument of the node, which is ignored in . The momentum is the radial velocity of the orbiter, whereas is the module and is the third component of its angular momentum. The variables , , and are the momenta for , , and , respectively. is the Legendre polynomial of degree , is the sine of the inclination of the orbit, is the maximum order of the zonal harmonic perturbation being considered, is the gravitational constant, is the equatorial radius of the planet, and are the zonal harmonic coefficients. More details about the Hamiltonian formulation of the artificial satellite problem can be found in [10, 17].
Figure 1 shows the combinations of Lie transforms and classical averaging techniques used in developing the two analytical theories, Delaunay normalization and Krylov-Bogoliubov-Mitropolsky method. Both theories begin by removing the long-period terms caused by the argument of the perigee through the combination of two Lie transforms, the elimination of the parallax, and the elimination of the perigee (Lie transforms and , resp.) from Hamiltonian (1). The transformed Hamiltonian has one degree of freedom but maintains the short-period terms caused by the mean anomaly. These two Lie transforms are carried out in a closed form of eccentricity and inclination, which leads us to an integrable problem in the variables .
The second order of the Hamiltonian, after the elimination of the perigee, takes the following form: where the small parameter represents , , are polynomials in eccentricity, and their coefficients are polynomials in .
Traditionally, to complete the theory and obtain the mean elements, as Brouwer [2] did, a further reduction is made through Delaunay variables by means of the Delaunay normalization [12] (Lie transform ), which averages the problem over the mean anomaly. However, if the time and variable are replaced by the perturbed true anomaly and , which is defined as the inverse of , in the Hamiltonian ( transform), the equations of motion become a one-dimensional perturbed harmonic oscillator in which the Krylov-Bogoliubov-Mitropolsky method [13, 14, 18] can also be used to integrate the equations of motion [19–24]. It is worth noting that Delaunay’s variables are singular for zero eccentricity and/or zero inclination.
In this study, the analytical theories are carried out maintaining the physical parameters in symbolic form, such that they can be valid for any gravity field model. On the other hand, no change of variables is introduced in these theories so as not to remove any problems generated by small eccentricities and inclinations [25], including critical inclination, which allows us to evaluate their real impact on the analytical expressions and so determine the region of the space where these problems may appear. However, in the case of critical inclination, the solution implemented in PPT2 is analyzed.
Our AOPP combines these two theories. The analytical expressions and the orbit propagator, coded in C programming language, are performed by the symbolic-numeric environment MathATESAT [23].
3. Methodology
In order to evaluate the sensitivity of the two analytical integration techniques in the zonal case for certain initial conditions, we use the following basic methodology. Each one of the two integration methods is used to propagate a TLE space catalogue over a time span of days. Both propagations are compared with the integration of the original problem using a highly accurate 8th order Runge-Kutta method [26]. Then, the distance, along-track, cross-track, and radial errors are calculated. Finally, an exploratory data analysis consisting of the study of these errors is made. The details of this methodology are given below.
3.1. Study Data
Data from a space catalogue with TLEs are considered so as to provide a comprehensive knowledge of the DN and KBM performances. Although it is well known that TLEs have been designed to be used in combination with the SGP4 orbit propagator, we consider that a TLE space catalogue contains a large and representative number of different types of orbits, which can be considered as a reliable and independent test for our study.
The orbit data can be separated from the space catalogue. Attending to the frequency distributions of , , and , we can conclude that, in general, they are distributed uniformly between and , whereas the distributions of , , and are not uniform. The values of the semimajor axis are between and Earth Radii, despite the fact that of the objects have Earth Radii. The eccentricity is between and , although almost have . The inclination is between and , including objects near critical inclination. Near-circular low-Earth, medium-Earth, eccentric, and low-Earth orbits represent of the orbit types belonging to the considered catalogue.
3.2. Exploratory Data Analysis (EDA)
The TLEs are propagated over a time span of days using the DN, KBM, and Runge-Kutta methods. Then, the distance, along-track, cross-track, and radial errors are calculated and analyzed over two time spans of and days. After that, we carry out an EDA consisting of the following.(i)A graphic study of DN and KBM errors uses the box and whisker plot. This plot of a data set consists of a box drawn around the median value, where the lower and upper box edges bound the first and third quartiles (Q1 and Q3), respectively. The lower and upper whiskers extending from each bound box are the minimum and maximum values without outliers (circles) and extreme values (asterisks), respectively, whereas and represent the upper and lower limits, respectively. We must remark that in this study the extreme values are also considered as outliers.(ii)Scatter plots aim to characterize the orbits in which the propagators have a poorer overall performance (greater presence of outliers—as above ). A study of contingency tables (considering the number of outliers and corresponding orbits as variables) and the calculation of Yule’s Q, a measure of the magnitude of association between the two rates, complete the study chart.(iii)Excluding the atypical values, the behaviour of the two propagators and the distance, along-track, cross-track, and radial errors are explored. In particular, TLEs whose distance errors are above 1000 m are considered atypical values for a time span of 30 days.(iv)Finally, we explore the results of the proposed improvements in the two estimates for the orbits near critical inclination.
It is worth noting that exploratory data analysis does not require the incorporation of any prior knowledge of this process.
4. DN versus KBM for Model
In this section, the force model taken into account in the analytical orbit propagator only contains the zonal harmonic coefficients , as in SGP4.
4.1. Time Span of 7 Days
The study begins by performing a box and whisker plot analysis in order to identify potential outliers for DN and KBM distance error values. The results of this analysis are shown in Figure 2. As can be seen, there are no outliers below the lower-whisker limit, but if there were, they would represent very good behaviour for the analytical methods in the corresponding TLE. Outliers larger than m have not been included in this graph. The values of the median, minimum, maximum, first and third quartiles, and the upper-whisker limit are shown in Table 1 including the data which generate overflows.
The number of outliers detected in the DN and KBM cases is and , respectively, of which only 55 are common. DN outliers represent of the TLEs in which the eccentricity is below 0.01. Thus, this strong association is confirmed by a highly significant Yule’s Q value of 0.97. When the distance error is obtained from TLEs with and near critical inclination (between and ), DN outliers represent of the corresponding TLEs. In this case, Yule’s Q of also indicates a very strong association. For KBM, outliers also represent of TLEs near critical inclination (between and ), with a Yule’s Q of , whilst of the remaining outliers correspond to TLEs with inclinations of below and , in which Yule’s Q is . It is worth noting that the maximum distance error in these cases is below m. The distributions of the eccentricity versus inclination of the full TLE catalogue and of the TLEs which correspond to DN and KBM distance error outliers are depicted in Figure 3.
(a) Full catalogue
(b) DN outliers
(c) KBM outliers
Next, the behaviour of the DN and KBM methods, starting from distance, along-track, cross-track, and radial errors, is analyzed.
The histogram of DN distance error is shown in Figure 4(a). It should be noted that the lack of robustness of estimates is due to the asymmetry on the right of the histogram. This asymmetry corresponds to TLEs with small eccentricity, , and near critical inclination. Figure 4(b) shows a box and whisker plot of the distance, along-track, cross-track, and radial errors without upper outliers. These data are classified into three sets: , , and , and any other data. The most influential errors for TLEs with are found in the along-track and radial components, whereas for TLEs with and the worst error behaviour is found to be in the cross-track component.
(a) DN distance error histogram
(b) Distance, along-track, cross-track, and radial errors
Figure 5(a) depicts the histogram of KBM distance error, which exhibits a fairly symmetric behaviour. It is necessary to bear in mind the robustness of estimates. Figure 5(b) shows a box and whisker plot of the distance, along-track, cross-track, and radial errors without upper outliers. These data are classified as mentioned in the DN case. However, in this case, the most influential errors for TLEs with are found in the along-track and cross-track components, whereas for TLEs with and the worst error behaviour is found to be in the cross-track component, as in the DN case. Note the scarce influence of radial error in the three cases.
(a) KBM distance error histogram
(b) Distance, along-track, cross-track, and radial errors
In general, the worst behaviour is found in the category of the Delaunay normalization. Figure 6 shows box and whisker plot analysis, without outliers, for semimajor axis, eccentricity, inclination, and argument of the perigee errors for both DN and KBM methods in this category. The graphics, on the one hand, of the argument of the perigee and mean anomaly errors and, on the other hand, of the inclination and argument of the node errors are similar for DN and KBM methods. It is worth noting that the semimajor axis error is only slightly better in DN than in KBM cases, that is, m versus m in mean values, respectively. Eccentricity, argument of the perigee, and mean anomaly errors are much better in KBM than in DN cases. The worst eccentricity and mean anomaly determinations using the Delaunay normalization method, in this case, explain the influence of the radial error over the cross-track error (Figures 6(b) and 6(d)).
(a) Semimajor axis errors
(b) Eccentricity errors
(c) Inclination errors
(d) Argument of the perigee errors
Finally, in this section some of the KBM outliers are analyzed, in particular, those corresponding to TLEs with and . Both DN and KBM methods exhibit similar behaviour, but when this set of outlier values is studied KBM is slightly better than DN, 0.12 m versus 0.14 m in mean values, respectively, as can be seen in Figure 7.
(a) DN versus KBM distance error
(b) Box and whisker plot
4.2. Time Span of 30 Days
The box and whisker plot analysis is used in order to identify potential outliers for DN and KBM distance error values. Results of this analysis are shown in Figure 8. Outliers larger than m are not included in this graph. As can be seen, there are no outliers below the lower-whisker limit. The values of the median, minimum, maximum, first and third quartiles, and the upper-whisker limit are shown in Table 2, including the data which generate overflows.
The number of outlier values detected in the DN and KBM cases is and , respectively, of which are common. DN outliers represent of the TLEs with and near critical inclination. Thus this association is confirmed by a Yule’s Q value of . When the distance error is obtained from TLEs with and near critical inclination, the outliers represent of the corresponding TLEs, in which case Yule’s Q of also indicates a significant association. For KBM, outliers also represent of TLEs near critical inclination with a significant Yule’s Q of . The distributions of the eccentricity versus inclination of the full TLE catalogue and of the TLEs which correspond to DN and KBM outliers with distance error values above m are depicted in Figure 9. The number of outliers for a distance error above m detected in the DN and KBM cases is and , respectively, of which are common.
(a) Full catalogue
(b) DN outliers
(c) KBM outliers
Next, the behaviour of the DN and KBM methods, starting from distance, along-track, cross-track, and radial errors, is analyzed for TLEs with distance errors below m.
The histogram of DN distance error is shown in Figure 10(a). Figure 10(b) shows a box and whisker plot of the distance, along-track, cross-track, and radial errors. These data are classified as mentioned in the previous section. The most influential error for TLEs is found in the along-track component in all categories. The cross-track component is greater than the radial component for TLEs with and , whereas in the other two categories the radial component is greater than the cross-track component.
(a) DN distance error histogram
(b) Distance, along-track, cross-track, and radial errors
Figure 11(a) depicts the histogram of KBM distance error. Figure 11(b) shows a box and whisker plot of the distance, along-track, cross-track, and radial errors. The behaviour of the three error components is the same as in the DN case.
(a) KBM distance error histogram
(b) Distance, along-track, cross-track, and radial errors
The worst behaviour of the two methods is found for eccentricities below , although the behaviour is more robust in the case of DN than in KBM, as the asymmetries on the right of the histograms in Figures 10(a) and 11(a) show. Then a detailed study is presented to determine the range of influence of the eccentricity over the distance error. Catalogue data are classified into five sets: , , , , and , which represent , , , , and of the catalogue, respectively. Figure 12 shows a box and whisker plot analysis of the DN and KBM distance error. As can be seen, DN is worse than KBM for in particular, m versus m in mean values, respectively, for , and m versus m in mean values, respectively, for .
5. Outliers Near Critical Inclination
The number of TLEs with in the space catalogue considered in this study is , of which have eccentricities below , and thus they are not taken into account with the Delaunay normalization method. Technically, the critical inclination singularity is directly related to the denominator , which appears in the direct and inverse transformations of the elimination of the perigee. We must point out that this Lie transform is not valid in the neighborhood of and .
In order to reduce the impact of the singularity over this transform, we replace the terms with where and . This solution is implemented in other analytical orbit propagator programmes (PPT2, HANDE, and AOP). DN2 and KBM2 will refer to the errors obtained with the modified analytical orbit propagator programme. It is worth noting that this modified propagator does not integrate Hamiltonian (1) exactly but is an approximation, which must be identified and its phase space needs to be characterized. The preliminary numerical analysis is here presented over a time span of days. The first subsection studies the behaviour of DN2 and KBM2 near critical inclination in the case of perturbations, while in the second subsection the Krylov-Bogoliubov-Mitropolsky behaviour is analyzed only for perturbations.
5.1. Perturbations
Data considered in this study correspond to TLEs with and . Several of these TLEs have distance errors above m, and overflow values are obtained in two cases. Table 3 shows through statistical parameters that the behaviour of DN and KBM distance errors is similar in the considered TLEs.
Table 4 shows the values of the median, minimum, maximum, first and third quartiles, and the upper-whisker limit in the cases of DN2 and KBM2 distance errors for all TLEs with and , including the two data which generate overflows. Figure 13 shows that the behaviour in both cases is very similar. The maximum error is m.
(a) DN2 versus KBM2 distance error
(b) Box and whisker plot of DN2 versus KBM2 distance error
Tables 5 and 6 show that robustness is not obtained from a general improvement on estimates but from a proper calculation of original outliers. In fact, in of the cases, the behaviour of both estimators is worse, but original outlier estimates have improved. estimates obtained by the modified models in the nonoutlier cases are worse, between m and m in mean values, while the remaining outliers and the two overflows have reduced to m in mean value.
Figure 14 shows the box and whisker plot of DN2 and KBM2 distance errors; in this case the cross-track error is also the most influential as was seen in the previous general study for TLEs with and .
5.2. Perturbations
As can be seen in the previous subsection, the behaviour of the two integration methods is similar for the near critical inclination case. Therefore, we will only consider the Krylov-Bogoliubov-Mitropolsky method to conduct the study when perturbations caused by zonal harmonic coefficients from to are taken into consideration. We will refer to this model as KBMJ12, and KBMJ12_2 will be the modified analytical orbit propagator programme when (3) is taken into account.
The same previously studied TLEs are considered here. Several of these TLEs have distance errors above m, and overflow values are obtained in ten cases. The values of the median, minimum, maximum, first and third quartiles, and the upper-whisker limit in the case of KBMJ12 distance error are , , , , , and , respectively. Figure 15(a) shows a box and whisker plot of the KBMJ12 distance error in which only out of outliers were included.
(a) KBMJ12 distance error
(b) KBMJ12_2 distance error
The values of the median, minimum, maximum, first and third quartiles, and the upper-whisker limit in the case of KBMJ12_2 distance error for all TLEs, including the ten data which generate overflows, are , , , , , and , respectively. Figure 15(b) shows a box and whisker plot for KBMJ12_2 distance error and as can be observed there are no outlier values.
Although the ephemeris generated by the second-order theory in the case provides good accuracy, this becomes worse when the number of zonal harmonics increases. The same accuracy is only achieved if third-order effects are included. The estimation in cases, including the ten overflow values, is better in KBMJ12_2, being m in mean value, although it is worse in cases, with mean values of m and m, as can be seen in Table 7. Another difference with respect to is that in the main error is obtained in the along-track component, as can be seen in Figure 16.
Finally, a detailed study is presented to determine the range of influence of critical inclination. As might be expected, the behaviour of the analytical theory becomes worse in the neighborhood of near inclinations. This study is restricted to the interval , which has been divided into five subintervals. Figure 17 shows the distribution of the KBMJ12 distance error in the five subintervals. This distribution is similar in the case of perturbations. As can be observed, the errors produced by the analytical theory are greater in the interval , whereas in the rest of the intervals other large errors can be found.
Figure 18(a) compares the KBMJ12 versus KBMJ12_2 distance errors, whilst Figure 18(b) shows the box and whisker plot analysis for KBMJ12_2 distance errors in each subinterval. It can be seen that the modified analytical method KBMJ12_2 clearly improves the behaviour in the critical inclination subinterval. The same conclusion can be drawn from the comparison between KBMJ12 and KBMJ12_2 distance error statistical parameters shown in Table 8.
(a) KBMJ12 versus KBMJ12_2 distance errors not showing KBMJ12 outliers above 10000 m
(b) KBMJ12_2 distance errors in each subinterval
6. Conclusion and Future Work
A scalable closed-form analytical orbit propagator programme is being developed. In this paper, only zonal harmonic perturbation is considered. This perturbation can be handled by means of two analytical theories. Both are derived after removing the long-period perturbation terms due to the argument of the perigee. The first removes the short-period terms due to the mean anomaly, using the classical Delaunay normalization (DN), and then the transformed Hamiltonian is integrated by quadratures. It should be noted that Delaunay variables are used here. The second removes the short-period terms by means of the Krylov-Bogoliubov-Mitropolsky method (KBM). This process is formulated in polar-nodal variables. In order to validate both models a TLE space catalogue is used. Although this input is only valid for SGP4, it can be considered a reliable set of initial conditions for validating our propagator. These data are classified into three sets: , with , and any other data. Then, an exploratory data analysis is applied to analyze the errors produced by the propagation of the TLE catalogue through both analytical theories with respect to the numerical integration of the original problem. The considered perturbations are the and effects for two time spans of 7 and 30 days.
When the TLEs are propagated over days, the number of outliers is and in the DN and KBM cases, respectively, of which are common. DN outliers concentrate in two groups: low eccentricities, , and near critical inclinations, with . Similarly, KBM outliers are also concentrated in two groups: near critical inclinations, and low inclinations, , with medium-high eccentricities, , which is a group of outliers with an error below m. The critical inclination is handled using the approximation implemented in PPT2, HANDE, and AOP in the and cases. Preliminary numerical analysis seems to be very promising in the inclination interval . In general, our conclusion is that KBM is clearly more robust than DN in of the study cases over days, and whenever DN is better than KBM, it is only slightly better. The behaviour of the polar-nodal variables is excellent in all cases and does not present any problem, not even when .
On the other hand, when TLEs are propagated over days, the behaviour of the two methods is pretty good. The number of outliers detected in the DN and KBM cases for a distance error above 1000 m is and , respectively, of which are common. DN is worse than KBM only when ; however, for the rest of the cases DN shows better behaviour than KBM, in contrast to the time span of days. Both methods show worse behaviour near critical inclination, producing very high distance errors and two overflow cases.
In the near future Delaunay variables will be replaced with nonsingular variables in our propagator so as to avoid the analyzed singularities, and third-order analytical theories for perturbation will be generated. We will also analyze how using polynomial , (3), instead of may modify the original Hamiltonian and what its corresponding phase space may be. Tesseral influence will be the next perturbation to be included in our propagator.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.