Application of Pulsar-Based Navigation for Deep-Space CubeSats

Abstract

This paper investigates the use of pulsar-based navigation for deep-space CubeSats. A novel approach for dealing with the onboard computation of navigational solutions and timekeeping capabilities of a spacecraft in a deep-space cruise is shown, and the related implementation and numerical simulations are discussed. The pulsar’s signal detection, processing, and exploitation are simulated for navigation onboard a spacecraft, thus showing the feasibility of autonomous state estimation in deep space even for miniaturized satellites.

1. Introduction

The vast majority of space missions have strongly relied—and still do—on Earth-based methods for determining the spacecraft’s position and velocity in deep space. Ground-based radiometric tracking, which is the state of the art for deep-space navigation, can reach a positioning accuracy of a meter at distances equivalent to Jupiter’s distance from Earth [1], yet requires dedicated ground infrastructure and the involvement of flight dynamics teams that lead to increased overall space-mission costs. This becomes unsustainable as the number of deep-space satellites increases. Autonomous navigation systems that do not need the constant aid of Earth-based infrastructure are desirable. This is to enable sustainable solar system exploration while allowing for a reduction in the space missions’ costs.

Autonomous navigation in deep space has been proposed through optical methods [2]. Celestial triangulation exploits the line-of-sight measurements to distant planets [3], which are used as navigation beacons in order to determine the observer’s position. This can be performed since the ephemeris of the planets are accurately known and the relative line of sight can be tracked while in deep space [4]. However, several planets need to be visible, which is not always the case for missions in deep space due to the presence of the Sun [5]. Other optical methods rely on close observations of moons [6] and asteroids [7], but in these cases, the navigation solution can only be computed when in proximity to these bodies. Therefore, optical navigation options are well suited for when geometries are favorable [8] and the motions of the observed bodies are well-characterized.

More recently, navigation methods exploiting the observation of celestial pulsed X-ray sources have been proposed [9,10,11,12,13,14,15,16,17]. These methods rely on the pulsar’s signal time-of-arrival acquisition to estimate a spacecraft’s position in deep space, comparing the acquired signal profile to an expected profile at the solar system barycenter. Within this context, the purpose of this work is to investigate the use of the X-ray pulsar-based navigation in a deep-space CubeSat mission scenario, assessing the feasibility and the related performances of autonomous navigation. This is performed by adopting performances and constraints of typical deep-space CubeSat missions.

The structure of the paper is as follows. Section 2 describes the methodology of the pulsar navigation concept and provides background regarding the physics of these stellar objects as well as signal modeling. Section 3 delves into the modeling of the measurements, reports the implementation procedures of the method, and presents the simulation performances in two mission scenario. Eventually, conclusions are drawn in Section 5.

2. Methodology

This section deals with the core methodology of the X-ray pulsar navigation, highlighting the mathematics behind the pulsar signals’ acquisition and exploitation for deep-space navigation.

2.1. Signal Time-of-Arrival

Pulsars are rapidly rotating, highly magnetized neutron stars consisting of extremely compressed matter. The formulated models regarding their nature predict a maximum mass up to twice the one of the Sun, possibly increasing when the magnetic field is particularly strong. Conversely, the radius is predicted to be around the size of a city. As they spin, charged particles are accelerated out along magnetic field lines in the magnetosphere. This acceleration drives the particle to emit electromagnetic radiation.

A beam of radiation is then emitted along the magnetic axis of the pulsar, which, in general, differs from its rotational axis. This misalignment causes the beam to be seen by an observer only once per rotation of the neutron star and leads to the “pulsed” nature of its appearance—hence the name pulsar. The repetition period of the pulses is simply the rotation period of the neutron star, and it is seen to be very regular and stable for most of these celestial objects. Therefore, the pulsar source emits a wave-front of photons in the directions of emission, hitting the observers in the solar system at precise epochs. Let us denote $\mathbf{n}$ as the direction of emission of a given pulsar. The arrival times of the pulsars’ wave-fronts at the solar system barycenter (SSB) are known by models developed by astronomers. The photons belonging to the front wave will arrive at an observer location with a given time that is different from the arrival at the SSB. Hence, the distance d between the observer location and the SSB is proportional, within the first order, to the time delay $t_d$ between the two timing locations, as shown in Figure 1.

The equation of the distance is:

$$d=c\left(t_{s/c}-t_{SSB}\right)=\mathbf{n}·\mathbf{r}+H$$

(1)

where c is the speed of light in vacuum; $t_{s/c}$ and $t_{SSB}$ are wave-front arrival epochs at the spacecraft and SSB locations, respectively; $\mathbf{r}$ is the position vector of the spacecraft with respect to the SSB; $\mathbf{n}$ the unit direction of the source of the photons expressed in the same inertial reference frame as the observer position; and H represents higher-order terms of various natures, which will be discussed later. By using three or more different pulsar sources, it is possible to estimate the spacecraft’s position vector in three dimensions in the same fashion as for a classical GNSS. Reasonably, some conditions have to be respected in order to apply this method. At first, the model for predicting the arrival times of the wave-front must be accurately known. This translates into a requirement for the source, which has to be stable in its pulsing so that its behavior can be safely predicted. In a second instance, the source must emit radiation that is strong enough to be detected by a hardware detector onboard a spacecraft. At the same time, this has to be above the background noise in order to allow for a correct detection.

2.2. Pulsar Phase

The pulsar rotation produces a precise interval between pulses that ranges from milliseconds to seconds for each individual star. Following this physical description, a mathematical timing model can be written as a Taylor expansion up to the third order to describe the signal phase evolution in time [18]

$$\varphi \left(t\right)=\varphi \left(t_0\right)+\left(t-t_0\right)f+{\left(t-t_0\right)}^2\frac{\dot{f}}{2}+{\left(t-t_0\right)}^3\frac{\ddot{f}}{6}$$

(2)

featuring the frequency of rotation and its derivatives as f, $\dot{f}$, $\ddot{f}$ with known initial phase $\varphi \left(t_0\right)$ at a reference time $t_0$. With these parameters known, the arrival time of the wave-front can be safely predicted.

In addition, since no pair of neutron stars are formed in exactly the same manner or have the same geometric orientation, the pulse frequency and signal shape produced are unique, identifying signatures for each pulsar. Thus, pulsars can act as natural beacons, or celestial lighthouses, on an intergalactic scale. In addition to that, the signal these stars produce can be observed in the radio, visible, X-ray, and Gamma-ray bands of the electromagnetic spectrum due to the wide energy range of the process leading to the emission of the particles. Observing them in the different bands might have advantages and disadvantages. Variable sources emitting signals in the radio band are certainly a potential candidate that can be used in a navigation system. However, at the radio frequencies that these sources emit (from 100 MHz to a few GHz), antennas with large diameters are required to detect their signals [19]. Due to neighboring sources that emit in radio bands and also due to the low signal intensity of radio pulsars, long integration times are needed to obtain a signal with an acceptable signal-to-noise ratio (SNR), suitable for use in a navigation system [20].

Similar issues exist for pulsars in the visible spectrum. In this case, less than one hundred of isolated pulsars are cataloged, and all of them are faint [21]. Few are also the pulsars discovered that emit in the Gamma-ray wavelengths; studies were conducted regarding this type of pulsar featuring their so-called Gamma-ray bursts (GRBs), which were found to have a very high signal to noise ratio, but due to their cataclysmic nature, they are non-repeating and non-periodic, thus presenting a very difficult challenge for navigation applications [22]. The lower amount of these specific pulsars and their issues related to faintness and instability make them less suitable for navigation purposes.

Luckily enough, pulsars exist emitting mainly in the X-ray band that would require relatively small-sized detectors to be mounted onboard spacecraft, therefore partially easing the system design.

2.3. Signal Model

Due to the poor strength level of the pulsar’s signal with respect to the background, photons arriving in a single period are generally not enough to fully characterize the signal and to precisely detect its peaks; for this reason, the pulsar has to be observed for several periods, time tagging all the incoming photons. Then, from the raw data obtained, the signal can be reconstructed by manipulating the photons’ time history.

Therefore, especially in presence of a low SNR, a fully deterministic representation of the process is not possible, but rather, a stochastic description can be more suitable for this kind of problem. In these conditions, one possible description of this stochastic process is to employ the Non-Homogeneous Poisson model, as it is extensively used for signal models and specifically for X-ray pulsars [23,24].

Given an observation time interval $T_{obs}$ bounded by an initial $t_0$ and a final $t_f$, let us denote $t_i$ the time corresponding to the arrival time of the i-th photon. Let us assume that the sequence ${\left\{t\right\}}_i^N$ is random, and it is expressed in incremental order within the observation window as

$$t_0\le t_1<t_2<\dots t_i\le t_f$$

where the variable N represents the total number of photons received and is itself random. The problem $\{N_{tot},t>0\}$ is said to be a point-process with $t_0=0$ and $N_0=0$ as initial conditions and with $N_{tot}$ as the total number of photons received in the interval $\left(0,t\right)$. A point-process is suitable for a Non-Homogeneous Poisson representation with a time-varying photon arrival rate $\lambda \left(t\right)$.

A process is said to belong to the Non-Homogeneous Poisson class when it satisfies the following three conditions:

(1)

The probability of detecting one photon in a time interval $\Delta t$ is given by:

$P\left(N_{t+\Delta t}-N_t=1\right)=\lambda \left(t\right)\Delta t$    with $\Delta t$ approaching zero.

(2)

The probability of detecting more than one photon in $\Delta t$ is given by:

$P\left(N_{t+\Delta t}-N_t>1\right)=0$    with $\Delta t$ approaching zero.

(3)

Non-overlapping increments are independent, with $N_t$ as the increment of the stochastic process:

$N_{l,q}=N_q-N_l$    with $q\ge l$.

In this framework, the number of received photons k equal to $N_q-N_l$ is a Poisson random variable in any fixed time interval $\left(l,q\right)$, featured by the integrated rate ${\int }_l^t\lambda \left(ϵ\right)dϵ$ and with probability [23]

$$P\left(k\right)=\frac{{\left({\int }_l^q\lambda \left(ϵ\right)\mathrm{d}ϵ\right)}^{k}exp\left(-{\int }_l^q\lambda \left(ϵ\right)\mathrm{d}ϵ\right)}{k!}$$

(3)

Equation (3) is the probability that k photons are received in a time frame from l to q. Now, the incoming rate of photons ($\lambda \left(t\right)$) is constituted by a background source and a signal source. Therefore, we can write

$$\lambda \left(t\right)={\lambda }_b+{\lambda }_{s}h\left({\varphi }_{det}\left(t\right)\right)$$

(4)

where ${\lambda }_b$ and ${\lambda }_s$ are the arrival rates of background and source photons, respectively, in terms of photons per second, and $h\left({\varphi }_{det}\left(t\right)\right)$ is the pulsar signal profile whose shape is different for each neutron star. It is periodic and dependent on the detected phase ${\varphi }_{det}\left(t\right)$, defined in the interval ${\varphi }_{det}\in \left(0,1\right]$, so that $h\left(\varphi +n\right)=h\left(\varphi \right)$, with n integer. It is also defined as normalized and non-negative, implying that ${\int }_0^{1}h\left(\varphi \right)\mathrm{d}\varphi =1$ and $mi{n}_{\varphi }h=0$. As an example, the normalized profile of the B0531+21 pulsar, also known as the Crab pulsar, is reported in Figure 2.

The phase with which the profile is expressed needs to be related to time to obtain an arrival function depending on the time. The phase will also depend on the observed pulsar frequency, which, including Doppler effects, can be written as:

$$f_0\left(t\right)=f_s+f_d\left(t\right)$$

(5)

where $f_s$ is the pulsar base frequency and $f_d$ the contribution due to the Doppler shift.

2.4. Epoch Folding

It is necessary to develop a method for the onboard processing of the information contained in the photons time tags to estimate signal phase at specific times, therefore producing a useful measurement to ingest in a Kalman filter. A feasible method that can be used directly on the onboard computer is the one generally employed by astrophysicists to generate the light curves of celestial objects starting from raw data, i.e., the so-called Epoch Folding. The idea is rather simple and consists of the following steps [25]:

• The photons’ time tags during the set observation window are collected.
• They are folded back into a single time interval equal to one pulse period.
• The period duration is divided into some equal-length bins.
• The number of photons in each bin is counted.
• The computed photon counts are normalized, and the empirical pulsar profile is derived.

Therefore, the estimated discrete rate function can be written as

$$\widehat{\lambda }\left(t_i\right)=\frac{1}{N_p{T}_b}\sum _{j=1}^{N_b}{c}_j\left(t_i\right)$$

(6)

where $t_i$ is the reference signal epoch, $T_b$ is the bin size, $N_b$ the number of bins the single period is divided into, and finally, $c_j\left(t_i\right)$ are Non-Homogeneous Poisson variables, whose mean and variance are described by:

$$E\left[c_j\left(t_i\right)\right]=var\left[c\left(t_i\right)\right]=\lambda \left(t_i\right)T_b$$

(7)

It must be remarked that the Epoch Folding technique relies on the knowledge of the spacecraft velocity with respect to the pulsar position, so that the Doppler frequency shift can be computed and taken into account in the signal reconstruction. However, the largest part of the Doppler shift is already taken into account from the solar system motion with respect to the pulsar, so that the relative motion of the spacecraft is negligible. This folding process is graphically depicted in Figure 3.

An alternative to these methods, which requires the knowledge of the rate function obtained out of the Epoch Folding, is the Maximum Likelihood Estimation (MLE) based on the probability density function of the photon time of arrivals [26]. As a reference, in [27], the performances of the MLE are studied with respect to the Cramer Rao Lower Bound, which measures the efficiency of the estimators based on their variance.

2.5. Phase Delay Estimation

The core of the navigation system is the measurement obtained through estimation of the time delay between two received signals. Therefore, a spacecraft in deep space is acquiring the pulsar signal, which is processed via epoch folding, and then the empirical rate function $\widehat{\lambda }\left(t_i\right)$ is obtained, and it will be used in the problem of phase estimation.

One possible approach for the estimation of the initial phase is the solution of a Non-linear Least Squares (NLS) optimization problem, constituted by the fit of the empirical rate function $\widehat{\lambda }\left(t_i\right)$ to the true known rate function $\lambda \left(t_i\right)$. The cost function can then be explicitly written as:

$$J_{\varphi }=\sum _{j=1}^{N_b}{\left(\widehat{\lambda }\left(t_i\right)-\lambda \left(t_i,\varphi \right)\right)}^2$$

(8)

for which the true known rate function is expressed in the unknown variable $\varphi$

$$\lambda \left(t,\varphi \right)={\lambda }_b+{\lambda }_{s}h\left(\varphi +\left(t-t_0)f_0\right)\right)$$

(9)

A solution is then sought by minimizing (8)

$$\widehat{\varphi }=\underset{\varphi \in \left(0,1\right)}{argmin}{J}_{\varphi }$$

(10)

It can be shown [23,25] that such an estimator is asymptotically unbiased, with a specified asymptotic variance.

3. Measurement Modeling

Accurately timing pulsar signals at the spacecraft’s location is crucial to enable their evaluation with respect to the SSB reference. As an example, to provide an accurate position information on the order of 300 m or lower, the system must accurately time pulses to at least $1\mu s$. Therefore, it is then mandatory to properly model the measurement model, which was previously introduced as the right-hand side of Equation (1).

3.1. Photon Path

The higher-order terms featured in Equation (1) come from the field of general relativity, which deals with how the signal travels when subjected to gravity fields. The path of a photon will follow the curvature of space-time when passing near a star; therefore, the photon will not necessarily travel in a straight line as (1) describes.

Several works described this process [28,29] and, by assuming a constant pulsar direction throughout the Solar System and not varying in time, a more complete expression than (1) can be found [28]:

$$\begin{array}{cc}\hfill t_{s/c}-t_{ssb}=& \frac{1}{c}\left\{\mathbf{n}·\mathbf{r}-\frac{r^2}{2D_0}+\frac{{\left(\mathbf{n}·\mathbf{r}\right)}^2}{2D_0}-\frac{\left(\mathbf{b}·\mathbf{r}\right)}{D_0}+\frac{\left(\mathbf{n}·\mathbf{b}\right)\left(\mathbf{n}·\mathbf{r}\right)}{D_0}\right\}+\hfill \\ \hfill & +\frac{2G{M}_{sun}}{c^3}ln\left|\frac{\mathbf{n}·{\mathbf{r}}_{\mathbf{sun}}+r_{sun}}{\mathbf{n}·\mathbf{b}+b}+1\right|\hfill \end{array}$$

(11)

for which $D_0$ represents the pulsar distance with respect to the Solar System, $\mathbf{r}$ is the spacecraft position vector, $\mathbf{b}$ is the position vector of the Sun with respect to the SSB, and ${\mathbf{r}}_{\mathbf{sun}}$ is the spacecraft position vector with respect the Sun.

3.2. Timing Model

For a correct measure, the time when the pulse is detected at the spacecraft’s location must be translated into the TCB frame, decoupling its flow from gravity and velocity feature. Additionally, this time is read from an on-board reference, which in general is not absolute, introducing one more issue in the navigation solution estimation.

Clocks measure time with respect to an oscillating phenomenon; specifically, atomic clocks observe the microwave signal emitted by electrons when they change energy levels for increased accuracy. If the frequency of the phenomenon is known, by observing its evolution the elapsed time can be inferred. In the ideal case, the time kept should be infinitely accurate, but both the physical phenomenon and the hardware built to observe it introduce several uncertainties which are usually collected under the definition of clock stability. The accuracy specification of a clock describes how much deviation there can be between the specified clock frequency and the actual one. The stability specification gives a measure of how much the frequency varies over time, introducing uncertainties in reading it.

3.3. Time Conversion

As stressed, the time conversion dynamics has to be considered to align the local time system of the spacecraft with respect to the proper one. The state space model consisting of the spacecraft’s three-dimensional position and velocity needs to be augmented with the time conversion equation, adding this new variable as a state. The time conversion follows the Lorentz transformation:

$$\mathrm{d}\tau =\left[1-\sum _{i=1}^{N_{ss}}\frac{G{M}_i}{|{\mathbf{r}}_i|}\frac{1}{c^2}-{\left(\frac{v_{ssb}}{c}\right)}^2\right]\mathrm{d}t$$

(12)

for which t represents the onboard coordinate time and $\tau$ the proper one, G is the gravitational constant, $M_i$ the mass of the i-th body, and ${\mathbf{r}}_i$ its relative distance from the spacecraft. Additionally, $N_{ss}$ is the number of planetary bodies considered, and finally, $v_{ssb}$ is the velocity with respect to an inertial SSB-centred frame. According to [29], (12) is accurate up to ${10}^{-12}$ s. Therefore, (12) can be rewritten in differential form

$$\delta \dot{\tau }=\frac{\mathrm{d}t-\mathrm{d}\tau }{\mathrm{d}t}={\left(\frac{v_{ssb}}{c}\right)}^2+\sum _{i=1}^{N_{ss}}\frac{G{M}_i}{|{\mathbf{r}}_i|}\frac{1}{c^2}$$

(13)

where $\delta {\dot{\tau }}_{TDB}$ is the time flow drift between the two time systems. Therefore, the state is then augmented as:

$$\mathbf{x}\left(t\right)={\left[\mathbf{r},\mathbf{v},\delta \tau \right]}^T$$

(14)

Now, the measurement model of Equation (11) can be rewritten as

$$\begin{array}{cc}\hfill h\left(\mathbf{x}\left(t\right),t\right)=& \frac{1}{c}\left\{\mathbf{n}·\mathbf{r}-\frac{r^2}{2D_0}+\frac{{\left(\mathbf{n}·\mathbf{r}\right)}^2}{2D_0}-\frac{\left(\mathbf{b}·\mathbf{r}\right)}{D_0}+\frac{\left(\mathbf{n}·\mathbf{b}\right)\left(\mathbf{n}·\mathbf{r}\right)}{D_0}\right\}+\hfill \\ & +\frac{2G{M}_{sun}}{c^3}ln|\frac{\mathbf{n}·{\mathbf{r}}_{\mathbf{sun}}+r_{sun}}{\mathbf{n}·\mathbf{b}+b}+1|+\delta \tau \hfill \end{array}$$

(15)

3.4. Atomic Clocks

The non-ideal timekeeping of the on-board clocks can be taken into account, where a stochastic description can be used for this problem [30]. The clock errors can be modelled by stochastic differential Equations (SDE). A two-state clock model can be written by introducing two new states, the clock bias $x_b$ and the clock drift $x_d$, which are related as:

$$\left\{\begin{array}{c}\hfill {\dot{x}}_b=x_d+{\sigma }_{b}d{W}_b\\ \hfill {\dot{x}}_d=a_d+{\sigma }_{d}d{W}_d\end{array}\right.$$

(16)

for which $d{W}_b$ and $d{W}_d$ are the Wiener processes associated with the drift and the bias, while ${\sigma }_b$ and ${\sigma }_d$ are their associated diffusion coefficients. Finally, $a_d$ is the frequency aging component. Apart from the Wiener processes which make the equations stochastic, the other features can be assumed constant and set according to the specific hardware considered. Now, the state (14) can be augmented once again as:

$$\mathbf{x}\left(t\right)={\left[\mathbf{r},\mathbf{v},\delta \tau ,x_b,x_d\right]}^T$$

(17)

and the measurement model updated accordingly:

$$\begin{array}{cc}\hfill h\left(\mathbf{x}\left(t\right),t\right)=& \frac{1}{c}\left\{\mathbf{n}·\mathbf{r}-\frac{r^2}{2D_0}+\frac{{\left(\mathbf{n}·\mathbf{r}\right)}^2}{2D_0}-\frac{\left(\mathbf{b}·\mathbf{r}\right)}{D_0}+\frac{\left(\mathbf{n}·\mathbf{b}\right)\left(\mathbf{n}·\mathbf{r}\right)}{D_0}\right\}+\hfill \\ & +\frac{2G{M}_{sun}}{c^3}ln|\frac{\mathbf{n}·{\mathbf{r}}_{\mathbf{sun}}+r_{sun}}{\mathbf{n}·\mathbf{b}+b}+1|+\delta {\tau }_{TDB}+x_b\hfill \end{array}$$

(18)

Equation (18) is the distance of the spacecraft with respect to the SSB along the pulsar direction $\mathbf{n}$, already accounting for relativistic effects and clock errors. This equation constitutes the measurement model that will be used in the navigation filter in Section 4. Including these aspects in the State Estimation becomes especially critical when CubeSat systems are considered as usually COTS components with predefined performances are selected for implementation. This augmented model provides the capability of simulating the impact of the clock selection on the final navigation results.

3.5. Signal-to-Noise Ratio

Additionally, this work includes the effect of the signal-to-noise ratio on the measurement accuracy. As shown in [22], the SNR measure of a pulsar signal is:

$$SNR=\frac{F_x{A}_{det}{p}_f{t}_{obs}}{\sqrt{\left[B_x+F_x\left(1-p_f\right)\right]\left(A_{det}{t}_{obs}\frac{W}{P}\right)+F_x{A}_{det}{p}_f{t}_{obs}}}$$

(19)

where $F_X$ is the X-ray photon flux, $A_{det}$ is the instrument collecting area, $p_f$ is the pulsed fraction of the source flux, $t_{obs}$ is the observation time, $B_x$ is the background photons flux, W is the pulse width, and P is the pulse period [18]. A direct first-order expression for the uncertainty of the pulsar measurements time-of-arrival, ${\sigma }_{TOA}$, can be obtained in relation to (19), it being [22]:

$${\sigma }_{TOA}=\frac{1}{2}\frac{W}{SNR}$$

(20)

which can be converted into a range uncertainty

$${\sigma }_{range}=c{\sigma }_{TOA}$$

(21)

The signal-to-noise modelling will be used in the two devised application scenarios, highlighting its impact on different spacecraft characteristics.

4. Simulation

The state of the spacecraft, during its deep-space travel, is known with some uncertainty and, if not updated or corrected continuously, tends to diverge from the reality. In the case of X-ray pulsars navigation, both the acquisition and processing of the data are affected by errors and uncertainties which introduce issues in the state estimation and possibly leading to performance degradation. This works considers the application of the unscented Kalman filter (UKF) [31] to retain all the non-linear features of the developed models. This is accomplished by representing the set of states by their mean and covariance, and these features are then propagated through the actual, non-linear model [31]. The devised methodology is applied to two different mission scenarios, one being the New Horizons mission case, and the second being the LUMIO mission case. The two missions differ mainly in their orbits and in the available platform budget characteristics, which in turn have an impact on the achievable SNR of the pulsar signals.

4.1. New Horizons Mission Case

The navigation model just built is first applied to the New Horizon mission as it represents one of the most iconic and successful deep space missions [32]. The orbit considered can be described as a very stretched arc of hyperbola since the spacecraft is moving at a velocity higher than the escape velocity of the Solar System. In particular, the trajectory leg just before the approach of Pluto is considered. In order to accomplish this, data were collected from the Jet Propulsion Laboratory’s (JPL) Horizon database (https://naif.jpl.nasa.gov/pub/naif/pds/data/nh-j_p_ss-spice-6-v1.0/nhsp_1000/data/spk/ (accessed on 12 January 2022)).

The kernel nh_recon_pluto_od122_v01 containing the data required for the simulation was then retrieved, and the initial state is the one at midday of 7 December 2014. This has been also reported in

Table 1. New Horizon initial state in J2000, Sun centered.

$x_0$	1,092,587,085.357 km	$v_{x,0}$	5.560 km/s
$y_0$	−4,206,812,531.514 km	$v_{y,0}$	−12.623 km/s
$z_0$	−1,655,791,253.618 km	$v_{z,0}$	−4.892 km/s

.

The simulation parameters are reported in

Table 2. Simulation parameters.

Parameter	Value	Unit	Parameter	Value	Unit
$T_{obs}$	3600	s	${\sigma }_{pos}$	${10}^2$	km
$T_{map}$	0.5	d	${\sigma }_{vel}$	${10}^{-4}$	km/s
$B_x$	0.005	ph/cm2/s	${\sigma }_{TDB}$	${10}^{-6}$	km
$N_{psr}$	4		${\sigma }_{bias}$	${10}^{-9}$	s
$A_{det}$	1	m2	${\sigma }_{bias}$	${10}^{-16}$	1/s
$\Delta t_{eul}$	10	s	$a_d$	${10}^{-14}$	s/s

.

In this scenario, four pulsars (see

Table 3. Selected pulsars’ parameters. Data are taken from [33] and updated where applicable with more recent data.

Pulsar	${\mathit{F}}_{\mathit{x}}$	Galactic Latitude	Galactic Longitude	Period
J1751−305	0.180000 ph/cm2/s	−0.0330 deg	6.27 deg	2.30 ms
B0531+21	1.540000 ph/cm2/s	−0.1000 deg	3.22 deg	33.4 ms
B1937+21	0.000050 ph/cm2/s	−0.0051 deg	1.00 deg	1.60 ms
B1821−24	0.000193 ph/cm2/s	−0.0970 deg	0.14 deg	3.10 ms

) are tracked sequentially, with measures acquired every $T_{map}$, which is half a day, with an observation window $T_{obs}$. The measurements fed to the UKF are generated by corrupting the true measure of (18) with Equation (21). The initial state of the filter is given by a distribution represented by the $\sigma$ in Table 2. The true initial state regarding position and velocity is the one in Table 1, the time conversion and clock dynamics are initialized as a vector of zeros.

With these settings, the UKF formulation can be run and the navigation accuracy can be evaluated with respect to the true states. The accuracy for the first six states, which are the position and velocity components, are presented in Figure 4, both for the residual and the 1$\sigma$ knowledge. As it can be seen, by tracking a different pulsar each half day, the estimation accuracy increases, and the filter manages to produce a suitable navigation solution even if starting with a high dispersion. The values of the position and velocity accuracy at the final time are reported in

Table 4. New Horizon final time navigation results.

${\sigma }_x$	4.532 km	${\sigma }_{v_x}$	0.8548 mm/s
${\sigma }_y$	0.751 km	${\sigma }_{v_y}$	0.1561 mm/s
${\sigma }_z$	2.022 km	${\sigma }_{v_z}$	0.4056 mm/s

.

The accuracies of the time conversion and clock states are displayed in Figure 5 and Figure 6, respectively. As it can be seen, the filter manages to keep the clock instability within bounds, in around 0.1 microseconds, showing that combined navigation and time keeping features can be obtained. The time conversion knowledge is also kept at acceptable levels with a distance-equivalent knowledge below 1 m; in this case, the good performances are the results of the combination of the applied filtering on the $\delta \dot{\tau }$ state and the fact that, it being highly correlated with the first orbital states, once they reach convergence, it tends not to diverge. The filter manages to reduce the knowledge to less than 10 km in 2–3 update steps.

The navigation results in Table 4 show that performances for position estimation with accuracy in the order of kilometers can be obtained with the setup of Table 2 which, taking into account the scale of the problem (more than 30 astronomical units), can be considered a very good result in view of similar mission scenario implementations for the future.

4.2. LUMIO Mission Case

The second case considered is the operative orbit of the LUMIO mission, which is a deep-space CubeSat under development by the European Space Agency (ESA) [34,35]. The scientific objectives of the LUMIO mission are to observe, quantify, and characterize meteoroid impacts on the lunar far-side by detecting their impact flashes, complementing Earth-based observations on the Lunar nearside, to provide global information on the Lunar Meteoroid Environment and contribute to Lunar Situational Awareness. LUMIO flies on a Halo orbit at the Earth–Moon L2 point. A window of 60 days is considered for this simulation, starting from the first kernel entry. The initial state is reported in

Table 5. LUMIO initial state vector in EME2000.

Start Date	29 August 2020
$x_0$	211,878.048 km	$v_{x,0}$	1.187 km/s
$y_0$	−345,283.782 km	$v_{y,0}$	0.691 km/s
$z_0$	−136,298.530 km	$v_{z,0}$	0.221 km/s

.

The simulation settings are shown in

Table 6. LUMIO simulation setup.

Parameter	Value	Unit	Parameter	Value	Unit
$T_{obs}$	300	s	${\sigma }_{pos}$	${10}^2$	km
$T_{map}$	1	d	${\sigma }_{vel}$	${10}^{-4}$	km/s
${\lambda }_{back}$	0.005	ph/cm2/s	$N_{bin}$	2048
$N_{psr}$	1		$A_{det}$	0.01	m2

. In this case, a single pulsar tracking will be adopted, with the Crab pulsar as the chosen option. Additionally, both the clock dynamics and the time conversion of the photon time tags are assumed to be perfectly known. Doppler effects are then included both in the photon generation and in their relative signal detection. Epoch Folding was performed with 2048 bins, which corresponds to a resolution of roughly 16 microseconds, a resolution that is compatible with a wide range of CMOS detectors. The collecting area is assumed to be 0.01 square meters in order to make it suitable for a CubeSat application (e.g., 12U application, with placement on one of the 4U faces, offering up to 0.04 square meters).

The accuracy of the X-ray pulsar navigation applied to the LUMIO mission scenario is shown in Figure 7. The covariance still achieves a good accuracy, and

Table 7. LUMIO final time navigation results.

${\sigma }_x$	2.595 km	${\sigma }_{v,x}$	1.484 $\times {10}^{-2}$ m/s
${\sigma }_y$	0.804 km	${\sigma }_{v,y}$	6.955$\times {10}^{-3}$ m/s
${\sigma }_z$	1.753 km	${\sigma }_{v,z}$	4.196$\times {10}^{-3}$ m/s

reports the final values in terms of $1\sigma$ knowledge for position and velocity components. The total uncertainty for positioning is in the order of kilometres, while the accuracy for the velocity is in the order of the millimetres per second. Both values are very promising considering deep-space satellite applications, where these uncertainties are completely acceptable with respect to the large distances travelled in deep space. Note that the convergence time is higher with respect to the New Horizon mission case; this is due to the smaller sensor and lower performances arising from the adoption of miniaturized components.

5. Conclusions

This work has investigated the application of pulsar-based navigation in two mission scenarios, highlighting its suitability for nanosatellite applications. A novel approach for dealing with the on-board timing features has been addressed, in conjunction with the pulsar signal modelling and onboard filtering implementations. All in all, reachable navigation accuracies are compliant with typical requirements of deep-space missions, making X-ray pulsar-based navigation promising for future onboard applications.