Robert Fonod
email: robert.fonod@ims-bordeaux.fr

David Henry

Catherine Charbonnel

Eric Bornschlegl

Ph D Candidate

Position and Attitude Model-Based Thruster Fault Diagnosis: A Comparison Study

This paper deals with performance and reliability evaluation of a fault diagnosis scheme based on two distinct models to detect and isolate a single thruster fault affecting a chasing spacecraft during rendezvous with a passive target in a circular orbit. The analysis is conducted in the frame of a terminal rendezvous sequence of the Mars Sample Return mission. A complete description of a robust residual generation design approach based on eigenstructure assignment is presented.Unknowntime-varying delays, induced by the thruster drive electronics and uncertainties on thruster rise times, are considered as unknown inputs. Particular novelty of the work is a new method for estimating the unknown input directions used to enhance the robustness properties of the diagnosis scheme. Monte Carlo results from a high-fidelity industrial simulator and carefully selected performance and reliability indices allows us to evaluate the effectiveness of both schemes. The obtained results reveal that the proposed fault diagnosis scheme based on a position model is a justified competitor to the conventionally used attitude model-based scheme.

I. Introduction

Many space exploration missions require critical autonomous proximity operations. Mission safety is usually guaranteed through a hierarchical implementation of fault/failure detection, isolation and recovery (FDIR) approach with several levels of fault containments defined from local component/equipment up to global system, i.e., through various equipments (sensors, thrusters, reaction wheels etc.) redundancy paths and ground intervention [START_REF] Osder | View of Redundancy Management Application and Theory[END_REF][START_REF] Olive | FDI (R) for satellites: How to deal with high availability and robustness in the space domain?[END_REF][START_REF] Zolghadri | Advanced Model-Based FDIR Techniques for Aerospace Systems: Today Challenges and Opportunities[END_REF]. In the case of the Mars Sample Return (MSR) mission, the hierarchical implementation of FDIR is concerned at three levels [START_REF] Peuvédic | Autonomous Rendezvous Control System: a High Fidelity Functional Engineering Simulator Developed for GNC/AMM/FDIR Validation[END_REF][START_REF] Strippoli | Integrated vision-based GNC for Autonomous Rendezvous and Capture around Mars[END_REF]: (i) based exclusively on sensor measurements in the fault detection and isolation (FDI) function which are mainly signal-based techniques, (ii) relying on both actuator commands and sensor measurements in the model-based FDI function, and (iii) based on navigation outputs for the monitoring of the trajectory in the safety monitoring function. Through this hierarchical strategy, it is assumed that sensor faults are diagnosed and recovered at the first level. This paper focuses on the second level of this hierarchy. The goal is to develop a robust model-based FDI unit for faults occurring in the spacecraft involved in the MSR mission.

Common FDIR techniques may be not sufficient in some cases, especially for faulty situations that cause quickly abnormal dynamics deviation in critical space operations. This is especially the case during rendezvous and docking/capture proximity operations. Advanced model-based FDI techniques should be particularly developed to safely accommodate on-board (and on-line) the necessary robustness/stability of the spacecraft control, the necessary trajectory dynamics and the vehicle nominal operation.

To ensure normal operation, real-time fault diagnosis is essential to provide information for the spacecraft to accommodate the fault in time. Numerous model-based FDI techniques are reported in the literature and have potential to be applied [START_REF] Isermann | Model-Based Fault Detection and Diagnosis -Status and Applications[END_REF][START_REF] Venkatasubramanian | A Review of Process Fault Detection and Diagnosis: Part I: Quantitative Model-Based Methods[END_REF][START_REF] Henry | Fault Detection and Diagnosis for Aeronautic and Aerospace Missions[END_REF][START_REF] Zolghadri | Fault Diagnosis and Fault-Tolerant Control and Guidance for Aerospace Vehicles[END_REF]. The problem of thruster faults is however less studied. Among other examples, methods based on the so-called unknown input observer (UIO) technique were applied to the Mars Express mission [START_REF] Patton | Reliable Fault Diagnosis Scheme for a Spacecraft Attitude Control System[END_REF][START_REF] Patton | Robust FDI Applied to Thruster Faults of a Satellite System[END_REF]; an iterative learning observer (ILO) was designed to achieve estimation of time-varying thruster faults [START_REF] Chen | Observer-Based Fault Diagnosis of Satellite Systems Subject to Time-Varying Thruster Faults[END_REF]; H ∞ /H -filters were used to addresses the problem of thruster fault diagnosis in the Microscope satellite [START_REF] Henry | Fault Diagnosis of Microscope Satellite Thrusters Using H∞/H-Filters[END_REF] and also the problem of faults affecting the micro-Newton colloidal thrust system of the LISA Pathfinder experiment [START_REF] Falcoz | Robust H∞/H-Thruster Failure Detection and Isolation with Application to the LISA Pathfinder Spacecraft[END_REF]. Both H ∞ /H -filters are based on a residual generator, that is robust against spatial disturbances (e.g., J 2 disturbances, atmospheric drag and solar radiation), measurement noises and sensor misalignment phenomena, whilst guaranteeing some fault sensitivity performance.

On the other hand, only limited results on FDI of time-delay systems have been developed in recent years in the literature. Among the contributions, an UIO was designed for fault detection of state-delayed systems with known delays [START_REF] Yang | Observer Design and Fault Diagnosis for State-Retarded Dynamical Systems[END_REF]. The well known parity space approach was extended for fault detection of retarded time-delay systems [START_REF] Kratz | Fault Detection for Time-Delay Systems: A Parity Space Approach[END_REF]. Two-objective optimization approach was considered for linear time-invariant (LTI) systems with constant time delays aiming at formulation of the optimization problem as: enhancing sensitivity of the residual to faults and at the same time suppressing the undesirable effects of unknown inputs and uncertainties in L 2 -gain sense [START_REF] Ding | An LMI Approach to the Design of Fault Detection Filter for Time-Delay LTI Systems with Unknown Inputs[END_REF][START_REF] Jiang | H∞ Fault Detection Filter Design for Linear Discrete-Time Systems with Multiple Time Delays[END_REF].

Robust fault diagnosis approach based on an adaptive observer was developed for uncertain continuous LTI systems with multiple discrete time-delays in both states and outputs [START_REF] Jiang | Fault Detection and Identification for Uncertain Linear Time-Delay Systems[END_REF]. Recently, a geometric approach for FDI of retarded and neutral time-delay systems was developed [START_REF] Meskin | Robust Fault Detection and Isolation of Time-Delay Systems Using a Geometric Approach[END_REF]. Problem of robust fault detector design for a class of LTI systems with some nonlinear perturbations and mixed neutral and discrete time-varying delays was investigated using a descriptor technique, Lyapunov-Krasovskii functional and a suitable change of variables [START_REF] Karimi | A Linear Matrix Inequality Approach to Robust Fault Detection Filter Design of Linear Systems with Mixed Time-Varying Delays and Nonlinear Perturbations[END_REF].

Often, uncertain time-varying delays occurs in the input channel of the system. Considering spacecraft's thrusters, such delays may be induced by the propulsion drive electronics. One of the main difficulty in fault diagnosis of systems with uncertain input delays lies in the fact that uncertainty caused by this delay is unstructured. Therefore, robustness cannot be achieved by applying existing robust unknown input (disturbance) decoupling approaches directly. There is an

important assumption for all such approaches that the disturbance distribution matrix must be known, but a generalized approach to obtain the matrix is still lacking [START_REF] Patton | Robust Fault Detection of Jet Engine Sensor Systems Using Eigenstructure Assignmen[END_REF][START_REF] Patton | Parameter-Insensitive Technique for Aircraft Sensor Fault Analysis[END_REF].

One possible solution to this problem is to transform the unstructured input uncertainty to unknown input by means of time-dependent Padé approximation [START_REF] Fonod | Robust Thruster Fault Diagnosis: Application to the Rendezvous Phase of the Mars Sample Return Mission[END_REF]. In this paper, we have taken an approach similar to [START_REF] Fonod | Robust Fault Detection for Systems with Electronic Induced Delays: Application to the Rendezvous Phase of the MSR Mission[END_REF], where by introducing a Cayley-Hamilton theorem-based and horder Taylor series expansion-based polytopic transformations, the influence/uncertainty of the timevarying delay on the system is summarized as an unknown input. By assigning the eigenstructure of the observer, the residuals can be decoupled from these unknown inputs [START_REF] Patton | Robust Fault Detection of Jet Engine Sensor Systems Using Eigenstructure Assignmen[END_REF]. The solutions investigated in this paper follow this strategy. Two FDI schemes are proposed based on two different models of the spacecraft dynamics, i.e., the so-called position model scheduled by judiciously chosen parameters, which exhibit the translation dynamic of the spacecraft in a judiciously chosen frame, and the attitude model.

The position model expressed in the local reference frame is well known and mastered for control, but rarely used for FDI of thrusters. The reason is quite evident, the attitude model seems to be more sensitive to thruster faults. The innovation that we pursue with this study is concerning the judiciously chosen position model for FDI and its comparison to the attitude model-based solution.

Simulation results from "high-fidelity" MSR industrial simulator, under realistic conditions and taking into account the effects that the GNC (Guidance, Navigation, Control) unit has on the FDI performances, demonstrate the efficiency and capabilities of the proposed methodology. These two FDI schemes have been successfully demonstrated as applicable for FDI of the chaser spacecraft involved in the MSR mission within the critical dynamics and operation constraints of the last terminal translation of the rendezvous/capture phase. The research work presented in this paper draws expertise from actions undertaken within the European Space Agency (ESA), Thales Alenia Space and the IMS laboratory (University of Bordeaux). The aim is to demonstrate the benefits of novel robust on-board FDI technologies, that may significantly enhance the spacecraft autonomy.

The paper is organized as follows: Section II briefly describes the MSR mission and introduces the modeling of the chaser spacecraft propulsion and control system. Necessities of modeling the spacecraft dynamics during the rendezvous phase are recalled in Sec. III. Section IV addresses the theoretical basements and computational procedure of the FDI scheme design. For a detailed assessment of the FDI performances, results of the Monte Carlo simulation campaign incorporating different fault scenarios are presented in Sec. V. Finally, some concluding remarks are given in Sec.

VI.

II. Background

A. Overview of the MSR mission

The MSR mission is one of the most exciting challenges in the international effort on the Solar System exploration. The mission concepts have been studied for years by NASA [START_REF] D'amario | Mars Orbit Rendezvous Strategy for the Mars 2003/2005 Sample Return Mission[END_REF], French National Space Agency (CNES) [START_REF] Cazaux | The NASA/CNES Mars Sample Return -A Status Report[END_REF], and ESA. Its main goal is to collect samples of Martian rocks, soils and atmosphere, and to return these samples safe and intact back to Earth for analysis. Finally, after successful capture, the sample container is inserted into an Earth re-entry capsule (ERC) inside the chaser vehicle and the chaser starts its interplanetary cruise towards the Earth. It should be noted, that this autonomous capture maneuver belongs to the non-cooperative rendezvous since the target is not equipped with any actuation system. Figure 1 provides an overview of the mission. The work addressed in this paper concerns the last 20m of the rendezvous phase.

B. Chaser Spacecraft Configuration

Figure 2 shows the general setup of the GNC system in the chaser vehicle (for confidential reasons, the numerical values with regards to the spacecraft geometry and characteristics are omitted). The selected attitude and orbital control system dedicated to spacecraft mode for the rendezvous phase corresponds to 6 degree of freedom (DOF) control. The 6 DOF control ensures the application of both commanded torque and force using thrusters only.

The control system relies on the following sensors:

• 2 3-axis inertial measurement unit (IMU) in hot redundancy;

• 2 star trackers (STR) in cold redundancy;

• 2 short-range rendezvous sensors with a functional hot redundancy, i.e., a light detection and ranging (LIDAR) sensor and a radio frequency sensor (RFS) as back-up.

The IMU unit is in charge of measuring the angular velocities ω = [p, q, r] T . The STR device gives the attitude measurement Θ = [ϕ, θ, ψ] T . The LIDAR unit is in charge of the measurement of the relative position ρ = [ξ, η, ζ] T between the chaser and target. The RFS sensor is used to monitor the chaser trajectory, and it can trigger a collision avoidance maneuver, if necessary.

The role of the navigation (NAV) unit is to perform an estimate Θ, ω, ρ of Θ, ω, ρ, respectively, by removing the misalignment phenomena, sensor bias and noises on these measurements. The NAV unit also provides an estimate of the target attitude quaternion Qtgt , that will be used later for the design of the FDI unit. However, we assume that the NAV unit is not perfect and, thus, that there still exists time delays and noises on Θ, ω, ρ, and Qtgt , respectively.

In terms of actuators, the chaser vehicle is equipped with a very precise chemical propulsion system composed of 16 thrusters (THR) in full redundancy. They are organized in two sets, the nominal set 'A' is used for the nominal vehicle control and the redundant set 'B' is reserved for the recovery actions after a fault has been detected in the active set. The thruster configuration is illustrated in Fig. 3.

Z X Y 1A 1B 2A 2B 4A 4B 5A 5B 6A 6B 7A 7B 8A 8B

C. Modeling the Control Loop and the Propulsion System

To carry out its mission and to ensure the mission performance, the chaser is equipped with 6 DOF control law, see Fig. 4 for an illustration. It consists of two dynamic controllers: K pos (s) that is in charge of controlling the position dynamics and K att (s) that aims to regulate the chaser's attitude. The controller outputs, the desired force F d and the desired torque C d , are sent to the thruster modulator unit (TMU) that integrates the small commanded pulses which are below the minimum impulse bit (MIB) and releases a pulse ( F d or/and C d ) when the total reaches a momentum threshold. The TMU is widely used to compensate the actuator nonlinearities and improve the control accuracy [START_REF] Peuvédic | Autonomous Rendezvous Control System: a High Fidelity Functional Engineering Simulator Developed for GNC/AMM/FDIR Validation[END_REF][START_REF] Peuvédic | Key Control Techniques for GNC Design of Martian Vehicles[END_REF]. The purpose of on-board thruster management function (TMF) is to select specific thrusters and compute their firing times T oni to realize force F d and torque C d command impulses coming from the TMU. The TMF algorithm relies on a simplified approach with respect to the Simplex and thrusters' non-linearities (minimum On/Off times) [START_REF] Peuvédic | Autonomous Rendezvous Control System: a High Fidelity Functional Engineering Simulator Developed for GNC/AMM/FDIR Validation[END_REF]. As mentioned earlier, the propulsion system uses only 8 of the 16 thrusters in the nominal case thus, in the following, we are focusing on these 8 thrusters. The thrusters have fixed directions and each one is able to produce a maximum thrust of 22N. The chemical propulsion drive electronics (CPDE), that drives the thrusting actuators, is initiating the opening of the thruster valve for the commanded duration 0 ≤ T oni ≤ 1, i = 1, . . . , 8. T oni are scaled on-times. The scaling is done versus the sampling period T of the control unit and is defined according to T oni = Toni /T , Toni being the actual/real firing time.

The propulsion system is a source of uncertainty in the system. The linear parameter-varying transfer function

H(s) = e -τ (t)s (1)
aims to model the effect of the unknown time-varying delays τ (t) ∈ R induced by the CPDE electronic device and the uncertainties on the thruster rise times [START_REF] Pettazzi | Design of Robust Drag-Free Controllers with Given Structure[END_REF]. In the following, it is assumed

that each thruster open duration T oni , i = 1, . . . , 8 is delayed with the same delay τ (t). This is a reasonable assumption from the practical point of view, since for all nominal thrusters the same CPDE device is used to control the opening of the thruster valve.

Let T oni t -τ (t) be the commanded open duration of the ith thruster delayed by τ (t), then the force vector

F thr = [F ξ , F η , F ζ ] T and the torque vector C thr = [C ϕ , C θ , C ψ ]
T , generated by the thrusters, are given by

    F thr (t) C thr (t)     = M T on t -τ (t) (2)
where T on = [T on1 , . . . , T on8 ] T and M ∈ R 6×8 is the thruster configuration matrix. The elements of M are the influence coefficients defining how each thruster affects each component of F thr and C thr , respectively.

D. Fault Considerations

It is obvious, that if a thruster fault occurs, for instance a hardover-type failure (thruster stuck open), it could lead to a drastic increase of the propellant consumption which is already very constrained by the travel to Mars. Dramatic consequences can occur, e.g., already in-placed GNC may not compensate such faults, possibly leading the chaser to lose the attitude and/or the position of the sample container. The problem becomes highly critical during the last 20 meters of the rendezvous phase when the chaser shall be correctly positioned in the rendezvous corridor (see Fig. 5 for illustration) in order to successfully capture the sample container as well as the chaser's attitude need to be maintained in the rendezvous sensor field of view.

Such faulty situations obviously cannot be diagnosed by ground support using telemetry information, due to the potential lack of communication between the chaser and the ground stations or due to significant communication delay. This motivates ESA to manage studies for the development of on-board fully autonomous FDI solutions that shall cope with all the failures which may occur and endanger the mission. A quick detection and isolation of the fault is the first step towards an efficient recovery action that has no impact on the mission success.

Through monitoring the sensor outputs, most of the faults in the chaser spacecraft can be detected at sensor data processing level (health status, operational mode, continuity checks, etc.) and recovery from these faults is usually quick enough, i.e., switching to a redundant sensor in hot redundancy (IMU case), or propagation of last valid measurement during the switch-on of the redundant equipment in cold redundancy (STR case). Most space agencies (NASA, ESA, CNES)

have already identified LIDAR as a preferred candidate instrument for autonomous rendezvous. It seems to be a justified and robust device for the mission success [START_REF] Pelletier | Lidar-Based Rendezvous Navigation for MSR[END_REF]. Moreover, the problem of fault diagnosis of multiple sensor faults occurring in the chaser spacecraft has been already addressed in the literature [START_REF] Issury | A Boolean Algebraic-Based Solution for Multiple Fault Diagnosis: Application to a Spatial Mission[END_REF]. Thus, the focus of this study concerns only thruster fault diagnosis.

Obviously, the mission can be endangered if a thruster fault occurs since it may have critical impact on the GNC system, i.e., it may lead to a degraded GNC performances or even the GNC become unstable. More precisely, we consider the following thruster fault scenarios:

• Case 1: stuck open valve -provides maximum thrust regardless of the demand;

• Case 2: thruster closing itself (blocked-closed) -thruster does not generate any thrust regardless of the demanded command by the TMF;

• Case 3: bi-propellant leakage -residual propellant leakage of size m leak (t), starting from 0 and reaching the maximum leakage size mleak > 0 with a given slope m s > 0, i.e., m leak (t) = min{m s (t -t f ), mleak }, where t f denotes the time of fault occurrence;

• Case 4: thrust loss -loss of efficiency of a particular thruster by a value mloss > 0;

Assuming no simultaneous faults, the considered thruster fault scenarios can be modeled in a multiplicative way according to (index "f is used to outline the faulty case)

    F f thr (t) C f thr (t)     = M (I 8 -Ψ(t)) T on (t -τ (t)) (3)
with Ψ(t) = diag ψ 1 (t), . . . , ψ 8 (t) , where 0 ≤ ψ i (t) ≤ 1, i = 1, ..., 8 are unknown. The health status of the ith thruster is modeled by ψ i (t) as follows

ψ i (t) =        0 if healthy 1 -ϕ i (t)/T oni (t) if faulty (4)
where ϕ i (t) allows to consider all four fault cases, mentioned earlier, as follows:

ϕ i (t) =        max{T oni (t), m leak (t)} if case 1 or 3 (1 -mloss )T oni (t) if case 2 or 4 (5)
where 0 < m leak (t) < 1 is the leakage size and 0 < mloss < 1 is the efficiency loss size. It is obvious that m leak (t) = 1, ∀t models a fully open (stuck open valve) and mloss = 1 a blocked-closed thruster fault, respectively.

III. Modeling the Chaser Dynamics During the Rendezvous Phase

In this section, a linear position model and an attitude model of the chaser spacecraft dynamics are introduced for FDI purpose. The two proposed models are able to describe the dynamics of the chaser in both, fault-free and faulty situations. For the sake of brevity, we recall only the necessary developments about modeling the spacecraft's dynamics, available in the extensive space literature [START_REF] Wie | Space vehicle dynamics and control[END_REF][START_REF] Schaub | Analytical Mechanics of Space Systems[END_REF].

A. Position Model

The translation motion of the chaser is derived from the 2nd Newton law. To proceed, let a, m, G and m M denote the radius of the circular orbit of the target, the mass of the chaser during rendezvous, the universal gravitational constant and the mass of Mars. Then, the orbit of the rendezvous being circular, the velocity of the target is given by the relation

µ a ( 6
)
target O T , with its Z l axis be perpendicular to the X l and Y l axis and oriented as shown in Fig. 6.

Fig. 6 The Mars rendezvous orbit with the associated frames.

The linear velocity of the target is given by the relation a.n, where n = ν stands for the uniform angular speed of the target. This velocity is given in the inertial frame

R i = {O M ; X i , Y i , Z i },
which is attached to the center of Mars O M . From the Kepler's third law it follows:

a.n = µ a ⇒ n = µ a 3 (7)
During the rendezvous phase, it is assumed that the chaser motion is due to the four following forces, all given in R l :

• the Mars attraction force • the centripetal force

F a = -m µ ((a+ξ) 2 +η 2 +ζ 2 ) 3/2 (a + ξ) X l + η Y l + ζ Z l ,
F e = m n 2 (a + ξ) X l + n 2 η Y l + 0 Z l ; • the Coriolis force F c = m 2n η X l -2n ξ Y l + 0 Z l ;
• the non-gravitational (chemical thrust, perturbations) force

F d = F dξ X l + F dη Y l + F dζ Z l .
Then, from the 2nd Newton law, it follows

                   ξ = n 2 (a + ξ) + 2n η - µ (a+ξ) 2 +η 2 +ζ 2 3/2 (a + ξ) + F dξ m η = n 2 η -2n ξ - µ (a+ξ) 2 +η 2 +ζ 2 3/2 η + F dη m ζ = - µ (a+ξ) 2 +η 2 +ζ 2 3/2 ζ + F dζ m (8)
Because the distance between the target and the chaser, during the rendezvous, is much smaller than the orbit, i.e., ρ a, it is possible to derive the so called Hill-Clohessy-Wiltshire (HCW) equations from Eq. ( 8) by means of a first order approximation [START_REF] Clohessy | Terminal Guidance System for Satellite Rendezvous[END_REF]. Thus, the motion of the chaser can be modeled in the target (local) frame R l , in both fault free (i.e., Ψ(t) = 0) and faulty (i.e., Ψ(t) = 0) situations, according to a linear 6th order state space model with state vector

x = [ξ η ζ ξ η ζ] T , i.e.
, from Eq. ( 8) it follows

       ẋ(t) = A p x(t) + B p R Qtgt (t), Qchs (t) F f thr (t) + E pw w(t) y(t) = C p x(t) (9)
where Qtgt ∈ R Considering Eq. ( 2) and Eq. ( 3), a new input vector u ∈ R 3 is defined according to

u(t) = R Qtgt (t), Qchs (t) F thr (t) (10)
and the fault model is approximated in terms of an additive fault vector f ∈ R 3 as follows:

f (t) = -R Qtgt (t), Qchs (t) M F Ψ(t)T on t -τ (t) (11)
where M F ∈ R 3×8 is the upper block of the thruster configuration matrix

M T = [M T F M T C ],
related to the forces F , see Fig. 4 if necessary. This type of approximation is widely used in the literature [START_REF] Isermann | Model-Based Fault Detection and Diagnosis -Status and Applications[END_REF][START_REF] Frank | Modelling for Fault Detection and Isolation versus Modelling for Control[END_REF].

Finally, the overall model of the position dynamics that takes into account both, the attitude Q chs , and the relative position ρ, can be derived from Eq. ( 9) using Eqs. ( 10) and ( 11) as follows [START_REF] Patton | Robust FDI Applied to Thruster Faults of a Satellite System[END_REF] where E pf = B p . This model is now suitable for the FDI filter design proposed in the next section.

       ẋ(t) = A p x(t) + B p u(t) + E pf f (t) + E pw w(t) y(t) = C p x(t)

B. Attitude Dynamics

The attitude control system works in a target pointing mode, which means that the chaser keeps one face of the spacecraft pointed to the target and maintains it during the whole rendezvous phase

ω(t) = I -1 chs C thr (t) + C dis (t) -ω(t) × I chs ω(t) ( 13
)
where × denotes the cross product of vectors and I chs ∈ R 3×3 is the inertia matrix. C thr and C dis is the external torque about the center of mass due to the propulsion and disturbances, respectively. ω = [p, q, r] T is the chaser rotational velocity vector, expressed in body frame R b .

Using the individual rotation matrices from Euler(3,2,1) rotation we can express the relation between the rotational velocities ω and the rate of the Euler angles Θ = [ϕ, θ, ψ] T . The inverse relationship becomes [START_REF] Wie | Space vehicle dynamics and control[END_REF]

Θ(t) = 1 cos(θ)         cos(θ) sin(ϕ) sin(θ) cos(ϕ) sin(θ) 0 cos(ϕ) cos(θ) -sin(ϕ) cos(θ) 0 sin(ϕ) cos(ϕ)         ω(t) (14)
Finally, noting that the chaser is controlled around the equilibrium point Θ = 0 and ω = 0, one can derive from Eqs. [START_REF] Henry | Fault Diagnosis of Microscope Satellite Thrusters Using H∞/H-Filters[END_REF][START_REF] Falcoz | Robust H∞/H-Thruster Failure Detection and Isolation with Application to the LISA Pathfinder Spacecraft[END_REF] a linear model by means of a first-order approximation of the nonlinear equations around the equilibrium point. Let the system's state vector be x = [ϕ, θ, ψ, p, q, r] T and the measurement vector y = [ϕ, θ, ψ] T , it boils down to a 6 order linear, pure attitude, model having the same structure as the position model given in Eq. ( 12), i.e.,

       ẋ(t) = A a x(t) + B a u(t) + E af f (t) + E aw w(t) y(t) = C a x(t) (15)
where w ∈ R 3 denotes the spatial disturbances. The input u and the fault f vector are given by

u(t) = C thr (t) (16
)

f (t) = -M C Ψ(t)T on t -τ (t) (17)
with M C ∈ R 3×8 being the lower block of the thruster configuration matrix M . The rest of the parameters A a , B a , E af and E aw are the result of the linearization procedure.

IV. Design of the FDI Scheme

This section is dedicated to the FDI task of thruster faults despite the presence of the delay τ (t) in the actuation system. Two fault diagnosis units based on the eigenstructure assignment (EA) approach of an observer-based fault detector are designed. The first is based on the position model ( 12) and the second one uses the attitude model [START_REF] Yang | Observer Design and Fault Diagnosis for State-Retarded Dynamical Systems[END_REF]. The isolation strategy is based on the evaluation of a cross-correlation like criterion between the component of the residuals and the thruster opening times.

A. Robust Residual Generator Design

Consider the following general description of the continuous LTI system

       ẋ(t) = Ax(t) + Bu(t) + E f f (t) y(t) = Cx(t) (18)
where x ∈ R nx , u ∈ R nu , y ∈ R ny , and f ∈ R n f is state, input, measurement, and unknown fault vector, respectively. The quadruplet {A, B, C, E f } represents the state-space matrices either of the position Eq. ( 12), or the attitude Eq. ( 15) model. It is assumed that all considered faults f are detectable (see [START_REF] Massoumnia | A Geometric Approach to the Synthesis of Failure Detection Filters[END_REF] for more details on fault delectability) and that the pair (A, C) is observable.

The spatial disturbances w are omitted from Eq. ( 18), because they have the same directional properties as those of faults, i.e., E f = E w , and thus, exact spatial disturbance decoupling cannot be achieved.

The TMF generates the thrusters' opening times T on equidistantly with a fixed sampling interval T > 0. The control signal T on (k), generated at time t = kT, k ∈ Z + , arrives at the actuator at time instant t = kT + τ (k). Assuming that the time-varying delay τ (k) is unknown but upper bounded, i.e., τ (k) ≤ τ , ∀k ∈ Z + , then the system's input, affected by delays, is given by

u(t) =          u c (k -1), t ∈ kT, kT + τ (k) u c (k), t ∈ kT + τ (k), (k + 1)T ( 19
)
where u c depends on the selected model according to

u c (k) =          R Qtgt (k), Qchs (k) M F T on (k), if position model is used M C T on (k), if attitude model is used ( 20
)
Our objective is to design discrete-time residual generator r(z) = H y (z)y(z) + H u (z)u(z), so that r is robust against the uncertain delay τ (k), H y & H u being observer-based filter transfer functions. To achieve this goal, the influence of the uncertainty τ (k) is first transformed to unknown input acting on the system model and then decoupled by means of eigenstructure assignment technique.

To proceed, assume τ (k) can be expressed as: τ (k) = lT + δ(k) ≤ τ , where l is a known integer, and δ(k) ∈ R is the unknown varying part of τ (k), bounded by 0 ≤ δ(k) < mT , with m being a known integer. Following, assume m = 1, which means that the variation part of the delay is less than one sampling interval. The case when m > 1 is discussed in [START_REF] Wang | A New Fault Detection Scheme for Networked Control Systems Subject to Uncertain Time-Varying Delay[END_REF].

If we assume that f is constant during each sampling interval T , what is a reasonable assumption from a practical point of view (the counterpart of this assumption is that the sampling interval T is chosen adequately so that Eq. ( 21) holds), then the discrete representation of Eqs. ( 18) and ( 19) is

                 x(k + 1) = Āx(k) + Γ 0 (δ(k))u c (k -l) + Γ 1 (δ(k))u c (k -l -1) + Ēf f (k) y(k) = Cx(k) (21)
where

Ā = e AT , Γ 0 (δ(k)) = T -δ(k) 0 e At dtB, Ēf = T 0 e At dtE f C = C, Γ 1 (δ(k)) = T T -δ(k)
e At dtB It's obvious that the following holds

B = Γ 0 (δ(k)) + Γ 1 (δ(k)) = T 0 e At dtB (22)
Using ( 21) and ( 22), and introducing a new augmented state vector z

T (k) = [ x T (k) u T c (k -l -1) ], we obtain                z(k + 1) = Â0 + Â(δ(k)) z(k) + B0 + B(δ(k)) u c (k -l) + Êf f (k) y(k) = Ĉz(k) (23)
where

Â0 =     Ā 0 0 0     , B0 =     B I     , Â(δ(k)) =     0 Γ 1 (δ(k)) 0 0     , B(δ(k)) =     -Γ 1 (δ(k)) 0     , Ĉ = C 0 Êf =     Ēf 0    
The system given in Eq. ( 23) is a time-varying system with unstructured uncertainty, where Γ 1 (δ(k))

is strongly dependent on the uncertain term δ(k). The remaining task is to transform this model to an uncertain polytopic system for which structured properties can be extracted in terms of unknown inputs. The polytopic system is then rewritten as a LTI system subject to an unknown input with a suitable distribution matrix.

Two polytopic transformations are used for this purpose. The first uses a Cayley-Hamilton theorem based transformation [START_REF] Wang | A New Fault Detection Scheme for Networked Control Systems Subject to Uncertain Time-Varying Delay[END_REF] and is introduced in Proposition 1.

Proposition 1. The Cayley-Hamilton theorem based transformation of Γ 1 (δ(k)) can be expressed as the convex matrix polytope

Γ a 1 (δ(k)) = 2nx i=1 µ a i (k)U a i ( 24
)
where µ a i (k) > 0, i = 1, ..., 2n x , ∀k ∈ Z + are uncertain scale factors satisfying

2nx i=1 µ a i (k) = 1, ∀k ∈ Z + ,
and U a i , i = 1, . . . , 2n x are known constant matrices given in Appendix A.

The second transformation is based on the h-order Taylor series expansion [START_REF] Hetel | Stabilization of Arbitrary Switched Linear Systems with Unknown Time-Varying Delays[END_REF] and is given in Proposition 2.

Proposition 2. The h-order Taylor series approximation of Γ 1 (δ(k)) can be expressed as the convex matrix polytope

Γ b 1 (δ(k)) = h+1 i=1 µ b i (k)U b i ( 25
)
where Taking into account the structure of the uncertain matrices Â(δ) and B(δ) in Eq. ( 23) and the two transformations of Γ 1 (δ) introduced in Eq. ( 24) and Eq. ( 25), the influence of the uncertain scalar factors µ a i and µ b i on the state x can be approximated in terms of unknown inputs as

µ b i (k) > 0, i =
2nx i=1 µ a i (k)U a i u c (k -l -1) -u c (k -l) = E a d d a (k) ( 26
) h+1 i=1 µ b i (k)U b i u c (k -l -1) -u c (k -l) = E b d d b (k) (27)
where

d a (k) = µ a 1 (k) u T c (k -l -1) -u T c (k -l) , ..., µ a 2nx (k) u T c (k -l -1) -u T c (k -l) T d b (k) = µ b 1 (k) u T c (k -l -1) -u T c (k -l) , ..., µ b h+1 (k) u T c (k -l -1) -u T c (k -l) T E a d = U a 1 , ..., U a 2nx , E b d = U b 1 , ..., U b h+1
Now, the mean value of the two unknown inputs, d a and d b , is considered. The augmented distribution matrix Êd and the augmented unknown input d take the following forms

Êd =     Ēd 0     =     E a d E b d 0     , d(k) = 1 2     d a (k) d b (k)     (28)
The elements (columns) Êd define the directions how each component of d affects the augmented state. This kind of approach gains advantage of combining two techniques to model the effect of the complex uncertainty δ on the state.

Finally, the augmented model with lumped unknown inputs can be expressed as

       z(k + 1) = Â0 z k + B0 u c (k -l) + Êf f (k) + Êd d(k) y(k) = Ĉz(k) (29)
This model is a quasi-equivalent representation of the augmented system given in Eq. ( 23), in other words, using two polytopic transformations, the influence of the uncertainty Γ 1 (δ) on the augmented state z is approximated in terms of unknown input d.

By closer examining the structure of Â0 , B0 , Êf , Êd in Eq. ( 29), one can see that only the upper state of z, i.e., the system state x, is influenced by f and d and, that there is no coupling between the lower and upper state. This allows us to consider only the upper state x in Eq. ( 29) for residual generator design. The following observer-based residual generator with matrices Ā, B, C, Ēf , and Ēd is thus considered:

       x(k + 1) = ( Ā -L C)x(k) + Bu c (k -l) + Ly(k) r(k) = Q y(k) -C x(k) (30)
where r ∈ R nr , x ∈ R nx is the residual, and the state estimation vector, respectively. The matrix Q ∈ R nr×ny is the residual weighting matrix.

The Z-transformed residual response to faults and unknown inputs is

r(z) = G rf (z)f (z) + G rd (z)d(z) (31)
where G rf (z) and G rd (z) denote the transfers between f (z) and r(z), and d(z) and r(z), respectively. Once Ēd is known, the remaining problem is to find matrices L and Q so that ( Ā -L C) is stable, and G rd (z) = 0 holds. The assignment of the observer eigenvectors and eigenvalues is a direct way to solve this design problem. Note that, because this technique does not consider a sensitivity constraint in the design procedure, therefore the fault sensitivity performance of the proposed FDI scheme can only be verified a posteriori. Especially, the subspace of considered faults should not intersect the subspace of decoupled disturbances [START_REF] Ding | A Unified Approach to the Optimization of Fault Detection Systems[END_REF], i.e., Im( Ēf ) ⊂ Im( Ēd ).

The following lemmas and theorems give the solution to the design problem.

Lemma 1. The transfer function G rd (z) can be expanded in terms of the eigenstructure as

G rd (z) = H(zI -Āc ) -1 Ēd = nx i=1 Hv i l T i Ēd z -λ i ( 32
)
where H = Q C, v i and l T i are the right and left eigenvectors of Āc = Ā -L C associated with eigenvalue λ i .

Lemma 2. A given left eigenvector l T i of Āc is always orthogonal to the right eigenvectors v j corresponding to the remaining (n x -1) eigenvalues λ j of Āc , where λ i = λ j .

Theorem 1 (Unknown input decoupling using left EA). If the necessary condition

Q C Ēd = H Ēd = 0 (33)
holds and all rows of the matrix H are left eigenvectors of Āc corresponding to n r eigenvalues of Āc , then G rd (z) = 0 is satisfied. The proof of this theorem can be found in [START_REF] Patton | Robust model-based fault diagnosis for dynamic systems[END_REF].

Following the above lemmas and theorem the first step for the design of the unknown input decoupled residual generator [START_REF] Pettazzi | Design of Robust Drag-Free Controllers with Given Structure[END_REF] is to compute the weighting matrix Q such that it satisfies Eq. ( 33).

The maximum row rank of Q is n y -rank( C Ēd ), thus the residual signal dimension should be chosen according to

n r ≤ n y -rank( C Ēd ) (34)
The second step is to determine the eigenstructure of the observer. All rows of H must be the n r left eigenvectors of Āc . The remaining n -n r left eigenvectors can be chosen without restraint.

For the given (stable) eigenvalue spectrum Λ( Āc ) = {λ i , i = 1, . . . , n x }, the following relation holds

l T i (λ i I -Ā) = -l T i L C = -m T i C, i = 1, . . . , n x ( 35
)
where m T i = l T i L. The assignability condition says that for each λ i , the corresponding left eigenvector l T i should lie in the column subspace spanned by { C(λ i I -Ā) -1 }, i.e., a vector m T i exists 21 such that

l T i = m T i K i , i = 1, . . . , n r ( 36
)
where

K i = -C(λ i I -Ā) -1 , i = 1, . . . , n r .
The projection of l i in the subspace span{K i } is denoted by

l •T i = m •T i K i , i = 1, . . . , n r ( 37
)
where

m •T i = l T i K T i (K i K T i ) -1 , i = 1, . . . , n r . If l T i = l •T i , l T i is in span{K i },
the required observer eigenstructure is assignable and perfect decoupling can be achieved. Otherwise, the eigenvectors must be chosen to be close, e.g., in a least-square sense l T i -l •T i , to the desired eigenvectors, i.e., an approximative procedure must be considered in order to replace l T i by its projection l •T i . In this situation, the residuals have low sensitivity to unknown inputs due to approximate decoupling [START_REF] Patton | Robust Fault Detection of Jet Engine Sensor Systems Using Eigenstructure Assignmen[END_REF].

The remaining n x -n r eigenvalues and corresponding eigenvectors can be chosen freely from the assignable subspace, e.g., using singular value decomposition (SVD). Then, the observer matrix L can be computed as follows

L = P -1 S ( 38
)
where

S = m • 1 . . . m • nr m nr+1 . . . m nx T P = l • 1 . . . l • nr l nr+1 . . . l nx T
It is obvious, that the first n r eigenvalues corresponding to the required eigenvectors l T i , i = 1, . . . , n r must be real because all these eigenvectors are real-valued.

Remark 1. Without loss of generality, it can be assumed that Ēd has a full column rank. When this is not the case, the following decomposition can be applied:

Ēd d(k) = Ēd1 Ēd2 d(k),
where Ēd1 is a full column rank matrix and Ēd2 d(k) can now be considered as a new unknown input.

Computational Results

First, the position model (Eq. ( 12)) and the attitude model (Eq. ( 15)) are transformed into the discrete form [START_REF] Karimi | A Linear Matrix Inequality Approach to Robust Fault Detection Filter Design of Linear Systems with Mixed Time-Varying Delays and Nonlinear Perturbations[END_REF] with l = 0 and m = 1. It practically means, that the unknown time-varying delay τ (k) is assumed to be in the closed interval [0, T ). Using the Cayley-Hamilton theorem given in proposition 1 and the 2nd order (h = 2) Taylor series expansion given in proposition 2, the uncertain time-varying delay τ (k) is modeled as an unknown input as in Eq. ( 28).

The resulting Ēd matrix has a large number of columns and the rank condition given by Eq. ( 34) cannot be explicitly satisfied. Choosing the desired residual dimension be equal to one, i.e., n r = 1, the following low rank factorization is performed for both models:

Ē * d = arg min Ēd - Ē * d 2 F , s.t. n r = n y -rank( C Ē * d ) (39)
By this factorization, the most significant directions are kept. Finally, a full column rank decomposition is performed on Ē * d using SVD decomposition.

The obtained distribution matrix is used for the residual generator design given by Eq. ( 30). The desired left eigenvectors of the observer are the rows of the matrix H = Q C where the weighting matrix Q is determined such that Eq. ( 33) is satisfied. In order to compare the FDI performances of both (position and attitude) models, the assigned eigenvalues (dynamics of the observer) were selected to be exactly the same for both models, i.e., Λ( Āc ) = {0.85, 0.87, 0.89, 0.91, 0.93, 0.95}.

Comments on Implementation Issues

In order to avoid using an optimization procedure to determine s max i and s min i in Eq. (49), the solutions s i (t), i = 1, . . . , n x of the differential equation (48) were found numerically, and therefore, s max i and s min i can be found using a simple iterative method. It is worth noting that other exact unknown input (disturbance) decoupling methods exist [START_REF] Patton | Robust model-based fault diagnosis for dynamic systems[END_REF]. In our particular case only the left EA technique appeared to be a viable candidate. Other methods, such as UIO or right EA technique, violated some necessary conditions of the solution existence.

B. Residual Evaluation and Fault Isolation

Once the residual generation problem is solved, the problem is to make a decision about the fault presence. The generalized likelihood ratio (GLR) test is used here [START_REF] Ding | Model-Based Fault Diagnosis Techniques: Design Schemes, Algorithms, and Tools[END_REF]. The decision is made based on two hypotheses: H 0 , the null hypothesis means no fault is present, while H 1 , the alternative hypothesis, indicates some anomaly in the system considered to be due to the thruster fault. In this case, the decision test J th is defined by

J th =        S N d (k) ≤ J th , H 0 is accepted S N d (k) > J th , H 1 is accepted (40)
where J th is a fixed threshold selected by the designer and S N d (k) is given by

S N d (k) = N d ln(σ 0 ) - N d 2 1 + ln(σ 2 1 (k)) - σ2 1 (k) σ 2 0 , σ2 1 (k) = 1 N d k j=k-N d +1 r 2 (j) ( 41
)
where r is the residual signal, σ 0 is the standard deviation of r in fault free situation, and

N d > 1
represents the detection sliding window due to on-line realization aspects. The interested reader can refer to the monograph of Basseville and Nikiforov [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF] for details on the threshold determination.

In this paper, a cross-correlation test between the residual r and the associated thruster open duration T oni is considered for fault isolation purposes. In the thruster configuration of Fig. 3, each thruster has its partner which provides same torque but force in exactly opposite direction.

Therefore, a full coverage of the isolation problem cannot be solved based on the so-called "sub-space isolation approach" because this test cannot distinguish between faults in either thruster, only in the thruster pair [START_REF] Massoumnia | A Geometric Approach to the Synthesis of Failure Detection Filters[END_REF]. Structured observer schemes could be possible candidates for thruster fault isolation [START_REF] Frank | Fault Diagnosis in Dynamic Systems Using Analytical and Knowledge-Based Redundancy: A Survey and Some New Results[END_REF]. They are based on making each residual signal sensitive to a subset of faults while being insensitive to another subset. This, however, requires a bank of observers to be designed and run in parallel. Due to the on-board computational limitations, the following considers an alternative solution that requires only one observer.

The proposed isolation strategy is based on minimum σ N or maximum σN cross-correlation criterion between the residual signal r and T oni , i.e.,

σ N (k) =        σ N (k) if fault case 1 or 3 σN (k) if fault case 2 or 4 (42)
where

σN (k) = arg i max 1 N + 1 k j=k-N +1 r(j)T oni (j) , i = 1...8, ∀k ∈ Z + (43) σ N (k) = arg i min 1 N + Ω N i (k) k j=k-N +1 r(j)T oni (j) , i = 1...8, ∀k ∈ Z + (44)
Ω N i (k) = 1 - k j=k-N +1 ϕ i (j), i = 1, . . . , 8, ϕ i (j) =        0 if T oni (j) = 0 1 if T oni (j) = 0
These cross-correlation functions are statistical quantities that try to find the associated thruster index that has the smallest/greatest impact on the resulting residual signal. For real-time reason, these criteria are computed on a N -length sliding-window. An increase in the value of N results in elongated isolation time delay. Hence again, an optimal value of N has to be selected. The resulting index σ N (k) ∈ {1, 2, . . . , 8} refers to the identified faulty thruster at time instance k.

Finally, the resulting thruster index is confirmed at time instant k, if the following holds:

σ N (k) = σ N (k -1) = . . . = σ N (k -N c + 1) (45)
where a confirmation window of length N c > 1 is introduced in order to avoid initial transition phenomena and to ensure robustness.

A key feature of this isolation strategy is that it is static, and thus, it has a low computational burden. Note that if the ith thruster is not used by the TMF, i.e., T oni = 0, the minimum crosscorrelation function will possibly result in σ(k) = i. This fact is taken into account by introducing a penalty function Ω N i (k) in [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF].

Remark 2. It should be noted that an event resolution algorithm must be implemented for the decision given by Eq. ( 42), since it is required to distinguish the fault cases 1-3 from the case 2-4.

Fortunately, because cases 1-3 imply a propellant overconsumption and since the chaser is equipped with a dedicated sensor that monitors the overall propellant consumption, this problem can be easily solved.

C. Computational Procedure

In order to ensure robustness, whilst being sensitive to faults, the threshold J th has to be selected carefully. A higher value of J th will increase the non-detection rate while a lower threshold will increase the false alarm rate. The optimal value of J th can be selected through Monte Carlo (MC) simulation. This approach is widely used in the FDI community to analyze the efficiency and performance of the designed algorithm [START_REF] Patton | Reliable Fault Diagnosis Scheme for a Spacecraft Attitude Control System[END_REF]. Here a set of 200 fault free MC simulations were performed using N d = 10 in Eq. 41, i.e., the notation S 10 is used. It can be clearly seen from Fig. 8, that for both models signal S 10 (k) does not exceed 20 in all the cases. Therefore, for the simulations in the next section, threshold of J th = 20 is chosen to ensure (ideally) a zero false alarm rate. For the isolation function σ N (k), a sliding window of N = 10 and confirmation window of N c = 15 samples was considered.

V. Simulation Results

The two FDI schemes described in the previous section are next implemented within the MSR "high-fidelity" industrial simulator, provided by Thales Alenia Space. All simulations are carried out under realistic conditions, i.e., the NAV unit is not considered to deliver "perfect" state estimates. We also assume constant time delay between the NAV and the control block, time-varying delays induced by the electronic devices, spatial disturbances (e.g., solar radiation pressure, gravity gradient, atmospheric drag) and parameter uncertainties. The simulations are all carried out during the last 20m of the rendezvous phase.

A set of n mc = 1600 Monte Carlo simulations, for each faulty case (4 × n mc ), have been run in order to assess the performances of the two proposed FDI schemes. For each run, model parameters, e.g., mass, center of mass (CoM), etc., were altered within a specific limit, see Table 1 for details.

The mass, CoM and inertia were scattered according to the normal distribution and truncated to the corresponding 3σ value. Uncertainties of thruster rise times and thruster misalignment phenomena are modeled by 1% uncertainty on thruster forces. Thruster faults were uniformly distributed among all of the 8 thrusters. Correspondingly, the leakage and the thrust loss size were drawn from the uniform distribution with the following range: mleak ∈< 10%, 30% > and mloss ∈< 40%, 90% >. The leakage is implemented as a dynamic lower saturation to the commanded thruster open rate, where this saturation starts at value 0 and ends at mleak with a slope of m s = 0.1. In all cases, fault occurs at time 1000s and is maintained. chosen faulty situations selected from the scattered parameter space listed in Table 1:

• A fault that corresponds to a stuck open valve, occurs in the thruster No.8., i.e., the thruster No. 8 is fully opened so that it provides a maximum thrust;

• A fully blocked-closed fault occurs in thruster No.3. In this case, the thruster does not generate any thrust regardless of the command by the TMF;

• The third faulty situation corresponds to thruster No.2 suffering from a leakage of size 19.2%;

• The fourth faulty situation corresponds to the case when the thruster No.7 loses its thrust level by a value of 54.7%. Clearly, the nonlinear simulations show that faults are detected and isolated by the proposed FDI units within a reasonable time. Moreover, Fig. 10 shows the ability to detect and isolate small (thrust loss) and incipient (residual leakage) thruster faults.

To evaluate the performance and reliability of the two FDI schemes, some statistical indices have been used, i.e., the detection time t d (time from fault occurrence to fault detection) and isolation time t i (time from fault occurrence to fault isolation) are computed. The considered indices are listed below.

• mean(t d )/mean(t i ) -mean detection/isolation time in seconds;

• std(t d )/std(t i ) -standard deviation of the detection/isolation time in seconds;

• nondetection rate P nd -number of non-detections divided by n mc ;

• true isolation rate P i -number of correctly isolated thrusters divided by n mc . These performance indices are calculated for each fault scenario and model separately. since no false alarms P f = 0 and non-detections P nd = 0 have been revealed. Furthermore, all thruster faults were correctly isolated at the end of the campaign, i.e., P i = 1. These achieved results also demonstrate that each FDI scheme is able to successfully detect and isolate the chaser thruster faults.

P nd /P f /Pi 0/0/1 0/0/1 0/0/1 0/0/1 0/0/1 0/0/1 0/0/1 0/0/1
In order to better appreciate the results presented in Table 2, a histogram plot is performed to graphically represent the distribution of the detection t d and isolation t i times, respectively.

Figure 11 shows a comparison of the obtained results in terms of detection times and Fig. 12 in terms of isolation times. This visual representation allows to evaluate the FDI performances in terms of minimum and maximum detection/isolation times, as well as to observe the median values.

It can be seen from Fig. 11, that (as expected) the FDI unit based on the attitude model presents a greater sensitivity towards all faulty situations. This can be easily explained by the fact that attitude is more sensitive to small thruster faults.

Note that the occurrence of incipient or small size thruster faults (e.g., small propellant leakage or thrust loss) may be covered by (robust) control actions, and the early detection of them is clearly The linearity of the attitude model during the fault presence is questionable.

VI. Conclusion

Design and implementation of two distinct model-based FDI strategies for thruster fault diagnosis of an autonomous spacecraft involved in the rendezvous phase of the MSR mission has been presented. Effects of unknown time-varying delays induced by the propulsion drive electronics and uncertainties on thruster rise times has been transformed into unknown input vector using Cayley-Hamilton theorem and h-order Taylor expansion. The estimation of the complex unknown input distribution matrix can be considered as a contribution to the theory. The robust (in the sense of unknown input decoupling) residual generation for FDI has been achieved using observer eigenstructure assignment technique, i.e., some left eigenvectors of the observer have been assigned to be orthogonal to the unknown input directions (columns of the distribution matrix). Thruster fault isolation has been achieved by cross-correlation test between the residual signal and the commanded thruster open durations. The core element of this study is the judiciously chosen position model. This model has been used to design the first FDI scheme and was compared (in terms of well established FDI performance indices) to the second, attitude model-based, scheme. A Monte Carlo simulation campaign, using a "high-fidelity" industrial simulator, has been performed under realistic conditions, considering imperfect navigation, delays, spatial disturbances and parameter uncertainties. Four different fault scenarios were injected throughout simulations. Obtained results indicate, that the proposed FDI strategies are effective and applicable for on-board implementation. Moreover, the selected performance indices reveal, that the position model-based scheme tends to achieve very similar FDI performances as the scheme based on the pure attitude model. The proof of this theorem can be found in [START_REF] Leonard | The Matrix Exponential[END_REF].

Using (47), we have 25) given in proposition 2. This proposition is proved in [START_REF] Hetel | Stabilization of Arbitrary Switched Linear Systems with Unknown Time-Varying Delays[END_REF].

Γ a 1 (δ(k)) =

  null matrix with appropriate dimensions In = n × n identity matrix I chs = inertia matrix (kg•m 2 ) M = thruster configuration matrix m = mass (kg) mloss , mleak = thrust loss size and maximum leakage size n = orbital rate of the target (rad•s -1 ) p, q, r = roll, pitch, and yaw rate (rad•s -1 ) Qtgt, Qchs = quaternion estimates describing the orientation of the target and chaser with respect to the Mars-centered inertial frame r, r = vector of residuals and scalar residual T = sampling interval (s) Ton i = scaled open duration of the ith thruster X, Y , Z = mutually perpendicular (orthogonal) unit vectors x (bold lowercase letters) = vectors ξ, η, ζ = Cartesian component of the relative position (m) ϕ, θ, ψ = roll, pitch, and yaw angle (rad) Λ(A) = set of all eigenvalues of a matrix A τ (t) = time-varying delay (s) R, Z + = field of real numbers and set of non-negative integers R l , Ri, R b = local (target), Mars-centered inertial, and chaser body-fixed reference frame • , • F = Euclidean norm (vectors) and Frobenius norm (matrices)

Fig. 1

 1 Fig. 1 Illustration of the MSR mission.

Fig. 2

 2 Fig. 2 General setup of the chaser's GNC system during the rendezvous.

Fig. 3

 3 Fig. 3 Thruster configuration of the chaser spacecraft.

Fig. 4

 4 Fig. 4 Control architecture of the thruster-based propulsion system of the chaser vehicle.

Fig. 5

 5 Fig. 5 The MSR rendezvous corridor.

  where ξ, η, ζ denote the elements of the three dimensional relative position vector ρ = [ξ, η, ζ] T of the chaser from the origin of the target frame O T , expressed in R l ;

  4 and Qchs ∈ R 4 denote the attitude quaternion estimate of the target, and the chaser, respectively. The attitude estimate Qtgt and Qchs describe the orientation of the target body frame and the chaser body frame with respect to R i , respectively. These estimates are provided by the NAV unit. The quaternions dependent rotation matrix R(•) performs the projection of the threedimensional force vector F f thr from the chaser's frame on to the target frame R l . The measurementvector y = [ξ η ζ] T is expressed in R l andmeasured by LIDAR. Spatial disturbances (solar radiation pressure, gravity gradient and atmospheric drag) are denoted by w ∈ R 3 . A p , B p , C p and E pw are matrices of appropriate dimension.

(

  see Fig 7 for illustration).

Fig. 7

 7 Fig. 7 Chaser target pointing.

1 ,

 1 ..., h + 1, ∀k ∈ Z + are uncertain scale factors satisfying h+1 i=1 µ b i (k) = 1, ∀k ∈ Z + , and U b i , i = 1, . . . , h + 1 are known constant matrices given in Appendix B.

Fig. 8 A

 8 Fig. 8 A set of 200 Monte Carlo simulations for fault-free case.

Figures 9 and 10 Fig. 9

 109 Figures 9 and 10 illustrate the behaviour of the most important characteristics of the FDI units

Figure 9 7

 97 Figure 9 is concerned with the two first situations whereas Fig. 10 considers the two last cases.

Fig. 11

 11 Fig. 11 Histograms of the detection times for the considered fault cases.

Fig. 12

 12 Fig. 12 Histograms of the isolation times for the considered fault cases.

Appendix A : i

 : Uncertainty Transformation Using Cayley-Hamilton TheoremLet us first consider the following theorem[START_REF] Leonard | The Matrix Exponential[END_REF]:Theorem 2. The characteristic polynomial of matrix A is p(λ) = det(λI -A) = λ nx + c n-1 λ nx-1 + . . . + c 1 λ + c 0 (46)then e At can be written ase At = s 1 (t)I + s 2 (t)A + . . . + s nx (t)A nx-1 (47)where s i (t), 1 ≤ i ≤ n x are solutions to the nth order homogenous scalar differential equation s (nx) (t) + c nx-1 s (nx-1) (t) + . . . + c 1 s (t) + c 0 s(t) (0) = 0 for j = i -1, 0 ≤ j ≤ n x -1

TΦ 1 =Φ 2 =

 12 -δ(k) s i (t)dt, i = 1, 2, . . . , n x s i (t)dt, i = 1, 2, . . . , n x then (49) can be rewritten asΓ a 1 (δ(k)) = nx i=1 α i,0 (k)s min i + α i,1 (k)s max i A i-1 B (50)where α i,0 (k) and α i,1 (k) are two time-varying unknown parameters satisfying 0≤ α i,0 (k) ≤ 1, 0 ≤ α i,1 (k) ≤ 1, and α i,0 (k) + α i,1 (k) = 1 for ∀k ∈ Z + . It can be verified that T T -δ(k) s i (t)dt, i = 1, 2, . . . , n x are Lipschitz-continuous on 0 ≤ δ(k) ≤ T , i.e., they satisfy T T -δ1(k) s i (t)dt -T T -δ2(k) s i (t)dt ≤ κ i |δ 1 (k) -δ 2 (k)| ∀δ 1 (k), δ 2 (k) ∈ [0, T ],where κ i , i = 1, 2, ..., n x are the Lipschitz constants.Settingµ a 2i-1 (k) = α i,0 (k)/n x , µ a 2i (k) = α i,1 (k)/n x , U a 2i-1 = n x s min i A i-1 B U a 2i = n x s max iwhereG i = (-1) i+1 A i-1 i! e AT B, i = 1, . . . , h ε h I ε h-1 I , ..., ε 2 I εI T ε h I ε h-1 I , ..., ε 2 I εI T . . . Φ h+1 = ε h I ε h-1 I , ..., ε 2 I εI T then using h-order Taylor expansion, Eq. (52) yields Eq. (

Table 1 Parameter uncertainties of the chaser spacecraft (3 sigma)

 1

	Parameter	Variation range	Unit
	Mass	±10%	[kg]
	Inertia	±20%	[kg.m 2 ]
	Thrusters forces	±1%	[N ]
	CoM	±3	[cm]

Table 2 FDI performances based on 4x1600 Monte Carlo runs FDI based on position model FDI based on attitude model Criterion CASE 1 CASE 2 CASE 3 CASE 4 CASE 1 CASE 2 CASE 3 CASE 4

 2

	mean(t d )	1.6	16.3	3.6	18.8	1.4	11.7	2.9	12.2
	std(t d )	0.017	3.803	0.167	6.594	0.048	2.895	0.076	3.064
	mean(ti)	3.4	17.9	5.5	20.4	3.6	15.0	5.5	15.6
	std(ti)	0.207	3.799	0.399	6.592	0.237	3.231	0.459	3.456

Table 2

 2 presents complete results obtained from the Monte Carlo simulation campaign. As it can be seen from this table, the two proposed FDI schemes present good reliability characteristics

Acknowledgments

This research work was supported by the European Space Agency (ESA) and Thales Alenia Space France in the frame of the ESA Networking/Partnering Initiative (NPI) program.

then Eq. ( 24) in proposition 1 yields Eq. (50).

Appendix B: Uncertainty Transformation Using Taylor Series Expansion

Theorem 3. Taylor series expansion of the uncertainty Γ 1 (δ(k)) is given by

The h-order approximation of Taylor expansion for the uncertainty Γ 1 (δ(k)) can be expressed as a finite sum of the first h elements


M Bichon

J Le Tanou

M Ropert

W Hamidouche

L Morin

L Zhang

INTER-BLOCK DEPENDENCIES CONSIDERATION FOR INTRA CODING IN H.264/AVC AND HEVC STANDARDS

Keywords: Inter-Block Dependencies, Joint Rate-Distortion Optimization, Intra Coding, H.264/AVC, HEVC

Recent MPEG video compression standards are still block-based: blocks of pixels are sequentially coded using spatial or temporal prediction schemes. For each block, a vector of coding parameters has to be selected. In order to limit the complexity of this decision, independence between blocks is assumed, and coding parameters are locally optimized to maximize the coding efficiency. Few studies have investigated the benefits of inter-block dependencies consideration using Joint Rate-Distortion Optimization (JRDO), especially in Intra coding. To the best of our knowledge, maximum achievable gains of such approaches have never been exhibited. In this paper, we propose two JRDO models performing joint optimization of multiple blocks applied to intra prediction mode decision. The proposed models have been evaluated in both H.264/AVC and HEVC standards. These two models enables a bitrate saving with respect to the classical RDO model up to -3.10% and -2.31% in H.264/AVC and HEVC, respectively.

INTRODUCTION

Rate Distortion Optimization (RDO) [START_REF] Sullivan | Rate-distortion optimization for video compression[END_REF] is a well-known technique used to decide the best set of coding parameters which maximizes coding efficiency. It consists in minimizing the distortion (D) under a bitrate (R) constraint. This minimization under constraint can be rewritten using the Lagrange multiplier method [START_REF] Wiegand | Lagrange multiplier selection in hybrid video coder control[END_REF]. Knowing λ the Lagrangian Multiplier (LM), it leads to minimization of the R-D cost function:

J = D + λ × R (1)
In block-based video compression standards like H.264/AVC [START_REF] Wiegand | Overview of the H.264/AVC video coding standard[END_REF] and High Efficiency Video Coding (HEVC) [START_REF] Sullivan | Overview of the High Efficiency Video Coding (HEVC) Standard[END_REF], the RDO technique is applied sequentially to groups of pixels, defined as Coding Units (CU)s or blocks. One denotes CUi the i th coded CU, pi the vector defining the set of coding parameters to estimate for CUi and Ji the local R-D cost of CUi. Then, the global solution is obtained by

min { p k } N -1 k=0 N -1 i=0 Ji | { p k } N -1 k=0 (2)
with N the number of CU within the video sequence, and { p k } N -1 k=0 representing all the set of coding parameters from CU0 to CUN-1. Without loss of generality, (2) can be rewritten as (3) by using the causality of the block-based compression scheme.

min { p k } N -1 k=0 N -1 i=0 Ji | { p k } i k=0 (3)
In this paper, one denotes (3) the global-RDO. The complexity of solving (3) exhaustively is tremendous and unrealistic. By writing K the vector space dimension for pi, K N would be the number of configurations to be tested.

Equation ( 3) is usually simplified based on the common assumption of independence between CUs: it assumes that Ji is only dependent on pi. Under this hypothesis, one defines the classical-RDO minimization problem as

min { p k } N -1 k=0 N -1 i=0 Ji | p i = min p i Ji | p i N -1 i=0 (4)
Classical-RDO reduces the number of tested configurations to K × N . Although the independence assumption significantly simplifies the computational complexity, we have not found any evaluation of its impact on the coding efficiency in the literature. Some studies [5∼10] investigated joint optimization over multiple CUs, considering spatial or temporal dependencies inherent to the prediction scheme. These approaches are noted Joint Rate-Distortion Optimizations (JRDOs), which differs in this paper from the joint optimization of pi components.

Regarding Inter coding, Yang and al. [START_REF] Yang | Source Distortion Temporal Propagation Model for Motion Compensated Video Coding Optimization[END_REF] proposed an inter-frame distortion propagation model using source motion estimation which, along with local delta-quantization, significantly increases the coding efficiency. Li and al. [START_REF] Li | Inter-frame dependent rate-distortion optimization using lagrangian multiplier adaption[END_REF] formalized the temporal dependency between two blocks as a rate function, dependent of distortion made on the reference signal used for prediction, noted as Ri+1(Di).

The JRDO in intra coding has also been studied. In [START_REF] Pang | Optimal distortion redistribution in block-based image coding using successive convex optimization[END_REF][START_REF] Pang | Optimal dependent bit allocation for AVS intra-frame coding via successive convex approximation[END_REF], Pang and al. used successive convex optimizations in order to reduce block boundaries distortions and thus to improve intra coding efficiency. Sun and al. [START_REF] Sun | Modified distortion redistribution problem for High Efficiency Video Coding (HEVC)[END_REF] achieved similar objective by providing an optimized quantization matrix, based on the assumption that coefficients distortion does not follow a uniform distribution in intra coded blocks. In [START_REF] Qingbo | A Novel Joint Rate Distortion Optimization Scheme for Intra Prediction Coding in H. 264/AVC[END_REF], Qingbo and al. experimentally estimated a linear distortion propagation model used for λ computation, coupled with off-line learning and multiple LM framework.

Nevertheless, these studies are using coarse assumptions resulting in simplified dependency models. In addition, they often use a two-in-one algorithm which makes difficult to exhibit gains brought by a JRDO standalone strategy. Finally, to the best of our knowledge, there is no reference proving the maximum achievable gain of a coding decision model that considers intrinsic inter-block dependencies.

In this paper, we propose to evaluate the maximum achievable gain of exhaustive joint optimization of multiple CUs applied to intra prediction mode decision. A brief review of intra coding is first given in Section 2. Section 3 identifies dependencies inherent to intra coding scheme and introduces two JRDO models: Dual-JRDO and Quad-JRDO. Experimental results and bitrate savings of the proposed JRDO approaches are presented and discussed in Section 4 for both H.264/AVC and HEVC. Finally, Section 5 concludes this paper.

INTRA CODING REVIEW

MPEG intra coding scheme can be decomposed into two main steps: mode prediction decision and block partitioning decision. The two decisions steps, for HEVC, are summarized in the following subsections. For a complete explanation on intra coding in HEVC standard, the reader is referred to [START_REF] Lainema | Intra Coding of the HEVC Standard[END_REF].

Intra Mode Prediction

Intra prediction consists in sequentially predicting the source signal from the previously reconstructed pixels within the same frame, used as reference. Prediction is built by copying or interpolating reference pixels onto target pixels, according to a rule specified by the predictor. In the case of intra coding, reference pixels used to predict the current CU are depicted in Figure 1.

In HEVC, 35 possible predictors for intra coding are available and given in Figure 1. The DC-mode uses average value of reference pixels and the Planar-mode is a bilinear interpolation designed to preserve continuities along block boundaries. The 33 angular modes represent a direction of projection for reference pixels.

Fig. 1. Intra Prediction in HEVC

The H.264/AVC standard uses the same prediction process, but the number of spatial predictors is limited to 9, or 4 in case of 16x16 size blocks.

Block Partitioning

In HEVC, each frame is uniformly partitioned in Coding Tree Units (CTU), equivalent to Macroblocks (MBs) in H.264/AVC. CTUs are sequentially compressed in a raster scan order. Then, each CTU can recursively be further sub-divided in multiple CUs, following a QuadTree structure. Figure 2 (a) shows an example of the partitioning of a 64x64 CTU in HEVC.

In HEVC, each CU at a given depth of the QuadTree can be compressed as a 2Nx2N or NxN partition, with N ∈ {32, 16, 8, 4}. The block partitioning decision can be summarized to decide for each CU whether to split or not, based on its R-D cost. The QuadTree structure implies a Z-scan order to process successive CUs, an example is shown in Figure 2 The next Section describes dependencies inherent to the intra coding scheme and the two proposed JRDO models.

PROPOSED JRDO MODELS

We define the two Inter-block dependencies inherent to the intra coding scheme. In H.264/AVC and HEVC, both predictive and contextual entropy coding are used.

The compression scheme implies a first dependency between the prediction efficiency on the current CUi and the distortion made on reference samples D ref i , i.e. spatial reference neighbors previously coded.

The entropy coding step is usually based on Context Adaptive Binary Arithmetic Coding (CABAC). Each Syntax Element (SE) is binarized and lossless coded, based on context adaptation and occurrence probability updates of each SE bin. Hence, by definition, the syntax coding efficiency of the current CUi is directly impacted by previous decisions and dependent of SE contexts updates after processing CUi-1. The syntax context state before to code CUi is noted SCi.

Once the inter-block dependencies SCi and D ref i are defined, we consider Ji | p i (0, 0), the R-D cost of CUi knowing pi, if no dependencies interfere, either in term of distortion or CABAC. Then, one writes the JRDO equation as follows:

Ji | p i = Ji | p i (0, 0) + ∆Ji | p i (SCi, D ref i ) (5)
where ∆Ji(SCi, D ref i

) is the intra dependency propagation cost that [START_REF] Pang | Optimal distortion redistribution in block-based image coding using successive convex optimization[END_REF][START_REF] Pang | Optimal dependent bit allocation for AVS intra-frame coding via successive convex approximation[END_REF][START_REF] Qingbo | A Novel Joint Rate Distortion Optimization Scheme for Intra Prediction Coding in H. 264/AVC[END_REF] have tried to model. In order to simplify implementations and to keep the computational complexity reasonably low for the study, the joint optimization is limited to intra prediction modes. Hence, the vector of coding parameters pi becomes the scalar pi.

For a CU to code with a JRDO strategy, the difficulty consists in how to consider neighboring CUs: either in terms of spatial distance, coding order, or both? For example, we can refer to Figure 2 (b), we see that CU9 has strong spatial dependencies with CU1. Optimizing both of them jointly should significantly improve the coding efficiency of both CUs. Nevertheless, compression is applied sequentially and the CUs 2 to 8 have to be compressed first, in order to have the correct coding context SC9 when deciding and coding CU9.

This observation highlights the difficulty of implementing JRDO approaches in intra coding: on the one hand jointly optimizing many CUs is too complex to compute, and on the other hand ignoring seven CUs leads to strong approximations on coding context. Models proposed in Sections 3.1 and 3.2 avoid these difficulties and apply an exhaustive joint optimization of the current CU with neighbor CUs taking into account both spatial and coding order distances.

Dual-JRDO

Z-scan is the coding order used to sequentially encode the four sub-CUs resulting from split. Based on Z-scan and local dependencies, SCi and D ref i , one predicts that CUi+1 is highly dependent of CUi if CUi+1 is the spatial right neighbor of CUi.

To confirm this, we experimentally verified that the right neighbor CU is often the most dependent on the current CU to code. In consequence, we propose the Dual-JRDO model that jointly optimizes the intra prediction mode of each CU with the prediction mode of its right neighbor. In order to avoid wrong syntax context states, Dual-JRDO handles two cases:

• If CUi right neighbor is CUi+1, apply the Dual-JRDO.

• Else, classical-RDO is applied on the current CUi.

With p *

i the chosen predictor to encode CUi and p i+1 the estimated optimal coding mode for CUi+1, Dual-JRDO solution is defined by

p * i , p i+1 = arg min {pi,pi+1} Ji |p i + Ji+1 |p i ,p i+1 (6)
In HEVC, the neighboring CUs can be further split, leading to p i+1 = p * i+1 . To overcome this problem, one considers that two CUs coming from the same split process have a high probability to have the same final partition size. Statistically, we can notice that the probability of this assumption to be true increases as CU size decreases. An example of Dual-JRDO is shown in Figure 3 (a), with dotted lines delimiting the optimization area and dark gray area refers to block coded independently using classical-RDO.

In Dual-JRDO, K 2 possibilities are explored for half of the CUs, whereas in classical-RDO, only K possibilities are explored for all CUs. We deduce that in this particular case Dual-JRDO multiplies the complexity of classical-RDO by K/2.

In H.264/AVC, block partitioning necessarily splits CUs into partitions of the same size, thus we have p i+1 = p * i+1 leading to complexity reduction. Results of Dual-JRDO for both H.

Quad-JRDO

We can see from Figure 1 that numerous spatial predictors exploit vertical spatial correlations. By definition, Dual-JRDO does not consider distortion propagated vertically, e.g. distortion on bottom samples of CUi will not affect CUi+1 if it is the right neighbor. Quad-JRDO proposes to include vertical neighbors of CUi in the joint optimization process.

Quad-JRDO optimizes all sub-CUs coming from the same split operation. At CTU level, raster scan order imposes to code the whole line of CUs before reaching bottom neighbors of CUi. This results in unachievable computational complexity or wrong syntax context states, reason why Quad-JRDO is not applied at CTU nor MB level. Equation ( 7) depicts the optimization formulation.

{p * k } i+3 k=i = arg min {p k } i+3 k=i i+3 k=i J k | {p l } k l=i (7)
Quad-JRDO supposes that all sub-CUs will not be further split, which is not matching the QuadTree structure in HEVC. To overcome this issue, ( 7) is only applied to the special case of NxN mode, other cases use the classical-RDO. An example of Quad-JRDO is shown in Figure 3 (b). The use of Quad-JRDO for NxN analysis multiplies the complexity of classical-RDO by 35 3 /4.

In H.264/AVC, MBs are always split in sub-partitions of same size, either 8x8 or 4x4 partitions. In the case of 8x8 partition mode, coding parameters of the four 8x8 blocks are optimized jointly. In the case of 4x4 partition mode, each 8x8 block optimizes jointly the four 4x4 sub-blocks. Experimental configurations and R-D results for both standards are presented and discussed in Section 4.2.

EXPERIMENTS

The two methods exposed in Section 3 have been implemented into HEVC and H.264/AVC reference test models, HM16.6 [START_REF] Mccann | JCTVC-R1002: High Efficiency Video Coding[END_REF] and JM19.0 [START_REF] Tourapis | JVT-AE010: H.264/14496-10 AVC Reference Software Manual[END_REF] respectively. The set of sequences utilized is picked among JCT-VC test set [START_REF] Bossen | Common test conditions and software reference configurations[END_REF]. Furthermore, because of the computational complexity of our methods, one encodes one frame of each sequence, in I-only configuration. Indeed, this study aims to estimate maximum gains of JRDO strategies, and not to propose low-complexity solutions.

Tables 1 and2 present coding efficiency improvements of Dual-JRDO and Quad-JRDO against classical-RDO. Results use Bjøntegaard metric [START_REF] Bjontegaard | Calculation of average PSNR differences between RD-curves[END_REF] and are expressed in BD-BR, i.e. the percentage of bitrate savings to achieve similar distortion, measured as frame PSNR. Even if initial metric is expressed using 4 different QP values, we use it with 5 QP values (QP ∈ (22, 27, 32, 37, 42)) to cover a larger range of bitrates. Besides, since the proposed solutions are mainly used to optimize luminance (Y) encoding, we focuses on Y BD-BR, nevertheless, similar gains have been obtained in YUV 4:2:0.

Results on Dual-JRDO

Results of Dual-JRDO are depicted in table 1. We observe constant gains against classical-RDO. Average bitrate savings are of -0.77% and of -0.71% in JM19.0 and HM16.6, respectively. Dual-JRDO outperforms classical-RDO up to more than -1.3% in both reference softwares. However, one observes that the BasketballPass sequence in JM19.0 is the only one to present negligible losses. Dual-JRDO slightly favors horizontal predictions. In few cases where vertical prediction is better than horizontal, Dual-JRDO can slightly penalize coding efficiency.

Table 2. Y BD-Rate of Quad-JRDO in JM19.0 and HM16.6

CONCLUSION

In this paper, one demonstrates the benefits of considering interblock dependencies in intra coding, for both H.264/AVC and HEVC standards. For intra coding, we identify two inter-block dependencies coming respectively from the distortion of reference pixels D ref i and the CABAC coding context SCi. Then, we propose an exhaustive JRDO scheme for intra prediction of multiple CUs in order to exhibit maximum achievable gains of dependencies consideration. For a first model, focused on dependencies with the right-CU only, one observes on average -0.77% (H.264/AVC) and -0.71% (HEVC) improvements in BD-rate against classical-RDO. In a second model, one includes up to four CUs into the joint optimization process. This last model outperforms on average classical-RDO by -1.78% (H.264/AVC) and -1.47% (HEVC). Despite complexity of proposed JRDO methods, we demonstrate that significant gains can be achieved with such a strategy.

Our future work will focus on the modelization of the two identified dependencies in order to propose a low-complexity JRDO scheme for intra prediction of multiple CUs, in the particular context of HEVC. Results of this study will be used for comparison. The QuadTree structure and the number of possible partitioning allowed within a CTU will be the main difficulties to overcome in HEVC. Besides, extension to this work to temporal dependencies across frame will be also envisaged.

  (b).

Fig. 2 .

 2 Fig. 2. Intra Partitioning in HEVC

Fig. 3 .

 3 Fig. 3. Example of Dual-JRDO and Quad-JRDO

The results presented in [START_REF] Qingbo | A Novel Joint Rate Distortion Optimization Scheme for Intra Prediction Coding in H. 264/AVC[END_REF] are coming from two separate solutions, the first contribution is related to a JRDO approach and the second contribution to a multiple LM framework. In their JRDO approach, similar dependencies as Dual-JRDO are considered and gains announced for video of 1920x1080 resolution are about -0.13%. Our study on identical test set shows that achievable gains are on average -0.80% for H.264/AVC. They estimate the distortion dependency with an off-line linear distortion propagation model, and analytically deduct the related optimal λ. Our exhaustive joint prediction optimization demonstrates there is room for improvement in modelization of dependencies or cost propagation.

Results on Quad-JRDO

In the case of the Quad-JRDO model presented in Section 3.2, the optimization is applied only to 4x4 and 8x8 blocks in JM19.0, and NxN case in HM16.6. The remaining decisions are based on classical-RDO. Results for both implementations are presented in Table 2.

As expected, much higher gains are observed with this second model, which is also much more complex. In average, bitrate savings over classical-RDO are -1.78% in JM19.0 and -1.47% in HM16.6. BasketballPass sequence, for which negligible losses were observed in Dual-JRDO, now outperforms classical-RDO from -1.09% in Quad-JRDO. Besides, we must note that the Quad JRDO is less efficient than Dual JRDO on some high resolution sequences such as Kimono and BasketballDrive. One explanation is that some HD sequences may have more homogeneous areas, where larger partitions are preferred for the prediction; the joint optimization of NxN (i.e. 4x4) partition mode is then of limited interest for these particular cases.

The significant coding efficiency improvement between Dual and Quad-JRDO mostly comes from the consideration of vertical predictions and 2-D spatial dependency. Based on these results, it seems relevant to assume that adding more CUs in the proposed joint optimization process, would bring much more gain. One could expect to tend towards global-RDO efficiency. In practice, the complexity of such process would lead to computationally intractable simulations.