Abstract
Autonomous cislunar spacecraft navigation is critical to mission success as communication to ground stations and access to global positioning system (GPS) signals could be lost. However, if the satellite has a camera of sufficient quality, geometric line-of-sight (unit vector) measurements can be made to known lunar landmarks (e.g., Tycho Crater) to provide observations that enable autonomous estimation of the position and velocity of the spacecraft. In this study, an improved batch gaussian initial orbit determination (IOD) differential corrector (DC) algorithm, based on the approximated values of the two-body f and g series, is applied to initialize a (non-conic based) circular restricted three body problem (CR3BP) extended Kalman Filter (EKF) navigator. This navigator collects geometric line-of-sight unit vector (angle only) measurements to a known location on the Moon to sequentially estimate the position and velocity of an observer spacecraft flying on an approximate southern L1 Halo orbit.
In this study, it was found that the best approach is to initialize the CR3BP EKF (navigator) using the solution from the batch DC filter with at least 10 measurements taken against the perceived centroid of Tycho Crater. Thereafter, it is best to continue the navigator with subsequent measurements taken against the same center coordinates of the Tycho Crater, where these coordinates are now expressed in the CR3BP rotating frame. For successful conic-based batch filter initialization and long-term CR3BP EKF convergence, it was found that the cadence for all optical measurements should be taken at 10 minutes for a simulated measurement noise of 0.1° one sigma uncertainty about the line-of-sight measurement unit vector.
1 INTRODUCTION
Recent years have seen a renewed interest in the development of cislunar space. As efforts are made to increase the number of satellites in both near-Earth and cislunar space, ground-based navigation assets are becoming increasingly taxed (Bradley et al., 2020). This necessitates the development of autonomous navigation strategies that will allow satellites in these regions to operate effectively when ground-based supporting assets are unavailable. One possible framework for such autonomous navigation is to use optical measurements to known objects in the Earth-Moon system, such as satellites, asteroids flying near Earth, and the Moon. Bradley et al. (2020) consider the feasibility of this navigation strategy and conclude that artificial satellites, the Moon, and lunar landmarks provide effective targets for navigation in the cislunar region. Their study addresses practical navigation concerns, such as visibility and geometry constraints, illumination conditions, and location uncertainty of potential targets. The study also demonstrates the performance of an optical AutoNav framework for trajectories used in the Lunar Reconnaissance Orbiter (LRO) and Gravity Recovery and Interior Laboratory (GRAIL) missions. Their treatment of the problem is quite general and thorough.
The purpose of this study is to investigate the performance of a custom-built angles-only extended Kalman filter (EKF) navigator that estimates the state (position and velocity) of a spacecraft flying in a halo orbit of the Earth-Moon L1 libration point. An approximate Southern L1 Halo orbit, based on computed initial conditions from Zimovan (2017), is chosen as the baseline orbit for this study. We choose to simulate the “lost in space” problem by initializing the states with solutions from a two-body f and g series batch differential correlator (DC) filter, which requires no a priori information to perform the initial orbit determination (IOD).
The reference L1 Halo orbit, along with a simulation and navigation result, are visualized in an Earth-Moon rotating frame in Figure 1 and a notional Earth-centered inertial (ECI) frame in Figure 2.
(Approx.) Southern L1 Halo orbit CR3BP frame
(Approx.) Southern L1 Halo Orbit ECI Frame
This study aims to identify promising targets for use by the navigator and offer some explanation of why such targets may be advantageous. We aim to identify which navigational targets are preferable for EKF initialization via the batch DC filter and which are best for long-term EKF state estimates. A recommendation is also given regarding the required frequency of measurements.
Among all target cases, it was found that the best candidate DC batch filter solution to use in initializing the circular restricted three body problem (CR3BP) EKF was based on the perceived center coordinates of Tycho Crater. Once initialized, the best EKF navigation performance was obtained by continuing to use Tycho Crater as the observation target (as opposed to switching to a different target).
In all cases, it is recommended that at least 10 measurements be provided to the DC batch filter to ensure a reasonable initialization of the EKF navigator state to prevent navigator divergence. It was found that the best cadence for all measurements should be 10 minutes for a simulated measurement noise of 0.1° one σ in the observed right ascension and declination of the navigation target.
The noise value of 0.1° one σ was chosen for the following reason: the accuracy performance for a novel and recent low-cost multi camera star tracker system for small satellites has been determined by Zhao (2020) to be significantly better than its stated accuracy objective of 0.1° one sigma error. For existing miniature star trackers, such as those used for the AeroCube-OCSD CubeSat, their pointing accuracy is 0.02° one sigma error, (Janson et al., 2016). Even though this uncertainty is smaller than that of the newer star trackers, we chose to be conservative and apply the larger, more challenging of the two values in this study.
The estimated states in the DC batch filter are the position and velocity in a notional inertial frame. All geometric line-of-sight unit vectors are equally weighted in the batch filter, and there is no a priori information or a priori uncertainty. The only assumption made by the batch filter is two-body motion about the Earth. Indeed, this may be a poor assumption for the CR3BP dynamics of a Halo orbit, but it could be the only assumption available for the lost in space scenario. The estimated states of the CR3BP EKF are position and velocity expressed in the CR3BP rotating frame. The EKF is initialized with the state and covariance solution from the DC batch filter rotated into the CR3BP frame.
Figure 3 illustrates an example of the observation model, which is that of the geometric line-of-sight unit vectors pointing from the observing spacecraft to a target asteroid (example target body), shown by the green arrows. Note that each geometric line-of-sight unit vector can be represented by a right ascension and declination from the point of view of the observer. Within Figure 3, the vectors r (red) represent the unknown coordinates of the spacecraft and can be expressed in either the inertial or rotating coordinate frame. The vectors R (in black) are the known position vectors of the navigation target object, which are known from ephemerides (e.g., Jet Propulsion Laboratory (JPL) DE441 (Park et al., 2021)). Finally, the vectors ρ (in blue) are the unknown relative range vectors between the asteroid and the spacecraft. The relative ranges are not known and are not directly measured in this observation model. By multiplying the green geometric line-of-sight unit vectors by minus one, we obtain the corresponding unit vectors associated with the unknown blue ρ relative range vectors. It is only the magnitudes of these relative blue range vectors that are unknown and are presented in Figure 3 for clarification. Within the batch IOD algorithm, these unknown relative (scalar) ranges (corresponding to the measured geometric line-of-sight unit vectors) are initialized to zero and are used to estimate the initial state of the spacecraft.
Observation model for both the batch and EKF filters.
Note that the direction of the vectors in Figure 3 is very important. If an arrow points from the spacecraft to the asteroid, this means that an actual geometric line-of-sight unit vector measurement is being physically measured from the spacecraft. In the computation of the spacecraft position, this unit vector is flipped (multiplied by minus one) and, thus, points from the asteroid to the spacecraft. In other words, we are computing the position of the spacecraft with respect to the asteroid, i.e., we form the measurement as if the observer is virtually located on the asteroid, not the spacecraft. If we knew the position of the spacecraft and were interested in computing the position and velocity of some unknown asteroid, then the system of arrows would be reversed, as seen in the study by Hinga (2018).
Within this study, several target bodies for navigation were considered, including a geostationary (GEO) and low Earth orbiting (LEO) satellite, the coordinates of the Moon’s center, and the center of Tycho Crater (−43.31° latitude and −11.36° longitude). In all cases, it was assumed that each target body was visible at all times.
1.1 Initial Orbit Determination - Inertial Reference Frame (Two-Body Problem - Phase One)
The problem of determining the orbit of an unknown object began with the advent of celestial mechanics seen in the works of Laplace (1780) and Gauss (1857) concerning the motion of heavenly bodies about the Sun. Their angle-only techniques utilized three observations to compute a position of a celestial object without the knowledge of the range, which had to be guessed with the help of the roots of an eighth-order polynomial. In 1889, Gibbs (Bate et al., 1971) developed his own technique enhancing the Gauss method of position estimation to include the determination of the velocity. Knowing both position and velocity, thus, defines an orbit in space. Herrick improved on Gibbs’ technique (for short arcs) with the use of a Taylor series to compute the velocity at the middle position vector. Clearly, these techniques were developed for celestial applications well before the beginning of the space age and availability of the computer Herrick (1971). However, astronomer Paul Herget (1965) introduced an algorithm that uses more than three angle-only measurements to estimate an orbit in which an iterative approach is applied through the variation of guessed geocentric distances to minimize a set of residuals in a least squares approach, using as many observations as are available, performed on an IBM 1620 computer.
Over the last several decades, many iterative methods to estimate the orbit of an unknown object (natural or artificial) using angle-only measurements have been developed. The Double r-iteration technique by Escobal (1965) iterates on an initial guess of the range between the observer and a target object via the numerical partial derivatives and a Newton–Raphson iteration to converge on the true range. The Gooding (1996) method, using a minimum of three measurements, requires an initial good guess of the first and third ranges and whether the orbit is pro or retrograde. Common to all of these methods are the assumptions they make about a target satellite to converge to a solution for its orbit.
However, instead of trying to guess “good” starting distances for scalar ranges, ρk, initializing all of them to zero provides for excellent starting values and have always converged to the nontrivial solution (in this study) for orbiting platform observations. A similar result was also found to be true by Karimi & Mortari (2011). A more detailed description of this iterative algorithm can be found in the study by Hinga (2018) and is also discussed in Sections 5.1 and 5.2.
1.2 Navigator - An EKF based in the Circular Restricted Three-Body Problem Reference Frame
Solving the IOD problem in the CR3BP introduces the issue of trying to determine an orbit in a rotating frame where conic sections no longer govern the motion of the spacecraft. The familiar f and g series are not applicable in the CR3BP because they are formulated based on the assumptions of two-body conic motion. In this study, an initial guess for the state of an unknown non-Keplerian orbit is first computed using a two-body approximated f and g series Gaussian-batch IOD. This estimated state is expressed in a notional inertial frame. It is then transformed into the rotating CR3BP frame. If the time span over which the measurements are taken is reasonable, then the IOD solution that uses the two-body motion can be sufficiently close to truth to prevent CR3BP EKF divergence. In this study, it was found that upwards of 90 minutes of measurements taken at a cadence of no more than 10 minutes apart can provide a good enough initial state estimate to enable convergence of the CR3BP EKF.
We find that if the spacecraft is in a region of the orbit that is fairly conic with respect to the Earth (i.e., near apolune), the two-body f and g batch filter will generally succeed in providing a reasonable IOD estimate. It is beyond the scope of this investigation to examine this issue further, although we might speculate that this leverages some of the underlying dynamics in this particular halo orbit. Once the CR3BP EKF orbit estimator is initialized, it continues to track the state using subsequent geometric line-of-sight unit vectors to the observed target body in the CR3BP reference frame.
CR3BP EKF navigation error is defined as the difference between the estimated positions and velocities from the known true states of the spacecraft within the simulation. The Jacobi Constant (the only known integral of motion in the CR3BP) is computed in the simulation to assess the quality of numerical integration. For every 50-day simulation and navigation solution, the difference between the very first and last Jacobi constant computed (fortunately) remained on the order of 10−15. This constancy of the first integral of motion satisfies only the necessary condition on the accuracy of the numerical integration but does suggest that numerical propagation is not obviously wrong.
Finally, it is necessary to mention that all satellite navigation and measurement simulations, algorithm software development, and verification were carried out inside the custom-built and validated Hinga SpaceNavSim (Hinga, 2023).
1.3 Observation Model G and the H Matrix
Within the EKF, the observation model is denoted 
. Intuitively, the observation model is the relative direction of the spacecraft (T) with respect to the observing platform (P), expressed as a unit vector. Recall that in our formulation of the problem, the observing platform (P) is the object we are using to navigate autonomously (e.g., Tycho Crater), while the satellite (T) is the object that is attempting to navigate. Letting 
denote the position of the spacecraft (T) and 
the position of the observing platform (P), the observation model is given by Equation (1).
1
Expressing 
and 
in terms of their components in some orthonormal coordinate system (rotating or inertial), 
we have:
2
3
We can then write:
4
The components lx, ly, and lz can be expressed in terms of a right ascension (α) and declination (δ) in the chosen coordinate frame, such that:
5
To relate the elements of the state vector to the target satellite X(t)satellite, we take the partial of G with respect to the elements of the state vector:
6
7
8
9
10
11
Now, we can define the H matrix and relate the change in a measurement to a change in the state vector.
12
Notice that the last three columns are all zero. This shows that velocity is not directly sensed by the geometric line-of-sight unit vector and there is no relation between the change in the elements of G, lx, ly, and lz to the change in velocity.
13
2 EXPERIMENT
In preparation for a potential flight experiment, simulations are carried out in custom-built astrodynamics software, and performance of the navigator is assessed. The formulation of the CR3BP EKF navigator is described in any standard textbook covering the EKF (see algorithm in the Appendix Section 5.4 for convenience).
Throughout our experiments, the CR3BP EKF navigator is intentionally given complete knowledge of the system dynamics to assess its ability to correct for simulated errors in measurement. It implements a high-quality variable step size numerical integrator found in Numerical Recipes in C, 2nd Edition (Vetterling et al., 1992). The fixed step Runge Kutta (RK4) is used by the simulation to propagate the true states.
To incorporate realistic measurement error, we must first clarify what the navigator is actually measuring. An angles-only optical navigator makes use of measurements of right ascension (α) and declination (δ) angles to compute geometric line-of-sight unit vectors (to the perceived center of the target object) in a chosen coordinate system (for our purposes, a CR3BP rotating coordinate system and an inertial). Due to inaccuracies in ephemeris data, knowledge of spacecraft orientation, and a variety of other error sources, measured values of α and δ are never exactly correct. In the effort to incorporate all real-world errors into this simulation, we perturb the true values of α and δ by small values (dα dδ) to form a corrupted measurement. These values of dα and dδ are sampled from a normal distribution about 0° with a standard deviation σ of 0.1° for reasons found by Zhao (2020) above. It is true that systematic and persistent errors from ephemeris of target objects or solar radiation pressure (and other unmodeled effects) cannot be adequately represented by random noise. However, they are assumed to be insignificant so that with appropriate process noise in the navigator, the state solution will not diverge but can capture dominant effects, thus, enabling coherent analyses. Moreover, one can increase the fidelity of the dynamical model to include known persistent effects, such as solar radiation pressure. Indeed, a higher fidelity ephemeris navigator would be necessary for mission applications.
As previously mentioned, our experiments aim to assess the performance of an optical navigator for a spacecraft flying in an L1 southern halo orbit. In our investigation, we aim to determine which objects the navigator might observe that offer the best navigation results. We consider an asteroid flying near Earth (EC2014), the perceived center of the Moon, the perceived center of the Tycho Crater, a LEO and a GEO satellite. (For these lunar targets, the perceived centers could be derived from a centroid of brightness or centroid of geometry.) For simplicity, we do not consider lighting conditions in this study. Neither do we consider geometric line-of-sight unit vector observability (or opportunity) to LEO and GEO satellites, with the understanding that there are sufficiently many of these satellites to allow a target switch when one passes out of view. Error introduced by switching between targets is not simulated or addressed in this investigation.
Simulations are carried out for a 50-day period beginning on March 5, 2014. This date was chosen for the close approach of the asteroid EC 2014 during this time. This nominal trajectory is illustrated in Figure 2 and was generated by propagating the initial conditions (provided by JPL’s Horizon Webpage) forward in time and transforming them to the CR3BP frame. Finally, we consider a nominal time interval between noise-free measurements (dt obs) of 60 seconds to examine which of the five target bodies provides the best initial IOD and long-term navigation. All other target body motions are simulated in the same manner. After the best target body is determined, a variation of noisy measurement cadences for the winning target body is investigated to examine how the state estimate is affected by increasing time spans between observations. The best target candidate for the two-body inertial batch filter solution is also determined.
2.1 Phase One - Spacecraft Inertial f and g Series Batch Filter IOD - Results
The f and g series batch filter computes a snapshot in time of the lost spacecraft’s position and velocity (see Sections 5.1 and 5.2). Although it provides an imperfect initial guess of the state in CR3BP, it is good enough so that the CR3BP EKF eventually converges close to truth for all cases with no measurement noise, as shown in Table 1. For the cases with measurement noise, as shown in Table 2, the successful scenarios are when the target bodies are the (perceived) Moon center or the (perceived) center coordinates of Tycho Crater. For nonlunar targets, our experiments suggest that the relative motion between the lost spacecraft and the target body is small enough such that the noise (0.1° one sigma) in the angles of the geometric line-of-sight unit vectors becomes significant enough to degrade state initialization to the point of causing eventual divergence in the Kalman filter. It is interesting to note that when the 10- or 30-minute measurement cadence is selected, only measurements made against the LEO satellite, Moon center, and the Tycho Crater are good enough to allow CR3BP EKF convergence. Although the batch filter (by itself) still converges, any solutions for measurement cadences above 60 minutes will cause divergences in all CR3BP EKF instances. For all measurement cadences performed, the batch filter solution uncertainties in position are on the order of 300–1000 km, and those for velocity are about 6 km/sec. Obviously, compared with the actual errors illustrated in Tables 1 and 2, the batch filter is overconfident in position and underconfident in velocity. Nevertheless, convergence of the CR3BP EKF (during noise conditions) using the Moon targets was obtained.
No Measurement Noise (1 min Interval): Inertial Batch Solution Difference to Truth for Lost Spacecraft
With Measurement Noise 0.1° one Sigma (1-min Interval): Inertial Batch Soln Difference to Truth for Lost Spacecraft
2.2 Phase Two - Spacecraft CR3BP EKF - Results
Performance of the navigator is assessed on the basis of the accurate recovery of both position and velocity. It is assumed that the target body trajectories in this study, e.g., an asteroid or LEO spacecraft, would be very well known at any time, as those of interest are available in the JPL binary ephemeris file DE441, (Park et al., 2021). This assumes that the ephemeris error of a target body would be small if the time on board the observing spacecraft is well known. Figure 4 illustrates errors in position using each navigation object (target body). It can be seen that all navigation examples require about 5 days to converge. This is because the navigator warms-up (covariance becomes confident) using the perceived center of a particular target to compute the geometric line-of-sight unit vector measurements.
Position error - all targets (no measurement noise)
From Figure 5, we see that the Earth satellite target navigators struggle for about 12 days until the actual error settles down to a maximum of about 15 km and 80 m/sec in position and velocity, respectively. The Moon target navigators require about 4 days of simulation time (the warm-up period), until their state covariance remains inside a consistent bound of values. For the Moon targets, this time until convergence is directly related to the time it takes the spacecraft to arrive at lunar periapsis (more on this later). Over more time, however, as the navigators continue, their performance becomes more apparent, thus allowing for a good comparison of the candidate targets. We note that navigation using the asteroid target body is significantly less accurate after 10 days than other observation candidates as the simulation progresses. After 15 days, the asteroid navigator begins to diverge. As the asteroid reaches a great enough distance from the observer, the relative change in the angle between subsequent geometric line-of-sight unit vectors becomes very small.
Position error - all targets except Asteroid (no measurement noise)
This causes each new measurement to be essentially the same (or contain little or no new information) as the previous measurement because the relative change in dynamic motion from the point of view of the observer is insignificant. To make matters worse, for the asteroid target, the noise present in the measurement becomes relatively significant compared with the low information content in the observation and ultimately contributes to filter divergence. See Figure 4.
Noting the significantly worse performance of the asteroid in the recovering position, subsequent figures demonstrate the other possible navigation targets. Figure 5 illustrates position error results for a subset of possible targets and provides more insight into navigator performance in recovering position. It can be seen that it requires about 35 days for the LEO and GEO navigators to converge to position errors that are already obtained by the Moon targets after 5 days. Also in Figure 5, the results of the two Moon targets are plotted. They demonstrate better navigator performance in recovering position, with the Tycho Crater performing slightly better than the perceived Moon center.
In comparing the navigator recovery of velocity, we consider both velocity magnitude, Figure 6, and direction, Figure 7. Due to the relatively poor performance of navigation based upon an asteroid, we consider only the Earth-orbiting artificial satellites and the Moon coordinates. In all cases, we see several spikes in velocity magnitude and direction, at times very close to the arrival time at lunar periapsis, which is about 12, 20, 28, and 36 days. For the Moon targets, their performances reach lower errors much sooner and have smaller maximums throughout the simulation. This is because the measurement information content (contained in the H matrices of the Moon targets) is much larger in magnitude. (See Section 1.3 for the definition of the H matrix and the model of observation G, both of which are standard elements of an EKF algorithm described in any standard textbook.) This allows the corresponding Kalman gains to better reject errors in the measurements and reduce the size of the spikes relative to those seen in other target bodies of navigation.
Velocity magnitude error (no measurement noise)
Velocity direction error (no measurement noise)
Having identified the Tycho Crater as the winner of the candidate target bodies, we next consider the effect of varying the time between subsequent observations while adding noise to its measurements. Figure 8 and Figure 9 illustrate the performance of position error for the measurement time intervals of 1, 10, 30, and 60 minutes, respectively. We can see that the 60-minute measurement time step causes a very large error compared with the smaller intervals. This can be interpreted to mean that too much linearization error is accumulating between Kalman filter measurement updates at this cadence. Because the uncertainty of the position solution is a factor of 100 times smaller than the actual error, this filter solution is unacceptable. See Figure 10. Discarding this 60-minute case and looking at Figures 9 and 11, we see the 1-, 10-, and 30-minute intervals for the error and the uncertainty of the position solution, respectively. As expected, the 10-minute cadence performs better than the 30-minute cadence. While the confidence converges after 5 days, it is not until 20 days that the uncertainty in the 10- and 30-minute solution reflects an appropriate magnitude of confidence. That is, the error is bounded by the uncertainty. It may be surprising that the 1-minute solution does not always perform better than the other two, in both error and uncertainty. This is because more measurement noise is getting through to the state updates and is not being filtered out successfully until after day 35. It may be that the measurement angles subtended during the 1-minute cadences at large distances from the Moon are not large enough compared with the noisy angle of 0.1°, one sigma. However, during close approaches to the Moon (perilune), the change in angles of measurement (δ and α) are more significant between each update and significantly larger than the assumed noise value. This allows for the overall error and uncertainty to go down over time. It is encouraging that these remaining covariance cases do indeed bound the actual error before and after the perilune approaches at every 8 days after the 5th day and thereafter for the remainder of the simulation. While the 1-minute measurement cadence covariance case of uncertainty is still acceptable, it does not perform as well as the 10-and 30-minute cases.
Position Error - Tycho Crater - All Measurement Intervals - (with measurement noise)
Position Error - Tycho Crater – 1-, 10-, 30-Minute Measurement Intervals (with measurement noise)
Position uncertainty - Tycho Crater - all measurement intervals - (with measurement noise)
Position uncertainty - Tycho Crater – 1-, 10-, 30-min measurement intervals (with measurement noise)
Figures 12, 13, 14 illustrate the corresponding performances in velocity magnitude recovery. It is not surprising to see that the 60-minute case shows poor performance in error and uncertainty, and by discarding it, we expose something more encouraging. Filter velocity error performance is again quite good at 20 days (for all remaining candences) and beyond, while the uncertainties in velocity, Figure 15, reach a conclusion at 5 days and merely oscillate around 10, 5, and 4 m/sec for the remainder of the simulation (1, 10, and 30 min, respectively), with a slight downward secular trend. It can be seen that the 10-min case performs better than the 30-min case. The explanation as to why the 1-minute measurement step (for velocity) underperforms compared with the other step size cases is the same as that for position recovery. That is, the magnitude of the noise in the measurement errors is more significant than the information content (in the component angles) of the geometric line-of-sight unit vectors.
Velocity magnitude error - Tycho Crater - all measurement intervals (with measurement noise)
Velocity magnitude error - Tycho Crater – 1-, 10-, 30-min measurement intervals (with measurement noise)
Velocity magnitude uncertainty - Tycho Crater - all measurement intervals (with measurement noise)
Velocity magnitude uncertainty - Tycho Crater – 1-, 10-, 30-min measurement intervals (with measurement noise)
Figures 16 and 17 illustrate the corresponding performances in velocity direction recovery. After we discard the 60-minute case, we observe similar behavior in the occurrences of peak amplitudes in the estimated velocity direction error versus time, just before perilune passages every 8 days starting at the 12th day. As we might intuitively expect, we see that navigator performance generally improves as the time between subsequent observations is reduced, but only down to a lower limit.
Velocity direction error - Tycho Crater - all measurement intervals (with measurement noise)
Velocity direction error - Tycho Crater – 1-, 10-, 30-min measurement intervals (with measurement noise)
For measurement cadences below 1 minute, filter performance degrades as the measurement noise in the measurement angle becomes larger than or at least more significant in size than the relative change in the measurement angles themselves. It was noticed that the filter struggles to discern between the actual change in the measurement angles and the noise. Therefore, a so-called tuning effort (in-process noise (Q) and measurement noise (R)) was focused over this measurement time frame so that the accumulation of linearization error in all time frames could be mitigated and the noise filtered. This allowed the filter under noisy measurement conditions to achieve similar residual errors, after 12 days, as those converged to at 5 days by the same filter with zero noise measurements. See Figures 5 and 9. At the other end of the scale of measurement cadences, the 120-minute measurement interval step, or greater, caused eventual filter divergence, no matter how the filter was “tuned”. The tuning process of this study is described as follows.
Under noisy measurement conditions, it was possible to improve the navigator performance by tuning or adjusting the (constant) diagonal elements of the process noise matrix Q (that compensates for mismodeled dynamics) of the standard EKF, which appears in the predict phase of the Kalman filter (see Section 5.4). In this study, the process noise was approximated as very small constant values on the diagonals of the Q matrix. The best constant diagonal elements of Q were obtained by varying them through a range of values over multiple performances of the filter, until a local minimum in the residual errors (with respect to the simulation truth) was found. Those elements discovered to give an optimum difference to truth were Qxx = Qyy = Qzz = 6.0908 × 10−22 and 
, where the values are nondimensional and unitless in position and velocity. All nondiagonal elements of Q are zero. See Table 3. The best values of the diagonals in the R matrix were found to be 1 × 10−10 in units of squared radians.
A Monte Carlo approach could have been employed to continue the search for still better values based on perturbations of the initial conditions of the state. However, the belief is that the tuning performed in this CR3BP study was thorough enough to compensate for noise at all measurement cadences and for the error due to mismodeling of the dynamics. Indeed, performing a Monte Carlo in a future study that involves a higher fidelity simulation and a naive navigator would be worthwhile.
3 THE FROBENIUS NORM
To understand why the Moon might provide better navigation results than other objects considered, we compute the Frobenius norm of the navigator’s sensitivity matrix, H, as a function of time. The role and definition of the H matrix in the formulation of an angles-only EKF navigator are discussed in Section 1.3. Roughly speaking, the sensitivity matrix relates a change (or error) in state to a change in observation, G. The Frobenius norm gives an idea of the size of the elements in a matrix and of the magnitude of the matrix itself, represented as a single scalar value.
From Figures 19 and 20 we clearly observe that the Frobenius norm (size) of H is significantly larger for Moon targets than any other object for the vast majority of the simulation. This implies that a small change in state of the spacecraft produces a relatively large change in the associated observation when compared with navigation targets closer to Earth (but very far away from the spacecraft).
Frobenius norm of H vs. time - all targets (no measurement noise)
Frobenius norm of H vs. time - ignoring Moon targets (no measurement noise)
Geometrically, this is intuitive to understand; as the distance between the spacecraft and navigation target is reduced, a small change in the position of the spacecraft can result in a large change in the geometric line-of-sight direction (the measurement). As such, for a halo orbit that remains relatively closer to the Moon than to Earth, we would expect that a small change in the spacecraft’s position would result in a larger observational change for a lunar navigation target than for a navigation target closer to Earth. This is reflected in the magnitudes of the elements of the corresponding H matrix, or similarly, the Frobenius norm thereof.
The Frobenius norm of the H matrix becomes a maximum precisely at the moment in time of lunar periapsis (perilune) when the spacecraft arrives at its closest approach to the Moon. A close examination of Figure 19 revealed the perilune passage times to be: 4.03327, 12.1002, 20.1667, 28.2336, 36.3001, and 44.3667 days. Figure 18 illustrates the plot of spacecraft position dotted with velocity versus time. The times at which these values pass through zero (at downward crossings) were found to agree exactly with the times corresponding to the maximum values of the Frobenius norm of H.
Perilune at Downward Crossings (at zero - green stars) of r dot v
The verification of the perilune passage times is important because they serve as reference times at which the behavior (spikes) of the Kalman gain and state error can be understood in regard to the times of the H matrix peak values. This allows us to further understand the behavior/performance of the navigator in its recovery of position and velocity. This will be discussed in the next section.
4 CONCLUSION AND FUTURE WORK
The results of this study indicate that optical angles-only navigation for a satellite flying in this particular L1 halo orbit can be accomplished by initializing the state of a CR3BP EKF with a solution (both state and covariance) from an f and g series batch filter using 10 measurements, taken at a cadence of 10 minutes, against the perceived center of the Tycho Crater. The same target (Tycho) should then be used for further sequential observations at the same 10-minute measurement cadence. This measurement cadence is meant to be applied during all times in the orbit regardless of the observability of the Tycho Crater. Cadences larger than 1 hour are not desirable as they cause the filter to fall behind the fast-changing geometric conditions, when the spacecraft repeatedly approaches and departs through perilune with every orbit, see Figure 21. Notice the 60-minute spike in the induced norm of the Kalman gain leading or lagging behind the other measurement cadences by approximately 2 or 3 days. This leads to eventual divergence. The 1-minute measurement cadence does provide an adequate state estimate, but it has larger uncertainties and residuals.
Induced norm of Kalman gain vs. time - Tycho Crater (with measurement noise - all cadences)
Induced norm of Kalman gain vs. time - Tycho Crater (with measurement noise – 1-, 10-, 30-min cadences)
The results of the approach considered in this investigation suggest that a lost in space satellite may be able to autonomously navigate without a priori knowledge of its state. In other words, a conic IOD method can sometimes provide a good-enough estimate to initialize a nonconic, sequential orbit estimator.
These results also suggest that a greater observation frequency tends to increase navigator performance, but only to a lower limit. As the navigation step size is reduced beyond this point (below 1 minute), errors in measurement are as, or more, significant than the actual observed changes in the angle. This leads to degraded navigator performance.
A word about the 4-day warm-up period that the Moon target navigator requires to become confident in its estimate is warranted. An additional simulation was conducted in which this spacecraft navigator was initialized at 1 day prior to perilune. It was discovered that the warm-up period was precisely 1 day, the time required to reach perilune. The time until the actual error was bounded by the converged uncertainties in state was (again) equal to the amount of time to complete one full revolution of the orbit, i.e., about 8 days. All behaviors in error, Kalman gain, H matrix magnitudes, etc., were repeated as seen in the original simulation, but merely offset by 3 days. This suggests that the navigator’s performance is highly dependent on the perilune approach.
Navigation results for LEO and GEO satellites and the asteroid suggest that these target bodies are not useful. Although the LEO satellite can provide for a good batch IOD initialization of the CR3BP EKF, an eventual switch to the Tycho Crater would be necessary to ensure stability in the presence of noise or larger measurement cadences.
Considering the times of perilune passage, we observe that the spikes in both the position and velocity estimate errors occur around 110 minutes prior to perilune arrival. By the time perilune is reached, these errors have been reduced (via the Kalman gain in successive state estimates). This can be observed through careful inspection of Figures 9, 13, and 17. These spikes illustrate the effect on navigator performance of both rapid dynamical changes near perilune as well as the effects of range and relative angle changes on the values of the partials of the sensitivity matrix H.
Larger errors in position and velocity estimates near perilune are not surprising from a strictly dynamical perspective. Approaching perilune, velocity and acceleration are large, and the spacecraft state evolves very rapidly. This implies that even a small initial error in the navigation state with respect to the true spacecraft state can be magnified as the spacecraft passes through perilune. The maximization of the H matrix also occurs, as mentioned previously, precisely at perilune. The rapid growth of the H matrix as perilune is approached, coupled with an increased covariance prediction (due to the sensitive dynamics as represented by the state transition matrix Φ), results in a spike in Kalman gain (K) slightly before perilune.
However, the Kalman gain is observed to decrease rapidly to very small values by the time perilune is reached. The covariance in state also behaves similarly, as the uncertainties in position and velocity collapse to smaller values just prior to perilune arrival. This illustrates the effectiveness of the update phase of the navigator, in which the Kalman gain is able to apply significant corrective updates to our posterior state and state covariance estimates. With successive passes at perilune, this behavior and performance of the Tycho Crater navigator continue with success for the entire duration of the 50-day simulation.
Further investigation is needed to determine if this overall assessment/understanding of the CR3BP EKF navigator is generally reasonable and can be extended to a higher fidelity simulation and navigator and in other orbits. It is expected that, because this Batch-CR3BP-EKF algorithm (and winning candidate Tycho Crater coordinates for navigation) captures the majority of the dynamics (the Earth–Moon perturbations) of the real cislunar environment, its elevation into an ephemeris-based navigator would be successful.
An investigation should be carried out that examines where in this halo orbit it is best to apply this two-body f and g series batch filter (as an initializer) and whether the results are improved if the batch IOD uses the Moon (as opposed to Earth) as the central body. A subsequent investigation could also consider the performance of this approach in a variety of other cislunar orbits. It is speculated that orbits with regions of “predominantly conic” motion with respect to either the Moon or Earth are likely to result in greater success of the batch IOD and subsequent EKF convergence.
Finally, additional work should also consider ways to reduce the warm-up period of this navigator even further. For example, initializing a few hours or minutes prior to perilune arrival (or also initializing the state covariance matrix off diagonals with some a priori information) would be interesting cases to examine. It might prove to be worthwhile to find the absolute minimum time required for convergence. However, it is suspected that a full revolution would still be necessary to bring the actual error down to within the estimated uncertainties.
HOW TO CITE THIS ARTICLE
Hinga, M B., & Williams, D A. (2023) Autonomous lunar L1 halo orbit navigation using optical measurements to a lunar landmark. NAVIGATION, 70(3). https://doi.org/10.33012/navi.586
AUTHOR CONTRIBUTIONS
All authors contributed equally to this work.
CONFLICTS OF INTEREST
The authors declare no potential conflicts of interest.
DISCLOSURES
Approved for public release; distribution is unlimited. Public Affairs release approval No. AFRL-2021-4143
ACKNOWLEDGMENTS
The authors would like to thank the AFRL Scholars’ Program for funding Mr. Williams’ work on this project. We thank Robert Morris (of Wabtec) for advice on Kalman filtering.
5 APPENDIX
5.1 Development of the Coplanar System of Equations in the Inertial Frame (Two-Body Problem (f and g Series) - Phase One)
We define the kth inertial position vector of a target satellite for n observation times, t1, t2, …, tn as the following:
14
where 
is the known observer position, ρk the scalar range from the observer to the target satellite, and 
is the corresponding “line-of-sight” unit vector, which can be expressed in terms of two angles representing right ascension and declination. The observer can be on the Earth surface or on an orbiting satellite somewhere in space above the Earth.
For three observations, the following standard Gaussian equation expresses the relation of three vectors in inertial space as a linear combination summing to zero with three distinct coefficients, c1, c2, and c3 as
15
This was the equation used by Gauss when he predicted the position of the first minor planet after its conjunction with the Sun, (Herget, 1948). If we substitute Equation (14) into Equation (15) and separate all known quantities on the right side, we end up with
16
After setting c2 = −1 and after many algebraic manipulations, it is found that
17
18
where f1, f3, g1, and g3, are the so-called “Lagrange f and g coefficients” (Curtis, 2013). For this three-observation example, the approximation (without the velocity term) of the f and g coefficients are
19
20
21
22
where τ1 and τ3 are the time intervals between successive measurements of 
, and 
.
Just like Gauss, we put Equation (16) into matrix format. Because the unknown ranges appear on both sides of the equation, no closed form solution exists, forcing us to solve the system of coplanar equations through an iterative procedure. By starting with an initial guess of zero for the three unknown ranges we avoid forming an eighth order polynomial and instead, iteratively solve for the scalar ranges of ρ1, ρ2, and ρ3.
23
We do this by inverting the 3 × 3 matrix of known 
unit vectors, premultiplying against the right-hand side by this inverse, avoid forming and solving the traditional eighth-order polynomial for the scalar magnitude of the middle inertial position vector r2 (Vallado, 2001), and solve for ρ3, ρ2, and ρ3.
24
where
25
is defined as the state vector to be estimated.
At each iteration step, both the scalar ranges, ρi, to the target and the corresponding magnitudes, ri, of the inertial position vectors, 
, are also computed. Both are required at each iteration step (after inversion for the state vector) for computing the f and g coefficients. ρi is given by the relation, which describes the simple geometry of the measurement scenario.
26
Carrying this out to n observations, we still group the observations into sets of three where the kth relation is defined as
27
Here, the coefficients of ck and dk are obtained similarly as above (in the three-observation example) by expressing the vectors 
and 
in terms of position and (including) velocity vectors at time tk, rk, and vk using the Lagrange coefficients f and g in the following (Curtis, 2013)/(Karimi & Mortari, 2011) format.
28
29
By inserting these two expressions into Equation (27), vector vk can be eliminated and a relationship between 
, and 
is completely defined giving the expression for ck and dk (Karimi & Mortari, 2011) as
30
Notice that by eliminating the velocity term in the Lagrange equations above we are constrained to a “triplet” of three position vectors to compute the f and g coefficients. We can still increase the number of observations above three, however we will group the measurements into sets of three as we take on more observations to compute an orbit. This is discussed below in Section 5.2. The Lagrange coefficients fk and gk, appearing in Equation (30), can be expanded in series with a time difference Δtk = tk – tk−1 up to a fourth-order series expansion approximation as in (Curtis, 2013)
31
32
33
34
where 
is the Earth’s gravitational parameter.
For time intervals Δtk that are small in comparison with the orbital period, these coefficients f and g can be well approximated using the first two terms of the series expansion, which yield the approximate expressions (Curtis, 2013). These expressions for ck and dk will actually be used in experimentation because the velocity magnitude at the middle position vector, vk, is unknown.
35
36
where k = 2, …, n −1. For equally space measured times (Δt = constant), it is very easy to see that
37
5.2 Multiple Observations in the Inertial Frame (Two- Body Problem - Batch IOD - Phase One)
To extend the number of measurements to four or more using the Lagrange f and g coefficients in the coplanar system of equations, we must arrange them in groups of three so that the coefficients are still based on f and g being expressed as a combination of one initial position and one initial velocity vector. Since velocity was solved for in Equation (28) and inserted into Equation (29), we ended up with three position vectors, resulting in the relation in which the middle position vector is a linear combination of the first and third, as seen in Equation (27). Repeating this relation into a series of triplet measurements for a number n > 3 observations, the corresponding indices are seen as:
38
As we increase the number of observations, from k = 2,3,…, n –1 each triplet relation is preserved by considering the additional geometric line-of-sight unit vector as the “third” observation to form rk+1, complementing the previous two. As this additional position vector is included, the corresponding kth “right-hand side,” similar to that seen in Equation (16), called the “residual” ξk, is formed. The following logic illustrates this principle.
In Equation (27), we substitute for each of the three inertial position vectors with the position vector given by Equation (14), 
, expand the terms and organize the geometric line-of-sight unit vectors on the left-hand side and the position vectors of the observer on the right. (Note: the observer could be a terrestrial or spaceborne platform.)
39
On the left side Equation (39), the measurement “triplet” is shown as three known unit vectors pointing to the target satellite multiplied by unknown scalar ranges and Lagrange coefficients. The right-hand side (labeled as ξk) contains the corresponding three known inertial position vectors of the observing (satellite) platform (or the known position of the asteroid in the case of autonomous navigation) with the same unknown Lagrange coefficients. (The Lagrange coefficients are unknown because they are a function of the scalar magnitude of the inertial middle position vector of the target body.) Notice that the Lagrange coefficient of the “middle” geometric line-of-sight unit vector remains as negative one for subsequent sets of triplet observations. Starting with the very first set of n = 3 observations, k = 2, we have
40
Putting Equation (40) into matrix format leads to
41
Adding a fourth observation (n = 4, k= 3), we group the measurements into sets of three, yielding two equations:
42
To prepare this set of equations for a least squares solution (Herget, 1965), we leave them set equal to their own residual and “stagger” them into matrix form in Equation (43),
43
where 0 represents a 3 × 1 vector of zeros.
Writing Equation (43) in compact matrix notation, we have
44
Expanding this to n observations, the system of equations is given as
Notice that H* (for this batch DC filter, not to be confused with the H of the EKF) is a 3(n – 2) × n matrix with 
and 
, having the dimensions of n × 1 and 3(n – 2) × 1, respectively.
Instead of inverting matrix H* on the left side and iteratively solving for the state vector of unknown scalar ranges 
, as seen in the similar example of Equation (24), here we are able form the normal equation, Equation (45), treating 
as the residual with equal weights for all measurements in matrix H*, and solve for 
in Equation (46). This iterative approach is similar to that of Herget (1965). However, instead of trying to guess “good” starting distances for scalar ranges, ρk, initializing all of them to zero provides for excellent starting values and have always converged to the nontrivial solution (in this study) for orbiting platform observations. This result was also found to be true by Karimi & Mortari (2011).
45
46
5.3 Solving for Velocity: Lambert and the Fixed Time-of-Arrival Solution
After the vector of unknown ranges have been estimated, the series of inertial position vectors of the target spacecraft are calculated. To evaluate the spacecraft’s velocity, a Lambert solver is applied to the first and last inertial positions with the known fixed time of flight between them. The solver used in this study is that developed by the European Space Agency (ESA) for their 10 year “Rosetta” mission (ESA, 2004–2015). However robust and dependable this Lambert solver is, an improvement to its solution is obtained by implementing a two point boundary value problem (TPBVP) shooting method to “fine tune” the departure and arrival velocity at initial and final time, respectively. The guidance or “shooting” algorithm is based on the Linear Perturbation Theory (Battin, 1960) developed for the America’s Program for Orbiting Lunar and Landing Operations (APOLLO) program during the 1960s. This TPBVP guidance algorithm can be applied equally well in either the rotating or inertial frames, where the Lambert method serves as a guess for the initial/final conditions in the inertial frame.
Using the notation of Battin (1960), the relation of the final state (position and velocity) error of the target ballistic spacecraft to the deviation of the current state from the nominal at time t is given by Equation (48). The term Φ(tf, t0) is known as the state transition matrix and is formed by taking the partial of the final state with respect to the initial state, (Gelb et al., 1974). This derivative is a matrix of partials, commonly known as the “Jacobian”, shown as matrix A in Equation (47). In this study, the Earth’s J2 gravitation model is used to define the inertial accelerations in this derivative and for the accelerations in the rotating frame derivative, the CR3BP accelerations are used. Integrating this equation produces the state transition matrix, which relates the change in state from a certain time t0 to another time tf. The homogenous solution, for a fixed time of integration, is given in Equation (48) and yields a fixed time-of-arrival solution.
47
48
49
Let us expand Equation (48) in terms of position “r” and velocity “v” components and assign convenient labels to the portions of the Φ matrix. Because we know the position of where the spacecraft starts and do not want to vary it, we set δr (t0) = 0. Then, the expression for the perturbation to the position and velocity at final time tf is,
50
51
Thus, the equation that defines how to vary (or perturb) the spacecraft velocity at time t0, based on the missed distance at the target impact, is
52
where R is the upper right 3× 3 matrix of the state transition matrix Φ(tf, t0).
The velocity correction defined in Equation (52) is used as the iterative correction term at initial time t0 in the search for the optimal improvement of the initial velocity to minimize the miss distance at tf. Upon convergence of this shooting/guidance method, the Lambert velocity solution has been improved, allowing for more accurate orbital elements to be evaluated in the inertial frame and for the improvement of the initial state velocity in the synodic frame of the CR3BP problem.
5.4 Extended Kalman Filter Formulation - Phase Two
The Extended Kalman Filter (EKF), a workhorse of real-time spacecraft state estimation (Crassidis et al., 2007), is a recursive nonlinear estimator (perturbed by Gaussian noise) that is discretized in the time domain by linearizing the physical dynamics of the current best estimate of the parameters of interest. It is not a requirement that the model of the state dynamics 
and observation model G [ X(t), t] be linear, only that they are differentiable. This means that the formulation of the state transition matrix Φ(t, t0) and the matrix H[X(t), t], which is the partials matrix needed to compute the predicted measurement from the predicted state (discussed below), is defined. In a non-Extended Kalman Filter, the same is true, but the linearization occurs for some precomputed nominal trajectory of the state. In the EKF, the current best estimate comes from an optimal combination of the state from the previous time step (that propagated to the current time) and the current measurement. The definition of this combination is determined by the so-called Kalman gain K, which is defined below.
Kalman filters are unusual in that most filters (i.e., Butterworth filter) are formulated in the frequency domain, then transformed back into the time domain for application. The EKF can be considered an adaptive low-pass infinite impulse response (IIR) digital filter, meaning that its response to an impulsive input is nonzero for infinite time (Gelb et al., 1974). The frequency response of the EKF in this study is of no interest.
Ideally, if the model of the state and measurements are complete and accurate and perpetrate no acts of error omission or commission, then the covariance P(t) of the estimate state will accurately reflect the confidence of the estimated state vector, and those parameters will have a mean error of zero. Stated differently, the variance and covariance of the estimated state parameters will have a distribution about the true state. Invoking the expectation operator 
, where f(ξ) is the function of interest, if there are no biases in the estimate 
and residual 
, then:
53
and the covariance matrices for the state estimate and residual, defined as,
54
will have zero bias. However, since the filter of this study is intentionally mechanized as a suboptimal filter, small biases will be present and the Equations of 53 will be close to zero.
The state of the filter is represented by 
, the estimated state at time k, and the error covariance matrix Pk, which is a measure of the confidence in that state estimate. The EKF has two separate phases, which are called prediction and update. In the prediction phase, the estimate from the previous time step (k–1), both the state and covariance matrix, are propagated forward to the current time step (k). Then, during the update phase, the state is refined with measurement information from the current time step. It is intended that after this refinement, the new estimate of the state is more accurate, i.e., closer to the truth. In this study (note that there is neither a control model nor a control input vector), the equations for these two phases are as follows (Kalman, 1960).
Predict Phase
where Φk,k−1 and Qk are the state transition and process noise matrices, respectively. In this investigation a Runge-Kutta 4 integration scheme (and when needed, a variable step size integrator with an appropriate tolerance) is used to propagate both Φk,k−1 and 
forward one time step to give 
, using the Jacobian matrix (a matrix of partial derivatives).
Update Phase
The terms 
and Hk are the measurement model and the partials of the measurement model with respect to the state, and are defined in Section 1.3. Yk is the actual measurement taken by a camera to produce a geometric line-of-sight unit vector and is discussed in Section 5.1.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.