1. Introduction
Nowadays, small bodies are interesting targets of in-suite explorations. Up to now, many missions have visited these objects. Determining the spacecraft’s orbit and the asteroid’s gravity field is one of the goals for these missions. To this end, the current technology heavily depends on observations from ground stations. There are several types of observation data often used in asteroid missions. In the earliest missions, Doppler and range-measurement are main types of data used in missions, such as the NEAR-Shoemaker and the DAWN mission [
1,
2,
3]. VLBI data is also widely used. For example, during the Chang’e-2 flyby of the asteroid Toutatis, its orbit was determined by the VLBI data together with the USB (Unified S-Band) data [
4]. Delta Differential One-way Ranging (DDOR) is a highly accurate measurement based on the VLBI technology [
5] and was used in the OSIRIS-Rex mission [
6,
7].
Most asteroids are irregularly shaped. Many studies have been carried out on asteroid shapes and gravity-field models. For example, the polyhedron model is one of the most widely used asteroid shape models [
8], but its associated form of gravity field is inconvenient to use in the gravity field recovery from observations. For elongated asteroids, a simple rotating mass dipole model is used to rapidly model the gravity field [
9]. Despite of these gravity models, the most commonly used one in the asteroid’s gravity field recovery is still the well-known spherical harmonics. By recovering the gravity harmonic coefficients, density and internal structures of an asteroid can be inferred [
10,
11,
12]. The gravity coefficients are usually measured simultaneously with the spacecraft’s orbit. For the first asteroid exploration mission, NEAR-Shoemaker, the gravity field up to the fourth order and degree was determined [
1]. Using Doppler and range data of the Dawn spacecraft, Park and Konopliv et al. obtained a gravity field up to 20th order and degree for Vesta and 18th order and degree for Ceres, the two most massive asteroids in the asteroid belt [
2,
3,
10,
13]. For asteroids of hundreds of meters in diameter or even smaller, it is usually difficult to obtain a high order and degree gravity field using traditional methods. One exception is the asteroid Bennu which is only about 450 m in diameter. In the OSIRIS-Rex mission, angle-measurement data between the spacecraft and particles around the asteroid made it possible to determine the gravity field to a high precision of 10th order and degree [
14,
15].
Experience from the OSIRIS-Rex mission indicates that accurate relative measurement around the asteroid may help improve the accuracy of the gravity-field recovery. We do not always have the chance to obtain relative angle measurement data as we did in the OSIRIS-Rex mission because not all targets have particles ejected from their surface, but we can use artificial relative-measurement data such as the inter-satellite range measurement, which is already extensively used in many applications around the Earth. For example, in the global navigation satellite system (GNSS), the addition of inter-satellite range data to the ground-observation data can significantly improve the accuracy of orbit determination [
16,
17,
18]. In addition, in the GRACE-FO mission, the precise inter-satellite range data are used to obtain the high-order gravity field of the Earth [
19].
An intuitive idea is to borrow this technology for asteroid missions. That is to say, we build accurate inter-satellite links between the spacecraft formation with one chief spacecraft and some cheap deputies which can move closer to the asteroid’s surface. For spacecraft formations in the proximity of asteroids, autonomous navigation using multiple sensors was already studied [
20]. For Earth satellites, it is impossible to autonomously determine the satellites’ orbit using only the relative inter-satellite range data, due to the well-known rank-deficiency problem. This is no longer the case for missions around irregular asteroids. The highly non-spherical gravity obviously disturbs the spacecraft’s orbit and efficiently removes the rank-deficiency problem, a fact which will be shown in our studies.
Taking an audacious step further, we attempt to study the feasibility of simultaneously determining the asteroid’s gravity field and the spacecraft’s orbit using only the inter-satellite range data. This is definitely impossible for Earth satellites, because determining the spacecraft’s orbit alone already requires the support from ground stations. But for asteroid missions it is possible due to the strong non-spherical gravity perturbation. To this end, we assume that the asteroid’s rotation is already well established, i.e, rotational parameters such as the pole angles and angular speed are already known. In the current study we assume that the asteroid is uniformly rotating along its largest principal axis in space. As a very preliminary setting, we use a gravity field up to fourth order and degree to simulate the orbits and the inter-satellite range data. However, we only try to determine the gravity field up to the second order and degree. By doing so, we simulate the force model errors in the orbit determination and gravity field recovery. The gravity field of the asteroid 433 is taken as an example [
1]. To simplify the notations, we call the problem of simultaneously determining the asteroid’s gravity field and the spacecraft the complete-inverse problem. We divide the complete problem into two sub-problems. Sub-problem 1 is the orbit determination with a prescribed gravity field, and sub-problem 2 is the gravity-field recovery with prescribed orbit elements.
In this paper, the feasibility of solving the complete problem using only the inter-satellite range data is analyzed. Our studies indicate that the complete-inverse problem is sensitive to errors in the initial estimation. It converges only when the initial estimation is close to the real value. When the initial estimation is poor, the iteration process usually fails. In our test, the accuracy of our best estimation is 3% for the term and 8% for the term. To help enhance the estimation accuracy, we propose two improved strategies. One is to extend the single link to multiple links by adding more deputies, and the other is to add the angle-measurement data between the chief and the deputy. Both strategies work fine in the sense that the convergence region is much wider. For the example case studied, the accuracy is 3% for the term and 0.15% for the term for the first improved strategy, and is 3% for the term and 0.1% for the term for the second improved strategy. At the end, we also discuss the influence of SRP and demonstrate the other perturbations only have little influence in our case.
This paper is organized as follows: The coordinate system, dynamic model and measurement model are introduced in
Section 2. Sub-problem 1 and sub-problem 2 are studied in
Section 3 and
Section 4, respectively. The complete inverse problem is studied in
Section 5. Two improved strategies and the influence of SRP are discussed in
Section 6.
Section 7 concludes this study.
2. Preliminaries
2.1. Coordinate Systems
The asteroid (433) Eros is a highly elongated asteroid, so it has a highly irregular gravity field. In our study, the gravity truncated at the fourth order and degree is chosen as an example. As a conceptual study, only two coordinate systems are involved—the asteroid-centered inertial frame and the asteroid’s body-fixed frame. Assuming that Eros is uniformly rotating in space, the transformation between the two coordinate systems is simply a rotation matrix along the asteroid’s rotational axis.
2.2. Force Model
For orbital motions around the asteroids, the two strongest perturbations are the non-spherical gravity and the SRP. For large-size asteroids such as Eros, the non-spherical gravity is much stronger than the SRP. For Eros, non-spherical gravity is several orders of magnitude larger than the SRP for low-altitude orbits within about ten kilometers above the asteroid’s surface. Therefore, it is appropriate for us to use the simplified force model.
The gravity field in the asteroid’s body-fixed frame can be described as
where
is the reference radius,
is the associated Legendre polynomials,
is latitude, and
is longitude.
and
are spherical harmonic coefficients. For the example case used in our study, the coefficients are adapted from Ref. [
1] and displayed in
Table 1.
Physical parameters of Eros are shown in
Table 2 [
1]. In our study, normalized units are used. The mass unit and the length unit take values in
Table 2. The time unit (T) is defined as
, and its value is also displayed in
Table 2. In this paper, all the units are normalized unless specified.
2.3. Measurement Model
Denote the positions of the two spacecrafts in the inertial coordinate system as
and
, and denote their relative position vector as
The inter-satellite range can be simply expressed as
In the discussion section, the inter-satellite angle measurement data is also used. It is a simple fact that
where
is the right ascension, and
is the declination.
To simulate the observation errors, we use the White Gaussian Noise as the measurement error. The standard deviation is set as 0.05 m for the range data and 5 arc-seconds for the angle data. One remark is made here. The two satellites are very close to each other. This fact leads to a negligible light time delay in the current study.
2.4. Orbit Types
Before we start studying the orbit determination problem, it is better to choose suitable orbits. There are two rules we should follow when we make our choice. First, the chief spacecraft should move on a practically stable orbit. Second, the orbits should be chosen such that the signal from the asteroid’s non-spherical gravity is strong.
In order to characterize the strength of the signal in the inter-satellite range data caused by the asteroid’s non-spherical gravity, we introduce the index
in Equation (
5). For the fourth order and degree non-spherical gravity field listed in
Table 1, we integrated two orbits and recorded the mutual distance
after an integration time. By adding a small error to the coefficients listed in
Table 1, we integrated two orbits starting from the same initial conditions and recorded the mutual distance
after the same integration time. The index
was, thus, defined as
. It is an obvious fact that the influence is stronger for a larger
value. Our studies find that the tesseral term has a more obvious influence on the mutual distance on short time scales. In addition, as we already mentioned in the introduction section, only the second order and degree gravity field will be recovered. Therefore, only the error to
is considered.
There are six orbital elements (
a for the semi-major axis,
e for orbit eccentricity,
i for orbit inclination,
for the longitude of the ascending node,
for the argument of perigee, and
M for the mean anomaly) to be considered for one spacecraft. To limit the dimension of the problem, we should introduce some constraints. The most important factor should be the semi-major axis which directly determines the spacecraft’s altitude. In order for the signal from the non-spherical gravity to be strong, the altitude should be sufficiently low. On the other hand, low-altitude orbits are generally unstable, which may lead to the spacecrafts’ collision. Based on this consideration and the highly elongated shape of Eros, the axes of both spacecrafts were chosen as
, which is a sufficiently low orbit and can avoid the spacecrafts’ collision with the asteroid within a moderate time. Due to the low-altitude orbit, we also set the initial orbit eccentricity as
. Moreover, considering the orbital stability of the chief spacecraft, the terminator orbit is a good choice [
21]. In our work, since the SRP is not considered, we approximated the terminator orbit as a near-circular polar orbit. As a result, we set the orbital element of the chief spacecraft as
.
For the deputy, we set
,
,
and
, and we changed the values of its longitude of ascending node
and orbit inclination
i to study the sensitivity of the inter-satellite range data to its orbital plane. Left frame of
Figure 1 displays the contour maps of measurement sensitivity when different values of orbit inclinations and longitude of ascending node are chosen.The abscissa is the orbit inclination of the deputy, and the ordinate is the difference in the longitude of the ascending node between the chief and the deputy. The sensitivity index is calculated using the following equation
where
n is the number of measurements.
is the inter-satellite range data when the gravity filed in
Table 1 was adopted, and
is the inter-satellite range data with a small error (+
) added to the
term. A measurement interval of 0.05(T) and an arc length of 100(T) are taken for
Figure 1, where (T) is the Unit time shown in
Table 2.
From the left frame of
Figure 1, we find that with the increase in the orbit inclinations, the measurement sensitivity
decreases. In addition, there is a periodic patten in the contour map when different values of the longitude of the ascending node are taken. It seems that the maximum sensitivity happens when the deputy’s
is close to the
value of the chief. Only judging from
Figure 1, we should choose prograde orbits (
) for the deputy. However, due to the instability of the prograde orbit, the deputy may quickly collide with the asteroid. In such a case, we are unable to obtain measurement data long enough. Therefore, in our study, we choose the orbit of the deputy as a polar orbit. This type of orbit is relatively stable, and receives enough influence from the asteroid’s non-spherical gravity. In this case, we set the deputy’s longitude of ascending node as
, shown as the red point in
Figure 1. In addition, the orbit of the chief and the deputy are shown in the left frame of
Figure 1.
5. The Complete-Inverse Problem
In
Section 3 and
Section 4.2, we showed the feasibility of the two sub-problems, respectively. However, the feasibility of the two sub-problems cannot guarantee the feasibility of the complete-inverse problem. In this case, the gravity field and orbit are coupled, and we need to estimate the orbit and gravity field together.
In our work, the mutual influence between the two sets of parameters (the gravity field and the orbits) is so strong that it is hard to converge. When the error of the initial estimate is too large, the complete-inverse problem is hard to converge. To solve this problem, we use the adaptive step size to control in the iteration process. As shown in Equation (
19).
is the calculated correction value by the iteration process. It is quite possible that
is too large. Directly correcting the estimate with this calculated value may lead to a worse estimate and eventually the failure of the iteration process. To overcome this problem, we deliberatively decrease the size of the correction value by introducing the step-size control parameter
k.
is the actual correction value taken and
k is an empirical constant which usually takes the value of 1. In such a case, when
is large (i.e., the estimate is far from the true value), the size of the actual correction
is always smaller than 1. When
is small (i.e., the estimate is close to the true value),
.
Test 3 shows two examples with good initial estimates. The simulation observation data of them are the same: The initial position is the same as
Section 3.2. But the initial estimate of initial position is different. Random error (with a threshold of 10 m) along all three directions is added differently in the two examples. The initial estimate of the gravity field is the same, where
and
. For an arc length of 50(T), the result is shown in
Figure 6 and
Table 3. And
Table 4 reflects the change of RMS for two examples before and after the calculation.
Obviously, the results after iteration are much better than the initial estimate. For test 3a, the residuals of observation are less than 0.5 m, and the position errors of estimated orbits are less than 0.8 m. For test 3b, the residuals are less than 0.3 m and the position errors are less than 2 m. In addition, the difference in the results of these two tests is small. The accuracy of the estimated gravity-field coefficients is about for test 3a and about for test 3b. The RMS of them are similar. Thus, we can say that the orbit determination process has converged.
However, when the initial estimate is poor, the iteration process converges harder, and the accuracy of the results obviously decreases. In test 4, we randomly added a position error of 50 m in all the three directions of both the spacecraft. The initial estimate of the gravity-field coefficients are
and
. Setting the arc length as 50(T), the result is shown in
Figure 7 and
Table 5.
Figure 7 and
Table 5 indicate that the iteration process has not converged, because some of the final result is even worse than the initial estimates.
Studies in this section demonstrate the feasibility of the complete-inverse problem. However, a good initial guess is necessary for the convergence of the problem. In order to obtain a good initial estimate, we tried to use the particle swarm optimization method (PSO) to find the global optimal solution. However, most of the time, the PSO fails to converge or converges to a local minimum far away from the true value. Thus, this attempt failed. In the following, we propose two improved algorithms to help solve the narrow convergence-region problem.
7. Conclusions
In this work, we mainly demonstrated the feasibility of orbit determination and gravity-field recovery using only inter-satellite range data in three steps. First, sub-problem 1, in which the gravity field is known a priori and the orbits of the chief and the deputy are autonomously determined, is studied. Second, sub-problem 2, in which the orbits of the chief and the deputy are known a priori and the gravity field is recovered, is studied. Then, the complete-inverse problem in which both the orbits of the chief and the deputy and the gravity field are determined is studied. Our study shows that the complete-inverse problem is invertible, but the convergence region is very narrow, i.e., an initial estimate close to the real one is necessary. Two improved algorithms are proposed to solve this problem: using multiple links by adding more deputies and adding the angle measurement data between the chief and deputies. By our calculation, both algorithms can increase the convergence region and improve the accuracy. Our studies find that the accuracy of the zonal term is worse than the accuracy of the tesseral term. In addition, we also analyzed the influence of different orbit types on the inter-satellite range data and the influence of different arc lengths on the orbit determination.
We need to stress that this work is only a feasibility study. Our conclusions are based on some assumptions. For asteroids of hundreds of meters in size or smaller, it is possible that SRP is one order of magnitude larger than non-spherical gravity which may seriously challenge the assumptions in the current study. The study indicates that it is possible to solve the complete inverse problem using only the inter-satellite range data, but a good initial estimation is necessary and the accuracy of the results is limited. A more realistic mission scenario is to combine the accurate inter-satellite range data along with ground observation data to do the complete inverse problem. We are currently working on this project.