Anthea Comellini
email: anthea.comellini@isae.fr

Davide Casu
email: davide.casu@thalesaleniaspace.com

Emmanuel Zenou
email: emmanuel.zenou@isae.fr

Vincent Dubanchet
email: vincent.dubanchet@thalesaleniaspace.com

Christine Espinosa
email: christine.espinosa@isae.fr

Incorporating delayed and multi-rate measurements in navigation filter for autonomous space rendezvous

I. Introduction

Autonomous rendezvous (RDV) and capture are key capabilities to answer main challenges in space engineering, such as In-Orbit-Servicing and Active Debris Removal. When the target does not assist the chaser in acquisition, track and rendezvous operations, it is referred to as non-cooperative [START_REF] Confers | Satellite servicing safety framework technical and operational guidance document[END_REF], meaning that the chaser has to estimate on board the target state for sake of autonomy. Inexpensive camera sensors, coupled with image processing (IP) and Computer Vision (CV) algorithms can provide cost effective and accurate measurements of relative pose (i.e., position and attitude) of the target. These tracking algorithms can have a relatively high latency time. This results in a delay between the time instant of data acquisition and the time instant when the processed measurements are available and ready to be fused into the navigation filter. The navigation filter will therefore need to merge infrequent and delayed measurements. Since the tracking can be provided by additional sensors and algorithms with different latency time, the filter must be able to fuse multi-rate measurements. While slow measurements are available after a certain delay, fast measurements (i.e., referred to as interim measurements, IM) are available at a higher rate and processed almost instantaneously and have to be fused within the delay period.

The problem of delay management in space applications and more precisely in space RDV scenarios has been sporadically assessed [START_REF] Zhou | A new approach for teleoperation rendezvous and docking with time delay[END_REF][START_REF] Rems | Rendezvous GNC-system for autonomous orbital servicing of uncooperative targets[END_REF][START_REF] Benninghoff | Development and hardware-in-the-loop test of a guidance, navigation and control system for on-orbit servicing[END_REF]. These works propose delay management techniques for the estimation of the chaser-target relative translational dynamics -which is described by the Clohessy-Wiltshire-Hill (CWH, [START_REF] Fehse | Automated rendezvous and docking of spacecraft[END_REF][START_REF] Curtis | Orbital mechanics for engineering students[END_REF]) equations-and are suitable only for linear systems where no IM are considered. No solutions have already been presented for the estimation of target rotational dynamics, which is, in addition, represented by a non-linear model. The current recommended best practice for S/C attitude estimation is the M-EKF (Multiplicative-Extended Kalman Filter, [START_REF] Markley | Attitude Error Representations for Kalman Filtering[END_REF]), where both attitude quaternion and rotation rate are measured and no consistent delay affects the measurements. However vision-based navigation algorithms usually provide only pose measurements, and these measurements can be affected by a substantial delay. Some RDV operations require the synchronization of chaser motion with target motion, implying the need of knowing also target velocity and rotation rate. In the case of high rotation rates typical of a tumbling object, the rotation rate can be estimated using a dynamic filter including the object angular momentum equation instead of a kinematic one. For this reason, the target rotational state estimation will be formulated as an Additive Extended Kalman Filter (A-EKF).

In fields such as automation industry, there is a vast literature on Kalman Filter (KF) with delayed and multi-rate measurements. In [START_REF] Gopalakrishnan | Incorporating delayed and infrequent measurements in Extended Kalman Filter based nonlinear state estimation[END_REF] methods are classified into two main families: "state augmentation approaches" and methods which fuse the measurements on arrival. State augmentation methods rely on augmenting the current state with appropriate past information required to fuse the delayed measurements. These methods, such as the fixed-lag smoothing [START_REF] Simon | Optimal state estimation: Kalman, H infinity, and nonlinear approaches[END_REF][START_REF] Gelb | Applied optimal estimation[END_REF], provide optimal estimates and can be extended to other filters (e.g., Particle Filter, Unscented KF), but are suitable only for fixed-delay measurements. Moreover the size of the system increases as the delay increases, leading thus to a proportional increase of the computational load. These methods are therefore suitable for applications where delays consist in small number of samples (e.g., [START_REF] Hsiao | Robust Kalman Filter Synthesis for Uncertain Multiple Time-Delay Stochastic Systems[END_REF][START_REF] Kaszkurewicz | Discrete-time state estimation with two counters and measurement delay[END_REF]) or applications where the computational burden is not an issue, such as industrial process control, but appear to be inapplicable to space RDV due to the limited on-board computational resources of S/C. On the other side, methods that fuse the delayed measurements on arrival can handle large and variable delays with a reasonably low computational load, while granting the optimality of the estimation under certain intervals and conditions. These methods can incorporate interim measurements while merging the slowest ones, which makes them suitable for autonomous navigation. Within this family of delay management techniques, one method always provides an optimal estimate even in presence of IM: the Filter Recalculation method (FR method) [START_REF] Prasad | Product property and production rate control of styrene polymerization[END_REF][START_REF] Li | On-line state and parameter estimation of EPDM polymerization reactors using a hierarchical extended Kalman filter[END_REF]. It can be applied also to non-linear systems without any loss of optimality at the expense of a high computational load. On the other side, the method that provides the best trade-off between optimality and computational burden is the Larsen's method [START_REF] Larsen | Incorporation of time delayed measurements in a discrete-time Kalman filter[END_REF]. This method has been theorized for linear systems but can be extended also to non-linear systems with the introduction of certain approximations (Sec. III.B.2). This Note offers a comparison of these two delay management techniques, which have never been applied to the vision-based autonomous RDV navigation problem, and especially to the estimation of the highly non-linear target relative rotational dynamics. The Note is structured as follows: in Sec.II the techniques are formalized, in Sec.III the implementation of these techniques is described towards the vision-based autonomous RDV problem, in Sec.IV the performance of the filters is investigated under different sources of uncertainties, and in Sec.V the conclusions are drown.

II. Filter Equations

In this section the methods are implemented on a linear time-discrete system for sake of clarity, but will be extended to a non-linear continuous system in Sec.III. A linear discrete system observed by non-delayed measurements, where the process noise and the measurement noise are zero-mean Gaussian white-noise processes, can be put in state-space form as follows [START_REF] Gibbs | Advanced Kalman filtering, least-squares and modeling: a practical handbook[END_REF]:

       = -1 + + = + , with =        0 ≠ = , =        0 ≠ = (1)
The associated KF is divided in: prediction of the a priori estimate of the state and the state error covariance matrix (Eq.( 2)); computation of the optimal gain minimizing the a posteriori estimate of the state error covariance (Eq.( 3)); update of state and covariance matrix (Eq.( 4)).

prediction

       ˆ | -1 = ˆ -1| -1 + | -1 = -1| -1 + (2)
gain computation = | -1 ( | -1 + ) -1 (3) update        ˆ | = ˆ | -1 + ( -ˆ | -1 ) | = ( - ) -1| -1 (4)
When delayed measurements are presents, at instant the system in Eq.( 1) receives a delayed measurement corresponding to time instant ( = -, number of delay samples), such that:

* = * + , with =        0 ≠ * = (5)
In such a case, Eq.( 3) is no more optimal and a new solution has to be found in order to compute the best estimates ˆ | , * and | , * which take into account the contribution of * .

A. Filter Recalculation method

The FR method consists of going back to the time step when the delayed measurement was taken, incorporating the measurement and recomputing the entire trajectory of the state until the current step. In so doing, the whole estimate time history will be optimal. The estimation is made as if two filters were employed: a principal one, which operates at constant rate by processing fast measurements , and a second one, which is activated any time a delayed (i.e., slow and infrequent) measurement * arrives. The filters operate as follows. At instant a slow measurement is taken and ˆ | -1 , | -1 are stored. For all the time steps

+ ( ∈ [1, - 1]
) the principal filter processes fast measurements + as in a KF (Eqs. ( 2),(3),( 4)), and measurements + and inputs + are stored.

At = + the slow measurement * and its corresponding covariance * become available: the secondary filter is activated, goes back to instant and computes the optimal update using the full measurements vector ˜ = [ , * ] , where are the fast measurements and * the slow ones:

               ˜ = | -1 ˜ ( ˜ | -1 ˜ + ˜ ) -1 ˆ | = ˆ | -1 + ˜ ( ˜ -˜ ˆ | -1 ) | = ( -˜ ˜ ) | -1 , with ˜ =         *         , ˜ =         ∅ ∅ *         (6)
The optimal estimates ˆ | and | are then propagated by the secondary filter from instant + 1 to instant = + according to Eqs.( 2),(3),(4), thus clarifying the need of storing the values of IM and inputs. Once the loop has reached instant = + , the filter has provided an optimal estimate of the current state and of the state error covariance matrix. The FR method can be extended to variable delays and provides an optimal estimate in presence of IM even for non-linear systems.

More precisely, in the case of non-linear systems, the loss of optimality is introduced by the use of the Extended KF, which is intrinsically sub-optimal due to the linearization of state and measurement equations, and not by the use of the FR method.

B. Larsen's method

The Larsen's (or extrapolation) method was proposed in [START_REF] Larsen | Incorporation of time delayed measurements in a discrete-time Kalman filter[END_REF] as an improvement of Alexander's method ( [START_REF] Alexander | State estimation for distributed systems with sensing delay[END_REF]) for delay management in discrete linear system. These methods rely on the computation, throughout the delay period, of a correction term to add to the filter estimate when the delayed measurement becomes available. The main difference between Alexander's and Larsen's methods is that the latter does not need to know, at instant , the covariance matrix * nor the measurement sensitivity matrix * : these matrices are supposed to become available at instant together with the delayed measurement * . Larsen's method is therefore suitable for systems relying on measurements processed by IP-CV algorithms, as these algorithms usually process * together with * . As for the FR method, at time , when a new slow measurement is acquired, ˆ | -1 and | -1 are stored. At

each instant + ( ∈ [1, - 1]
) the classic KF structure in Eqs. ( 2),( 3),( 4) is applied using fast measurements. Moreover, the term + = ( -+ + ) + + -1 is computed, with = . At time instant = + , * , * and * become available. The filter firstly computes the gain and the updates as in Eqs.( 3),(4) using fast measurements . The final correction term * will be [START_REF] Larsen | Incorporation of time delayed measurements in a discrete-time Kalman filter[END_REF]:

* = + = =1 ( -+ + ) + (7)
Then an extrapolated measurement is computed to derive a representation of * at instant :

= * - * ˆ | -1 + * ˆ | -1 (8)
In [START_REF] Larsen | Incorporation of time delayed measurements in a discrete-time Kalman filter[END_REF] Larsen provides the demonstration of the computation of the optimal gain * and the resulting state and covariance updates, which are:

               * = * | -1 * ( * | -1 * + * ) -1 ˆ | , * = ˆ | + * ( - * ˆ | ) | , * = | - * * | -1 * (9)
where the Kalman gain * is actually the Kalman gain * (i.e., the gain that would have been computed if the measurement * had become available at instant ) pre-multiplied by the Larsen correction term * . In the presence of IM this method performs sub-optimally: at each interim step, the gain + is computed using a covariance matrix

+ | + -1 that is not optimal because it
has not yet taken into account the contribution of * . In any case, Larsen's method always requires only two matrix multiplications at each time instant and the storage of two variables any time a slow measurement is acquired, with no need of storing IM and inputs. As the FR method, Larsen's method can be extended to variable delays which are not known a priori.

C. No interim measurements case

In the absence of IM and for a linear system, FR and Larsen's methods give the same estimates.

Let us assume that the delayed measurement * and its corresponding sensitive matrix and noise covariance matrix are available at instant . The update at instant is then given by:

         ˆ | = ˆ | -1 + * ( * - * ˆ | -1 ) | = ( - * * ) | -1 with * = | -1 * ( * | -1 * + * ) -1 (10)
For all the time steps = + ( ∈ [1, -1]) state and error covariance matrix evolve in open-loop:

         ˆ | = ˆ | -1 = ˆ -1| -1 + | = | -1 = -1| -1 + (11)
After loops, at = + , the estimated state ˆ + | + will be:

ˆ + | + = =1 + ˆ | + =1 = +1 + + + (12)
Recalling the expression of ˆ | in Eq.( 10) and calling * the residual * - * ˆ | -1 , Eq.( 12) becomes:

ˆ + | + =        =1 + ˆ | -1 +       =1 = +1 + + +              + =1 + * * (13)
where the expression inside the curly brackets corresponds to the open-loop evolution of the predicted state ˆ | -1 from to + , and the term =1 + is the correction term * of Eq.( 7) in the absence of IM, and therefore in the absence of interim Kalman gains. The same demonstration can be done for the state error covariance matrix. At = + , the estimate will be:

+ | + = =1 + | =1 + + =1 = +1 + + =1 + =        =1 + | -1 =1 + +       =1 = +1 + + =1 +              - =1 + * * | -1 =1 + ( 14
)
where the term in the curly brackets is the open loop evolution of | -1 and the term outside the curly brackets is equal to Larsen covariance update in Eq.( 9). These demonstrations show how, in the particular case of the absence of IM, using the delayed measurement at instant and then propagating the a posteriori estimate ˆ | (i.e., FR method) is equivalent to propagate the a priori estimate (prediction) ˆ | -1 and to perform the update at = + by pre-multiplying the update * * by a correction factor, which is Larsen's correction factor. Larsen's method exploits the superposition property of linear systems to project in the "future" the update * * , which can be seen as a Δ that is propagated through the same transformation of ˆ | -1 . In such a case, Larsen's method is preferable since it always requires a lower amount of computation. FR method should be selected only if the whole optimal time history of the estimate from to + needs to be known. These considerations are valid only for linear systems: non-linear systems cannot exploit the superposition property and the approximated transition matrix has to be used instead of .

III. Application to the space rendezvous problem

When expressed at its Center-of-Mass (CoM), the motion of a S/C can be decoupled between the translational motion of its CoM and its rotational motion. Therefore the navigation filters for translational and rotational dynamics will be completely decoupled. Actually, when the relative translational dynamics is modeled according to the CWH equations, a small coupling between translational and rotational motion exist, but can be neglected as long as the chaser is controlled with respect to a reference frame that originates in its CoM [START_REF] Fehse | Automated rendezvous and docking of spacecraft[END_REF].

A. Translational dynamics in space rendezvous

In a RDV the relative distance chaser-target is much smaller than the distance planet-target, allowing the introduction of some simplifications to derive the dynamic model which describes the relative motion of the chaser with respect to the target. This model, in the case of a circular orbit, can be put under the form of a system of linear differential equations called Hill's equations [START_REF] Fehse | Automated rendezvous and docking of spacecraft[END_REF]:

               -3 2 -2 = / + 2 = / + 2 = / (15)
where , , are the control forces acting on the chaser CoM, is the chaser's mass and is the target's orbit angular rate. The relative position is expressed in target (Local Orbital Frame), according to the convention used in [START_REF] Curtis | Orbital mechanics for engineering students[END_REF] ( axis along the radial Earth-target, axis along target orbit angular momentum and axis completing the right-handed trihedron). The homogeneous solution of the system in Eq.( 15), which leads to the analytical computation of the transition matrix Φ( ) of the time-continuous linear state-space representation corresponding to Eq.( 15), is known under the name of Clohessy-Wiltshire equations. The CWH transition matrix is valid only for circular orbits, even if the analytical computation of a transition matrix valid for elliptical orbits is described in [START_REF] Yamanaka | New state transition matrix for relative motion on an arbitrary elliptical orbit[END_REF]. Given the transition matrix Φ( ) of the time-continuous linear state-space representation corresponding to Eq.( 15), the time-discrete state matrices , can be computed according to:

= = Φ( ) , = ∫ 0 ( 16
)
where , are the time-continuous state matrices associated to the system in Eq.( 15) and is the time step size of the time-discrete system. The result of the computation can be found in [START_REF] Fehse | Automated rendezvous and docking of spacecraft[END_REF]. The system can be therefore written in the form of Eq.( 1) and both Larsen's and FR methods can be easily implemented. In the case of vision-based navigation, the CV algorithms compute a relative target position in camera reference frame. Assuming that the rotation quaternion from camera to chaser reference frame is known by on-ground calibration, it will be necessary to know at each instant the quaternion ℎ- , in order to rotate the measurements (as well as chaser control forces) from chaser ℎ reference frame into target coordinates before using them in the KF. ℎ- can be computed knowing chaser absolute attitude quaternion, denoted -ℎ with the inertial reference frame, and the rotation quaternion from inertial reference frame to target -. The chaser inertial attitude -ℎ is estimated by the usual Attitude and Orbit Control System. Besides, the is related to target orbital parameters, and therefore to its absolute velocity and position. At each time step, the KF-estimated relative position will be added to the absolute chaser position (i.e., whose estimation can rely on GNSS and accelerometers measurements), resulting in an absolute translational target state from which an updated estimation of - can be derived.

B. Rotational dynamics in space rendezvous

As anticipated, certain close proximity operations require the knowledge of the complete rotational state of the target (i.e., attitude quaternion -and rotation rate -). The estimation of the absolute rotational dynamics of the chaser is not considered in this work. In this Note rotations are described using quaternions according to Hamilton convention [START_REF] Sola | Quaternion kinematics for the error-state KF[END_REF]. The Target rotational dynamics will be modeled according to the following prediction model:

               - = 1 2 -⊗       0 -       - = --1 -× - (17)
where is the inertia matrix of the target at its CoM. The orbital disturbance torques that affects the target dynamics will be modeled in the system as process noises. This second order system is formulated as an A-EKF, in order to ensure observability of the rotation rate from attitude measurements. Moreover, the application of Larsen's method to the M-EKF is not straightforward due to the presence of a multiplicative update. The M-EKF formulation was developed in order to avoid the ill-conditioning problems that could appear in the state error covariance matrix when the quaternion normalization is forced in the A-EKF update step [START_REF] Markley | Fundamentals of spacecraft attitude determination and control[END_REF]. An in-depth study is carried out in [START_REF] Carmi | On the covariance singularity of quaternion estimators[END_REF] to understand whether or not such a constraint leads to an ill-conditioned matrix. The work provides the mathematical demonstration that the quaternion estimation error covariance matrix does lead to ill-conditioning problem as its trace tends to zero, but also that this property is not inherited by its corresponding second order Taylor approximation, as the one computed by the A-EKF. This theoretical demonstration explains how several works relying on an A-EKF attitude estimation (e.g., such as [START_REF] Bar-Itzhack | Attitude determination from vector observations: Quaternion estimation[END_REF]) have never shown ill-conditioning problems. Moreover [START_REF] Carmi | On the covariance singularity of quaternion estimators[END_REF] has provided a practical demonstration that the covariance matrix computed by the A-EKF is well conditioned even in static and noise-free attitude estimation problem, contrary to predictions made in literature (e.g., [START_REF] Zanetti | Quaternion estimation and norm constrained Kalman filtering[END_REF][START_REF] Pittelkau | An analysis of the quaternion attitude determination filter[END_REF]).

This A-EKF estimation filter has been used in the sequel, and the tests in Sec.III brought numerically stable results.

The time-discrete measurement equation must express the relative chaser-target attitude quaternion ( = ℎ-) as a function of the state vector -. ℎ-also depends on the absolute chaser quaternion -ℎ , which is supposed to be previously estimated by a classic M-EKF. Exploiting the quaternion properties it is possible to write the quaternion product in matrix form:

ℎ-= * -ℎ ⊗ -= Σ( * -ℎ ) - , with Σ( ) =           0 -1 -2 -3 1 0 -3 2 2 3 0 -1 3 -2 1 0           (18)
The time-discrete output equation = + is therefore linear in the variable -, with:

= Σ( * -ℎ ) ∅ 4x3 , =        - -        (19)
Since -ℎ is the result of an estimation process, the covariance associated to ℎ-should take into account also the uncertainties introduced by Σ( * ℎ ) through the computation of the composed variance of the function ℎ-= ( -ℎ , -). Since this Note focuses on the characterization of the intrinsic performance of the delay management techniques, the true -ℎ will be used during the performance analysis not to introduce coupling between the covariances of chaser and target states.

Filter Recalculation method implementation

Since the system dynamics is a continuous process (Eq.( 17)), while the measurement is a discrete process (Eq.( 18)), the CD-EKF (Continuous Discrete EKF, [START_REF] Frogerais | Various ways to compute the continuous-discrete extended Kalman filter[END_REF][START_REF] Kulikov | Accurate numerical implementation of the continuous-discrete extended Kalman filter[END_REF]) structure will be exploited.

Φ = + + 2 2 , with          = ( ˆ -1| -1 ) ( ˆ | -1 ) = ( ˆ -1| -1 ) + ( ˆ | -1 ) /2 (21)
where ( ˆ -1| -1 ) and ( ˆ | -1 ) are the Jacobian of the state equation , evaluated respectively in the estimated state at -1 and the predicted state at . The prediction of the covariance matrix can be done using the discrete Lyaponouv equation in Eq.( 2), adapted for a non-linear system:

| -1 = Φ -1| -1 Φ + (22)
In fact, as discussed in [START_REF] Carpenter | Navigation filter best practices[END_REF], both methods in Eq.( 20) and Eq.( 22) can be used to compute the prediction of the state error covariance matrix for non-linear systems. The choice of using the discrete Lyaponouv equation is done in order to save computational resources. Indeed, the matrix Φ needs to be computed to implement Larsen's method, and is computed as follows [START_REF] Gibbs | Advanced Kalman filtering, least-squares and modeling: a practical handbook[END_REF] :

= ∫ -1 Φ( , ) ( ) ( ) ( )Φ ( , ) = ∫ -1 + ( -) + ( -) 2 2 ( ) ( ) ( ) + ( -) + ( -) 2 2 ( 23
)
where ( ) and ( ) have been considered constant along the interval (which is a good assumption since matrix depends on target inertia matrix and disturbance torques have a slow dynamics) and Φ( , ) has been substituted by the second order approximation of Eq.( 23). The expression of the resulting can be truncated at lower orders if needed, namely for small values of .

IV. Simulations and Performance analysis

The

A. Performance of the translational dynamics estimation

The simulated scenario is the following: the target is on a circular orbit at an altitude of 765 km and the chaser is approaching along the -R-Bar side. The initial relative target-chaser position in is: 0 = [-50, 0, 0 ] m. The chaser is subjected to a continuous profile of thrust in order to perform an R-bar maneuver and intercepts the target after 500 seconds at a relative speed of 10 cm/s.

Four MC scenarios have been selected in order to test the performance of the filtering techniques under uncertainties in the knowledge of the true covariance of the measurements and the knowledge of the chaser applied acceleration (i.e., which is the input of CWH state-space representation).

The latter may be due to an error in the knowledge of chaser mass or to a difference between the commanded and the true thrust. All the simulations have a relative position initialization error uniformly distributed in the interval [-10 , +10 ] for component and [-5 , +5 ] for and components. The relative velocity estimate is always initialized to be 0 m/s along each direction.

Table 1 summarizes the different conditions tested in each MC scenario. In cases . and .

the standard deviation of the generated measurements noise is fixed to a constant value of 2 m for the component and 1 m for and components. These values of are quite representative

T.A Δ ∈ [-10 , +10 ] Δ , Δ ∈ [-5 , +5 ] Δ = 0 = 2 , = 1 , = 1 Case T.B Δ ∈ [-10 , +10 ] Δ , Δ ∈ [-5 , +5 ] Δ ∈ [-0.25 , 0.25 ] = 2 , = 1 , = 1 Case T.C Δ ∈ [-10 , +10 ] Δ , Δ ∈ [-5 , +5 ] Δ = 0 = 0 + Δ , Δ ∈ [-0.8 0 , 0.8 0 ], 0 = [2, 1, 1] Case T.D Δ ∈ [-10 , +10 ] Δ , Δ ∈ [-5 , +5 ] Δ ∈ [-0.25 , 0.25 ] = 0 + Δ , Δ ∈ [-0.8 0 , 0.8 0 ] 0 = [2, 1, 1]
for a distance around 50 m, but are overestimated for shorter distances (i.e., the second half of the simulations), since camera sensors are characterized by an increase of measurement accuracy with decreasing range [START_REF] Fehse | Automated rendezvous and docking of spacecraft[END_REF]. The measurement component along direction is generated with a higher standard deviation with respect to the other components since during an R-bar maneuver the target axis corresponds to the optical axis of the camera. Cases . and . test the sensitivity of filter performance under the presence of random variations in the standard deviation of the measurement noise, in order to assess the robustness to a variable covariance . The standard deviation will be equal to 0 + Δ , where Δ is a uniformly distributed variable that varies at any instant in the interval [-0.8 0 , +0.8 0 ] and 0 corresponds to the values defined for cases . and . . The filter is unaware of the measurement noise variation, therefore both cases . and .

will have the nominal tuning respectively of . and . . In cases . and . , chaser acceleration is supposed to be known and only the sensitivity to uncertainties on the initials condition is tested.

In cases . and . , an uniformly distributed uncertainty in the interval of [-25%, +25%] is added to the knowledge of chaser control accelerations. In order to grant the convergence of the filter in these cases, the diagonal terms of matrix Q are increased of two order of magnitude with respect to case . . Tables 2 show the steady-state performance of the filter. The reported values,

], which could be the convergence interval of a classical model-based recursive tracking algorithm [START_REF] Lepetit | Monocular Model-Based 3D Tracking of Rigid Objects: A Survey[END_REF]. The uncertainty on the target inertia diagonal terms is uniformly distributed in the interval [-0.2 , +0.2 ] which is quite probable for RDV where the industrial model of the target is supposed to be known but still there could be uncertainties on the amount of remaining propellant or on the degradation of the S/C. Table 4 shows the steady-state performance of the filters for each one of the scenarios described, averaged on 200 MC runs over a simulation of 500s. The performance of the two delay management techniques is compared to the performance of a classic CD-EKF processing the infrequent measurements without delay. For the attitude estimation, the estimation error is given in the axis-angle representation, which provides a scalar representation of the error. The reported values in Table 4, averaged on 200 MC runs, are: (i.e., the root mean square of the generated measurement noise), (i.e., the root mean square of the AKE for the filter with infrequent non-delayed measurements), (i.e., the root mean square of the AKE for the FR method), and (i.e., the root mean square of the AKE for the Larsen's method). As for the translational dynamics, also the attenuation index for each filter is shown. For the rotation rate components, the root mean square corresponds to the standard deviation of the estimation since the mean of the estimate is zero. The results confirm the expectations: the performance of the attitude estimation is in all the cases better for the filter without delay, followed by FR method and then by Larsen's method. On the other side the performance of Larsen's method in the cases without IM ( . and . ), but is more performing in the case with IM ( . and . ), showing much lower error overshoots and a comparable convergence time. This degradation of Larsen's method performance can be explained as follows: in the case without IM Larsen's method was optimal (i.e., to the extent that it was linearizing the propagation of the update); in the case with IM Larsen's method becomes sub-optimal also with respect to Kalman theory, since the correction term * is computed using Kalman gains that do not take into account the contribution of the delayed measurement. The angle error of the no-IM cases shows the peculiarity of appearing piece-wise linear: this is due to the fact that the estimate evolves in open loop as long as a new measurement arrives and only every steps the state is corrected by the Kalman update.

C. Execution time and needed storage

Both Larsen's and FR methods have almost the same latency time for the time steps in which no delayed measurements arrive, while, at the arrival of a delayed measurement, FR method has to completely recompute the estimate through the delay period. Another important aspect is the amount of data that need to be stored for each method, with respect to a classical KF without delayed measurements. Table 7 shows the amount of "double" to be stored within a delay period (from to + ) for each method, under the following hypothesis: = 10 (i.e., the number of delay samples of the delayed measurements), = 10 Hz (i.e., the navigation filter run frequency), / = 1 Hz (i.e., the frequency of the delayed measurements). is the size of the IM vector . Larsen's method, due to the need of propagating through each time instant the matrix * , requires to allocate a higher memory with respect to a classical KF, but the required space does not depend on the number of the delay samples nor on the size of the IM vector. On the other side, the FR method, which for small values of delay and small size of IM vector requires to store a lower amount of data, rapidly increases its storage burden as or increase. Figure 2 shows the required number of double to be stored as a function of the number of delay samples , for different values of . In the case with no IM ( = 0), FR requires a higher storage than Larsen's methods for > 13 (rotational dynamics) and for > 12 (translational dynamics). Therefore Larsen's method needs to store a rather high amount of data, which is almost comparable to the one needed by the FR method for values of ∼ 10 (which is reasonable value taking into account latency time of IP-CV algorithms and the typical navigation filter run frequency). However, the main advantage of Larsen's method concerns the computational burden, which is equally distributed over every time step.

V. Conclusion

The problem of incorporating delayed and multi-rate measurements in a navigation filter for the estimation of the dynamics of a non-cooperative target has been assessed. A dynamic filter for the estimation of full target rotational and translational state exploiting relative pose measurement has been formalized. Two delay management techniques have been compared: Larsen's method, which provides a fast but sub-optimal solution, and Filter Recalculation method, which always provides the optimal estimate but has a higher computational load. The Monte Carlo validation campaign has shown that Larsen's method performance is comparable to Filter Recalculation method performance. The latter shows remarkably better performance only in the transient phase of simulations exploiting interim measurements but at the expense of a higher computational and storage need. When a delayed measurement arrives, Filter Recalculation method computational load is multiplied by a factor equal to the number of delay samples, which might be incompatible with the critical applications run by the on-board computer for this particular time step. This suggests that, in applications where the on-board resources are limited, Larsen's method is preferable since it provides a faster estimation without any significant degradation of the steady-state performance.

  Kalman Filters have been implemented in a full RDV simulator developed in Simulink and have been tested in a Monte Carlo (MC) Campaign. The generated true S/C dynamics takes into account all the principal sources of disturbance in LEO environment, as well as chaser thrusters acceleration and control torques. The filters are tested with simulated measurements that are generated by adding a Gaussian noise to the true relative state. It is necessary to go through this stage in order to test the performance of the filters under Kalman optimality hypothesis before coupling it with IP-CV measurements. Measurements resulting from IP-CV algorithms are affected by noise which depends on many factors, such as the intrinsic noise of the sensor, the relative distance camera-target, the relative rotation rate and velocity, the camera capture rate, the illumination conditions and even the target relative pose itself. This makes it very difficult to compute on-line in real time a representative model of the covariance matrix which is valid in any condition. The delay management techniques must therefore demonstrate robustness with respect to uncertainties in the knowledge of matrix . In all the simulated scenarios, the relative pose measurements are acquired at a rate of 1 Hz and become available for the filter after a delay of 1 second, which corresponds to = 10 assuming that the navigation filter operates at 10 Hz. These are all reasonable values taking into account the typical latency time of an IP-CV algorithm and navigation filter run frequency with typical space processing capabilities. All the presented scenarios are tested over 200 MC runs on a 500 seconds simulation. The state error covariance matrix is initialized as the identity matrix.

  D, FR Method Case R.D, Larsen's Method Case RI.D, FR Method Case RI.D, Larsen's Method (b) Cases . and .

Fig. 1

 1 Fig. 1 Transient phase estimation error in a single run

Fig. 2

 2 Fig. 2 Amount of double to be stored as a function of the delay samples

Table 1 Definition of the MC scenarios for the translational dynamics

 1

	Position initialization error	Control thrust knowledge error	Simulated measurement noise
	Case

Table 2 Performance of the translational dynamics estimation

 2 , the standard deviation of the estimation error -Absolute Knowledge Error AKE measurement noise variation, both . and . have the nominal tuning for a constant = 2 deg. Case . tests the performance of the methods under uncertainties in the initialization of the state of the filter. Target attitude quaternion estimation is initialized by adding a random error Δ ∈ [-40 , 40 ] expressed in Euler angles to the true state. Case . tests the filters performance under the presence of uncertainties in the knowledge of target inertia matrix and no error in the state initialization. The inertia matrix used in the filter is obtained by adding to each diagonal term of the true inertia matrix a value , where is a uniformly distributed variable in the interval [-0.5 , +0.5 ]. Cases . and . reproduce scenarios with uncertainties values closer to the ones encountered in a real vision-based RDV in space. The initialization error Euler angles in both cases is in the interval [-20 ,

			Case T.A			Case T.B			Case T.C			Case T.D
			(1 -/ )		(1 -/ )		(1 -/ )		(1 -/ )
	[m]	1.999	0.053	97.33%	1.999	0.145	92.77%	1.924	0.052	97.29 %	1.921	0.138	92.82%
	[m]	1.000	0.033	96.73%	1.000	0.099	90.22%	0.985	0.032	96.71%	0.967	0.095	90.22 %
	[m]	1.000	0.025	97.47%	1.000	0.070	93.03 %	0.965	0.024	97.46%	0.989	0.069	93.01 %
	[m/s]	-	0.232•10 -3	-	-	0.182•10 -2	-	-	0.225•10 -3	-	-	0.173•10 -2	-
	[m/s]	-	0.203•10 -3	-	-	0.173•10 -2	-	-	0.201•10 -3	-	-	0.168•10 -2	-
	[m/s]	-	0.103•10 -3	-	-	0.088•10 -2	-	-	0.099•10 -3	-	-	0.087•10 -2	-
	averaged over 200 MC runs, are:	(i.e., the standard deviation of the generated measurement
	noise),	(i.e.

Table 6 Performance of the rotational dynamics estimation with interim measurements

 6 , within a delay interval, the state does not evolve in open loop but it continues being corrected by the fast measurements. The steady-state performance of the FR method is only slightly better than Larsen's one, and it is only during the transient phase that FR method shows a remarkably better performance. The performance of the filters in the transient phase can be observed in Fig.1, which shows the first 200 seconds of a single-run in cases . , . (Fig.1a), and . , . (Fig.1b), with the same filter state initialization (i.e., an error on Euler attitude angles of

											Case RI.C	Case RI.D
								[	]		0.667 0.690	0.736 0.790
									[	/ ] 0.0077 0.0077	0.0137 0.0140
									[	/ ] 0.0101 0.0101	0.0200 0.0212
									[	/ ] 0.0113 0.0114	0.0281 0.0293
	without IM since[+20 , -10	, +15
		40						Case R.C, FR Method
		35						Case R.C, Larsen's Method
	error angle [deg]	15 20 25 30						Case RI.C, FR Method Case RI.C, Larsen's Method
		10
		5
		0 0	20	40	60	80	100	120	140	160	180	200
							time [s]

]) and the same uncertainty on the diagonal terms of target inertia matrix (i.e., 20%). In the figures it can be noticed how the filters exploiting IM converge faster than

Table 7 Amount of double to be stored within a delay period by both estimation filters

 7

		Rotational Dynamics	Translational Dynamics
		Recalculation	Larsen	Recalculation	Larsen
	ˆ | -1	7x1		7x1	6x1		6x1
	| -1	7x7		7x7	6x6		6x6
	-ℎ	4x		4x1 ( -ℎ )	-		-
		x		-	x		-
		-		-	3x		-
	*	-		7x7	-		6x6
		56+(4 +	)	109	42+(3+	)	78

Both state and covariance prediction equations can be written in the following time-continuous form:

where the second equation is a differential Lyaponouv equation derived using the approximated state equation ∼ ( , ) + ( ) ( ), with ( ) the input matrix of the process noise [START_REF] Crassidis | Optimal estimation of dynamic systems[END_REF]. The prediction step is computed by numerical integration from -1 to of the ordinary differential equations in Eq.( 20) using the explicit fourth order Runge-Kutta (RK) method. Numerical integration procedures may provide solution of that are not necessarily positive semi-definite matrix, which is in contradiction to the intrinsic properties of the state error covariance matrix. The numerical integration of the class of coupled differential equations in Eq.( 20) is investigated in [START_REF] Mazzoni | Computational aspects of continuous-discrete extended Kalman-filtering[END_REF], which suggests procedures that ensure stable solutions and guarantee positive semi-definite covariance matrices . Nevertheless, for the rotational dynamics estimation problem, it has not been experienced any issues either with the stability of the solution or with the properties of the covariance matrix.

With a filter run frequency of 10 Hz, a single sub-step for the fourth order RK integration is necessary.

The so computed ( ) and ( ) are the prediction of the state ( ˆ | -1 ) and of the covariance matrix ( | -1 ). From this moment on, the CD-EKF will follow the steps of a classical Discrete KF to compute the gain , the covariance update and the state update using the discrete measurement . When the delayed measurement * arrives, the update at time is computed using the stored ˆ | -1 , | -1 , and -ℎ (that is needed to compute matrix * ). No inputs need to be stored since the disturbance torques acting on the target affect the dynamics as process noises. Then the filter implementation follows the steps explained in Sec.II.A. If IM are present, -ℎ has to be stored for any time step going from to + .

Larsen's method implementation

Larsen's method has been expressly designed for linear system: the correction term * requires the knowledge of the transition matrix Φ = , implying that, for a non-linear system, an approximation must be computed. Here a second order RK approximation will be used [START_REF] Carpenter | Navigation filter best practices[END_REF]:

according to [START_REF] Ott | Esa pointing error engineering handbook[END_REF] standard-), and (1 -/ ) (i.e. an index in percentage of the noise attenuation introduced by the filter. As no IM are present, both FR method and Larsen's method provide the same results, so there is no distinction in the performance of the estimation. Comparing . and

. to the correspondent cased with known control acceleration . and . , it is possible to see a slight performance degradation. Indeed, the calibration of used in cases . and . provides very high level of attenuation (around 97%) that is paid by a loss of robustness with respect to uncertainties in the knowledge of the control accelerations, which could even make the filter diverge.

The attenuation index is slightly degraded from case . to case . , while it is almost the same for cases . and . , where the uncertainty on the process noise is dominant with respect to the uncertainty on the measurement noise.

B. Performance of the rotational dynamics estimation

The simulated rotational dynamics is the following: the chaser is rotating very slowly around the axis in order to ensure target pointing during the R-bar maneuver, while the target is rotating under the effect of its initial conditions and of the orbital disturbances. Four scenarios have been analyzed, case . , case . , case . and case . , whose different conditions are summarized in Table 3. The initial target rotation rate is equal to 1 deg/s around each body axis for cases . , . , and . , and equal to 3 deg/s around each body axis for case . . These last rotational rates are representative of drifting S/C rotation rates, and remain particularly challenging for an IP-CV algorithm running on a space processor. In all the four scenarios, the estimated rotation rate is initialized at 0 deg/s. In cases . and . the standard deviation of the measurements noise variable that varies at any instant in the interval [-0.8 0 , +0.8 0 ]. Since the filter is unaware of the the estimation of the angular rate is almost comparable for all the filters. In general the filter shows better performance under uncertainties in the initialization ( . ) with respect to uncertainties in the knowledge of the inertia matrix ( . ) -which actually corresponds to uncertainties in the knowledge of the prediction model. The performance degrades as the rotation rate increases (case . is the only one having an attenuation index lower than 70%). In any case, it is worth highlighting that the attenuation index of Larsen's method is only one percent below the one of FR method. Also in the transient phase the FR method shows very slightly better performance than the Larsen's method. The scenarios described by case . . and . . have been tested also adding a set of IM of ℎ-. This could be the case of measurements coming from marker-based methods, which have a relatively low latency time. The measurements are generated with a rate equal to the filter run frequency. They are affected by a Gaussian noise having standard deviation equal to 0 + Δ , with Δ uniformly distributed in the interval [-0.8 0 , +0.8 0 ] and 0 = 4deg. The characteristics of the simulated scenarios, named . . and . . are summarized in Table 5. Table 6 shows the performance of FR and Larsen's method. The root mean square of the averaged estimation errors of each case can be directly compared to the values of the corresponding case in Table 4, taking into account that the only difference is the presence of high frequency attitude measurements with higher standard deviation. The performance of both methods is highly increased with respect to the cases

Table 4 Performance of the rotational dynamics estimation