Long Li

Franc ¸ois-Xavier X Le Dimet

Jianwei Ma

Arthur Vidard
email: arthur.vidard@inria.fr

A level-set based image assimilation method: applications for predicting the movement of oil spills

Keywords: Image assimilation, Structural information, Level set method, Oil spills. I

In this paper, we present a novel method for assimilating geometric information from observed images. Image assimilation technology fully utilizes structural information from the dynamics of the images to retrieve the state of a system, and thus to better predict its evolution. Additionally, the level set method, which describes the evolution of the geometry shapes of a given system, is taken into account to include the dynamics of the images. This method differs from previous methods of image assimilation in that it takes advantage of Lagrangian information in an Eulerian numerical framework. In our numerical experiments, we apply this technique of image assimilation based on the level set method to an oil pollution problem, to calibrate the initial contours of oil pollutants and to identify diffusion coefficients of the model. Topological merging and breaking of oil slicks are well defined and easily performed by this proposed approach. The results show good agreement between simulated values and observed images.

This information is possibly recorded in the form of pixels representing static images or sequences of such images. For example, these images record evolutions of large weather bodies, such as hurricanes, which can travel great distances and cause large amounts of damages. It is therefore vital to couple these images with numerical models to be able to predict the evolution of events.

First, it was necessary to define what is recorded in these images. It is clear that a large part of the dynamical information is contained in the discontinuities or edges of the images, such as fronts, vortexes, and filaments. These discrepancies are linked to physical properties within a threshold, such as the condensation of water vapor, or the salinity level in an estuary. It is important to use images to reveal the dynamics of a system and identify the underlying physics. For example, in the field of meteorology, it is common to consider small cumulus clouds, mainly at tropical latitudes, as Lagrangian tracers, and to use estimates of their drift rates to in turn estimate wind speed. This can actually be misleading, for instance "Lenticularis Clouds" (as shown in Fig. 1) which are common on mountains, characteristically occur in strong winds, but look almost steady state because they are created by a change of phase in the water particles. Therefore if the drift of these clouds were used as Lagrangian markers to evaluate winds, the results would be in full disagreement with the true situation. Before coupling images and numerical models, it is necessary to identify what physical properties are produced within the images and to ensure that they are included in the numerical model.

Second, it was necessary to mathematically define the images in order to make them compatible with the physical model. This task has two constraints:

1) The discretized images must be compatible with the discretization of the model. In the process of data assimilation, one should compare recorded images with artificial images created through the solution of the numerical model. Therefore, we must construct mappings from the model toward the images.

2) To extract information from the evolutions of these images, one needs to include information about the system's dynamics.

In recent years, several methods of assimilation making full use of structural information contained within images have been proposed. These methods can be divided into two basic categories. The first method estimates motion velocity fields by assimilating pseudo-observations of velocity fields produced from image sequences [START_REF] Herlin | Coupling models and data: which possibilities for remotely-sensed images?[END_REF], [START_REF] Horn | Determining optical flow[END_REF]. The drawback to this approach is that it only takes into account very little physical information about the underlying processes controlling features in the image.

The second method directly assimilates data from an image into a dynamic numerical model in a process called direct image assimilation (DIA). In DIA, image's pixels can de assimilated as a tracer concentration measurement [START_REF] Corpetti | Pressure image assimilation for atmospheric motion estimation[END_REF], alternatively, the structural information can be extracted from the image for instance by a multiscale sparse transform [START_REF] Candès | Fast discrete curvelet transforms[END_REF]- [START_REF] Ma | Multiscale geometric analysis of turbulence by curvelets[END_REF], and assimilated through a model-to-structure operator [START_REF] Le Dimet | Toward the assimilation of images[END_REF]- [START_REF] Chabot | Accounting for observation errors in image data assimilation[END_REF].

Going back to the specific topic of this paper, it is well established that geometry is an important aspect of prediction models because it can be used to identify polluted areas. This geometry has both a Lagrangian character because the pollution spot is advected with the flow and an Eulerian character because the proper physics of the pollutant is taken into account and interact with the environment. The geometry of an oil slick can be added to state variables of the model (including velocity, sea surface temperature, and elevation of the surface) as well as the concentration of oil. However, images of the ocean's surface can often be obfuscated (by a cloud for instance), in that case one would need to employ an inpainting method to reconstruct the shape of oil slick. The geometry of oil slicks is complex. For instance a unique source slick can be divided into several slicks, due to the effects of wind and local turbulence on both the atmosphere and the ocean. Lagrangian methods, which define an individual interface as an ensemble of particles following their own trajectories, do not adequately describe the complex topology of many oil slicks. Currently, the level set method is more suitable for this purpose, as it uses Lagrangian information within an Eulerian numerical framework. With this method, topological merging and breaking are well defined and easily performed. The level set method was introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF] for the purposes of computing and analyzing the motion of an interface in two or three dimensions. Since then, the level set method has been successfully used to implicitly track moving interfaces and automatically detect changes in topology [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF]- [START_REF] Sethian | Level set methods for fluid interfaces[END_REF].

In this paper, we attempt to use the level set method to assimilate recorded images. To study the feasibility of this method, we have chosen to model oil pollution in the ocean. This is an ideal site of study because, to a first order, oil remains on the surface of the ocean, and thus a 2D model can be used to model oil spots. For other oceanographic or meteorological applications, what is seen is the integral of the concentration of the contaminant and a 3D information will be necessary.

The paper is organized as follows. Section 2 describes the theory behind the method, including the application of a conventional formulation of variational data assimilation to a pollution transport problem, as well as the level set method and the framework of the level-set based DIA method. Section 3 details our use of the proposed method to address a global oil spill problem, including the details of our numerical implementation and results. Finally, Section 4 summarizes the main results and outlines possible directions of future research.

II. THEORY

Since the aim of this paper is to combine variational data assimilation and level set methods, both are presented separately in the first two parts of this section. The last part proposes a combined approach tailored for oil spill monitoring.

A. Formulation of VDA to address a pollution transport problem

As mentioned in the introduction, VDA is a class of techniques that combines the mathematical information provided by equations with the physical information yielded by observational data in an optimal way to retrieve the state of a system. To use a VDA method for solving a pollution transport problem, the following components are required [START_REF] Le Dimet | General sensitivity analysis in data assimilation[END_REF], [START_REF] Le Dimet | Sensitivity analysis applied to a variational data assimilation of a simulated pollution transport problem[END_REF]:

A model describing the evolution of the state variable X, which governs flow, and of the concentration C of a given contaminant which is defined as:

                 dX dt = F(X, V ), t ∈ [t 0 , t f ] X(t 0 ) = U, dC dt = G(X, C, Z), t ∈ [t 0 , t f ] C(t 0 ) = W. (1)
A set of control variables: U , V , W , Z. U , W represent the initial conditions, and V , Z are parameters of the models.

A set of observations X obs , C obs ∈ O obs . For the sake of simplicity, we will assume that our observations are continuous in time. They can be obtained by physical measurements of the system state. However, these observations typically comprise partial or indirect measurements of the state variables. For example, The observation space O is not necessarily the same as the state space X , and the observation operator H(•) is defined as the mapping operator from space X onto space O. Similarly, the operator of observation D(•) maps the space of the concentration variable C onto the other observation space. Note that in the specific case of image observation, H(•) and D(•) may include structure extraction operators and become what is called model-to-structure operator in the introduction.

Background estimations U b , V b , W b , Z b of the initial states of model parameters U , V , W , Z. For meteorological applications, this background estimation can be deduced from previous forecasts.

Statistical information. Some estimates of the statistics of errors take the form of a covariance matrix, as with the error covariance matrices Q X , Q C of the observation error and the covariance matrices B U , B V , B W , B Z of the background estimation.

A cost function, J. The cost function contains two terms: the first one, J o , measures the discrepancy between observations and the solution of the model associated with U , V , W , Z. The second term, J b , is the background term, which requires the solution to be located in the vicinity of U b , V b , W b , Z b . This term can be referred to as the Tikhonov regularization, which is applied to remedy ill-posed inverse problems [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF]. The cost function to be minimized is defined in the continuous case as:

J(U, V, W, Z) = J o + J b , (2)
J o = 1 2 t f t0 H(X, t) -X obs (t) 2 Q -1 X dt + 1 2 t f t0 D(C, t) -C obs (t) 2 Q -1 C dt, (3)
J b = 1 2 U -U b 2 B -1 U + 1 2 V -V b 2 B -1 V + 1 2 W -W b 2 B -1 W + 1 2 Z -Z b 2 B -1 Z , (4)
where the unknown parameters V , Z are space-dependent control variables. The norm takes the form of X 2 M = M X, X , where •, • is an inner product operator in Hilbert space.

Variational data assimilation defines the optimal initial conditions U a , W a and model parameters V a , Z a which best fit observations acquired over the time range [t 0 , t f ]. To minimize the equation (2), A necessary condition for this optimality is expressed by the Lagrange-Euler equation, which requests the gradient of J with respect to U , V , W , Z to be equal to 0.

However, evaluating this condition (i.e. computing ∇J) is made difficult by the facts that J can be an implicit function of these variables and that the size of the problem can be large. Therefore, during the process of variational data assimilation, ∇J is often computed using the so-called adjoint model:

                 dP dt + ∂F ∂X * P + ∂G ∂X * R = ∂H ∂X * Q -1 X (H(X) -X obs ), t ∈ [t0, t f ] P (t f ) = 0, dR dt + ∂G ∂C * R = ∂D ∂C * Q -1 C (D(C) -C obs ), t ∈ [t0, t f ] R(t f ) = 0, (5)
where * denotes the transpose, and P , R represent the so-called adjoint variables. The gradient of the cost function with respect to these control variables is then given by the following equation:

        ∇J(U ) ∇J(V ) ∇J(W ) ∇J(Z)         =          -P (t 0 ) + B -1 U (U -U b ) - t f t0 [ ∂F ∂V ] * P dt + B -1 V (V -V b ) -R(t 0 ) + B -1 W (W -W b ) - t f t0 [ ∂G ∂Z ] * Rdt + B -1 Z (Z -Z b )          . (6)
Finally, U a , V a , W a , Z a are determined by running a descent-type optimization algorithm of the Newton-type or L-BFGS variety [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF].

B. Basic level set method

Mathematically, a given shape (e.g., an oil slick) can be represented by a subdomain Ω(t) ⊂ R d . Let to define the boundary ∂Ω(t) as the zero level set of a mapping φ : R d → R:

∀x ∈ R d ,          φ(x) < 0, x ∈ Ω φ(x) = 0, x ∈ ∂Ω φ(x) > 0, x / ∈ Ω. (7)
The function φ(t, x) could be the signed distance from (t, x) to ∂Ω(t, x). It defines the evolution of the shape whose movement is governed by some velocity field V (t, x) : R d → R d . The motion of Ω(t, x) follows the advection equation:

∂φ ∂t (t, x) + V (t, x) • ∇φ(t, x) = 0, on [t 0 , t f ] × R d . (8)
In fact, the motion of Ω(t, x) depends only on v(t, x), the component of V (t, x) that is normal to the boundary of Ω(t, x).

Therefore, the motion of the boundary verifies the Hamilton-Jacobi equation:

∂φ ∂t (t, x) + v(t, x) ∇φ(t, x) = 0, on [t 0 , t f ] × R d . (9)
for which we need to provide an initial condition:

φ(t 0 , x) = φ 0 (x). (10)
A more general form of this equation will be introduced in Section II-C2. However, during the evolution of the level set equation, the function φ(t) can become very flat or very steep, both of which are undesirable effects that can lead to high numerical errors either in the reconstruction of the zero level set or elsewhere in the numerical schemes. An optimal control based approach of this important reinitialization aspect is discussed in Appendix B.

C. Level-set based direct image assimilation 1) Edge capturing and image observation operators: For considering the problem of pollution within a fluid, the variable Ω can be used to represent one or more patches of the pollutant. For adapting our theoretical methodology to analyze pollution, we express the size of the visible pollutant spot in terms of the pollutant's concentration, in which case Ω and the threshold function θ(•) are defined as:

Ω(t, x) = {x(t) ∈ Ω|C(t, x) ≥ }, (11)
θ(φ obs (t)) =    φ obs (t), C(t, x) ≥ 0, C(t, x) < ( 12
)
where C is the concentration of the pollutant, Ω ⊂ R 2 is the domain taken into account, φ obs is the observation image and represents the concentration threshold above which the spots of pollutant are visible.

In this paper, we use the wavelet modulus maximum method [START_REF] Canny | A computational approach to edge detection[END_REF]- [START_REF] Mallat | Zero-crossings of a wavelet transform[END_REF] to exploit image edge detection after undergoing thresholding. The result obtained after edge detection is expressed as W(θ(φ obs )), from which we can distinguish the inside and outside edges clearly. Some edge detection methods based on wavelet can refer to [START_REF] Alonso | Edge enhancement algorithm based on the wavelet transform for automatic edge detection in sar images[END_REF]- [START_REF] Niedermeier | Detection of coastlines in sar images using wavelet methods[END_REF]. Following this procedure, we apply an edge-based signum operator H W→S (•) to the edge-detected images W(θ(φ obs )). The procedure is as follows:

H O→S (φ obs ) = H W→S (W(θ(φ obs ))) =          -1, x ∈ Ω 0, x ∈ ∂Ω 1, x ∈ Ω\ Ω (13)
where H O→S (•) is the image-to-structure operator, and the structure image is expressed as φ S (t) = H O→S (φ obs (t)). This construction procedure is presented in Fig. 4a-4d. Additionally, we constructed another adapted observation operator, called model-to-structure operator. This operator maps the signed distance function toward the signed function, which links variables of the level set equation with the structural information of images. This model-to-structure operator is defined as:

H φ→S (φ) =          -1, φ < 0 0, φ = 0 1. φ > 0 (14)
This procedure is shown in Fig. 4e-4f.

2) General level set model: In order to apply the model-to-structure operator one has to know the temporal evolution of the edge mapping φ. This can be obtain by application of level set techniques. In general, the motion of the interface ∂Ω(t)

is analyzed by convecting the level set function φ(t, x, y) with the external velocity field V = (u, v) and diffusing along both the normal direction by a constant value, a, as well as proportionally to the mean curvature, κ. This elementary equation can be expressed as a fully nonlinear second-order partial differential equations:

   ∂φ ∂t + V • ∇φ + a ∇φ -σκ ∇φ = 0, on [t 0 , t f ] × Ω φ(t 0 ) = φ 0 , (15)
κ = div( ∇φ ∇φ ) = ( ∂φ ∂x ) 2 ∂ 2 φ ∂y 2 -2 ∂φ ∂x ∂φ ∂y ∂ 2 φ ∂x∂y + ( ∂φ ∂y ) 2 ∂ 2 φ ∂x 2 ∇φ 3 . ( 16
)
If the term contained curvature σ is present, closed curves of arbitrary shape will collapse to a circle and then disappear. Unlike the curvature, as it constantly evolves under a, a relatively closed smooth curve will contain sharp edges and the topological structures of the curve may change. In practical cases, the values of the diffusion coefficients σ, a depend on the physical characteristics of the pollutant. In this model, we adopt a periodic boundary condition.

3) Continuous variational model and the adjoint model: We use level set equation [START_REF] Osher | Level set methods: An overview and some recent results[END_REF] as the governing equation to predict the evolution of oil slicks over time. The method used to determine the structure of the images, φ S (t), is described in [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF].

The control variables are the initial condition and diffusion coefficients σ, a. In this case, the total cost function associated with the observed images J o , the backgrounds for initial conditions J b φ(t0) , and their coefficients J b σ,a , can be expressed as:

J = J o + J b φ(t0) + J b σ,a , (17)
J o + J b φ(t0) = 1 2 t f t0 H φ→S [K t (φ(t 0 ), σ, a)] -H O→S [φ obs (t)] 2 L 2 ( Ω) dt + α 2 φ(t 0 ) -φ(t 0 ) b 2 L 2 ( Ω) , (18)
J b σ,a = β 2 (R p1,q1 (σ) + R p2,q2 (a)), (19)
where K t (•) presents the forward operator of the differential equation ( 15) from t 0 to t, and the assimilation window is set as

[t 0 , t f ].
The variables α, β denote regularization parameters, and the expression φ(t 0 ) b represents the background fields. When constructing an adjoint model, an additional difficulty is that of the non-differentiability inherent in mapping H φ→S (•). Here, we modify the edge-based signum function to its differentiable form:

H φ→S (φ) = φ φ 2 + ε , ( 20
)
where ε is a small constant, and Section III-A details the method for determining it. The R p,q (x) is a regularization term for σ and a. If the coefficients are defined as being constant values, the regularization term takes the following form:

R p,q (x) =          1 2 (x -q) 2 , x > q 0, p ≤ x ≤ q 1 2 (x -p) 2 , x < p (21)
where p, q represent the lower and upper bounds of coefficients, respectively.

The new cost functions ( 17)-( 19) are the main contribution of this paper. In the following, we will focus on how to compute and numerically implement above optimization problem of the cost functions, and apply the new method to oil spill VDA.

According to the definition of the adjoint operator in Hilbert space and the method of integration by parts, we can obtain a continuous adjoint model from the continuous forward model [START_REF] Osher | Level set methods: An overview and some recent results[END_REF]:

               ∂R ∂t = ∂(A-u-F )R ∂x + ∂(B-v-G)R ∂y -∂ 2 CR ∂x 2 -∂ 2 DR ∂y 2 + ∂ 2 ER ∂x∂y +[ ∂H φ→S ∂φ ] * (H φ→S (φ) -H O→S (φ obs )), on [t 0 , t f ] × Ω R(t f ) = 0, R| ∂ Ω = 0, (22)
where

A = σ 2 ∂φ ∂y [( ∂φ ∂x ) 2 ∂ 2 φ ∂x∂y -( ∂φ ∂y ) 2 ∂ 2 φ ∂x∂y + ∂φ ∂x ∂φ ∂y ∂ 2 φ ∂y 2 -∂φ ∂x ∂φ ∂y ∂ 2 φ ∂x 2 ] ∇φ 4 , (23)
B = σ 2 ∂φ ∂x [-( ∂φ ∂x ) 2 ∂ 2 φ ∂x∂y + ( ∂φ ∂y ) 2 ∂ 2 φ ∂x∂y -∂φ ∂x ∂φ ∂y ∂ 2 φ ∂y 2 + ∂φ ∂x ∂φ ∂y ∂ 2 φ ∂x 2 ] ∇φ 4 , (24)
C = σ ( ∂φ ∂y ) 2 ∇φ 2 , D = σ ( ∂φ ∂x ) 2 ∇φ 2 , E = σ 2 ∂φ ∂x ∂φ ∂y ∇φ 2 , F = a ∂φ ∂x ∇φ , G = a ∂φ ∂y ∇φ , (25)
[ ∂H φ→S ∂φ ] * (•) = ( 1
φ 2 + ε - φ 2 (φ 2 + ε) 3 2
)(•).

Then, the gradient can be expressed as:

     ∇J(φ 0 ) ∇J(σ) ∇J(a)      =      -R(t 0 ) + (φ(t 0 ) -φ(t 0 ) b ) - t f t0 Ω κRdΩdt + ∇R p1,q1 (σ) t f t0 Ω (( ∂φ ∂x ) 2 + ( ∂φ ∂y ) 2 ) 1 2 RdΩdt + ∇R p2,q2 (a)      . (27)

III. APPLICATIONS TO EVOLUTION OF OIL SPILLS

A. Numerical implementation

The evolution of the fluid can be modeled using the state variable X = (u, v, h), whose components verify the shallow water equations in the conservative form expressed as:

             ∂h ∂t + ∂uh ∂x + ∂vh ∂y = 0, on [t 0 , t f ] × Ω ∂uh ∂t + ∂ ∂x (u 2 h + gh 2 2 ) + ∂uvh ∂y = h(f v -g ∂H ∂x ), on [t 0 , t f ] × Ω ∂vh ∂t + ∂uvh ∂x + ∂ ∂y (v 2 h + gh 2 2 ) = h(-f u -g ∂H ∂y ), on [t 0 , t f ] × Ω (28)
where a numerical simulation is performed on a rectangular domain

Ω = [0, 2] × [0, 2].
Additionally, H is the height of the orography, h = z -H is the water elevation, and z is the free surface elevation. The variables u(t, x) and v(t, x) represent the zonal and meridional components of the current velocity, respectively, and g represents the acceleration due to gravity. In this model, the Coriolis parameter on the β-plane is defined as f = f 0 + βy. The following numerical values are used for these experiments: f 0 = 0.25s -1 , β = 0.0406m -1 s -1 , g = 9.81ms -2 and H = 0. We integrate this scheme forwards-in-time using the Lax-Wendroff scheme [START_REF] Liska | 2d shallow water equations by composite schemes[END_REF], which is accurate to a second-order.

To discretize the level set equation, we use an upwind scheme for convection term discretization, and employ a central difference method for diffusion terms. For the sake of simplicity, time integration is performed using the Euler scheme.

In applied cases, a third-order accurate TVD Runge-Kutta scheme could also be used. Detailed descriptions of high-order discretization schemes for level set equations can be found in the literature [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF], [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF], [START_REF] Jiang | Weighted eno schemes for hamilton-jacobi equations[END_REF], [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. In our model, the domain is discretized on a N ×M = 81×81 uniform square grid. There are a total of 201 time steps. Each time step ∆t is set to 0.0005s, and the assimilation window is [0, 0.1s]. The variable ε, as expressed in [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF], is usually defined as (∆x) 2 . Additionally, in the assimilation experiments, we assume that the current velocity fields u, v are known.

Because the forward model and the observations are both discrete, we can use the discrete adjoint model in this test to obtain the gradient of the control variables. The validation of both the tangent linear and the adjoint is of the utmost importance for nonlinear models. The methods of producing the discrete adjoint model and testing the gradient are shown in detail in the Appendix A.

In this section, we have contributed directly assimilate the images of the structures of oil pollutants φ S (t) into numerical models to correct the contours of initial oil pollutants and diffusion coefficients of models. True structure images φ S (t) are taken at each time step and at every spatial grid point using the image-to-structure operator H O→S (•). In this test, the sequences of structure images proceeds according to H φ→S (φ), where the level set function φ is produced from the forward model. The sequences of synthetic structure images can thus be expressed as:

φ S = H φ→S [K t (φ(0) t , σ t , a t )], (29)
where φ(0) t represents the true initial state of the system. To test the efficiency of this proposed method, we have contributed a missing data point at t = 0 to generate the initial image of the structure of the oil slicks, φ S (0), which is shown in Fig. 5a.

To assess the background fields of initial conditions, we assumed the missing oil slick began as a rectangular domain whose zero level set is shown as the solid red line in Fig. 6a. The true model parameters are σ t = 0.1 and a t = 0.5, and weights are set at α = 1 and β = 2.

To assess the performance of the numerical test, we computed the root mean square error (RM SE) of the solution estimated relative to the true value. The RM SE is defined as

RM SE φ(0) a = 1 | Ω| x∈ Ω (φ(0, x) a -φ(0, x) t ) 2 , ( 30
)
RM SE σ a = |σ a -σ t |, ( 31
)
RM SE a a = |a a -a t |, (32)
where φ(0) a , σ a , a a represent optimal solutions.

B. Numerical results

To illustrate the level set method, we first show the numerical results, the level set function φ(t), produced from the forward model ( 15) at three different moments in Fig. 2. The level set function at the first moment derives from a synthetic oil spot image, which is obtained by constructing a signed distance function based on the oil spot edges. In this paper, we only take into account the affect of the shear flow fields provided by the shallow water equation [START_REF] Niedermeier | Detection of coastlines in sar images using wavelet methods[END_REF]. Fig. 3 presents the zero level set of φ(t) at corresponding moments of Fig. 2, representing the point set (x(t), y(t)) at which φ(x(t), y(t), t) = 0. The zero level sets express the edges of oil slicks. The changes of the contours reveal the evolution of the shapes of the polluted areas, which drifted on the surface of the water at a velocity determined by the current and the wind in practical cases.

In the framework of level-set based DIA, we need two observation operators for the purpose of mapping the numerical solutions and the observations into the same structure space. One is the so-called image-to-structure operator H O→S (•),

the other is the model-to-structure operator H φ→S (•). This construction procedure is presented in Fig. 4. First, a synthetic observation image φ obs with noise is given displayed in Fig. 4a. The light gray part in the image indicates the oil pollution.

Edge detection result W(φ obs ) by the method of wavelet modulus maximum is shown in Fig. 4b. Indeed, this method can balance noises and edges. Fig. 5 displays portions of the structure image sequences φ S (t) constructed by the method defined in [START_REF] Liska | 2d shallow water equations by composite schemes[END_REF], which are used as true observations in numerical implementation. In order to test the efficiency of the proposed method, the structure image sequences include a missing image at t = 0. In practical cases, it can be considered as an image obfuscated by clouds.

To clearly compare numerical results, we can examine the zero level set of the analysis value φ(0) a (shown as a solid red line in Fig. 6b) and the first guess or background fields φ(0) b (shown as a solid red line in Fig. 6a), with respect to the corresponding true initial states φ(0) t (shown as dashed black line in Fig. 6a and6b). Fig. 6b suggests that the optimal solution will eventually converge to true values with increasing numbers of iterations. A part of the first observation is missing, the proposed method is capable of recovering it. Additionally, the evolution of the l 2 norm of the gradient and its cost function are shown in Fig. 6c, in which the l 2 norm of the gradient gradually moves closer to zero. Fig. 7 displays the evolution of the estimated coefficients σ, a, as well as their gradients and cost function. This shows that the optimal coefficients will be close to the true values with increasing numbers of iterations.

Finally, Table I and II list the results of our numerical tests, including the RM SE, gradients and cost function values, which change with increasing numbers of iterations. These results suggest that our proposed method performs well.

IV. CONCLUSIONS AND FUTURE WORK

A. Conclusions

In this paper, we propose a level-set based image assimilation method for assimilating geometric information from the observed images into numerical models. We present numerical experiments which can be applied to better predict oceanic oil spill dispersion. In general, our method takes into account both the complex topological structures of oil slicks and the geometric shape of the polluted areas to avoid delaying cleanup efforts by making inaccurate predictions from observations. This method represents a key step to better predict the evolutions of oil spill. Furthermore, the adapted edge-based signum observation operator could eventually link the data from the level set function to the structural information contained in satellite images. Our assimilation test yielded promising results.

B. Future work

In this paper it is assumed that both the external velocity and the concentration field are known. In practice they have to be estimated as well. This can be done sequentially with the proposed approach, but for the sake of consistency one may prefer to use a more integrated approach called splitting method for global oil spill problem. We first need to define the framework of optimal control.

1) Model and control variables:

• A model of flow: F(X, U ) = 0, where X is the state variable describing the evolution of the flow and initial conditions X(t 0 ) = U is a control variable.

• An equation of concentration: G(X, C, W )= 0, where C is the concentration of pollutant and C(t 0 ) = W is a control variable.

• An equation for the evolution of an image: K(X, C, I, φ 0 , λ) = 0, where I is the image and φ 0 is a control variable defining the initial level set associated with the spot of pollution. λ is an additional parameter (e.g., the viscosity of oil on the ocean's surface). In this paper, λ represents (σ, a).

2) Optimal control: The goal of this problem is to determine (U a , W a , φ a 0 , λ a ) while minimizing the cost function; this is expressed as:

J(U, W, φ 0 , λ) = J F + J G + J I = 1 2 t f t0 Γ X X -X obs 2 + Γ C C -C obs 2 + Γ I I -I obs 2 dt, (33)
where Γ X , Γ C and Γ I are observation operators. Throughout the process of quantifying minimum values, U n ,W n , (φ 0 ) n can be obtained from the iterative formula:

• U n+1 minimizes J F with the constraint F(X, U ) = 0,

• W n+1 minimizes J G with the constraint G(X n+1 , C, W ) = 0,
• ((φ 0 ) n+1 , λ n+1 ) minimizes J I with the constraint K(X n+1 , C n+1 , I, φ 0 , λ) = 0.

The two first formulae are well known in the VDA community [START_REF] Le Dimet | General sensitivity analysis in data assimilation[END_REF], [START_REF] Le Dimet | Sensitivity analysis applied to a variational data assimilation of a simulated pollution transport problem[END_REF]. This splitting method solves large numerical problems using a sequence and iteration of smaller, simpler problems. In the numerical experiments of this paper, we focus on the third part of this algorithm.

Additionally, in the future we would like not only to further consider the physical properties controlling the movements of oil spills but also to utilize real satellite images of oil slicks.

APPENDIX A THE DISCRETE ADJOINT MODEL AND THE GRADIENT OF THE CONTROL VARIABLES

We consider the cost function J as follows:

J(u 0 , ν) = 1 2 ∆t N n=1 (H(u n ) -u n obs ), (H(u n ) -u n obs ) + 1 2 (u 0 -u b ), (u 0 -u b ) , (34)
where the discrete model is defined as follows:

u n+1 = M(u n , ν), (35)
where u 0 represents the initial condition, having been arranged in a column vector, and ν is a parameter of the model. The The Gatéaux derivative of J in the direction of perturbation δu 0 , δν is defined as:

Ĵ(u 0 , ν, δu 0 , δν) = lim α→0 J(u 0 + αδu 0 , ν + αδν) -J(u 0 , ν) α

= δu 0 , (u 0 -u b ) + ∆t N n=1 û, ( ∂H ∂u (u n )) * (H(u n ) -u n obs ) , (36)
where * denotes the transpose. We can obtain the discrete tangent linear equation (DTLE) from the discrete model given in (35):

ûn+1 = ∂M ∂u (u n , ν)û n + ∂M ∂ν (u n , ν)ν. (37)
Then, we can take the scalar products of (37) with the adjoint variable p n+1 and sum them all in time:

N n=1 ûn+1 , p n+1 = N n=1 ( ∂M ∂u (u n , ν)û n ), p n+1 + N n=1 ( ∂M ∂ν (u n , ν)ν), p n+1 , (38)
which is the same as:

N n=1 ûn , p n -û1 , p 1 + ûN+1 , p N +1 = N n=1 ûn , ( ∂M ∂u (u n , ν)) * p n+1 + N n=1 ν, ( ∂M ∂ν (u n , ν)) * p n+1 . (39)
Next, we can assumes that p N +1 is 0 and use the fact that û1 = δu 0 . This yields the following equation:

N n=1 ûn , (p n -( ∂M ∂u (u n , ν)) * p n+1 ) = δu 0 , p 1 + N n=1 ν, ( ∂M ∂ν (u n , ν)) * p n+1 . ( 40
)
The Gatéaux derivative of J can be expressed as:

Ĵ(u 0 , ν, δu 0 , δν) = δu 0 , ∇ u0 J + ν, ∇ ν J = δu 0 , (u 0 -u b ) + ∆t N n=1 û, ( ∂H ∂u (u n )) * (H(u n ) -u n obs ) . (41)
If the so-called variable p is the solution of the following equation,

p n -( ∂M ∂u (u n , ν)) * p n+1 = ∆t( ∂H ∂u (u n )) * (H(u n ) -u n obs ), (42)
then

δu 0 , ∇ u0 J + ν, ∇ ν J = δu 0 , (u 0 -u b ) + δu 0 , p 1 + N n=1 ν, ( ∂M ∂ν (u n , ν)) * p n+1 = δu 0 , ((u 0 -u b ) + p 1 ) + N n=1 ν, ( ∂M ∂ν (u n , ν)) * p n+1 . (43)
Finally, we can obtain the gradient:

∇ u0 J = (u 0 -u b ) + p 1 , (44)
∇ ν J = N n=1 ( ∂M ∂ν (u n , ν)) * p n+1 . (45)
Next, we can use the test of the gradient to validate both the tangent linear and the adjoint models. The test of the gradient is based on the Taylor expansion of the cost function, which is written as follows:

J(x + αδx) = J(x) + α ∇ x J(x), δx + O(α δx ), (46)
where J is the cost function and x is a symbolic representation of the initial condition or the parameters of the models.

Rearranging (46) leads to the following expression:

lim α→0 J(x + αδx) -J(x) α ∇ x J(x), δx = 1. (47)

APPENDIX B THE REINITIALIZATION BASED OPTIMAL CONTROL

In this section, we take into account the process of reinitialization mentioned in Section II-B in the framework of level-set based DIA. In order to maintain the numerical stability of ( 9), one needs to reinitialize the level set function so that, after a period of time, it moves closer to the signed distance function of the actual shape. This reinitialization procedure is expressed as [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]:

   ∂ψ ∂τ = sgn(φ)(1 -∇ψ ), on [τ 0 , τ f ] × R d ψ(τ 0 ) = φ(t), (48)
where sgn(•) is the signum function and the variable ψ is called renormalized level set function. This makes possible the computation of ψ by setting ∇ψ(τ f ) = 1 as the long-time limit τ f → ∞, related to the corresponding evolution equation.

Then φ(t) will be substituted by ψ(τ f ).

For numerical implementation, we denote the discrete equation ( 15) and (48) respectively as:

φ n+1 = A(φ n ), (49)
ψ m+1 = B(ψ m ), ψ 1 = φ n+1 , (50)
where A(•) and B(•) are discrete forward operators from the moment t = n∆t to t = (n+1)∆t and t = m∆t to t = (m+1)∆t.

R d in (48) is set to Ω. And the process of initialization is used for each time step. We can obtain the DTLE from the discrete model given in (49) and (50):

φn+1 = ∂A ∂φ (φ n ) φn , (51) ψm+1
= ∂B ∂ψ (ψ m ) ψm , ψ1 = φn+1 . (52)
Then we have:

ψm+1 = ∂B ∂ψ (ψ m ) • • • • • ∂B ∂ψ (ψ 1 ) ψ1 = ∂B ∂ψ (ψ m ) • • • • • ∂B ∂ψ (ψ 1 ) φn+1 = ∂B ∂ψ (ψ m ) • • • • • ∂B ∂ψ (ψ 1 ) ∂A ∂φ (φ n ) φn . (53)
The variable ψm+1 represents the updated φn+1 . In addition, the adjoint of the operator

B = ∂B ∂ψ (ψ m ) • • • • • ∂B ∂ψ (ψ 1 ) ∂A ∂φ (φ n )
used in the adjoint model (42) will be expressed as follows:

(B) * = [ ∂B ∂ψ (ψ m ) • • • • • ∂B ∂ψ (ψ 1 ) ∂A ∂φ (φ n )] * = [ ∂A ∂φ (φ n )] * [ ∂B ∂ψ (ψ 1 )] * • • • • • [ ∂B ∂ψ (ψ m )] * . (54)
Because of the linearity of the Jacobian matrix, the adjoint of the operator is equivalent to its transpose. In practical cases, the number m is fixed in each iteration, which represents the duration of the reinitialization. For the sake of simplicity, we will not take this procedure into account in this paper.

Fig. 1 :

 1 Fig. 1: Lenticular clouds image of Eulerian type (Photo by Dahlia Rudolph at Mt. Rainer, October 5, 2011).

  Fig. 4c-4d depict the structure and modified structure image (i.e. the edge-based signum function and smoothing signum function image ϕ S = H O→S (φ obs )). This can show the structure information of edge clearly. From the aspects of variables from the forward model, the level set function φ and the structure image ϕ S = H φ→S (φ) are shown in Fig. 4e-4f. As a result, Fig. 4d approximately equals to Fig. 4f.

variables u n obs , n = 1 ,

 1 • • • , N represent the observations sampled at t = n∆t, and u b is the background field or the first guess. H(•) represents a nonlinear observation operator which combines the model variables with observations. M(•) is the discrete operator of the forward model from the moment t = n∆t to t = (n + 1)∆t.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Level set function φ(t) at three moments produced by forward model. (a) Initial moment t = 0. (b) After 100 time steps. (c) Final moment t = 0.1s.

Fig. 4 :

 4 Fig. 4: (a) Synthetic observation image φ obs . (b) Edge detection result W(φ obs ). (c) Structure image. (d) Modified structure image ϕ S = H O→S (φ obs ). (e) Signed distance function image φ. (f) Modified structure image ϕ S = H φ→S (φ).

Fig. 5 :Fig. 6 :

 56 Fig. 5: Portions of the sequences of structure images φ S (t) including a missing image at t = 0. (a) Initial moment t = 0. (b) t = 0.05s. (c) Final moment t = 0.1s.

Fig. 7 :

 7 Fig. 7: Assimilation results of coefficients. (a) Evolution of coefficients σ, a estimated with the iterations (σ t = 0.1, a t = 0.5). (b) Evolution of the l 2 norm of the gradients and cost function.

TABLE I :

 I ASSIMILATION RESULTS OF INITIAL CONDITIONS

	Iterations RM SE φ(0)	∇J(φ(0))	J(φ(0))
	70	3.5748e-02	1.3687e-01	1.0920

TABLE II

 II

		: ASSIMILATION RESULTS OF COEFFICIENTS
	Iterations	RM SE σ	RM SE a	∇J(σ)	∇J(a)	J(a, σ)
	30	1.7118e-02 9.5625e-02 6.9408e-04 3.9667e-03 3.5445e-01
	60	4.9719e-03 2.7232e-02 1.9990e-04 1.1045e-03 2.2405e-01

ACKNOWLEDGMENT

This work is supported by NSFC (grant number: 91330108, 41374121, 61327013), the Fundamental Research Funds for the Central Universities (grant number: HIT.PIRS.A201501).