Robust Onboard Orbit Determination Through Error Kalman Filtering

Abstract

Accurate and robust on-board orbit determination is essential for enabling autonomous spacecraft operations, particularly in scenarios where ground control is limited or unavailable. This paper presents a novel method for achieving robust on-board orbit determination by integrating a loosely coupled GNSS/INS architecture with an on-board orbit propagator through error Kalman filtering. This method is designed to continuously estimate and propagate a spacecraft’s orbital state, leveraging real-time sensor measurements from a global navigation satellite system (GNSS) receiver and an inertial navigation system (INS). The key advantage of the proposed approach lies in its ability to maintain orbit determination integrity even during GNSS signal outages or sensor failures. During such events, the on-board orbit propagator seamlessly continues to predict the spacecraft’s trajectory using the last known state information and the error estimates from the Kalman filter, which were adapted here to handle synthetic propagated measurements. The effectiveness and robustness of the method are demonstrated through comprehensive simulation studies under various operational scenarios, including simulated GNSS signal interruptions and sensor anomalies.

1. Introduction

In the realm of space engineering and satellite missions, the drive towards spacecraft autonomy has become increasingly vital. Autonomous operations enable spacecraft to perform critical tasks and decision making independently, reducing reliance on continuous communication with Earth and enabling missions with challenging goals where real-time ground control is limited or impractical. One fundamental aspect of spacecraft autonomy is the ability to accurately determine and predict the spacecraft’s orbital trajectory onboard, particularly in scenarios where the precise estimation of the orbital states is mandatory to achieve mission requirements, such as in proximity operations, on-orbit inspection and servicing, sub-millimetre Earth monitoring, large constellation control, precise station keeping and continuous orbit maintenance. The motivation behind enhancing spacecraft autonomy stems from the complexity and risks associated with relying solely on ground-based operations. Communication delays, limited coverage in certain regions of space, and high costs in operating large ground infrastructure highlight the necessity of developing self-sufficient spacecraft capable of navigating and adapting to dynamic environments autonomously.

This paper addresses the critical need for robust onboard orbit determination by proposing a novel approach which integrates a global navigation satellite system (GNSS) and an inertial navigation system (INS) in a loosely coupled architecture with an onboard orbit propagator, utilizing error-state Kalman filtering techniques. This integrated method leverages real-time sensor measurements from the GNSS receiver and the INS to accurately estimate the spacecraft’s orbital state and the continuously computed dynamical solutions from the orbital propagator to feed the error-state Kalman filter and guarantee the robustness of the proposed orbit determination system (ODS). In fact, by integrating error-state Kalman filtering within the GNSS/INS and orbit propagation framework, accurate and reliable orbit determination is possible, even during GNSS signal outages or sensor failures. This capability is essential for maintaining mission continuity and enabling spacecraft to autonomously adapt to unforeseen circumstances without immediate ground intervention. In case the failure is permanent or extremely prolonged in time, the unavoidable drift in the estimated orbital states can be reduced or fixed with periodic updates from ground in the form of two-line element (TLE) states.

The effectiveness and robustness of the proposed method are validated through comprehensive simulation studies under various operational scenarios, including simulated GNSS signal interruptions and sensor anomalies. The model-in-the-loop (MIL) simulations exploit a validated functional engineering simulator (FES) with hardware-calibrated sensor models [1]. Their results demonstrate the method’s capability to enhance mission reliability by enabling spacecraft to autonomously determine and navigate their orbital states with a high degree of accuracy and resilience. Existing literature studies have addressed the problem of onboard orbit determination for many years, exploiting various architectures and proposing multiple alternative methods.

The integration of a GNSS and INS is classically exploited as an integrated navigation system designed for onboard space applications. For many years, the market’s available systems have exploited these integrated instruments to autonomously provide an orbit determination solution. Initially deployed on the International Space Station (ISS) for the Crew Return Vehicle (CRV), the onboard orbit determination systems had to meet navigation needs across all flight phases and adhere to size, weight, computing and power constraints. Notably, these integrated systems commonly utilize a loosely coupled architecture for integrating a GNSS and INS, and thus the performance of the system is inherently associated with those of the GNSS receiver. A tightly coupled GNSS/INS integration could overcome this constraint by directly accessing the GNSS’s raw observables, even if the complexity of the overall navigation system increases. In these regards, significant advancements have been made in integrated GNSS/INS systems for aircraft and missiles in particular, with many devices currently in operational use. For instance, the Honeywell H764G Embedded GNSS/INS, widely used in various military aircraft, has demonstrated its precision in approach and landing operations, as highlighted by Moya and Elchynski [2].

In the domain of spacecraft navigation, Ebinuma [3] has demonstrated the application of a GNSS for rendezvous missions. He developed an extended Kalman filter to facilitate real-time relative navigation, utilizing a hardware-in-the-loop test facility which integrates navigation and guidance systems for rendezvousing. However, his models presupposed the continuous availability of GNSS signals without accounting for multipath signals or blockages [3]. Um’s recent research has focused on navigation during proximity operations, particularly examining spacecraft relative navigation using an integrated GNSS/INS system near the ISS [4]. The system employed a loosely coupled architecture with two Kalman filters: one for processing GNSS data and another for merging INS data with the GNSS filter outputs. The INS component was simulated via software, while the GNSS data were sourced from a Mitel Architect GNSS receiver stimulated by a GSSI STR4760 GNSS simulator (Spirent Communications: Crawley, West Sussex, UK).

Considering recent studies enhancing the precision and reliability of on-orbit orbit determination techniques, Mao demonstrated a highly effective method for GPS-based low Earth orbit (LEO) determination, achieving 1 cm precision using advanced GNSS processing techniques which significantly enhance orbit calculation accuracy for LEO satellites [5]. Extending this work, a later work of Mao explored GNSS-based precise orbit determination for maneuvering LEO satellites, showing that these techniques can adapt to dynamic conditions and maneuvers, thus improving the accuracy of orbit determination for satellites in complex trajectories [6]. These studies established a strong benchmark for accuracy relative to the current state of the art, as comprehensively reviewed by Selvan, who listed precise orbit determination techniques applied to various LEO satellites, emphasizing both progress and challenges in achieving high accuracy [7].

Different research works also considered the possibility of integrating ground-based and onboard navigation solutions. Cho focused on NORAD TLE sets for orbit determination, integrating GNSS navigation solutions to enhance accuracy, although traditional TLE methods face limitations in terms of update frequency and precision compared with modern GNSS approaches [8]. Similarly, Coffee examined the accuracy and update frequency of NORAD TLE sets for spacecraft in LEO, noting that while TLE sets provide useful data, they often lack the precision required for high-fidelity onboard orbit propagation [9]. Conversely, Park verified the navigation performance of space GNSS receivers using full-orbit solutions for satellites, demonstrating the reliability of modern GNSS technology across entire orbital periods [10]. Thus, a positive integration of onboard GNSS/INS navigation methods with onboard orbit solutions and ground-based TLE can be sought in the context of non-nominal degraded scenarios to guarantee continuity of operations rather than precision enhancement only. Indeed, in case of GNSS outage periods, a GNSS/INS system integrated with an onboard orbital propagator could enhance the reliability of the navigation solution in both nominal and failure scenarios. The last valid GNSS fix can be used to propagate the onboard orbit determination with acceptable accuracy, keeping its integration with the INS and the navigation functions. In case prolonged failures occur, the ground solutions in the form of TLE data can be exploited to keep the system away from large drifts and errors.

An interesting research point deals with the performance of GNSS solutions in various orbital regimes, including medium Earth orbit (MEO), geostationary orbit (GEO), and highly elliptical orbit (HEO). Capuano and Liu provided insights into GNSS performance across these different regimes, highlighting the need for tailored approaches to meet the specific requirements of each orbital environment [11,12]. A later Capuano work further explored high-accuracy GNSS navigation in GEO, showcasing advancements which enable precise navigation even in geostationary orbits [13]. Additionally, Guan investigated the use of multi-GNSS sidelobe signals for navigation in GEO, HEO, and lunar trajectories, indicating that leveraging sidelobe signals can improve navigation performance in challenging orbital environments [14]. This study highlights ongoing innovations in GNSS technologies and their applications to diverse space missions with different possible accuracy levels. In MEO, GNSSs can achieve accuracies ranging from a few meters to sub-meter levels. In HEO, accuracy can be less stable due to varying satellite visibility and signal strength, generally ranging from several meters to tens of meters. For GEO, the accuracy is more challenging to maintain because of the limited visibility of GNSS satellites, typically resulting in accuracies from tens of meters to several hundred meters. The precision in GEO is affected by factors such as geometric dilution of precision (GDOP) and signal attenuation, though combining a GNSS with other navigation aids can improve performance. Finally, in lunar trajectories, the accuracy level range us in the order of kilometers, but there exist possibilities for improvement.

In this paper, a positive integration of a loosely coupled GNSS/INS onboard orbit navigation system with an onboard propagator capable of being updated and resynchronized by means of ground solutions in TLE format is proposed. The accuracy of the navigation solution is limited by the one for the GNSS receiver, as is typical for such loosely coupled architectures. However, the system has proved to be robust to short GNSS data outages. Moreover, the capability to accept ground-based navigation solutions allows keeping reliable orbit knowledge onboard, even after a complete GNSS receiver failure. In this non-nominal case, the accuracy is limited by the ground data but can be acceptable for absolute orbit navigation requirements.

2. Mathematical Background

The integration of a GNSS with an INS enhances the quality and integrity of both navigation systems. The use of a GNSS allows for the calibration of inertial instrument biases, while the INS can improve the tracking and reacquisition performance of the GNSS receiver. Two error calibration techniques can be applied in an integrated GNSS/INS system—the feedforward (or open-loop) method and the feedback (or closed-loop) method—as illustrated in Figure 1. The feedforward method calculates estimated navigation errors and directly applies them to the navigation output. However, this approach has a significant drawback; the INS remains unaware of the external aid, which can lead to unbounded error growth. If this error becomes too large, it can push the system beyond its linear operating range, compromising the effectiveness and accuracy of the linear Kalman filter. In contrast, the feedback method addresses this issue by feeding the estimated errors back into the INS mechanization to correct the equipment’s calculations.

Furthermore, there are two primary methods for integrating GNSS and INS data: loosely coupled and tightly coupled. In a loosely coupled integration [15], a navigation processor within the GNSS receiver calculates the position and velocity using GNSS observables alone. An external navigation filter then computes the position, velocity and attitude from the raw inertial sensor measurements, using the GNSS-derived position and velocity to calibrate INS errors. One advantage of a loosely coupled system, depicted in Figure 2, is that the GNSS receiver can be treated as a stand-alone processed unit. The blended navigation filter is simpler when utilizing GNSS preprocessed position and velocity measurements. However, during a GNSS outage, the GNSS ceases to provide processed measurements, and the calibration of inertial sensors from the GNSS/INS filter halts.

A more complex GNSS/INS filter design can mitigate issues related to GNSS satellite blockage. Figure 2 illustrates a tightly coupled GNSS/INS integration. In this configuration, the external navigation filter receives raw GNSS measurements of a pseudo-range and Doppler or delta range. The tightly coupled GNSS/INS filter can utilize GNSS measurement updates even when fewer than four satellites are available, which is insufficient for a complete GNSS navigation solution. In practice, a third technique often called deep or ultra-tight coupling exists, in which the INS and GNSS receiver are combined at the signal processing level [16]. In particular, the filter pseudo-range and Doppler predictions are injected into the GNSS receiver carrier phase and code tracking algorithms to aid the processing of the incoming signals, allowing for acquisition of weaker signals and a faster position fix.

An INS can enhance the performance of a GNSS receiver on multiple levels. By providing outputs such as the position, velocity and attitude, the INS can serve as external inputs to a GNSS receiver, facilitating faster signal acquisition through pre-positioning calculations and aiding in interference rejection during signal tracking. Navigation systems are inherently nonlinear dynamic systems, and a fundamental aspect of navigation system development is the estimation of the states of such systems. In the context of state estimation for nonlinear systems, no single approach universally outperforms all others.

2.1. GNSS Receiver Model

The GNSS satellite ephemerides broadcasted in the navigation message express the satellite antenna phase center in a dedicated Earth-centred, Earth-fixed (ECEF) reference frame. For example, the United States GPS uses the World Geodetic System WGS-84, Galileo satellites broadcast their ephemerides in the Galileo Reference Frame (GTRF), BeiDou adopts the CGCS2000, and Glonass uses the PZ-90 reference frame. Nevertheless, all of the latest realisations of these reference frames are either closely aligned to the International Terrestrial Reference Frame (ITRF) or provide corrections in the form of a set of translation and rotation parameters, with maximum deviations from the ITRF to the order of a few centimeters. For this reason, throughout this work, these minor residual differences are neglected. However, since the orbital state estimated by the filter is expressed in an inertial frame, the navigation solution provided by the GNSS receiver must undergo a post-processing step which rotates it from the ITRF to the Geocentric Celestial Reference Frame (GCRF).

The latest transformation procedures between the ICRF and the GCRF are defined by the International Earth Rotation and Reference System Service (IERS) in the 2010 conventions [17]. In particular, according to the IAU 2006/2000A CIO-based procedure, the rotation matrix from the ITRF to the GCRF and its derivative are expressed as a product of three fundamental components:

$$\begin{array}{cc}\hfill {\mathbf{r}}_{GCRF}& =\mathbf{Q}(t)\mathbf{R}(t)\mathbf{W}(t){\mathbf{r}}_{ITRF}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\mathbf{v}}_{GCRF}& =\mathbf{Q}(t)\mathbf{R}(t)\left[\mathbf{W}(t){\mathbf{v}}_{ITRF}+{\omega }_e\times \mathbf{W}(t){\mathbf{r}}_{ITRF}\right]\hfill \end{array}$$

(2)

where $\mathbf{Q}(t)$, $\mathbf{R}(t)$ and $\mathbf{W}(t)$ are the rotation matrices which account for the motion of the Celestial Intermediate Pole (CIP) in the GCRF, the rotation of the Earth around the CIP and the polar motion, respectively. In the velocity transformation, ${\omega }_e$ is the rotation rate of the Earth, and the time derivatives of the long-term nutation and precession effects are neglected. Meanwhile, $\mathbf{R}(t)$ and $\mathbf{W}(t)$ have relatively simple expressions. The precise computation of the transformation matrix for the celestial motion of the CIP $\mathbf{Q}(t)$ is quite cumbersome as it requires the evaluation of lengthy Poisson series which express the $X,Y$ coordinates of the CIP in the GCRS with microarcsecond accuracies. Additionally, since detailed VLBI observations have shown that there are deficiencies to the order of 0.2 mas in the IAU 2006/2000A precession-nutation model, further time-dependent offsets $\delta X$, $\delta Y$ must be added to the nominal X and Y values. These corrections, together with the $\Delta UT1$ offset, the excess length of day (LOD) and the polar coordinates of the CIP in the ITRF required for the computation of $\mathbf{W}(t)$ (${\omega }_e$ and $\mathbf{R}(t)$, respectively) form the Earth orientation parameters (EOPs) and are periodically published online in the IERS Bulletins.

In practice, simplified expressions with different degrees of accuracy have been derived to reduce the computational time required for the computation of the GCRF-to-ITRF rotation matrix [18]:

• IAU2000B: This retains about 80 coefficients in the nutation series and provides the CIP coordinates with a worst-case accuracy of 1 mas with respect to the IAU 2006/2000A model during the 1995–2050 period.
• CPN_c: This delivers a worst-case accuracy of about 16 mas throughout 1995–2050, leaving only 42 coefficients in the $X,Y$ series.
• CPN_d: By keeping only four coefficients of the entire $X,Y$ series, a 1 arcsec worst-case accuracy can be achieved between 1995 and 2050. Additionally, at this level of accuracy, both the CIO locator s as well as the polar motion matrix can be neglected.

Given the resulting accuracies, none of these simplified models require the EOP corrections to the CIP coordinates ($\delta X$, $\delta Y$). The resulting rotation errors with respect to the IAU 2006/2000A model for a spacecraft in LEO and GEO are reported in Figure 3.

In general, despite the used formulation, each contribution in Equation (1) plays a crucial role in accurately performing the post-processing step to convert GNSS readings from the ITRF to the GCRF. Neglecting any reference frame rotation contribution can lead to significant errors in the position and velocity estimates. The most significant term is obviously the one due to the Earth’s rotation around its axis, expressed by the matrix $\mathbf{R}(t)$. This rotation is so fast that an error of a few seconds can lead to intrinsic conversion mistakes which may completely compromise the whole navigation solution. For instance, neglecting the $\Delta UT1$ term in computing the current sidereal time completely jeopardizes the onboard orbit estimation. However, performing the ITRF-to-GCRF conversion by simply accounting for the Earth’s rotation, even with the $\Delta UT1$ knowledge, is not enough for an accurate onboard navigation system. Indeed, the simple reference frame rotation around the Earth’s axis without any other EOP corrections can introduce errors of several kilometers and tens of meters per second in LEO. Incorporating at least precession and nutation corrections, represented by the matrix $\mathbf{Q}(t)$, significantly improves accuracy, reducing errors to a few meters and a few tenths of a meter per second. For instance, accounting only for precession can lead to positional errors of 200 m and velocity errors of 0.25 m per second in LEO. Incorporating both precession and nutation further reduces these errors in LEO to a few meters and 0.1 m per second, as also evident from the analyses reported in Figure 3, where the CPN_d formulation does not consider polar motion corrections. These numerical values are presented as relative errors in comparison to the exact full mathematical conversion. Thus, the overall navigation error must also account for the inherent errors in the rotation formulation. If continuous updates of the EOP parameters from the IERS Bulletin are not feasible onboard, then they can be linearly interpolated over time. In this case, the relative error of the complete mathematical rotation with interpolated data is approximately 2 m and 0.002 m per second compared with using the actual tabular data.

While onboard updates of tabular EOP can further enhance accuracy, the interpolation error is often comparable to the GNSS receiver accuracy, making interpolated EOP a practical choice. Linear interpolation is particularly effective for parameters like $X,Y$, the LOD and the nutation coefficients. However, for $\Delta UT1$, which exhibits a sawtooth profile, predictions based on linear trends between discontinuities are necessary. The optimal frequency for updating $\Delta UT1$ remains related to the natural evolution of the universal time. In conclusion, accurate EOP corrections are essential for ensuring the reliability of GNSS-based applications and maintaining the integrity of reference frame conversions.

However, these are not the only factors which require attention. Time and data precision also impact the intrinsic error in GNSS measurement post-processing. Reference frame conversions depend on accurate timing, as well as the use of EOP data, whether interpolated or tabular. For example, time conversion errors of less than 1 s can cause positional errors of up to 100 m. If GPS time is confused with Coordinated Universal Time (e.g., a 20 s error at time of writing), then positional errors can increase to tens of kilometers. Similarly, using single-precision calculations for reference frame conversions can introduce numerical errors of up to tens of kilometers in position and tens of meters per second in velocity. Therefore, GNSS data post-processing requires at least double-precision floating point arithmetic to ensure accurate results.

With this in mind, given the results presented in Figure 3, we opted to leverage the simplified CPNc model to convert the GNSS state to the GCRF within our GNSS post-processing algorithm. Indeed, throughout the 2000–2050 timespan, the maximum error caused by the truncation of the CIP series is limited to a few meters in the worst-case scenario for a GEO satellite and to less than 0.5 m in LEO.

In real-world applications, a finite amount of time will pass between the time when the GNSS solution is measured and the time at which such a solution becomes available to the onboard computer. Moreover, the velocity solution can be further affected by an additional latency, which introduces desynchronization between the times at which the output position and velocity are associated. Depending on the field of application, these delays can severely affect the accuracy of the navigation solution. For example, as Figure 4 highlights, an LEO satellite at an altitude of 550 km has an orbital speed of about 7.5 km/s such that even a delay of 0.1 s can result in errors of up to 750 m.

In this work, a dedicated post-processing algorithm is introduced to mitigate the impact of these delays on the navigation solution. In particular, the GNSS latency between the position and velocity solutions is neglected, because as Figure 4 shows, the velocity errors are negligible with respect to the position ones. Then, using a simple two-body model, the state provided by the GNSS is propagated with a forward Euler scheme to compensate for the average GNSS time delay, allowing reducing the position error from hundreds of meters down to a few tenths of meters. In practice, the propagation time can be selected by statistically analysing the inference time between the real hardware components.

When considering a realistic GNSS receiver for onorbit applications [1], position fix accuracies can be to the order of 1.5 m and 0.03 m per second in the ITRF. Additionally, it provides time accuracies of 20 nanoseconds. While cable and antenna delays can introduce errors into time measurements, their impact is typically negligible in spacecraft due to the short lengths of cables and antennas. Communication delays and velocity fix latencies can be to the order of 15–100 ms, also depending on the communication bus employed onboard.

With this in mind, Figure 5 illustrates the receiver model which we developed for this work. In this implementation, the true GCRF position and velocity are affected by bias and white or brown noise terms which simulate the internal processing and tracking of the GNSS signals. Using the EOP data, these are then converted to the ITRF using the exact IAU 2006/2000A formulation. All of the results are then delayed by a variable time latency and transformed into a digital signal. An additional latency is added to the velocity measurement to mimic the behaviour of real-world GNSS receivers. An internal clock model is also implemented to provide an absolute measurement of time to the onboard computer. The position and velocity outputs of the GNSS sensor are then post-processed and rotated to the GCRF according to the methodology explained above before being exploited by the error-state Kalman filter.

2.2. Inertial Navigation System Model

An INS typically consists of three orthogonal accelerometers and three orthogonal gyroscopes, providing a full six-axis measurement of linear accelerations and angular rates in the body-fixed reference frame. The accelerometers measure the specific force applied to the vehicle, while the gyroscopes detect the angular velocity. Together, these sensors enable the determination of the vehicle’s position, velocity and orientation relative to an inertial reference frame.

In the common strapdown configuration, the sensors are rigidly attached to the vehicle’s body, meaning they are directly subjected to the vehicle’s full motion, including angular rotations and vibrations. Unlike gimballed systems, which mechanically isolate the sensors, a strapdown INS requires real-time computational algorithms to transform body-fixed measurements into the inertial frame using the vehicle’s current attitude. This attitude is computed by integrating the angular rates measured by the gyroscopes over time.

Strapdown INS technology is widely used in modern space applications, where real-time processing and advanced algorithms enable accurate navigation without the mechanical limitations of traditional gimballed systems. The accelerometer output in the body frame ${\mathbf{a}}^b$ is evaluated using the direction cosine matrix, rotating a vector from the inertial to the body reference frame ${}_b\mathbf{A}_i$, as computed from the gyroscope, and correcting for the gravitational acceleration $\mathbf{g}(\mathbf{r})$ to model the specific non-gravitational force sensed by the sensor:

$${\mathbf{a}}^b={}_b\mathbf{A}_i\left(\ddot{\mathbf{r}}-\mathbf{g}(\mathbf{r})\right).$$

(3)

It is evident how in ideal free-fall orbiting scenarios, where the only force acting on the spacecraft is the gravitational pull of celestial bodies, an ideal accelerometer located at the spacecraft’s center of gravity would measure zero acceleration. This is because the accelerometer senses non-gravitational forces, and in such ideal conditions, no other external forces are present. To accurately compute the total acceleration acting on the spacecraft, which is necessary for propagating its motion, a gravitational field model must be employed. This model is used to correct the accelerometer’s readings by accounting for gravitational accelerations.

The navigation functions exploiting the inertial sensors must incorporate several elements: body rotation estimation, reference frame transformation, gravitational acceleration evaluation, correction for translational and rotational coupling and the integration of the system’s dynamics. In this work, the INS is developed following the scheme reported in Figure 6.

Accelerometers and gyroscope sensors face inherent errors which limit their measurement accuracy. They need to be carefully considered when used in the proposed navigation functions, For instance, there are deterministic errors, such as mounting errors or scale factors, which are repeatable and can be corrected through a fixed calibration; temperature-induced errors, which are correctable with thermal calibration; switch-on to switch-on bias errors, which are stochastic but constant within a session and correctable with onboard calibration; and in-run stochastic variations, which vary throughout the measurement run and can be corrected if their dynamics are slow enough. In fact, fast dynamics errors like random noise are difficult to compensate for.

All of these error terms affect the overall performance of the onboard orbit determination function, and a credible design should consider them in order to verify the achievable performance in a realistic scenario. The classical errors considered in the sensor models used in this work are described in the scheme in Figure 7. They compose a classic model for a sensor processing pipeline, illustrating the various sources of errors and distortions which can occur during the acquisition, processing and output of a physical signal.

Note that in addition to these sensor errors, extrinsic errors from estimation processes, like imprecise gravitational field models or numerical computation errors, can further affect the accuracy of position, velocity or attitude measurements in the orbit determination process.

3. Onboard Orbit Determination Algorithm

INSs are essential for spacecraft navigation, especially in scenarios where external references such as a GNSS are unreliable or unavailable. An INS estimates the position, velocity and orientation of a spacecraft using internal measurements of the acceleration and angular velocity. The core components of an INS include accelerometers, gyroscopes and sometimes additional sensors like magnetometers.

Inertial navigation relies on the principles of Newtonian mechanics and measures acceleration and angular rates to estimate the spacecraft’s trajectory and orientation over time. The state variables in an INS typically include the position $\mathbf{r}(t)$, velocity $\mathbf{v}(t)$ and attitude $\mathbf{q}(t)$. The position vector $\mathbf{r}(t)$ represents the spacecraft’s location, the velocity $\mathbf{v}(t)$ indicates the rate of change of position, and the attitude $\mathbf{q}(t)$ describes the spacecraft’s orientation, often represented using quaternions $\mathbf{q}(t)$ or Euler angles. The fundamental equations of motion in an INS are given by

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\mathbf{r}(t)}{dt}}& =\mathbf{v}(t)\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\mathbf{v}(t)}{dt}}& =\mathbf{a}(t)-\mathbf{g}(t)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\mathbf{q}(t)}{dt}}& ={\displaystyle \frac{1}{2}}\mathbf{\Omega }(t)\mathbf{q}(t)\hfill \end{array}$$

(6)

where $\mathbf{a}(t)$ is the specific acceleration measured by the accelerometers, $\mathbf{g}(t)$ is the gravitational acceleration which can be compensated for and $\mathbf{\Omega }(t)$ is a 4 × 4 matrix defined as follows:

$$\mathbf{\Omega }(t)=\left[\begin{array}{cc}0& -\mathit{\omega }{(t)}^{\mathrm{T}}\\ \mathit{\omega }(t)& -{\lfloor \mathit{\omega }(t)\rfloor }_{\mathrm{x}}\end{array}\right]$$

(7)

in which $\mathit{\omega }(t)$ is the angular velocity vector measured by the gyroscopes and ${\lfloor ·\rfloor }_{\mathrm{x}}$ denotes a skew-symmetric matrix.

The INS’s system-level error model is designed to account for various errors which can affect the system’s performance. These errors include sensor biases, scale factors, misalignment and random noise. Sensor biases are constant or slowly varying offsets in accelerometer and gyroscope measurements, which can cause drift in the position and velocity estimates. Scale factors introduce inaccuracies in sensor measurements, affecting the magnitude of the acceleration and angular rate. Misalignment errors arise from deviations in the sensor alignment with the coordinate system, impacting attitude estimates. Random noise represents high-frequency disturbances in the sensor measurements, typically modelled as zero-mean Gaussian noise.

To address these errors, the INS’s system-level error model includes an error state vector $\mathit{\delta }{\mathit{x}}_k$, which comprises

$$\mathit{\delta }{\mathit{x}}_k=\left[\begin{array}{c}\mathit{\delta }{\mathit{r}}_k\\ \mathit{\delta }{\mathit{v}}_k\\ \mathit{\delta }{\mathit{q}}_k\\ \mathit{\delta }{\mathit{b}}_{\omega }\\ \mathit{\delta }{\mathit{b}}_a\end{array}\right]$$

(8)

where $\mathit{\delta }{\mathit{r}}_k$ represents the position errors, $\mathit{\delta }{\mathit{v}}_k$ represents the velocity errors, $\mathit{\delta }{\mathit{q}}_k$ represents the attitude errors, $\mathit{\delta }{\mathit{b}}_{\omega }$ represents the gyroscope bias errors and $\mathit{\delta }{\mathit{b}}_a$ represents the accelerometer bias errors. The dynamics of these errors are governed by the discrete-time error model:

$$\mathit{\delta }{\mathit{x}}_k={\mathbf{\Phi }}_k\mathit{\delta }{\mathit{x}}_{k-1}+{\mathbf{w}}_{k-1}$$

(9)

where ${\mathbf{\Phi }}_k$ is the state transition matrix which describes the evolution of errors over time and ${\mathbf{w}}_{k-1}$ is the process noise vector, typically modelled as zero-mean Gaussian noise with covariance ${\mathbf{Q}}_k$. The process noise covariance ${\mathbf{Q}}_k$ accounts for uncertainties in sensor errors and other model inaccuracies:

$${\mathbf{Q}}_k=\left[\begin{array}{ccccc}{\mathbf{Q}}_p& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}_v& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{Q}}_q& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{Q}}_{b_{\omega }}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{Q}}_{b_a}\end{array}\right]$$

(10)

where ${\mathbf{Q}}_p$, ${\mathbf{Q}}_v$, ${\mathbf{Q}}_q$, ${\mathbf{Q}}_{b_{\omega }}$ and ${\mathbf{Q}}_{b_a}$ are the covariances associated with the position, velocity, attitude, accelerometer bias and gyroscope bias errors, respectively.

The measurement model relates the observed measurements to the true state and includes measurement noise:

$${\mathbf{z}}_k={\mathbf{H}}_k\mathit{\delta }{\mathit{x}}_k+{\mathbf{v}}_k$$

(11)

where ${\mathbf{z}}_k$ is the measurement vector, ${\mathbf{H}}_k$ is the measurement matrix which maps the error state to measurements and ${\mathbf{v}}_k$ is the measurement noise, assumed to be zero-mean Gaussian noise, with covariance ${\mathbf{R}}_k$. Incorporating these models into the Kalman filter framework allows the INS to estimate the spacecraft’s position, velocity and attitude while compensating for various errors and uncertainties. This comprehensive approach ensures improved accuracy and reliability in spacecraft navigation. In case the system implements a stand-alone attitude determination algorithm, the above derivation can be reduced to include only uncertainties coming from the position and velocity estimation. This means that the error state vector is reduced to

$$\mathit{\delta }{\mathit{x}}_k=\left[\begin{array}{c}\mathit{\delta }{\mathit{r}}_k\\ \mathit{\delta }{\mathit{v}}_k\\ \delta {\mathit{b}}_a\end{array}\right]$$

(12)

and consequently the process noise covariance ${\mathbf{Q}}_k$ accounts for uncertainties in sensor errors and other model inaccuracies:

$${\mathbf{Q}}_k=\left[\begin{array}{ccc}{\mathbf{Q}}_p& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}_v& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{Q}}_{b_a}\end{array}\right]$$

(13)

3.1. Derivation of the Error-State Kalman Filter

The error-state Kalman filter is designed to estimate the errors in state variables such as the position, velocity and attitude rather than the state variables themselves. This is particularly beneficial for spacecraft navigation near the International Space Station (ISS), where precise error estimation can significantly enhance mission performance.

Consider the continuous-time state-space model

$$\dot{\mathbf{x}}(t)=\mathbf{F}(t)\mathbf{x}(t)+\mathbf{G}(t)\mathbf{u}(t)+\mathbf{w}(t)$$

(14)

where $\mathbf{x}(t)$ is the state vector, which may include the position, velocity and attitude, $\mathbf{F}(t)$ is the state matrix, $\mathbf{G}(t)$ is the control input matrix, $\mathbf{u}(t)$ is the control input vector and $\mathbf{w}(t)$ is the process noise, which is assumed to be zero-mean Gaussian noise, with covariance $\mathbf{Q}(t)$.

The measurement model is given by

$$\mathbf{z}(t)=\mathbf{H}(t)\mathbf{x}(t)+\mathbf{v}(t)$$

(15)

where $\mathbf{z}(t)$ is the measurement vector, which may include GNSS measurements and INS outputs, $\mathbf{H}(t)$ is the measurement matrix and $\mathbf{v}(t)$ is the measurement noise, assumed to be zero-mean Gaussian noise, with covariance $\mathbf{R}(t)$.

Let $\widehat{\mathbf{x}}(t)$ be the estimate of the state $\mathbf{x}(t)$. The error state $\delta \mathbf{x}(t)$ is defined as follows:

$$\delta \mathbf{x}(t)=\mathbf{x}(t)-\widehat{\mathbf{x}}(t)$$

(16)

The error dynamics can be derived by differentiating $\delta \mathbf{x}(t)$:

$${\displaystyle \frac{d}{dt}}\delta \mathbf{x}(t)=\dot{\mathbf{x}}(t)-\dot{\widehat{\mathbf{x}}}(t)$$

(17)

Using the state-space model, we obtain

$${\displaystyle \frac{d}{dt}}\delta \mathbf{x}(t)=\mathbf{F}(t)\mathbf{x}(t)+\mathbf{G}(t)\mathbf{u}(t)+\mathbf{w}(t)-\left[\mathbf{F}(t)\widehat{\mathbf{x}}(t)+\mathbf{G}(t)\mathbf{u}(t)\right]$$

(18)

By simplifying this, we obtain

$${\displaystyle \frac{d}{dt}}\delta \mathbf{x}(t)=\mathbf{F}(t)\delta \mathbf{x}(t)+\mathbf{w}(t)$$

(19)

For practical implementation, we often use the discrete-time form. Assuming a sampling interval $\Delta t$, the discrete-time error dynamics are

$$\delta {\mathbf{x}}_k={\mathbf{\Phi }}_k\delta {\mathbf{x}}_{k-1}+{\mathbf{w}}_k$$

(20)

where $\delta {\mathbf{x}}_k$ is the error state at time step k, ${\mathbf{\Phi }}_k$ is the state transition matrix over $\Delta t$ and ${\mathbf{w}}_k$ is the discrete process noise. The discrete-time measurement model is

$${\mathbf{z}}_k={\mathbf{H}}_k\delta {\mathbf{x}}_k+{\mathbf{v}}_k$$

(21)

The Kalman filter algorithm consists of two main steps: predict and update. In the prediction step, we propagate the error state and error covariance forward in time. The predicted error covariance is

$${\mathbf{P}}_{k|k-1}={\mathbf{\Phi }}_{k-1}{\mathbf{P}}_{k-1}{\mathbf{\Phi }}_{k-1}^T+{\mathbf{Q}}_{k-1}$$

(22)

where ${\mathbf{P}}_{k-1}$ is the error covariance at the previous time step. To calculate ${\mathbf{\Phi }}_k$, we need the Jacobian matrix of the state transition function. Assuming that the state transition function $\mathbf{\Phi }(t,t_{k-1})$ can be approximated linearly, we have

$${\mathbf{\Phi }}_k\approx \mathbf{I}+{\mathbf{F}}_k\Delta t$$

(23)

where ${\mathbf{F}}_k$ is the state jacobian matrix evaluated at time $t_k$.

In the update step, we incorporate the new measurements to correct the error state and covariance. The Kalman gain ${\mathbf{K}}_k$ is computed as follows:

$${\mathbf{K}}_k={\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T{\left({\mathbf{H}}_k{\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T+{\mathbf{R}}_k\right)}^{-1}$$

(24)

The updated error state estimate is

$$\delta {\widehat{\mathbf{x}}}_k=\delta {\widehat{\mathbf{x}}}_{k|k-1}+{\mathbf{K}}_k\left({\mathbf{z}}_k-{\mathbf{H}}_k\delta {\widehat{\mathbf{x}}}_{k|k-1}\right)$$

(25)

where $\delta {\widehat{\mathbf{x}}}_{k|k-1}$ is the predicted error state. The updated error covariance uses the Joseph formula:

$${\mathit{P}}_k^{+}=(\mathit{I}-{\mathit{K}}_k{\mathit{H}}_k){\mathit{P}}_k^{-}{(\mathit{I}-{\mathit{K}}_k{\mathit{H}}_k)}^T+{\mathit{K}}_k\mathit{R}{\mathit{K}}_k^T$$

(26)

The process noise covariance ${\mathbf{Q}}_k$ and the measurement noise covariance ${\mathbf{R}}_k$ are generally provided by the system’s specifications or can be empirically determined. On the other hand, similar to the state transition matrix computation, the measurement matrix ${\mathbf{H}}_k$ may also need to be linearized if the measurement function is nonlinear:

$${\mathbf{H}}_k={\left.{\displaystyle \frac{\partial \mathbf{h}}{\partial \mathbf{x}}}\right|}_{\mathbf{x}={\widehat{\mathbf{x}}}_{k|k-1}}$$

(27)

where $\mathbf{h}$ represents the nonlinear measurement function.

3.2. Derivation of the Jacobian Matrices

In the following section, we derive the linearized transition and measurement matrices for the case of complete orbit and attitude reconstruction. The equations can be greatly simplified if the attitude determination is performed separately. Recalling Equation (4), we can derive the error differential equations as follows:

$$\left[\begin{array}{c}\dot{\delta \mathbf{r}}\\ \dot{\delta \mathbf{v}}\\ \dot{\delta \mathbf{q}}\end{array}\right]=\left[\begin{array}{c}\delta \mathbf{v}\\ {\mathbf{R}}_B^I(\widehat{\mathbf{q}})\left[2{\mathbf{A}}^T{\mathbf{C}}^T\delta \mathbf{q}+\delta \mathbf{f}\right]+{\displaystyle \frac{\partial \mathbf{g}(\mathbf{r})}{\partial \mathbf{r}}}\delta \mathbf{r}\\ \mathbf{Q}(\widehat{\mathbf{q}})\delta \mathit{\omega }+\mathbf{\Omega }(\widehat{\mathit{\omega }})\delta \mathbf{q}\end{array}\right]$$

(28)

in which ${\mathbf{R}}_B^I$ is the rotation of the INS body with respect to inertial space while $\delta \mathit{\omega }$ and $\delta \mathbf{f}$ express the errors in the accelerometer force and gyro angular velocity measurements, respectively. Following [19], we introduce the following matrices:

$$\begin{array}{cc}\hfill \mathbf{\Omega }(\mathit{\omega })& ={\displaystyle \frac{1}{2}}\left[\begin{array}{cccc}0& {\omega }_3& -{\omega }_2& {\omega }_1\\ -{\omega }_3& 0& {\omega }_1& {\omega }_2\\ {\omega }_2& -{\omega }_1& 0& {\omega }_3\\ -{\omega }_1& -{\omega }_2& -{\omega }_3& 0\end{array}\right]\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \mathbf{Q}(\mathbf{q})& ={\displaystyle \frac{1}{2}}\left[\begin{array}{ccc}{q}_4& -q_3& q_2\\ q_3& q_4& -q_1\\ -q_2& q_1& q_4\\ -q_1& -q_2& -q_3\end{array}\right]\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \mathbf{C}& =\left[\begin{array}{cc}2\mathbf{Q}& -\mathbf{q}\end{array}\right]\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \mathbf{A}& =\left[\begin{array}{ccc}0& -f_3& f_2\\ f_3& 0& -f_1\\ -f_2& f_1& 0\\ -f_1& -f_2& -f_3\end{array}\right]\hfill \end{array}$$

(32)

where ${\omega }_1$, ${\omega }_2$ and ${\omega }_3$ are the body frame components of the angular velocity. The gravity gradient matrix $\mathbf{G}={\displaystyle \frac{\partial \mathbf{g}(\mathbf{r})}{\partial \mathbf{r}}}$ for a two-body gravity field is defined as follows:

$${\mathbf{G}}_{3\times 3}=\left[\begin{array}{ccc}{\displaystyle \frac{3\mu x^2}{r^5}}-{\displaystyle \frac{\mu }{r^3}}& {\displaystyle \frac{3\mu xy}{r^5}}& {\displaystyle \frac{3\mu xz}{r^5}}\\ {\displaystyle \frac{3\mu xy}{r^5}}& {\displaystyle \frac{3\mu y^2}{r^5}}-{\displaystyle \frac{\mu }{r^3}}& {\displaystyle \frac{3\mu yz}{r^5}}\\ {\displaystyle \frac{3\mu xz}{r^5}}& {\displaystyle \frac{3\mu yz}{r^5}}& {\displaystyle \frac{3\mu z^2}{r^5}}-{\displaystyle \frac{\mu }{r^3}}\end{array}\right]$$

(33)

Thus, the INS error quantities in Equation (28) can be rearranged as follows:

$$\left[\begin{array}{c}\dot{\delta \mathbf{r}}\\ \dot{\delta \mathbf{v}}\\ \dot{\delta \mathbf{q}}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{0}}_{3\times 3}& {\mathbf{I}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{G}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{D}}_{3\times 3}\\ {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& \mathbf{\Omega }(\widehat{\mathit{\omega }})\end{array}\right]\left[\begin{array}{c}\delta \mathbf{r}\\ \delta \mathbf{v}\\ \delta \mathbf{q}\end{array}\right]+\left[\begin{array}{cc}{\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{R}}_B^I(\widehat{\mathbf{q}})\\ \mathbf{Q}(\widehat{\mathbf{q}})& {\mathbf{0}}_{4\times 3}\end{array}\right]\left[\begin{array}{c}\delta \mathit{\omega }\\ \delta \mathbf{f}\end{array}\right]$$

(34)

where $\mathbf{D}=2{\mathbf{R}}_B^I{\mathbf{A}}^T{\mathbf{C}}^T$. Finally, given the error-state vector of Equation (8), the full Jacobian matrix of the coupled orbital-attitude error dynamics is:

$$\mathbf{F}=\left[\begin{array}{ccccc}{\mathbf{0}}_{3\times 3}& {\mathbf{I}}_{3\times 3}& {\mathbf{0}}_{3\times 4}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{G}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{D}}_{3\times 4}& {\mathbf{0}}_{3\times 3}& {\mathbf{R}}_B^I(\widehat{\mathbf{q}})\\ {\mathbf{0}}_{4\times 3}& {\mathbf{0}}_{4\times 3}& \mathbf{\Omega }(\widehat{\mathit{\omega }})& \mathbf{Q}(\widehat{\mathbf{q}})& {\mathbf{0}}_{4\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{F}}_{b_{\omega }}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{F}}_{b_a}\end{array}\right]$$

(35)

in which it is assumed that the gyroscope and accelerometer error models used within the filter, which provide the $\delta \mathit{\omega }$ and $\delta \mathbf{f}$ terms, respectively, only include the bias contributions. Furthermore, if these biases are modelled as first-order Gauss–Markov processes, then

$${\mathbf{F}}_{b_a}=-{\displaystyle \frac{1}{{\tau }_a}}{\mathbf{I}}_{3\times 3}$$

(36)

$${\mathbf{F}}_{b_{\omega }}=-{\displaystyle \frac{1}{{\tau }_{\omega }}}{\mathbf{I}}_{3\times 3}$$

(37)

where ${\tau }_a$ and ${\tau }_{\omega }$ are the correlation times of the accelerometer and gyroscope biases.

In a loosely coupled INS/GNSS configuration using an error-state Kalman filter, the linearized measurement matrix $\mathbf{H}$ relates the error-state vector to the GNSS measurement residuals $\delta \mathbf{r}$ and $\delta \mathbf{v}$ (i.e., the differences between the state provided by the GNSS and that predicted by the INS):

$$\delta \mathbf{r}={\mathbf{r}}_{\mathrm{GNSS}}-{\mathbf{r}}_{\mathrm{INS}},\delta \mathbf{v}={\mathbf{v}}_{\mathrm{GNSS}}-{\mathbf{v}}_{\mathrm{INS}}.$$

(38)

Thus, the linearized measurement matrix $\mathbf{H}$ takes the form of

$$\mathbf{H}=\left[\begin{array}{ccccc}{\mathbf{I}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{I}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\end{array}\right].$$

(39)

in which the two row blocks map the position and velocity errors to the GNSS position and velocity residuals, respectively. Please notice that in this configuration, the orientation errors ($\delta \mathit{q}$) and bias terms ($\delta {\mathbf{b}}_a,\delta {\mathbf{b}}_{\omega }$) do not directly affect the GNSS position and velocity measurement residuals, and thus their corresponding blocks are null.

3.3. Onboard Orbit Determination During GNSS Outages

The loosely coupled GNSS/INS integration relies on the continuous availability of a GNSS fix. However, in certain situations, such as in orbit, the GNSS receiver may temporarily lose the signal from multiple GNSS satellites. This is particularly likely in high-altitude orbits or during mission phases where the spacecraft’s attitude angular rates exceed a specific threshold. To ensure uninterrupted orbit determination onboard, dead reckoning navigation can be supported by integrating an onboard orbital propagator with the navigation filters.

In the event of a GNSS signal loss, when position and velocity measurements are unavailable, the filter architecture is seamlessly adapted. Instead of relying on GNSS data, the system switches to the output of the orbital propagator, which is initialized based on the last GNSS-supported navigation solution. The propagated orbital states serve as artificial measurements, fed into the navigation filters to maintain position estimates. During this phase, the filter parameters are adjusted to account for changes in state covariances and measurement models. In particular, the covariance matrix is reset, enlarging the measurement uncertainties based on the known mechanization uncertainties (i.e., the INS contribution together with onboard numerical propagation.) Furthermore, the bias estimates are no longer updated while the system operates under these conditions.

The navigation filter performs dynamical mechanization using an orbital propagator, where the selection of the appropriate method must strike a balance between accuracy and computational efficiency, especially for real-time onboard applications. Several well-established numerical orbital propagators, such as simplified general perturbations (SGP), SGP4 and the perturbed Kepler method, can be employed in the proposed GNSS/INS onboard orbit navigation system. These propagators have a rich history, having been developed to meet the growing need for accurate satellite tracking and prediction in both scientific and defense applications.

The SGP and SGP4 models were initially developed by the United States Air Force for satellite tracking purposes in the 1960s. SGP4, a more refined version of the SGP model, became the standard for propagating TLE sets, which provide orbital parameters for thousands of satellites. SGP4 is particularly efficient in modelling LEO satellites, taking into account perturbations such as Earth’s oblateness (J2 effect), atmospheric drag and resonance effects. Its widespread adoption is due to its balance between accuracy and speed, making it suitable for operational satellite tracking, where real-time performance is critical.

On the other hand, the perturbed Kepler method, while based on classical two-body dynamics, is enhanced by incorporating perturbative forces which significantly affect satellite orbits. These forces include the Earth’s J2 oblateness effect, which causes precession of the orbit, and atmospheric drag, which gradually reduces the altitude of satellites in LEO. The perturbed Kepler approach is more flexible than SGP4 in that it can be adapted to include additional perturbations, such as solar radiation pressure or third-body gravitational effects, depending on the mission requirements. It is also versatile in its ability to accept both classical orbital elements and state vectors in the form of position and velocity, making it ideal for integration with a GNSS/INS system.

For the proposed GNSS/INS navigation system, any of these orbital propagators can be applied, provided they can accurately propagate the satellite’s state vector. Whether initialized with orbital elements or directly with position and velocity data, the output must be in position and velocity form to ensure compatibility and interchangeability with GNSS measurements. The propagation formulation based on the numerical integration of position and velocity vectors has demonstrated greater robustness in the face of numerical degradation, with respect to orbital element methods. In fact, it is less prone to the accumulation of errors which can affect individual orbital elements, such as eccentricity or inclination.

The results presented in this paper specifically use a perturbed Kepler orbital propagator, where the position and velocity are propagated while accounting for the $J_2$ effect and atmospheric drag using an exponential atmosphere model. Specifically, the spacecraft acceleration is computed as follows:

$$\mathbf{a}(t)={\mathbf{a}}_{2B}(\mathbf{r})+{\mathbf{a}}_{J_2}(\mathbf{r})+{\mathbf{a}}_d(\mathbf{r},\mathbf{v})+{\tilde{\mathbf{a}}}_{\mathrm{IMU}}(t)$$

(40)

where ${\mathbf{a}}_{2B}$ and ${\mathbf{a}}_{J_2}$ are Earth’s central gravity and $J_2$ contributions, respectively, ${\tilde{\mathbf{a}}}_{\mathrm{IMU}}$ represents the non-gravitational accelerations estimated by the IMU and ${\mathbf{a}}_d$ is the perturbing acceleration due to the atmospheric drag, computed with a standard cannonball model according to Equation (41):

$${\mathbf{a}}_d=-{\displaystyle \frac{1}{2}}\rho (\mathbf{r})v_r{\displaystyle \frac{C_{d}A}{m}}{\mathbf{v}}_r$$

(41)

in which $\rho$ is the atmospheric air density, ${\mathbf{v}}_r$ is the spacecraft velocity relative to the atmosphere, $C_d$ is the drag coefficient, A is the cross-sectional area (i.e., the area normal to the relative velocity) and m is the spacecraft mass. This configuration was chosen to simulate a worst-case scenario, with limited the fidelity models to stress test the system’s robustness. The inclusion of only two perturbations ensures that the system remains computationally efficient while still addressing the dominant forces acting on most satellites in Earth’s orbit. Please notice that when simulations are performed in GEO, the drag contribution is also removed because of its negligible impact. By demonstrating the system’s ability to maintain accurate orbit determination under these challenging conditions, we show that the GNSS/INS navigation system is broadly applicable to a wide range of satellite missions, from small nano-satellites to larger Earth observation satellites operating in various orbital regimes.

Note that the performance of the proposed orbit determination method in case of GNSS outages, namely coupling the orbital propagator with the error-state Kalman filter fed INS measurements, is not remarkably higher than a simpler isolated orbital propagator initialized with the last valid GNSS-supported orbit determination solution. However, the proposed formulation guarantees a more robust solution and is capable of supporting maneuvered phases and thrusting arcs. Moreover, by including the control and actuation effects in the estimation filter, it is possible to identify differences between the computed control commands and the actuated ones. This is feasible for mission phases with or without the GNSS measurements available.

The scheme of the GNSS/INS orbit determination filter with the integrated orbit propagation section is shown in Figure 8.

Ground Updates of the Onboard Orbit Determination

In typical mission scenarios, GNSS outages may last for several tens of minutes. However, in the event of extended outages or a complete failure of the GNSS receiver, the proposed GNSS/INS system can continue to operate by leveraging orbit determination solutions from ground-based data, such as in TLE format. While this approach allows the system to remain operational, it is important to note that the navigation accuracy in this case will be limited to the precision of the ground-based TLE data, which are generally less accurate than GNSS-derived solutions, particularly in LEO. Nevertheless, the onboard orbital propagation and the INS-coupled filter can continue functioning using these data.

TLEs are widely used for satellite tracking and provide a compact representation of a satellite’s orbital parameters. These elements are updated periodically by organizations, such as NORAD, and are often accurate to within 1 km in LEO [20]. The TLE consists of various orbital elements which describe the satellite’s orbit at a given epoch. However, since they are based upon the SGP4 propagation model, their accuracy rapidly degrades over time due to the action of unmodelled perturbations. In particular, the omission of the $J_{22}$ effect, approximations in the solar and lunar orbit computation, and uncertainties over the atmospheric drag introduce systematic periodic errors as large as 1 km for LEO satellites [21,22], making them unsuitable for missions requiring high-precision orbit determination.

In the proposed onboard orbit determination architecture, the uploaded TLE must be synchronized with the actual universal time. This synchronization is achieved by updating the satellite’s orbital elements, exploiting the mean motion and its rates of change (i.e., mean motion variation) according to the time difference between the epoch of TLE generation and the onboard update. Ensuring the time differential remains within a specific threshold is crucial, as outdated TLE data can significantly impact navigation accuracy. As part of this process, the system must also convert the TLE data from the true equator mean equinox (TEME) reference frame, in which TLEs are defined, into the inertial GCRF frame used onboard for orbital propagation. This frame transformation ensures compatibility between the TLE data and the onboard propagator, allowing for accurate propagation of the satellite’s position and velocity.

While the use of TLEs enhances system robustness in the face of GNSS failures or extended outages, it should be recognized that the accuracy provided by ground-based TLE solutions is not sufficient for certain high-precision operations. For example, proximity operations or missions requiring stringent orbital navigation—such as docking, formation flying or autonomous rendezvousing—would suffer from the lower fidelity of TLE-derived solutions. In such cases, more accurate onboard redundant solutions would be necessary.

On the other hand, for short-duration GNSS outages, the coupled INS and onboard orbital propagator filter can reliably maintain orbit determination, ensuring continuity of navigation even for missions where relative accuracy is critical. Thus, while the use of TLE data offers a fallback solution during extended outages or definitive failures, the real strength of the proposed GNSS/INS system lies in its ability to autonomously navigate during short GNSS outages, maintaining a reliable and accurate orbit determination solution.

3.4. Delay Handling in Real-Time Implementation

Real-world estimation processes are affected by several delays which are typically neglected in early GNC modelling phases. However, depending on the type of application, these delays could severely affect the navigation performance if they are not properly accounted for. In this regard, Figure 9 reports an overview of the operational timeline of a traditional predict-update Kalman filter. In particular, sensor measurements are shifted forward with respect to $t_k$, depending on the order in which they are retrieved by the filter and the latencies over the communication buses. When all of these delays are taken into account, it becomes clear that the filter will estimate a state at $t_x$, lying between $t_k$ and $t_y$. Then, if these delays remain approximately constant, after initialisation, the prediction step will start propagating the state at $t_x$ to one which should align closer to the one associated to the measurements. Therefore, the actual time delay between the state provided by the filter and the real one at $t_s$ will be $\Delta KG+\Delta T_x$, where $\Delta KG$ is the time required for computation of the Kalman gain together with the state and covariance update equations, while $\Delta T_x=t_y-t_x$ is the time offset between the actual state estimated by the filter and the time at which the filter has received all of the measurements. Within this work, we only model the GPS sensor latencies, as reported in

Table 1. Summary of the main simulation settings for the LEO scenario.

Parameter	Value
Simulation start date	1 April 2020, 12:30:00 p.m. UTC
Simulation duration	95 min
Initial state distribution ($1\sigma$)	Semi-major axis: 6921 km ± 1 km
	Eccentricity: 0.001 ± 0.0
	Inclination: 98.88 deg ± 0.4 mdeg
	RAAN: 324.12 deg ± 0.4 mdeg
	Argument of perigee: 337.85 deg ± 0.4 mdeg
	True anomaly: 17.80 deg ± 0.4 mdeg
Satellite mass	17.5 kg
Satellite dimensions	Spacecraft body: 0.1 m × 0.1 m × 0.3 m
	Solar arrays: 0.3 m × 0.2 m
Drag coefficient	2.2 ± 0.18
Optical coefficients (1$\sigma$)	Diffusion ${\rho }_d$: 0.1 ± 0.017
	Specular ${\rho }_s$: 0.1 ± 0.017
GNSS bias ($1\sigma$)	10 cm
GNSS noise ($1\sigma$)	1.5 m
GNSS delay ($1\sigma$)	15 ms ± 7.5 ms
GNSS delay noise ($1\sigma$)	1 ms ± 0.5 ms
IMU static bias ($1\sigma$)	270 μg
IMU velocity random walk ($1\sigma$)	0.008 m/s/$\sqrt{(}h)$
IMU bias stability ($1\sigma$)	3.7 μg
IMU scale factors ($1\sigma$)	Linear: 34.34 ppm
	Asimmetry: 1 ppm
	Nonlinearity: 3000 ppm
Initial filter uncertainty ($1\sigma$)	Position ${\sigma }_r$ (x, y, z): 8000 km
	Velocity ${\sigma }_v$ (x, y, z): 1000 m/s
Process noise covariance ($1\sigma$)	Position $q_p$ (x, y, z): 2 m
	Velocity $q_v$ (x, y, z): 0.1 m/s
	Bias $q_b$ (x, y, z): 0.001 m/s2
Measurement noise covariance ($1\sigma$)	GNSS position $R_p$ (x, y, z): 5 m
	GNSS velocity $R_v$ (x, y, z): 0.5 m/s
Filter rate	10 Hz

and

Table 2. Summary of the main simulation settings for the GEO scenario.

Parameter	Value
Simulation duration	12 h
Initial state distribution ($1\sigma$)	Semi-major axis: 42,164 km ± 0.1 km
	Eccentricity: 0 ± 0.0001
	Inclination: 0 deg ± 0.4 mdeg
	RAAN: 0 deg ± 0.4 mdeg
	Argument of perigee: 0 deg ± 0.4 mdeg
	True anomaly: 68.97 deg ± 0.4 mdeg
GNSS bias ($1\sigma$)	1 m
GNSS noise ($1\sigma$)	250 m
GNSS delay ($1\sigma$)	50 ms ± 25 ms
GNSS delay noise ($1\sigma$)	20 ms ± 10 ms
Initial filter uncertainty ($1\sigma$)	Position ${\sigma }_r$ (x, y, z): 300 m
	Velocity ${\sigma }_v$ (x, y, z): 10 m/s
Process noise covariance ($1\sigma$)	Position $q_p$ (x, y, z): 25 m
	Velocity $q_v$ (x, y, z): 0.2 m/s
	Bias $q_{b_a}$ (x, y, z): 0.002 m/s2
Measurement noise covariance ($1\sigma$)	GNSS position $R_p$ (x, y, z): 20 m
	GNSS velocity $R_v$ (x, y, z): 1 m/s
Filter rate	10 Hz

, and assume that the others can be neglected at the current stage.

4. Orbit Determination Performance

A set of Montecarlo campaigns was performed to validate the proposed navigation architecture and test its robustness across different operational scenarios in both nominal and off-nominal conditions (i.e., with failures of the GNSS receiver). A validated functional engineering simulator (FES) was used for the simulations [1]. It featured a high-accuracy orbital propagator which included Earth’s gravity harmonics (EGM2008 up to 60 × 60) as well as the Sun and Moon’s third-body effects. The solar radiation pressure (SRP) and drag contributions were evaluated using the standard plate model, in which the spacecraft body was modelled as a 6U CubeSat with deployable solar panels. The IMU and GNSS sensors were implemented according to the formulation reported in Section 2, while their parameters were tuned to match the Allan variance of real-world commercial sensors. Time delays with both deterministic and stochastic components were added between the sensors and the navigation block to simulate the real sampling behaviour of standard hardware for CubeSat missions.

4.1. LEO Nominal Scenario

The performance of the proposed orbit determination algorithm was initially tested through a Monte Carlo simulation, with 100 runs in a nominal LEO scenario featuring a spacecraft in a 550 km sun-synchronous orbit. The initial conditions for the true spacecraft trajectory were generated from a distribution around the nominal Keplerian parameters which was in line with traditional TLE errors: 1 km for the semi-major axis and 0.4 mdeg for the angular measures. Similarly, the filter’s initial state was generated from the same distribution, but with the true anomaly set to a random value between 0 and 360 deg. This was intended to mimick a real operational scenario, in which the actual true anomaly of the spacecraft at injection may differ from the one embedded within the flight software. For this reason, the initial filter covariance was increased to account for all possible orbital positions. A recap of the main simulation settings is reported in Table 1. Additional parameters are available upon request from the authors.

Figure 10 highlights that the average errors for the position components were limited to a few tenths of meters in the worst case. In particular, while the average error for the Z component was within 10 m (1$\sigma$), the X and Y components exhibited larger oscillations (up to 40 m $1\sigma$) because the CPNc model used for the ITRF-to-GCRF conversion led to greater approximation errors in the horizontal components. Within these simulations, the high accuracy of the GNSS solution implies that when the receiver operates in nominal conditions, the filter continuously corrects the drift of the INS such that the solution provided by the navigation block strictly follows the GNSS output. Additionally, despite the large initialization errors (due to the random position along the orbit), the filter required less than a minute to make the INS solution converge towards the correct state. Concerning the velocity, Figure 11 shows that the worst estimation error reached values up to 1 m/s for the X and Z components, while it remained near 0.5 m/s for the Y component.

The benefits of the post-processing propagation step applied to the original GNSS output to compensate for the delay between the time the GNSS solution were measured, and the time at which they became available to the onboard computer are shown in Figure 12. Specifically, when the time correction was not applied, the average error on each position component increased from a few tenths of meters up to 150 m, confirming that even a correction of a few milliseconds had a key impact on the navigation accuracy.

4.2. GEO Orbit Nominal Scenario

To further stress the capabilities of the architecture, we tested it on a spacecraft placed in GEO. The major differences of this scenario are related to the accuracy of the GNSS solution, which is inherently affected by the lower number of signals available. Indeed, since GNSS constellations are located at lower altitudes, most of the available signals are those from the side lobes of the GNSS antennas of the satellites which are not obstructed by Earth. Thus, the combination of a reduced number of satellites, lower signal-to-noise ratio (SNR) and bad constellation geometry (as seen from a user in GEO) caused a drop in the GNSS accuracy. In practice, this behaviour was simulated by updating the GNSS receiver settings as reported in Table 2. For example, an increased delay simulated the longer processing time due to the weaker GNSS signals, while greater delay noise accounted for the greater variability in the signal strength. On the other hand, the IMU settings and all other non-repeated parameters remained unchanged with respect to Table 1.

The resulting position estimation error is reported in Figure 13. In this scenario, the lower accuracy of the GNSS solution was reflected in average estimation errors of up to 400 m (3$\sigma$) for the horizontal components, while it remained bounded to less than 5 m (3$\sigma$) in the vertical direction.

4.3. GNSS Outage Scenarios

As highlighted in Section 4.1, when the GNSS state was available, the filter strictly followed the receiver output, with a maximum error to the order of a few tenths of meters. However, the novelty of this loosely coupled architecture is its ability to cope with GNSS outages or failures by replacing the GNSS output with the state provided by an onboard dynamic propagator. To simulate such a scenario, we performed another Monte Carlo simulation of 50 runs, in which a GNSS receiver failure was introduced starting from a random time sampled from a uniform distribution between 8 and 16 min and lasted a duration which randomly varied between 5 and 20 min.

The results in terms of position estimation error in LEO are reported in Figure 14, which exhibited a standard diverging pattern when the GNSS measures were removed. In particular, when the GNSS unavailability time was below 10 min, the navigation error remained below 1500 m (3$\sigma$), while for an average outage duration of 20 min, the maximum navigation error increased to about 3000 m. On the other hand, in the GEO scenario reported in Figure 15, the larger baseline estimation error caused worst-case peaks of about 11 km at the end of the outage period.

It is important to highlight that although the resulting accuracies did not differ much from those obtained by an isolated propagator which used as the initial state the latest available GNSS solution, the advantage of this approach is its capability to account for the real control action and unmodelled accelerations via the IMU measurements. This statement is proven in Figure 16, which displays two simulations run in an LEO scenario, with orbital control maneuvers performed during the GNSS outage time. In particular, in each sample run, a finite-time maneuver was performed through the onboard thrusters with a random duration between 60 and 480 s and with a start time which varied randomly within the GNSS outage time. The acceleration magnitude was fixed to 0.03 m/s², while its direction was also selected at random. The results display that by taking into account the real control action through the IMU feedback, the propagation error during the outage time was reduced to 5000 m (3$\sigma$), whereas a simple propagation of the latest available GNSS state resulted in a maximum range error of 9000 m (3$\sigma$).

Moreover, this position error was strictly related to the accuracy of the dedicated onboard propagator and, depending on the capabilities of the onboard computer, could be further improved either by adopting more accurate integration schemes or dynamic perturbation models. Finally, once the GNSS returned to nominal conditions, the filter only required a few steps to converge towards the correct solution, mitigating the entire drift accumulated during the propagation.

5. Final Remarks

This research contributes to advancing spacecraft autonomy by providing a practical and efficient solution for onboard orbit determination. The integration of error Kalman filtering within a GNSS/INS and orbit propagation architecture offers a promising pathway towards enhancing spacecraft operational autonomy, reducing dependence on ground support and enabling missions which require a high level of self-sufficiency in navigation and orbital determination. The nominal system provided an accuracy for each axis better than 40 m and 1.5 m/s at $3\sigma$ in LEO. The accuracy was worsened up to a few hundred meters in GEO. In the case of GNSS data outages and failures up to 30 min, the system was able to satisfy relative navigation requirements with accuracies below 3 km at 3$\sigma$ in LEO and to the order of 10 km at $3\sigma$ in GEO orbit.