Claire Thomas

Thomas Corpetti
email: tcorpetti@gmail.com

Etienne Mémin

Data Assimilation for Convective Cells Tracking on Meteorological Image Sequences

Keywords: Convective cells, data assimilation, levelsets, tracking

HAL is

Introduction

Convective cells are atmospheric events that are known to be associated with hazardous consequences, such as strong wind drafts, lightnings, heavy rainfalls, hails or even tornadoes. Over tropical areas such as central Africa, convective cells produce most of the rain during the monsoon period. They are also indirectly linked to droughts and floods, which might afflict this area. Their analysis and forecasting are thus of the utmost interest for meteorologists and forecasters.

The measurement of the convective cells parameters can be done either by the use of conventional probes or through satellite data. However, over remote areas (like the Amazonian forest for instance), conventional sensors such as radiosonde, rain gauges, radars or lightning detectors have a limited coverage. Satellite information constitutes therefore an appealing alternative for the study of convective atmospheric activity over these particular regions.

Except for large convective systems, the study of the convective cells development is particularly intricate since they are generally very sudden, short-lived and constitute highly deformable objects. As a consequence, the study of the convective activity requires high spatial and temporal resolution for precise diagnosis and forecasts. The new generation of geostationary space-borne sensors provides a wider variety of image channels with narrower bandwidths and an increased spatial and temporal resolution suitable for the analysis of synoptic to mesoscale phenomena. In this study we will rely on data of the Infrared (IR) and Water Vapor (WV) channels of the Meteosat Second Generation satellite (MSG-2). Compared to the first version, this sensor has been improved both on spatial resolution (1 pixel corresponding to 3km instead of 5km in the first generation) and temporal sampling rate (15min between two images instead of 30min in the previous system), in response to users requirements.

As introduced previously, this paper aims at studying methods for the reconstruction of a convective clouds system trajectory. Despite of the spatial and temporal resolution of satellite images, this task remains challenging, as images constitute only a 2D discrete sequence along time of a 3D phenomenon. Besides, as convective cloud systems exhibit a complex dynamics in which cells may split or merge, the 2D shape of such clouds system may change drastically between two successive frame instants. The reconstruction of a continuous trajectory to infer the whole history of a convective system in term of cells birth/death time, or of their activity phases constitutes hence an intricate issue that can hardly be handled by running "decoupled" segmentation processes on successive images.

In the past, several attempts have been proposed in the meteorological community for an automatic detection and tracking of these phenomena from satellite image data [1, 2, 17, 29-31, 33, 34, 39]. The first detection approaches were based on temperature threshold of IR images. These simple techniques permit to detect the cold core of convective cells. They, however, do not enable a reliable extraction of the cells contours, which are very sensitive to the threshold value. A prediction of the cells location roughly identified to the cold blob gravity centers has been then carried out through correlation techniques. Assuming an effective overlapping between the predicted cells and the detected ones at the currant instant, such techniques provide successive estimations of the cells core average location. Velocity and sometimes others parameters (such as the ellipticity factor or the distribution of the temperature gradient within the cells cores) are also provided.

In computer vision, a primal strategy to delineate the contours of complex objects relies on the implementation of partial differential equations encoding the evolution of parametric or non parametric curves toward the minimizers of an energy functional [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Chan | Active contours whithout edges[END_REF][START_REF] Kass | Snakes: Active contour models[END_REF]. The corresponding functionals include generally a data term representative of a photometric distribution characterizing the object of interest and some regularity constraints on the curves. These non linear deterministic minimization strategies are in practice very sensitive to the initial curve location and topology. Among the different extensions of those methods, the Eulerian "level-set" approach [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] has been specifically proposed to alleviate such shortcomings. In this framework, the contour's shape is represented as the zero level-set of an implicit higher-dimensional scalar function. The evolution of this implicit surface describes the contour's evolution and naturally enables the handling of topology changes. Such techniques have been already proposed in a meteorological context [START_REF] Cohen | Tracking meteorological structures through curve(s) matching using geodesic paths[END_REF][START_REF] Papin | Tracking and characterization of highly deformable cloud structures[END_REF][START_REF] Yahia | Segmentation of deformable templates with level sets characterized by particle systems[END_REF]. However they generally remain static detection/segmentation approaches applied (quasi-)independently on the images of the sequence with no temporal consistency guaranty on the recovered shape sequences. Time consistency is brought considering the estimation of the whole curves trajectory on the basis of a system of equations that couples an evolution model representative of the target dynamics and the complete image sequence. This constitutes a data assimilation problem that can be handled either relying on stochastic filtering or variational assimilation techniques. Stochastic filtering techniques implemented through particle filters have been proposed for the tracking of parametric curves [START_REF] Blake | Active contours[END_REF] or level set [START_REF] Rathi | Tracking deforming objects using particle filtering for geometric active contours[END_REF]. These techniques, which have the advantage of being recursive are nevertheless limited by difficult sampling issues in high dimensional state spaces. This limitation refered in the litterature as the curse of dimensionality, requires in practice the use of simple parametric descriptors for the curve (which prevents any topological changes) or the definition of simplified linear motion models. This latter simplification makes impossible the reconstitution of complex deformations and can hardly describe accurately the dynamics of convective cells. Compared to stochastic filtering techniques, the variational data assimilation framework [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF], provides a convenient procedure to handle high dimensional data assimilation problems. This technique belongs to the class of smoothing approaches, that operates on the whole sequence through forward and backward integrations of two coupled systems of pde's (the dynamical model and its adjoint). This process supplies a single trajectory, solution of the specified dynamics that fits at best a sequence of observations related to the state variable of interest.

In this work, we aim at specifying an image based variational data assimilation procedure [START_REF] Corpetti | Pressure image assimilation for atmospheric motion estimation[END_REF][START_REF] Papadakis | Dynamically consistent optical flow estimation[END_REF][START_REF] Papadakis | Variational optimal control technique for the tracking of deformable objects[END_REF][START_REF] Papadakis | A variational technique for time consistent tracking of curves and motion[END_REF] for the tracking of a convective clouds system. The definition of such a scheme requires to establish an observation operator linking the unknown curve with the available data and to define a dynamics that i) represents, as accurately as possible, the phenomenon of interest and ii) allows an efficient estimation of the unknown state variables. In our case this specification should describe the evolution along time of the 2D projection of the external frontiers of a convective cell and of its core. An accurate dynamics defined at the image grid scale that accounts for the different physical processes involved in those meteorological phenomenon is far to complex to be included in a practical image based assimilation process. A simplified "imperfect" model has then to be considered. Such assimilation scheme referred as weak constrained variational assimilation relies generally on a physical model simplification defined up to an uncertainty variable. In this work, we will comply with the same principle. However instead of sticking to a simplification of the cells dynamics, which would be in its own a difficult issue (as to our knowledge no satisfying evolution models are for the moment available), we will rely on a generic dynamics driven from image based motion measurements. This generic model will allow us to describe the evolution of two nested curves representing the cell cores and their external contours.

The simplicity of this model will be paid in return by the fact that as we do not have any dynamical models for the flow within or in vicinity of the cells, it will be impossible to set up a denoising data assimilation process of the velocity measurement as proposed in [START_REF] Corpetti | Pressure image assimilation for atmospheric motion estimation[END_REF][START_REF] Papadakis | Dynamically consistent optical flow estimation[END_REF]. Nevertheless, since the model considered incorporates an uncertainty related to the curve dynamics, we will show that the proposed model is robust to inaccurate motion measurements and does not require any specific assimilation of the flow velocity fields.

The structure of the paper is organized as follows: section 2 briefly describes the data used and justifies the choice of specific MSG channels; section 3 outlines the main principles of variational assimilation. In section 4, we describe the proposed data assimilation procedure for the tracking of two nested level sets delineating respectively the core and the boundary of the convective cell. The section 5 presents some experimental results and discusses the tracking robustness. Section 6 draws some conclusions on this work. Finally, the appendix opens new perspectives on the use of the estimated level set for the temporal analysis of convective cells.

Dataset

Channel for the convective cell detection

The first issue consists in describing properly the clouds designated as "convective cell" in order to define an adequate geometrical description together with photometric signatures that characterize them. In the literature the term "convective clouds" carries some ambiguities as it is employed to designate a multitude of very different meteorological phenomena such as: derecho, thunderstorm, multicell complex, supercell, squall line or convective cell. As explained by Morel and Senezi in [START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF], convective cells, or Mesoscale Convective Systems (MCS), are convective phenomena that take place at meso-α scales (250-2500 km) or upper meso-β scales (25-250 km). They have been studied for decades. In our study, we focus on the meso-β cells, and rely on the spatial description proposed for those cells by Houze in [START_REF] Houze | Mesoscale convective systems[END_REF]. A convective cell is a cloud producing low temperatures, and composed of two distinct regions: the actual convective part which consists of intense, coldest, vertically extended cores, and the stratiform region characterized by a more uniform texture and lighter precipitation. The stratiform area is partly produced by the dissipation of older convective cells and partly produced by broader sloping mesoscale layer ascent. The horizontal pattern formed by the convective core and the stratiform cloud part exhibit large variations along time and constitute both important indicators of the internal cloud dynamics. In this work, we will therefore describe a convective cell as two nested curves in order to distinguish the stratiform area from the actual convective part of the cloud.

The discrimination of convective cells is usually performed using a temperature threshold on IR imagery, since those events constitute the coldest events observed by satellite. In order to isolate these cloud systems, the image luminance needs to be transformed into temperatures. In general, physical temperatures of the sky can not be directly inferred from satellite imagery [START_REF] Szantai | Using multiple channels from msg to improve atmospheric motion wind selection and quality[END_REF]. However relying on a black body assumption an equivalent blackbody brightness temperature can be extracted for top of clouds observed through IR channels.

This assumption is valid only for thick clouds like the convective cells. However, other elements such as gases and aerosols located between the cloud of interest and the satellite sensor may absorb or diffuse a portion of the radiation emitted by the cloud. As a consequence, due to strong water vapor absorption in the WV 6.2 and WV 7.3 channels, these channels cannot be used for temperature determination, as opposed to IR bands, which are thus commonly adopted for the detection of convective system.

Relying on previous studies [START_REF] Arnaud | Automatic tracking and characterization of convective systems on Meteosat pictures[END_REF][START_REF] Houze | Mesoscale convective systems[END_REF][START_REF] Mathon | Life cycle of the Sahelian mesoscale convective cloud systems[END_REF][START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF][START_REF] Szantai | Using multiple channels from msg to improve atmospheric motion wind selection and quality[END_REF], we have selected the IR 10. Let us recall that, following the formation process image, the brightness temperatures are computed from the IR luminance values using:

T b = C 2 v c Aln C1υ 3 c R + 1 -B (1)
where C 1 and C 2 are constants respectively equal to 1.1910410 -5 mW m -2 sr -1 cm -4 and 1.43877Kcm.

Parameters A and B differ for each satellite/channel and their values are available on the EUMETSAT website http : //www.eumetsat.int for infrared and water vapor channels (IR and WV). The radiance R (in mW m -2 sr -1 cm) is received by the detector and is linearly related to the luminance value I by:

R = R 0 + αI, (2)
where α and R 0 are constants that might slightly change from day to day and whose values are provided by EUMETSAT. Let us note however that some very localized clouds associated to shorter lifetime (between 30 minutes and one hour) and a weaker convection activity exhibit brightness temperature below the threshold of -20 o C. Even if these clouds reach also the tropopause, they do not belong to the class of Mesoscale Convective Systems. To discard those isolated cloud storms, it is possible to discriminate them based on the lightnings occurences [START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF]. However, this requires additional data. In our study, this discremination will be performed naturally through the proposed assimilation process. Let us now present the data used for the motion estimation technique.

Channel for the optic flow estimation

As further presented in section 4, the assimilation process we built requires the estimation of velocity fields describing the top of the convective motion. On the basis of comparisons studies leaded on atmospherical wind field estimation from meteorological satellite images [START_REF] Endlich | Automatic cloud tracking applied to GOES and METEOSAT observations[END_REF][START_REF] Laurent | Measurement and validation of atmospheric motions detected on water vapor Meteosat imagery[END_REF][START_REF] Szantai | Using multiple channels from msg to improve atmospheric motion wind selection and quality[END_REF], we selected the WV channels to proceed the wind estimation. Let us point out that an important characteristic of pure water vapor structures is that they represent a mixture of radiation originating from different altitudes, located between pressure levels 700 and 200 hpa, (i.e. altitudes between 3 and 12 km). Consequently, WV channel provides on pure water vapor area an integration of structures at different heights. In the case of the convective cells, thunderstorms have a thickness of several kilometers and their wind motion information might strongly differ at the cloud's top and at its base. The use of WV channels for the determination of an average wind field related to the convective cells is hence pertinent. The wind estimation may however not be relevant in areas without clouds and with a poor water vapor content [START_REF] Laurent | Measurement and validation of atmospheric motions detected on water vapor Meteosat imagery[END_REF]. This has to be taken into account in the assimilation process.

Data Assimilation

In this section we will describe the general principles of the assimilation scheme we devised in this study.

Definition

Variational data assimilation is a technique derived from optimal control theory [START_REF] Lions | Optimal control of systems governed by PDEs[END_REF] to recover a state variable's trajectory from a sequence of measurements. Opposite to sequential Bayesian filters, which share the same aim, this framework allows to handle high dimensional systems and is thus intensively used in environmental sciences [START_REF] Bennett | Inverse Methods in Physical Oceanography[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF]. In this study, we rely on this framework to estimate the complete trajectories along an image sequence of closed curves representing convective clouds systems.

The data assimilation setup is basically composed of an ideal dynamics of the system variables, an initialization of the system's state and an observation equation that relates the system variables to some measurements:

               ∂ t X + M(X) = ν(x, t) X(x, t o ) = X o (x) + η(x) Y (x, t) = H(X(x, t)) + ǫ(x, t). (3)
The right hand side of the first equation describes, through a differential operator M, the evolution of the state function X(x, t) defined over the image plane Ω and the whole image sequence time range [t o ; t f ]. In our case, the components of this function will be composed of two implicit scalar surface functions representing two nested closed curves. This will be further detailed in the next section, but for sake of generality, let us first consider a general state function. The second equation sets up an initial condition for the state function through a given initialization X o (x). The last equation links an observation function Y (x, t) to the state function . This relation is formalized through a differential operator H. The observation function is usually composed of noisy and sparse measurements. In these three equations ν, η and ǫ are time varying Gaussian noise functions defined on the whole image plane. Since in our case the observation depends directly on the luminance function, I(x, t), whereas the dynamics depends on motion measurements estimated from the luminance spatial and temporal gradients, (∇I, ∂ t I), we assume that ν and ǫ are conditionally independent with respect to the image sequence. Besides, in order to simplify the estimation scheme and in order to fit to simplified variational assimilation schemes we will assume these noises are correlated in space but uncorrelated in time. This restriction could easily be dropped in the following but at the price of a much higher computational cost. Noises ν, η and ǫ are thus respectively associated to definite positive endomorphisms Q(x, x ′ )δ(t -t ′ ), B(x, x ′ ) and R(x, x ′ )δ(t -t ′ ) refered as covariance tensors. These tensors represent the errors involved in the different components of the system (i.e dynamical model errors, initialization errors and measurement errors).

Penalty function

The system of equations ( 3) could as well be specified through three conditional probability densities p(X(., t)|X(., t o ), I), p(X(., t o )|X o (.)) and p(Y (., t)|X(., t), I) related respectively to i) the distribution of the state trajectory X(., t) given the initial state value X(., t o ) and the whole image

sequence I = {I(., t), t ∈ [t o , t f ]};
ii) the distribution of initial state value given the initial condition X o (.) and the image sequence, and

iii) the distribution of Y (., t), the complete set of measurements, given the state X and the image sequence.

For the system we are considering the second distribution is Gaussian whereas the first and the third ones are Gaussians only if their associated operator (observation and dynamic) are linear. In practice the dynamical operator is highly non linear and due to for instance 2D-3D projections or indirect measurement of the quantity of interest, the associated observation operators we have to cope with in image are often non-linear.

Generally, it is considered that the less linear these operator, the more their distributions depart from Gaussian distributions.

As in any stochastic filtering problem, we aim here at estimating the conditional expectation of the state trajectory given the whole set of available measurements and the complete image sequence. A maximum a posteriori estimation consists in estimating the mode of the distribution p(X|Y, X o , I). Assuming Gaussian approximations of the pdf involved and uncorrelated errors in time, this leads to the minimization of a quadratic functional defined as:

J(X, ν, η) = 1 2 t f to Y -H(X) 2 R -1 dt + 1 2 X(t o ) -X o 2 B -1 + 1 2 t f to ∂X ∂t + M(X) 2 Q -1 dt, (4)
where we have introduced the norm X 2 A -1 associated to the scalar product in L 2 (Ω):

X, A -1 Y = Ω Ω X(x)A(x, x ′ )Y (x ′ )dxdx ′ .
When the assimilation relies on an exact version of the dynamics, the last term of this functional (the model part) disappears. In that case, the functional then depends only on the initial condition, and comes to an initial value control problem [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF]. This scheme is commonly referred in the literature as variational data assimilation and abbreviated as "4Dvar" whereas associated to dynamics with uncertainty variables, it is refered as "weak-constraint" variational data assimilation methodology [START_REF] Bennet | Generalized inversion of a global numerical weather prediction model[END_REF][START_REF] Bennet | Generalized inversion of a global numerical weather prediction model (ii). analysis and implementation[END_REF][START_REF] Egbert | Topex/Poseidon tides estimated using a global inverse model[END_REF][START_REF] Uboldi | Time space weak-constraint data assimilation for non linear models[END_REF][START_REF] Vidard | 4D-variational data analysis with imperfect model[END_REF]. Both schemes rely on the same adjoint optimization principles [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF] to minimize the functional. In the following we give a presentation of this principle.

A minimizer X of functional J is also a minimizer of a cost function J(X + βθ(x, t)), where θ(x, t))

belongs to a space of admissible functions and β is a positive parameter. In other words, X must cancel out the directional derivative:

∂ X J(θ) = lim β→0 dJ(X + βθ(x, t)) dβ = 0.
The computation of this gradient can be realized using an "adjoint formulation", as explained in the next section.

Functional gradient

Introducing an "adjoint variable" λ :

λ(x, t) = Ω Q -1 (∂ t X + M(X)) dx ′ , (5)
as well as the Gâteaux derivatives at point X of the operators M and H (called the linear tangent operators), denoted ∂ X M and ∂ X H respectively and defined by:

lim β→0 dM(X + βθ) dβ = ∂ X Mθ, (6)
the directional derivative of functional J(X, ν, η) reads:

lim β→0 dJ dβ = θ(t o ), B -1 (X(x ′ , t o ) -X o ) - t ∂ X Hθ, R -1 (Y -H(X)) dt + t ∂ t θ + ∂ X Mθ, λ dt (7)
Considering the three following integrations by parts, one can get rid of the partial derivatives of the admissible function θ in expression [START_REF] Blake | Active contours[END_REF], i.e.

t ∂ t θ, λ dt = θ(t f ), λ(t f ) -θ(t 0 ), λ(x, t 0 ) - t θ, ∂ t λ dt, (8) t
∂ X Mθ, λ dt = t θ, ∂ X M * λ dt, (9)
t ∂ X Hθ, R -1 (Y -H(X)) dt = t θ, ∂ X H * R -1 (Y -H(X)) dt. (10)
In the two last equations, we have introduced adjoint operators ∂ X M * and ∂ X H * as compact notations of the integration by parts of the associated linear tangent operators. Gathering all the elements we have so far, equation ( 7) can be rewritten as:

∂ X J(θ) = θ(t f ), λ(t f ) + θ(t o ), B -1 (X(t o ) -X o ) -λ(t 0 ) dt+ t θ, (-∂ t λ + ∂ X M * λ) -∂ X H * R -1 (Y -H(X)) dt. ( 11
)
This functional gradient must cancel out to satisfy to optimality conditions. This yields a "forward/backward" integration of the linear tangent dynamical system and its adjoint.

Forward/backward equations

Since the functional derivative must be null for any arbitrary independent admissible functions, the right members of the three scalar products of expression ( 11) must be identically null. It follows a coupled system of forward and backward PDE's with initial and final conditions:

λ(x, t f ) = 0, (12)
λ(x, t 0 ) = Ω B -1 (x, x ′ )(X(x ′ , t o ) -X o (x ′ ))dx ′ , (13)
∂ t X(x, t) + M(X(x, t)) = Ω Q(x, x ′ )λ(x ′ , t)dx ′ , (14)
∂ t λ -∂ X M * λ = - Ω ∂ X H * R -1 (x, x ′ )(Y -H(X))dx ′ . (15)
The forward equation ( 14) corresponds to the definition of the adjoint variable (5) and has been obtained introducing Q, the pseudo-inverse of Q -1 , defined as [START_REF] Bennett | Inverse Methods in Physical Oceanography[END_REF]:

Ω Q(x, x ′ )Q -1 (x ′ , x ′′ )dx ′ = δ(x -x ′′ )δ(t -t ′′ ).
Let us remark that for Q = 0 we retrieve the case of a perfect dynamics. Otherwise, equation ( 12) constitutes an explicit final condition for the adjoint evolution model equation [START_REF] Egbert | Topex/Poseidon tides estimated using a global inverse model[END_REF]. This adjoint evolution model can be integrated backward from this final condition assuming the knowledge of an initial guess for X to compute the discrepancy Y -H(X). To perform this integration, an expression of the adjoint evolution operator is required. Let us recall that this operator is defined from an integration by part of the linear tangent operator associated to the evolution law operator. The expression of such an operator can be obtained analytically or from the discrete scheme associated to the linear tangent operator. As a matter of fact, the adjoint of the linear tangent operator discretized as a matrix simply consists in the transpose of that matrix.

Knowing a first solution for the adjoint variable, an initial update condition for the state variable can be obtained from ( 13) and a pseudo inverse expression of the covariance tensor B. From this initial condition, equation ( 14) can be finally integrated forward.

It should be noted that equation ( 14) provides at convergence an expression of the model error:

ν(x, t) = Ω Q(x, x ′ )λ(x ′ , t)dx ′ . ( 16
)
The knowledge of the adjoint variable enables thus to estimate the error associated to the state variable evolution law. This may be sometimes useful to validate or invalidate a tracking result, or at least to infer some confidence measures on the result obtained.

This scheme somewhat differs from traditional weak constraint data assimilation. As a matter of fact, in geophysical applications, for computational reasons it is generally impossible to consider an uncertainty variable of the same dimension as the state variable. To control this variable a dimension reduction is generally applied or an explicit deterministic evolution model is imposed to the uncertainty model in order to come back to an initial value control problem [START_REF] Vidard | 4D-variational data analysis with imperfect model[END_REF]. In our case, the state space is much smaller than in those applications and a full uncertainty model can be applied.

Incremental state function

The previous system relies on a Gaussian assumption of p(X(., t)|X(., t o )) and p((Y., t)|X(., t)) (and implicitly thus on linear assumption of the model and observation operators). Besides, as due to the non linearity of the operators involved in the minimization described previously, this procedure constitutes a direct non-linear least squares estimation of the state variable. This process is likely to converge toward a bad local minimum if the initial guess is too far from the sought solution. The previous minimization can be slightly modified to follow a Gauss-Newton incremental methodology [START_REF] Courtier | Dual formulation of four-dimensional variational assimilation[END_REF][START_REF] Uboldi | Time space weak-constraint data assimilation for non linear models[END_REF], which relies on successive linearization of the operators (and therefore better fits to the Gaussianity assumption). This scheme is built introducing a function of state increments linking the state variable and a background solution: δX = X -X o . The linearization of the model operator M around the background X 0 :

M(X) = M(X o ) + ∂ Xo M(δX),
enables to split equation ( 14) into two pde's with an explicit initial condition performing a model integration from a first guess solution X o (x):

X(x, t o ) = X o (x), (17)
∂ t X o + M(X o ) = 0, (18)
∂ t δX + ∂ X0 MδX = Ω Q(x, x ′ )λ(x ′ , t)dx ′ . (19)
Combining equations [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF][START_REF] Corpetti | Dense estimation of fluid flows[END_REF] and [START_REF] Feidas | A software tool for monitoring the features of convective cloud systems with the use of Meteosat images[END_REF][START_REF] Feidas | Monitoring mesoscale convective cloud systems associated with heavy storms with the use of Meteosat imagery[END_REF][START_REF] Heas | 3D motion estimation of atmospheric layers from image sequences[END_REF] leads to the final assimilation algorithm. The method first consists of a forward integration of the background X o with a perfect system dynamics [START_REF] Feidas | Monitoring mesoscale convective cloud systems associated with heavy storms with the use of Meteosat imagery[END_REF]. The current solution is then corrected by performing a backward integration [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF][START_REF] Egbert | Topex/Poseidon tides estimated using a global inverse model[END_REF] of the adjoint variable. The evolution of λ is guided by a discrepancy measure between the observation and the estimate: Y -H(X). The initial condition is then updated through equation ( 13) and a forward integration of the increment δX is realized through the equation [START_REF] Heas | 3D motion estimation of atmospheric layers from image sequences[END_REF]. The overall process is iteratively repeated until convergence. A sketch of the whole process is summarized in figure [START_REF] Adler | A satellite infrared technique to estimate tropical convective and stratiform rainfall[END_REF]. As a Gauss-Newton optimization process this non-linear minimization may fail to converge if the initial condition of the incremental function is too far from the sought solution. In our case, we will show experimentally that even for an incorrect initial trajectory the system still show good convergence property. This nice behavior is due to the observations we considered, and that enables to correct efficiently the initial trajectory. These observations are spatially highly resolved and allow getting noisy but good representations of the clouds system at discrete instants. The assimilation process manage to removes the

Let X o (t 0 ) = X o . (i) From X o (t o ), compute X o (t), ∀t ∈ ]t o , t f [
with a forward integration of system [START_REF] Feidas | Monitoring mesoscale convective cloud systems associated with heavy storms with the use of Meteosat imagery[END_REF].

(ii) X o (t) being given, realize a backward integration of the adjoint variable with the system [START_REF] Egbert | Topex/Poseidon tides estimated using a global inverse model[END_REF].

(iii) Compute the initial value of the incremental function [START_REF] Corpetti | Dense estimation of fluid flows[END_REF].

(iv) From δX(t o ), compute δX(t), ∀t ∈ ]t o , t f [ with a forward integration of system [START_REF] Heas | 3D motion estimation of atmospheric layers from image sequences[END_REF].

(v) Update X = X + δX.

(vi) Return to (ii) and repeat until convergence.

Figure 1: Incremental variationnal data assimilation process noise or the false alarms associated to these observations and provides a continuous trajectory of the convective clouds system.

15

This section presents the principles of the variational data assimilation scheme we propose for the tracking of a convective cells system represented by a set of nested curves. In is organized as follows:

• the section 4.1 introduces the system state to assimilate;

• the section 4.2 is devoted to the dynamical model;

• the section 4.3 deals with the initial condition;

• the section 4.4 presents the observations and

• the section 4.5 is concerned with the definition of the error covariance matrices.

We now turn to the definition of the system state.

System state: convective cells representation

Among the different possibilities available to represent closed curves, the level sets formalism [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Theory, Algorithms, and Applications of Level Set Methods for Propagating Interfaces[END_REF] has proves its ability to manage in an elegant way the changes of topology (split, merge, ...) that are likely to appear. As the convective cells often undergo rapid changes in their topology and shape, such a formalism is well adapted to our problem. Within such framework, the curve at time t,

C t (p) = C(p, t) : [a, b] × R+ → R 2 ,
is implicitly described by the zero level set of a scalar function φ(x, t) : Ω × R+ → R:

C t (.) = C(., t) = {x ∈ Ω | φ(x, t) = 0},
where Ω stands for the image spatial domain. The surface is driven by the contours dynamics and is chosen so as to have for instance positive values inside the curve and negative values outside. A common choice for the implicit function is the signed distance function but any other surfaces whose level sets fits the curves of interest is possible. This representation has the great advantage to allow describing through a single implicit surface a set of non intersecting closed curves. Besides, geometric properties of the curve can be directly obtained from the implicit surface. The inward unit normal and the mean curvature are for instance respectively given by:

n = ∇φ |∇φ| and κ = div ∇φ |∇φ| . (20)
In order to encode the description of a convective cloud system encapsulating a set of convective cells, the system state X we used is then composed of two implicit surfaces:

• The first one, φ e , is devoted to the description of the external delimitations of the convective cells, C e (t).

These curves will be associated to a first level of brightness temperature threshold T e .

• The second one, φ c represents the cold core region of each cells, C c (t). It is related to a colder threshold value T c of the brightness temperature.

The first level set is defined on the whole image domain Ω:

φ e (x, t) : Ω × R+ → R,

Dynamical model

In this section, we first introduce (in 4.2.1) a common way to represent the evolution of implicit surfaces, where these latter are assumed to be driven by a smooth external motion field without any uncertainties. As will be shown, this assumption needs to be relaxed and we then propose a new and original dynamical model for the two nested curves where the motion uncertainties are explicitly taken into account (in 4.2.2). The resulting adjoint and linear tangent operators, some associated numerical aspects and the way to obtain the external motion fields are successively presented in sections 4.2.3, 4.2.4 and 4.2.5.

Implicit surface evolution without motion uncertainties

As commonly done when one deals with implicit representations, we assume that the curves of interest C i t , i = {e, c} are driven by a smooth external motion field w(x, t) along the normal n i and diffuse according to mean curvature flow [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Theory, Algorithms, and Applications of Level Set Methods for Propagating Interfaces[END_REF]:

∂C i t ∂t = (w i • n i + σκ i )n i , i = {e, c} (21)
where σ and κ respectively represent a diffusion coeficient and the curvature of the closed curve. As φ i (C i t (p), t) = 0, the implicit surface dynamics reads:

∂ t φ i + ∇φ i ∂C i t ∂t = 0. (22)
Replacing in this equation the expression of the curve velocity [START_REF] Houze | Mesoscale convective systems[END_REF] and of the normal [START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF], and extending this evolution law to the whole plane, we obtain finally the following transport equation:

∂ t φ i + ∇φ i • w i + σκ|∇φ i | = 0, (23)
where we have considered a velocity field w i (x, t) defined as an extension to the whole plane of the contour C i (., t) velocity. This representation provides an Eulerian representation of the evolving contours. As such, it allows to avoid the ad-hoc regridding processes of the different control points required in any explicit Lagrangian -spline based -curve description.

However, one can remark that this dynamics relies on (i) an ideal noise free external motion field and requires (ii) a smoothing parameter on the curves to guarantee a viscosity solution [START_REF] Sethian | Theory, Algorithms, and Applications of Level Set Methods for Propagating Interfaces[END_REF]. This parameter is not easy to fix, as it is not related to any physical quantity of the object of interest. We will see in the following section how this dynamics can be adapted to include some uncertainties on the motion field and how this provides us as a by product a more meaningful expression of the curve smoothing term.

Implicit surface evolution with motion uncertainties

In this section we introduce uncertainties on the motion field w of equation ( 21) that drives the two curves. To that end, we will rely on a stochastic formalism. For an introduction to stochastic calculus, we refer the reader to [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF].

We assume that the curves C i t are two stochastic processes driven by an Ito diffusion:

dC i t = (w i • n i )n i dt + σ n n i dB n t + σ τ n i⊥ dB τ t . (24)
Here n i dB n and n i⊥ dB τ are two independent Brownian motions directed along and perpendicularly to the curves normals and σ n , σ τ are the uncertainties respectively related to the normal and the tangent of the curve.

In this model, the curve is transported by a deterministic drift associated to the velocity field w (first term of relation [START_REF] Kass | Snakes: Active contour models[END_REF]) that is mitigated with isotropic Gaussian incertitudes -along the normal and tangent-whose covariances grow linearly in time (second and third terms of ( 24)). As previously, the curve location is defined from the evolution of an implicit surface φ i (x, t). As this surface depends on the stochastic process C t , its differential must be calculated using stochastic calculus differentiation rules, the so called Ito formula [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF].

We get the following stochastic partial differential equation (in the following, for sake of clarity, we will drop the curve's index i = {e, c} unless explicitly needed):

(∂ t φ + w • ∇φ + 1 2 {xi,xj }={x,y} ∂ 2 φ ∂x i ∂x j d C xi t , C xj t )dt + σ n |∇φ|d B n t = 0, (25)
where C x t and C y t correspond to the coordinates of the curve-points in the Euclidian plan at time t. As previously, we have introduced velocity drifts w(x) and Brownian motions Bn t (x) which are extensions on the whole plane of the curves drift and noise. For the drifts, we considered a unique smooth velocity field defined on the whole plane for both curves. This velocity field will be further detailed subsequently in section 4.2.5.

As for the noise, we considered a noise function defined from a set of independent Brownian noise variables localized on the points of a fixed discrete grid and a set of gaussian smoothing functions, ϕ ν , with variance ν and centered on the same points:

Bt (x) = i ϕ ν (x -x i )B n t (x i ). (26)
For a number of point tending to infinity, this function admits a limiting covariance tensor defined as:

C(x, y, t, s) = ϕ 2ν (x -y)δ(t -s). (27)
The notation

X t , X t = lim ∆ti ti≤t ((X ti+1 -X ti ) T (X ti+1 -X ti ))
involved in (eq. 25) denotes the quadratic variation of X t . Following stochastic calculus properties, this quantity is null for continuously differentiable functions and the following relation holds:

σ t dB t , σ t dB t = t 0 σ 2 t dt.
The term d X t , X t = dX T t dX t is computed according to the rule dtdB t = dtdt = 0 and dB t dB t = dt. In our case, we get:

d C x t , C x t = σ 2 n φ 2 x + σ 2 τ φ 2 y |∇φ| 2 dt, (28)
d C y t , C y t = σ 2 n φ 2 y + σ τ φ 2 x |∇φ| 2 dt, (29)
d C x t , C y t = (σ 2 n -σ 2 τ )φ x φ y |∇φ| 2 dt. (30)
Introducing the surface normal expression, the diffusion driving the implicit surface evolution reads finally:

(∂ t φ + w • ∇φ + 1 2|∇φ| 2 (φ xx (σ 2 1 φ 2 x + σ 2 2 φ 2 y ) + φ yy (σ 2 1 φ 2 y + σ 2 2 φ 2 x ) + 2(σ 2 n -σ 2 τ )φ x φ y φ xy ))dt + σ n |∇φ|d B = 0. ( 31
)
With the expression of the mean curvature:

κ = curv(φ) = 1 ∇φ (∆φ -∇φ T ∇ 2 φ ∇φ).
we get,

(∂ t φ + M(φ))dt + σ n |∇φ|d B = 0, (32)
with

M(φ) = w • ∇φ + σ 2 τ 2 κ|∇φ| + σ 2 n 2|∇φ| 2 ∇φ T ∇ 2 φ ∇φ, (33)
where ∇ 2 φ denotes the Hessian matrix. It can be observed from (33) that if both incertitudes have the same strength (i.e. σ n = σ τ ) this model takes a particular simple form, which reads:

(∂ t φ + w • ∇φ + σ 2 τ 2 κ∆φ)dt + σ n |∇φ|d B = 0, (34)
The dynamical model (32-33) constitutes a general stochastic process allowing to guide a curve through an implicit surface. Nevertheless, as in our work, we wish to stay within a deterministic framework, we will assume a Gaussian approximation of the state variable probability distribution. As the expectation of the dynamical model is null we get the following approximation for the model pdf:

p(φ t |φ 0 ) ∝ exp - t ∂ t φ + M(φ) 2 Q -1 (35)
where from [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF] the covariance tensor is

Q(x, y, t, s) = σ 2 n |∇φ(x, t)| × |∇φ(y, s)|C(x, y)δ(t -s). (36)
This prior defines our complete evolution model. Let us note that the inverse of this covariance tensor is not required in the assimilation process of relation [START_REF] Heas | 3D motion estimation of atmospheric layers from image sequences[END_REF] and does not have to be computed. Compared to the model in [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF] where no motion uncertainties were taken into account, the new formulation in [START_REF] Mcanelly | The precipitation life cycle of mesoscale convective complexes over the central United States[END_REF][START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF] is more generic since it includes both uncertainties along the normal and the tangent of the curve. It is interesting to observe that on the deterministic part of the model (first term of relation ( 32)), the motion incertitude along the curves tangent corresponds to a smoothing of the curve and leads to the mean curvature flow introduced in ( 23). The uncertainty along the curves normal generates also a smoothing on the basis of the Hessian ∇ 2 φ of the levelset.

In practice, the constant σ n and σ τ are fixed from an external observed displacements field between two consecutive images. More precisely, this field d(x, t) : R 2 × R → R 2 , given in practice by a procedure described in section 4.2.5 is assumed to be constituted by noisy measurements on the normal and the tangent of the photometric level-lines and of the form:

d(x, k) = w(x, k) + σ n ( B(x, k + 1) -B(x, k)) ∇I |∇I| + σ τ ( B(x, k + 1) -B(x, k)) ∇ ⊥ I |∇ ⊥ I| . ( 37
)
We are therefore making here the assumption that the observed velocity noises and the curve noises are collinear and have the same variances. Furthermore, we assume that the smooth drift velocity field, used to drive the curves, is such that w(x, t) = E(d(x, k))/dt and we approach this expectation through a local spatial mean over a small neighborhood V (x) of point x:

w(x, k) = 1 |V (x)| xi∈V (x) (d(x i , k)). ( 38
)
The variances of this observed displacement projected along the normalized curve normal and tangent provide an estimation of the noise variances σ 2 n and σ 2 τ . Let us now turn to the constitution of the tangent linear operator and of its adjoint necessary to the assimilation process application.

Adjoint and linear tangent operator of the evolution model

The tangent linear operator of the curve evolution model around a previous estimate φ required within the incremental Gauss-Newton assimilation process reads:

∂M ∂ φi (δφ i ) = ∇δφ T i w -σ t ∆δφ i + (σ n -σ t )    ∇ φT i ∇ 2 δφ i ∇ φi ∇ φi 2 -2 ∇ φT i ∇ 2 φi ∇ φi 2    ∇ φi ∇ φT i ∇ φi 2 -Id    ∇δφ i    . (39)
The corresponding adjoint can be computed from the adjoint of the discrete scheme used for the tangent linear operator. This strategy is straightforward and exact but nevertheless requires storing a lot of intermediate information at each grid point when advanced numerical schemes are used. A simpler scheme consists to operate a discretization of the adjoint of the tangent linear operator with the same discrete scheme. In our case this operator is defined as:

22 ! "! #!! #"! ! $! %! &! '! #!! !$! ! $! %! &! '!
∂M ⋆ ∂ φi (δφ i ) = -∇δφ T i w-σ t ∆δφ i +(σ n -σ t )    ∇ φT i ∇ 2 δφ i ∇ φi ∇ φi 2 + 2 ∇ φT i ∇ 2 φi ∇ φi 2    ∇ φi ∇ φT i ∇ φi 2 -Id    ∇δφ i    . ( 40
)
Both strategies are equivalent when the discrete schemes used to implement the adjoint operator constitutes the adjoints of the discrete schemes used for the direct model. This is the strategy we adopted in our implementation. Exact discrete adjoints operators of the direct model have been used for the adjoint evolution law. Some hints on the numerical schemes on which we relied are given in the following.

Numerical schemes

In practice, the two level sets of the system state and the velocity field w are defined on the image grid. The direct dynamical system [START_REF] Mcanelly | The precipitation life cycle of mesoscale convective complexes over the central United States[END_REF][START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF], the backward adjoint equation and the forward incremental model associated respectively to operators [START_REF] Papin | Tracking and characterization of highly deformable cloud structures[END_REF] and ( 40) are integrated using an explicit discretization scheme with a third order Runge-Kutta scheme that respects a Total Variation Diminishing property [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. The advective term of the level set dynamics has been also implemented through a TVD numerical scheme in order to achieve an accurate and stable discretization. We employed a semidiscrete central scheme with a numerical convection flux derived from the monotone Lax-Friedrichs flux [START_REF] Kurganov | A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations[END_REF][START_REF] Xu | Anti-diffusive finite difference weno methods for shallow water with transport of pollutant[END_REF]. The second order operators have been implemented through traditional centered discrete schemes.

Observed displacements field

The dynamical system [START_REF] Papin | Tracking and characterization of highly deformable cloud structures[END_REF] we consider depends on the local mean of an external displacements field (cf relation [START_REF] Papadakis | A variational technique for time consistent tracking of curves and motion[END_REF]). This field can be estimated in practice from the image sequence using a fluid motion estimator.

The estimator on which we relied in this study is formalized as the following minimization problem [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF][START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]:

arg min d f ( ∂I ∂t + ∇I • d + Iζ) + α(f (|∇ζ|) + f (|∇ξ|)) dx, (41)
where I = I(x, t) stands for the image luminance function, d = d(x, t) = (u(x, t), v(x, t)) T is the velocity field to estimate, α a smoothing parameter and ζ = ∂u/∂x + ∂v/∂y (resp. ξ = ∂v/∂x -∂u/∂y) is the divergence (resp. the vorticity) of the flow. The function f (•) is often the quadratic norm but it can either be a "robust" or "softer" penalty function that enables to deal with corrupted data (see [START_REF] Huber | Robust Statistics[END_REF] for some introduction to robust statistics). Roughly, this states that the velocity field to estimate should satisfies 1. the mass conservation law applied to the image luminance: ∂I/∂t + ∇I • d + Iζ ≈ 0 that assumes that the variation on the density is related to the diverging motions. It has been proved in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF] that this assumption holds in IR imagery and ;

2. a spatial consistency that allows the estimation of highly diverging and rotating motion fields, unlike usual optical flow techniques based on a first order smoothing of the components u and v.

Variations of such schemes has been successfully used for the estimation of atmospheric 2D or 3D wind fields [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Heas | 3D motion estimation of atmospheric layers from image sequences[END_REF][START_REF] Heas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF]. An example of a motion fields estimated through such an estimator on IR images is shown on figure 3.

Initial surface

The initial surface is in practice often chosen as a signed distance function to the contour. Since a distance function is characterized by a unity gradient modulus ( ∇φ = 1), this strategy leads, following the same strategy than [START_REF] Papadakis | Dynamically consistent optical flow estimation[END_REF], to a constant covariance at the initial time and hence does not appear to be an efficient choice in our case.

As in our application a simple thresholding technique provides the initial contours, it appears more natural to infer the error covariance between the initial state and the initial contour (the background error) from the data, and besides to guarantee some consistency with the model error covariance at the initial time. This can be simply done fixing the initial surface as a smoothed temperature surface observed at the same moment.

This smoothing is operated in practice through a Gaussian convolution with a small standard deviation and guaranties differentiability properties of this surface. Furthermore, this initialization has the great advantage to include in the implicit surface structural informations of the temperature map before an eventual cell formation.

This avoids in empty cells areas imposing level sets values that depend only on the geometry of the other detected active cells -or worst on an artificial boundary when no cell has been detected. Figure 2 shows an example of the benefit provided by this initialization.

Observations

A common way to track objects in images consists in relying on a reference histogram. However in our application, such a generic description of the regions of interest is difficult to specify as the spatial organization of the temperature inside a convective cell strongly evolves along time. This evolution is depicted in figure 4.

This figure illustrates a typical evolution of the temperature histograms observed along time for a given of a convective cell. From this illustration, it is obvious that a fixed reference object histogram would be strongly sensitive to the instant of acquisition of the first image of the sequence. We thus propose an observation operator based on the temperature threshold related respectively to the convective cell contour and the cell's core. Assuming, as it is commonly done in convective cells detection issues, that the core's activity (resp. its external area) lies in temperatures under a given threshold T c (resp. above T e ), the function φ, whose zero-level corresponds to the border to extract, should be such as:

H(I(x) -T • )(1 -H(φ(x))) + H(T • -I(x))H(φ(x)) = ǫ, (42)
where • = {e, c} and H is the Heaviside function (H(x) = 0 if x < 0 and H(x) = 1 elsewhere). This relation is zero if I(x) > T • and φ(x) > 0 or if I(x) < T • and φ(x) < 0. This measurement model favors thus the curve interior to be above of a given threshold and the exterior to be lower of the same threshold. Otherwise this function is always positive. Hence, this relation enables to define the observation system Y = H(φ) as:

       Y = 0 H(φ) = H(I(x) -T • )(1 -H(φ(x))) + H(T • -I(x))H(φ(x)) (43)
The operator H is non-linear with respect to φ. Its associated tangent linear operator yields:

∂H ∂φ (δφ) = [-H(I(x) -T • ) + H(T • -I(x))] δ(φ), (44)
where δ(φ) is the dirac function applied to φ.

Error covariance tensors

The dynamical model error covariance matrices Q and the background error covariance B have been defined in section 4.2 and 4.3 respectively. They are account for the quality of both the measured velocities from the image sequences and the input luminance data. As for the observation error, a similar procedure has been applied: since the observations are extracted from an image threshold independently at each point of the image grid, the error observation covariance tensor has been defined as a diagonal matrix whose values are directly linked to the surface gradient magnitude (i.e a thresholding measure is associated to high covariances value and hence low uncertainties -resp. a low value and high uncertainty-for high magnitude of the temperature gradients -resp. low temperature variations-). The proportionality coefficient used in this covariance matrix constitutes the single parameter of our assimilation process. In practice this coefficient has been set to 1 in all the experiments on which the technique has been run.

Let us now present some experimental results.

Experimental results

This section describes some experiments that have been carried out on real data of the Meteosat Second Generation satellite. Before showing the results, let us first describe the datasets that have been used as

Data and geographical area

Our different image sequences were extracted from August 8th 2005 (time 00h00) to September 1st (time 23h59). The upper left hand corner coordinates of the geographical area are (0,-20) degrees in latitude and longitude, and [START_REF] Kurganov | A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations[END_REF][START_REF] Rathi | Tracking deforming objects using particle filtering for geometric active contours[END_REF] at the opposite point. The region comprizes a great part of central Africa, which extends from the Sahel in the northern part of the image to the equatorial region in the southern part. These data occurred during the Sahel rainy season, spreading from June to September [START_REF] Arnaud | Automatic tracking and characterization of convective systems on Meteosat pictures[END_REF], and during the monsoon rains, which takes place in the southern part. The following bands were used: WV 6.2, IR 9.7, IR 10.8, and IR 13.4.

IR images were transformed into brightness temperatures using equation (1).

We have selected two sequences, respectively named "single cell" and "multiple cells". The figure 5 after the assimilation process. The contour of the cell is represented in black while the core is in white.

Experiments on the single cell sequence

As explained previously, in order to ensure a non overlapping between the two curves represented by φ e and φ c (respectively related to the external contour and the core of the cell), the surface φ c is defined over a support corresponding to the interior of the domain delineated by zero level set of φ e . Its assimilation depends therefore on the estimation of surface φ e . The estimation of the cell boundary surface φ e does not depend however on the core surface and can be performed first. The assimilation of the core surface is then handled afterward.

The assimilation results are shown in Figure 6. The first line represents five images of the sequence where the propagation of the initial contours through the dynamical model are superimposed (the contour in black represents the external cell and the one in white is related to the core). The second line depicts the results after the assimilation process. On the first row of this figure, one can observe that the dynamical model provides already a coarse approximation of the final contour location for the two initial cells (that are merging later in the sequence). During the initial integration, we see that the both cells evolve independantly from each other along the sequence (first line of Figure 6). The benefit of the assimilation process then clearly appears on the second line of the figure 6 since the two initial cells merge together in a unique cell. It can be observed that the estimated curves related to the contour and the core of the cell are well in line with the visualized sequence. They also are in accordance with the usual behavior of these convective systems. Let us now study the relevance of out approach when several cells split and merge.

Experiments on the multiple cell sequence

This second experiment represents a more complex scene since it contains several cells at different degrees of evolution. In addition, as few convective phenomena are initiated in the first image (i.e. only the biggest cells are detected through the threshold technique on the first image), the first integration of the dynamical system has provided only a rough initial system state for the most proeminent cells: the very first forecast of the initial flat temperature surface (first line of Figure 7) corresponding to the initial condition indeed obviously failed in predicting a complete series of contours for the different cells observed. The results after the assimilation process are shown in Figure 7. Despite an incomplete and unsatisfactory initial solution, we can observe that the process enables to estimate accurately the contours and the cores of the different convective cells. The different split and merge were correctly managed (the complete sequence can be seen at http://www.sites.univ-rennes2.fr/costel/liama/corpetti/Thomas Corpetti eng/ANR MSG.html).

These two experiments on real meteorological data demonstrate the ability of the proposed assimilation process to recover reliable curves describing the contour and the core of the observed convective cell.

To highlight the robustness of the assimilation procedure, let us now study the relevance of the proposed technique in presence of highly corrupted data.

Robustness in presence of noisy and corrupted data

The Meteosat Second Generation data are received, analyzed and stored by the EUMETSAT consortium. For some specific applications, it is also possible to install a simple reception station using a conventional TV antenna. The MSG signal is captured and transform into a luminance brightness. In such simple reception stations, it is frequent that the signal is corrupted due for instance to interference with aerosols, clouds, ... Therefore, many images are likely to be unusable. To analyze the stability of our approach in presence of outliers, we have performed the assimilation on a corrupeted version of the single cell sequence where the following transformations have been applied:

• the fourth image has been half bottom cut;

• a Gaussian-noise with a known standard deviation has been introduced in the seventh;

• the same noise has been also introduced in the tenth, and the top left-hand quarter has been discarded;

• the last image has been suppressed.

On this sequence, the motion field have been estimated using the cost function of relation [START_REF] Scofield | The NESDIS operational convective precipitation-estimation technique[END_REF]. Thanks to the robust norm f (.) used, the outliers areas are "rejected" of the minimization procedure and the solution consists in a regularization of the motion field by propagating its divergence and vorticity on the "noisy" data. The Figure 8 presents the results after the assimilation superimposed to the original images. The method turns out to be very efficient since the estimated curves are still very closed from the ones obtained from the uncorrupted data. On all the noisy areas, the associated covariance matrices are very high since the observation provides no relevant information. As mentioned in previous sections, these latter are indeed related to the magnitude of the gradient of the image luminance. Therefore, on such outliers regions, the observation adds no information and the system state results in a propagation of the dynamical system. This experience emphasizes well the robustness of the presented approach with respect to different potential sources of degradations of MSG data.

Conclusion

In this paper we have proposed an original and efficient variational assimilation framework for the tracking of convective cells from both IR and WV MSG imagery. The technique relies on an optimal control formulation which consists to monitor the trajectory of the curves descriptors from the complete sequence of data. Each step of the process is composed of two iterative phases: the forward linear development of the state variable from a dynamical model, and the backward integration of an adjoint dynamics incorporating a data measurement model. A time consistent family of state variables is by this way generated as a satisfying trade-off with respect to uncertainties of the system dynamics and of its observation. Excepted the proportionality coefficient used in the background error covariance, all the other parameters involved in the proposed data assimilation process are estimated from the image data and do not require any ad-hoc tuning.

The cells are characterized using a state variable composed of two levelsets, related respectively to their boundary and their core. The level set formulation enables a natural topology adaptability that fits well to the possible splitting and merging phases underwent by a convective cell during its life cycle. This approach is also relevant to properly process such highly deformable cloud structures. The use of two level sets allows identifying the cold clouds from the other surrounding and warmer ones, but also enables to focus on the actual convective activity part of the convective cell and to distinguish it from the whole cloud.

The dynamical model associated to the level sets assumes that the curves are propagated by an external velocity field (in practice estimated on the image but that can be provided by any measurement techniques)

where the uncertainties along the normal and the tangent of the curve are taken into account. This provides a generic stochastic process allowing to guide the curves through an implicit surface. Of course, the quality of this dynamical model is related to a global accuracy of the velocity fields but local velocity errors are handled in our scheme.

We have outlined the strength of the variational assimilation formulation in its ability to deal with missing or altered data. Furthermore, it provides a continuous framework for the state variable tracking, and provides therefore a continuous trajectory of each cell. This enables us to provide a cell localization even in case of missing or altered images and to oversample the final trajectory (with regard to the image temporal lapse rate).

In the intricate situation where a cell undergoes multiple merges and splits simultaneously, this decomposition of the trajectory into shorter intervals of time provides a clearer temporal interpretation of the diverse successive topological changes.

We believe that such a tool is very relevant for the climatologist and meteorology community since the analysis of convective cells is still a very intricate problem.

In order to explore further the potentialities of the level sets assimilation, we present in the appendix a tool that labels and estimates several of the geometric and radiative features of the cells that characterizes their development's phase. This is based on several automatic and semi-automatic tools encountered in the literature. The association of the assimilation results to the convective cell activity characterization constitutes • the averaged wind speed vector of the center of gravity, weighted by the temperatures of the cell.

In addition to the labelisation and characterization of each cell, a "contextual" parameter has been added

to precise the origin the cell:

• Normal birth: the cell has just been created and has no antecedent;

• Normal evolution: the cell has a unique antecedent;

• Merge: the cell has several antecedents;

• Split: several cells have the same antecedent;

These labels allow us memorizing the cell history. For each cell, we are able to identify its date of birth and death or its type of evolution (merge, split, ancestor, filiation, ...). To illustrate such labelisation, the figure 9 represents for the second sequence a picture of the labeled cells at different stages of development.

The following structure is an example of a cell identity card extracted on this sequence:

• Label number: 2

• Born at image #: 4

• Dead at image #: 7

• Type of born: NATURAL

• Futur merge: YES

• Futur split: NO

• Merge with cell(s): 13

• Overlapping (in pix): 22

• Average temperature along life: [347.5, 343.5, 340.7, 339.0]

• Surface along life: [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Lions | Optimal control of systems governed by PDEs[END_REF]53,76] • Perimeter along life: [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Papin | Tracking and characterization of highly deformable cloud structures[END_REF] • Average divergence along life: [0.1, 0. In addition to such an "identity card" for each convective cell, it is also possible to characterize the activity of a convective system in the sens of the definition given in [START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF] that links to the same convective cells future merge cells. As already mentioned, we have assumed that a convective cell can be characterized into three phases : growing, maturity and decay. Following [START_REF] Papin | Tracking and characterization of highly deformable cloud structures[END_REF], this convective activity can be characterized considering the joint evaluation of the local divergence contained in the 2D motion fields and the temporal variations of the brightness temperature. For instance, the growing activity corresponds to the time interval in which the averaged temperature of the cell decreases while the averaged divergence is positive; the maturity is associated to a calming down activity where the top cloud stops cooling down and its averaged temperature can even start rising up slightly whereas the decay is characterized by a period where the divergence is negative in conjunction with an increasing temperature. Figure 10 illustrates the phase characterization for the cell labeled #2. We have also added both indicators (evolution of the mean divergence and temperature inside the cell) that enables to characterize its activity.

With such available data, it is possible to print out other useful characterization indicators. For instance, the figure 11 exhibits some of common descriptors used in different previous studies for the characterization of convective cells [START_REF] Arnaud | Automatic tracking and characterization of convective systems on Meteosat pictures[END_REF][START_REF] Feidas | A software tool for monitoring the features of convective cloud systems with the use of Meteosat images[END_REF][START_REF] Feidas | Monitoring mesoscale convective cloud systems associated with heavy storms with the use of Meteosat imagery[END_REF][START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF]. From these plots, one can see that even if some indicators are correlated in some phases (for instance, a positive divergence is related to a positive gradient surface evolution), they contain all a specific information. The figure 11 (c) exhibits the quantity 1 A ∂A ∂tdiv(v i ) where A stands for the area of the cell. This index, which links the area evolution to the average divergence has been introduced in [START_REF] Machado | Seasonal and diurnal variability of convection over the Amazonia: a comparison of different vegetation types and large scale forcing[END_REF], is null during the growing phase. The others indexes, related to the surface, perimeter, volume and temperature, are estimated from the estimated curves and the associated motion fields. They may provide precious information that can be used for specialists to better analyze and understand the convective cell's behavior. Such summary cannot be robustly extracted without the tracking and denoising process performed through the curve assimilation.

  8 channel for brightness temper-ature estimation. Two thresholds (T 1 = -20 o C = 253K and T 2 = -60 o C = 213K) are used to separate respectively the convective cells and the core from the other clouds.

  whereas the second one is defined inside the domain encapsulted by the convective cells contours, I = {x | φ e (x, t) > 0} (i.e. the domain delineated by the zero level set of φ e ):φ c (x, t) : I × R+ → R.The contours of the convective cores, represented by the zero level set of the second implicit surface, are thus constrained to lie inside the domain delineated by the external frontiers of the convective cells set. Both implicit surfaces constitute the components of the state variable of our data assimilation problem. Let us now define the dynamical model.

Figure 2 :

 2 Figure 2: Differences between two strategies of initialization of the level-sets. (a): a meteosat image including a convective cell on the right part; (b) the corresponding initial level-set obtained with a signed distance and (c) the corresponding initial level-set obtained based on the temperature. One can see that the apparition of a new cell in the left part of the image will be less favored with the signed distance function (b).

Figure 3 :

 3 Figure 3: Illustration of the estimation of the velocity fields. (a) a convective cell image superimposed with the estimated velocity and (b) the velocity field.

Figure 4 :

 4 Figure 4: Evolution of the non-normalized temperature histogram of a convective cell sequence. Each line corresponds to a representation of the temperature distribution at a given instant of the sequence. The

Figure 5 :

 5 Figure 5: Examples of IR 10.8 images on the two different sequences of the benchmark.

  depicts one image of two different selected sequences. In the single cell sequence (figure5(a)), a single large convective cell is visible. With a diameter of approximately 240 kilometers at its mature stage and a temperature below -60C, this convective cell belongs to the category of mesoscale events also called super cells. The multiple cells sequence of figure5 (b) is more complex. It exhibits numerous diagonally organized smaller cells, and might be more considered as a multicellular storm than as a true convective cell. Indeed, the biggest cell's diameter does not spread more than one hundred kilometers, and the temperature does not dive much below the core temperature threhold.

Figure 6 :

 6 Figure 6: Results on the single cell sequence. The first line corresponds to the contours obtained after the propagation of the initial contour with the proposed level set's dynamics; the second line shows the final results

Figure 7 :

 7 Figure 7: Results on the multiple cell sequence. The first line corresponds to the contours obtained after the propagation of the initial contour with the proposed level set's dynamics; the second line shows the final results after the assimilation process. The contour of the cell is represented in black while the core is in white.

Figure 8 :

 8 Figure 8: Illustration of the assimilation results on corrupted data; the half bottom part of the fourth image has been cut out, the seventh has been blurred by a known Gaussian noise, the tenth has a quarter missing and is

Figure 9 :

 9 Figure 9: Example of labeled convective cells computed from the assimilation results along the sequence.

•

  2, 0.1, 0.3] • Volume index along life: [100, 463, 1021, 1593] • (Volume index)/(Surface) along life: [12.5, 16.5, 19.3, 21.0] Gravity center along life: [(140,46), (140,45), (140,44), (140,42)]This expresses that the cell number #2 triggers on image # 4 with a normal birth (it does not result from any merge or split of previous cells) and died on image # 7 of the sequence by merging with another cell to create the cell # 13. The overlapping between cell # 2 and cell # 13 is 22 pixels. As the cell # 2 is visible on four images, four values for each parameter field have been extracted.

Figure 10 :

 10 Figure 10: Characterization a convective cell's activity. For a given cell (cell labeled #2 in this figure), an image corresponding to each phase and key events is represented. The mean displacement of the cell is also

Figure 11 :

 11 Figure 11: Area and gradient of the area of the cell 2 of the second sequence along the 23 images of its

Acknowledgment

This work was supported by the French National Agency of Research through the ANR-06-JCJC-0083 STANDS-MSG project.

APPENDIX APPLICATION TO THE ESTIMATION OF CHARACTERISTIC PARAMETERS OF THE CONVECTIVE CELLS

Many studies have demonstrated the interest of characterizing each convective cell individually, as well as its phase of development. It is common to characterize the evolution of a convective cell by three states: growing, maturation and decay. According to Barrett and Martin [START_REF] Barrett | The Use of Satellite Data in Rainfall Monitoring[END_REF], each step of the evolution has its own particularity. For instance, a cloud produces less rain in its decreasing stage than during its growth. Thus, the characterization of the phase should improve the rainfall estimation [START_REF] Mcanelly | The precipitation life cycle of mesoscale convective complexes over the central United States[END_REF]. In addition, several indicators related to the stage of development are very useful for specialists. As an example, in [START_REF] Scofield | The NESDIS operational convective precipitation-estimation technique[END_REF], the authors have presented a technique to estimate precipitations using satellite images based on the "convective cloud area evolution" index. More generally, several attempts have been carried out in order to provide to the community automatic or partially-automatic detection and tracking tools [START_REF] Arnaud | Automatic tracking and characterization of convective systems on Meteosat pictures[END_REF][START_REF] Feidas | A software tool for monitoring the features of convective cloud systems with the use of Meteosat images[END_REF][START_REF] Feidas | Monitoring mesoscale convective cloud systems associated with heavy storms with the use of Meteosat imagery[END_REF][START_REF] Morel | A climatology of mesoscale convective systems over europe using satellite infrared imagery. i: Methodology[END_REF]. In this section, we aim at showing the interest of the assimilation for the characterization of the life cycle of a given convective cell. From the different curves extracted with the proposed assimilation process, we present a tool that allows us to label each convective cell, to memorize its history (origin and new cells issued from the current), and to compute several significant parameters for the characterization of the convective phenomenon. These parameters are:

• the averaged brightness temperature;

• the surface;

• the perimeter;

• the averaged divergence;

• the volume index [START_REF] Szejwach | Dynamic classification of mesoscale cloud patterns[END_REF]: n i (I(x) -T ) where n i is the number of pixels with the same brightness temperature;

• the gravity center coordinates;