A Distributed Radio Beacon/IMU/Altimeter Integrated Localization Scheme with Uncertain Initial Beacon Locations for Lunar Pinpoint Landing

As a growing number of exploration missions have successfully landed on the Moon in recent decades, ground infrastructures, such as radio beacons, have attracted a great deal of attention in the design of navigation systems. None of the available studies regarding integrating beacon measurements for pinpoint landing have considered uncertain initial beacon locations, which are quite common in practice. In this paper, we propose a radio beacon/inertial measurement unit (IMU)/altimeter localization scheme that is sufficiently robust regarding uncertain initial beacon locations. This scheme was designed based on the sparse extended information filter (SEIF) to locate the lander and update the beacon configuration at the same time. Then, an adaptive iterated sparse extended hybrid filter (AISEHF) was devised by modifying the prediction and update stage of SEIF with a hybrid-form propagation and a damping iteration algorithm, respectively. The simulation results indicated that the proposed method effectively reduced the error in the position estimations caused by uncertain beacon locations and made an effective trade-off between the estimation accuracy and the computational efficiency. Thus, this method is a potential candidate for future lunar exploration activities.


Introduction
Safe and soft pinpoint landing (within 100 m at 3σ from the target site [1]) on an extraterrestrial body has been a central objective since the beginning of human space exploration missions. To date, many manned and unmanned landers have successfully landed on moons (Apollo and Chang'E), planets (Curiosity and Opportunity) and asteroids (Rosetta and Hayabusa-2), with landing footprints in the scale of kilometers [2,3]. Among all these celestial bodies, a growing number of exploration missions will land on the Moon in the future, as it is the most suitable outpost for deep space exploration. A robotic lunar base followed by a human base will likely be constructed during these missions, which will be the first potential implementation of co-located pinpoint landing [4]. With the existence of these infrastructures, the natural consequence is to make use of more ground sources to enhance the current onboard navigation capability for lunar landing missions using, for instance, radio beacons.
Radiometric localization is an extensively used technology on Earth. Although there are many specific kinds of signals and modulations for different applications, the information provided by radio beacons can be divided into three categories [5]: (a) range measurements, (b) range rate measurements and (c) bearing measurements. All of these have drawn attention to improving the accuracy of pinpoint landing. In 2008, NASA proposed a radio measurement-enhanced navigation architecture based on lunar relay satellites (LRS) and lunar communication terminals (LCT) for the Autonomous Landing and Hazard Avoidance Technology (ALHAT) project [6]. By combining range and range rate measurements from LRS and LCT with other basic measurements of the ALHAT sensor suite, the landing accuracy could reach a level of 10 m at 3σ [7].
Several studies regarding the optimization of beacon configurations have been presented, both in terms of the filter accuracy [7,8] and observabilities [9,10]. On this basis, Theil et al. from the German Space Center (Deutsches Zentrum für Luft-und Raumfahrt, DLR) investigated the impact of one to four ground beacons for the small integrated navigation for planetary exploration (SINPLEX) project of ESA [4,11]. In this research, bearing measurements, as well as range and range rate measurements, were utilized by the navigation filter to augment the navigation baseline of SINPLEX. Different cases of beacon locations and setups were compared in the Autonomous Terrain-Based Optical Navigation (ATON) project of DLR, which indicated that the impact of additional bearing measurements was quite small.
All studies listed above were based on a strong assumption that uncertainties only existed in the sensor measurements, which means that all beacons were stationary and their locations were precisely pre-determined. In fact, most ground beacons can only be located either by orbiters or by the deep space network on the order of hundreds of meters [12]. Without the protection of the atmosphere, the diurnal temperature variation on the Moon is more than 300 • C. If radio beacons are stationary, meaning that they will be exposed to a harsh environment for a long time, the onboard electronic equipment will be fatally damaged. One possible solution is to maintain radio beacons in the lunar base and redeploy them before each landing mission. Such scenarios serve as a strong motivation to investigate a novel navigation scheme that can locate the lander and update all beacon locations at the same time, which is similar to the range-only simultaneous localization and mapping (RO-SLAM) approach in the field of robotics.
Almost all existing radio beacons/inertial measurement unit (IMU)-integrated navigation algorithms in pinpoint landing are formulated using the state by the covariance matrix Σ and the mean vector µ of the multivariate Gaussian distribution and tracked via the extended Kalman filter (EKF) [4,11,[13][14][15] or unscented Kalman filter (UKF) [16,17]. However, the computational complexity and memory space of these algorithms are both quadratic to the dimension of state vectors [18], which makes them unsuitable for dealing with this SLAM-like localization problem.
A large amount of strategies for RO-SLAM [19][20][21][22][23][24] have been developed in recent decades, in which the sparse extended information filter (SEIF) proposed by Thrun et al. [18] has been widely used. Unlike Kalman filters, SEIF represents the state by the information matrix Λ = Σ −1 and the information vector η = Λµ, which makes the update stage more efficient due to the natural sparseness of the information matrix [25]. A conditional sparsification operation is implemented before the prediction stage to maintain the sparse representation of Λ, which leads to a constant execution time online.
While SEIF has been demonstrated to be more efficient and scalable in many applications [26][27][28], the consistency and accuracy of SEIF may be worse than EKF due to three approximations [29]: (a) the linearization of the motion and measurement models, (b) the recovery of the mean vector by iterations and (c) the sparsification via approximated conditional independence. Several methods have been designed to refine the performance of SEIF, such as exactly SEIF (ESEIF) [30], exactly sparse delayed-state filter (ESDSF) [31] and hybrid SEIF (HSEIF) [32]; each method only addresses one aspect of the approximation errors.
Iterated versions of EKF and UKF can increase their consistency and robustness against approximation errors compared to the original versions [33]; however, the rate of convergence is frustrated due to the dense covariance matrix. Due to the sparseness of the information filter, the computational burden is no longer a problem for iterated SEIF. He et al. [34] proposed the first iterated SEIF (ISEIF) for autonomous underwater vehicles, which utilized the Gauss-Newton algorithm for measurement updates rather than the traditional linear combination of predictions and innovations. However, divergence may occur in the Gauss-Newton algorithm when the initial values are far from the optimal values or the coefficient matrix in the normal equation is approximately singular [35].
In addition to the efficiency, the update stage in SEIF is also additive, which makes it capable of deployment in distributed systems. Torres et al. [23] first proposed a distributed RO-SLAM scheme based on SEIF for an autonomous ground vehicle (AGV). The prediction stage and update stage are implemented on the AGV and surrounding beacons, respectively. Except for the traditional direct robot-beacon measurements, additional inter-beacon measurements are also integrated into the SEIF, which reduces uncertainty both in the map estimation and localization accuracy. This work was then extended to a 3D resource-constrained operation with an auxiliary selection tool which only integrates the most informative measurements into the SEIF [36]. Further analysis regarding the effect of inter-beacon measurements on the building of an information matrix followed in [37], validating the preservation of the main properties of SEIF.
Beacon initialization is another challenge in RO-SLAM due to the partial observability of range measurements. As a multi-mode problem, the initial estimation of beacon locations can be obtained either by delayed or undelayed initialization methods. The trilateration method [38][39][40] is the most simple and efficient delayed initialization method; however, it suffers from measurement noise and the odometry error of robots. Olson et al. [41] attempted to address this problem with a 2D probability grid, with a performance that was limited by the size and resolution of the grid.
The shortcomings of these two methods resulted in them being quickly replaced by another widely used delayed initialization method: particle filters (PFs) [23,36,37,42,43]. PFs model the probability density of the beacon location as a uniform distribution around the measurement circle, and then the weight of each particle is updated and finally converges to a single solution by accumulating successive range measurements. Even though PFs can provide more accurate results, their significant computational burden cannot be neglected, which leads to several corresponding investigations [44,45].
In contrast to the delay initialization method, the initial beacon location is modeled with a multi-hypothesis distribution, which makes it available to integrate range measurements from the start. Blanco et al. [46] first proposed a Gaussian mixture model (GMM)-based initialization strategy, which provides accurate results but makes the integration of inter-beacon measurements impossible as each beacon is inserted in an independent EKF. On the basis of this research, Felipe et al. [47] proposed a centralized EKF-based initialization framework, which allowed the integration of inter-beacon measurements, and they enhanced it with a outlier rejection filter [48,49]. Geneve et al. [50] proposed a short-delayed composite initialization method based on Euclidean parameterization and a two-hypothesis GMM, which showed a lower computational cost.
Based on all the existing studies listed above, a novel distributed radio beacon/IMU-integrated localization scheme based on an adaptive iterated sparse extended hybrid filter (AISEHF) is presented in this paper. To the best of our knowledge, this is the first time that beacon location errors have been taken into consideration in the integrated navigation architecture for lunar pinpoint landing. This scheme is inspired by the SEIF-based distributed framework sketched in [37] and the ISEIF proposed in [34], while several aspects have been modified according to the characteristics of lunar pinpoint landing: • Only lander-beacon range measurements are utilized as inter-beacon measurements and are easily blocked by obstacles, such as craters on the lunar surface. • A batch-least-squares trilateration algorithm similar to [39] was adopted for the initialization of beacon locations, as their prior rough estimations are available and the inertial measurement unit (IMU) used in our application is much more accurate than that used in the majority of RO-SLAM applications.

•
The navigation state between two successive measurements was propagated by a hybrid form of the "mean vector + information matrix" to avoid the frequent conversion from the information vector to the mean vector, which can improve the computational efficiency of the prediction stage.

•
An adaptive iteration algorithm with the damping factor was adopted in the update stage, which can both improve the accuracy and efficiency of the estimator.
The remainder of this paper is structured as follows. The problem formulation is given in Section 2. In Section 3, we detail the implementations of the proposed distributed localization scheme. In Section 4, we verify the performance of the proposed algorithm by simulations, which is followed by further discussions. Eventually, our conclusions are presented in Section 5.

Problem Formulation
As illustrated in Figure 1, the target application scenario of this paper is the final phase of the power descent, which is short both in terms of the flight time and downrange distance (commonly around 200 s and 5-10 km). Without generality, the landing frame (L) can be chosen as the navigation frame. The origin of the L frame o l is located at the targeted landing site, and three coordinate axes (x l , y l , and z l ) are aligned with the geographic directions east, north and up (ENU), respectively. Another significant coordinate system is the body frame B, whose origin lies in the mass center of the lander, and the directions of the three axes are front-left-up (FLU); i.e., the x b axis points to the front along the body symmetric axis, the y b axis points to the left and is orthogonal to x b and the z b axis completes the right-handed orthogonal coordinate system. In addition, measurements of the laser altimeter are present in the sensor frame (S), whose origin and x s axis are the same as the B frame but whose y s axis and z s axis are opposite to y b and z b .  Overview of the proposed distributed radio beacon/inertial measurement unit (IMU)-integrated localization scheme and definitions of the landing frame L (east, north and up (ENU)) and the body frame B (front-left-up (FLU)).

Lunar Surface
The kinematic equations of the lunar lander, defined with respect to the landing frame L, are given as follows: [51]: where [·] ∧ is the skew symmetric operator; p L is the lander position; v L is the lander velocity; C L B is the rotation matrix from B to L; f B i and ω B ib , respectively, denote the specific force and the angular rate sensed by IMU, whose further interpretations can be found in [52]; and g L is the gravity vector, whose magnitude is assumed to vary with the altitude as follows [51]: where g 0 = 1.622 m/s 2 is the gravity constant at the lunar surface, h is the current altitude of the lander and R m is the radius of the Moon. The superscripts that denote the corresponding frame are omitted for notational convenience in the following parts if there is no confusion. The beacon receiver is attached on the bottom of the lander, which makes it directly visible from the ground beacons. Several radio beacons are mounted on rovers and pre-deployed along the designed ground track of the landing mission, while their prior rough locations are determined by the onboard navigation systems of the rover. Due to the limited effect of bearing measurements, each beacon is assumed to be equipped with a range-only sensor whose measurements r li are affected by independent Gaussian white noise v b : , denote the positions of the lander and the ith beacon. In addition, only the lander-beacon range measurements are available, as there are too many obstacles (such as craters or huge rocks) on the lunar surface, which could heavily degrade the inter-beacon measurements. Considering the terrain accessibility of the rover, all ground beacons are assumed to be on the same horizontal plane. This coplanarity will lead to uncertainty in altitude estimation; thus, a laser altimeter is mounted next to the beacon receiver to offer additional range measurements r a : where θ is the pitch angle, γ is the roll angle and v a is the Gaussian white noise. The flow diagram of the proposed distributed localization scheme is shown in Figure 2, which can be divided into two parts: the lander and the beacon. The prediction stage is continuously conducted onboard based on the IMU measurements. After the visible beacons have been detected, a message that contains the current estimated lander positionp l will be transmitted to all visible beacons. When this message has been received by beacon b i , it will switch from the "Standby" mode to the "Initialization" mode, in which the beacon location is initialized by minimizing the following batch-least-squares cost function with N range measurements: where N can be pre-determined offline and Equation (5) can be solved by the nonlinear optimization methods given in [53] with prior rough beacon locations as initial values. When the initialization of b i is convergent, it starts to work in the "Localization" mode. Its contribution (local information vectorη i t and information matrixΛ i t ) to the update stage is calculated and sent back to the lander. Then, the locations of all visible beacons are inserted into the state vector, and their contributions are summed together as follows: where( ·) and( ·) denote the variables obtained by the state prediction and measurement update, respectively. S i = 0 · · · I · · · 0 T is a projection matrix that allocatesη i t andΛ i t at the corresponding position ofη t andΛ t : Finally, an iterated update was performed based on the global information vectorη t and information matrixΛ t , which was followed by the sparsification step of SEIF when the number of visible beacons were out of bounds or some links in the information matrix were too weak [34].

State Definition and Jacobians
Like EKF, the first-order Taylor series expansion was applied to linearize the state space model as follows in Equation (8) to deal with nonlinear propagation and measurement functions, where w t and v t+1 are the white process noise and measurement noise, respectively. , · · · , p T b n ] T . As a navigation (or localization) filter, it is common to choose the specific force f i and gravity vector g as the control input u t = [ f T i , g T ] T . As B is kept unchanged and is independent of x t , we compute the Jacobians F and G from the kinematic model (Equation (1) with respect to the x t , which is more efficient. The results are given as follows: Similarly, based on the measurement function of Equation (3) and (4), the measurement Jacobian H is given in Equation (11)∼(13), where all entries that are not given are zero.

Temporal Alignment of Asynchronous Measurements
Signal transmission or triggering delays will lead to a temporal misalignment between measurements from different sensors, which cannot be ignored in the integrated localization scheme. As depicted in Figure 3, a virtual measurement method based on linear interpolations was adopted in our scheme to eliminated the temporal misalignment. Assuming there are j IMU measurements between every two sequential beacon measurements and the two closest IMU measurements to t k are denoted as m i t k−1/j and m i t k/1 , the virtual measurement at t k can be obtained as follows: When all measurements are temporally aligned, the state prediction and measurement update stage of SEIF are then executed separately.  Figure 3. An example of processing asynchronous inertial measurement unit (IMU) and beacon measurements.

Hybrid State Prediction
The typical prediction stage in SEIF can be separated into two main processes [31]-state augmentation and state marginalization-assuming the state vector ξ t is subjected to the conditional distribution given, z t and u t−1 , as expressed in Equation (15): At time t + 1, a new variable x t+1 , which represents the newest lander state, is first inserted into the state vector ξ t , and the augmented distribution p (x t , x t+1 , B) ∼ N −1 η t+1 , Λ t+1 is given as follows: The propagated distribution p (x t+1 , B|z t , u t ) can then be recovered by marginalizing out x t from the joint distribution p (x t , x t+1 , B|z t , u t ): From Equations (16) and (17), the recovery of µ t from η t must be performed before the state augmentation to re-linearize F and G [18]. This is not a problem for a synchronous application, as η t will be utilized in the following update stage, whereas in an asynchronous application, the prediction stage will be followed by several other prediction stages, meaning that the state recovery will be implemented several times. As the state recovery is commonly performed based on iterated algorithms, as with a pre-conditioned gradient, these extra steps will lead to an unnecessary computational burden.
Based on the analysis above, we proposed a hybrid form of the "mean vector + information matrix" for the prediction stage, which will be called the sparse extended hybrid filter (SEHF) in the following. SEHF propagates the uncertainty of ξ t using Λ t as with the typical SEIF, while the expectation of ξ t is propagated by µ t instead of η t as follows: Figure 4 shows a comparison of the typical form and our hybrid form prediction when there are two prediction stages between successive update stages. Clearly, in the hybrid form, the state recovery only needs to be performed once before the filter switches from the update stage to the prediction stage, regardless of the number of extra prediction stages, whereas it will increase with the number of extra prediction stages in the typical form. Although there is an extra operation in the hybrid form that constructs η t for the update stage, the corresponding computational cost is quite small, as this is a simple matrix multiplication.

State Recovery
Update Stage Prediction Stage Figure 4. Comparison of the typical form and our hybrid form prediction. Different processing steps between these two forms are marked with different colors depending on their computational burden.

Iterated Measurement Update
When the current measurement z t+1 is achieved, the SEIF typically updates the new posterior p (ξ t+1 |z t+1 , u t ) according to the chain rule based on existing distributions p (ξ t+1 |z t , u t ) ∼ N −1 η t+1 ,Λ t+1 and p (z t+1 |ξ t+1 ) ∼ N (h(ξ t+1 ), R t+1 ): Varying from solving the measurement update via a linear combination of the prediction and innovation in the traditional filter, the iterated filter models it as a nonlinear least squares (NLS) problem as follows: where χ 2 (ξ t+1 ) is the objective function of the NLS problem and ξ * t+1 is the optimal value when χ 2 (ξ t+1 ) takes the minimum value, defined by the following matrices: The objective function χ 2 (ξ t+1 ) can be rewritten into the following a matrix form [33]: Then, the iterative sequences of ξ i t+1 and Λ i t+1 can be achieved as follows [34]: where H i is the Jacobian of the measurement function at ξ i t+1 .

Iterated Measurement Update with Damping Factor
The iterative form update equations given above are obtained by solving the following normal equation of the Gauss-Newton algorithm: where J i represents the Jacobian of e(ξ t+1 ) with respect to ξ t+1 evaluated at ξ i t+1 . As the coefficient matrix J T i J i is only positive semidefinite, there is no guarantee that ∆ξ is always along the descent direction. ∆ξ would become unstable when the initial value was far away from the solution as the second derivative (Hessian) of the object function was approximated by the Jacobian. To address these issues, a damped Gauss-Newton algorithm was proposed by Levenberg and Marquardt, successivel [53]. In this method, a positive damping factor λ was added to the coefficient matrix of the normal equation, which is known as the damping factor: Clearly, J T i J i + λI is always positive and the method is equivalent to the steepest descent algorithm and the Gauss-Newton algorithm when µI J T i J i and µI J T i J i , respectively. With this modification, the iterative sequence of ξ t is modified as follows (see Appendix A for more details): As the uncertainty of ξ t will only be updated when the iteration converges [33], the update of Λ t is the same as Equation (28). As a hybrid form is adopted for the following prediction stages in our scheme, there is no need to convert µ t to η t , which further reduces the computational burden.

Damping Factor and Stopping Criteria
A proper value of λ is vital to the iteration. One possible solution is to choose the initial value of λ depending on the maximum value of the coefficient matrix: and then update λ based on the gain ratio ρ at each iteration [54]: where ρ is the ratio of the actual and approximated decrease of the objective function: Another parameter of concern is the number of iterations, which balances the accuracy and efficiency of the iteration process. Commonly, this parameter is difficult to predetermine offline; thus, several adaptive stopping criteria are applied in the proposed algorithm. First, the iteration process should be terminated when the difference between two successive residuals is smaller than a threshold 1 , which reflects the global minimum: Another stopping criterion is related to the local minimum; i.e., the increment ∆ξ i t+1 is smaller than the threshold 2 : Lastly, the maximum iterations k max should be settled in the case of a dead loop. With all notations and definitions above, the complete working flow of the adaptive iterative update is summarized in Figure 5.

Simulation Scenario and Parameters
The simulation scenario is a final landing phase performed after the main break phase of the powered descent, as shown in Figure 6.  Figure 7 presents the profile of the specific force and angular rate in detail. Ten beacons were deployed uniformly along both sides of the ground track in the simulation, and their true locations are listed in Table 1. An initial position error of 200 m was added to their true positions to simulate the prior rough locations obtained by rovers.
With respect to the navigation filter, the number of measurements K in the beacon initialization in Equation (5) was chosen to be 50. The attitude was propagated based on raw gyroscope measurements and periodically calibrated by the star tracker independently with a residual around 9.1" [4]. The residuals of the initial position and velocity were set to 300 m and 30 m/s, respectively. The 1σ parameters of IMU, the laser altimeters and ground beacons are listed in Table 2, and their update frequencies are, respectively, settled to 200 Hz, 100 Hz, and 20 Hz.   Other parameters of the localization filter-i.e., the initial covariance matrix P 0 , the process noise covariance matrices Q t and the measurement noise covariance matrices R t+1 -are given as follows: P 0 = diag 1 × 10 4 , 1 × 10 4 , 1 × 10 4 , 1 × 10 2 , 1 × 10 2 , 1 × 10 2 , 1 × 10 4 , · · · , 1 × 10 4 R k+1 = diag 1 × 10 4 , · · · , 1 × 10 4 , 25 where n is the number of available beacons.
The simulation system was implemented in C++ using the Robot Operating System (ROS) on a laptop with an Intel Core TM i7-4790 @ 3.68 GHz and 8 GB RAM. Each sensor runs an independent ROS node to simulate the real interactions.

Simulation Results and Analyses
The estimation errors of the navigation state obtained by the proposed AISEHF-based distributed localization scheme are shown in Figure 8. It can be seen from the results that the residuals converged into the 1σ uncertainty bound rapidly both in the triaxial position and velocity estimations with steady values less than 50 m and 5 m/s, which meet the demands of pinpoint landing.  Then, the proposed scheme (Scheme 3) was further compared with another two schemes with the same sensor configurations: (a) the SEIF-based scheme proposed in [37] without inter-beacon measurements (Scheme 1) and (b) the ISEIF scheme proposed in [34] (Scheme 2). The performance of all three schemes was evaluated in 100 independent Monte Carlo runs. In each of the 100 runs, all random seeds (e.g., beacon locations and sensor noise) were changed randomly, while all parameters of estimator were kept unchanged. Two indicators-root-mean-square errors (RMSEs) and averaged RMSEs (ARMSEs), which are defined in Equation (42)-were adopted to acquire a fair result.
where K is the number of Monte Carlo runs and T is the total time from the convergence to the end.
x i t and x i t , respectively, denote the estimated and true state at time t in the ith run. Figure 10 gives the position RMSEs and the average computational cost of the three schemes mentioned above from 100 independent Monte Carlo runs. The improvement of Scheme 3 in the position estimation was 4.55% w.r.t. Scheme 2 and 8.11% w.r.t. Scheme 3. The velocity RMSEs are not given here due to the fact that the iteration algorithm mainly focused on addressing the linearization error caused by the linear approximation of the measurement function, and it can be seen from Equations (12) and (13) that triaxial velocities make no contribution to the measurement function in this case.
The computational cost of Scheme 3 was 3.64% more w.r.t. Scheme 1 and 18.96% less w.r.t. Scheme 2. The performance of the three schemes is summarized in Table 3, indicating that Scheme 3 results in a good trade-off between position estimation accuracy and computational efficiency.  The influence of the beacon number n is given in Figure 11. From the results, the position ARMSE continued to decrease when n was lower than 10 and increased to a large error in the order of kilometers, while the computational cost continued to increase with n. We defined an indicator call cost-ARMSE percentage, which was equal to 100/(Cost * ARMSE) to further illustrate the results and indicated that n = 10 was the most suitable configuration in our case. Figure 12 gives the comparison of the proposed scheme under different initial location errors. We can infer from the results that the batch-least-square trilateration initialization algorithm achieved a consistent residual (around 40 m) regardless of the beacon initial location error in our application, which also led to a good robustness of the proposed localization scheme to the beacon initial noise. Extensive results with the initialization procedure disabled are further illustrated, which further prove the conclusion. Finally, another 100 independent Monte Carlo runs were performed to evaluate the final drops of the proposed method, and the results are plotted in Figure 13, showing a circular error probable (CEP) of 37.52 m and a mean value of [35.66 m, −12.17 m]. Another interesting point is that the distribution of the final drops in the north (11.60 m) was much larger than that in the east (1.58 m). One possible reason for this phenomenon may lie in the geometrical configuration of all ground beacons, which has a significant effect on the observability of the whole system.

Conclusions
In this article, a distributed localization scheme was proposed to handle the situation of having no accurate prior knowledge regarding the beacon configuration in lunar pinpoint landing. This scheme integrated onboard IMU/altimeter measurements together with lander-beacon range measurements within a SEIF-based framework. The prediction and update stages of SEIF were replaced by a hybrid-form propagation and a damping iteration algorithm, respectively, which led to the so-called AISEHF. The simulation results indicated that the proposed scheme could satisfy the navigation demands of lunar pinpoint landing and showed a robustness to uncertain initial beacon locations. The Monte Carlo simulation indicated that our scheme makes a good trade-off between estimation accuracy and computational efficiency compared with the existing algorithms.
Although the research in this paper has proved to be significant, there are still some questions that require further investigation, such as the augmentation of the current range-only measurements with the range-rate and the optimization of the geometrical pattern of all ground beacons. Several realistic problems, such as transmission loss and glint noise, should be taken into consideration and verified by experiments. Our work will be expanded to include these issues in the future.

Conflicts of Interest:
The authors declare no conflict of interests.

Appendix A. Adaptive Iterated Update Algorithm
The detailed matrix form of J i is given as Substituting Equation (A1) into Equation (30) leads to where S =Λ t+1 + λI. Expanding the above equation with the matrix inversion lemmas given in [33], we can achieve