I T Castro

K T Huynh

A Barros

C Bérenguer

A predictive maintenance strategy based on mean residual life for systems subject to competing failures due to degradation and shocks

Keywords: Mean residual life, Condition-based maintenance, Degradation-threshold-shock model, Stochastic processes

In condition-based maintenance area, mean residual life function is recognized as a promising reliability indicator for maintenance decision-making. However, very few maintenance strategies based on the mean residual life have been developed. This paper therefore analyzes a predictive maintenance strategy based on the mean residual life given the degradation level collected periodically in different inspection times. This strategy is developed for systems which are subject to competing and dependent failures due to degradation and traumatic shocks. Numerical examples are given to illustrate the analytical results of the strategy.

INTRODUCTION

Condition based maintenance is a maintenance program that utilizes condition monitoring technologies in order to construct a reliability indicator representing the health state of a system and, through that knowledge, set up proper maintenance decisions. Synthesizing, the literature shows that a robust reliability indicator has to characterize precisely not only the current, but also the future state of systems. The system degradation level is a traditional indicator, based on which numerous maintenance models have been analyzed and their performance has been proved through practical applications (see e.g., [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF] and [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF]). However, since the degradation level can only reflect the current system state, it does not suffice to be a robust indicator. From the degradation level, one can estimate the mean residual life (MRL) of system. Since the MRL can provide an idea about how long a system at any particular age and deterioration is expected to survive, it can also reflect the system future state. The MRL is therefore considered as a very promising indicator for maintenance decision-making. This is why we seek in the present paper a maintenance strategy based on this reliability indicator.

It is well known that most of the industrial systems degrade physically over time, and a stochastic process is often used to describe this phenomenon [START_REF] Abdel-Hameed | Degradation processes: An overview[END_REF]. But, considering a degradation process only seems to be unsatisfactory for modeling dynamic systems, because they can fail due to another causes rather than the degradation [START_REF] Singpurwalla | Survival in dynamic environments[END_REF]. When the system is subject to the competing risks of degradation and shocks, one can show that Degradation-Threshold-Shock (DTS) models are relevant for the modeling purpose [START_REF] Lehmann | Joint modeling of degradation and failure time data[END_REF]. For these models, the system is regarded as failed when its degradation reaches a critical threshold or when a catastrophic shock occurs although the degradation process has not reached the threshold. The "classical" DTS models assume that the degradation process and the shock process are independent [START_REF] Li | An inspection-maintenance model for systems with multiple competing processes[END_REF][START_REF] Van Noortwijk | Gamma processes and peaks-over-threshold distributions for time-dependent reliability[END_REF]. However, in many practical situations, the dependence between them is of importance and should not be neglected in the modeling [START_REF] Wang | Dependent competing-risk degradation systems[END_REF]. Recently, Huynh et al. provide in [START_REF] Huynh | A periodic inspection and replacement policy for systems subject to competing failure modes due to degradation and traumatic events[END_REF] and [START_REF] Huynh | Modeling age-based maintenance strategies with minimal repairs for systems subject to competing failure modes due to degradation and shocks[END_REF] an extension of DTS models by introducing the dependence of the shock process on the degradation process. The authors show that such a model can be seen as a combination and a more versatile and hopefully realistic extension of many classical failure models based either only on degradation or only on parametric lifetime distributions. Hence, we will develop our MRL-based maintenance strategy under the assumption of this kind of model.

The performance of the considered MRL-based maintenance strategy is assessed by a cost criterion which is the long-run expected maintenance cost rate C ∞ . We assume that a maintenance action, either preventive or corrective, restores the system to an as-good-as-new condition. The behavior of the maintained system is then assimilated to a regenerative process. As such, by applying the wellknown renewal-reward theorem [START_REF] Tijms | A first course in stochastic models[END_REF], the long-run expected maintenance cost rate can be calculated as the expected cost in a renewal cycle divided by the expected length of the cycle

C ∞ = lim t→∞ E [C (t)] t = E [C (S)] E [S] , (1)
where C (t) is the cumulated maintenance cost at time t and S is the length of the first renewal cycle.

The remainder of this paper is organized as follows. The detailed description of the maintained system according to the MRL-based maintenance strategy is introduced in section 2. Section 3 is devoted to computing the mean residual life function and to analyzing its monotony. In section 4, we formulate the cost model of the MRL-based maintenance strategy for the considered system. Section 5 shows the maintenance strategy optimization and some numerical examples. Finally, section 6 concludes the paper.

SYSTEM DEFINITION AND MAINTENANCE STRATEGY

We consider a system subject to the competing and dependent failure modes due to shocks and degradation. Such a system, under the proposed predictive maintenance strategy based on MRL, is described as follows.

1. The system is subject to an accumulative degradation and starts working at time 0. Let X(t) be the degradation level of the system at time t with X (0) = 0. We assume that {X (t) , t ≥ 0} follows a homogeneous gamma process whose increment X (t) -X (s) between two times s and t, 0 ≤ s < t, follows a Gamma density function with shape parameter α • (ts) and scale parameter β

f α•(t-s),β (x) = 1 Γ (α • (t -s)) β α•(t-s) x α•(t-s)-1 e -βx , x ≥ 0,
where Γ denotes the gamma function given by

Γ(α) = ∞ 0 u α-1 e -u du.
We assume that the system fails when its degradation level exceeds the value L. Let σ L be the hitting time of the threshold L (i.e., failure time due to degradation), Fσ L and f σ L be its survival function and its density function respectively, then

Fσ L (t) = P (σ L ≥ t) = P (X (t) ≤ L) = 1 - Γ(αt, Lβ) Γ(αt) , t ≥ 0, (2)
f σ L (t) = α Γ(αt) ∞ Lβ {log (z) -ψ (αt)} z αt-1 e -z dz, (3)
where Γ (α, x) = ∞ x z α-1 e -z dz and ψ (α) = Γ ′ (α) Γ(α) are respectively the Gamma function and Digamma function.

2. The system also fails whenever a shock occurs. We assume that the shock arrival times are modeled by a non-homogeneous Poisson process {N s (t), t ≥ 0} with stochastic increasing intensity depending on both degradation level X (t) and working time t

r (t, X (t)) = r 1 (t) 1 {X(t)≤Ms} + r 2 (t) 1 {X(t)>Ms} , (4)
where r 1 (t) and r 2 (t) denote two non-decreasing failure rates and r 1 (t) ≤ r 2 (t) for t ≥ 0. The quantity M s represents a fixed deterioration threshold. The expression 4 means that the degradation evolution affects to the occurrence of traumatic shocks in the sense that the system is more prone to shocks when the degradation increases and exceeds a given level. Let F1 (t) and F2 (t) be the survival functions associated to the failure rates r 1 and r 2 , that is,

Fi (t) = exp - t 0 r i (z)dz , t ≥ 0, i = 1, 2. (5)
3. In order to monitor the system health, a periodic inspection scheme with period T is planned at a cost rate C i > 0 to detect the degradation level. At each inspection, if the system is still running, we estimate the MRL from the detected degradation level, and we make a preventive maintenance decision based on this indicator. If the value of MRL is less than a threshold m, the system is preventively replaced with a cost C p > C i ; otherwise, nothing is done. But if the system is detected to be in a failure state (i.e., due to degradation or shock), a corrective replacement is performed with a cost C c > C p . In this case, because of the system inactivity after failure, an additional cost is incurred from the failure time until the replacement time at a cost rate C d > 0. We assume that all the maintenance operations are perfect and take negligible time.

According to the above descriptions, constructing the MRL and choosing the optimal values of T and m are very important deciding the performance of the maintenance strategy for the considered system. Two next sections deal with the mathematical models to construct and optimize these quantities.

MEAN RESIDUAL LIFE FUNCTION

Generally, the mean residual life function at time t ≥ 0 can be defined as [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF] E

[σ f -t | σ f > t, Z (t)] ,
where σ f denotes the time to failure and Z (t) denotes the covariates vector related to the system condition obtained at time t. In the model presented in this paper, we assume that the data collected in the inspection times can perfectly reflect the degradation level of the system at time t, Z (t) = X (t). Thus, given the degradation level at time t X (t) = z and known the system is still working at time t, σ f > t, the MRL function is given by

m z (t) = E [σ f -t | σ f > t, X (t) = z] = ∞ t R (u | X (t) = z) du, (6)
where R (u | X (t) = z) denotes the conditional reliability of the system at time u, u ≥ t, given

X (t) = z.
For the considered system, where the system fails whenever degradation level exceeds a fixed threshold

L or a shock occurs, R (u | X (t) = z) is hence computed by: R (u | X (t) = z) = P (X (u) < L, N s (u) = 0 | X (t) = z) = R 1 1 {z≤Ms} + R 2 1 {z>Ms} , (7)
where 1 denotes the indicator function which equals 1 if the argument is true and 0 otherwise, R 1 and R 2 are given respectively by:

R 1 = Fσ Ms-z (u -t) F1 (u) F1 (t) + u t Fσ L -σ Ms (u -v) F1 (v) F2 (u) F1 (t) F2 (v) f σ Ms-z (v -t) dv, R 2 = Fσ L-z (u -t) F2 (u) F2 (t) ,
where Fσ Ms-z , Fσ L-z , f σ Ms-z , F1 and F2 is computed from equations 2, 3 and 5. Fσ L -σ Ms denotes the survival function of the random variable σ Lσ Ms given by [START_REF] Bérenguer | Maintenance policy for a continuously monitored deteriorating system[END_REF] Fσ L -σ

Ms (t) = F(σ Ms ,σ L ) (t, x, y)1 {y>x>0} + Fd (t, x, y)1 {y=x>0} , (8)
where Consequently, the MRL function of the considered system is computed by equations 6 and 7.

F(σ Ms ,σ L ) (t, x, y) = - Ms<x<L,0<x+y<L,0<y ∞ 0 f αu,β (x)du ∂f αt,β (
Analyzing the monotonicity of the MRL function is important to compute the probabilities of the different maintenance actions. Lemma 1 shows the monotonicity of the MRL function. In order to fulfill the space requirements of the conference, we omit its proof in the present paper; interested readers can refer to [START_REF] Huynh | Quantification de l'apport de l'information de surveillance dans la prise de décision en maintenance[END_REF].

Lemma 1 For a fixed degradation level z, 0 ≤ z ≤ L, the MRL function, m z (t), is a non-decreasing function of the time t. For a fixed time t, 0 ≤ t, the MRL function, m z (t), is a non-decreasing function of the degradation level z.

To illustrate the property of monotonicity of the MRL function, Figure 1 shows the MRL surface as a function of t and z when the shocks intensities r 1 (t) and r 2 (t) are linear in time t. the MRL function, m z (t), depends on both z and t, so it contains more information on the system compared to the degradation level. This is why the MRL function is more significant for maintenance decision-making than degradation level, and it is then interesting to estimate and integrate the MRL in maintenance planning.

MAINTENANCE COST MODEL

In this section, the maintenance cost model of the considered strategy based on the MRL is analyzed. By applying equation 1, the long-run expected maintenance cost rate is computed by the following formula:

C ∞ (T, m) = C p P p (T, m) + C c P c (T, m) + C i E [N i (T, m)] + C d E [W d (T, m)] E [S (T, m)] , (9)
where E [S (T, m)] is the expected length of a replacement cycle under the considered policy, P p (T, m), P c (T, m), E [N i (T, m)] and E [W d (T, m)] are, respectively, the probability of a preventive replacement, the probability of a corrective replacement, the expected number of inspections and the expected down time of the system in a replacement cycle. Computing these quantities is complicated but quite classical. The idea is, firstly, to represent all the possible events of MRL related to each quantity, then, to transform these events into the ones of degradation levels thanks to the monotony property of the MRL function, and finally, to compute the associated probabilities. In this way, we must define a degradation level corresponding to a given MRL value. For a fixed working time t ≥ 0, let z t,m be the minimum value of degradation level z such that the mean residual life is less than m ≥ 0, that is,

z t,m = inf {z ≥ 0, m z (t) ≤ m} (10)
(if z t,m > L, then we impose z t,m ≡ L). Furthermore, using the monotonicity of the MRL function, for t 1 ≤ t 2 , m z (t 1 ) ≥ m z (t 2 ); hence z t 1 ,m ≥ z t 2 ,m for t 1 ≤ t 2 . The remainder of this section reports the detailed mathematical results of the quantities of interest.

Preventive maintenance probability

Under the considered MRL-based maintenance policy, a preventive maintenance is performed at inspection time if the system is still working and the MRL calculated from the detected degradation level of the system is less than the value m. As such the probability of the preventive replacement in a renewal cycle is calculated by

P p (T, m) = ∞ k=0 P p ((k + 1) T ) , (11)
where, P p ((k + 1) T ) is the probability of a preventive replacement at time (k + 1) T , k = 0, 1, 2, . . ., given by:

P p ((k + 1) T ) = P m z kT (kT ) > m ≥ m z (k+1)T ((k + 1) T ) , (k + 1) T < σ L , σ s > (k + 1) T = P X(kT ) < z kT,m , z (k+1)T,m ≤ X((k + 1) T ), (k + 1) T < σ L , σ s > (k + 1) T = P kT < σ z kT,m , σ z (k+1)T,m ≤ (k + 1) T < σ L , σ s > (k + 1) T ,
where σ zt,m is the hitting time of the degradation level z t,m and σ s is the time to the first shock.

To explicite the formula of P p ((k + 1) T ), we consider three exclusive events: M s ≤ z (k+1)T,m , z (k+1)T,m < M s ≤ z kT,m and {z kT,m < M s }. Before going ahead with the calculation, we define, for 0

≤ u ≤ v ≤ w and 0 ≤ z 1 ≤ z 2 ≤ z 3 (i.e., 0 = σ 0 ≤ σ z 1 ≤ σ z 2 ≤ σ z 3 ), the functions G (u, v) and H σz 1 ,σz 2 ,σz 3 (u, v, w) by G (u, v) = P (N s (v) = 0 | σ Ms = u) = F1 (u) F2 (v) F2 (u) , (12)
H σz 1 ,σz 2 ,σz 3 (u, v, w) = f σz 2 -σz 1 (v -u) Fσz 3 -σz 2 (w -v) , (13)
where F1 and F2 is given from equation 5, Fσz 3 -σz 2 and f σz 2 -σz 1 is computed and derived respectively from equation 8.

When M s ≤ z (k+1)T,m , the probability P p ((k + 1) ∆T ) is computed as

P p,1 = kT 0 kT u G (u, (k + 1) T ) H σ Ms ,σz (k+1)T,m ,σz kT,m (u, v, (k + 1) T ) f σ Ms (u) dvdu + kT 0 (k+1)T kT G (u, (k + 1) T ) H σ Ms ,σz (k+1)T,m ,σ L (u, v, (k + 1) T ) f σ Ms (u) dvdu + (k+1)T kT (k+1)T u G (u, (k + 1) T ) H σ Ms ,σz (k+1)T,m ,σ L (u, v, (k + 1) T ) f σ Ms (u) dvdu + kT 0 kT u (k+1)T kT -G (u, (k + 1) T ) ∂H σ 0 ,σ Ms ,σz (k+1)∆T,m (0, u, v) ∂v × H σz (k+1)T,m ,σz kT,m ,σ L (v, w, (k + 1) T ) dwdvdu. ( 14
)
When z (k+1)T,m < M s ≤ z kT,m , the probability P p ((k + 1) T ) is rewritten as

P p,2 = kT 0 kT u G (v, (k + 1) T ) H σz (k+1)T,m ,σ Ms ,σz kT,m (u, v, (k + 1) T ) f σz (k+1)T,m (u) dvdu + kT 0 (k+1)T kT G (v, (k + 1) T ) H σz (k+1)T,m ,σ Ms ,σ L (u, v, (k + 1) T ) f σz (k+1)T,m (u) dvdu + F1 ((k + 1) T ) (k+1)T 0 H σ 0 ,σz (k+1)T,m ,σ Ms (0, u, (k + 1) T ) du + (k+1)T kT (k+1)T u G (v, (k + 1) T ) H σz (k+1)T,m ,σ Ms ,σ L (u, v, (k + 1) T ) f σz (k+1)T,m (u) dvdu + kT 0 kT u (k+1)T kT -G (v, (k + 1) T ) ∂H σ 0 ,σz (k+1)T,m ,σ Ms (0, u, v) ∂v × H σ Ms ,σz kT,m ,σ L (v, w, (k + 1) T ) dwdvdu. ( 15
)
And when z kT,m < M s , the probability P p ((k + 1) T ) is given by

P p,3 = F1 ((k + 1) T ) kT 0 (k+1)T kT H σz (k+1)T,m ,σz kT,m ,σ Ms (u, v, (k + 1) T ) f σz (k+1)T,m (u) dudv + F1 ((k + 1) T ) kT 0 H σ 0 ,σz (k+1)T,m ,σz kT,m (0, u, (k + 1) T ) du + F1 ((k + 1) T ) (k+1)T kT H 0,σz (k+1)T,m ,σ Ms (0, u, (k + 1) T ) du + (k+1)T kT (k+1)T u G (v, (k + 1) T ) H σz (k+1)T,m ,σ Ms ,σ L (u, v, (k + 1) T ) f σz (k+1)T,m (u) dudv + kT 0 (k+1)T kT (k+1)
T v -G (w, (k + 1) T ) ∂H σ 0 ,σz (k+1)T,m ,σz kT,m (0, u, v) ∂v ×
H σz kT,m ,σ Ms ,σ L (v, w, (k + 1) T ) dwdvdu.

Consequently, the probability P p ((k + 1) T ) is given by: P p ((k + 1) T ) = P p,1 1 {Ms≤z(k+1)T,m} + P p,2 1 {z(k+1)T,m<Ms≤zkT,m} + P p,3 1 {zkT,m<Ms} .

where P p,1 , P p,2 and P p,3 are given by 14, 15 and 16 respectively.

Corrective maintenance probability

Similarly to the above analysis, the probability of a corrective replacement under the considered MRL-based maintenance policy in a renewal cycle is calculated by

P c (T, m) = ∞ k=0 P c ((k + 1) T ) , (18)
where, P c ((k + 1) T ) is the probability of a corrective replacement at time (k + 1) T , k = 0, 1, 2, . . ., given by:

P c ((k + 1) T ) = P kT < σ z kT,m ≤ σ L ≤ (k + 1) T, σ L < σ s + P kT < σ s ≤ (k + 1) T, σ s < σ z kT,m + P kT < σ z kT,m ≤ σ s ≤ (k + 1) T, σ s < σ L .
When M s ≤ z kT,m , the probability P c ((k + 1) T ) is given by

P c,1 = kT 0 f σ Ms (u) (k+1)T kT (k+1)T v - ∂ ∂w H σ Ms ,σz kT,m ,σ L (u, v, w) G (u, w) dwdvdu + (k+1)T kT (k+1)T u - ∂ ∂v Fσ L -σ Ms (v -u) G (u, v) f σ Ms (u) dvdu + kT 0 (k+1)T kT - ∂ ∂v (G (u, v)) H σ 0 ,σ Ms ,σz kT,m (0, u, v) dvdu + (k+1)T kT f 1 (u) Fσ Ms (u) du. (19)
When z kT,m < M s , the probability P c ((k + 1) T ) is computed by

P c,2 = (k+1)T kT (k+1)T u (k+1)T v - ∂ ∂w H σz kT,m ,σ Ms ,σ L (u, v, w) G (u, w) f σz kT,m (u) dwdvdu + (k+1)T kT (k+1)T u H σ 0 ,σz kT,m ,σ Ms (0, u, v) f 1 (v) dvdu + (k+1)T kT f 1 (u) Fσz kT,m (u) du. (20)
Hence, the probability P c ((k + 1) T ) is given by

P c ((k + 1) T ) = P c,1 1 {Ms≤zkT,m} + P c,2 1 {Ms>zkT,m} , (21)
where P c,1 and P c,2 are given by 19 and 20 respectively.

Expected downtime in a replacement cycle

The system downtime corresponds to the time interval from a failure to the next replacement. Assuming (k + 1) T is the failure time, the downtime of the considered system is expressed by

W d (T, m) = ∞ k=0 ((k + 1) T -σ L ) 1 kT ≤σz kT,m ,kT ≤σ L ≤(k+1)T,σ L <σs + ∞ k=0 ((k + 1) T -σ s ) 1 kT ≤σz kT,m ,kT ≤σs≤(k+1)T,σs<σ L . ( 22
)
As such, the expected downtime is given by

E [W d (T, m)] = ∞ k=0 d 1 (k, T ) 1 {Ms<zkT,m} + ∞ k=0 d 2 (k, T ) 1 {Ms≥zkT,m} , (23)
where d 1 (k, T ) is given by

d 1 (k, T ) = kT 0 (k+1)T kT (k+1)T v -((k + 1) T -w) ∂ ∂w H σ Ms ,σz kT,m ,σ L (u, v, w) G (u, w) dwdvdu + (k+1)T kT (k+1)T u -((k + 1) T -v) ∂ ∂v Fσ L -σ Ms (v -u) G (u, v) f σ Ms (u) dvdu + kT 0 (k+1)T kT -((k + 1) T -v) ∂ ∂v (G(u, v)) H σ 0 ,σ Ms ,σz kT,m (0, u, v) dvdu + (k+1)T kT ((k + 1) T -u) f 1 (u) Fσ Ms (u) du, (24)
and d 2 (k, T ) is given by

d 2 (k, T ) = (k+1)T kT (k+1)T u (k+1)T v -((k + 1) T -w) ∂ ∂w H σz kT,m ,σ Ms ,σ L (u, v, w) G (v, w) × f σz kT,m (u)dwdvdu + (k+1)T kT (k+1)T u ((k + 1) T -v) H σ 0 ,σz kT,m ,σ Ms (0, u, v) f 1 (v) dvdu + (k+1)T kT ((k + 1) T -u) f 1 (u) Fσz kT,c (u) du. (25)

Expected replacement time and expected number of inspections

Under the considered MRL-based maintenance policy, the system can be replaced either preventively or correctively. So the expected time to a system replacement is given by

E [S (T, m)] = ∞ k=0 (k + 1) T (P p ((k + 1) T ) + P c ((k + 1) T )) , (26)
where P p ((k + 1) T ) and P c ((k + 1) T ) are given by 17 and 21 respectively. And the expected number of inspections in a replacement cycle is given by

E [N i (T, m)] = E [S (T, m)] T . (27)
where E [S (T, m)] is given by 26.

MAINTENANCE OPTIMIZATION AND NUMERICAL EXAMPLE

As indicated in the previous sections, the inspection period T and the preventive replacement threshold m decide the performance of the considered maintenance strategy. Thus, optimizing the strategy is reduced to find the optimal values T opt and m opt that minimize the long-run expected maintenance cost rate C ∞ (T, m) given by equation 9, that is

C ∞ (T opt , m opt ) = min T,m {C ∞ (T, m) , T > 0, 0 ≤ m ≤ M T T F } , (28)
where M T T F denotes the mean time to failure of system, which is computed by

M T T F = ∞ 0 F1 (u) Fσ Ms (u) du + ∞ 0 ∞ 0 Fσ L -σ Ms (v) F1 (u) F2 (v + u) F2 (u) f σ Ms (u) dvdu. (29)
Due to the complexity of the cost function C ∞ (T, m), obtaining the analytical solution of T opt and m opt is almost impossible; numerical methods (e.g., policy iteration algorithm [START_REF] Tijms | A first course in stochastic models[END_REF]) are used to overcome the problem. To illustrate the cost surface and the optimum solution of the considered MRL-based maintenance strategy, we use the following data set: α = 0.5, β = 0.5, M s = 15, L = 30, r 1 (t) = 0.0025t + 0.01, r 2 (t) = 0.0025t + 0.1, C i = 5, C p = 50, C c = 100 and C d = 25. The numerical example of the strategy is shown in Figure 2. Sub-figures 2a and 2b illustrate respectively the shape and the iso-level cost curves of the expected maintenance cost rate when the inspection period T and the preventive replacement threshold m vary. The convexity of the cost surface shows the existence of an optimal setting of parameters T and m. In fact, for the chosen data set, optimal values are T opt = 4.5 and m opt = 8.5, which correspond to the optimal cost rate given by C ∞ (4.5, 8.5) = 7.8416.

CONCLUSION

The present paper introduces the cost model of a predictive maintenance strategy for a system subject to the computing failure modes due to degradation and shocks. The system is described by an extension of DTS model where the shock intensity is assumed to be dependent on the system degradation level. We have shown that, under the assumptions of the model, the MRL function of the system is a decreasing function of the degradation level and working time. The preventive maintenance decisions of the proposed strategy is based on the MRL function calculated at different inspection times. Due to the monotonicity of the MRL function, the degradation levels at different inspection times for the same value of MRL are different. Hence, the considered MRL-based maintenance strategy is equivalent to a degradation-based maintenance strategy with multiple preventive degradation thresholds. Since these thresholds can finely adapt to the health state of system, this strategy is better than the one based on single degradation threshold. Furthermore, for our new strategy, multiple preventive thresholds are integrated in an unique MRL indicator, so it is generally more simple and more complete than degradation-based maintenance strategies.

Figure 1 Figure 1 :

 11 Figure 1: MRL functionα = 0.1, β = 0.1, L = 30, M s = 15, r 1 (t) = 0.0025t + 0.01 and r 2 (t) = 0.0025t + 0.1.

  Mesh of C∞ (T, m)m T C ∞ (T opt ,m opt ) = 7.8416, m opt = 8.5, T opt = 4Iso-level of C∞ (∆T, m)

Figure 2 :

 2 Figure 2: Long-run expected maintenance cost rate C ∞ (T, m)

ACKNOWLEDGEMENTS

For I.T. Castro, this work was supported by the Ministerio de Educación y Ciencia (Spain) under grant MTM2009-07634


Sebastien Beyou

Anne Cuzol

Sai Subrahmanyam Gorthi

Etienne Mémin
email: etienne.memin@inria.fr

Weighted Ensemble Transform Kalman Filter for Image Assimilation

This paper proposes an extension of the Weighted Ensemble Kalman filter (WEnKF) proposed by Papadakis et al. (2010) for the assimilation of image observations. The main contribution of this paper consists in a novel formulation of the Weighted filter with the Ensemble Transform Kalman filter (WETKF) incorporating directly as a measurement model a nonlinear image reconstruction criterion. This technique has been compared to the original WEnKF on numerical and real world data of 2D turbulence observed through the transport of a passive scalar. It has been in particular applied for the reconstruction of oceanic surface current vorticity fields from Sea Surface Temperature satellite data. This latter technique enables a consistent recovery of oceanic surface currents, vorticity maps along time in presence of large missing data areas and strong noise.

Introduction

The analysis of geophysical fluid flows is of crucial interest in Oceanography, Meteorology, or Hydrology for forecasting applications, or for the monitoring of hazards. In all those domains a great number of orbital or geostationary satellites provides a huge amount of image data with still increasing spatial and temporal resolutions. The combination of such fine scale observations with larger scales dynamical models in data assimilation schemes might be of great interest to infer with a greater accuracy the evolution of the observed geophysical fluid flows. Compared to in situ measurement techniques supplied by dedicated probes or Lagrangian drifters, satellite images give access to a much denser observation field but provide unfortunately only an indirect access to the physical quantities of interest. The direct estimation from a satellite image sequence of the flow state variables leads hence to difficult inverse problems. This is particularly true for the recovery of fluid flow motion fields from successive pairs of images. Such a direct estimation is intrinsically ill-posed and requires the introduction of additional smoothing constraints imposed either globally on the image domain in variational approaches (see [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF] for a review) or within a local window in correlation based techniques [START_REF] Adrian | Particle imaging techniques for experimental fluid mechanics[END_REF][START_REF] Emery | An objective method for computing advective surface velocities from sequential infrared satellite images[END_REF][START_REF] Ottenbacher | Low-level cloud motion winds from Meteosat high-resolution visible imagery[END_REF][START_REF] Szantai | A method for the construction of cloud trajectories from series of satellite images[END_REF]. Those techniques rely heavily on this implicit or explicit smoothing assumption and are prone to strong errors in presence of noise, missing data or regions of low photometric variability. Consequently, there is a strong difficulty in imposing a trajectorial consistency to the sequence of (generally independently) estimated motion fields. In the best scenario case these methods provide accurate instantaneous displacements only for a limited range of scales and experience severe difficulties for medium to small scales measurements due to the smoothing prior used. Their accuracy is also restricted to regions that exhibit sufficient luminance variations. Although these techniques are very useful to measure in a fast and tireless way motion fields from an image sequence, all these difficulties generate in large areas inconsistent measurements along time in real world atmospheric or oceanic satellite images.

One way to enforce such a dynamical consistency of the velocity measurements consists in embedding the estimation problem within an image based assimilation process. Variational assimilation of image information has been recently considered by several authors for the tracking along time of convective cells [START_REF] Thomas | Data assimilation for convective cells tracking on meteorological image sequences[END_REF] or for the estimation of consistent motion fields [START_REF] Corpetti | Pressure image asimilation for atmospheric motion estimation[END_REF][START_REF] Titaud | Assimilation of image sequences in numerical models[END_REF][START_REF] Papadakis | An optimal control technique for fluid motion estimation[END_REF]. All these techniques have been assessed on simplified dynamics: level set transport equation associated to a transportation component defined from motion measurements with Brownian uncertainties for [START_REF] Thomas | Data assimilation for convective cells tracking on meteorological image sequences[END_REF], different forms of the Shallow-water equations for [START_REF] Corpetti | Pressure image asimilation for atmospheric motion estimation[END_REF] and [START_REF] Titaud | Assimilation of image sequences in numerical models[END_REF] or 2D velocity-vorticity Navier-Stokes formulation in [START_REF] Papadakis | An optimal control technique for fluid motion estimation[END_REF]. The techniques differ mainly in the way they are dealing with the relation linking the measurements and the state variables (the observation operator). In [START_REF] Corpetti | Pressure image asimilation for atmospheric motion estimation[END_REF]; [START_REF] Thomas | Data assimilation for convective cells tracking on meteorological image sequences[END_REF] linear observation operators derived from the normal equations associated to a local least square motion estimation formulation are considered. The problem is addressed in a slightly different way in Titaud c 0000 Tellus, 000, 000-000 et al. (2010), simulating from an initial image a sequence transported by the state variable dynamics. The observation operator is then defined as a relation between this simulated sequence and the observed image sequence. This approach relies on the hypothesis that the actual image sequence is close to a reference sequence transported by the model.

Let us note that all these techniques may incorporate either preconditioning stages in their internal iterative optimization or covariance specifications that include smoothing functional similar as those employed for the motion estimation problems [START_REF] Souopgui | Vector field regularization by generalized diffusion[END_REF].

Another family of data assimilation procedures is defined as stochastic filters implemented through Monte Carlo samples called ensemble members or particles whether the filtering is implemented through an ensemble Kalman filter or with a particle filter. Contrary to the previous deterministic techniques, those assimilation schemes aim at recovering a description of the probability density function of the state variable trajectory based on the complete history of measurements up to the current instant. These schemes have the advantage to deal with the whole nonlinearity of the dynamics and to include directly a representation of the uncertainty associated to the considered dynamics. Obviously the building of such a stochastic model is debatable on the precise form it may take. However, this ability is very interesting in cases for which the dynamical model is only imprecisely known. Let us note that second order variational methods allows also computing the error covariance matrix through second order adjoint systems (Le [START_REF] Dimet | Second-order information in data assimilation[END_REF]. This is however at the price of a significant computational cost increase.

Recently a data assimilation procedure embedding an Ensemble Kalman filter into the particle filter framework, entitled Weighted Ensemble Kalman filter, has been proposed [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF]. This filter combines all the nice properties of the ensemble Kalman mechanics with the correction scheme of the particle filter. Conceptually, this filter is more adapted to a non Gaussian distribution of the ensemble members, which is a situation that may arise at very short time horizon with nonlinear stochastic dynamics 1 . It has been also

shown experimentally that this modified filter exhibits a faster convergence. This has been however assessed only on synthetic 2D turbulence models in the simple case of measurements artificially built from a noisy version of the true simulated data [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF].

The main contribution of the paper consists thus first in extending the study of WEnKF 1 Other solutions have been proposed to deal with non-Gaussianity in the framework of variational methods or stochastic filters (interested readers may refer to the review of [START_REF] Bocquet | Beyond gaussian statistical modeling in geophysical data assimilation[END_REF] or [START_REF] Yang | Handling nonlinearity and non-Gaussianity in ensemble Kalman filter[END_REF] c 0000 Tellus, 000, 000-000 to the assimilation of observations supplied by image data. We will mainly focus on the incorporation of a direct image reconstruction error measurement. Here the great advantage is to avoid the use of external techniques to supply so called pseudo observations at the price of a nonlinear observation operator. The second purpose of the paper lies in extending the weighted ensemble framework based on randomly sampling of the measurement noise (referred usually as ensemble schemes with perturbed observation [START_REF] Evensen | Sequential data assimilation with a non linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF][START_REF] Houtekamer | Data assimilation using an ensemble Kalman filter technique[END_REF] to a weighted version that uses the ensemble transform Kalman filter (ETKF), proposed by [START_REF] Bishop | Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects[END_REF] and developed by [START_REF] Tippett | Ensemble square root filters[END_REF] and [START_REF] Wang | Which is better, an ensemble of positive-negative pairs or a centered simplex ensemble?[END_REF]. This transform is built from a (mean-preserving) square-root matrix to update the mean and covariances, in accordance with the Kalman update equations. This extension was suggested in [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF] but was not implemented. We will see that this transform is well adapted to the case of nonlinear direct image measurements and enables processing the fluid flow image sequences directly, without the need of any intermediate velocity fields observations. The proposed WETKF makes use of the ETKF mechanism as the proposal distribution of the particle filter and corrects the different members according to the particle filtering selection step. This selection scheme associates to each ensemble member a filtering weight determined in our case through a mean square error of displaced frame differences, and performs a resampling favoring particles corresponding to high values of the filtering probability distribution.

We organize the paper as follows: Section 2 recalls the principles governing the construction of the Weighted Ensemble Kalman filter proposed in [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF]. After a brief discussion of ensemble transform Kalman filter in Section 2.2, we describe its particle filter extension. Its application to fluid flow analysis through the direct observation of an image sequence depicting the transportation of a passive scalar is presented in Section 3. Section 4 presents the experimental validation of these approaches and brings out some elements of comparison. Real results on 2D sequences of 2D turbulence and satellite oceanic data are presented. Finally, Section 5 concludes the paper with a discussion to future directions.

Assimilation with the weighted ensemble transform Kalman filter

We describe in this section how the assimilation problem can be formulated sequentially by a filtering approach adapted to a nonlinear and high-dimensional framework.

c 0000 Tellus, 000, 000-000

Filtering in a nonlinear and high-dimensional setting

Stochastic filters aim at estimating the posterior probability distribution p(x 0:k |y 1:k ) of a state variable trajectory x 0:k starting from an initial state x 0 up to the state x k ∈ R n at current time k. The state variable trajectory is obtained through the integration of a dynamical system:

dx t = M(x t )dt + ηdB t , (1)
where M denotes a deterministic linear/nonlinear dynamical operator, corresponding to a physical conservation law describing the state evolution. The random part dB t accounts for the uncertainties in the deterministic state model. It is usually a Gaussian noise, white in time but correlated in space with covariance Q = ηη T . It is assumed that the true state is unknown and observed through y k ∈ R m at discrete time instants. These observations are assumed to be linked to the state variable through the following measurement equation:

y k = H(x k ) + γ k , (2)
where the observation noise γ k is a white Gaussian noise with covariance matrix R and H stands for the linear/nonlinear mapping from the state variable space to the observation space. We note that the (integration) time step used for the state variable dynamics δt is usually much smaller (about 10-100 times), than the time interval δk between two subsequent measurements. A sequence of measurements or observations from time 1 to k will be denoted by a set of vectors of dimension m as: y 1:k = {y i , i = 1, . . . , k}, where the time between two successive measurements is arbitrarily set to δk = 1.

In a nonlinear context, it is known that the filtering equations can no longer be solved by the Kalman algorithm and its direct variants. Ensemble Kalman filters [START_REF] Anderson | A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts[END_REF][START_REF] Bishop | Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects[END_REF][START_REF] Evensen | Sequential data assimilation with a non linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF][START_REF] Houtekamer | Data assimilation using an ensemble Kalman filter technique[END_REF][START_REF] Ott | A local ensemble Kalman filter for atmospheric data assimilation[END_REF] have been developed to tackle the filtering problem in nonlinear and highdimensional systems. However, these methods still rely on a Gaussian assumption. On the other hand, particle filters [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF][START_REF] Del Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF][START_REF] Gordon | Novel approach to non-linear/non-gaussian bayesian state estimation[END_REF] are able to solve exactly the filtering equations (up to the Monte Carlo approximation). But, in their simplest formulation, they are not suitable to high-dimensional systems [START_REF] Van Leeuwen | Particle filtering in geophysical systems[END_REF][START_REF] Snyder | Obstacles to high-dimensional particle filtering[END_REF].

Particle filtering techniques implement an approximation of the state posterior density c 0000 Tellus, 000, 000-000 p(x 0:k |y 1:k ) (called filtering law), using a sum of N weighted Diracs:

p(x 0:k |y 1:k ) ≈ N i=1 w (i) k δ x (i) 0:k (x 0:k ) , (3)
centered on hypothesized elements (called particles) of the state space sampled from a proposal distribution π(x 0:k |y 1:k ). This distribution, called the importance distribution, approximates the true filtering distribution. Each sample is then weighted by a weight, w

k , accounting for the ratio between the two distributions. Any importance function can be chosen (with the only restriction that its support contains the filtering distribution one).

Under weak hypotheses the importance ratio can be recursively defined as:

w (i) k ∝ w (i) k-1 p(y k |x (i) k )p(x (i) k |x (i) k-1 ) π(x (i) k |x (i) 0:k-1 , y 1:k ) . ( 4
)
By propagating the particles {x

(i) k , i = 1, . . . , N } from time k-1 through the proposal density π(x k |x (i) 0:k-1 , y 1:k ),
and by weighting the sampled states with w (i) k , a recursive sampling of the filtering law is obtained. The standard choice when applying the particle filter is to fix the proposal distribution to the transition law p(x k |x k-1 ). In that case, it is easy to sample from the proposal (this requires only to sample from the dynamical model ( 1)). However, this choice has the major drawback of sampling the particles without taking into account the observation y k . This makes the particle filter not suitable in a high-dimensional context. As a matter of fact, particles are led to explore the state space blindly following the dynamical model only and are most of the time far away from the observation at correction time, yielding hence very quickly to a filter divergence. It is thus necessary to introduce the observation within the proposal distribution in order to guide the particles, and then make the particle filtering technique suitable to high-dimensional problems [START_REF] Van Leeuwen | Particle filtering in geophysical systems[END_REF]. Such a strategy has been used in [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF] through the ensemble Kalman filter mechanism. Relying on the usual assumption of the ensemble techniques (i.e. considering the dynamics as a discrete Gaussian system), the proposal can be set to a Gaussian conditional

distribution π(x k |x (i) k-1 , y k ).
In order to make the estimation of the filtering distribution exact (up to the sampling), each member of the ensemble must be weighted at each instant k, with appropriate weights w (i) k defined from (4). With a systematic resampling scheme and for high dimensional systems represented on the basis of a very small number of ensemble members it can be shown [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF] that the weights reduces to the computation of the data likelihood:

w (i) k ∝ p(y k |x (i) k ), (5)
c 0000 Tellus, 000, 000-000 k and provide for the filtering objective a natural quality ranking of a given sample. In the following we describe how this ensemble Kalman proposal step can be performed with a particular scheme called the Ensemble Transform Kalman Filter.

Ensemble transform Kalman filter

Several forms of ensemble Kalman filters have been successfully proposed since the seminal work of Evensen [START_REF] Evensen | Sequential data assimilation with a non linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF]. Those filters which correspond to the ensemble or particle implementation of the Kalman Filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] have demonstrated their performances in geophysical sciences as they allow propagating very high-dimensional state vectors. This advantage springs from the ability of the EnKF to approximate or compute the forecast and analysis two first moments (needed in KF recursion) directly from the ensemble, while propagating the ensemble members through the exact dynamics and observation models, without an explicit linearization -as it is done in Extended Kalman filter [START_REF] Anderson | Optimal Filtering[END_REF]. The EnKF avoids explicit computation of the covariance matrix, which may not even be affordable in very high-dimensional states, and further the algorithms of [START_REF] Evensen | The ensemble Kalman filter, theoretical formulation and practical implementation[END_REF] or [START_REF] Bishop | Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects[END_REF] present computationally efficient schemes for the inversion of the covariance matrix required in the Kalman gain computation, through a low dimensional singular value decomposition, when the number of ensemble members is much lower than the state dimension, which is usually the case. Two kinds of ensemble filters have been proposed so far in the literature. The first kind of method [START_REF] Evensen | The ensemble Kalman filter, theoretical formulation and practical implementation[END_REF] introduces artificially perturbed observations in order to enable a complete Monte Carlo estimation of the mean and covariance of the filtering distribution.

This technique has been reported to yield good Gaussian approximation in a nonlinear regime but may show some bias for a too small ensemble. The second kind of techniques relies on so-called deterministic square-root filters. Among the different variations of squareroot filters, the Ensemble Transform Kalman filters, proposed by [START_REF] Bishop | Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects[END_REF] and developed by [START_REF] Tippett | Ensemble square root filters[END_REF] and [START_REF] Wang | Which is better, an ensemble of positive-negative pairs or a centered simplex ensemble?[END_REF], employ a (mean-preserving) square-root matrix to update the mean and covariances, in accordance with the Kalman update equations. We present now the essential ideas of this scheme.

Considering first a linear observation operator H, the covariance and mean update equac 0000 Tellus, 000, 000-000 tions in the Kalman filter at instant k are respectively given by:

P a k = P f k -P f k H T (HP f k H T + R) -1 HP f k , (6)
xa k = xf k + P a k H T R -1 (y k -Hx f k ), (7)
where superscript "f " stands for the forecast step, whereas "a" corresponds to the correction step (also called analysis) and the overline symbol denotes the empirical average. Using a forecasted ensemble members perturbation centered around the mean,

X f k = {x f,(i) k - xf k , i = 1, .
. . , N }, an empirical approximation of the forecast covariance is obtained as:

P f k = (N -1) -1 X f k X f k
T , and we get:

(N -1)P a k = X f k X f k T -X f k X f k T H T (HX f k X f k T H T + R) -1 HX f k X f k T , with R = (N -1)R.
Let us note that in general the ensemble dimension N is much lower than the state space dimension (N << n). A nonlinear observation operator H(X f k ), can be incorporated by replacing HX f k with a centered matrix H(X f k ), defined as:

H(X f k ) = {H(x f,(i) k ) -H(x f k ), i = 1, . . . , N }. (8)
The same expression of the nonlinear observation is used in [START_REF] Houtekamer | A sequential ensemble Kalman filter for atmospheric data assimilation[END_REF]; [START_REF] Hunt | Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter[END_REF]. This leads to the following approximate covariance update with a nonlinear observation model:

(N -1)P a k = X f k X f k T -X f k H(X f k ) T (H(X f k )H(X f k ) T + R) -1 H(X f k )X f k T ,
which can be written under the form (N -1)

P a k = X f k D k X f k T with: D k = I -H(X f k ) T (H(X f k )H(X f k ) T + R -1 H(X f k ), (9)
and where I is the

N × N identity matrix. Writing (N -1)P a k = X a k X a k
T where X a k gathers desired centered corrected ensemble members, and expressing

D k = A k A T k through its ma- trix square root, we get X a k X a k T = X f k A k A T k X f k T , or X a k = X f k A k .
In order to compute the matrix square-root A k , we rely on the Sherman-Morison-Woodbury formula as suggested by [START_REF] Tippett | Ensemble square root filters[END_REF] and which leads to:

D k = (I + H(X f k ) T R-1 H(X f k )) -1 . ( 10
)
From the eigenvalue decomposition:

I + H(X f k ) T R-1 H(X f k ) = U k (I + Λ k )U T k , (11)
where U k is an orthogonal matrix the square root A k can be written as:

A k = U k (I + Λ k ) -1/2 . ( 12
)
Note however that this square-root is not unique and can be replaced by A k V, where V is c 0000 Tellus, 000, 000-000 any arbitrary orthonormal matrix (VV T = I) of suitable dimension. In particular, one can look for a mean preserving square-root such that:

X a k 1 = X f k A k 1 = 0, ( 13
)
where 1 is a vector of ones. As the forecast ensemble is centered by definition, this condition is fulfilled for an ensemble matrix X f k of rank N -1 if and only if A k 1 = α1. It can be shown [START_REF] Ott | A local ensemble Kalman filter for atmospheric data assimilation[END_REF][START_REF] Sakov | Implications of the form of the ensemble transformation in the ensemble square root filter[END_REF][START_REF] Wang | Which is better, an ensemble of positive-negative pairs or a centered simplex ensemble?[END_REF] , that taking:

A k = U k (I + Λ k ) -1/2 U T k . ( 14
)
satisfies the above mean preserving condition for the transform. The improvement brought by the use of this mean preserving transform over non mean preserving square root filters or ensemble filters with perturbed observations has been experimentally shown and analyzed in [START_REF] Sakov | Implications of the form of the ensemble transformation in the ensemble square root filter[END_REF].

As a final step a mean analysis term xa k has to be added to each element x a,(i) k

. The mean analysis vector reads:

xa k = xf k + (N -1) -1 X f k D k H(X f k ) T R -1 (y k -H(x f k )). (15)
Note that computing (10) only requires the inversion of a N × N matrix, if the inverse of the measurement error covariance matrix R is available. Usually R is defined as a diagonal or band diagonal (very sparse) matrix and is easy to invert. The analysis step can thus be accomplished efficiently. In our case R will be assumed to be a diagonal matrix and is thus trivially inverted. Let us note that in the case of image data, it would be very interresting to investigate the use of sparse banded information matrices (R -1 ) associated to statistical graphical models [START_REF] Lauritzen | Graphical models[END_REF].

Ensemble filters may be also accompanied with so-called localization procedures introduced to discard spurious long range correlations that appear in empirical covariance matrices built from a small number of samples. Such a localization procedure is carried out either explicitely by considering only the observation within a window centered around a given point grid, or implicitely by introducing in the forecast covariance matrix a factor decaying gradually to zero as the distance between two grid points increases.

Weighted ensemble transform Kalman filter

We describe here how to perform the assimilation of observation y k at time k, i.e. how to update sequentially the filtering distribution from p(x k-1 |y

1:k-1 ) = N i=1 w (i) k-1 δ x (i) k-1 (x k-1 )
c 0000 Tellus, 000, 000-000

to p(x k |y 1:k ) = N i=1 w (i) k δ x (i) k (x k ).
The extension to the ensemble transform Kalman filter (ETKF) of the WEnKF described in Papadakis et al. ( 2010) is straightforward. It consists in sampling the ensemble members from the ETKF procedure and then to weight these members according to their likelihood (5). A final resampling is then performed, drawing a new set of N members from the ensemble with probabilities proportional to their weights.

Those steps are summarized below.

At time k, starting from a set of N particles {x

(i)
k-1 , i = 1, . . . , N }, one iteration of the WETKF algorithm splits into the following steps:

• Get x a,(i) k
from the forecast and analysis steps of ETKF:

-Forecast: For all i = 1, . . . , N , simulate x f,(i) k from model (1), with initial condition x a,(i)
k-1 and a given discrete scheme with time step δt; -Analysis: Compute x a,(i) k from Section 2.2.

• Compute particle filter weights w

(i) k ∝ p(y k |x a,(i) k
).

• Resample particles to obtain {x

(i) k , i = 1, . . . , N } and set w (i) k = 1 N for all i = 1, . . . , N .
Let us note that the final resampling step is implemented with a small Gaussian perturbation of the ensemble members. This random perturbation avoids a strict duplication of some members and implements an additive inflation procedure as it is usually done in ensemble Kalman filtering implementation. There is no particular constraint on the proposal ETKF step, it may thus includes any additional localization or inflation procedure. This extension could be for instance applied to localized ETKF version [START_REF] Hunt | Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter[END_REF]; [START_REF] Ott | A local ensemble Kalman filter for atmospheric data assimilation[END_REF].

Fluid motion estimation from images

In this section, we describe the application of the proposed WETKF method to the problem of fluid motion estimation from a sequence of images.

Dynamical model

In this work we will rely on the 2D vorticity-velocity formulation of the Navier-Stokes equation with a stochastic forcing to encode the dynamics uncertainty terms. Let ξ(x)

denote the scalar vorticity at point x = (x, y) T , associated to the 2D velocity w(x) = (w x (x), w y (x)) T through ξ(x) = ∂wy ∂x -∂wx ∂y . Let ξ ∈ R |Ω ξ | be the state vector describing the vorticity over a spatial domain Ω ξ of size |Ω ξ |, and w ∈ R 2|Ω ξ | the associated velocity field c 0000 Tellus, 000, 000-000 over the same domain. For incompressible flows, the stochastic dynamics is then written:

dξ t = -∇ξ t • wdt + ν∆ξ t dt + ηdB t ( 16
)
where ν is the viscosity coefficient and dB t is random forcing term. This a white in time noise model is an isotropic homogeneous random field corresponding to idealistic turbulence toy model [START_REF] Kraichnan | Small-scale structure of a randomly advected passive scalar[END_REF][START_REF] Majda | Simplified models for turbulent diffusion:theory, numerical modelling, and physical phenomena[END_REF]. It is characterized spatially by a covariance model Q exhibiting a power-law function spatial within an inertial range of scale.

These random fields are in practice simulated in the Fourier domain [START_REF] Elliot | A Fourier-wavelet Monte Carlo method for fractal random fields[END_REF][START_REF] Evensen | The ensemble Kalman filter, theoretical formulation and practical implementation[END_REF].

Note that the velocity field can be recovered from its vorticity, for incompressible fluids, using the Biot-Savart integral as:

w = ∇ ⊥ G ξ, where G(.) = ln(|.|) 2π
is the Green kernel associated to Laplacian operator and ∇ ⊥ = (-∂ y , ∂ x ) T represents the orthogonal gradiant.

This transformation can be efficiently implemented in Fourier domain, noting that ∆Ψ = ξ and w = ∇ ⊥ Ψ, where Ψ represents the stream function.

Nonlinear image observation model

There are two ways to assimilate the image data within the system. The first solution consists in computing pseudo-observations (motion estimates) from two consecutive images.

This leads to a linear observation model. Another approach consists in observing directly the image data, which makes the problem harder since the observed quantity is then a highly nonlinear function of the system state. As shown in the previous section, such a nonlinearity in the observation model can be handled within the WETKF framework.

In order to assimilate directly the image data, a nonlinear image observation model that relies on the brightness consistency assumption can be formulated as follows:

I(x, k) = I(x k+1 , k + 1) + γ k (x) ∀x ∈ Ω I , ( 17
)
where I is the luminance function, Ω I is the spatial image domain,

x k+1 = x + d(x) is the displaced point at time k + 1 and d(x) = k+1 k
w(x t )dt denotes the displacement between time k and time k + 1, with the initial point location at time k fixed to the grid point x (i.e.

x k = x).

Let I k be the vector gathering values I(x, k) for all x in the image domain Ω I . This vector constitutes the observation at time k. It is nonlinearly related to the state vector ξ k through the luminance function of the image data at time k + 1. The perturbation γ k c 0000 Tellus, 000, 000-000 is a Gaussian random field with covariance R. This measurement noise is assumed to be independent of the dynamics noise.

In the following, since this nonlinear image model will be used within the WETKF framework, the covariance R will be assumed to be diagonal in order to compute the analysis step more efficiently (see section 2.2). The uncertainty γ k (x) is then a Gaussian random variable of variance R(x, x) = σ 2 k (x). As the displaced image difference depends on the unknown displacement, this variance depends necessarily on the state variable (which is not contradictory with the independence property between the noise dynamics and the measurement noise). We propose to estimate empirically the variance from the ensemble members:

σ 2 k (x) = 1 N -1 N i=1 (I(x + d (i) (x), k + 1) -I d (x, k + 1)) 2 + , ( 18
)
where I d denotes the empirical mean of the displaced image. The parameter corresponds to a minimal variance value associated to homogeneous photometric regions. Let us note that in order to enforce some homogeneity between neighbor perturbations, a small Gaussian smoothing is applied to this variance.

The likelihood of observations (used to compute the particles weights) is then written:

p(I k |ξ (i) k ) ∝ exp   - |Ω| x=1 I(x, k) -I(x + d (i) (x), k + 1) 2 R(x, x)   . ( 19
)

Image assimilation scheme

We summarize in this section the WETKF image assimilation technique for fluid flow estimation.

At time k = 0, the vorticity ξ0 is initialized from a given motion estimator applied on the two first images I 0 and I 1 .

For k = 1, . . . , K, the following steps are performed to estimate ξk :

• Obtain {ξ (i)
k , i = 1, . . . , N } and associated weights w

(i)
k , i = 1, . . . , N from the WETKF algorithm described in section 2.3, with the following practical details:

-In order to compute the forecast step, estimated vorticity maps {ξ (i) k-1 , i = 1, . . . , N } at previous step have to be propagated from the dynamical model ( 16). This supplies the forecast ensemble members ξ f,(i) k

, at time k, as the Itô integral:

ξ f,(i) k = ξ (i) k-1 + k k-1 (ν∆ξ f,(i) t -∇ξ f,(i) t • w(ξ f,(i) t ))dt + η k k-1 dB t , (20)
c 0000 Tellus, 000, 000-000

where the velocity w(ξ f,(i) t

) is recovered from the Biot-Savart integral (see Section 3.1).

In practice, equation ( 20) is discretized with the Euler-Maruyama method [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] associated to a small time increment δt. A conservative numerical scheme [START_REF] Shu | Advanced Numerical Approximation of Nonlinear Hyperbolic Equations[END_REF] is used for the advective term ∇ξ • w, which respects the (flux) conservation law by integrating the flux values at cell boundaries. This total variation diminishing (TVD) scheme (monotonically preserving flux) prevent from an increase in oscillations over time and enable to transport shocks, as detailed in [START_REF] Shu | Advanced Numerical Approximation of Nonlinear Hyperbolic Equations[END_REF]. The Laplacian is implemented through a centered finite difference scheme.

-The analysis step of ETKF is performed with observation model ( 17), where the image data I k = {I(x, k) ∀x ∈ Ω I } constitutes the observation at time k, nonlinearly related to the vorticity state vector.

-Weights are computed from the likelihood ( 19).

-The variance of observation noise defined by ( 18) is not valid in missing data areas.

Such regions appear frequently in oceanic data due to the cloud cover (see the experiments on oceanic data in Section 4.3) or correspond to areas occluded by layers of higher altitudes than the layer of interest in meteorological images. The absence of measurement can in general be easily identified. For such regions we propose to impose a high uncertainty.

The idea here is to rely either on dynamics or on the neighborhood estimates rather than on false interpolated measurements. The strength of these uncertainties is set to a value proportional to the number of missing data over a centered local window. The magnitude of the proportionality factor is scaled to the highest uncertainty of the data (excluding the missing data). Moreover, for coastal regions, a mask corresponding to the land region is created with zero velocity values (Dirichlet boundary conditions). The dynamical update and measurement correction are frozen on such locations.

• Compute the vorticity estimate ξk as the weighted average ξk =

N i=1 w (i) k ξ (i)
k , and associated mean velocity fields estimates ŵk through Biot-Savart integration.

Experimental results

In this section, the proposed image assimilation method will be tested on synthetic and real sequences describing nonlinear and high-dimensional fluid phenomena. Results will be compared to another data assimilation method which can be used within this framework: the weighted ensemble Kalman filter (WEnKF, [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF]) which also relies on c 0000 Tellus, 000, 000-000 an ensemble Kalman filtering step, integrated within a particle filtering framework. In Papadakis et al. ( 2010), the WEnKF was used to assimilate linear observations. In case of image data that is nonlinearly related to the state of the system, these linear observations consist in pseudo-observations, i.e. velocity fields (and associated vorticity maps) computed from a given motion estimation technique. In the experiments below, these pseudo-observations are computed from a stochastic version of the well-known Lucas and Kanade motion estimator [START_REF] Corpetti | Stochastic uncertainty models for the luminance consistency assumption[END_REF], applied on each pair of the image sequence.

Let wk be the velocity field estimated by the stochastic Lucas-Kanade method, called SLK in the following, and ξk its corresponding vorticity. The observation model is written than as:

ξk (x) = ξ k (x) + γ k (x), ( 21
) with γ k (x) = σ k (x)η k (x)
where η k is a Gaussian random field with given covariance g λ (x-y), and σ k (x) is defined by the uncertainties of the motion measurements at each point x (obtained as an output of the SLK method). The additive uncertainty γ k is then a Gaussian random field with non-stationary variance R(x, x) = σ 2 k (x), and covariance R(x, y) = σ k (x)σ k (y)g λ (x -y).

The likelihood of observations (used to compute the particle weights) reads then:

p( ξk |ξ (i) k ) ∝ exp - 1 2 ( ξk -ξ (i) k ) T R -1 ( ξk -ξ (i) k ) . ( 22
)
As noted in [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF], a singular value decomposition of the |Ω I | × N observation perturbations matrix can be performed in order to compute efficiently this likelihood.

Results on a synthetic 2D image sequence

Our first set of experiments concerns a synthetic image sequence with images of size 256 × 256 showing the transport of a passive scalar by a forced 2D turbulent flow. This sequence is, to some extent, representative of typical satellite images depicting transport processes by oceanic streams, such as sea surface temperature or sea color images. This scalar exhibits a small diffusion, and does not respect strictly a luminance conservation assumption. Furthermore, due to low spatial gradients of this scalar in large areas (see experiments. For the nonlinear image measurement model, the observation noise covariance has been first fixed to a constant diagonal matrix with a unity standard deviation. Some trials with an image dependant observation covariance are shown at the end of this section.

Let us outline that the vorticity dynamics used for the filtering does not corresponds here exactly to a perfect dynamical model as the true deterministic forcing of the 2D turbulence is unknown and represented only by a zero mean random field.

The results are first compared in terms of global RMSE between the ground truth vorticity and the estimated vorticity at each time step, and similarly for the motion fields.

The WEnKF and WETKF results have been obtained with N = 700 particles. As may be observed on Figure 1, the improvement obtained with the WEnKF filtering technique over SLK measurements is minor, and even deteriorates the results in terms of velocity fields reconstruction. The proposed WETKF technique, directly based on image luminance observations, leads to much better results for the estimation of vorticity and velocity fields.

In addition to global error comparisons, the analysis of spectra allows to observe more precisely the accuracy of the different techniques at different scales. As a matter of fact the RMSE constitutes only a performance measure at large scales. It is indeed interesting to note that all methods give relatively coherent results only from the large scales up to the beginning of the inertial range. On the other hand, the result obtained with the WETKF assimilation scheme is closer to the ground truth over the whole scale range (from the largest scales up to the small dissipative scales). Again, this comes from the fact that this technique can take benefit of all available information in the image data, while pseudo-observations given by a local motion estimator will not be able to improve the estimation at all scales.

Filters based on pseudo-observations are unable to correct the loss of energy caused by the smoothing operator used in the external estimation procedure.

For a qualitative visual comparison, estimated vorticity maps are plotted on Figure 3 clearly lacks of energy compared to the true vorticity. The direct assimilation of image data through the WETKF scheme leads to a better estimation of vorticity structures, and in particular of small scales structures as discussed previously with the comparison of spectra.

To conclude, note that the improvement obtained with the WETKF technique comes not only from the direct introduction of image data within the assimilation system, but also from the chosen ETKF strategy for the ensemble Kalman step. We highlight this improvement on Figure 4, comparing the image assimilation technique with the ETKF scheme and EnKF scheme where the linear observation operator Hx k is replaced by H(x k ) as described in section 2.2 equ. ( 8) for the ETKF approach. Figure 4 shows that the best results are obtained with the ETKF scheme, for different sample sizes.

On this benchmark, we also compared in terms of vorticity the results obtained by WETKF and a similar ETKF version (non weighted implementation of the filter). The results corresponding to different amounts of ensemble members are gathered in tables 6 (a-b). These results correspond to the average mean square error obtained over the last 20 images of the sequence and 10 realizations. The first table shows the root mean square vorticity error associated to the empiral mean, whereas the second table gives the mean square dispersion of the ensemble around the ground truth (i.e. (N -1)

-1 i ξ a,(i) k -ξ ref 2 ).
It can be observed that the results for ETKF and WETKF remain of the same order. The first moment and the ensemble dispersion are sligtly lower for the weighted version. It is however important to outline that RMSE constitutes only a large scale criterion. To draw finer conclusions it is interesting to observe the error power spectrum (spectrum of the difference between the reference and the results averaged on the last 20 images and 10 realizations).

These curves are plotted in figure 5(a) in log-k×spectrum for different number of particles.

This graph depicts the amount of energy of the velocity error (w err = ww ref ) for each frequency band. For comparison purpose the ground truth spectrum and these curves are reproduced in the inside box graph. The error energy curves obtained by ETKF are plotted in dotted lines whereas those corresponding to WETKF are in plain line. It can be noticed that for both methods and for the different number of particles the error energy is mainly concentrated within a frequency band corresponding to 80 to 8 pixels and the error energy bump inside this range decreases significantly when the number of particles increases. The results obtained with 10 particles have much higher error energy than for the other amount of particles and are not acceptable. There is clearly a limit number of particles below which errors are too important. For 50 particles WETKF provides lower errors on almost all the c 0000 Tellus, 000, 000-000 whole frequency range that concentrates the most part of the energy. Another point of view can be obtained from error power spectrum in log-log scale normalized by the ground truth power spectrum. This is shown in figure 5(b). As can be observed for 10 particles there is no scale range for which the normalized error energy drops below 10%. For 50 particles this 10% error limit scale corresponds to 40 pixels. For 400 particles this scale reaches 25 pixels. Hence for both methods, low energy error levels are only obtained at coarse scales and accurate results at fine resolution requires to significantly augment the particles. For a weak number of particles WETKF seems to perform slighly better than ETKF whereas for a higher number of particles both filter provides very similar results.

The figure 8 depicts the spatial maps of the second order moment. For both methods the evolution of this moment is plotted. We can observe that WETKF provides in a faster way a smoother variance map and ETKF shows at some localized places a higher variance. The ensemble covariance values for the center line is also pictured (for 50 ensemble members)

for WETKF and ETKF. Both techniques provide comparable results in terms of covariance length scale (around 15 pixels).

As a last result on this sequence, we plot in figure 7 an example of the image-based observation variance (18). It is interesting to see that high variance value concentrates on thin structures corresponding to high photometric gradients. For this observation covariance and 400 particles we plot the evolution of the root mean square error. We notice this parameter free covariance observation model yields a better result than the previous constant variance.

The results described in this section comfort us on the use of the WETKF. The weight updating step and the resampling are not computationally expensive and strengthen the robustness of ensemble Kalman filter for relative low number of particles. In the following we will rely uniquely on such filters.

Results on real 2D turbulence

The next set of experiments corresponds to real experimental images of a passive scalar transported by a turbulent flow. This image sequence corresponds to the experiments described in [START_REF] Jullien | Experimental observation of batchelor dispersion of passive tracers[END_REF]. The size of images is 512 × 512, and WEnKF and WETKF results have been obtained with N = 300 particles. Vorticity results at a given time instant are plotted on Figure 9 for the SLK technique, and both WEnKF (assimilating SLK results) and WETKF (assimilating image data). We note that the SLK motion estimation c 0000 Tellus, 000, 000-000 the local motion estimation approach (SLK). We see that this initialization provides only a rough large scale motion field. This estimate is refined afterward by the filtering process.

We can note that the motion fields estimated along the sequence stick quite well to the big image structures observed on the SST images. The sequence of motion fields does not seem to be perturbed by the big missing data regions observed intermittently. Finally, we can note that the result obtained shows well an increase of the turbulent agitation with an intensification of the El-Niño phenomenon.

Let us note there seems to be a good adequation between the 2D vorticity dynamics used here and the data model based on (local in time) conservation of the observed quantity (the SST). This coupling enables extracting a 2D component explaining the image observation deformation and which is consistent with a 2D flow. The incorporation of more realistic 3D oceanic models would require more advanced data model allowing explaining luminance variations in a consistent way with respect to the dynamical models. In that prospect, the elaboration of data models relying on physical conservation laws might be of interest [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF].

Conclusion and discussion

In the context of image data assimilation, this paper extends to mean preserving square root filters [START_REF] Wang | Which is better, an ensemble of positive-negative pairs or a centered simplex ensemble?[END_REF] the weighted ensemble strategy proposed in [START_REF] Papadakis | Data assimilation with the weighted ensemble kalman filter[END_REF]. This technique is defined as a particle filter in which the state variable trajectories sampling mechanism (i.e. the proposal distribution) is defined from an ensemble Kalman filter. This type of proposal function revealed to be experimentally much more efficient than the traditional sequential importance sampling based on the dynamics alone. Furthermore compared to an importance sampling based on perturbed observations, a mean preserving square root filter sampling has shown to be more efficient for high-dimensional state space and a highly nonlinear observation equation. A significant improvement was brought by Mean square estimation errors with respect to the ground truth are presented on the first table, the ensemble dispersion around the ground truth is indicated on the second table. Note the ground truth is associated to a zero variance as it is a deterministic system. The mean has been computed through an empirical average over the 20 last image frames of 10 realizations c 0000 Tellus, 000, 000-000 9 Real 2D turbulent flow [START_REF] Jullien | Experimental observation of batchelor dispersion of passive tracers[END_REF]. c 0000 Tellus, 000, 000-000 c 0000 Tellus, 000, 000-000

  Those weights may be thus interpreted through (3) as the a posteriori probability associated to the sample trajectory x (i)

Fig. 3 -

 3 Fig. 3-(a)), most of the usual motion estimators perform quite poorly on this sequence.The parameters of the power law function involved in the model noise covariance are here fixed from the ground truth data. The noise samples have thus in this case the same energy spectrum as the reference. The model noise standard deviation has been set to 0.01 for all the c 0000 Tellus, 000, 000-000

  at a given time instant. The scalar image observation is first presented, together with the ground truth vorticity on Figure 3(b).Figure 3(c) shows the vorticity estimated with the local SLK technique. The result obtained with the WEnKF assimilation scheme based on these SLK observations is presented on Figure 3(d), while Figure 3(e) shows the result obtained with the proposed WETKF technique, assimilating image data directly. As can be seen on Figure 3(c), the vorticity estimated by the local motion estimation technique is far from the ground truth. As a consequence, the WEnKF assimilation based on these measurements only brings a small improvement since these observations do not carry enough information. The solution c 0000 Tellus, 000, 000-000

  technique captures only coarse-scale vortices. The scales are successfully refined with the WEnKF assimilation relying on the same measurement and with the WETKF assimilation incorporating directly the nonlinear image reconstruction error. The WETKF image assimilation reveals small motion scales exhibiting vortex filaments or stretching areas, which are smoothed out by SLK and WEnKF methods. This can be observed in particular inside the white delimited area of Figure9(d-e-f), zoomed in on Figure9(g-h-i). It must be pointed out that for the three techniques the estimated motion are quite close at large scale, so the results are consistent although based on different measurement operators. This might help validating these results, even if the ground truth is unknown here.4.3 Results on oceanic dataOur final set of results concerns the application of the proposed WETKF technique on satellite images of Sea Surface Temperature (SST), with large areas of missing data due to the cloud cover and presence of coastal regions.First of all, in order to highlight the capability of the method to cope with missing data, it is tested on a synthetic oceanic sequence. This sequence has been built from an initial real world 200×200 SST image transported by a synthetic flow, and respects the brightness consistency assumption given by (17). The synthetic flow is the same 2D turbulent flow as the one used in section 4.1. Along the whole sequence, four square holes with side 20-45 pixels have been added randomly on each image. Figure10shows an example at a given time k of the sequence. Figures10(a) and 10(b) show two consecutive images of the sequence, with the corresponding simulated ground truth velocity and vorticity on Figure 10(c). The result obtained with the WETKF assimilation scheme is presented on Figure 10(d) from the images without missing data, and on Figure 10(e) for the images with holes. It can be observed that the vorticity and velocity fields are estimated accurately, even in the presence of missing data. As a matter of fact, since for such regions the method relies on the dynamical model or on neighborhood estimates, the spatiotemporal consistency is maintained. The results can be compared on Figure 6 in terms of global RMSE between the mean estimates and the ground truth at each time step. The accuracy of the estimation deteriorates in presence of missing data, but the error level remains low and stable with time.The method is then tested on a real oceanic sequence of SST images. The sequence consists of 48 images of size 256×256. The images have a spatial resolution of 0.1 degree (10 c 0000 Tellus, 000, 000-000 km) and a temporal period of 24 hours. The sequence is centered on an area of the Pacific ocean off the Panama isthmus and shot during an El Niño episode. Representative results at different time instants of the sequence are shown in Figure12. The first column of Figure12shows the estimated velocity fields at times k = 1, 10, 24 and k = 39, superimposed on the corresponding SST images. The second column shows the estimated vorticity maps and velocity fields. The initialization for image 1 is based on the estimation provided by

Figure 1 .

 1 Figure 1. Synthetic 2D turbulent flow. (a) RMSE between mean estimate of vorticity and ground truth; (b) RMSE between mean estimate velocity and ground truth.

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Synthetic 2D turbulent flow. (a) Energy spectra; (b) Error energy spectra.

Figure 6 .Figure 8 .Figure 9 .Figure 10 .

 68910 Figure 6. Visualization(with 50 members) of the covariance matrix values corresponding to the central line at t = 99 left: ETKF, right: WETKF . The significative covariance length-scale is about 15 × ∆x large for both filters

Figure 11 .Figure 12 .

 1112 Figure 11. Synthetic oceanic sequence with missing data. (a) RMSE between mean estimate of vorticity and ground truth; (b) RMSE between mean estimate of velocity and ground truth.

Table 1 .

 1 Experiments on 2D turbulence. Results obtained from WEKF and ETKF are compared. For both filters the dynamics noise has a standard deviation of 0.01 and the parameters of the self-similar power law are fixed from the ground truth data.

	Results (mean estimation)		Results (ensemble dispersion)
	Number of members WETKF ETKF	Number of members WETKF ETKF
	10	0.145	0.146	10	0.146	0.147
	50	0.079	0.083	50	0.080	0.084
	150	0.051	0.052	150	0.053	0.053
	400	0.043	0.044	400	0.045	0.046

c 0000 Tellus, 000, 000-000

Acknowledgments

The authors thank the INRIA-Microsoft joint centre, and the ANR projects PRE-VASSEMBLE (ANR-08-COSI-012) and Geo-Fluids (ANR-09-SYSC-005) for their financial support. They acknowledge also the Ifremer CERSAT for providing them the oceanic SST image sequence.

this nonlinear observation over velocity and uncertainty measurements provided by a local motion estimation techniques [START_REF] Corpetti | Stochastic uncertainty models for the luminance consistency assumption[END_REF]. Pseudo-observation of flow velocities constructed from local motion estimator techniques such as the Lucas and Kanade estimator [START_REF] Lucas | An iterative image registration technique with an application to stereovision[END_REF] corresponds to measurements at too coarse scales.

These measurements should not be used in image based data assimilation. As those estimators have generally better performances than correlation techniques on scalar images (such as top of cloud pressure, sea surface temperature, sea color, etc.) the same conclusion applies also to maximum cross correlation or particle image velocimetry techniques. The nonlinear reconstruction error based on the displaced frame difference criterion has demonstrated to provide much better results. This criterion relies on a noisy transport assumption from frame to frame. This transport assumption is however only local since it applies only on the time interval in between two images. It is thus a weaker assumption than a strong transport assumption over the whole assimilation window used in [START_REF] Titaud | Assimilation of image sequences in numerical models[END_REF]. Besides it enables when coupled with ensemble techniques to access to empirical observation error covariance that can easily be adapted to missing information regions. This ability has allowed us to apply such an assimilation on a real noisy sequence of sea surface temperature. The results obtained for a simple 2D velocity-vorticity formulation of Navier-Stokes (with stochastic forcing) are very promising and demonstrate the possibility to estimate a coherent velocity fields sequence from difficult images exhibiting poor photometric contrast and large regions with missing data. In future works, we plan to extend such an image based assimilation strategy to more realistic geophysical models in oceanography or meteorology.

LIST OF TABLES

1

Experiments on 2D turbulence. Results obtained from WEKF and ETKF are compared. For both filters the dynamics noise has a standard deviation of 0.01 and the parameters of the self-similar power law are fixed from the ground truth data.

Mean square estimation errors with respect to the ground truth are presented on the first table, the ensemble dispersion around the ground truth is indicated on the second table. Note the ground truth is associated to a zero variance as it is a deterministic system. The mean has been computed through an empirical average over the 20 last image frames of 10 realizations c 0000 Tellus, 000, 000-000