Autonomous Crosslink Radionavigation for a Lunar CubeSat Mission

This study presents an autonomous orbit determination system based on crosslink radiometric measurements applied to a future lunar CubeSat mission to clearly highlight its advantages with respect to existing ground-based navigation strategies. This work is based on the Linked Autonomous Interplanetary Satellite Orbit Navigation (LiAISON) method which provides an autonomous navigation solution solely using satellite-to-satellite measurements, such as range and/or range-rate, to estimate absolute spacecraft states when at least one of the involved spacecraft has an orbit with a unique size, shape, and orientation. The lunar vicinity is a perfect candidate for this type of application due to the asymmetrical gravity field: the selected lunar mission, an Earth-Moon L2 (EML2) Halo orbiter, has an inter-satellite link between a lunar elliptical frozen orbiter. Simulation results show that, even in case of high-measurement errors (in the order of 100 m, 1σ), the navigation filter estimates the true states of spacecraft at EML2 with an error in the order of 500 m for position, and 2 mm/s for velocity, respectively and the elliptical lunar frozen orbiter states can be estimated in the order of 100 m for position and 1 cm/s for velocity, respectively. This study shows that range-only measurements provide better state estimation than range-rate-only measurements for this specific situation. Different bias handling strategies are also investigated. It has been found that even a less accurate ranging method, such as data-aided ranging, provides a sufficient orbit determination solution. This would simplify the communication system design for the selected CubeSat mission. The most observable states are found to be position states of the lunar orbiter via the observability analysis. In addition, the best tracking windows are also investigated for the selected mission scenario.


Introduction
There has been an increasing interest in small satellites for lunar missions lately, forming almost 40% of all planned deep space small satellite missions [21].In these missions, the baseline option for orbit determination is, in general, ground-based radiometric navigation.However, ground-based tracking could be expensive and limited considering crowded ground networks.In addition, small satellite missions are expected to be low-cost and there are also limitations from small satellite themselves such as on-board power for communication.Autonomous navigation, on the other hand, could be a possible approach for these lunar missions.Until now, various autonomous navigation methods have been proposed and implemented.
One of them, called LiAISON, uses solely satellite-to-satellite observations, such as range or range-rate, to estimate the absolute states of the involved satellites when at least one of them has an orbit with unique size, shape, and orientation [1,9,10].The characteristics of the acceleration function determine whether inter-satellite range or range-rate measurements can be used alone to estimate the absolute and relative Spacecraft (S/C) states.Considering the asymmetric gravity field in the cislunar vicinity, it is possible to build such an autonomous navigation system.Up to now, various studies have presented the capabilities of LiAISON over the past decade in lunar and deep space mission studies [11,23,12,7,22].
This navigation method will be tested via the link between Lunar Reconnaissance Orbiter (LRO) and the CAPSTONE CubeSat soon [4].However, many more missions targeted to the EML 2 could benefit from such technique.One of these, Lunar Meteoroid Impact Observer (LUMIO) is a CubeSat at EML 2 designed to observe, quantify, and characterize the meteoroid impacts by detecting their flashes on the Lunar farside, to provide global information on the Lunar Meteoroid Environment [3,20,18,19].The baseline navigation strategy, as for almost all small satellites targeting this orbit, is ground-based radiometric.Beside this, LUMIO also includes an inter-satellite link to a larger Lunar orbiter for telecommand purposes but not for navigation.Having an inter-satellite link provides an opportunity to investigate the performances of autonomous navigation via crosslink radiometric measurements, potentially extending the mission possibilities with this new technique.
This study presents the autonomous navigation performances of the selected mission scenario: LUMIO, via the link between Lunar Pathfinder (LPF) and the LUMIO CubeSat.A simulation-based analysis will determine the achievable orbit determination accuracy considering realistic radio frequency measurement errors derived from the Phase-A study inter-satellite link design, which was not designed to perform navigation.
In the following sections, at first, the selected mission scenario is presented and then, dynamical models are provided.Orbit determination models are introduced including the observability analysis.Thereafter, the navigation simulation setup and results are presented.Finally, conclusions are drawn.

A Lunar CubeSat scenario
For this study, the LUMIO mission has been selected.It features a 12U CubeSat in a halo orbit at the EML 2 to observe, quantify, and characterize meteoroid impacts on the Lunar farside by detecting their flashes [20,3].The mission aims at determining the spatial and temporal characteristics of meteoroids impacting the Lunar surface to characterize their flux.The operative orbit for LUMIO has been selected as a quasi-periodic halo orbit with a Jacobi constant C j = 3.09 [3].The LUMIO Phase-A study has been completed in March 2021 while the mission's Phase B is expected to start soon [19].The radiometric navigation system uses existing hardware of the communication system featuring a combination of Inter-Satellite Link (ISL) and Direct-to-Earth Link (DTE) links.The latter has been designed to provide payload data downlink, ranging and tracking in nominal conditions.Based on the orbit determination (OD) analysis given in [17], ground-based radiometric navigation via Cebreros, ESTRACK or the Sardinia Deep Space Antenna for 3 hours per track following a 7 + 7 + 14 days scheme meets the OD requirements of 1 km and 1 cm/s position and velocity accuracy, respectively.In this study the same OD requirements have been considered for an autonomous navigation scenario.The ISL, on the other hand, has been designed to provide a redundant commanding link without involving a dedicated deep-space class ground station but reusing commercial resources.Such link provides optimal performances in terms of visibility, despite the fact that data rates are quite limited.The SSTL LPF spacecraft has been considered as relay satellite.Depending on the relative distance, between a minimum of 31 000 km and a maximum of 89 000 km, data rates are expected in the order of 0.5 kbps-4 kbps based on S-band, 9 dBW Equivalent Isotropic Radiated Power (EIRP) link (also including a 3 dB safety margin [19]).

Autonomous Radiometric Navigation
This study investigates the LiAISON orbit determination performances for the LUMIO mission: autonomous orbit determination requires absolute position and velocity estimates without any ground-based observation.This requires that the S/C states, obtained from inter-satellite observations, must be observable.
In the two body problem, relative range or range-rate measurements do not provide full state estimation due to the rank deficiency: absolute orientations of the orbital planes are not observable, which is related to the symmetrical gravity field.However, in an asymmetrical gravity field, spacecraft-to-spacecraft tracking provides absolute state estimation.The LiAISON orbit determination method uses solely inter-satellite measurements to estimate absolute S/C states of the involved S/C when at least one of them has an orbit with a unique size, shape and orientation [1,9,10].Basically, if one of the S/C has an unique orbit, it is possible to estimate absolute states of all involved spacecrafts.This study derives inter-satellite observations, namely range and range-rate, from radio frequency measurements.
Inter-satellite radiometric measurements can be collected with various methods.Range observations, for example, can be collected via phase or time-based measurements: the phase shift on a ranging signal at the receiver with respect to a transmitter provides an accurate ranging solution (in the order of 1 m).This still requires the phase ambiguity to be solved in case the ranging signal period exceeds the round-trip light time.
There are lots of common standards used for this purpose such as pseudo noise/sequential tone-based ranging [21].However, small satellites, in general, have limited on-board power available for data transfer.If a small satellite requires ranging for its navigation, such signal reduces the power available for telemetry in case both signals are modulated at the same time window.On the other hand, timing exchange between satellites provides also a ranging solution.Basically, in order to deal with on-board power limitations, time-derived ranging methods are a possible options for small satellites.These methods are not quite accurate (in the order of 150 m at 10 kbps) but they can provide sufficient ranging solutions to meet the navigation requirements.
Because time-derived ranging methods are data-rate dependent, high-data transfer between satellites is beneficial.Another measurement type is range-rate and this can be derived from the Doppler shift between satellites but, again, this measurement type suffers from on-board limitations.In this study, high and low measurement ranging errors derived from the conventional (radiometric) and time-derived ranging methods are considered for the selected mission scenario.

Orbit Determination Models
This section presents the orbit determination models used in this study.Dynamical, measurement and estimation models are given in the following subsections.

Dynamical Model
The dynamical model used in this study is formulated as Circular Restricted Three-body Problem (CRTBP).This is simple but accurate enough for this type of application, considering the mission required position accuracy.In [1], the LiAISON OD performances remain in the same order of magnitude for various force models.In addition to CRTBP, in the end, high-fidelity dynamic simulation models results will also be given for the selected measurement configuration to have a more realistic analysis and to compare results obtained from CRTBP dynamics.
The CRTBP assumes there are two massive bodies, Earth (P 1 ) with mass m 1 and Moon (P 2 ) with mass m 2 in this case, moving under their mutual gravitation in a circular orbit around each other with a radius r 12 , [5].Considering a non-inertial, co-moving reference frame (see Figure 1) with its origin at the barycenter of the two bodies, the positive x−axis points from the barycenter to P 2 .The positive y−axis is parallel to the velocity vector of P 2 and the z−axis is perpendicular to the orbital plane. 2 where For the Earth-Moon system, the gravitational parameter µ is 0.01215, the normalized time t * 4.343 days, and the normalized length l * 384 747.96 km, respectively.The transformation from the rotating frame to the inertial frame can be done as follows [8]: where and where X rot represents the non-dimensional Earth-Moon barycentric coordinates, X inert represents the inertial state centered on the barycenter.It is assumed that two frames coincide at the initial time t 0 = 0.A rotating position vector of the S/C in the Moon-centered frame is as follows: and in the Earth-centered frame: Based on the equations given above, states can be transformed from the Earth-Moon barycenteric rotating frame to the primary-centered, either Earth or Moon, inertial frame.In the study, the Lunar orbiter's initial states are expressed in the Moon-Centered Inertial (MCI) frame, a coordinate transformation is needed to convert to the non-dimensional Earth-Moon barycentric frame.Basically, these steps are followed in reverse order.
In high-fidelity dynamical simulations, in addition to gravitational acceleration due to Earth and Moon (treated as point masses in the Earth-Centered Inertial (ECI) J2000 frame), the gravitational acceleration due to Sun treated as a point mass and acceleration due to Solar Radiation Pressure (SRP) with a spherical model have been used.

Measurement Model
In this study, the autonomous navigation method uses solely inter-satellite measurements and the observables are collected via radiometric measurements.The first radiometric data type used is range, which can be derived from either conventional ranging methods or data-aided ranging.A conventional ranging method using is pseudo-random noise has been selected and, considering a non-coherent transponder with a pseudo-noise square wave shaped ranging signal, and a chip tracking loop, the following one-way ranging error (1σ) would be expected [2]: and the range bias due to a chip rate mismatch: with c the speed of light, f rc the frequency of the ranging clock component, B L one-sided loop noise bandwidth, P RC power of the ranging clock component, T integration time, N 0 one-side noise power spectral density, ∆f chip the difference in frequency between the received chip rate and the local chip rate.
The round-trip light time can also be measured via time transfer between satellites.This process requires four successive time stamps to be obtained which represent the time of transmission and reception of both S/C (see Figure 2).If the timing is measured in units of telemetry/telecommand symbols, instead of directly in seconds, the performance of time-derived ranging (one-way) would be: where T sd is the symbol duration, T l the correlator integration time and E S /N 0 the symbol-to-noise ratio.
Note that both equations 9 and 11 represent one-way ranging performance and this requires calculations for both uplink and downlink.The measurement error for two-way Doppler due to thermal noise, influencing range-rate observations, can be given as follows [6]: where f c the downlink carrier frequency, P C /N 0 uplink carrier power to noise spectral density ratio, ρ L the downlink carrier loop signal-to-noise ratio, G the turn-around ratio.Doppler data noise can be expressed by the phase noise, σ ϕ , in radians and converted to range-rate noise as follows [13]: In this study, the formation includes two S/C, a lunar orbiter, LPF, and EML 2 Halo orbiter, LUMIO.
The state vector being estimated consists of the position and velocity components of both S/C is as follows: The measurement model in this paper, referred as the pseudo-range, involves the geometric range, the overall clock bias, and other error sources.The two-way ranging measurement concept can be seen in Figure 2. The geometric range is given as follows: By ignoring the light-time correction, and by assuming the speed of light is greater than the S/C relative velocity, c v, the geometric range can be modeles as where x i , y i and z i represents the position components of S/C, i = 1, 2 states and the pseudo range observations can be modeled as where ψ t4 and ψ t1 are the clock states at t 4 and t 1 respectively.The transponder transmit and receive line delays are ∆ tx and ∆ rx , respectively and ∆ trx is the line delay on the S/C transponding the ranging signal.
All these terms are combined as ρ bias and ρ noise representing the un-modelled statistical error sources.
The range rate measurements, ρ, can be modelled as: ρ noise and ρnoise are calculated based on equations given in 9, 11, and 12.

Estimation Model
In this study, an Extended Kalman Filter (EKF) is adopted as a common method used in real-time navigation.The EKF consists of a prediction and a correction step: in the former, predicted state and error covariance P are [14]: where Φ(t k , t k−1 ) is the state transition matrix from t k−1 to t k and Q is the process noise matrix.The correction step is: where X is the state estimate, K is the Kalman gain, H is the measurement sensitivity, P is the error covariance estimate, and W is the state noise compensation matrix.
In this study, the state noise compensation is introduced by Q which can be constructed for each S/C as follows [10]: The measurement bias is expected to affect the the navigation system performances: in general, this can either be neglected or estimated by including dynamic or measurement model parameters into the state vector.Another approach would be to assume its a priori estimate and associated covariance matrix are known.In this study, all three cases (neglected bias, estimated bias and considered bias) were investigated.
In case of bias estimation, the estimated state vector must be expanded with a bias component, ρ bias .This requires H and Φ to expand as follows: In a similar way, the clock drift can be estimated by expanding H and Φ with (t − t 0 ) representing the step size and 1 representing unchanged clock parameter during the time update process, respectively.However, expanding them into the state vector may affect the performances by reducing observability of the navigation system.That's why this study also investigated the considered bias case by implementing a sequential consider filter (known as Consider-Kalman Filter (CKF), Schmidt-Kalman filter (SKF)).Here, the specific time invariant measurement bias case is studied assuming that consider parameter, b k , is constant for all k and the bias covariance matrix, B k , is time invariant.CKF is slightly different than EKF and requires to implement the following equation into the time update: and changes in measurement update: where where C k is the cross-covariance matrix, B 0 is the bias covariance, b 0 a priori the measurement biases estimate, and N k the measurement bias vector sensitivity matrix.Note that the CKF turns into EKF in case of zero bias.

Observability
The observability analysis is used to relate OD performance and observation data: in this study, the degree of the system observability is used to evaluate the estimation performances.For this purpose, the observability Gramian is used as follows: The observability can be assessed via Singular Value Decomposition (SVD) of the observability Gramian.
In general, two metrics are used for this purpose: the condition number, which is the ratio of the largest singular value to the smallest one, and unobservability index, which is the reciprocal of the smallest local singular value.Most and least observable states can also be derived from unitary matrices via SVD.

Navigation Simulations
This section presents the orbit determination results for the selected mission scenario.The simulation setup is given first and results are given thereafter in the corresponding subsections.

Simulation setup
The selected mission scenario consists of two S/C, the LPF and LUMIO, at the Elliptical Lunar Frozen orbit and EML 2 Halo orbit [15,3].The simulation duration is set to be 14 days.For the the N-body orbital dynamics based analysis, the simulation start time and duration are set to be 18 April 2024, 21:00:00 UTC.
The initial states of LPF expressed in MCI are listed in Table 1.LUMIO has a quasi-halo orbit at EML 2 with a Jacobi energy of C j = 3.09.For this study, due to the duration of the simulation, a similar southern Halo orbit (with the same Jacobi energy) has been used in CRTBP.The initial conditions can be found in Table 2 including the LPF converted states from MCI to the non-dimensional Earth-Moon barycentric frame.
Corresponding trajectories of both S/C for 14 days can be seen in Figure 3.In the mission scenario, two-different measurement error cases have been simulated which were representing the high-accuracy conventional pseudo-noise ranging method and the time-derived ranging method.Intersatellite ranging is relaying on the communication system and the link budget for LUMIO is presented in  3.In this table, different from the existing link budget, the uplink data-rate and thus symbol duration is increased in consideration of a higher gain antenna configuration or higher transmit power.This is reasonable because the existing S/C configuration did not take into account the inter-satellite link based autonomous navigation (the existing configuration would give 2700 m 1σ ranging accuracy which is not sufficient for the mission).Based on these assumptions, range and range-rate measurement performances can be seen in Table 3.
During the simulations, in each kth time step, the Root Mean Square (RMS) error for the N th case of the Monte Carlo simulation has been calculated by using following where x i,k and xi,k are ith component of state vector and its estimate respectively.Parameters used in simulations can be found in Table 4. Regarding filter uncertainty and initial errors, a more detailed analysis has been done during the LUMIO Phase-A design study, previously.A ground-based tracking for 7 hours between 18 April 2021 14:00:00 and 21:00:00 UTC from the Sardinia Deep Space Antenna based on range and range-rate measurements (measurement errors of 1 m range and 0.33 mm/s range-rate with a measurement bias of 2.5 m) would give 0.11 km and 0.95 cm/s position and velocity errors with 1σ uncertainty of 1.65 km and 4.7 cm/s, respectively (Estimation time Epoch: 18 April 2021 21:00:00 UTC with Batch-Least squares).
Basically, this single tracking session estimation can be used to initialize the autonomous navigation system on-board.A ground-based state estimation has not been done for LPF.However, it is assumed that groundbased state estimation results for LPF is the same order of magnitude.In brief, initial parameters given in  As already mentioned, OD requirement of the EML 2 orbiter is 1 km for position and 1 cm/s for velocity, respectively.This has been taken as baseline goal for this study.The trajectory used in the study is also considered as a true reference, so there is no error in the dynamics.In addition to CRTBP dynamical models, the N-body orbital dynamics simulations are also performed.In such model (based on the JPL DE405 ephemeris model), Earth, Moon, and Sun are treated as point masses.Regarding the SRP model, SRP areas are set to 3 m 2 and 0.41 m 2 , and reflectivities are set to 1.8 and 1.08, for LPF and LUMIO, respectively [15,16].The other settings for the high-fidelity analysis are the same.

Results
This section presents the performance of radiometric autonomous navigation for the selected mission scenario.The effects of measurement accuracy, precision, data type and navigation filter on the OD performance have been investigated.
At first, the baseline case has been presented to show the autonomous navigation method works in the lunar vicinity.This case is based on the inter-satellite ranging derived from the conventional Pseudo-Noise (PN) method.Based on the settings given in the previous section, the navigation filter estimates the true states of LUMIO in the order of 100 m for position and 1 mm/s for velocity, respectively.LPF states are estimated within the order of 10 m position and 1 cm/s velocity, respectively.Estimation results can be seen in Figure 4 including RMS error and covariance values.As it can be seen, the filter converged after day 6.
This is due to the fact that halo orbit has a period of 14 days and the half-orbit is sufficient to fit for the EML 2 orbiter.Fluctuations in the LPF state estimation are related to the relative geometry between S/C.The position estimation converges when LPF approaches the periselene, and diverges when LPF approaches the aposelene.It is beneficial for LPF position estimation to be performed when the LPF is in the high velocity region.In case range-rate measurements are used, instead of range, the filter estimates are not improved (see Figure 5).Monte-Carlo simulation results can be seen in Table 5.
As part of the analysis, it was determined that the most observable states are in order: z 2 , x 1 , y   and 4.324 × 10 12 for range-only and range-rate only case, respectively.This shows that s range-only system has a higher observability than a range-rate only system.However, the range-rate only case provides higher information to the filter on the least observable state: ż via the lower unobservability index (0.5829 for anytime via any type of data transfer between S/C.However, this method is not as accurate and this would affect the autonomous navigation performances.In brief, this method would simplify the communication system design and reduce the on-board power required, making this method a perfect candidate for this mission.The time-derived ranging method simulation results are given in Figure 9.As it can be noticed, the filter can estimates the true states of LUMIO in the order of 140 m for position and 1 mm/s velocity, respectively.Errors are almost 4 times higher than the conventional ranging method in simulation.However, this still meets the mission navigation requirements (1 km).This shows, from the navigation performance perspective, that collecting measurements at different time intervals (and thus orbit geometry) has more importance than the absolute measurement errors.As a side note, even in the high measurement error case (time-derived ranging), the system is observable enough to estimate the systematic bias.This only requires more time than in the PN case.
In the last simulation of the paper, the ephemeris model based simulation results are given.This shows the effects of additional perturbations in the system on the navigation performances.Based on the ephemeris model (including SRP and the gravitational acceleration due to Sun), the estimation results are shown in Figure 10: position and velocity estimation uncertainty increased up to three times in this case with respect to the CRTBP case.However, the LUMIO estimation errors are lower than 500 m 1σ for position and 2 mm/s for velocity, respectively, and for LPF, 1σ uncertainties are less than 100 m for position and 1 cm/s velocity.
This would meet the OD requirements of 1 km position and 1 cm/s velocity, respectively.

Conclusion
This study showed that the LiAISON orbit determination technique could be a possible navigation approach for a proposed mission, LUMIO, based on the existing inter-satellite link between LPF and the LUMIO CubeSat without using any ground based measurements.Simulation results show that the navigation filter estimates the true states of LUMIO in the order of 500 m for position, and 2 mm/s for velocity, respectively.LPF states can be estimated in the order of 100 m for position and 1 cm/s velocity, respectively.Considering the range-only case with 1σ error of 2.98 m provides better states estimation than the range-rate only case with 1σ measurement error of 0.97 mm/s.The observability analysis showed that the system is observable and found that LPF's position and velocity components are the most and the least observable states, respectively.The best tracking windows are also presented by means of the observation effectiveness analysis.It has been showed that bias would affect the performances, and it can be estimated along with the dynamical states, thanks to the high observability.In addition, the time-derived ranging method provides sufficient information to the filter in order to meet the navigation requirements.In these simulations, only the initial phase of the operative orbit has been considered.Basically, after one orbital period, it is quite expected to improve the estimation due to relative position change between S/C.This would bring additional information to the filter and decrease the uncertainty.This study considered only the first 14 days of the mission.Even though this is sufficient to see the expected performances, further work would be needed to investigate what would happen for the full operative mission lifetime (1-year).This study also didn't consider the dynamic model errors and the effects of clock drift on the performances and they are considered as topics for future research.

Figure 1 :
Figure 1: Circular Restricted Three-body Problem: Earth-Moon co-moving reference frame

Figure 2 :
Figure 2: Time exchange between satellites for the purpose of time-derived ranging.t i represents time, ψ i represents onboard clock states at t i

Figure 3 :
Figure 3: Trajectories of both S/C for 14 days in the non-dimensional Earth-Moon barycentric frame

Figure 4 :Figure 5 :
Figure 4: State estimation results based on range-only measurements

Figure 9 :
Figure 9: State estimation results based on the time-derived ranging method (Monte-Carlo simulation 100-execution)

Figure 10 :
Figure 10: State estimation results based on the ephemeris model with time-derived ranging method)

Table 2 :
Initial states used in the simulations (expressed in non-dimensional Earth-Moon barycentric frame)

Table 3 :
[19]mptions for LUMIO (the link budget is given based on the maximum inter-satellite distance[19].Range measurement parameters are for conventional and time-derived methods and measurement errors are 1σ, two-way)

Table 3 .
Inter-satellite range measurement parameters for both conventional and time-derived methods are given in Table

Table 4
are sensible.

Table 5 :
100-execution Monte-Carlo simulation results (values in parenthesis represents the results after day 6, and values are averaged over two S/C)