O Titaud

A Vidard

I Souopgui

F.-X Le Dimet

Assimilation of Image Sequences in Numerical Models

Keywords: Data assimilation, Inverse problem, Image sequence, Curvelets

Understanding and forecasting the evolution of geophysical fluids is a major scientific and societal challenge. Forecasting algorithms should take into account all the available informations on the considered dynamical system. The Variational Data Assimilation (VDA) technique combines in a consistent way all these informations in an Optimality System in order to reconstruct the model inputs. VDA is currently used by the major meteorological centres. During the last two decades about thirty satellites were launched to improve the knowledge of the atmosphere and of the oceans. They continuously provide a huge amount of data that are still underused by numerical forecast systems. In particular, the dynamical evolution of some meteorological or oceanic features (such as eddies, fronts, . . . ) that a human vision may easily detect is not optimally taken into account in realistic applications of VDA. Image Assimilation in VDA framework can be performed using pseudo-observation techniques : they provide some apparent velocity fields which are assimilated as classical observations. These measurements are obtained by some external procedures which are decoupled with the considered dynamical system. In this paper, we suggest a more consistent approach which directly incorporates image sequences into the Optimality System.

Assimilation de séquence d'image dans les modèles numériques

Résumé :

Comprendre et prévoir l'évolution des fluides géophysiques est un grand challenge scientifique à fort impact societal. Les méthodes de prévision doivent prendre en compte toutes les informations disponibles sur le système dynamique considéré. Les techniques d'Assimilation Variationnelle de Données (AVD) combines de façon cohérente toutes ces informations au sein d'un système d'optimalité afin de reconstruire les entrées du modèle. l'AVD est actuellement utilisée par les grands centres de prévision météorologiques . Durant les deux dernières décennies une trentaine de satellites ont été lancés afin d'améliorer notre connaissance de l'atmosphère et de l'océan. Ils fournissent de façon continue une grande quantité de données qui est toujours sous-utilisée par les systèmes de prévision numériques. En particulier l'évolution dynamique de certains "objets" atmosphériques ou océaniques (tourbillons, fronts, . . . ) que l'oeil humain peut détecter facilement, n'est pas prise en compte de façon optimale en AVD. L'assimilatio de séquence d'image en AVD peut être effectuée en utilisant une approche de type pseudo-observation : elles fournissent des champs de vitesse apparente qui peuvent être assimilé comme des observations classiques. Ces pseudo-mesures sont obtenues par des procédures externe qui sont découplées du système dynamique considéré. Dans ce rapport, on propose un approche plus cohérente qui incorpore les séquences d'image directement dans le système d'optimalité.

Mots-clés : Assimilation de données, problèmes inverse, séquences d'image, courbelettes 1 Introduction: images and assimilation Understanding and forecasting the evolution of geophysical fluids (ocean, atmosphere, continental water) is a major scientific and societal challenge with important applications (prediction of extreme meteorological events, estimation of climate change, droughts and floods forecast, etc.). It is currently the subject of an intensive research effort by the international scientific community. In order to achieve this goal scientists need to take into account all the available informations on the studied dynamical system. These informations can take several forms:

• Numerical models: they attempt to simulate the system dynamics on a wide range of time and space scales. Based on the Navier-Stokes equations, they are complicated because they must deal with the irregular shape of the domain, with the parameterisation of unresolved processes, and with highly nonlinear processes at the mesoscale and below. Numerical models represents mathematical and physical informations under the form of a set of coupled sets of non linear PDEs.

• Observations: By essence models simulate simplified real physical processes by a limited number of variables. Moreover, despite their regular improvement, some model parameters are poorly known (e.g. initial and boundary conditions). This is why it is important to also observe the system. Observations are now increasingly numerous as result of both satellite and in situ observing systems. However, their accuracy is not always satisfactory, and the processing of such a large quantity of heterogeneous data is difficult. Moreover, observations provide only a partial view of reality, localised in time and space, and some of them are only indirectly related to model variables.

• Statistical informations: On the first hand, there is an a priori knowledge about the error of the two above sources of information: some physical processes are simplified or even neglected and errors are introduced through numerical treatments (discretisation, computing errors, etc.). Most of the time there are known deficiencies of the model and therefore information about the structure of the needed correction needed is available. On the other hand observation errors are mainly driven by so called representativity errors i.e observation can represent phenomena that are not simulated by the model (subgrid processes). Most of the time the statistics of the representativity error can be estimated. All these informations are important and should be taken into account.

• Qualitative informations: forecasters use their knowledge on the system and additional materials such as satellite images to interpret the output of numerical prediction systems and provide corrected forecasts.

Since models and observations taken separately do not allow for a deterministic reconstruction of real geophysical flows, it is necessary to use simultaneously these heterogeneous but complementary sources of information through so called Data Assimilation methods. These are inverse modelling techniques based on the mathematical theory of statistical estimation. Their aim is to combine observations and prior model estimates, taking into account their respective statistical accuracies, in such a way that the combined estimate is more accurate than either source of information taken individually. In the early 80's Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations : Theoretical aspects[END_REF] proposed to use optimal control techniques in order to perform this task. This method, commonly referred as Variational Data Assimilation (VDA), has now been adopted by most of the main meteorological operational centres (European Centre for Medium Range Weather Forecast, Japan Meteorological Agency, Météo France . . . ). In variational data assimilation, the analysis problem is defined by the minimisation of a cost function that measures the statistically weighted squared differences between the observational information and their model counterpart. The cost function is minimised with respect to the control variables and this is done iteratively using a gradient-type descent method.

Models are based on fluid dynamics equations and on thermodynamic laws. After discretisation in space (finite difference and spectral methods) the model is represented by a system of PDE's with about 500 millions of variables, and this number will increase further in the future. Data are from numerous sources (in situ, radiosondes, aircraft, geostationary and polar orbiting satellites, etc.) and they are heterogeneous in density, quality and nature.

During the last two decades many satellites have been launched for the observation of the Earth in order to improve the knowledge of the atmosphere and of the oceans. They provide, among other data, photographic images of the earth system. Image-type data could be used in a classical data assimilation scheme, using the grey levels as observations of radiances emitted or reflected by the system. But connecting model variables with this kind of measurements is often complex and finally pixel levels are not assimilated in practice. This approcah is also too numerically expensive. However, satellites provide sequences of images that show the dynamical evolution of some meteorological or oceanic "objects": fronts, clouds, eddies, vortices, etc. The dynamics that the human vision easely detects in this kind of image sequences has clearely a strong predictive potential. Unfortunately this information is not currently optimally used in conjunction with numerical models. Furthermore, in some other domains, such as astrophysics, images are the only source of observation.

The purpose of this paper is to present an introduction to image sequences assimilation where two classes of techniques can be considered:

• Pseudo-Observation: an apparent velocity field is estimated from an image sequence using some image processing techniques. This estimated field can then be used as pseudo-observations in a classical variational assimilation scheme.

• Direct Image Sequence Assimilation (DISA): image sequences are assimilated through an appropriate observation operator directly into the Optimality System. This can

INRIA

be done on a pixel by pixel basis, alternatively some features can be extracted from the image sequence and then assimilated.

Using image sequences together with numerical models may present several difficulties:

• The state variables of the numerical models (for instance: wind, temperature, pressure, humidity in atmosphere modelling and current velocities, temperature, surface elevation, salinity in ocean modelling) are not directly measured by satellites. What is observed is more or less linked to those variables by diagnostic equations which will have to be included in the assimilation procedure.

• The physical processes that are observed are not always taken into account in the model (e.g. local convection in atmospheric modelling, ocean colour in ocean modelling). This may happen for other kind of observation as well but are likely to be systematic in the case of image sequences.

• Images are bidimensionnal informations whereas physical processes of geophysical fluids are three dimensional. From this point of view a major difference between the ocean and the atmosphere comes from their radiative properties. In the first case images give informations about the ocean surface whereas in the second case images integrate all the radiative profile of the observed atmosphere column.

• As for human vision, we can presume that the pertinent informations coming from an image is mainly brought by its discontinuities or high gradients. Unfortunately, numerical models have a tendency to smooth these discontinuities out.

The paper is organised as follow: basic principles of variational data assimilation techniques are described in section 2. Section 3 briefly presents some mathematical considerations about images. Section 4 browses the different techniques currently available to assimilate -in a regular assimilation scheme -velocity fields derived from image sequences (pseudo-observations). A new method for direct variational image sequences assimilation is described in section 5: this forms the original part of this paper. Some preliminary results are presented in section 6.

Variational Data Assimilation

Principle of VDA

Let us briefly present the basic principle of Variational Data Assimilation (VDA). The state of the system is described by a state variable x depending on time and space ; we have x(t) ∈ X . It represents the variables of the model (e.g. velocity, temperature, elevation of the free surface, salinity, concentration of biological or chemical species. . . ). The evolution of the state variable x is governed by the differential system:

dx dt = M(x, U ), t ∈ [0, T ], x(0) = V, (1)
where U ∈ P is some unknown parameters of the model: boundary conditions, model errors, parametrisation of subgrid effects. . . U may depend on space and time. The initial condition V ∈ X is unknown and depends on space. Let us assume that U and V being given, the problem (1) has a unique solution x(U, V, t), t ∈ [0, T ]. y(t) ∈ O are given observations of the system available between 0 and T . For the sake of simplicity of this description we suppose the continuity in time of the observations and that U depends on space only. The discrepancy between the observations and the state variable intervenes in the following cost function:

J(U, V ) = 1 2 T 0 H[x(U, V, t)] -y(t) 2 O dt + 1 2 U -U 0 2 P + 1 2 V -V 0 2 X , (2)
where H is the so-called Observation operator: it is a mapping from the space of the state variable toward the space of observations O where the comparison is carried out (for more details see section 2.2). The second and the third terms are regularisation terms in Tikhonov's sense [START_REF] Tikhonov | Regularization of ill-posed problems[END_REF]. They also allow to introduce some a priori informations. It is important to point out that J involves three different norms that define the topological structure of the functional spaces X , P and O; in practice norms come from an inner product and can take into account the statistical information by introducing the error covariance matrices. For instance a discretised in space version of the cost function (2) may involve

V -V 0 2 X := (V -V 0 ) * B -1 (V -V 0 ) ( 3
)
where B is the background error covariance matrix.

The problem of VDA can be considered as the determination of the analysed state (U a ,V a ) minimising J, that is verifying (U a , V a ) = argmin J(U, V ).

(4)

A necessary condition for (U a , V a ) to be a solution of (4) is that

∇J(U a , V a ) = 0 (5)
which becomes a sufficient condition if J is convex. As a first approximation one has to solve a problem of unconstrained optimisation despite the fact that, for instance, some quantities are known to be non-negative (such as humidity, salinity or concentration). From a practical and numerical point of view (U a , V a ) is estimated by a descent-type algorithm i.e by finding the limit of the sequence

U k+1 V k+1 = U k V k + λ k D k , (6)
INRIA where λ k is the step size realizing the minimum of J along the direction of descent D k deduced from the gradient of J. For computing the gradient ∇J with respect to the control variables (U, V ), one introduces the adjoint variable p as the solution of the adjoint model:

   dp dt + ∂M ∂x * • p = ∂H ∂x * .(H(x) -y), t ∈ [T, 0], p(T ) = 0. ( 7
)
After a backward integration of the adjoint model, the gradient of J is given by:

∇J = ∇ U J ∇ V J =   - T 0 ∂M ∂U * .p -p(0)   . ( 8
)
The derivation of the system (7) can be found in Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations : Theoretical aspects[END_REF]. The model ( 1) plus the adjoint model ( 7) along with the cost function (2) form the Optimality System (O.S.). The adjoint variable has the same dimension as the state variable of the direct model.

The observation space and operator

Measures are partial -in space and time -estimations of some continuous physical quantities. Moreover observations that are actually used in data assimilation schemes are in general a sampling of filtered measures. The space of observations O is a subset of R p where p is the number of observations. Observations y and state variables x are linked through the observation operator

H : X → O x(t) → y(t) = H(x(t)) (9)
which maps the state variable space X onto the observation space O. The observation space is equipped with an inner product whose definition takes into account some uncertainties on the measurements and some representativeness errors. We have

H(x) -y 2 O = (H(x) -y) * R -1 (H(x) -y), ( 10
)
where R is the observation error covariance matrix.

We talk about direct observations when the observed variables belong to the set of model state variables, and about indirect observations otherwise. In the first case, H may define a projection or an interpolation from X onto O. In the second case, evaluating H(x) may be difficult from a computational and mathematical point of view because y and x may be connected through complex physical laws. Observation operator should be regular enough in order to ensures that the partial derivative ∂H ∂x involved in the adjoint model ( 7) exists. In the case of direct observations, this property is generally satisfied which is could be far from the case for indirect observations. For instance, let us take the case of an ocean model with state variables velocity current, surface elevation, temperature and salinity. Direct observations could be some in situ measures of temperature and salinity in some specific locations of the studied domain. Indirect observations could be measures of surface salinity inverted from microwave radiances coming from remote sensing captors. In that case H should represent the relation between these microwave radiances and surface salinity.

Interpretation levels and mathematical definitions of an image

Our goal is to define a methodology based on the theoretical background of variational data assimilation (see section 2) to combine in a consistent way the informations coming from image sequences, and the ones coming from numerical model, a-priori knowledges and classical observations. A first step is to find a suitable level of interpretation of the image and a suitable mathematical framework to process it: this problematic is often called Image modelling.

Numerical raw images

The acquisition device converts the radiations coming from the observed scene into an electric signal which is then digitalised to form the raw image.

The device contains a sensor that is composed of an array of photo-sites that accumulate photons. A numerical raw image is an array of pixels whose values represent the quantity of photons the corresponding photo-site received during the acquisition procedure. Pixels of raw images can then be considered as radiance measurements. Pixels of black and white images are described by one single channel called grey level whereas pixels of coloured images are described with three channels (Red, Green, Blue). The intensity of a pixel is usually described by a number between 0 (black) and 1 (white) which is stored in a binary format by a string of digits (bits) whose length, the pixel depth, depends on the acquisition device and digitalisation procedure. For instance, METEOSAT images has a 8 bits pixel depth (256 grey levels) as it is the case for common grey-scale images. Some other remote sensing devices manage 10 or 12 bits pixel depth images. Note that computers usually store single precision floating point numbers with 32 bits (64 bits in double precision).

If we consider that a pixel intensity is a measurement, then it is represented in a lower precision than single precision floating point numbers. Finally, the spatial resolution of a satellite image refers to the quantity of earth's surface a single pixel covers.

However, from a mathematical point of view, an image is commonly considered as a real valued function of two continuous real variables. This allows the use of powerful math-INRIA ematical tools developed by the image processing scientific community. A numerical image is the discrete version of a mathematical one that is finally and practically manipulated by computers.

Levels of interpretations

At the pixel level an image of a fluid in motion is a digital instantaneous expression of the state of the fluid through the radiation that it reflects (visible light) or emits (infrared light). It is then a snapshot of an apparent state which may involve other physical quantities than the model simulates, like tracers (e.g. nebulosity and aerosols in the atmosphere, chlorophyll in the ocean). The radiances are generally not part of the studied system and the dependence between this radiation and the state variables is often non trivial. For instance, nebulosity is an intricate function of the state variable of meteorological models: clouds reflectance and intrinsic infrared emission depend both on thermodynamics (temperature and humidity) and on the micro-physical processes that occur inside the clouds (involving water, ice, snow and the size of the particles). Satellite images that come from infrared sensors (see Fig 1 right) give some indirect informations about the temperature of the observed system. At this level, information contained in images are not structured enough and represent a huge amount of data: typical satellite images are about 5000 × 5000 size that is 25 millions of pixels. This number should be multiplied by the number of frames in the observed sequence. For this reason, this level will probably not be efficient in realistic image sequences assimilation. Note that low level vision problematics (denoising, deblurring, inpainting, segmentation, etc.) are usually performed at the pixel level.

At the analysis level, an image is symbolically described by its contents, typically by the objects and their characterisation (geometry, layout, etc.) that appear in the observed scene. This kind of interpretation is used for high level vision problematics as motion estimation, edge detection, pattern recognition, etc. The advantage of this level is that it may less depend on the effective radiance emission as suggests Figure 1. In these images, visible channel (left) corresponds to a chlorophyll concentration map whereas the infrared channel (right) corresponds to a surface temperature map. But in both of them, one may detect the same structures as eddies and filaments. Currently, going from pixel level to analysis level in fluid flow images processing is a hard and still challenging task but this is beyond the topic of this paper. In the following we make use of structure as a generic word for the different features that may compose a geophysical image at the analysis level of interpretation.

Mathematical definition There exists several mathematical definitions for images and it is not the aim of this section to describe them all. They depend essentially on the way the image need to be processed. One can cite the multi-scale decomposition approach which comes from signal processing theory (wavelet [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and the more recent curvelet [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF]) transforms). Another approach uses Partial Differential Equation framework [START_REF] Aubert | Mathematical problems in image processing[END_REF] where the image is considered as a function of two variables verifying a PDE equation depending on the considered image processing problem. Some applications use stochastic approach as well to define images. For a good review of image analysis and the underlying mathematical definition see [START_REF] Chan | Image processing and analysis[END_REF]. At the end of this paper, we present an example where the anisotropic multi-scale decomposition approach by curvelet transform is used.

In the VDA framework, the level of interpretation together with the mathematical definition of the image may be represented by the image-type observation space and its corresponding observation operator (see section 5).

Sequence of images With a static image, one has only access to Eulerian informations. Sequences of images provide dynamical informations which can be interpreted as Lagrangian observations. The use of this type of information in variational data assimilation schemes was recently studied for ocean and hydrology problems (see [START_REF] Nodet | Modélisation mathématique et assimilation de données lagrangiennes pour l'océanographie[END_REF] and [START_REF] Honnorat | Assimilation de données Lagrangiennes pour la simulation numérique en hydraulique fluviale[END_REF]). Fig 2 is an image sequence taken over the Atlantic Ocean. On this example, one can see the evolution of a front over western Europe and a cyclonic depression over the ocean. This illustrates the human vision capability for estimating the evolution of the atmosphere through the features that it is available to detect and track. Because it takes into account all the observations at once in the assimilation window, the VDA framework seems to be naturally well adapted to consistently process the dynamical information the image sequences provide. Processing and analysing image sequences may be done by extending the concepts listed above for static images using two different ap-INRIA proaches: first, the sequence may be assimilated as a set of single two dimensional images (i.e. the frames) which are processed and analysed individually. An other approach is to consider the sequence as a three dimensional data : two spatial dimensions and time: the use of some extension of image processing techniques to the 3D case may reinforce the temporal consistency of the direct image sequences assimilation. Sequence of images are widely available in meteorology, oceanography, hydrology, astrophysics, and even medicine. Historically they have been mainly used for a "by eye" analysis from experts in the field. However since the mid 80' researches have been carried out to derive velocity fields from the sequences, with applications mainly for fluid dynamics (and very recently for movie compression and medical applications). One can use the velocity field derived from the image processing techniques as pseudo-observation of the studied system, for instance in a regular VDA scheme [START_REF] Herlin | Coupling models and data : which possibilities for remotely-sensed images[END_REF].

There are several ways to extract a velocity field from a sequence of image. They can be divided into two categories: the frame-to-frame estimators and the so-called Image Model (see section 4.2).

Frame-to-frame Motion Estimator

The first class of Motion Estimator, directly based on image processing techniques, aims at estimating the velocity field that would 'transport' one frame to another.

Among them one can cite the statistical methods (mostly based on correlations between 2 images). They are mainly used for PIV (Particle Image Velocities) and PTV (Particle Tracking Velocities) experiments i.e. for lab experiments in fluid mechanics where particles are added to the fluid in order to follow the flow [START_REF] Adrian | Particle imaging techniques for experimental fluid mechanics[END_REF]. They are easy to implement, but they can be quite expensive when the number of particles increases. Applications of such methods are not really in the scope of what we are studying here.

Variational methods are more adapted to assimilation of image sequences in geophysics. The most common one is usually called optical-flow; This classical approach in computer vision is based on the conservation of the global luminance between two images [START_REF] Horn | Determining optical flow[END_REF]: let I : Ω × R → R be the luminance of the pixel, the optical flow is the vector v(x, y) verifying

dI dt = ∂I ∂t + ∇I • v = 0. ( 11
)
This is not the only law of conservation which could be considered according to the nature of the image: with an image of the colour of the ocean an equation of conservation of chlorophyll (with source and sink) terms could be considered, with an image of Sea Surface Temperature (SST), the Boussinesq approximation could be used. . . Problem ( 11) is an ill-posed problem therefore one need to add a regularisation constraint; the choice of the regularisation term depends on the application field [START_REF] Vigan | Sea surface velocity from sea surface temperature image sequences[END_REF] for oceanography, [START_REF] Amodei | A vector spline approximation[END_REF] for wind field). This problem can be summarised by finding v ⋆ that minimises F defined by

F (v) = 1 2 Ω ∂I ∂t + ∇I • v 2 dxdy + 1 2 λ S(v) 2 , (12)
where S : R 2 → R 2 is a spatial regularisation function (see [START_REF] Auroux | Identification of velocity fields for geophysical fluids from a sequence of images[END_REF]) and λ ∈ R the regularisation factor. The minimisation problem ( 12) is treated either by optimal control or by vector spline [START_REF] Amodei | Étude d'une classe de fonctions splines vectorielles en vue de l'approximation d'un champ de vitesse[END_REF], [START_REF] Suter | Motion estimation and vector splines[END_REF], [START_REF] Isambert | Fast and stable vector spline method for fluid apparent motion estimation[END_REF]).

Image Model

Recently [START_REF] Herlin | Estimation of a motien field on satellite images from a simplified ocean circulation model[END_REF], [START_REF] Huot | Data assimilation of satellite images within an oceanographic circulation model[END_REF] and [START_REF] Korotaev | Analysis of the black sea surface currents retrieved from space imagery[END_REF] proposed an extension of the optical flow methods called "Image Model". This is a more advanced INRIA motion estimator using techniques based on VDA that allows the use of several successive frames at once by adding a governing equation for v to the optimality system that generates the pseudo observations. The aim is to retrieve a 'continuous' in time motion vector field from a sequence of image I obs (t), t ∈ [0, T ]. As in 4.1 it assumes the conservation of light and the image pixels are transported from one image to the other by the velocity field to be estimated. In order to improve the temporal consistency of the estimated velocities, this method requires a simple temporal evolution model N , along with the transport of the pixels it forms the so-called Image Model:

         ∂I ∂t + ∇I • v = 0, ∂v ∂t = N (v), I(0) = I 0 , v(0) = v 0 . ( 13
)
The pixels value of the images I obs is then assimilated in the previous model using variational data assimilation (see 2.1) where the control vector is the initial condition v 0 of N .

F (v 0 ) = 1 2 T 0 I -I obs 2 dt + λ 2 S(v 0 ) 2 . ( 14
)
The advantage of this method compared to the one described in 4.1 is that it takes into account all the images of the sequence at once and therefore can cope with missing data on some of the frames.

Limitations of the pseudo-observation technique

The velocity field produced by the motion estimation methods proposed in the two previous sections could be used as pseudo observation in a classical VDA scheme such as presented in section 2 [START_REF] Herlin | Coupling models and data : which possibilities for remotely-sensed images[END_REF]. However, despite their relative simplicity for implementation and their rapidity, pseudo-observations techniques based on motion estimation may suffer from some limitations. First, frame to frame motion estimator does not take into account a sequence of more than two images. Due to the lack of consistency in time, it cannot capture the dynamical evolution of the image during the entire assimilation window. One of the consequences is that it cannot deal with missing data which can be quite frequent in some of the applications targeted here. They can be due to a glitch in the observing system or an obstacle (e.g. clouds).

Image Model technique gives an answer to these two problems. However the model N of the temporal evolution of the velocity only depends on v and therefore is likely to be a somewhat poor representation of the real behaviour of the velocity.

Neither of these two techniques takes into account physical informations about the underlying processes observed in the images, moreover some assumptions like the conservation of the luminance are clearly not valid for fluid flow images showing break up structures. Finally, velocity fields obtained through the pseudo-observation approach are apparent velocities that can be quite different from the actual velocities. For instance a tracer absence in a subset of an image sequence will lead to null apparent velocities in that area. Such techniques then need to provide error maps in order to account for this kind of error in the interpretation of the observation. Therefore the definition of the observation operators and corresponding observation error matrices required by the assimilation of such data is likely to be a very difficult task.

Direct Assimilation of Image Sequences

General Framework

In the pseudo observation approach, the information content is extracted from the images externally from the Optimality System equations (1, 2, 7). An other way to proceed is a direct assimilation of the image sequence in the O.S. This is more consistent with the VDA philosophy, in particular the information content of the images is extracted consistently with the physics of the model. This approach intrinsically avoids the problem of observation interpretation error mentioned above in the absence of tracer.

The principle is to add to the cost function a term directly linked to the images and their dynamics. Let us denote by M t (U, V ) the state variable x(t) when initial condition is x(0) = V and the model parameters are U . Then J would take the form:

J(U, V ) = T 0 y -H[M t (U, V )] 2 O dt Classical term Jo + T 0 f (t) -H X →F [M t (U, V )] Model to Image Operator 2 F dt. (15)
In this expression f (t) ∈ F is a frame of an image sequence of the observed dynamical system at time t. It belongs to the space of images F equipped with a norm • F . The operator H X →F is a mapping from the space of the state variable toward the space of images: in a sense, it generates an image of the dynamical system at time t from the values of its state variable at this date: in the following H X →F is called Synthetic Image Operator. Note that by the sake of simplicity, we omit in (15) the regularisation term which appears in the expression (2) of the classical VDA cost function. Synthetic image operator works at the pixel level interpretation of the image (see section 3). [START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF] showed on a simple case that assimilating directly pixel levels gives a better result than using an optical flow-based pseudo-observation approach. However, due to the size of satellite images (typically 5000 × 5000) such formalism would require the processing of a huge amount of data and would prove not feasible for realistic applications in geophysics.

INRIA

At the analysis level (see section 3), an image is characterised by coherent structures (eddies, fronts,. . . ). If we suppose that these structures can be approximated in a functional space S equipped with a suitable norm • S , then the comparison between image-type observation and model output can then be done in that space and the cost function becomes: With the usual variational approach described in section 2 a term in the form:

J(U, V ) = T 0 y -H[M t (U, V )] 2 O dt + T 0 H F →S [f ] Image to Structure Operator -H X →S [M t (U, V )]
∂H X →S ∂x * . (H X →S [x] -H F →S [f ]) (17)
is added in the right hand side of the adjoint model ( 7). In practical applications, we will then prefer working with differentiable Model to Structure operators. Note finally that one may consider assimilating several image sequences by adding corresponding terms in the DASI cost function (16).

Image to Structure operator H F →S

As we previously said, the Image to Structure operator represents the level of analysis we consider to manipulate images. It could represent some transformations we apply to the image. The image space (in the mathematical sense) of this transformation is a space of another data type. The constraints which are imposed to choose S is that it should have a structure of normed space in order to use simple rules for the differentiation of J. Moreover, its dimensionality must remain reasonable to store the extracted data. It may exist many different approaches to define the space S. For instance, one can consider 1. Frequency characteristics: in that case S is the image space of a multi-scale transformation such as wavelet and curvelet transforms. This approach is interesting for its well known capability in data compression, denoising and edge extraction for example.

2. Geometric characteristics: some image analysis techniques can detect some geometrical feature in the image (e.g. by the mean of active contours). This features may take the form of a parametrised curve in R 2 which could be assimilated as Lagrangian observation. This approach will necessitate some high level image analysis.

3. Qualitative characteristics: in geophysics there are several recognisable structures connected to physical phenomena (e.g. cyclone, front, etc.). One could define S as a dictionary of such couples observable structure / physical phenomenon.

Model to Structure operator H X →S

Features in geophysical images as eddies or filaments may correspond to some tracers in the observed system. Tracers can be defined as any quantities that are transported by the fluid flow. For instance, potential vorticity, temperature, salinity can be considered as atmospheric and ocean tracers. Modelling and simulate their evolution from the model state variables may help to define the Model to Structure operator. Indeed, one may use their description maps (e.g. concentration maps) as a Synthetic Image Sequence from which structures can be extracted using the Image to Structure operator. More precise correspondence between tracer description and measured radiances (i.e. pixel intensity) may also be considered by modelling the acquisition procedure. For instance, this may connect properly the three dimensional informations of the physical model with the two dimensional informations of an image. However, as suggests Figure 1, we may think that such precise modelling is not always useful. Using a tracer evolution modelling leads to a Model to Structure operator of the form

H X →S = H F →S • H X →F , (18)
where H X →F is the Synthetic Image Sequence operator. Note that if tracers take part in the model the Synthetic Image operator is just some kind of projection (to be composed with an image acquisition model if necessary). However some features of interest in geophysical fluid images do not correspond to any tracers (such as waves or dry intrusions in atmosphere images). In that case, other techniques for extracting structures from model outputs avoiding synthetic image sequence production should be developed.

6 Numerical example

Experimental framework

Laboratory experiment can be a good testbed for the new technique suggested in 5. Indeed it offers the complexity of true physical processes while keeping it manageable, and more importantly it opens the way to realistic applications with images.

As a first step to study the numerical feasibility of this new approach, one will focus on so called twin-experiments where observations do not come from real world but are INRIA generated by the numerical model. In this configuration, the "true" state is known and this allows us to evaluate and understand the behaviour of a new method.

The experimental framework is that of the study of a drift of a vortex on a turntable. The evolution of a vortex in the atmosphere is simulated at the CORIOLIS experimental turntable (Grenoble, France) which re-create the effect of the Coriolis force on a thin layer of water. A complete rotation of the tank takes 60 seconds which corresponds to one Earth rotation. The vortex is created by stirring the water and made visible thanks to the addition of a passive tracer (fluorescein). Photographs of the vortex are taken from above the turntable which constitute the observed image sequence. For more details on the statements ant the motivations of this experiment see [START_REF] Flór | Dynamcis of monopolar vortices on a topographic beta-plane[END_REF].

In this configuration, the evolution of the fluid can be modelled by the state variable x = (u, v, h) whose components verify the Shallow-Water equations

   ∂ t u -(f + ζ)v + ∂ x B = -r * u + ν∆u, ∂ t v + (f + ζ)u + ∂ y B = -r * v + ν∆v, ∂ t h + ∂ x (hu) + ∂ y (hv) = 0. (19)
Unknowns are the zonal component u(t, x, y) and meridional component v(t, x, y) of the current velocity and the water elevation h(t, x, y). They depend on time t and the two horizontal directions x and y. One defines the relative vorticity = ∂ x v -∂ y u and the Bernoulli's potential

B = g * h + 1 2 (u 2 + v 2 )
, where g * is the reduced gravity. The Coriolis parameter on the β-plane is given by f = f 0 + βy, ν is the diffusion coefficient and r * the bottom friction coefficient. In this paper the following numerical values are used for the parameters: r * = 0.9, 10 -7 , ν = 0, f 0 = 0.25, g * = 0.02 and β = 0.0406.

The simulation is performed on a rectangular domain Ω =]0, L[×]0, H[ representing a sub-domain of the turntable with L = H = 3m. The domain is discretised on a 128 × 128 uniform Arakawa C-type square grid. A finite differences scheme is used for space discretisation. Time integration is performed using Leap-Frog scheme. Time step is set to 0.05s of the turnable experiment which may correspond to 72s in the atmosphere. Image sequence acquisition is simulated using synthetic images whose grey levels correspond to the concentration of a passive tracer advected by the velocity field (see section 6.2.1).

Assimilation procedure

In our test, we consider the problem of recovering the initial state of the fluid V = x 0 = (u, v, h)(0, •, •) which constitute our control variable. The observations are images and we do not assimilate any measurements of some components of the state variable. We use an Images to Structure operator of the type "frequency characteristics" by the mean of a multi-scale and multi-orientation transformation of the images: the only source of information we assimilate is a subset of the coefficients of the decomposition (see below for more details). We consider an assimilation window of 100s (2000 time steps) which corresponds approximatively to 40 hours in the Earth atmosphere. Images are captured every 5s (100 time steps) which then simulates a real capture frequency of one image per two hours. Observations then consist of a set of 20 images of 128 × 128 resolution.

The application of VDA to the above test case using section 5 formalism relies on the minimisation of the following cost function:

J(x 0 ) = α o T 0 H F →S [f ] -H X →S [M t (x 0 )] 2 S dt Jo +α b J b x 0 -x b 2 X , (20)
where observation and regularisation terms J o and J b are balanced by the mean of two adjustable coefficients whose ratio α b /α o may represent the confidence in the background x b . The following subsections are devoted to a detailed description of these terms, for this experiment.

Observation term J o

Synthetic and true images generation As mentioned in 5.3, one way to define the Model to Structure operator H X →S is to first create an image from the model outputs using a Synthetic Images Operator H X →F and then apply the Image to Structure Operator H F →S . In the current experimental framework H X →F is obtained by simulating the advection of the monochromatic passive tracer into the fluid. To that purpose, we define

H X →F [u, v, h] = q, (21)
where q = q(t, x, y) is the passive tracer concentration verifying the following conservation equation posed in Ω and for t ≥ 0

∂ t q + u∂ x q + v∂ y q = 0, (22)
where (u, v) is given by ( 19). One obtains the synthetic image at a given date t using the concentration map (x, y) → q(t, x, y) as grey levels. True images are generated from a given true state (u t , v t , h t ) using the same technique, since we are in a twin experiments framework.

Image to Structure operator: Curvelet transform and Thresholding As mentioned in 5.2, one way to extract consistent structures from the images is to work in the frequency domain by applying a multi-scale transformation. In this paper we use a curvelet transform [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF]): An image f can be decomposed as

f = j,k,l ϕ j,k,l , f ϕ j,k,l , (23)
INRIA where (ϕ j,k,l ) j,kl are the elements of the curvelet frame. The inner product ., . is the usual L 2 (R 2 ) one. Curvelet transform is a linear isometry operator (see Appendix) and then its adjoint needed in ( 7) is given by the inverse transformation (reconstruction). This is a convenient feature for practical use in the VDA framework. For more details on curvelet transform see the Appendix. In order to extract the structures from the images, we need to apply a threshold on the curvelet coefficients ϕ j,k,l , f . We then define the Image to Structure operator by the threshold of some coefficients of the curvelet transform of the image:

H F →S [f ] = T (DCT(f )), (24)
where T denotes a threshold operator applied on the curvelet coefficients of the image f computed by the Discrete Curvelet Transform (DCT) algorithm (see Appendix). By the sake of simplicity, we use in our experiments the following hard threshold function

τ (x) = x if |x| ≥ σ, 0 if |x| < σ, (25)
where x represents a coefficient in the curvelet series ( 23) and σ > 0 is a given threshold value. The space S is then the image space of the curvelet transform. We assimilate a fixed subset of the curvelet coefficients of the image decomposition. In our experiments, we keep only 10 percents of these coefficients. Using hard thresholding is more or less equivalent to assimilate images at the pixel level with this difference that it uses the efficient compression ability of this transformation (see Appendix): this feature is critical for future realistic applications. A large choice of threshold function is available depending on the applications. Some classical threshold functions are listed in Antoniadis (2007). More sophisticated thresholding using a rule of Total Variation (TV) minimisation can reduce the pseudo-Gibbs and element-like artifacts that appears with classical threshold functions [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF], and may be used in this framework. DASI may give a starting point to very interesting investigations in the area of the application of curvelet and other multi-resolution transformations [START_REF] Ma | Spatiotemporal structure extraction and denoising of geophysical fluid image sequences using 3d curvelet transforms[END_REF].

Model to Structure operator Finally our Model to Structure operator is the composition of the previous operator (24) with the Synthetic Image operator ( 21):

H X →S = H F →S • H X →F = T • DCT. (26)

Regularisation term J b

The regularisation term

J b (x 0 ) = x 0 -x b 2 X (27)
involves the background state x b and some nondimensionnalisation of the equations through the definition of the norm • X . This norm is also used to define a pre-conditioning of the minimisation problem.

Background term In our experiments, the background is set to the constant value x b = (u b , v b , h b ) = (0, 0, h m ) where h m is the mean water elevation in the tank. This approximatively corresponds to a fluid at rest when the turntable speed is stationary. We then set the ratio α b /α o to a small value (about 10 -6 ).

Preconditioning The control state variable space X is equipped with the following norm

x X = x T B -1 x ( 28
)
where B is the so-called background error co-variances matrix. Its diagonal coefficients represent the standard deviations error on each variable with respect to itself: they are a-priori estimated from model outputs. The velocity field is regularised using a method based on the Gradient Vector Flow technique [START_REF] Souopgui | Vector field regularization by gradient vector flow[END_REF]. The off-diagonal blocs are the cross-covariances between variables and come from balance relations between the observed variables and the others: this is a classical technique in VDA. Our observation operator connects images with the velocity field and that is why we decided to balance the water elevation h with respect to (u, v) by "inverting" the geostrophic equilibrium relation

u = -g f ∂h ∂y , v = g f ∂h ∂x
. Finally, we use the classical change of variable x 0 ′ = B -1/2 x 0 for preconditionning. For more details about this strategy see [START_REF] Vidard | Variational data analysis with control of the forecast bias[END_REF].

Optimisation

The cost function J is successively minimised on assimilation sub-windows of increasing sizes using the strategy described in [START_REF] Luong | A variational method for the resolution of a data assimilation problem in oceanography[END_REF]. This technique avoids some stagnation of the optimisation algorithm around a local minima which may occur when J is not convex. The minimisation procedure is initialised with the previous background. In each successive minimisation, we use the N1QN3 algorithm of the MODULOPT library [START_REF] Gilbert | The modules m1qn3 and n1qn3. version 2.0[END_REF].

Numerical Results

Following the experimental set-up described above, one try to reconstruct the true initial state assimilating only the image sequence and starting from a background at rest. The true state represents the evolution over the time window of a single vortex subject to the Coriolis force. A passive tracer is added into the vortex and its concentration is used to generate the images, only one snapshot every 100 time steps is retained as observation. Figure 3(left) shows a subset of the assimilated image sequence (i.e. the evolution of the passive tracer). The background state being the system at rest it does not allow for the transport of the tracer, as shown in Figure 3 and initial water elevation (bottom). Background state variables correspond to a system at rest.

INRIA

that the estimated vortex is not sufficient to properly generate the images at the end of the sequence. Figure 4 shows the true velocity field (top left) that was used to generate the observation (left sequence of images of Figure 3) and the estimated velocity field (top right) that is retrieved thanks to the minimisation of J. It also shows the true water elevation (bottom left) and the estimated elevation retrieved at the same time (bottom right). The estimated vortex is at the right position, but its intensity is a bit weak compared to the true one, this explain why the reconstructed images at the end of the sequence of figure 3(right) are not perfect. Note that, in this example, the background x b (i.e. the starting point of the minimisation) was taken to be the system at rest. Additionally, the assimilated image sequence (figure 3(left)) does not show very contrasted images, these are smoother than real ones. Moreover, in realistic applications, one would be able to assimilate other kind of data. This is not a very favourable case and therefore it is encouraging to be able to recreate a vortex with about the same characteristics. A more realistic background, more contrasted images and/or a better threshold, should improve the efficiency of this test-case.

Conclusions

Predicting the evolution of natural dynamical systems requires to link together heterogeneous sources of information. There is a strong social demand for the improvement of the prediction of geophysical fluids evolution and VDA is a powerful tool to achieve this goal. The observation of Earth by satellites is very expensive but unfortunately dynamical information that image sequences contain are still underused by forecasting numerical systems. However they represent a large amount of structured data and it is then of great interest to suggest some way to process them in conjunction with numerical forecast systems. Some previous and still in development works could lead to the assimilation of data coming from image sequences into a VDA scheme: but this so-called pseudo-observation technique suffer from the lack of physical consistency with the observed system. The Direct Image Sequences Assimilation approach presented in this paper aims at processing the image sequences within the Optimality System of the VDA scheme: it then combines image-type information consistently with the underlying physical model. In practice, DASI is based on structures extraction from both image sequences and model outputs. One may use some already existing processing and analysing image techniques to perform the structures extraction. Features in geophysical images usually corresponds to some tracers presence in the observed system. A natural way to extract structure from the model output is then to simulate the evolution of these tracers. This approach builds a synthetic image sequence from which the structures are effectively extracted. Other techniques for extracting structures from model outputs avoiding synthetic image sequence production should be developed. Our first result shows that direct image sequence assimilation approach using structure extraction with anisotropic mutli-scale transformation of the images is achievable.

There is no doubt that investigating the way to optimally integrate image sequence in VDA framework will bring up some very interesting questions and gives some improvements in both areas of VDA and image analysis. For instance, because resolution of a satellite images is often higher than the spatial resolution of numerical models, we can expect that satellite image assimilation gives some improvements in the subgrid parametrisation problematic in VDA. DASI may also initiate some new definition of image interpretations and improve special features tracking in geophysical fluid image sequences.

Appendix: Curvelet transform

In 1999, an anisotropic geometric wavelet called ridgelet was proposed by [START_REF] Candès | Ridgelets: A key to higher-dimensional intermittency ?[END_REF]. It is optimal at representing straight-line singularities. However its application is limited to objects with global straight-line singularities. In order to take into account images with local line or curve singularities, the so-called first generation of Curvelet transform [START_REF] Candès | A surprisingly effective nonadaptive representation for objects with edges[END_REF] was developed using ridgelet transforms applied on sub-partitions of the image. A second generation of curvelet transform based on a frequency partition technique was introduced by the same authors in 2004 [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF]. A 2D discrete curvelet transform was implemented in 2006: the Matlab and C++ source codes are freely available for non-commercial research purposes at http://www.curvelet.org. It also exists a parallel version of this algorithm. The following succinct description of the continuous curvelet transform is mainly inspired from [START_REF] Candès | Fast discrete curvelet transforms[END_REF] and [START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF]. We first fix some notations: Let W : [1/2, 2] → R and V : [-1, 1] → R be respectively a radial and an angular windows -in the spatial domain -, supported on [1/2, 2] and [-1, 1] respectively and obeying the admissibility condition

+∞ l=-∞ V 2 (t -l) = 1, t ∈ R and +∞ j=-∞ W 2 (2 -j r) = 1, r > 0.
Examples of such functions are given in [START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF]. Note that curvelet transform efficiency few depends on the choice of these functions. As in wavelet theory, we define the low-pass window W 0 (r) = 1 -j≥0 W 2 (2 -j r), for all r > 0. Finally, let us denote by f the Fourier transform of an element f of L 2 (R 2 ): it is a function of the frequency variables ξ = (ξ 1 , ξ 2 ) whose polar coordinates are denoted by (r, θ). Curvelets are multi-scale, multi-orientation transformations with atoms ϕ j,k,l indexed by three integers j ≥ -1, k ∈ Z and l ∈ N which respectively stand for scale, location and orientation. At the coarsest scale (j = -1) curvelets ϕ -1,k are translations of a "mother" curvelet ϕ -1 :

ϕ -1,k (x) = ϕ -1 (x -k), x ∈ R 2 , k = Z 2 . ( 29
)
They are then isotropic. The mother curvelet ϕ -1 is the inverse Fourier transform of the above low-pass window W 0 . We have φ-1 (ξ) = W 0 (|ξ|).

(30)

INRIA

For each scale 2 -j , j ≥ 0, (finest scales), curvelets are defined through the inverse Fourier transform of compactly supported windows: we define a window U j in the frequency domain for all ξ = (ξ 1 , ξ 2 ) with polar coordinates (r, θ) by

U j (ξ) = 2 -3j/4 W (2 -j r)V 2 ⌊j/2⌋ θ 2π , (31)
where ⌊•⌋ denotes the integer part operator. The support of U j is a polar 'wedge' determined by supp(r → W (2 -j r)) = [2 j-1 , 2 j+1 ] and supp(θ → V (2 ⌊j/2⌋ θ)) = [-2 -⌊j/2⌋ , 2 -⌊j/2⌋ ].

Figure 5 (left) shows the support of the first window U 0 in the frequency domain (see also [START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF]). We consider a uniform discretisation of the angle coordinate θ ∈ [0, 2π[

θ l = 2π × l × 2 -⌊j/2⌋ , 0 ≤ l < 2 ⌊j/2⌋ , (32)
and the rotation by θ radians matrix

R θ = cos θ sin θ -sin θ cos θ . ( 33
)
Note that the number of discrete orientations θ l depends on the scale index j: this number is doubled when scale change from an odd index 2j +1 to the larger next even one 2j +2. It is unchanged when scale index is changing from 2j to 2j + 1 (see Figure 5 (right)). Finally, for j ≥ 0, curvelets are defined for all x ∈ R 2 by ϕ j,k,l (x) = ϕ j R θ (xx

(j,l) k ) , (34)
where the waveform ϕ j is the inverse Fourier transform of the frequency window U j (i.e. φj = U j ) and the translation vector x (j,l) k

, k = (k 1 , k 2 ) ∈ Z 2 is defined by

x (j,l) k = R -1 θ l (2 -j k 1 , 2 -j/2 k 2 ). ( 35
)
It is proved in [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] that the above family of curvelets defines a tight frame (see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] for a detailed definition of this term): for all f ∈ L 2 (R 2 ) we have (Parseval equality)

j,k,l | f, ϕ j,k,l | 2 = f 2 L 2 (R 2 ) , (36)
and the reconstruction formula

f = j,k,l ϕ j,k,l , f ϕ j,k,l (37)
where •, • is the usual inner product in L 2 (R 2 ): f, g = R 2 f g.

Remarks:

RR n°6701

• Curvelet coefficients in (37) are determined by c j,k,l := ϕ j,k,l , f = 1 (2π) 2 f (ξ)U j (R θ l ξ)e i x (j,l) k ,ξ dξ (38)

• In practice we use Discrete Curvelet Transform (DCT) which is designed for digital inputs f defined on a Cartesian grid {(n 1 , n 2 ) ∈ N, 0 ≤ n 1 , n 2 < n ∈ N}. Digital coefficients are defined

c D j,k,l = n 1 ,n 2 f (n 1 , n 2 )ϕ D j,k,l (n 1 , n 2 ) ( 39
)
where the construction of the digital curvelet waveforms ϕ D j,kl is an adaptation of the above continuous transform for Cartesian grids. This lead to a non uniform discretisation of the orientations. For more details and C++ implementation see [START_REF] Candès | Fast discrete curvelet transforms[END_REF].

• Let f F n (resp. f W n , f C n ) denotes the partial reconstruction of a C 2 image away from C 2 edges (i.e. containing curve-line singularities) obtained after a hard thresholding by selecting the n largest terms in the Fourier series (resp. in an efficient wavelet basis, in the above curvelet decomposition), then we have [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF]

f -f F n 2 2 = O(n -1/2 ), (40)
f -f W n 2 2 = O(n -1 ), (41)
f -f C n 2 2 = O((log n) 3 .n -2 ), (42)
where . 2 denotes the L 2 (R 2 ) norm. This can be explained by the oriented tight shape of the curvelets supports. Curvelet support tiling is then more efficient for curve-line singularities than wavelet supports do: this efficiency is illustrated in Figure 6. Right: partition of the frequency plane: the shaded area shows an example of a "wedge"on which a curvelet is supported in the frequency plane. At scale j, this wedge has a width of order 2 j and a length of order 2 j/2 . Figure 6: Reconstruction of a circle (solid) at a given scale. Curvelet supports have a tight shape and can be oriented. They fit better with curve-line singularities (left) than wavelet supports do (right). That is why curvelet transform is more efficient in terms of compression and edge enhancement for this kind of images than wavelet transform.

Figure 1 :

 1 Figure 1: MODIS satellite images of the Black Sea in visible (left) and infrared (right) channels. In both image, one detects easily the same structures as eddies and filaments even though they are not linked to the same physical tracer radiance: chlorophyll concentration (left) and surface temperature (right). (Courtesy : MHI)

Figure 2 :

 2 Figure 2: sequence of cloud coverage images above Europe from Meteosat (source Météo-France) from 28/04/2008 to 29/04/2008

Figure 3 :Figure 4 :

 34 Figure 3: True images (true tracer concentration) (left) and concentration of the tracer advected by the background (centre) and by the retrieved velocity field (right) at initial time (bottom), t=33s and t=66.5s (centre) and t=100s (top).

Figure 5 :

 5 Figure5: Left: support of the frequency domain window U 0 . Right: partition of the frequency plane: the shaded area shows an example of a "wedge"on which a curvelet is supported in the frequency plane. At scale j, this wedge has a width of order 2 j and a length of order 2 j/2 .

  F →S extracts structures from images (see next subsection): it represents the analysis level of the image. The operator H X →S extracts the same type of structures from model state variables (model outputs). The norm • S depends on the mathematical definition of the image at the level S. By choosing S = F and H F →S = f one go back to formulation (15).

	2 S dt (16)
	Model to
	Structure
	Operator
	where H

RR n°6701


Mengchen Wang

Nicolas Férey

Frédéric Magoulès

Patrick Bourdot

Interactive simulation for easy decision-making in fluid dynamics

Keywords:

published or not. The documents may come

Interactive simulation for easy decision-making in fluid dynamics

Introduction

Interactive simulations have been proved to be useful to solve many scientific questions. For instance, Fold'It [KDC * 11] helps users to solve the structure of an HIV protein by using an interactive molecular simulation approach. We believe that interactive simulations applied to fluid mechanics can be an interesting tool to support decision making, problem solving and optimization in this field, especially when humans are needed in the loop. Using Computational Fluid Dynamics (CFD) requires a large set of expertise and tools to design and perform a simulation, then to visualize, explore and analyze the simulation results. Different applications are used sequentially, including 3D modeling software, simulation and visualization software. This process is timeconsuming because it involves a lot of technical and unproductive tasks such as data processing and wrapping that require different tools. Moreover, to get accurate relevant physical results, realistic 3D fluid simulations need expansive computational resources. Therefore, it is crucial to detect early errors on the simulation model or the parameters used before the end of the computation. Consequently, it is necessary to regularly extract and analyze simulation snapshots. By taking advantage of advances in calculation, algorithms, and simulation methodology, an interactive simulation approach was proposed to address these issues, especially to reduce the time lost during this complex process. . These methods target extreme performance, high level rendering, and stability, but at the expense of the physical relevance, which is required by decision-making, physics education or research. Using more relevant methods requires computing centers that are often isolated from the visualization resources which are usually undersized. Even if a lot of technical results were published to support the interactive fluid simulation approach, this methodology remains rarely used. Moreover, there are only a few studies that deal with the problem of the usefulness of interactive simulations in terms of performance and user experience. We propose in this paper a work-in-progress to address this issue by designing an interactive fluid simulation platform based on Unity 3D and evaluating the benefit of the interactive simulation approach on decision making from fluid simulation on a simple but realistic use case.

A platform coupling Unity 3D with classical simulation tools

To provide interactive fluid simulation features, the visualization tool needs to be graphically synchronized with the initial model used as an input to the simulation software. Usually, the user input includes physical properties of the fluid (e.g., speed, pressure, temperature, fluid type), and boundary conditions (e.g., wall, floor, obstacles). During the simulation, the visualization application needs to get the state of the model sent at each time-step from the simulation to update the rendering. To take into account interactive modification on the fluid or its boundary conditions, the interactive visualization tool sends to the simulation the modified fluid states or boundary conditions. According to the rendering and simulation frame rates, it is possible to tune the frequency of data exchange.

To realize the communication between two systems, we chose to use a network protocol.

A dedicated network protocol and API for interactive fluid simulation

To be able to integrate any simulation tools daily used by experts, and to ensure the quality of the simulation results, we have developed an API called CFDriver [WFBM19] embedding a network protocol to perform the data communication between the simulation tool and visualization one. As a cross-platform API, it works with different programming languages and on different operating systems. It allows an expert to connect our interactive simulation platform to any simulation tools with only a few lines of code. This API provides a distributed architecture to enable the use of computational centers for large fluid simulations. This network protocol allows us to receive the simulation results and also sent the updates from the users to the simulation server during the simulation.

Unity 3D as a user interface for interactive simulation

For the interactive visualization tools, we decided to use Unity 3D for our platform. This software was initially designed for game programming but is a good compromise between speed of development and rendering performance. We provide classical camera manipulation tools for rotating, panning, and zooming. For the user's input, we chose to use ray casting to select and change the state of each voxel corresponding to a part of the simulation grid. We use the same cursor to perform multiple tasks such as setting a wall, filling a voxel with water, or emptying a voxel. Besides, it is possible to stop, start or restart, pause or continue the simulation after updating its content such as boundary conditions, grid cell states.

Lattice Boltzmann as simulation approach

We chose to use Lattice Boltzmann Method (LBM) for the simulation tool because it is a good compromise between performance and physical relevance. New solvers such as adaptive relaxation method can simulate more and more conditions like fluid-solid two way coupling simulation [LCD * 20]. The LBM has become very popular for CFD simulations as it can provide fully parallel algorithms [Yeo06]. Since macroscopic variables are translated into particle distribution functions [NCGB95], the simulation is less sensitive to parameter changes during the simulation. Furthermore, this method uses a regular grid which makes the user's input easier to process. Therefore, we can change the parameters cell by cell during the simulation. This method can also provide high performance to achieve interactive rendering times. We use Palabos [LMK * 20] as our LBM software coupled with Unity 3D.

A user study to compare interactive vs. classical simulations

The goal of this experiment is to compare the performance between the interactive simulation condition and the non-interactive one. We aimed to measure the added value of interactive simulation in terms of user experience and performance for decision making in a simple water dam design scenario.

A water dam design as a use case

In this scenario, a wall is fixed on the left to stop a water flow coming from the right. In front of the wall, there are different kinds of obstacles to slow down the water wave. Each scene has a certain type of obstacle and different water quantity, corresponding to a different optimal solution in terms of wall height. These scenes allow us to test our hypotheses on different levels of difficulty. The task of the participants was to find the minimum height of the wall to stop the water wave, and optimal solutions were computed by offline simulations performed before the experiment.

During the simulation, the system detects if the water overflows the wall or the wall is built too high for the water wave. In this case, the participant has to perform another attempt. Once the system detects that the correct height has been achieved after fluid stabilization, the next scene is loaded (see Fig. 1).

In the interactive simulation condition, boundary height can be edited during a simulation in progress to perform the targeted task. In the non-interactive condition, the simulation must be restarted from the initial conditions when modifications are performed.

The simulation grid is composed of 30 × 30 × 72 cells (86400 cells in total), simulating a box whose size is of 1m × 1m × 3.2m. The water surface is estimated and rendered using the water fraction method for each voxel/cell of the grid.

Figure 1: Interactive simulation of water dam conception: 1 -The simulation starts, 2 -The system detects that water is higher than the wall, 3 -User modifies the wall height, 4 -User finds the optimal height.

Hypothesis

We expected that the interactive simulation would have better user performance regarding task completion time. We also assumed that the task was easier in the interactive simulation condition. Indeed, this condition allows users to make more attempts with different setups to achieve the task, reducing the time to make a decision, and shortening the loop between simulation, visualization and analysis.

• H1: Less time is required to find the correct wall height in the interactive simulation condition.

• H2: Interactive simulation allows users to spend less time to observe and analyze the scene and result. • H3: Overall workload is lower with interactive simulation • H4: Interactive simulation allows users to perform more attempts to explore larger solution space.

Experimental process

Because of the Covid-19 sanitary condition, an experiment in the lab was not permitted. Therefore, before the experiment, an email containing an information notice and a consent form was sent to the subjects. Once they signed the document, they were contacted by phone and invited to connect to our computer remotely via Teamviewer. Then, they had to perform a 5-minute familiarisation task with a training scene to familiarize themselves with interactive features and goals. After this training stage, the first session of the experiment started with each mode of simulation (interactive mode or non-interactive condition). After 3 scenes with the same mode of simulation, they were invited to fill the NASA TLX questionnaire [Har06]. Then the subject entered the second session: another training part with the other mode of simulation started followed by 3 other scenes of the experiment. After the second part, the experiment ended with the same questionnaire for the second session.

Participants: 24 participants (13 males and 11 females) aged from 17 to 36 (µ = 27, σ = 4.83) Figure 2: 24 participants were divided in two groups. 12 of them started with interactive simulation then finished with noninteractive simulation. The other 12 started with non-interactive simulation then finished with interactive simulation. The Latinsquare design was constructed for counterbalancing the condition and scene.

Data collection: For each session of the experiment, we logged all the user activities (changing viewpoint, interaction events, editing boundary condition, etc.), the results received from simulation, and completion time. From this raw data, the following measures were extracted:

• Task completion time (TCT): The duration for the user to complete the task. The measurement started when the scene was loaded until the user found the optimal height of the wall. • Number of failures: The number of failures before achieving the task for each scene. • Observation time: The average observation time between each interactive attempt including camera manipulation.

Results

The results presented in this section were considered statistically significant when p < 0.05. All the analyses were performed using R. In the bar plots, error bars show the standard deviation (SD).

Quantitative results: We found the data was not all normally distributed in a Shapiro test. In the following part of this section, we performed Wilcoxon Signed Rank test. We registered 144 trials: 2 conditions × 3 scenes × 24 participants.

For TCT, the Wilcoxon Signed Rank test showed that the task was significantly quicker (p < 0.0001) to achieve with interactive simulation (Mean = 230.4, SD = 50.0) than with non-interactive simulation (Mean = 141.625, SD = 23.4).

For the number of failures, the Wilcoxon Signed Rank test showed that the participants tend to try more times (p = 0.0004) with interactive simulation (Mean = 7.58, SD = 1.32) than with non-interactive simulation (Mean = 9.79, SD = 2.41).

For the observation time, the Wilcoxon Signed Rank test showed that the participants spent more time (p = 0.002) with noninteractive simulation (Mean = 8.26, SD = 3.89) than interactive simulation (Mean = 5.66, SD = 3.21).

Subjective questionnaire: For the NASA-TLX questionnaire, the ratings are averaged to calculate the overall workload. We ran the Wilcoxon Signed Rank test for the overall workload, but we did not found a significant difference (p = 0.346). While the Wilcoxon Signed Rank test showed that participants made less effort (p = 0.02) with interactive simulation (Mean = 3.71, SD = 1.83) than non-interactive simulation (Mean = 4.38, SD = 1.93).

Discussion

From the results, we can first highlight that the participants had 39% better times completion performances when in the interactive simulation condition (H1). The interactive simulation allows the user to modify the wall height more frequently without restarting the simulation. It is much more efficient for interactive simulation to show the result with the modified wall, as it skipped most parts of the water simulation.

We also found that users made more attempts (H4) and spent less time observing the scene (H2) to achieve the task in the interactive simulation. We interpret this as that the interactive simulation costs less time to make more trials. In non-interactive mode, the participants tried to analyze the result after having the result of the first simulation as they would need to wait longer time to get the results of the modified simulation. Interactive mode let users to make more attempts in limited time to achieve the goal or get a better solution.

On the contrary to our expectations, we did not find significant differences in the overall workload, which rejects H3. From the discussion with the participants, we observed that as the participants were not experts in fluid dynamics, they had difficulties in predicting the movement of the water. We also assumed that the time needed to wait for the end of the simulation in the non-interactive condition might decrease the whole cognitive load. In the meanwhile, we did find that there were significant differences in effort. We interpret this as that since users tended to make more attempts to achieve the task, they did not make efforts to try to understand the simulation. Furthermore, as participants spent less time to achieve the task, it appeared much less difficult to them in performing it.

Conclusion and perspectives

In this work we present a platform that couples fluid simulation tools with Unity 3D to realize the visualization and interaction with a fluid simulation. We address the complex issue of modifying boundary conditions during a simulation in progress by using the LBM. We conducted a first user study to evaluate the added value of interactive simulations compared to the conventional approach in a water dam design scenario. The first results show that users have better performance in achieving the task with an interactive simulation in task completion time. As interactive simulation is less time consuming, more attempts of simulation can replace some work of analyzing the simulation results, which also makes users feel easier to achieve the task. In the meanwhile, the interactive simulation does not reduce the overall users' workload in achieving the task.

Our platform coupling Unity 3D and classical simulation tools are presently mainly used to conduct user studies to evaluate the added value of interactive fluid simulation and to compare different contexts such as classical and immersive environments, augmented reality context, or collaborative ones. We plan in short term to conduct further works to confirm the benefit of interactive simulation to larger fluid simulations, including other interaction contexts such as immersive and augmented reality to support collaboration. We will especially test if using ecological movement around the fluid phenomena using tracking and if collaboration features increase the performance and the quality of the decision-making process, and decrease the cognitive load. This platform also lays the first bricks for a serious game fluid simulation application. Serious games are used for nonentertainment purposes and are applied on different domains [START_REF] Backlund | Games and traffic safety-an experimental study in a gamebased simulation environment[END_REF], such as education, oil and gas application, health care, etc. Serious games enable learners to adapt learning to their cognitive needs and to increase intrinsic motivation of student [START_REF] Malone | Toward a theory of intrinsically motivating instruction[END_REF]. Fold'It is one of the most famous examples. In our case, our goals are to design scenarios and to work on gameplays to allow students to get a better understanding of the fluid physics phenomena. Inspired by Fold'It, we will use our platform as a serious game to illustrate fluid mechanics concepts through appropriate scenarios designed with experts, and to know if serious game and citizen science approaches could solve practical fluid mechanical problems.

  There are many works dedicated to interactive data visualization or real-time data visualization in the CFD domain. ViSTA FlowLib [SGvR * 03] uses haptic rendering techniques to give a better understanding of the unsteady fluid flows data. Another work [VMN * 09] provides and evaluates multimodal feedback such as sonification during interaction with fluid simulation, especially to address visual overload. Stam J. has pioneered the development of the interactive simulation using an unconditionally stable model [Sta99]. Until today, most applications based on interactive simulation approach use particle-based methods such as Translating Eulerian Grids [CTG10] or Smoothed Particle Hydrodynamics (SPH) [MCG03]

Figure 3 :

 3 Figure 3: Mean TCT (top-left), effort (top-right), number of failure (bottom-left) and observation time (bottom-right). Error bars show standard deviation.

This is the authors' version of the work. It is posted here for your personal use. Not for redistribution. The definitive version of record was published in Proceedings, https://doi.org/10.2312/egs.20211022