Precise Positioning and Timekeeping in a Lunar Orbit via Terrestrial GPS Time-Differenced Carrier-Phase Measurements

.

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.

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(Capuano et al., , 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(Capuano et al., , 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(Bhamidipati et al., , 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.

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.

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 FIGURE 1 Our proposed framework performs precise positioning and timekeeping onboard a lunar satellite by leveraging the intermittently available terrestrial GPS TDCP measurements.The figure is not to scale.
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.

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.

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.

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.

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: where γ denotes the ballistic coefficient of the solar radiation pressure, t k denotes the time based on GPS time, r u ∈ R 3 and v u ∈ R 3 denote the position and velocity of the lunar satellite in the Moon-centered inertial (MCI) frame, respectively, and  rv denotes the associated process noise.In this paper, we use a Moon-centered J2000 frame for the MCI frame.Additionally, 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): where U g and ∇U g denote the gravitational potential and its gradient, respectively, 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 p is the corresponding mass.At each time step, 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.Pnm is the normalized Legendre function, and Cnm and Snm are the normalized geopotential coefficients.For the Moon-fixed frame, we used the "IAU_MOON" frame provided in the NAIF SPICE library.N sph is the maximum degree and order of the higher-order gravity terms.
In Equation (1), the acceleration due to solar radiation pressure a  is calculated from the cannon ball model as follows: where Φ  denotes the solar flux at 1 astronomical unit (AU) (= 1360 W/m 2 ; 1 AU = 1.498 × 10 8 km) and A u and M u denote the surface area and mass of the lunar satellite, respectively.Additionally, d denotes the heliocentric distance to the lunar satellite, ˆSun e denotes the sun unit direction vector from the lunar satellite, and C R 1  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.458km/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 J  be the vector of states related to the orbital dynamics.To obtain the linearized dynamics, consider the linearized state s k around the reference state S k * as follows: 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 r v ( , , , ) t u u γ from Equation (1) as a S ( , ( )) t t , we obtain the linearized dynamics as follows: where: By propagating Equation (6) from t k−1 to t k , we have the following: R is the state transition matrix.The state transition matrix can be computed by propagating the following equation from t k−1 to t k :

Clock Dynamics
The clock bias W u R and clock drift  W u R are propagated by using the following discrete, two-state error model taken from Galleani (2008): where the noise vector ε ε ε is Gaussian noise sampled from the random-walk process of the clock.

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 I using carrier-phase values at time steps k and ( ) k − 1 , i.e., φ k i and I k i 1 , respectively, is given by the following: where, eliminating subscripts k and k − 1 for the sake of simplicity, τ u is the receiver's clock bias and  φ i is the carrier-phase measurement noise.For any terrestrial GPS satellite i i , r 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 τ i is the satellite clock error.Additionally, ∆ 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 ' T i  I .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).

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.

Augmented State
Let  X k u R 9 1 be the lunar satellite state vector defined in Equation ( 13): The lunar satellite state vector  X k includes the position r u , velocity v u , solar radiation pressure reflection coefficient C R , clock bias τ u , and clock drift  τ u .The position 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), e k u −1 and e k u are functions of  X k−1 and  X k , respectively.Therefore, the TDCP measurement Kaplan and Hegarty (2017)) lunar satellite states.This breaks one of the assumptions of the Kalman filter, which is the conditional independence of  X k−1 and  z k , given  X k , as follows: 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:

Time Update Step
We first define the linearized state vector  x around the reference lunar satellite state (i.e., current estimate of the state in the extended Kalman filter) vector  X * as follows: Based on Equations ( 16), ( 8), and ( 11), the linearized dynamics of the lunar satellite state vector  x k at any k time step are given by the following: ) where  k denotes zero-mean Gaussian process noise, whose covariance  Q k depends on the noise covariance in the lunar satellite position and velocity, which is given by  Q k rv , and that in the lunar satellite clock, which is given by  Q clk .Note that the state transition matrices ) , 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  Q clk is modeled as , are the diffusion components of the clock phase deviation and clock frequency deviation.
The process noise covariance  Q k rv is formulated as follows (Tapley et al., 2004): where Q a u R 3 3 is defined as the unmodeled acceleration, R 6 6 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): Next, we formulate the linearized augmented state vector x k given by , where  x k was introduced in Equation ( 16).In the time update of our proposed framework, we propagate the linearized augmented state x k and its covariance P k as follows: where x k and P k are the predicted augmented state vector and its covariance at time step k, respectively.Similarly,  0 0 I with the state transition matrix  )( , ) t t k k1 defined in Equation ( 18).In addition, the process noise covariance Q k is given by Q , where the process noise covariance matrix  Q k of the linearized state vector  x 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 t k−1 and thus remains unchanged.

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 ˆk x (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., ˆ( 1 )).
The full estimate ˆk x and its covariance ˆk P 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 I. Let z k be the stacked measurements obtained at time step k as follows: where m and ′ m are the number of available terrestrial GPS signals and TDCP signals, respectively.The total number of measurements will be M m m c 2 .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 2 m ( ) σ to the pseudorange and 0 01 2 .
/ ( ) m s σ to the pseudorange rate as a SISE.This value is based on the SISE values ( m % ≤ 7 0 95 .user range error and ≤ 0 006 95 ./ m s % 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 R 9 9 × , as follows: We then calculate the Kalman gain as follows: Here, H k denotes the linearized measurement model (defined in Equation ( 28)) that stacks the linearized measurement models associated with the pseudorange H k ρ, , pseudorange rate H k  ρ, , and TDCP H T k ' I, .The augmented state and its covariance will be reconstructed from these values following the procedure shown in Section 3.4.For any i-th satellite, the linearized measurement models of the pseudorange H k i ρ, , pseudorange rate H k i  ρ, , and TDCP H T k i ' I, are defined in Equations ( 30), (31), and (32), respectively: With our measurement vector z k from Equation ( 25), the modeled measurement covariance R k , and the updated Kalman gain K k from Equation ( 27 where 1 ˆM k y × ∈  denotes the measurement residual, i.e., the difference between the received measurements z k and the expected measurement for the predicted augmented state vector X k (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).

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: , , , ,

ˆˆˆ and ˆˆk
where , ˆk k P 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 ˆk x and its covariance ˆk P 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)): 1, 1 1, 1 1 9 9 9 9 1 9 9 9 9 9 9 9 9 9 9 9 9 9 9 1, 1 1, 1 which has the same computational cost as calculating because the other three sub-matrices of P k , i.e., ˆ, k k P − − can be obtained during the computation of P k k , .

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: ( ) Here, the state update ∆ k x is defined as ˆ, where ˆk y 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 Q a u R 3 3 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): where: Here, 6 p rv u R 6 6 is the first six rows and columns of the covariance of the state update 6 p p p p p K S K S u  R 9 9 ( 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.

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 m } c 1, , 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 ( ) .We discard measurement i if the normalized residual metric is greater than the pre-defined threshold γ cs 2 or if the ratio of the carrier wavelength λ c and the expected covariance of the residual is below a predefined threshold γ cs 2 .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 γ cs depends on the desired balance between sensitivity and false alarm rate as follows:

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.

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) ) 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: where ( , , , , , ) x y z x y z    are the true positions and velocities (in the MCI frame), ˆˆ( , , , , , ) x y z x y z    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.

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 u u u diag .([ . . .]) 1 0 10 1 0 10 1 0 10 14 14 14 When ASNC is used, the Q a 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.
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.

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 J 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%.

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 Navigation Results (24 Monte-Carlo Runs) When Cycle Slips are Added to the Simulator The filter crashed because of the ill-matrix condition when our outlier detector was not incorporated in the filter.In contrast, our proposed cycle slip detector avoided this degradation by successfully removing outliers.The percentage of corrupted measurements in any Monte-Carlo run denotes the percentage of measurements among the total that are induced with cycle slips.

Position and Clock
Bias Error (m)  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.

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.

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.

a c k n o w l e d g m e n t s
We would like to thank Daniel Neamati for reviewing this paper.

FIGURE 3
FIGURE 3 The lunar satellite's trajectory and skyplot and the statistics of terrestrial GPS signals tracked from an ELFO (a) ELFO shown in the Moon-centered J2000 (b) Skyplot at a chosen time epoch where 11 terrestrial GPS signals are tracked.The GPS satellites are clustered in the satellite's boresight direction.(c) Number of terrestrial GPS satellites tracked in the ELFO.The gray-shaded region in (c) illustrates the terrestrial GPS outage period (48 min per orbit, which corresponds to 6.1% of the orbital period).(d) DOP value of the terrestrial GPS signals.The DOP values are defined in the satellite body frame, where the z axis is the antenna boresight direction.The VDOP and GDOP values are over 100 times larger than the HDOP values.

FIGURE 4
FIGURE 4 Trace of the unmodeled acceleration matrix ( )Q a computed by the ASNC algorithm The unmodeled acceleration noise grows near the perilune, where the effects of unmodeled highorder gravity terms become more significant and a larger linearization error occurs because of the fast motion of the lunar satellite.

FIGURE 5
FIGURE 5 Estimation error (red line) and 3σ covariance (red shaded area) of the lunar satellite states for a Monte-Carlo simulation in an ELFO The position and velocity state errors are represented in the radial-tangential-normal (RTN) frame with respect to the true orbit.The estimation errors in the normal direction are larger compared with the radial and tangential directions, as they are more aligned with the GPS line-of-sight direction, where pseudorange error exists.The gray shaded area shows the terrestrial GPS outage period near the perilune.The duration corresponds to four orbital periods, and the simulation begins at the apolune.The state estimates other than the solar radiation pressure (SRP) coefficient γ converge after approximately one orbital period.Spikes in the state covariance estimate occur during the terrestrial GPS outage.

FIGURE 7
FIGURE 7 Residuals of the TDCP measurements, shown in meters The residuals from a single run are split into two figures: the left figure shows the residuals of the cycle-slip-corrupted measurements, and the right figure shows the residuals of the non-corrupted TDCP measurements.The black dashed horizontal lines show multiples of the carrier-phase wavelength.The outlier detector discards 100% of the cycle-slip-corrupted measurements and accepts 99.7% of the non-corrupted measurements.

FIGURE 8 FIGURE 9
FIGURE 8 Trajectory of the lunar satellite and the number of GPS satellites tracked in the QFLLO (a) QFLLO shown in the MCI J2000 frame (b) Number of terrestrial GPS satellites tracked in the QFLLO

TABLE 1
Initial Orbit Elements and Orbital Period of the ELFO The orbital elements are given in the J2000 frame.

TABLE 2
Dynamics Model Used in the Simulation

TABLE 3
Tracking Loop Parameters Used in the Simulation

TABLE 4
PNT Performance Obtained with Our Proposed Framework in an ELFO with an Onboard Mini-RAFS "Below Req." denotes the time ratio at which the PCBE/VCDE ratio falls below the SISE requirements of the LunaNet SRD (13.43 m for SISE of position and 1.2 mm/s for SISE of velocity).

TABLE 7
Initial Orbit Elements and Orbital Period of a QFLLO The orbital elements are given in the J2000 frame.

TABLE 6
Navigation Results (24 Monte-Carlo Runs) for Different Measurement Combinations Compared with the pseudorange + pseudorange rate solution, the addition of the TDCP measurement reduces both the PCBE and VCDE.

TABLE 8
Navigation Performance for a Satellite in a QFLLO Equipped with a Mini-RAFS (24 Monte-Carlo Runs) Each performance metric is calculated from the last two orbital periods of all of the Monte-Carlo runs for eight orbital periods each."PR" and "PRR" in the measurement column denote pseudorange and pseudorange rate measurements, respectively.