Abstract
This paper develops an asynchronous measurement processing technique for sequential filtering that effectively handles small errors in the measurement sampling epoch within a linearized framework. The derived method relaxes the assumption that sensing systems generate and communicate measurements instantaneously and suggests a linearized method for extracting information from latent measurements via a temporal measurement update that considers uncertainty in the measurement acquisition epoch. To investigate performance, numerical simulations are performed utilizing the consider/neglect extended Kalman filter framework applied to a lunar descent-to-landing scenario in which latent vision-based measurements with uncertain acquisition times are used to navigate the vehicle. Through Monte Carlo simulation and analysis, this paper shows that the presented approach can be used to maintain filter consistency for latent measurements with low measurement-time uncertainties. Furthermore, an error budget and sensitivity analysis are presented to provide insight into the impact of the measurement-time uncertainty on navigation performance.
1 INTRODUCTION
In the design of navigation and estimation systems, it is commonplace to assume that measurements are collected coincident with their processing time. Unfortunately, this assumption unravels when applied to a real-world system. In such a setting, the measurements must be sampled, processed, and communicated to the navigation computer from independent and asynchronously operating sensing hardware. For example, vision-based measurements utilized for terrain relative navigation (TRN) often require demanding image-processing techniques that can lead to significant delays in data delivery, dependent upon the selected image-processing algorithms (Woicke et al., 2018). Uncertainty in the sampling time of the measurement may result from many factors, including computational load, measurement preprocessing, and communications buffering. For the remainder of this paper, uncertainty in the measurement sampling epoch is referred to as jitter. While latency and jitter are the phenomena of interest in this paper, not all sensor systems are subject to significant delays or sampling-time uncertainty. However, the adoption of camera-based systems that rely upon image capture and processing makes sensor systems particularly susceptible to latency and time variability. If handled naively within a highly dynamic environment, latency can cause an otherwise accurate measurement to be rejected by residual editing. If latent measurements are repeatedly processed, however, errors in the estimate can accumulate such that navigation performance is degraded, potentially leading to filter divergence and navigation system failure.
A layer of complexity is added by realizing that many navigation applications consider multi-sensor fusion, which may cause measurements to be processed in an order that is not reflective of their true sampling time. If this occurs, the measurements are then referred to as out-of-sequence measurements (OOSMs). A depiction of a scenario with latency and OOSMs is shown in Figure 1. Another factor that complicates the study and implementation of latent measurement processing is the reliance upon high-frequency, discrete inertial measurements from an inertial measurement unit (IMU) of the vehicle’s angular velocity and non-gravitational acceleration for state propagation as a model replacement. In general, inertial measurements do not directly align with collected measurements, requiring that some assumptions be made to obtain a state realization between discrete inertial measurement epochs. Latent measurement processing capabilities are of paramount importance in a real-time environment (Carpenter & D’Souza, 2018), although latency and the associated uncertainty are often not primary concerns in the development of a navigation system.
Timeline illustration for measurement sampling times and associated arrival times for navigation processing, depicting latent measurements (z1, z2, z4, z5, z6, dashed) and OOSMs (z3, solid)
Conventional navigation architectures, typically based upon adaptations to the Kalman filter (Kalman, 1960; Kalman & Bucy, 1961), cannot directly address measurements that do not coincide with the navigation epoch or contain timing uncertainties without further alteration. A straightforward and basic approach known as in-sequence processing maintains a history of measurements, states, and estimation error covariances and chronologically reprocesses all measurements after an OOSM is obtained. Bar-Shalom (2002) formulated an exact update for the case in which an OOSM occurs within the current sampling interval; to extend this method to older measurements, the application of a smoothing-like algorithm is necessary, as measurements will need to be reprocessed. Both approaches can become computationally burdensome because of the necessary reprocessing of measurements. An augmented form of the fixed-lag smoother (Simon, 2006) presented by Yoon et al. (2016) addresses the problem of OOSMs by adaptively generating the necessary nodes for OOSM processing. However, both of these approaches estimate the desired quantity using an augmented system to reprocess measurements and produce an optimal estimate that lags behind the current state. Unfortunately, this construction makes both methods nonideal for high-stakes navigation scenarios, as they require significant computational overhead and produce optimal estimates that lag behind real time. Another approach for fusing latent measurements into a general nonlinear system includes the state augmentation methods introduced by van der Merwe et al. (2004). To address uncertainty in the measurement epoch, Julier and Uhlmann (2005b) proposed a covariance union, whereas Lee and Johnson (2017) utilized multiple-model adaptive estimation to characterize uncertain measurement latencies. To produce a method that is computationally efficient and accurate, the approach developed in this paper instead treats the latency as a filter state, so that it can be directly accounted for within a typical navigation framework.
Many adaptations to the Kalman filter have been developed to extend its applicability to nonlinear systems and address complications that arise. In particular, the use of considered parameters, as initially posed in the Schmidt–Kalman filter (Schmidt, 1966), is of particular interest within this work. Under this formulation, the concept of considering (or neglecting) weakly observable parametric states is introduced: when a state is deemed as considered, its impact upon the system is still incorporated, without a direct estimation of the state. Alternatively, when a state is weakly observable and has a low impact on the system, it may be useful to instead neglect the state altogether. The use of considered states and neglected states is common for real-time navigation development, where computational resources are often significantly limited. The consider/neglect architecture outlined by DeMars and Ward (2020) is employed for this work, as it facilitates an effective method for examining of the dependence of navigation performance upon measurement timing errors. This paper examines the proposed methods by assessing the impact of jitter in terrain camera measurements on navigation performance when it is either neglected or considered.
In the remainder of this paper, a linearized approach for latent measurement processing with measurement epoch uncertainty (or jitter) is developed and examined via numerical simulation. A compact formulation is presented that makes implementation of the developed approach easier and consistent with existing navigation architectures. The extension to IMU-based navigation is then discussed, where specific methods must be developed for the necessary backpropagation or retrodiction procedure. Finally, the developed methods are examined in a lunar descent-to-landing scenario via Monte Carlo analysis for a variety of configurations.
2 METHODOLOGY
In this section, a latent measurement processing scheme is derived and presented alongside a review of the consider–neglect filtering architecture within which the approach is employed. For latent measurements, the derivation is presented in a manner such that it is easily implementable within the linearized, extended Kalman filter architecture with a generalized measurement model, allowing application within existing navigation architectures for a wide range of sensing systems.
2.1 Nonlinear Systems with Latent Measurement Processing
Let a system with a vector state x(t) evolve according to nonlinear dynamics, such that
where f(·) is a nonlinear dynamics function, dependent upon both the system state x and zero-mean white process noise v with a constant power spectral density, Q. The discrete-time evolution of the nonlinear system, describing the evolution from some time tk−1 to another tk is given by the integration
Defining the state transition matrix for x as Fx(tk, tk−1) and the estimation error as ea = a – ma, where 
is the estimate of the state a, the linearized, discrete evolution of estimation errors is
1
where Fw(tk, tk−1) is the discrete mapping of the process noise into the state and w is a zero-mean white process noise sequence with covariance Pww that serves as the discrete-time analog to the continuous-time process noise.
Suppose measurements of the state at time tℓ are given by
where h(·) is the function mapping the states into the measurement space and vℓ is a white-noise sequence with covariance Pvv,ℓ. Let the measurement time tℓ be
described in terms of the navigation time tk, latency λk ≥ 0, and jitter jk. Let us assume that λk is known and accompanies the measurement zℓ. We let the jitter instead be described as a white-noise timing error with mean mj,k and variance Pjj,k. Note that while no assumption is made here about the relative magnitude of the jitter and latency, in practice, it is nominally assumed that the latency will generally be much larger than the jitter.
Through the extended Kalman filter approximation, the expected measurement can be expressed as
2
Note that the expected measurement in this case is dependent upon the predicted state 
at the expected time point 
, where the superscript “–” denotes a priori quantities, i.e., those that exist immediately preceding a measurement update. From Equation (2), the residual 
is expressed as
3
Utilizing a first-order expansion about the prior state estimate 
at 
, Equation (3) can be approximated as
where the Jacobian of the measurement function for quantity q is as follows:
Recall that the latency is considered as known, such that eλ,k = 0. Additionally, the jitter Jacobian can be decomposed via the chain rule, such that
giving
4
where 0 is a zero vector. Note that the linearization employed in Equation (4) relies upon a good characterization of 
, which can be difficult to achieve for some systems. Therefore, the residual can be expressed as
5
The first term can be shifted to time tk via the error dynamics in Equation (1), where
6
Next, recall that the state transition matrix composes via
7
With this decomposition, Fx(tℓ, tk) can be expressed as
8
Based on a first-order truncation of the Peano–Baker series (Rugh, 1996), the state transition matrix 
can be approximated by
9
where Ax = ∂f(x, w)/∂x is the continuous-time dynamics Jacobian. Additionally, by expanding Fw(tk, tℓ) about the expected measurement-time 
via a first-order Taylor series, Fw(tk, tℓ) can be expressed as
10
With the use of Equations (8)–(10), Equation (6) can be restated by
11
Truncating Equation (11) to first order, the state error at 
can be expressed as
12
Substituting Equation (12) into Equation (5), the residual is given as
13
It is worth noting that, even for a traditionally linear system, timing errors cause nonlinear error terms to arise. Thus, the presence of timing errors elevates the estimation of a linear system such that tools and approaches from nonlinear estimation are necessary, i.e., it can be concluded that the linear system with timing errors can be considered nonlinear. Consequently, even for a linear Gaussian system, processing measurements with sampling-time uncertainty with a linear update will produce a suboptimal measurement update. Additionally, it is important to recognize that the error expressions presented in Equations (12) and (13) are used in the remainder of this section to generate a realization of uncertainty at a time where tℓ < tk; such a realization is better classified as a rewinding of time rather than a backpropagation. With a true backpropagation, uncertainty would tend to grow because of the integration of process noise, w. However, these results describe a reversal of an already forward-moving process, rather than a backward integration. Throughout the remainder of this manuscript, the process of predicting a previous state is called a “rewinding” or “retrodiction” in place of a backpropagation. Through the forward integration of process noise, the state estimation error and the process noise become correlated.
Let the update to the prior state 
with measurement zℓ be an approximate linear minimum mean-square error update of the form
14
where the expected measurement is dependent upon the prior state estimate at tk and the superscript “+” denotes a posteriori quantities, i.e., those that result from the measurement update. Notice that the update in Equation (14) occurs “through time;” the linear gain maps the innovation through the temporal displacement between the measurement and state estimate, producing an estimate at time tk. Based on the linear update defined in Equation (14) and the residual expression in Equation (13), the posterior estimation error (defined by 
) is
15
Let the posterior estimation error covariance be 
. Utilizing the posterior estimation error from Equation (15), the posterior estimation error covariance can be expressed as
16
where 
is the covariance between the state estimation error and process noise. Notice that in the development of Equation (16), it is assumed that the jitter is not correlated with any other quantities and that the measurement noise v and dynamic process noise w are uncorrelated.
When the measurement time and filter epoch coincide (tℓ = tk), as is generally accepted in idealized filtering applications, it is common to assume that the state estimation error and process noise are uncorrelated, i.e., Pxw,kℓ = 0. This assumption is easy to satisfy, as the discrete noise at that time point has not yet been integrated into the state estimate. However, if the measurement time is such that tℓ < tk, the state at tk is dependent upon the process noise sequence from tℓ that was integrated into the state moving forward through the dynamics. Through this integration, the process noise and state have necessarily become correlated for any time where tℓ < tk. From Equation (1), the correlations are
17
which allows the posterior estimation error covariance in Equation (16) to be written as
18
Notice that the expression in Equation (18) can be reduced to a form that parallels that of the standard Kalman filter with some minor modifications for measurement latency and jitter. Thus, the same optimization can be applied to obtain the linear gain that minimizes the posterior mean-square error of the estimated state, as follows:
19
where the residual covariance and cross covariance are respectively defined by
Utilizing the residual expression in Equation (13) and making use of Equation (17), the residual covariance is
20
With the residual error described in Equation (13) and assuming the state estimation error to be uncorrelated to the measurement noise, the cross covariance can be expressed as
21
where 
describes the cross-correlation between jitter and state estimation errors. It is likely that 
will initially be zero, although cross-correlations may build as measurements are processed. Given Equations (19)–(21), the linear update of Equations (14) and (18) can be applied for processing latent measurements with sampling-time uncertainty due to jitter.
2.2 A Compact Formulation
The derivation presented above can also be used to arrive at another representation that is more compact and consistent with the concatenated consider–neglect formulation of DeMars and Ward (2020). To develop this alternate definition, the estimated state is instead defined such that 
. With this concatenated representation, the residual in Equation (13) can be expressed as
where the measurement Jacobian is now expressed such that
22
and the process noise mapping takes the form
23
Following the same approach as before, with the alternatively defined state, the residual covariance of Equation (20) can be expressed as
24
With this result, the posterior estimation error is given by
with a posterior estimation error covariance stated as
25
Note that the cross covariance and linear gain are similarly given by
26
and
27
respectively. Using this formulation, correlations that occur between measurement jitter and the dynamic states of interest are no longer neglected, although the jitter and dynamic states are likely to be uncorrelated initially. This more compact formulation appears much more similar to the expected Kalman filter, although a few necessary modifications are incurred by the presence of process noise. The general procedure for processing measurements with latency and measurement jitter is summarized in Algorithm 1, where it is assumed that the rewound state estimate 
, dynamics Jacobian 
, time rate of change of the estimate 
, and process noise Jacobian 
have been pre-computed for the system of interest.
Extended Kalman Filter Update with Measurement Latency and Jitter
3 MODELING CONSIDERATIONS FOR DESCENT-TO-LANDING NAVIGATION
Within this section, a system definition is presented for a spacecraft undergoing a descent-to-landing scenario. The dynamics for a vehicle aided by a strapdown IMU are briefly introduced alongside the utilized IMU model, comprised of three linear accelerometers and three gyroscopes. Furthermore, modifications to the approaches for latent measurement processing, as introduced in Section 2, are discussed for IMU-based navigation.
3.1 Dynamics Modeling for Inertial Navigation
Let the vehicle state xk be defined as
where 
and 
are the inertial position and velocity of the IMU case frame origin, respectively, while 
is the vehicle attitude quaternion and describes the transformation from an inertially fixed frame to the IMU case frame orientation, here taken to be right-handed and vector-first. Additional parametric states of interest, such as sensor biases, are estimated and given by p. Throughout this paper, it is assumed that all of the parametric states are constant in time, i.e., ṗ = 0. This assumption serves to simplify the model and is by no means a requirement. The continuous-time dynamics for the vehicle position, velocity, and attitude are given by (DeMars, 2007)
28a
28b
28c
The non-gravitational acceleration of the vehicle is given by ang, whereas the gravitational acceleration is ag, which depends upon the position of the vehicle’s center of gravity, given by 
. The relative position of the center of gravity with respect to the IMU case frame origin is given by d, and 
is the transformation matrix representing the transformation from the IMU case to the inertially fixed frame. The pure quaternion representation of the angular velocity of the case frame with respect to the inertial frame, denoted by ωc/i, is 
. The superscripts i and c are used to denote quantities expressed in the inertial and IMU case frames, respectively. Additionally, note that ⊗ indicates quaternion multiplication and composes in the same order as transformation matrices.
By integrating Equations (28) from tk−1 to tk and linearizing about the states at tk−1, it can be shown that the forward evolution of the discrete-time dynamics are
29a
29b
29c
29d
where Δrng,k|k−1 and Δvng,k|k−1 describe the change in position and velocity, respectively, from tk−1 to tk due to non-gravitational accelerations and the total body rotation Δϕk|k−1. Additionally, Δrg,k−1 and Δvg,k−1 are the change in position and velocity due to gravity, linearized about the states at tk−1.
In navigation applications, it is commonplace to approximate the non-gravitational accelerations and the total angular velocity of the vehicle using a strapdown IMU, which generates inertial measurements as
where 
and 
are the incremental velocity and incremental angle, respectively. These quantities can be influenced by several measurement corruption phenomena resulting from a variety of sources, such as manufacturing tolerances, sensor installation, and unit degradation. Thus, a typical model used to represent inertial measurements ym of quantity y is given by
30
where b and s are bias and scale factor parameters, m and n are parameters characterizing the misalignment and nonorthogonality of the sensing axes, and w is a zero-mean, time-wise uncorrelated noise in the measurement. Note that y can be taken as either Δθc or Δvc, such that the model in Equation (30) is appropriate for both quantities. The operators used in Equation (30) are defined as follows for 
:
By applying the model given in Equation (30) and invoking the matrix inversion lemma (Golub & Van Loan, 1996), we can show that the true incremental angle and velocity can be approximated as
31a
31b
where the subscripts a and g denote quantities specific to the accelerometer and gyroscopes, separately. As a result of Equation (31), an estimate of the true incremental angle and velocity can be determined, given estimates for each of the parametric error source terms.
IMUs are typically capable of operating at much higher frequencies than those at which the navigation computer can truly process the data; thus, the non-gravitational changes in position and velocity, alongside the total incremental rotation of the body, are described instead by the result of an additional preprocessing algorithm (Zanetti et al., 2017), often housed in a separate computer. Oftentimes, the high-frequency data are pared down into lower-frequency increments via some form of coning, sculling, and scrolling correction algorithms, such as those presented by Savage (2000) and Miller (1983). The increments resulting from the preprocessing architecture are then used to propagate the states, as shown in Equations (29). For the purposes of this paper, the coning, sculling, and scrolling algorithms applied to down-sample the IMU data are those presented by Savage (2000), whereas the covariance propagation utilizes the first-order expansion error propagation described by Brouk (2019).
3.2 Latent Measurement Processing with Discrete-Time Dynamics
The development and operation of navigation algorithms within a real-time environment require care and forethought to prevent system failure and to achieve the desired precision. Figure 2 portrays some of the challenges that occur within a real-time navigation system as internal and external measurements are processed. As discussed in Section 3.1, discrete IMU measurements are typically fed at a very high frequency into a preprocessor, which then delivers a pared-down version of the measurements for state estimate propagation. In the traditional filtering sense, propagation need only occur when an external measurement is ready for processing. However, because the state estimate is used by the other vehicle systems, propagation occurs immediately after the preprocessor output is received. This allows dependent systems to obtain an up-to-date estimate rather than waiting for measurements that may never arrive. Once external measurements are received, the predicted state can be updated. Ideally, this process is synchronous with the filter cycle, at some nominal frequency. Unfortunately, this is an unrealistic scenario, as each measurement contains some time-varying latency due to several factors, including preprocessing techniques, communication complexity, and the physical state of the sensor hardware; the distribution of possible measurement times is indicated by the probability density functions in Figure 2. To accurately incorporate latent measurements with uncertainty in their sampling time, one must make adaptations to the typical filtering architecture. By utilizing the method proposed in Section 2.1, the processing of latent measurements in the real-time navigation framework is possible. However, because the state transition matrix, the time rate of the change vector, and the rewind process rely upon discrete IMU measurements, appropriate approximations for each must be determined; a method for each is proposed in the remainder of this section.
Individual IMU measurements (diamonds) are pared down into their accumulated incremental angle θi and accumulated incremental velocity vi to propagate mx,i between navigation epochs (triangles); latent measurements zi are generated at a true time (squares) but processed according to the sensor-reported time (circles).
To begin, it is first necessary to develop a method for retrodiction, rewinding the current prior 
to 
. Consider the continuous-time dynamics provided in Equation (28); by integrating from tℓ to tk with tℓ ≤ tk and linearizing about the state estimate at tk, the discrete rewinding of the position, velocity, and attitude of the vehicle is described as
32a
32b
32c
Equations (32) can also be described as an inversion of Equations (29), where the linearization has instead been performed at tk. Note that, because the parametric states have been assumed constant, pℓ = pk. It is worth pointing out that if a measurement update occurred at tk−1 with 
, one could instead process the measurement by discarding 
and propagating the previous posterior 
to 
. Thus, given low-latency measurements and the previous state estimate, rewinding the state is unnecessary, although propagation will need to occur twice more: once to tℓ and again to tk. However, if the measurement arrives at 
, the use of earlier estimates would require the reprocessing of all measurements that occurred after that estimate was generated. By rewinding from 
and processing the latent measurement through time, the need to reprocess data and buffer previous state estimates is eliminated; moreover, the predicted measurement is conditioned on the state, which contains information extracted from the set of already processed measurements.
To accomplish the retrodiction as posed in Equations (32), the internal measurements must be reprocessed to generate Δrng,k|ℓ, Δvng,k|ℓ, and Δϕk|ℓ. For this, a queue or buffer of inertial measurements must be maintained so that the latent measurement processing can occur; this is already seen in practice (Zanetti et al., 2017). Temporary memory on a flight processor is limited, requiring the specification of a buffer limit. The maximum size of the buffer N can be declared such that either the maximum expected or allowed latency can be processed; once the buffer meets this limit, the buffer will constantly replace older IMU measurements with the latest measurements. In this discussion, it is important to note that only uncorrected IMU measurements are buffered and that no parametric states have been accounted for, as described in Equation (31b); otherwise, unintended cross-correlations may go unaccounted. Once the latent measurement has been received and is ready to be processed, one can obtain the estimated incremental angles 
and velocities 
by correcting the raw quantities with the most recent a priori parametric state estimates 
; because the parametric states are assumed to be constant, the dependence upon every previously processed measurement is maintained when the most recent parametric state estimate 
is used. However, if the parametric states were not constant, they could be rewound to the desired time, starting from the most recent knowledge.
To approximate the incremental angle and velocity at a time between discrete IMU measurements, it is assumed that the angular velocity and non-gravitational acceleration vary linearly between measurements, such that linear interpolation may be applied. Additionally, to approximate Δrng,k|ℓ, Δvng,k|ℓ, and Δϕk|ℓ, an accumulation is defined such that
33
where ykℓ can be the accumulation of either error-corrected incremental angles 
or velocities 
. The n and n +1 subscripts denote measurements obtained at tn and tn+1, respectively, and are here defined to be the time steps that immediately precede and succeed the expected latent measurement time, i.e., 
. Note that both the accumulation and linear interpolation are built into Equation (33). The expression in Equation (33) relies heavily on several assumptions that are easily violated, i.e., that the rotation vectors describing the incremental angles are commutative and that the incremental velocities are expressed within the same IMU case frame realization. These assumptions can be relaxed by performing a single-measurement backpropagation or by introducing coning, sculling, and scrolling corrections, although this would introduce computational complexity into the filtering algorithms that is usually offloaded to a separate processor. Alternatively, the impact of the approximation might also be reduced through the application of a recursive rewinding, where the relinearization and accumulation are performed over several smaller intervals. Finally, it is worth noting that this assumption only affects the expected measurement calculation and comes into the posterior via the residual in Equation (14), i.e., if a higher-fidelity propagation method has been utilized to obtain the prior state estimate, that information is still contained within the posterior.
If it is assumed that the accumulation is an adequate approximation of the rotation and non-gravitational acceleration between 
and tk, a single retrodiction can be applied from tk to tℓ, as described in Equations (32). Estimates for the non-gravitational position and velocity increments, as well as the incremental rotation, from 
to tk are calculated by
34a
34b
34c
when this approximation is used. However, this can introduce significant errors as the retrodiction time interval grows because of the non-commutativity of the incremental angles and the different frames describing each incremental velocity. Thus, it is necessary to instead break this interval into separate segments and back-propagate iteratively, relinearizing at each iteration.
Recall that the discrete-time Jacobian property described in Equation (1) is equivalent to the state transition matrix of the system. Thus, both Jacobians must satisfy the composition property presented in Equation (7), such that
35
Because the state transition matrix is used to propagate the covariance forward in time from tk−1 to tk, the same approach can be used to generate the mapping from tℓ to tk; the result can then be inverted to obtain the desired state transition matrix. However, because determining the matrix inverse is often computationally expensive for large matrices and can cause issues with numerical accuracy, a different approach is sought to generate the discrete-time Jacobian. Notice that by deriving linearized error mapping for Equations (32), one can obtain an analytical description for the inverted state transition matrix. Thus, by following the approach outlined by Zanetti (2007), the discrete-time Jacobian for retrodiction can be constructed, although this is not explicitly described within this paper.
The final topic of interest within this section is the approximation of the time rate of change for the vehicle states, as described in Equations (28). From the buffered IMU measurements, the non-gravitational acceleration and angular velocity of the case frame can be approximated. Assuming that the angular velocity of the vehicle can be modeled as a linear variation between discrete IMU measurements, the linear interpolation presented earlier for Equation (33) can be applied and allows the estimated angular velocity at to be approximated as
36
where Δt = tn+1 – tn and ζ is defined in Equation (33). The non-gravitational acceleration can similarly be approximated by applying a linear interpolation to the accelerations observed via the incremental velocities, such that
37
From these definitions, approximate temporal derivatives for the vehicle states can be obtained via Equations (28), where the gravitational acceleration is evaluated at the estimated position, 
. Note, however, that because of timing uncertainties, this approach can be flawed when the vehicle rapidly performs maneuvers because of its dependence upon only two sets of inertial measurements. For larger jitter distributions with agile vehicles, a better approximation would account for the entire interval over which the measurement could have been sampled. The implementation for latent measurement processing with discrete-time dynamics is summarized in Algorithm 2, which can be used to generate the rewound state estimate 
, dynamics Jacobian 
, time rate of change of the estimate 
, and process noise mapping 
. The results of this algorithm are then used to perform the latent measurement update, as described in Section 2.
Descent-to-Landing State Retrodiction for Latent Measurement Processing
4 SIMULATION AND RESULTS
In this section, the simulated scenario is outlined, and results from several Monte Carlo simulations are examined. First, the configuration for the spacecraft of interest is presented, and the lunar descent-to-landing trajectory is discussed, including the utilized sensors and accompanying specifications. Several simulation configurations that facilitate an effective examination of the architecture proposed in Section 2.1 are identified for application to the outlined descent-to-landing scenario. Through these configurations, the impact of latent measurements with and without sampling epoch uncertainties on navigation performance is examined. The outlined configurations are as follows: configuration A serves as a baseline and examines the nominal simulation in which no latency or jitter is present, such that all measurements are synchronous to the navigation system; configuration B incorporates static latency into the terrain camera measurements without any measurement jitter; and configuration C introduces jitter into the measurement timing error. Configuration C also examines the impact when jitter is treated as a considered or neglected parameter, separately, to assess the importance of including jitter in the measurement model.
4.1 Trajectory, Sensor, and Vehicle Models
The simulated spacecraft utilizes a strapdown IMU for inertial navigation, whereas the equipped star camera, terrain camera, and three-beam slant-range/speed sensor provide external measurements; the operational specifics of each are provided in Table 1. The simulated IMU emulates the Northrop Grumman LN-200S inertial fiber-optic gyroscope,1 selected because of its medium-to high-end performance characteristics (Vandersteen et al., 2018). The operational characteristics of the simulated terrain camera mirror those of the Mars Descent Imager (MARDI) optics and detector2 utilized on the Mars Polar Lander. In an attempt to model the limitations of current image-processing techniques, we utilize a probability of 0.9 for feature detection, emulating the average true crater detection rate outlined by Woicke et al. (2018). It is also assumed that perfect feature/label associations are made between the true and estimated terrain camera measurements. While difficult to accomplish in practice, careful data curation and feature matching prior to filtering can limit the number of bad associations (Maass et al., 2020; Mourikis et al., 2009), although some recently proposed methods have attempted to solve the data association problem directly within the filter (McCabe & DeMars, 2020). Specifications for the star camera are taken to represent a combination of those for the HE-5AS Star Tracker3 and previous optical navigation literature (Zanetti, 2009). The slant-range/speed three-beam configuration and operation window are modeled after the NASA COBALT Navigation Doppler Lidar (Amzajerdian et al., 2016; Restrepo et al., 2018) with bias and noise values following the line-of-sight (LOS) model of the Mars Science Lander Terminal Descent Sensor (Pollard & Chen, 2009). Exact models for each of the sensors explored in this paper have been explicitly outlined by Ward et al. (2020).
Specifications for the Simulated External and Internal Sensors
Each trial of the Monte Carlo simulations examined in the following sections maintains the same true trajectory for the spacecraft position, velocity, and attitude. It is assumed that guidance and control inputs for the vehicle are fixed for the trajectory, resulting in the same angular velocity and non-gravitational accelerations experienced for each trial; the magnitude of each of these values is provided in Figure 3, where only active periods are shown. The resulting spacecraft trajectory has attitude and altitude profiles as shown in Figures 4 and 5. As shown in Figure 5, the vehicle begins at an initial altitude of 50 km and gradually descends to 50 m above the surface, near the landing target at Beer Crater. At approximately 24 min mission elapsed time (MET), the vehicle undergoes a brief attitude maneuver, triggering the star camera to shut down. At 25.5 min, another attitude maneuver takes place, coinciding with the beginning of the powered-descent phase of the mission. Throughout the powered descent, translational and rotational maneuvers are continuously performed for the remaining 7 min of the trajectory, and the powered descent cuts off just before landing. The operational window for each sensor is listed in Table 1 and provided as an overlay in Figures 4 and 5. It is worth pointing out that the star camera deactivation coincides with the vehicle maneuvers, whereas the terrain camera remains active until the vehicle achieves an altitude of 3 km, at which point the slant-range/speed sensors are within range to obtain measurements. Here, the simulations are specifically examined in the period following the terrain camera activation, when the impact of jitter and latency becomes relevant.
Magnitude of vehicle angular velocity (top) and acting non-gravitational accelerations (bottom)
Vehicle attitude in terms of Euler angles with sensor operational windows indicated (Ward et al., 2020)
Altitude profile with labeling for each sensor’s operational window
To enable latent measurement generation and processing, a continuous form of the discrete reference trajectory has been developed via a modified version of the interpolation method of Akima (1970), which favors data with a lower slope in the weight generation process.4 The interpolation method utilizes a piece-wise, cubic-polynomial representation of the discrete data that is restricted to have a continuous first derivative. The resulting polynomials have a smoother, more natural appearance and avoid overshooting. By applying this interpolation to the trajectory, one can extract non-gravitational accelerations, angular velocities, and the dynamic quantities describing the vehicle behavior at any point along the trajectory.
The lunar ground-track presented in Figure 6 displays the features contained within the onboard map (based upon the database constructed by Robbins (2019)), highlighting those visible to the terrain camera along the trajectory. Measurements are generated from the visible feature set, where the number of detected features is limited to 50 individual detections for any single measurement epoch. The elevation map utilized in this work is the SELENE digital elevation model (DEM) (Barker et al., 2016), which has precision down to 1/512-degree grid spacing between data points (i.e., 512 pixels per degree). To generate this DEM, images from the SELENE terrain camera were co-registered with altimetry data collected by the Lunar Orbiter Laser Altimeter (LOLA) Lunar Reconnaissance Orbiter (Smith et al., 2010), producing a DEM product with improved accuracy over those that use only LOLA data. In this paper, it is assumed that there is no mismatch between the onboard map and the map used to generate true measurements; a case of model mismatch in the terrain has previously been examined by Ward et al. (2020).
True vehicle ground-track with the terrain camera feature set
To initialize the filter states for each Monte Carlo trial, estimates for the dynamic states (position, velocity, and attitude) are sampled from a Gaussian distribution with the initial truth vector as the mean and 1σ uncertainties of 1,000 m for position states, 0.1 m/s for velocity states, and 100 arcsec for attitude states, where the filter states are assumed to be initially uncorrelated. The true parametric states are sampled from a Gaussian distribution with zero mean and covariance defined according to the sensor specifications listed in Table 1.
Many factors, such as linearization, contribute to errors in the filtering solution, even if care is taken to develop high-fidelity dynamics and sensor models. To mitigate the effects of these extraneous and sometimes unpredictable error sources on the navigation filter performance, several tools are employed within the filtering architecture. First and foremost, residual editing (Carpenter & D’Souza, 2018) is applied with a probability gate corresponding to 99.0% for each sensor, while measurement underweighting (Zanetti et al., 2010) is applied with a factor of 5/7; both of these techniques seek to maintain conservative estimates. Tuning noise is constantly injected into the filtering solution via the dynamics states and further promotes conservative performance: parametric states are only influenced through correlations. The injected uncertainty has statistics (1σ) given by 0.1 μm, 1 μm/s, and 0.1 arcsec in the position, velocity, and attitude channels, respectively, and is introduced during the filter propagation. The specific values have been roughly selected such that the standard deviations of the Monte Carlo errors are roughly encompassed by the filter standard deviation in the nominal filter configuration, which is examined in Section 4.2. The slightly conservative filter estimate provided through the tuning noise selection helps prevent filter performance degradation in the presence of other, unmodeled sources of error.
4.2 Configuration A: Nominal Simulation
The nominal simulation serves as a baseline against which to compare the other configurations. Within this simulation, all measurements are instantaneously sampled, communicated, and processed by the navigation algorithms, which run at a fixed frequency of 10 Hz. The IMU preprocessor obtains inertial measurements and produces a corrected output at 10 Hz for propagation, using coning, sculling, and scrolling algorithms presented by Savage (2000) and Jacobians described by Brouk (2019). The filter utilizes a consider–neglect formulation of the multiplicative extended Kalman filter (Crassidis & Junkins, 2004), as described in Section 2, with system dynamics and an IMU model as presented in Section 3.1. For the purpose of simplification, IMU scale factors, misalignments, and nonorthogonalities have been eliminated from both the true and simulated system models. Simulation results from 1,000 Monte Carlo trials are shown in Figure 7(a), where the root-sum-square (RSS) of the Monte Carlo error covariance is compared with the RSS of the averaged filter statistics for vehicle position, velocity, and attitude states. The results shown focus on the trajectory from the terrain camera initialization to the end of the powered descent because the analyses in the following sections will focus primarily on the impact of measurement latency and jitter in the terrain camera measurements. Note that there appears to be a large non-zero mean error in Figure 7(a). This is a result of the trajectory beginning with 10 min of propagation before TRN acquisition; thus, the uncertainty at this point is initially very large. With such a large uncertainty at acquisition, this seemingly large bias is comparatively zero-mean; this same observation can be made in the other simulation configurations as well. In Figure 8, the filter consistency can be examined by comparing the quantities directly as a ratio, i.e., the ratio shown compares the averaged filter covariance with that of the Monte Carlo error statistics as
38
where x is taken to be any state or subset of states of interest. A ratio of unity from Equation (38) indicates perfect agreement, whereas a ratio less than unity indicates a filter that is overconfident and a ratio greater than unity suggests that the filter is conservative in its estimate. It is, however, necessary to point out that only autocorrelation is considered for this consistency measure. Note that agreement is seen throughout the examined period: the Monte Carlo statistics approximately overlap with the averaged filter statistics, giving a consistency ratio of approximately unity for much of the trajectory. Additionally, following the initial acquisition, the mean error appears to be approximately zero, indicating that the filter estimate is unbiased.
Monte Carlo results showing sample mean error (RSS), sample covariance (RSS), and averaged filter covariance (RSS) for both the (a) nominal and (b) 2.5-s latent terrain camera measurement configurations
Ratio of the averaged filter covariance (RSS) to the Monte Carlo sample covariance (RSS), comparing the nominal simulation (blue, solid) and the outlined 2.5-s latent measurement processing (orange, dashed)
4.3 Configuration B: Latent Terrain Camera Measurements
In the second simulation considered for this paper, the simulation architecture and filter described in Section 4.2 are augmented with the developed latent measurement handling approach, applied to the terrain camera sensor. The terrain camera measurements are ingested by the filter via an accompanying time-tag as reported by the sensor. Terrain camera measurements arrive 2.5 s latent and are intended to model a sensor that requires a significant overhead, due to capture, processing, and data transport. In this specific scenario, the latency is static and known exactly without any measurement jitter. To process latent measurements, the methods described in Sections 2.1 and 3.2 are utilized, where multistep retrodiction is performed at the operational speed of the IMU. Simulation results from 1,000 Monte Carlo trials with latent terrain camera measurements are shown in Figures 7(b) and 8. When the results of the simulation with latent terrain camera measurements in Figure 7(b) are compared with those of the nominal simulation without latency in Figure 7(a), it is clear that some performance degradation has occurred; the filter’s velocity estimate is briefly overconfident, recovering consistency at approximately 13 min MET. Throughout the entire interval, the position uncertainty also appears to be slightly undercharacterized. In Figure 8, the consistency ratio of Equation (38) is examined. From the trends present in Figure 8, the earlier observations are confirmed: consistency in the velocity seems to degrade until approximately 13 min, at which point the estimate begins to recover consistency. The position estimate, however, appears to be consistently undercharacterized throughout the entire trajectory, with gradual degradation from 19 min until 27 min MET. Note that consistency for the attitude estimate begins to degrade following the attitude maneuver at 24 min MET, as shown in Figures 3 and 4, although this is likely due to the limited characterization of the retrodiction, which rewinds the state to the expected time of the measurement without considering any of the vehicle dynamics surrounding the expected measurement time.
4.4 Configuration C: Latent Terrain Camera Measurements with Uncertain Epochs
To introduce and examine the combined effects of measurement latency and jitter in a descent-to-landing scenario, the simulation examined in Section 4.3 is further expanded to include jitter, modeled as a Gaussian white-noise error in the measurement sampling time, tℓ. The terrain camera system described in the previous sections is retained, where measurements are received with a 2.5-s latency. Additionally, the impact of a model mismatch in the navigation system is examined through the use of considered and neglected parameter designations; when jitter is considered, the developed approach is employed, and when jitter is neglected, the associated uncertainty is eliminated entirely from the filter model. To examine the impact on navigation performance with variable jitter, several cases are considered. Configurations in which jitter is modeled as a considered parameter have standard deviations of 1 ms, 2 ms, 5 ms, and 10 ms. For the model-mismatch case, where jitter is treated as a neglected parameter, the Gaussian jitter distribution is varied with standard deviations of 2 ms, 5 ms, and 10 ms.
Results for each configuration are provided in Figures 9, whereas Figures 10 and 11 show a direct comparison of certain configurations of interest. By examining Figures 9, it can be deduced that neglecting jitter generally produces degraded navigation performance, tending to more significantly impact position and attitude estimates. Figure 10 shows a comparison of the filter’s consistency with Monte Carlo RSS results when the jitter is treated as a considered parameter with varying error standard deviation. In general, it is clear that increased jitter degrades estimation performance, likely resulting from the strong nonlinearities present in the descent-to-landing scenario combined with high-impact measurements such as the terrain camera measurements. Following the initial acquisition of terrain camera imagery at approximately 10 min, a reduction in consistency is observed until approximately 14 min. An examination of individual trial results shows that this first update is particularly sensitive to terrain camera jitter, often causing an increase in the estimation error with larger jitter errors. During the few minutes following the initial terrain camera measurements, this increased error is typically reduced, and the filter regains consistency with additional measurements. Fortunately, each configuration is capable of returning to approximately consistent or even conservative performance before landing occurs. A comparison of the results for the configuration with 10-ms jitter, shown in Figure 9(e), and the nominal simulation in Figure 7(a) shows that the final averaged filter covariance values (3σ, RSS) are 20.8 m and 19.8 m, respectively, suggesting that the use of slant-range and slant-speed sensors at the end of the simulation help to regain precision for final landing.
Monte Carlo results showing sample mean error (RSS), sample covariance (RSS), and averaged filter covariance (RSS) for 2.5-s latent terrain camera measurements with 2-ms jitter (top row), 5-ms jitter (middle row), and 10-ms jitter (bottom row), where jitter is treated as either a considered (left column) or neglected (right column) parameter (a) 2-ms jitter (considered) (b) 2-ms jitter (neglected) (c) 5-ms jitter (considered) (d) 5-ms jitter (neglected) (e) 10-ms jitter (considered) (f) 10-ms jitter (neglected)
Ratio of the averaged filter to the Monte Carlo sample covariance for cases without jitter and for configurations with 1-ms, 2-ms, 5-ms, and 10-ms jitter included for terrain camera measurements, treated as a considered parameter
Ratio of the averaged filter to the Monte Carlo sample covariance with 2-ms, 5-ms, and 10-ms jitter (1σ) included for terrain camera measurements, with jitter treated as a considered or neglected parameter
Figure 11 provides a comparison of filter consistency for cases in which jitter is present but treated as either a considered or neglected parameter. Given these results, it is clear that 2 ms of jitter (1σ) has little impact on navigation performance. However, for 5 ms of jitter (1σ), the difference between the considered and neglected cases no longer seems to be negligible. At 10 ms of jitter (1σ), navigation performance begins to suffer in both cases, although it is still better to consider jitter in the terrain camera measurements. In further experimentation beyond the results shown, consistency throughout the period during which the terrain camera is active continues to steadily decrease as the jitter magnitude is increased. Thus, new methods or approximations should be sought for sensors exhibiting sampling-time uncertainties larger than 10 ms (1σ). In summary, neglecting measurement jitter can significantly deteriorate a filter’s estimation capability, although by treating jitter as a considered parameter, the proposed methods can recover some performance. In general, it appears that the most significant impact on navigation performance due to jitter being neglected is realized in the attitude estimate after a large maneuver, while position and velocity estimation exhibit degraded performance throughout, with a somewhat less significant reduction in consistency immediately following maneuvers. However, this impact from maneuvers may be due to the limited characterization of the retrodiction, which rewinds the state to the expected time of the measurement without considering the vehicle dynamics surrounding the expected measurement time.
4.5 Error Budget and Sensitivity to Uncertain Measurement Epochs
Based upon linear covariance analysis, error budgets and sensitivity analyses are common tools for examining the impact of certain states on the navigation solution; the basic procedures for error budget construction and sensitivity analyses have been described by Gelb (1974). Following the development of an error budget, it is possible to make informed decisions about sensor requirements and identify which states are most sensitive to a given error source. To better characterize the impact of terrain camera measurement jitter on vehicle state uncertainty, this section examines an error budget focusing on the terrain camera and initial vehicle state error sources. In particular, the initial position, velocity, and attitude errors for the vehicle are examined alongside the terrain camera jitter, noise, and bias errors, while additional error sources are all grouped together (represented as a solid gray region in Figure 12).
Position (top), velocity (middle), and attitude (bottom) uncertainty contributions separated into initial vehicle state uncertainties and terrain camera errors
The position covariance breakdown in Figure 12 clearly shows that jitter does not dominate the uncertainty. However, during coast periods, the contribution of terrain camera jitter tends to increase, reaching a peak of 15.7% at 24 min. Following the attitude maneuver at 24 min, the impact of jitter relative to the other sources begins to grow again, reaching a peak of 22.8% at 29.2 min before decreasing again. The same type of behavior is observed in the velocity (also shown in Figure 12), although the growth following the attitude maneuver at 24 min peaks at approximately 28 min with a value of 27.9%. When the attitude trend is examined, it is clear that the gradual growth in contribution resulting from terrain camera jitter, shown in position and velocity, is not present. However, following the attitude maneuver at 24 min, the contribution of the terrain camera jitter to the attitude uncertainty grows, peaking at 24.4 min with a value of 13.8%. Note that this corresponds to the loss of star camera measurements at 24 min. Examination of each trend in Figure 12 together shows that the terrain camera jitter errors first impact the attitude at 24 min, followed by velocity at 25 min, and position at approximately 27 min, with the peaks occurring sequentially after 24 min. For brevity, further figures illustrating the impact on other estimated states are omitted, although many show little dependence on the terrain camera jitter. However, biases in the accelerometer and terrain camera show a maximum contribution from terrain camera jitter of 12.1% and 15.2%, respectively, between 27 and 32 min MET. Both peak at approximately 29.2 min MET and steadily decrease for the remainder of the trajectory.
For the sensitivity analysis shown in Figure 13, six different times are selected for examination: (1) 9.76 min, corresponding to immediately after the first set of terrain camera measurements, (2) 10.05 min, corresponding to 10 terrain camera updates after the first terrain camera measurements, while the filter is rapidly converging, (3) 14 min, corresponding to a few minutes after terrain camera initialization, where the filter convergence rate has slowed, (4) 23.73 min, corresponding to a few terrain camera measurements preceding the large attitude maneuver at approximately 24 min, (5) 30.63 min, which precedes the large translational maneuver defining the beginning of powered descent, and (6) 32.94 min, corresponding to the final simulation time. For this examination, the baseline terrain camera jitter is 5 ms, making the 10-ms case examined in previous sections the upper limit of the scaling shown. Figure 13 shows that the percentage contribution of jitter to position uncertainties is most significant, although this primarily occurs for times when errors are allowed to accumulate between maneuvers and during stretches of time when the terrain camera is the primary sensor. Examining the sensitivity at 23.73 min and 30.63 min clearly shows that position is especially sensitive, although the ratio appears to plateau under 30%; this result suggests that jitter under 2 ms can likely be ignored. Additionally, the rapid increase in position ratio above 30% indicates that terrain camera measurement jitter should be mitigated when possible. Examination of the velocity results in Figure 13 suggests that velocity errors are much less sensitive to the terrain camera jitter, although the velocity at 23.73 min seems to show a sensitivity similar to that of position at the same time. Finally, the attitude ratio shows very little variation with terrain camera jitter.
Sensitivity to terrain camera jitter in position, velocity, and attitude, generated with a 5-ms jitter (1σ) baseline
5 CONCLUSION
A method for handling latent measurements with sampling-time errors due to jitter has been presented and described within this paper. Adaptations to the method, such as extending the method for implementation within a system reliant upon inertial measurements, are also discussed. The implementation is applied to a descent-to-landing scenario and examined via Monte Carlo simulation and analysis. In the simulation, the impact of latency and jitter on the terrain camera system is examined. As shown in Section 4.3, the estimation error covariance can be reasonably recovered for large latencies. When examining the introduction of jitter in Section 4.4, it becomes apparent that degraded performance can be expected, although the developed methods help to mitigate the performance degradation due to measurement sampling-time errors. Following the sensitivity analysis and error budget, it is also clear that the vehicle state uncertainty is most strongly impacted by terrain camera jitter during periods of little dynamic activity where the error is allowed to accumulate. In systems with large jitter errors, the impact will be much more significant and could be a subject of further study. However, following these results, it is clear that investing in improved sensing and communication technologies with accurate timing characterizations is much more important than pursuing improved processing strategies.
HOW TO CITE THIS ARTICLE
Brouk, J. D., & DeMars, K. J. (2024). Kalman filtering with uncertain and asynchronous measurement epochs. NAVIGATION, 71(3). https://doi.org/10.33012/navi.652
AUTHOR CONTRIBUTIONS
Both authors worked to develop and refine the methodology presented in this paper. James Brouk performed the simulations and wrote the manuscript. Dr. Kyle DeMars provided guidance and feedback throughout the development of the methodology and manuscript.
CONFLICT OF INTEREST
The authors declare no potential conflicts of interest.
ACKNOWLEDGMENTS
This research was supported by the National Aeronautics and Space Administration (NASA) Space Technology Graduate Research Opportunity and NASA grant number 80NSSC19M0208. The authors would also like to thank Gunner Fritsch for his insights on error budget development and sensitivity analyses.
Footnotes
↵1 Northrop Grumman LN200S specifications are available from the manufacturer here.
↵2 MARDI specifications are available from Malin Space Science Systems here.
↵3 Characteristics of the HE-5AS Star Tracker are available here.
↵4 Documentation describing the MATLAB implementation of the modified Akima interpolation method implementation is available here.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.