INS/GNSS Integrated Rover Navigation Designed With MDPO-Based Dual-Satellite Lunar Global Navigation Systems

Aiming to provide navigational framework to upcoming early lunar exploration missions by lunar rovers, interests in robust and low-cost global satellite navigation systems (GNSS) around the Moon are growing more than ever before. In response, dual-satellite lunar global navigation systems (LGNS) based on Multi-epoch Double-differenced Pseudorange Observation (MDPO) algorithm was proposed in our earlier work. While the MDPO-based dual-satellite LGNS can provide reasonably high positioning accuracy at an order of several tens of meters within a one-minute observation, the positioning calculation is only available when the two navigational satellites are in user’s view. This limitation can be overcome by integrating other navigational sensors such as inertial navigation system (INS) and compensating for user’s positions in the absence of GNSS signals. The main objective of this research is to provide an integration model of INS and MDPO-based dual-satellite LGNS measurements and quantitatively show the benefits of the proposed integration by numerical simulations. More specifically, the major contributions of this paper are three-fold: 1) proposed a mathematical model of INS/GNSS integration adopted for the MDPO algorithm, 2) developed a numerical simulation that combines INS measurements and dual-satellite LGNS measurements, and 3) performed a quantitative comparison between the proposed INS/GNSS integration and raw dual-satellite LGNS measurements.


I. INTRODUCTION
A great deal of research conducted on the design of the Global Navigation Satellite Systems (GNSS) that use a fewer number of satellites in the applications for Earth GNSS [1]- [6] and applications for lunar GNSS [7]- [11]. The author's prior research collectively studied different algorithms for lunar GNSS that comprises only two satellites, i.e., dualsatellite, and made a comparative analysis by summarizing their pros and cons [12]. Dual-satellite lunar global navigation systems (LGNS) are a cost-efficient and useful platform to provide multiple users with reasonably accurate position measurements, such as a few tens of meters within one-minute observation as reported in [11]. One of the shortcomings of dual-satellite LGNS is lower availability (availability: a percentage of time at which the GNSS The associate editor coordinating the review of this manuscript and approving it for publication was Ehsan Asadi . system provides the number of satellites required for position fixing). Reportedly, in the case that two navigation satellites are injected into low lunar orbiters, availability is capped at around 15 %. Consequently, rover's traveling distance over the mission period will be limited by the availability of GNSS measurements in case the rover solely relies on GNSS measurements.
As an approach to extend the traveling distance, integration of inertial navigation system (INS) can be considered: when observation data from navigation satellites are not enough to pinpoint the user position, INS measurements are used to obtain rough user position estimation. During the last decade, different approaches for INS/GNSS integration have been adopted [13] and many of them have been investigated for different types of applications, for instance, urban canyon [14], [15], indoor [16], maritime [17], airborne [18] as well as lunar exploration [19]. In [19], an autonomous crewed lunar landing mission was simulated where range measurements, VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ an inertial measurement unit (IMU), and a star tracker were used to compensate for navigational measurements from single lunar navigation satellite. Furthermore, the complementary characteristics of GNSS and INS have made them well-suited to integration. In particular, INS obtains better short-term accuracy; however GNSS can provide navigation solution with long-term stability. In order to cope with demands from missions that requires higher accuracy than several tens of meters, INS/GNSS integration can be used to improve the user position estimation accuracy.
In this work, we provide an INS/GNSS integration designed to work with a dual-satellite LGNS that is based on the multi-epoch double-differenced pseudorange observation (MDPO) [11]. We present a mathematical model and numerical simulation to quantitatively show how the proposed integration will improve the aforementioned shortcomings. In the numerical simulation, rover trajectory and dual-satellite LGNS measurements were produced using our prior research work [11]. As for INS measurements, our analysis assumes that the rover is equipped with sensors that measure the rover's orientation and traveling distance. Those measurements are commonly used as navigational information in many planetary rovers. For example, rover's orientation can be derived from a gyro [20], a sun sensor [21], [22] or star tracker [23], [24], and traveling distance can be derived from a wheel odometry [25], [26], lidar [27], [28] or visual odometry [29]- [31].
To our best knowledge, this research provides an INS/GNSS integration for dual-satellite LGNS applications for the first time. Therefore, the obtained numerical simulation results are highly useful to the discussion of the future LGNS design. Moreover, the provided integration form is specialized for the MDPO algorithm, which has not been presented in any other work.
To summarize, our major contributions include: 1) A mathematical model of INS/GNSS integration adopted for the MDPO algorithm, 2) A numerical simulation that combines INS measurements and dual-satellite LGNS measurements, 3) A quantitative comparison between the proposed INS/GNSS integration and raw dual-satellite LGNS measurements under the simulated conditions. The rest of the paper is organized in four sections. In Section II, we review the mathematical model of the MDPO algorithm, and fuse it into INS/GNSS integration form. In Section III, we summarize a simulation scenario, simulation setting, and the simulation results. In Section IV, we presents discussion of the results. In Section V, we provides some concluding remarks.

II. METHODS
In Section II-A and II-B, we refer relating mathematical models from previous research with respect to the MDPO algorithm and INS/GNSS integration. In Section II-C, we develop a mathematical model of the INS/GNSS integration specialized for the MDPO algorithm to the purpose of this research.
In Section II-D, II-E, and II-F, we describes systematic error, uncertainty in modeling parameters, and other considerations related to algorithm implementation on real hardware.

A. MULTI-EPOCH DOUBLE-DIFFERENCED PSEUDORANGE OBSERVATION
This section reviews the multi-epoch double-differenced pseudorange observation (MDPO). In this paper, we only explain important equations which are finally fused into the INS/GNSS integration. The complete formulation derivation of the MDPO can be found in [11].
The MDPO reduces the number of navigation satellites required from four to two, while dealing with the clock instability of the space segment and user segment as well as satellite orbit determination error, as shown in Fig. 1. In the conventional Time Of Arrival (TOA) algorithm, the pseudorange (ρ) measurement between one user (user1) and one satellite (satellite1) is presented by the following equation: where X X X s (t s i ) is the satellite1 position at the time of signal transmission t s i ; X X X R (t i ) is the user1 position at the time of signal reception t i ; c is the speed of light; dτ R is the user clock bias; dT S is the satellite clock bias; dX X X Rsa corresponds to the user1 position transition due to the Sagnac effect; and ω r S R is the pseudorange receiver noise. In this study, we assume that the receiver noise ω r S R follows a white Gaussian distribution with the standard deviation of σ ωr .
MDPO uses double-differenced pseudorange observations to eliminate the clock bias of the space segment and user segment as well as satellite orbit determination error, by subtracting four pseudorange measurements between two users (user1 and user2) and two satellites (satellite1 and satellite2) as where ∇(•) is the double difference operator. In the double difference method, user2 is used as a reference station whose position is fixed and known, and the position of user1 is estimated in relation to the position of user2, i.e., user2's position is referenced as the origin of navigation. In a lunar navigation system, the lander can be used as a reference station (user2), and its geodetic position must be pre-known. Hereafter, the rover corresponds to user1, and the lander corresponds to user2.
In the MDPO method, multiple double-differenced pseudorange observations are obtained from multiple epochs, i.e., t i = t k , . . . , t k+N −1 , where N is the number of observation epochs, and k is the epoch number at which the estimation starts. It is important to note that the rover position must be fixed in place during all multi-epoch observations taken, in order to keep the number of estimation parameters lower than the number of observation equations. Otherwise, the rover position cannot be identified deterministically by the MDPO. We will estimate X X X R (t k ) = (x R (t k ), y R (t k ), z R (t k )) which represents a fixed rover position during multi-epoch The rover position can be estimated by solving the following Newton-Raphson equations iteratively. The following equations correspond to 'two dimensional(2D) MDPO' in [11], which calculates an estimated two-dimensional (X-Y) position by the Newton-Raphson equation whereas a rover altitude (Z) is calculated by using a lunar digital elevation model (DEM) as we discuss later.
In the 2D MDPO, the number of multi-epoch observations can be reduced to as low as 2 (N = 2). First, we define a new parameter R for t i = t k , . . . , t k+N −1 : where ∇ρ is the measured double-differenced pseudorange value and ∇ r 0 is an estimated double-differenced range based on an initial rover position estimate, i.e., X X X 0 R (t k ) (the superscript represents the number of iteration.). Then, the following equations can be derived: where G G G in (9) is called an observation matrix, which is equivalent to the Jacobian of R R R with regard to X X X . By solving the least-square problem that minimizes the residual error |R R R − G G GdX X X |, an estimated value of dX X X , defined as dX X X , is obtained: Then, a new estimated value X X X 1 R is given by (15), which provides a better fit to the observation.
This estimation process continues until the number of iterations reaches the designed value n, i.e., X X X 1 R ,X X X 2 R ,. . . ,X X X n R , and then the final estimated value X X X n R (t k ) is acquired. An expression for the quality of the estimates can be written asd Suppose w w w follows a white Gaussian distribution that has a mean value of zero and covariance matrix C C C: the covariance ofd X X X − dX X X , defined as P, is given by: Assuming the components of w w w are uncorrelated and have an identical variance, i.e., C C C = (σ ∇ω ) 2 I I I , the expression can be simplified as where σ ∇ωr is the standard deviation of the doubledifferenced receiver noises. In GNSS terminology, (G G G T G G G) −1 is known as the dilution of precision (DOP) matrix, which is used to specify error propagation as a mathematical effect of the navigation satellite geometry on positional measurement precision. We define the DOP matrix as where σ DOP is the elements of DOP. Using DOP, the root mean square (rms) of the user position accuracy at a time of t k can be obtained by where UPE represents user position error, which is the distance between the true rover position and an estimated rover position. VOLUME 10, 2022 It is important to highlight that the standard deviation of MDPO's double-differenced receiver noises, i.e., σ ∇ωr , is amplified from the standard deviation of the original receiver noises, i.e., σ ωr , as a result of the doubledifferencing process, in particular ω r (7), and becomes as large as σ ∇ωr = σ ωr 2 + σ ωr 2 + σ ωr 2 + σ ωr 2 = 2σ ωr assuming that the receiver noises follow a white Gaussian distribution.
Further, by defining GDOP as (20) can be written as (22) This can be decomposed into two directions as where UPE x and UPE y represents a user position error in X and Y direction respectively. In the 2D MDPO, the rover altitude z R is estimated by using a lunar digital elevation model (DEM): Here, z RDEM represents a lunar DEM model that is a function of X and Y positions. The estimation proceeds in the following sequence: First, X X X 0 R (t k ) is estimated using the rover position before its relocation, i.e., X X X 0 R (t k ) = X X X R (t k−1 ). Then, a new estimated rover position, i.e., X X X 1 (8)- (15): z R is not updated at this moment. After that, the altitude of the rover is updated to ) is acquired. The calculation continues until the number of iterations reaches the designed value, i.e., n.
According to (25), as z R changes along with x R and y R , errors in the X-Y position induce errors in the Z position, which ultimately induce errors in estimated x R and y R in turn. As a result, UPE deteriorates stochastically. In our research, we did not apply the case in which the rover altitude changes too rapidly, such as the rover dropping off the cliff or roving on steep slopes. In that case, the output of z RDEM does not change with x R and y R too significantly, and UPE would not deteriorate, which was confirmed by numerical simulations in our prior research [11].
It is important to note that GNSS measurements are also subject to systematic errors, such as the satellite orbit determination error, time tag error, and DEM information error as modeled in [11]. Systematic errors are implemented in the numerical simulation.

B. INS/GNSS INTEGRATION
This section provides a review of the INS/GNSS integration form that we use in our research. We fuse the selected INS/GNSS integration form with the mathematical model of the MDPO-based GNSS measurement in the next section.
In the INS/GNSS integrations forms, there are the three most common integration strategies: the looselycoupled [32], the tightly-coupled [32]- [34] and the ultra-tight integration [35], [36]. Since the ultra-tight integration involves the baseband signal processing of GNSS receivers (i.e., the digital tracking loops), that is typically not accessible using commercial products. The basic difference between tight-integration and loose-integration is the type of data shared by the GNSS receiver. In the loosely-coupled technique, the positions and velocities estimated by the GNSS receiver are blended with the INS navigation solution, whereas in the case of tightly-coupled method, GNSS raw measurements (i.e., pseudorange and Doppler observables) are processed through a unique Kalman filter with the measurements. Whereas tightly-coupled integration works more robustly with however many navigation satellites available, loosely-coupled integration is simpler in implementation and still beneficial in case the computing power of spacecraft (i.e., rover) is limited. In our study, we formulate the loose INS/GNSS integration in combination with the MDPO algorithm.
In the loose INS/GNSS integration format, positions derived by GNSS measurement are merged as updates of the INS estimates through a Kalman filter. The Kalman filter model assumes that the true state at time k is evolved from the state at k − 1 according to where x x x k denotes the state, u u u k denotes the control input, ω ω ω k denotes the process noise, F F F k denotes the state-transition model, B B B k denotes the state-transition model, G G G k denotes the noise-transition model. ω ω ω k is assumed to be drawn from a zero mean multivariate normal distribution N with covariance, i.e., Q Q Q k : ω ω ω k ∼ N (0, Q Q Q k ). At time k, an observation z k of the true state x x x k is made according to: where v v v k denotes the observation noise, H H H k is the observation model, which maps the true state space into the observed space. v v v k is the process noise and assumed to be drawn from a zero mean multivariate normal distribution N with covariance, i.e., To solve Kalman filter, the initial state, and the noise vectors at each step, i.e., x 0 , w 1 , . . . , w k , v 1 , . . . , v k , are all assumed to be mutually independent. The state of the filter is represented by two variables:x x x k|k which is the a posteriori state estimate at time k given observations up to and including at time k, and P P P k|k which is the a posteriori estimate covariance matrix. In Kalman filter, typically, the two phases alternate, with the prediction advancing the state until the next scheduled observation, and the update incorporating the observation in the following sequence: Predictx Update e e e k = z z z k − H H H kx x x k|k−1 (30)

S S S k = R R R k + H H H k P P P k|k−1 H H H
x x x k|k =x x x k|k−1 + K K K k e e e k (33) P P P k|k = (I I I − K K K k H H H k )P P P k|k−1 (34) where e e e k is the innovation pre-fit residual, S S S k is the innovation covariance, K K K k is the optimal Kalman gain.

C. INS/GNSS INTEGRATION DESIGN FOR THE MDPO ALGORITHM
In our study, we assume that the system is equipped with sensors that measure the rover's orientation and traveling distance as shown in Fig. 2: the states x x x k consists of a twodimensional rover's position relative to the reference station in a local topocentric coordinate (East-North-Up), i.e., x k and y k , and an orientation θ k . The control input u u u k consists of a traveling distance d k and reorientation angle θ k . The control input u u u k is known to the system through INS measurements but in reality u u u k has errors, which refers to ω ω ω k in the equations and are not known to the system. According to the 2D MDPO equations (16), GNSS provides two-dimensional user position observations, i.e., x ob and y ob with estimation errors in the form of (G G G T G G G) −1 G G G T w w w whose covariance are represented by XDOP(t k ) × 2σ ωr and YDOP(t k ) × 2σ ωr respectively. We assume that INS provides measurements of the rover's orientation, i.e., θ ob , with estimation errors of ω θ . These measurements are summarized in observation data matrix z z z k . To summarize, each matrix was constructed as where σ d represents the standard deviation of the traveling distance measurement, σ θ represents the standard deviation of the reorientation angle measurement, and σ θ represents the standard deviation of the orientation measurement.

D. SYSTEMATIC ERROR
Actual INS measurements contain systematic errors such as temperature coefficient error, bias, random walk, alignment error depending on their characteristics. Those errors can be effectively eliminated through offline calibration process [32] or online estimation [37]. However, some unremoved component of systematic errors can remain. Therefore, in this study, the following bias component were added to measurements data to evaluate their influences in the analysis: where δd is the systematic error of traveling distance measurement, δ θ is the systematic error of reorientation angle VOLUME 10, 2022 measurement, δθ ob is the systematic error of orientation measurement, and δx ob and δy ob are the systematic errors of GNSS measurements. It is important to note that the systematic errors of GNSS measurements are already implemented in our LGNS simulation [11], and not implemented additionally for this research.

E. UNCERTAINITY OF COVARIANCE MATRICES
When inferring parameters from a Gaussian-distributed data set by computing a likelihood, covariance matrices are needed that describes the data errors and their correlations. In our case, that refers to Q Q Q k and R R R k where Q Q Q k is the covariance matrix of the traveling distance and reorientation measurements and R R R k is the covariance matrix of the GNSS and orientation measurements. If the covariance matrices are not known perfectly in advance, estimator needs to rely on information known at that moment which may contain errors from the true covariance matrices. There are several techniques to infer the true covariance matrices out of measurements [38]. However, in some combinations of sensors, the true covariance matrices of measurements are not easy to estimate a priori or in real time (we don't specify what sensors to use in this study for generalization). Therefore, in this simulation, we added some intrinsic uncertainty to the true covariance matrices to create the estimator covariance matrices:

F. OTHER CONSIDERATION
For real-time implementation, the system needs to cope with the synchronization between measurements, in particular GNSS data latency for this application. Each time the rover uses the GNSS receiver, the rover gets observable delayed with respect to 1 pulse-per-second (PPS) signal. Furthermore, additional delay referring to the double-differenced calculation must be processed. In order to manage the GNSS data latency, INS measurements are buffered and processed only at proper timing.

A. SIMULATED SCENARIO
A scenario was developed to evaluate how the total traveling distance can be extended by the proposed INS/GNSS integration. In the developed scenario, the rover continues to relocate regardless of the visibility of the navigation satellites, whereas the rover moves only when both of two satellites are in view in the original scenario of the dual-satellite LGNS simulation analyzed in [12] (we call the original scenario 'GNSS-only' in this article). Each time the rover moves, the rover obtains INS measurements. The rover obtains GNSS measurements additionally when the two satellites were in view.

B. SIMULATION SETUP
The rover trajectory and GNSS measurements in this work were created with the same simulation parameters as our prior work [12]. Therefore, this article can help provide the comparison between the proposed INS/GNSS integration and the original raw dual-satellite LGNS measurements.
In the numerical simulation, the rover relocation and position estimation were repeated in turn: (step1) in the first epoch, for instance t = t k , the rover relocates its position and orientation, i.e., x x x k , according (46), (step2) in the next two epochs, for instance t = t k+1 and t = t k+2 , the rover obtains observation data z z z k according to (47) and applies the proposed INS/GNSS integration to update the estimated rover trajectoryx x x k according (28)- (34). Then, step1 and step2 continue over the course of the simulated mission period.
The rover trajectory, i.e., x x x k , was created dynamically by changing the rover position and orientation every three epochs according (46). The traveling distance, i.e., d, reorientation angle θ , process noise, i.e., σ d and σ θ , and systematic error, i.e., δd and δθ ob , are specified in Table. 1. The total simulated mission period was set 15,000 minutes and one epoch corresponded to 30 seconds. Accordingly, the total simulated mission period was equivalent to 30,000 epochs.
GNSS measurements, i.e., x ob and y ob , were provided to the generated rover trajectory using the LGNS simulation that was developed in our prior work. The rover obtained an updated GNSS measurement when the two navigation satellites are in view. With respect to the observation noise and systematic errors, the same parameters as Section 4.2 of [12] were used. In particular, the receiver noise σ ωr was set 0.2 m. Systematic error, i.e., δx ob and δy ob , such as satellite orbit determination error, time tag errors, DEM error were already included in the generated data.
These sensor error parameters, σ d , σ θ , σ θ , δd, δθ ob , and δθ ob , assume micro-sized rovers and were created by adding safety margins to the values reported in other micro-rover navigation studies, namely [29] for the traveling distance measurement error and [20] for the reorientation and orientation measurement error. It is important to note that these values are conservative and can be improved for larger rovers.
Finally, position estimates were provided according to (28)- (34) and compared with the true rover trajectory. As described in Section II.E, we assumed that the covariance matrices of measurements were not known perfectly to the estimator. We thus added intrinsic uncertainties: δσ d is 10 % of σ d , δσ θ is 10 % of σ θ , δσ θ is 10 % of σ θ , and δσ ωr is 20 % of σ ωr assuming that the GNSS measurements had more uncertainity than the INS measurements assuming that the performance of INS can be validated and calibrated on the ground in advance of the actual missions.  To the purpose of a comparative study, user trajectory estimated only with INS sensor was also created, so-called INS-only navigation scenario as shown in Fig. 3. For that case, Kalman filter equation were modified to exclude GNSS measurements: C. NUMERICAL SIMULATION RESULT Fig. 4 shows the true and estimated rover trajectories. Fig. 5 shows the user position error projected onto the local topocentric coordinate. Fig. 6 shows the time sequence of the user position error. In Fig. 4, GNSS measurements were observed in a locally clustered form. This is due to the fact that GNSS measurements were not available all the time and obtained in limited locations corresponding to the visibility of the two navigation satellites. Also, a Monte Carlo simulation was conducted 100 times, and averaged data were presented for each navigation algorithms. Table. 2 shows the traveling distance and twice the distance root mean square (2drms) of the user position error. The traveling distance of the INS-only and INS-GNSS integration were larger than the GNSS-only scenario. This is because in the INS-only and INS-GNSS integration scenarios, the rover continued to travel regardless of the visibility of the navigation satellites, whereas in the GNSS-only scenario the rover only traveled when both navigation satellites were in view.

IV. DISCUSSIONS
As seen in Fig. 6, it was confirmed that the user trajectory estimated by the INS-only navigation algorithm had offsets from the true trajectory, which increased over time due to the accumulated INS measurement errors. On the other hand, GNSS measurements had scattered errors rather than offset errors, which were derived from normally-distributed receiver noises at each observation. The characteristics of the two errors were complementary, and both errors were corrected and largely removed through Kalman filter in the proposed INS/GNSS integration.
Consequently, the proposed INS/GNSS integration successfully improved the positioning accuracy compared to our previous scenario [12], i.e., GNSS-only in this simulation. Under the simulated condition, we confirmed that the proposed INS/GNSS integration achieved beyond 25 m (2drms) user position accuracy, which is a significant improvement from the user position accuracy of the GNSS-only scenario that was about 60 m (2drms). Furthermore, the proposed INS/GNSS integration extended the rover's traveling distance almost 10 times longer by compensating for the vehicle position in the absence of GNSS measurements.
One limitation of the proposed method is derived from the availability of INS measurements. In shadowed regions such as permanently shadowed regions (PSRs) at poles, visual odometry is constrained, if not, invalidated. Under those circumstances, the rover needs to rely on GNSS measurements rather than INS measurements to maintain the position accuracy. In such a case, the rove may needs to lengthen the observation period of the MDPO to improve the accuracy of the GNSS measurements.
Another limitation may be imposed by the availability of the GNSS measurements. For instance, the orbital parameters of the two navigation satellites changes over time due to the lunar gravity perturbations. While navigation satellite commonly owns a propulsion system to maintain its orbital parameters, the rover may face a situation where the rover receives only one GNSS signal at a time instead of two satellites. In that case, tightly-coupled integration is a viable VOLUME 10, 2022  approach that can use single-satellite pseudorange observations in position estimation process. We will study the tightly-coupled integration in our future work.

V. CONCLUSION
In this study, we proposed an INS/GNSS integration model adopted for Multi-Epoch Double-Differenced Pseudorange (MDPO)-based dual-satellite lunar global navigation satellite systems. We also evaluated its performance in a simulated environment.
In the scenario where a rover travels on the lunar surface using GNSS measurements from two navigation satellites, we used the proposed INS/GNSS integration to 1) estimate rover positions when GNSS measurements are not available, and as a result, extend rover's traveling distance, 2) improve the user position accuracy by fusing two different sensor characteristics. Under the selected numerical simulation condition, it was confirmed that the proposed INS/GNSS algorithm successfully extended the total traveling distance by compensating for navigational measurements in the absence of GNSS measurements while it also improved the user position estimation accuracy beyond the accuracy of the raw GNSS measurements.
Our future work includes case studies with different LGNS parameters. In particular, the trajectory of the two navigation satellites is a key trade-off design parameter which affects the accuracy and availability of the GNSS measurements. To obtain further insights for the future dual-satellite LGNS design, other simulation cases will be studied.