## Abstract

There is a growing interest in the use of legacy terrestrial Global Positioning System (GPS) signals to determine the precise positioning and timing onboard a lunar satellite. Unlike prior works that utilize pseudoranges with meter-level accuracy, we propose a precise positioning and timekeeping technique that leverages carrier-phase measurements with millimeter-level accuracy (when integer ambiguities are correctly fixed). We design an extended Kalman filter framework that harnesses the intermittently available terrestrial GPS time-differenced carrier-phase (TDCP) values and gravitational accelerations predicted by the orbital filter. To estimate the process noise covariance, we implement an adaptive state noise compensation algorithm that adapts to the challenging lunar environment with weak gravity and strong third-body perturbations. Additionally, we perform measurement residual analysis to discard TDCP measurements corrupted by cycle slips and increased measurement noise. We present Monte-Carlo simulations of a lunar satellite in an elliptical lunar frozen orbit and quasi-frozen low lunar orbit, wherein we showcase higher positioning and timing accuracy as compared with the pseudorange-only navigation solution.

## 1 INTRODUCTION

After more than 50 years since the Apollo program, there is a growing interest in establishing a sustainable human presence on the Moon. To support increasing plans for lunar exploration, various missions are being planned in different lunar orbit regimes, e.g., the European Space Agency’s lunar Pathfinder (Giordano et al., 2022), which is to be deployed in an elliptical lunar frozen orbit (ELFO), and the Lunar Gateway (Winternitz et al., 2019), which is being spearheaded by the National Aeronautics and Space Administration (NASA), in a near-rectilinear halo orbit (NRHO). These missions require precise positioning and onboard timing to execute reliable communication; onboard guidance, navigation, and control; and mission activity planning. As per NASA’s Services Requirement Document (SRD) for lunar relay satellites (National Aeronautics and Space Administration, 2022), the required signal-in-space error (SISE) in position is 13.34 m (3*σ*, 99.7% of the time), and the SISE requirement in velocity is 1.2 mm/s (3*σ*) over any 10-s interval.

### 1.1 Prior Work on Positioning, Navigation, and Timing in Lunar Orbits

In the Lunar Reconnaissance Orbiter (LRO) mission, positioning accuracy better than 20 m is achieved by combining range and Doppler radiometric data tracked from the Universal Space Network and altimetric data from the lunar orbiter laser altimeter (Mazarico et al., 2012). However, in ground-based tracking, the position, navigation, and timing (PNT) information is estimated with a delay due to the time required for data post-processing and uplinking, thus restricting the real-time estimation of lunar satellite position. In addition, the scalability of these ground-station-based techniques is limited because they require dedicated deep-space antennas for recurring downlink/uplink tasks. Moreover, conducting orbit determination via ground station measurements requires human-in-the-loop operations, which could lead to an increase in mission operation costs. To this end, efforts have been made to achieve onboard navigation for lunar orbits. For example, the cislunar autonomous positioning system technology operations and navigation experiment (CAPSTONE) satellite, which entered an NRHO around the Moon, estimates its position in space without having to rely exclusively on Earth tracking stations (Cheetham et al., 2022). In particular, the CAPSTONE satellite communicates directly with the LRO and utilizes crosslink data to measure its distance and rate of distance change from the LRO, thus determining CAPSTONE’s position.

There has been growing interest in harnessing the signals being broadcast by the legacy terrestrial Global Positioning System (GPS) to achieve precise, standalone positioning in a lunar orbit. This research is motivated by the widespread use of terrestrial GPS for near-Earth space applications. Prior simulation works (Capuano et al., 2016, 2017; Schonfeldt et al., 2020) have demonstrated the presence of continual time segments at lunar distances during which terrestrial GPS signals are not occluded by Earth or the Moon, and their received signal power is sufficient for signal acquisition and tracking. These findings establish that, in a lunar orbit, it is possible to intermittently track a sufficient number of terrestrial GPS satellite signals required for estimating a full navigation solution (i.e., position, velocity, and timing). SpacePNT has developed a high-sensitivity, spaceborne global navigation satellite system (GNSS) receiver named NAVIMOON (Pultarova, 2021; Scotti et al., 2022), which will fly onboard the Lunar Pathfinder in 2026 and will perform a satellite navigation positioning fix in a lunar orbit. Additionally, the NASA Goddard Space Flight Center (GSFC) has pioneered work in developing the Navigator, a fully space-flight-qualified receiver with the ability to track very weak terrestrial GPS signals (Winternitz et al., 2004). GSFC is planning to test GPS navigation in the Earth-Moon transfer orbit in 2024 during the Lunar GNSS Receiver Experiment (Parker et al., 2022).

Particularly, achieving precise positioning and timing of satellites in lunar orbits poses significant challenges, including the following key issues:

Increased measurement noise due to the low carrier-to-noise ratio (

*C/N*_{0}) of weak terrestrial GPS signals;Weak central gravitational potential and significant gravitational anomaly of the Moon, stronger third-body effects due to Earth, and relatively stronger solar radiation perturbations compared with central gravity, especially in highly eccentric orbits such as ELFOs and NRHOs;

Intermittent outage of terrestrial GPS signals in low lunar orbits and near perilunes of eccentric orbits; and

Poor geometry of terrestrial GPS satellites as seen from lunar satellites, leading to large dilution of precision (DOP) values and low observability of terrestrial GPS measurements.

Several works have also evaluated terrestrial GPS-based positioning and timing performance in various lunar orbits, including Earth-Moon transfer orbits, ELFOs, NRHOs, and low lunar orbits (Capuano et al., 2016, 2017; Lopes et al., 2014; Small et al., 2022; Winternitz et al., 2019). Capuano et al. (2016) and Capuano et al. (2017) developed an orbital filter that fuses observations of terrestrial GPS pseudorange and pseudorange rate with an orbital force model through an adaptive extended Kalman filter framework. In our prior works (Bhamidipati et al., 2022a, 2022b), we successfully demonstrated the use of terrestrial GPS pseudoranges to achieve precise and standalone timing onboard small satellites (SmallSats) in various lunar orbits such as low lunar orbit, prograde circular orbit, ELFO, and NRHO. In another prior work (Iiyama et al., 2021), we developed a decentralized, autonomous state estimation algorithm for SmallSats that provides PNT service using terrestrial GPS pseudoranges and inter-satellite communication links. However, all of these works rely on pseudorange measurements with only meter-level accuracy, thus restricting the positioning and timing accuracy achieved in lunar orbits.

### 1.2 Opportunities and Challenges of Using Terrestrial GPS Carrier-Phase Measurements in a Lunar Orbit

To achieve precise navigation in a lunar orbit, it would be beneficial to explore leveraging carrier-phase measurements that provide millimeter-level ranging accuracy when integer ambiguities are correctly fixed, unlike the meter-level pseudorange measurements used in prior lunar PNT works.

While carrier-phase measurements could provide more precision compared with pseudorange, this precision comes at the cost of fixing the unknown integer ambiguity. Integer ambiguity is often resolved via a least-squares ambiguity decorrelation adjustment (Teunissen et al., 1997). For pairs of satellites in an Earth orbit, several missions have demonstrated that precise real-time relative navigation can be achieved by processing single-difference carrier-phase measurements (D’Amico et al., 2013; Kahr et al., 2018). However, it is difficult to resolve integer ambiguity for the absolute positioning of lunar satellites because of the following challenges:

Currently, base stations or reference satellites with well-known positions that could be used for differential positioning do not exist on the Moon or in a lunar orbit.

It is difficult to fix integer ambiguities when the terrestrial GPS geometry seen from the user is poor compared with Earth orbits.

It is challenging to keep track of the carrier-phase signal from the same terrestrial GPS satellite for a long time in a lunar orbit, because of signal blockage by the Earth and Moon.

The time-differenced carrier-phase (TDCP) formulation avoids the problem of integer ambiguity fixing by differentiating the carrier-phase measurement between two time steps (Guillard et al., 2022; Kim et al., 2019; Zhao, 2017). The TDCP appraoch utilizes the fact that the integer ambiguity is constant over time while the carrier phase is tracked, which is true except when a cycle slip occurs. The TDCP formulation is also beneficial for removing slow-varying signal errors caused by terrestrial GPS satellite clocks and relativistic effects. Moreover, the TDCP formulation is well suited for the rapid motion of a lunar satellite, as we can benefit from the geometrical diversity of the terrestrial GPS satellites between two carrier-phase values at different time steps.

### 1.3 Proposed Framework

We propose a precise positioning and timing technique (see Figure 1) for lunar satellites that harnesses the intermittently available terrestrial GPS TDCP measurements and the gravitational accelerations predicted by the orbital filter. Note that our proposed framework can improve state estimates only when terrestrial GPS satellites are available. When terrestrial GPS measurements are not visible because of occultations from Earth and the Moon, the orbital filter can mitigate the drift in the position estimate (see results in Section 4 for further details).

Here, we design a TDCP-based extended Kalman filter framework. To account for time correlations across TDCP measurements, we design an augmented state vector that incorporates consecutive lunar satellite states and performs fixed-point smoothing at each measurement update. The increase in the computational cost of the time and measurement update due to the increase in state size is minimized by effectively skipping unnecessary computations in the classical Kalman filter. The filter also executes cycle-slip detection, which discards corrupted TDCP measurements via measurement residual analysis. Moreover, we also implement an adaptive state noise compensation (ASNC) algorithm (Stacey & D’Amico, 2021) to adaptively adjust the process noise of the filter and to enhance the robustness of estimation in the challenging, dynamic lunar environment (i.e., with weak gravity and strong external perturbations). This approach is particularly beneficial for elliptical lunar orbits such as ELFOs, wherein the dynamic environment changes significantly between periapsis and apoapsis, thus also changing the process noise statistics. Our proposed framework will bridge the gap in lowering the current standards for positioning and timing accuracy achievable in a lunar orbit and thereby facilitate reliable operations during the entire mission life cycle of a lunar satellite. To the best of the authors’ knowledge, this is the first work to leverage TDCP measurements of terrestrial GPS satellites to navigate lunar satellites.

### 1.4 Key Contributions

Our key contributions are listed below. This work is based on our prior conference paper presented at the 2023 International Technical Meeting of the Institute of Navigation (Iiyama et al., 2023a).

We propose a precise PNT technique for lunar satellites that harnesses the intermittently available terrestrial GPS TDCP values and gravitational accelerations predicted by the orbital filter.

We address the problem of conditional dependence that occurs in the processing of GPS TDCP measurements by introducing augmented states constructed by concatenating states of two time steps.

We propose a time update and measurement update scheme that effectively reduces the computational cost compared with directly applying the Kalman filter update equations to the augmented state.

We detect outliers to discard TDCP measurements corrupted by increased measurement noise due to a low received (

*C/N*_{0}) and cycle slips.We execute an ASNC algorithm to optimally adjust the process noise of the filter based on the dynamic characteristics of a lunar satellite within orbit.

Through high-fidelity Monte-Carlo simulations of a lunar satellite in an ELFO and quasi-frozen low lunar orbit (QFLLO), we confirm that our proposed framework achieves improved positioning and onboard timing as compared with the pseudorange-only navigation technique.

The remainder of this paper is organized as follows. Section 2 describes the formulation of lunar satellite dynamics and terrestrial GPS TDCP measurements. Section 3 discusses our proposed TDCP-based extended Kalman filter framework. Section 4 explains our simulation setup and experimental results. Section 5 provides concluding remarks.

## 2 PRELIMINARIES

We now describe the preliminaries regarding the linearized dynamics associated with the nonlinear motion of lunar satellites. We also introduce the math related to terrestrial GPS-based TDCP measurement formulation.

### 2.1 Lunar Satellite Dynamics Formulation

The satellite dynamics in a lunar orbit can be explained in terms of the equations of motion associated with their a) position and velocity via orbital dynamics, and b) time via dynamics of the onboard clock.

#### 2.1.1 Orbital Dynamics

Based on the work by Montenbruck and Gill (2000), at any time step *k*, the satellite dynamics is a nonlinear function of time, position, velocity, and the solar radiation pressure coefficient and is given by the following:

1

where *γ* denotes the ballistic coefficient of the solar radiation pressure, *t _{k}* denotes the time based on GPS time, and denote the position and velocity of the lunar satellite in the Moon-centered inertial (MCI) frame, respectively, and

*ϵ*denotes the associated process noise. In this paper, we use a Moon-centered J2000 frame for the MCI frame. Additionally,

^{rv}**(**

*a**t*,

*r*_{u},

*v*_{u},

*γ*) is the total acceleration on the lunar satellite, which equals the sum of the gravitational acceleration

*a*_{g}(

*t*,

*r*_{u},

*v*_{u}) and the acceleration due to solar radiation pressure

*a*_{⊙}(

*γ*).

In Equation (1), the gravitational acceleration of the lunar satellite, which is calculated as the sum of the high-order gravitational terms from the Moon and third-body perturbations from the Earth, Sun, and Jupiter, can be expressed as follows (Montenbruck & Gill, 2000):

2

3

where *U _{g}* and ∇

*U*denote the gravitational potential and its gradient, respectively,

_{g}*G*is the universal gravitational constant, which equals 6.67430 × 10

^{−11}m

^{3}/kg s

^{2},

*M*

_{☾},

*R*

_{☾}are the mass and position of the Moon,

*r*^{p}is the position of the third-body planets (Earth, Sun, and Jupiter), and

*M*is the corresponding mass. At each time step,

_{p}

*r*^{p}is obtained from the DE440 ephemerides data available in the Spacecraft, Planet, Instrument, C-matrix, Events (SPICE) library of the Navigation and Ancillary Information Facility (NAIF) (Acton et al., 2018). Additionally,

*ϕ*and

*λ*are the latitude and longitude of the spacecraft position in the Moon-fixed frame, respectively. is the normalized Legendre function, and and are the normalized geopotential coefficients. For the Moon-fixed frame, we used the “IAU_MOON” frame provided in the NAIF SPICE library.

*N*is the maximum degree and order of the higher-order gravity terms.

_{sph}In Equation (1), the acceleration due to solar radiation pressure *α*_{⊙} is calculated from the cannon ball model as follows:

4

where Φ_{⊙} denotes the solar flux at 1 astronomical unit (AU) (=1360 W/m^{2}; 1 AU = 1.498 × 108 km) and *A ^{u}* and

*M*denote the surface area and mass of the lunar satellite, respectively. Additionally,

^{u}*d*denotes the heliocentric distance to the lunar satellite,

*ê*_{Sun}denotes the sun unit direction vector from the lunar satellite, and

*C*= 1 +

_{R}*ϵ*denotes the coefficient of reflection, where 0 ≤

*ϵ*≤ 1 is the radiation reflection ratio of the spacecraft body. Finally,

*c*denotes the speed of light, which equals 299792.458 km/s. To account for orbital perturbations due to solar radiation pressure with high precision, existing literature on orbit determination of satellites (Montenbruck & Gill, 2000) estimates as part of the state.

Let ** S**(

*t*) = [

*r*^{u}(

*t*)

*v*^{u}(

*t*)

*γ*]

^{⊤}be the vector of states related to the orbital dynamics. To obtain the linearized dynamics, consider the linearized state

*s*around the reference state as follows:

_{k}5

where the reference state is initialized at the current estimate of the states, (e.g., the predicted state estimated via a time update step, as explained in Section 3.2).

By re-defining ** a**(

*t*,

*r*^{u},

*v*^{u},

*γ*) from Equation (1) as

**(**

*a**t*,

**(**

*S**t*)), we obtain the linearized dynamics as follows:

6

where:

7

By propagating Equation (6) from *t*_{k–1} to *t _{k}*, we have the following:

8

where is the state transition matrix. The state transition matrix can be computed by propagating the following equation from *t*_{k–1} to *t _{k}*:

9

#### 2.1.2 Clock Dynamics

The clock bias and clock drift are propagated by using the following discrete, two-state error model taken from Galleani (2008):

10

11

where the noise vector is Gaussian noise sampled from the random-walk process of the clock.

### 2.2 Terrestrial GPS TDCP Formulation

It is difficult to resolve the integer ambiguity in a non-differential mode and poor GPS geometry scenario, as considered in this paper. However, as discussed in Section 1, under the condition that the integer ambiguity is constant over time, the integer ambiguity can be eliminated by differencing the carrier-phase measurements between two time steps.

The terrestrial GPS TDCP measurement of any *i*-th satellite Δ* _{T}ϕ^{i}* using carrier-phase values at time steps

*k*and (

*k*– 1), i.e., and , respectively, is given by the following:

12

where, eliminating subscripts *k* and *k* – 1 for the sake of simplicity, *τ ^{u}* is the receiver’s clock bias and is the carrier-phase measurement noise. For any terrestrial GPS satellite

*i*,

*r*^{i}is the satellite position in the MCI frame, (

*e*^{u})

^{i}is the unit line-of-sight vector from the lunar satellite to the terrestrial GPS satellite

*i*, and

*τ*is the satellite clock error. Additionally, Δ

^{i}_{T}is sufficiently small such that the time difference in the tropospheric, ionospheric, and multipath-related terms in the standard carrier-phase equation (Kaplan & Hegarty, 2017) is considered negligible and is thus absorbed into the residual noise term . Additionally, we utilize terrestrial GPS-based pseudorange and pseudorange rate measurements in our proposed filter, with standard equations taken from the work of Kaplan and Hegarty (2017). A detailed explanation regarding the terrestrial GPS measurement formulation can be found in our prior work (Iiyama et al., 2023b).

## 3 PROPOSED TDCP EXTENDED KALMAN FILTER

Figure 2 shows a high-level overview of our proposed framework, which is based on an extended Kalman filter. In this framework, we execute the time update step based on the orbital dynamics and the measurement update step via terrestrial GPS-based pseudorange, pseudorange rate, and TDCP measurements. We perform outlier detection in carrier-phase measurements during the measurement update step. We implement an ASNC algorithm to adjust the process noise to be used in the time update. The details of each process are provided in the following subsections.

### 3.1 Augmented State

Let be the lunar satellite state vector defined in Equation (13):

13

The lunar satellite state vector includes the position *r*^{u}, velocity *v*^{u}, solar radiation pressure reflection coefficient *C _{R}*, clock bias

*τ*, and clock drift . The position

^{u}

*r*_{u}and velocity

*v*_{u}are in the MCI frame, where the axes of the inertial frame are chosen to coincide with J2000.

In Equation (12), and are functions of and , respectively. Therefore, the TDCP measurement is a function of two consecutive lunar satellite states. This breaks one of the assumptions of the Kalman filter, which is the conditional independence of and , given , as follows:

14

Previous works have addressed this issue by modifying the extended Kalman filter update equations to account for the conditional dependance (Guillard et al., 2022; Kim et al., 2019). However, the framework in those previous works processes the TDCP measurement separately from the other measurements (e.g., pseudorange), which ignores the effect of future measurements on the past state estimate. Note that when TDCP measurements are present, future measurements will also provide information on the past states because the TDCP measurements will correlate the states at two different time steps.

In this paper, we propose an alternative approach, where we construct an augmented state:

15

by concatenating the two states at consecutive time steps at which the two carrier-phase measurements that construct the TDCP measurement are taken. By introducing this new state, the TDCP measurement becomes a function of *X _{k}* only; thus, conditional independence can be secured. One advantage of the augmented state approach is that it accounts for the aforementioned effect in the measurement update by processing other measurements together with the TDCP measurements. Another advantage of the augmented state approach is that it can handle an arbitrary number of intermediate measurements (e.g., inertial measurement unit [IMU] data, inter-satellite range) between two TDCP measurements. While intermediate measurements are not considered in this paper, we have presented a framework to jointly process IMU and TDCP measurements with the augmented state approach in our previous work (Iiyama et al., 2023b).

### 3.2 Time Update Step

We first define the linearized state vector around the reference lunar satellite state (i.e., current estimate of the state in the extended Kalman filter) vector as follows:

16

Based on Equations (16), (8), and (11), the linearized dynamics of the lunar satellite state vector at any *k* time step are given by the following:

17

18

19

where *ϵ _{k}* denotes zero-mean Gaussian process noise, whose covariance depends on the noise covariance in the lunar satellite position and velocity, which is given by , and that in the lunar satellite clock, which is given by . Note that the state transition matrices Φ

*(*

^{rvc}*t*,

_{k}*t*

_{k–1}) and Φ

^{clk}(

*t*,

_{k}*t*

_{k–1}) were defined in Equations (9) and (11), respectively. We set the process noise of the solar ballistic coefficient

*γ*to 0 because its variation will most likely be slow enough to justify the assumption of zero process noise.

In Equation (19), the process noise covariance is modeled as , where *σ*_{1}, *σ*_{2} are the diffusion components of the clock phase deviation and clock frequency deviation.

The process noise covariance is formulated as follows (Tapley et al., 2004):

20

21

where is defined as the unmodeled acceleration, is the position and velocity component of the state transition matrix, and *κ* is the variable of integration, which represents time. With the assumption that the process noise does not perturb the state enough to change the accelerations over the update interval, Equation (20) can be approximated as follows (Myers, 1974):

22

Next, we formulate the linearized augmented state vector *x _{k}* given by , where was introduced in Equation (16). In the time update of our proposed framework, we propagate the linearized augmented state

*x*and its covariance

_{k}*P*as follows:

_{k}23

24

where and are the predicted augmented state vector and its covariance at time step *k*, respectively. Similarly, and are the estimated augmented state vector and its covariance at time step (*k* – 1), respectively. The state transition matrix of the augmented state Φ(*t _{k}*,

*t*

_{k–1}) is given by with the state transition matrix defined in Equation (18). In addition, the process noise covariance

*Q*is given by , where the process noise covariance matrix of the linearized state vector was defined in Equation (19). Note that the state transition matrix and process noise covariance matrix corresponding to the second part of the state are identity and zero matrices, respectively, because the second half of the state represents the fixed past state

_{k}*t*

_{k–1}and thus remains unchanged.

### 3.3 Measurement Update Step: Modified for Computational Efficiency

Compared with simply applying the measurement update step equations performed in the extended Kalman filter to the full augmented state, we can reduce the computational cost by only calculating half of the Kalman gain matrix. This approach is possible because the second half of the updated linearized augmented state vector (i.e., after the measurement update step) will be the estimate of the lunar satellite state vector at time *t*_{k–1} based on the measurements up to *t _{k}* (i.e., ). The full estimate and its covariance need not be calculated unless we want to refine our past state estimate at which the last terrestrial GPS measurement was received. Because the augmented state will be re-initialized after the measurement update, as shown in Section 3.4, and will not be used in the later time steps, the Kalman gain only has to be calculated for the first half of the states (as shown later in Equation (26)).

In the measurement update step, our proposed filter processes a stacked measurement vector that combines three types of terrestrial GPS measurements, i.e., pseudorange *ρ*, pseudorange rate , and TDCP Δ* _{T}ϕ*. Let

*z*be the stacked measurements obtained at time step

_{k}*k*as follows:

25

where *m* and *m*′ are the number of available terrestrial GPS signals and TDCP signals, respectively. The total number of measurements will be *M* = 2*m* + *m*′. We additionally model the measurement covariance matrix *R _{k}* as a time-dependent diagonal matrix. We modeled the pseudorange error based on the thermal noise of the receiver delay lock loop (DLL), carrier tracking errors with phase lock loop (PLL) thermal noise, and pseudorange rate errors with frequency lock loop (FLL) tracking loop jitter. More details regarding this covariance modeling can be found in the work by Kaplan and Hegarty (2017). We also added a Gaussian noise term of 10 m (2

*σ*) to the pseudorange and 0.01 m/s (2

*σ*) to the pseudorange rate as a SISE. This value is based on the SISE values (≤7.0 m 95% user range error and ≤0.006 m/s 95% user range rate error) specified in the GPS performance standard (Department of Defense, 2020). We added some margins to the specified values to account for additional errors, such as group delays in the sidelobe signal and ionospheric delays. Time correlations in the SISE are not modeled in this work, but will be addressed in future work. For the TDCP measurements, we assumed that SISEs and propagation errors were completely removed via time differentiation.

As explained earlier, to improve computational efficiency, we compute the Kalman gain for only the first half of the states. To do so, we first divide the predicted state covariance (the covariance matrix after the time update in Equation (24)) into four sub-matrices of size , as follows:

26

We then calculate the Kalman gain as follows:

27

28

29

Here, *H _{k}* denotes the linearized measurement model (defined in Equation (28)) that stacks the linearized measurement models associated with the pseudorange

*H*, pseudorange rate , and TDCP . The augmented state and its covariance will be reconstructed from these values following the procedure shown in Section 3.4. For any

_{ρ,k}*i*-th satellite, the linearized measurement models of the pseudorange , pseudorange rate , and TDCP are defined in Equations (30), (31), and (32), respectively:

30

31

32

With our measurement vector *z _{k}* from Equation (25), the modeled measurement covariance

*R*, and the updated Kalman gain

_{k}*K*from Equation (27), the state estimate and the covariance estimate are updated via the sub-component of the Joseph form equation for the extended Kalman filter (Bucy & Joseph, 1987) as follows:

_{k}33

34

35

where denotes the measurement residual, i.e., the difference between the received measurements *z _{k}* and the expected measurement for the predicted augmented state vector (defined as a stack of the measurement equations related to pseudorange, pseudorange rate, and TDCP). The Joseph form is mathematically equivalent to the standard Kalman filter measurement update equations, but has better numerical robustness and prevents the covariance matrix from becoming non-positive definite (Tapley et al., 2004).

### 3.4 Augmented State Vector Reconstruction

We reconstruct the augmented state vector and covariance after the measurement update to process the next TDCP measurements. The reconstructed augmented state for processing the next TDCP measurement between time steps *k* and (*k* + 1) is as follows:

36

where is the obtained covariance in Equation (34). Note that the covariance matrix of the augmented state given by Equation (36) is positive semi-definite and therefore valid. See Appendix A for a corresponding proof. The augmented state and its covariance are propagated via time update steps until the next terrestrial GPS measurement at time step *k* + 1. As noted by Gopalakrishnan et al. (2011), this procedure can be viewed as a fixed-point smoothing of the lagged state *x _{k}* based on future predictions from the orbital filter and the terrestrial GPS measurements at time step

*k*+ 1.

Note that the time update of the augmented state covariance matrix of size 18 × 18 shown in Equation (24) can be executed with the same computational cost as updating the 9 × 9 covariance matrix, as shown below in Equation (37):

We expand Equation (24) using (which is similar to Equation (36)):

37

which has the same computational cost as calculating , because the other three sub-matrices of , i.e., , and , can be obtained during the computation of .

### 3.5 Adaptive State Noise Compensation

As explained in Section 1, the dynamic environment changes significantly between the periapsis and apoapsis for elliptical lunar orbits such as ELFOs and NRHOs. Therefore, it would be beneficial to estimate the process noise covariance online through adaptive filtering techniques to automatically tune the process noise covariance to improve performance and robustness.

In ASNC, the empirical process covariance of the position and velocity is calculated as an average over the sliding window of the past *N* time steps as follows:

38

Here, the state update is defined as , where denotes the measurement residual defined in Equation (33) and *K _{k}* denotes the Kalman gain defined in Equation (27). Correspondingly, note that the theoretical process noise covariance matrix and its approximation were defined in Section 3.2.

The objective of this ASNC step is to calculate the covariance of the unmodeled acceleration that minimizes the difference between Equation (38) and Equation (22), assuming that *Q ^{a}* is a diagonal matrix. This is a weighted least-squares problem, and the solution can be obtained analytically as follows (Stacey & D’Amico, 2021):

39

40

41

where:

42

Here, is the first six rows and columns of the covariance of the state update (*S _{p}* is defined in Equation (29)), and (∘) denotes Hadamard (element-wise) multiplication.

After *Q ^{a}* is computed from Equations (39)–(42), it is inserted into Equation (22) to perform the time update.

### 3.6 Outlier Detection

To correctly process TDCP measurements, outliers caused by different factors, such as increased measurement noise or cycle slips, should be detected and discarded. For each TDCP measurement *i* ∈ 1,…*m*′, as discussed for Equation (43), we perform outlier detection by comparing the normalized residual metric, which is formulated as the ratio of the squared measurement residual and its expected covariance . We discard measurement *i* if the normalized residual metric is greater than the pre-defined threshold or if the ratio of the carrier wavelength *λ _{c}* and the expected covariance of the residual is below a predefined threshold . The first criterion discards outliers, whereas the second criterion discards all TDCP measurements whenever we have large residual expectations and it is difficult to distinguish whether the residuals are due to cycle slips or measurement noise. The selection of

*γ*depends on the desired balance between sensitivity and false alarm rate as follows:

_{cs}43

## 4 EXPERIMENTAL SETUP AND RESULTS

We validate our proposed framework using a MATLAB simulation of lunar satellites equipped with a terrestrial GPS receiver and an onboard clock.

### 4.1 Experimental Setup

To verify the proposed algorithm, we simulate the PNT performance of a satellite operating in an ELFO. The starting epoch of the simulation is set to 2022 August 01, 01:00:00, and the simulation time duration is set to four orbits for the ELFO. The orbit is shown in Figure 3(a). The initial orbital elements and the orbital periods are shown in Table 1.

To model the transmit antenna onboard the terrestrial GPS satellites in our simulation setup, we utilize the transmit power and antenna gain patterns of the L1 coarse acquisition signals from the NASA antenna characterization experiment (Donaldson et al., 2020) for Block II-F and the original data from Lockheed Martin for Block IIR and IIR-M satellites (Marquis & Reigh, 2015). Because of a lack of data for Block III satellites, we assumed they have the same antenna patterns as Block II-F satellites. For the GPS receiver, we consider a high-gain steering antenna with 14 dBi at an off-boresight angle of 0° and a 3-dB beamwidth of 12.2°, representative of the current design of the Pretty CubeSat mission (Fragner et al., 2020).

The terrestrial GPS signal from each GPS satellite is considered to be visible when it is not blocked by the Earth and Moon and the *C/N*_{0} continuously exceeds the acquisition and tracking threshold of 15 dB-Hz. All terrestrial GPS signals that pass through a region within an altitude of 1000 km from the Earth’s surface are also discarded to remove signals subject to tropospheric or ionospheric delay. It is assumed that the demodulation of the almanac and ephemeris is not performed onboard but is provided from ground stations or other satellites, as the modulation requires a higher *C/N*_{0} threshold of 26.5 dB-Hz for legacy navigation (LNAV) (Delépaut et al., 2020). Based on the simulation results reported by Delépaut et al. (2020), the mean number of GPS satellites whose ephemeris can be demodulated is 3.08 for an NRHO. Computing the number of GPS satellites whose ephemeris can be demodulated onboard will be addressed as future work.

Figure 3(c) shows the number of terrestrial GPS signals that can be tracked. At some time points, a lunar receiver can track a maximum of 11 terrestrial GPS satellites. However, near the perilune, there exist complete GPS outages, lasting for an average of 48 min per orbit, which corresponds to 6.1% of the orbital period.

Figure 3(d) shows the geometric DOP (GDOP), time DOP (TDOP), horizontal DOP (HDOP), and vertical DOP (VDOP). Here, the HDOP and VDOP were calculated with respect to the body frame of the lunar satellite, where the +z axis points toward the center of the Earth. The VDOP and TDOP values are over 100 times larger than the HDOP values, which implies that we incur larger positioning errors in the line-of-sight direction of the signal. This result occurs because the lines of sight to the terrestrial GPS satellites are all clustered around the boresight direction and lack variation in elevation, as shown in Figure 3(b), because of the large distance between the lunar orbit and Earth compared with the size of the terrestrial GPS orbit.

We sample the position and velocity of the terrestrial GPS satellites at a rate of 1 Hz. Correspondingly, we design our filter with an update interval of Δ*t* = 1.0 s. To simulate the frequency random-walk noise and white frequency noise of the clock in Equation (10), parameters of the mini-rubidium atomic frequency standard (RAFS) (Orolia, 2021) (*σ*_{1} = 1.0 × 10^{−11}, *σ*_{2} = 1.1 × 10^{−15}) are used in modeling the on-board clock. The tracking loop parameters used to compute the thermal noise in the DLL, PLL, and FLL are summarized in Table 3.

The parameters shown in Table 2 are used for the truth dynamics model (used to generate measurement data) and the filter dynamics model (used for time updates). We assume that the filter model is only able to calculate the high-order gravity terms of the Moon up to a lower order and degree compared with the truth model because of the limited computation available onboard.

To evaluate the navigation performance of the filter and compare this performance with the LunaNet SRD specifications, we calculate the root mean square (RMS) and the 68, 95, and 99.7 percentile values of the two values shown below:

44

45

where are the true positions and velocities (in the MCI frame), are the estimated positions and velocities, are the true clock bias and clock drifts, and are the estimated clock bias and clock drifts. We assess the mean performance over 24 Monte-Carlo simulations with different clock motions and state estimate initializations. However, in this work, we do not consider other sources of SISE, such as orbit and clock fitting errors in broadcasted ephemeris parameters (Cortinovis & Iiyama, 2023) or uncalibrated or unknown group delays. Therefore, the position SISE and velocity SISE will be larger than the PCBE and VCDE values, respectively.

Assuming navigation aids from the ground station are provided for initial orbit determination, the initial state estimation error (1*σ*) was set to 100 m for each position axis in the MCI frame, 1.0 m/s for each velocity axis in the MCI frame, 100 m for the clock bias (converted to meters by multiplication with the speed of light), and 0.1 m/s for the clock drift. Each performance metric was calculated by using the last orbital period of all of the Monte-Carlo runs to evaluate the PNT performance after convergence.

### 4.2 PNT Performance

The PNT performance was simulated for the cases in which all three types of measurements (psueduorange, pseudorange rate, and TDCP) were processed with or without ASNC. For the case without ASNC, the unmodeled acceleration was fixed to *Q ^{a}* = diag([1.0 × 10

^{−14}1.0 × 10

^{−14}1.0 × 10

^{−14}]). When ASNC is used, the

*Q*matrix is adaptively updated; its trace value is shown in Figure 4. We can see that the unmodeled acceleration noise grows near the perilune. This increase occurs because the effects of unmodeled high-order gravity terms become more significant near the perilune; moreover, the linearization errors between fixed time steps become larger because of the faster motion of the satellite.

^{a}The obtained estimation error is summarized in Table 4. Both PCBE and VCDE estimates are improved by adding the ASNC compared with using a fixed unmodeled acceleration value. However, despite the improved performance obtained with ASNC, the filter was not able to meet both SISE requirements of NASA, as shown by the cumulative distribution function (CDF) histogram in Figure 6. Potential ways to further reduce the estimation error include incorporating additional sensor measurements, receiving signals from other GNSS constellations (e.g., Galileo), and tuning the ASNC window length, which is currently set to *N* = 10.

The time history of the estimation error and covariance of the lunar satellite states for the ASNC case is shown in Figure 5. From Figure 5, we observe that the position, velocity, and clock bias estimations converge after roughly two orbital periods. This result occurs because the satellite is able to receive terrestrial GPS measurements from various directions while circulating around the orbit. Recall that it is challenging to estimate a lunar satellite’s position using single or limited terrestrial GPS measurements because of the high DOP values in the line-of-sight directions, as explained in Section 4.1. The ballistic solar radiation pressure coefficient *γ* takes a longer time to converge than the other parameters, as it cannot be directly observed from the measurements, and the perturbations created by solar radiation pressure are relatively small compared with gravitational forces.

### 4.3 Outlier Detection

Our outlier detector (as explained in Section 3.6) was tested by randomly adding cycle slips to 30% of the simulated carrier-phase signals. The detection threshold was set to a constant value of *γ _{cs}* = 3. Each cycle slip was simulated as an instantaneous change in the integer ambiguity term from the previous time step. The number of cycles that change from the previous time step when a cycle slip occurs was limited to integers between –5 and +5.

The navigation performance for the scenario with simulated cycle slips is shown in Table 5. When our outlier detector was not incorporated into the filter, the state estimation diverged by unknowingly processing biased estimates. In contrast, our proposed outlier detector was able to avoid degradation by successfully removing outliers, as shown in Figure 7. When cycle slips were present, the implemented detector was able to detect and discard 100% of corrupted measurements. When cycle slips were not present, the detector identified 99.7% as non-corrupted. In simple terms, our detector had a misdetection rate of 0% and a false alarm rate of 0.3%.

### 4.4 Different Terrestrial GPS Measurement Combinations

Different measurement combinations were tested to investigate the contribution of each measurement to the navigation accuracy, as follows: 1) pseudorange only; 2) pseudorange + pseudorange rate; and 3) pseudorange + pseudorange rate + TDCP. The results are summarized in Table 6. As expected, adding the pseudorange rate and/or TDCP measurements to the pseudorange measurement reduces the state estimation errors. In particular, compared with the pseudorange + pseudorange rate solution, which has a 99.7% value of 25.4 m in PCBE and 8.76 mm/s in VCDE, the addition of the TDCP measurement reduces the 99.7% value to 14.0 m in PCBE (a 44.0% improvement) and 5.37 mm/s in VCDE (a 38.6% improvement). These improvements can be attributed to the precise position difference information provided by the TDCP measurement.

### 4.5 Performance in a QFLLO

The PNT performance that can be achieved by using terrestrial GPS signals depends on the geometrical features and dynamic environment of the lunar orbit. To investigate the achievable PNT performance of the proposed algorithm in a more challenging scenario, we also conducted a PNT simulation in a QFLLO. The initial orbital elements of the QFLLO were obtained via the method introduced by Singh et al. (2020), with values shown in Table 7. The generated orbit is shown in Figure 8(a), and the number of GPS signals that could be tracked is shown in Figure 8(b).

It is challenging to performed GPS-based PNT in a QFLLO for two main reasons. First, a satellite at a low lunar orbit experiences a proportionally longer GPS outage period compared with an ELFO, with an average GPS outage period of 32 min per orbit, which corresponds to 25% of the orbital period. Second, there exists a larger discrepancy between acceleration modeled by the truth and filter model, because the effects of unmodeled high-order gravity terms become more significant in low lunar orbits.

A simulation was conducted for eight orbital periods for a QFLLO, with the results summarized in Table 8. Compared with an ELFO, a satellite in an QFLLO experiences larger state estimation errors because of the higher ratio of the GPS outage period to the orbital period. As seen for the ELFO case, the positioning error increases during the outage period, as shown in Figure 9, but the error is still bounded and soon decreases again after the outage period has passed.

Similar to the ELFO, compared with the pseudorange + pseudorange rate solution, which shows a 99.7% PCBE of 121.2 m and VCDE of 115.2 mm/s, the addition of the TDCP measurement reduces the 99.7% PCBE to 48.5 m (a 75.9% improvement) and 3*σ* VCDE to 27.7 mm/s (a 77.1% improvement). The improvement in the PCBE is larger than that of the ELFO case because a satellite in a QFLLO experiences faster motion compared with an ELFO, and therefore, obtaining a precise time-differential position estimate via TDCP measurements more strongly improves the position estimate accuracy.

## 5 CONCLUSIONS

We designed an augmented state Kalman filter framework for lunar satellite PNT that utilizes terrestrial GPS pseudorange, pseudorange rate, and TDCP measurements with the orbital dynamics predicted by the filter. The developed filter accounts for the issue of conditional dependence that occurs in the processing of TDCP measurements by introducing an augmented state vector that incorporates consecutive satellite states. Our proposed framework discards corrupted TDCP measurements (outliers such as cycle slips) by comparing the post-measurement residuals and their estimated covariance. Additionally, our framework adaptively adjusts the process noise covariance by minimizing the difference between the theoretical and empirical unmodeled acceleration covariance.

The performance of the proposed filter was validated through Monte-Carlo simulations of a lunar satellite in an ELFO, equipped with a mini-RAFS. Our work successfully isolated the TDCP measurements suffering from cycle slips and demonstrated 99.7% PCBE and VCDE values below 20 m and 10 mm/s, respectively. The results of the navigation simulation with different measurement combinations indicate that we can reduce the 99.7% PCBE by 44.0% and 3*σ* VCDE by 38.6% by using the precise TDCP measurements, compared with using only pseudorange and pseudorange rate measurements. In addition, we conducted navigation simulations for a QFLLO and demonstrated the applicability of the proposed framework over different lunar orbit regimes. Overall, although the simulation results of our proposed TDCP-based extended Kalman filter framework did not meet the target set by NASA’s SRD, our findings show a promising improvement in position, velocity, and clock bias estimate when compared with the use of an extended Kalman filter with only pseudorange and pseudorange rate measurements.

Future research areas include higher-fidelity modeling of the carrier-phase tracking process to investigate the feasibility of carrier-phase tracking of low *C/N*_{0} signals, simulating the onboard ephemeris and almanac availability in ELFOs and QFLLOs, investigating the required uplink intervals from ground stations, and computing achievable performance when other GNSS signals (e.g., Galileo) are available. Moreover, we may consider simulating the achievable performance for lower-grade clocks, although a higher *C/N*_{0} threshold will likely be required to track the carrier with lower-grade clocks. Finally, we may assess the integration of other sensor measurements and/or satellite crosslinks to reduce positioning and timing errors and to meet the LunaNet SISE requirements.

## HOW TO CITE THIS ARTICLE

Iiyama, K., Bhamidipati, S., & Gao, G. (2024). Precise positioning and timekeeping in a lunar orbit via terrestrial GPS time-differenced carrier-phase measurements. *NAVIGATION, 71*(1). https://doi.org/10.33012/navi.635

## ACKNOWLEDGMENTS

We would like to thank Daniel Neamati for reviewing this paper.

## A PROOF THAT THE RECONSTRUCTED AUGMENTED STATE COVARIANCE MATRIX IS POSITIVE SEMI-DEFINITE

Here, we prove that the reconstructed augmented state covariance matrix given by (from Equation (36)) is positive semi-definite (p.s.d.) for an arbitrary time step *k* when the covariance estimate at the initial time step is p.s.d. We prove this by mathematical induction.

If is p.s.d., is p.s.d.

Proof: From the assumption, is p.s.d. Here, let , where . Then, we have the following:

46

Therefore, the augmented state covariance matrix is p.s.d.

If is p.s.d., is also p.s.d.

Proof: In the time update (as per Section 3.2) and the measurement update (as per Section 3.3), the state covariance matrix is updated via the Joseph form of the Kalman filter equation (Bucy & Joseph, 1987): Equations (24) and (34). In these two equations, is p.s.d. from the assumption, and to satisfy the setup conditions for an extended Kalman filter, the process noise and measurement noise matrices,

*Q*and_{k}*R*, are p.s.d._{k}For any p.s.d. matrix and an arbitrary matrix ,

*TAT*^{⊤}is also p.s.d. because of the following:47

Moreover, the sum of two p.s.d. matrices is also p.s.d. Therefore, from Equation (24), is p.s.d.; based on this finding and Equation (34), is p.s.d.

Finally, we can prove that the following matrix is also p.s.d.:

48

by performing operations analogous to those in Equation (46).

From mathematical induction, is p.s.d. for arbitrary *k*, if the covariance estimate at the initial time step is p.s.d.

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.