Long Li

Arthur Vidard

François-Xavier Le Dimet

Jianwei Ma

Topological data assimilation using Wasserstein distance

published or not. The documents may come

Introduction

Numerical prediction of geophysical fluids requires a good knowledge of initial and boundary conditions, external forces, and internal physical parameters. Such information is indirectly accessible from observation and previous or external numerical simulations; however, they are partial and uncertain. Data assimilation (DA) is an efficient method for calibrating the state of dynamical equations by combining all the available heterogeneous information (i.e., mathematical models based on physical laws, observation, and a priori knowledge). Variational data assimilation (VDA) [1, 2] and filtering DA [3, 4, 5] are the two most common methods. Their applications include numerical weather prediction (NWP), marine and pollution monitoring [START_REF] Nodet | Variational assimilation of Lagrangian data in oceanography Inverse problems[END_REF], hydrologic forecasting, soil moisture monitoring, land surface flux and wave speed estimation [START_REF] Fehrenbach | Variational data assimilation to estimate the velocity in the wave equation Inverse[END_REF].

The problem of determining the unknown variables called control variables is considered as a nonlinear least-squares problem in the framework of VDA and is solved using the adjoint method [2]. This is the preferred approach in numerical weather prediction, in its incremental formulation [START_REF] Courtier | A strategy for operational implementation of 4D-Var, using an incremental approach[END_REF].

For the pollution problem, reliable prediction from the concentration transport model is not possible when the initial concentration information is limited, observations are scarce, and velocity fields from the ocean and wind are inaccurate (e.g., oil spill).

In this scenario, the level-set method can be considered. It was pioneered by Osher and Sethian, aiming at computing and analyzing the motion of an interface in two or three dimensions [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Osher | Level set methods: an review and some recent results[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF][START_REF] Gibou | A review of level-set methods and some recent applications[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF]. Topological merging and breaking are well defined and easily performed by the method of implicit representation of a time-evolving contour. At present, the level-set method has been successfully applied to tumor growth [START_REF] Martins | Multiscale models for the growth of avascular tumors[END_REF][START_REF] Macklin | An improved geometry-aware curvature discretization for level set methods: application to tumor growth[END_REF], flame propagation [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], and wildland fire propagation [START_REF] Mallet | Modeling wildland fire propagation with level set methods Computers and Mathematics with[END_REF]. However, owing to the uncertainty in the model such as the parameters of the indirect measurement and inaccurate initial contour position, several data-driven methods have gained increased attention. The interface evolution can be corrected and become reliable by combining the model variables with observation and statistical methods. Among these, the level-set methods based on VDA [2] have been used for motion estimation and structure tracking [START_REF] Li | A level-set-based image assimilation method: potential applications for predicting the movement of oil spills[END_REF][START_REF] Papadakis | An variational technique for time consistent tracking of curves and motion[END_REF][START_REF] Lepoittevin | Continuous tracking of structures from an image sequence[END_REF][START_REF] Lepoittevin | Object's tracking by advection of a distance map[END_REF].

In [START_REF] Li | A level-set-based image assimilation method: potential applications for predicting the movement of oil spills[END_REF], the theoretical framework for assimilating the contour information in ocean image observation was illustrated and a shape-fitting functional used in image processing was considered. The continuous cost function and adjoint-based gradient were given, and the results showed its potential for oil spill and pollutant transport. In addition, there have been attempts to use the level-set structure to indirectly retrieve the flow fields in fluids [START_REF] Lepoittevin | Continuous tracking of structures from an image sequence[END_REF][START_REF] Lepoittevin | Object's tracking by advection of a distance map[END_REF]. Moreover, filtering-based methods have been proposed for front tracking the wildfire spread by assimilating observation [START_REF] Rochoux | Towards predictive data-driven simulations of wildfire spread-Part 1: reduced-cost ensemble Kalman filter based on a polynomial chaos surrogate model for parameter estimation Natural Hazards and Earth System[END_REF][START_REF] Rochoux | Towards predictive data-driven simulations of wildfire spread-Part 2: ensemble Kalman filter for the state estimation of a fronttracking simulator of wildfire spread Natural Hazards and[END_REF].

However in DA, because of the displacement and obstacle in the remote-sensing observer, position and shape errors usually arise in observation data (e.g., image-type sequence) compared with true state and background simulation. To avoid this effect, in traditional approaches, a step of data alignment is usually included before using DA for correcting the position error [START_REF] Hoffman | A technique for assimilating SSM/I observations of marine atmospheric storms: tests with ECMWF analyses[END_REF][START_REF] Ravela | Data assimilation by field alignment[END_REF][START_REF] Nehrkorn | Correcting for position errors in variational data assimilation[END_REF]. A shape similarity measure based on the Chan-Vese contour fitting functional was studied for DA [START_REF] Rochoux | Front shape similarity measure for shape-oriented sensitivity analysis and data assimilation for Eikonal equation[END_REF]. This novel method could more efficiently compare the position information of the level-set structure from observation and model variables for the front-tracking problem. However, the background information could not be made full use in these DA methods. One possible solution would be to use optimal transport-based metrics.

The optimal transport theory has been significantly developed and applied to wide fields since it was introduced by Monge [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais De l'Imprimerie Royale[END_REF]. In [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovitch mass transfer problem[END_REF], the transportation problem between two probability measures was described as the minimization of the kinetic energy. The optimal "geodesics" maintained the mass conservation over time, driven by the flow in a continuity equation. A series of numerical methods have been proposed for obtaining the optimal transportation such as the primal-dual method [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF], Fisher information regularization method [START_REF] Li | Computations of optimal transport distance with Fisher information regularization[END_REF] and parallel method [START_REF] Li | A parallel method for Earth Mover's distance[END_REF]. Recently, multi-physics optimal transportation, in which a constraint term is added to the cost function based on the physical criterion [START_REF] Bouharguane | Multiphysics optimal transportation and image analysis[END_REF], has been introduced. In addition, a generalized metric interpolating the optimal transport and Fisher-Rao was proposed and found to be beneficial for different mass transportations [START_REF] Chizat | An interpolating distance between optimal transport and Fisher-Rao metrics[END_REF]. Recent techniques have demonstrated impressive performance in geophysical inverse problems such as full waveform inversion (FWI) [START_REF] Métivier | An optimal transport approach for seismic tomography: application to 3D full waveform inversion Inverse Problems[END_REF] and in the field of image processing: image registration [START_REF] Haker | Optimal mass transport for registration and warping[END_REF], video reconstruction [START_REF] Delon | Movie and video scale-time equalization application to flicker reduction[END_REF], color transfer [START_REF] Papadakis | A variational model for histogram transfer of color images[END_REF], image segmentation [START_REF] Ni | Local histogram based segmentation using the Wasserstein distance[END_REF].

Motivated by the advantages of the Wasserstein metric in the optimal transport theory, in [START_REF] Farchi | Using the Wasserstein distance to compare fields of polltants: application to the radionuclide atmospheric dispersion of the Fukushima-Daiichi accident[END_REF], it was used to compare the air pollutant dispersion from observation and model simulation. This non-local metric can avoid the "double penalty" effect in classical point-wise indicators when there is a misplacement. At present, a learningtype Wasserstein metric [START_REF] Gulrajani | Improved training of Wasserstein gans[END_REF] has been introduced for noise reduction [START_REF] Yang | Low dose CT image denoising using a generative adversarial network with Wasserstein distance and perceptual loss[END_REF]. The model error in DA based on the Wasserstein measure was investigated in [START_REF] Ning | Coping with model error in variational data assimilation using optimal mass transport[END_REF]. Viewing the 2-Wasserstein metric as a displacement interpolation [START_REF] Mccann | A convexity principle for interacting gases[END_REF][START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF], recently, an optimal transport-based DA [START_REF] Feyeux | Transport optimal pour l'assimilation de données images[END_REF][START_REF] Feyeux | Optimal transport for variational data assimilation Nonlinear Processes in[END_REF] considering the position error in the observation and background was proposed. The information from the observation and background was sufficiently used and promising results for shape preservation and position correction were obtained.

Owing to its tremendous effect on real weather and pollution forecasting (e.g., the scope), generally, the position and observation error are considered for DA. Moreover, the concentration observation is limited or even impossible to obtain (e.g., oil spill). To deal with these two key problems, in this work, our main contribution is the utilization of the Wasserstein measure to misfit the discrepancy and a combination with levelset method in the DA framework. First, a distance regularization-based method is established. Then, based on the Wasserstein metric instead of the Euclidean measure, two novel methods called original and splitting topological data assimilation (OTDA and STDA) are proposed. The proposed methods are suitable for an event lacking concentration information and different from previous approaches of dealing with the position error such as performing a concentration information alignment step before DA [START_REF] Hoffman | A technique for assimilating SSM/I observations of marine atmospheric storms: tests with ECMWF analyses[END_REF][START_REF] Ravela | Data assimilation by field alignment[END_REF][START_REF] Nehrkorn | Correcting for position errors in variational data assimilation[END_REF]. In new methods, contours in observations can be assimilated for dispersion prediction and simultaneously, the shape and position information between contours from observations and background can be integrated when there are position and shape errors (e.g., partially missing observation because of cloud cover and a displaced remote-sensing observer). The proposed method reduces the non-Gaussian type of error, position error, to some extent and optimizes the shape.

The remainder of the article is organized as follows. Section 2 presents a brief introduction on the classical and distance regularized level-set-based DA, including the contour data-fitting cost function and gradient. Section 3 details the optimal transport theory, Wasserstein distance, and topological data assimilation (OTDA and STDA) using the Wasserstein distance. The numerical experiment framework and results of the application of the novel level-set-based DA approach to a model test problem are discussed in Section 4. Finally, the main results are summarized and future work is outlined in Section 5.

Level-set-based assimilation using Euclidean distance

Level-set function and model

A given shape (e.g., pollutant spot C) can be represented by sub-domain Ω c ⊂ R 2 , whose boundary ∂Ω c can be defined by the zero level set of a mapping called level-set

function φ : R 2 → R ∀x ∈ R 2 ,        φ(x ) < 0, x ∈ Ω c φ(x ) = 0, x ∈ ∂Ω c φ(x ) > 0, x / ∈ Ωc (2.1)
where function φ(x ) is generally the signed distance from x to ∂Ω c = {x ∈ Ω|φ(x ) = 0}.

The motion of interface ∂Ω c (t) over time is considered as the evolution of level-set function φ, which is modeled by a two-dimensional (2D) nonlinear partial differential equation called the level-set model expressed as

∂φ ∂t + V • ∇φ + a ∇φ -σκ ∇φ = 0, on Ω × [0, T ] φ(0) = φ 0 , (2.2) κ = div( ∇φ ∇φ ) = ( ∂φ ∂x ) 2 ∂ 2 φ ∂y 2 -2 ∂φ ∂x ∂φ ∂y ∂ 2 φ ∂x∂y + ( ∂φ ∂y ) 2 ∂ 2 φ ∂x 2 ∇φ 3 (2.3)
where φ 0 is the initial value, a signed distance function (SDF) of contour ∂Ω c (0). ∂Ω c denotes the boundary of the polluted domain, ∇ is the gradient operator in the spatial domain, V = (u, v) represents the external velocity field from the flow, constant values a and σ denote the diffusion coefficients, and κ is the mean curvature. The modeling domain is Ω and simulation time is equal to T .

In the level-set model, level-set function φ(t) is driven by external velocity field V, and over time, the diffusion is along the normal direction by constant value a and in proportion to mean curvature κ.

Regarding the physics in the level-set model, in the presence of term σκ ∇φ , the closed curves of an arbitrary shape will collapse into a circle and then disappear, whereas under term a ∇φ , a relatively closed smooth curve will contain the sharp edges, and the topological structures of the curve may change. Both diffusion coefficients σ and a depend on the physical characteristics of the pollutant. The numerical periodic or homogenous Neumann boundary condition can be used for this model. To maintain the stability of the level-set function in the numerical simulation, the re-initialization approach [START_REF] Sethian | Level set methods for fluid interfaces[END_REF] can be adopted. Fig. 1 gives an example of the level-set method applied to two spots of pollutant.

Contour-fitting observation operators

To assimilate the active contours from the observation, two observation operators suitable to both metrics, namely the Euclidean metric and Wasserstein distance, are introduced. These types of operators can combine level-set function φ with the corresponding contours so that the Euclidean or Wasserstein distance can be used to efficiently measure the discrepancy between them. Let us define the "observation to structure" operator as

t ∈ [0, T ], C S (x , t) = H O→S (C obs (x , t)) = 1, x ∈ Ωc (t) 0, x ∈ Ω\ Ωc (t) (2.4)
where C S (t) belongs to the new observation space.

Another observation operator called the "model to structure" operator is defined as where φ denotes the level set function. H φ→S (φ) is a Heaviside-like function, which approximates to the following differentiable form:

H φ→S (φ) = 1, φ ≤ 0 0, φ > 0 (2.5)
H φ→S (φ) = 1 2 (1 -sgn(φ)) ≈ 1 2 (1 - φ φ 2 + ε ) (2.6)
where ε is a small constant. This operator maps level-set function φ to the new observation space.

Contour-fitting cost function with Euclidean distance

Generally, the optimal control theory is involved in the framework of traditional fourdimensional variational data assimilation (4D-Var). The cost function estimates the residual between the control variable from the governing model and observation over time. For the level-set model, the control variable is initial condition φ 0 . In [START_REF] Li | A level-set-based image assimilation method: potential applications for predicting the movement of oil spills[END_REF], the optimal control problem aimed at obtaining the minimizer of the following cost function:

J(φ 0 ) = J o +J b = δ o 2 N i=1 sgn[G i (φ 0 )]-C S (t i ) 2 L 2 (Ω) + δ b 2 φ 0 -φ b 0 2 L 2 (Ω) (2.7)
where G i (•) denotes the forward operator of the level-set model from time 0 to t i and sgn is a smoothing sign operator. C S is the corresponding structural image from the observation, and φ b 0 represents the background, which is an SDF. The zero level set of φ b 0 is the contour of the first guess. When weight δ b is small, the initial value is corrected by the observation. Due to the inconsistent way of comparison, information contained in both the observation and background cannot be used efficiently. In this case, the cost function can be modified by a mass conservation-type functional,

J(φ 0 , a, σ) = J o + J b + J a reg = δo 2 N i=1 N (H φ→S [G i (φ 0 , a, σ)]) -N (C S (t i )) 2 L 2 (Ω) + δ b 2 N (H φ→S [φ 0 ]) -N (B S ) 2 L 2 (Ω) + J init (φ 0 ) (2.8)
where

B S = H φ→S (φ b 0 ) (2.9)
and

J init (φ 0 ) = λ Ω R(φ 0 ), R(φ 0 ) = 1 2 (|∇φ 0 | -1) 2 , |∇φ 0 | ≥ 1 1 (2π) 2 (1 -cos(2π|∇φ 0 |)), |∇φ 0 | < 1 (2.10) H φ→S [•] is defined as (2.6)
, and B S is the structure from the background. N (•) is the normalization operator defined as

N (f ) = f Ω f
. The minimum problem is ill-posed when only the contour information from the observation and background is assimilated.

Therefore, a stable function called distance-based regularization term J a reg used in image segmentation [START_REF] Li | Distance regularized level set evolution and its application to image segmentation[END_REF] is added. This term contains an SDF-fitting function making φ 0 an SDF. λ denotes the weight.

In the following, we will focus on the computation of the gradient of the above optimization problem. The gradient of cost function J o with respect to the initial value can be obtained by solving an adjoint equation backward. For the term related to the background and distance-based regularization, the gradient is, respectively, expressed as

gradJ b (φ 0 ) = δ b (H φ→S ) * (H φ→S (φ 0 ) -B S ) (2.11) and gradJ a reg (φ 0 ) = -λdiv((1 -1 |∇φ 0 | )∇φ 0 ), |∇φ 0 | ≥ 1 -λdiv(sinc(2|∇φ 0 |)∇φ 0 ), |∇φ 0 | < 1 (2.12)
where * denotes the adjoint of a linear operator and sinc(x) = sin(2πx) 2πx .

Topological data assimilation

Optimal transport theory and Wasserstein distance

In [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovitch mass transfer problem[END_REF], a fluid mechanics formulation of the optimal transport was introduced with the objective of obtaining optimal geodesic ρ(t, x) that changed over time between two densities. The mechanics formulation was described as determining the minimizer of the following kinetic energy:

[0,1]×Ω ρ(t, x)|ν(t, x)| 2 dtdx (3.13)
where densities ρ(t, x) should belong to mass function space P(Ω) defined as

P(Ω) := {ρ ≥ 0 : Ω ρ(x)dx = 1} (3.14)
and the density and velocity (ρ, ν) satisfy a continuity equation,

C(ρ 0 , ρ 1 ) := { ∂ t ρ + div(ρν) = 0, (ρ, ν) s.t. ρ(t = 0) = ρ 0 , ρ(t = 1) = ρ 1 , ρν • n = 0 on ∂Ω } (3.15)
ρ 0 and ρ 1 are given as the densities at the initial and final moment, respectively. This type of optimal transport satisfies mass conservation. At present, there are numerous numerical methods for solving the optimal transport problem [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF][START_REF] Li | Computations of optimal transport distance with Fisher information regularization[END_REF].

Wasserstein distance W between mass functions ρ 0 and ρ 1 is defined as

W(ρ 0 , ρ 1 ) = min (ρ,ν)∈C(ρ 0 ,ρ 1 ) [0,1]×Ω ρ(t, x)|ν(t, x)| 2 dtdx (3.16)
Wasserstein metric W can be considered as a displacement interpolation between two densities functions. Hence, in the next section, we will mainly focus on the level-setbased DA with the Wasserstein metric-type contour-fitting functional to tackle the effect of the position and shape errors in the data.

Topological data assimilation using Wasserstein distance (TDA)

Both position and shape errors in the observation and background affect the quality of the assimilation. Motivated by the optimal transport theory, the Wasserstein distance could be used to simulate the residual instead of the l 2 -norm to make full use of their information and optimize the shape. Thus, the contour-fitting functional can be modified by

J W (φ 0 ) = J o W + J b W + J a reg = δo 2 N i=1 W(N (H φ→S [G i (φ 0 )]), N (C S (t i ))) 2 + δ b 2 W(N (H φ→S [φ 0 ]), N (B S )) 2 + J init (φ 0 ) (3.17)
where W(•, •) represents the Wasserstein distance. The normalization operator makes the mass of transport equal to 1, i.e., their masses are the same.

N (H φ→S [G i (φ 0 )]), N (C S (t i )), N (H φ→S [φ 0 ]
) , and N (B S ) belong to mass function space P(Ω), which satisfies the optimal transport theory. Alternatively, the Fisher-Rao metric-based approach [START_REF] Chizat | An interpolating distance between optimal transport and Fisher-Rao metrics[END_REF] could be used to solve the problem of different masses in optimal transport.

δ o and δ b are the weights and N is the number of observations. J init (φ 0 ) is the SDFfitting term.

To obtain the minimizer of the variational model (3.17), we propose splitting method-based topological data assimilation (STDA), which is summarized in Algorithm 3.2. The cost function is split into two parts, with one term being associated with the observation and background (J o,b W ) and the other being SDF-fitting term J init (φ 0 ). First, the gradient of J o,b W is obtained and the solution is updated by the descent algorithm. Then, the updated solution is re-initialized to an SDF, taken as the optimal solution in the next iteration. To accelerate the convergence, in this work, the re-initialization process [START_REF] Sethian | Level set methods for fluid interfaces[END_REF] is performed by

∂ψ ∂τ = sgn(φ 0 )(1 -∇ψ ), on Ω × [0, τ 0 ] ψ(0) = φ 0 (3.18)
As time t → ∞, ψ gradually satisfies ∇ψ = 1, which approximates an SDF.

Alternatively, the TDA method based on the SDF-fitting function in (2.8) is called the original approach (OTDA), which is achieved by Algorithm 3.2. Next, we concentrate on the computation of the gradient of new cost function J o,b W with respect to initial value φ 0 .

In [START_REF] Feyeux | Transport optimal pour l'assimilation de données images[END_REF], the Wasserstein distance-based variational model proved to be differentiable, and by choosing different representations of the inner products and metric methods in the tangent space, the gradient of the Wasserstein-based cost function changed according to the definition of the Gâteaux derivative expressed as

∀η ∈ T φ 0 P, lim →0 J o,b W (φ 0 + η) -J o,b W (φ 0 ) = grad J o,b W (φ 0 ), η (3.19)
where •, • represents a type of inner product and grad represents the gradient of the cost function using such inner product. T φ 0 P is called the tangent space. The following theorem provides the gradient of the proposed cost function.

Theorem 3.1 If the Euclidean scalar product and l 2 -norm are used in tangent space

T φ 0 P, then the gradient of J o,b
W contained in the objective function (3.17) with respect to initial value φ 0 becomes

grad 2 J o,b W (φ 0 ) = δ o N i=1 (G i ) * (H φ→S ) * N * Ψ i + δ b (H φ→S ) * N * Ψ b (3.20)
where Then, nonlinear descent algorithms such as the gradient descent and conjugate gradientlike method can be used to obtain the optimal solution of the cost function. To maintain the stability of the optimal value, at each iteration, the re-initialization method should be adopted to make initial contour φ 0 an SDF. The re-initialization method can accelerate the convergence during the shape optimization (Remark 3.4). A schematic of proposed methods is presented in Fig. 2.

Ψ i , i = 1, • • • , N and Ψ b are the Kantorovich potentials of the transport between N (H φ→S [G i (φ 0 )]) and N (C S (t i )), N (H φ→S [φ 0 ]) and N (B S ), respectively. G i ,
Remark 3.2 For the optimal transport problem (3.13), the optimal velocity fields ν(t, x)

have the property [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovitch mass transfer problem[END_REF][START_REF] Feyeux | Transport optimal pour l'assimilation de données images[END_REF][START_REF] Feyeux | Optimal transport for variational data assimilation Nonlinear Processes in[END_REF],

ν(t, x) = ∇Φ(t, x)
where the function Φ(t, x) follows a Hamilton-Jacobi equation

∂Φ ∂t + |∇Φ| 2 2 = 0
In addition, the so-called Kantorovich potential Ψ(x) is defined by

Ψ(x) := -Φ(t = 0, x)
So the following representation can be obtained

ν(t = 0, x) = -∇Ψ(x)
where ν(0) is the initial velocity fields for the minimizer of (3.13).

Therefore, Kantorovich potentials Ψ i and Ψ b can be determined by solving a series of 2D Poisson equations

∆Ψ i = -div(ν(0) i ),
∂Ψ i ∂n = 0 on ∂Ω where ν(0) i and ν(0) b are the initial velocities in the optimal transport problems, in which the available mass functions are described by

ρ i 1 = N (H φ→S [G i (φ 0 )]) and ρ i 2 = N (C S (t i )), i = 1, • • • , N , ρ b 1 = N (H φ→S [φ 0 ]) and ρ b 2 = N (B S
). ν(0) i and ν(0) b can be obtained by the iterative algorithm in [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF].

Remark 3.3 For normalization operator N (•), the following condition is well satisfied:

N * (Ψ + C) = N * Ψ
where C belongs to the constant space. Thus, the effect of the non-unique solution of the above Poisson equation can be avoided. A brief illustration is given in the Appendix.

Remark 3.4 In [START_REF] Feyeux | Transport optimal pour l'assimilation de données images[END_REF], two types of gradients with respect to the initial value were defined by using different inner products in tangent space T φ 0 P, called as the L 2 gradient and W gradient. W gradient shows a better performance than L 2 gradient for convergence but under the assumption of a non-flux at the boundary of the L 2 gradient. Our proposed method also has a rapid convergence, mainly owing to the re-initialization process when the L 2 gradient is utilized. Algorithm 3.2 : Splitting method-based TDA Algorithm (STDA)

Input: Observations C obs (t i ), i = 1, • • • , N , background φ 0 0 = φ b 0 , weight parameters δ o and δ b , stopping criterion τ . Output: Analysis φ a 0 . 1: while |J(φ k 0 ) -J(φ k-1 0 )| > τ do 2: Compute the gradient of J o,b W with respect to φ k 0 : grad 2 J o,b W (φ k 0 ) = δ o N i=1 (G i ) * (H φ→S ) * N * Ψ i k + δ b (H φ→S ) * N * Ψ b k . 3: Iterate φ k
0 by the descent algorithm (e.g., steepest descent or conjugate gradient method)

φ k+1 0 = φ k 0 -µ k Θ(grad 2 J o,b W (φ k 0 )
) where µ k is the optimal step size in the k iteration and Θ(grad

2 J o,b W (φ k 0 )) represents a kind of combination of grad 2 J o,b W (φ k 0 ). 4: Re-initialize φ k+1
: grad 2 J W (φ k 0 ) = δ o N i=1 (G i ) * (H φ→S ) * N * Ψ i k + δ b (H φ→S ) * N * Ψ b k + gradJ init (φ k 0 )
, where gradJ init (•) is given in (2.12).

3: Iterate φ k

0 by the descent algorithm (e.g., steepest descent or conjugate gradient method)

φ k+1 0 = φ k 0 -µ k Θ(grad 2 J W (φ k 0 )
) where µ k is the optimal step size in the k iteration and Θ(grad 2 J W (φ k 0 )) represents a kind of combination of grad 2 J W (φ k 0 ). 4: k = k + 1.

5: end while

Numerical examples

Numerical setting

To perform the numerical simulation, the finite difference scheme is adopted for the forward level-set model (2.2). The upwind scheme is used for discretizing the convection term, and the central difference method is employed for the diffusion terms. For the The total simulation time is equal to 6s. The assimilation window is [0, T ] with T = 3s, and the moment of prediction is 6s. Heaviside-like parameters ε used in (2.8) and (3.17) are defined as (∆x) 2 and (2∆x) 2 , respectively. The model (2.2) has the following parameters: a = 0.05 and σ = 1e -5. The velocities are known in the test

with u = v = 1ms -1 .
The control variable is initial value φ 0 only. The observations are generated by the forward level-set model using a given initial value φ obs 0 , whose zero level set is shown as a pink solid line and named as Observation in Figs. 4a, 4d, and4g

C obs (t) = C(t), C(t, x ) > δ 0, C(t, x ) ≤ δ (4.24)
here, δ represents the threshold above which the spots of the pollutants are visible. In this test, δ is set to 0.2. The contour can be extracted from the observation by the wavelet-based method [START_REF] Mallat | Zero-crossings of a wavelet transform[END_REF].

The total simulation time is equal to 20s. The assimilation window is [0, T ], with T = 8.5s. Heaviside-like parameters ε used in (2.8) and (3.17 Remark 4.1 In a real event, first, unknown diffusion parameters a and σ in the levelset model have to be identified using available observation data, which highly depend on the characteristics of the corresponding pollutants. Then using the identified parameter couples, the first guess or background, and observed contours, the initial value can be optimized via our proposed method.

Numerical test and results

In this section, the results of the comparison between the ED method and Wasserstein distance-based method (WD method) are presented. By using normalization operator

N (f ) = f Ω f
, the mass of the observation in one moment and background will remain unchanged, satisfying the theory of optimal transport. Motivated by this, the proposed WD method can deal with the problem that the observation and background not only have position errors but also include different shapes. First, we show the results of test cases where the observation sequence is acquired from the level-set method and the parameters in the first group of DA test described in Section 4.1 are used. Let us consider the optimal solution of the following cost function (time-independent):

J W (φ) = δx 2 W(N (H φ→S [φ]), N (H φ→S [φ x ])) 2 + δy 2 W(N (H φ→S [φ]), N (H φ→S [φ y ])) 2 + J init (φ) (4.25)
where φ x and φ y are the given level-set functions that have the same shape but different positions (shown in Fig. 3a). δ x = δ y = 1e + 3. Zero level sets {x ∈ Ω|φ(x , t) = 0} of optimal value φ obtained by using the two methods are displayed in Fig. 3b. It can be seen that the WD method can maintain the shape, which is different from the results obtained by the ED method. The changes in the cost function and l 2 -norm of the gradient with increasing number of iterations are presented in Figs. 3c and3d. The fast convergence shows the efficiency of the proposed method.

For the time-dependent data assimilation problem (3.17 Here, we assume that the background is in the right position and observation provides the right shape information. The zero level set of the optimal value called analysis φ a 0 with weights couples a and b (Table 1) are shown in Figs. 4b and4c respectively (ED method shown as blue dashed lines and named as Analysis-E; WD method shown as red solid lines and named as Analysis-W). The WD method makes full use of the information between the observation and background, which can be illustrated as the optimal interpolation by using Wasserstein distance W(•, •) when the background and observation contain position errors. However, compared with the WD method, the ED method, only the linear interpolation is used in by the l 2 -norm. Once position errors occur, a larger shape will be obtained, the so-called "double penalty" effect, where one effect is from the observation and the other from the background. Therefore, compared with the ED method, the proposed WD method provides position and shape correction for the results only from observation and the relatively accurate shape information when weights δ o and δ b are uncertain. In this regard, the ED method will be ineffective. In addition, with weight δ o δ b , the results obtained by both the methods are similar and are closer to those obtained only from the observation (presented in Fig. 4c).

δ b N b 1e+3 30
25δ b N Test 2 c 1e+4 3e+2 2δ b N d 1e+3 30
25δ b N Test 3 e 1e+4 3e+2 2δ b N f 1e+3 30
25δ b N
In order to further verify the validity of the proposed approach, another two tests (Test 2 and 3) have been done in which both the observation and background simulation contain complex multiple structures of interest (two and three spots, respectively).

Figs. 4d and 4g display the boundaries of the initial observation and zero level sets of the background (shown as pink solid lines and black dashed lines). Similarly, some reasonable results of the test cases are presented in Figs. 4e and4f with weights couples c and d (Table 1), and Figs. 4h and 4i with weights couples e and f (Table 1). As shown, the WD method preserves the shape of the contours compared with the ED method. Because the weights δ o = 2δ b N used in the test is a little larger than that in test 1 (δ o = δ b N ), the position of the contours could be closer to the true state rather than in the middle (shown in Figs. 4e and4h). But from the results it can be seen that the proposed method will smooth the edge of the contours and if two contours are too close to each other, they will merge together easily (shown in Figs. 4h and4i). It may lie in the regularization method, SDF-fitting method and the inaccurate computation for the minimizer of the cost function. Nevertheless, it can still indicate that the novel method applies equally to multiple structures of interest problems. In addition, for the multiple structures of interest problem, numerical tests show that SDF-fitting function

J init (•) in (2.8
) is used for the stable convergence although it is slower. By using identified parameter a and a given background, the tests for assimilating the contours of the pollutants from the full concentration observation images and background with position and shape errors and three weights couples g-i (Table 2) are implemented, the results of which are presented in Figs. 6b-6d. For comparison, Fig. 6a displays the contours of the pollutants in the first observation image (shown as a pink solid line), which can be missing in reality. Similar to the test above, under the assumption that the background is in the right position and observation provides the right shape information, an obvious modification in the position and optimal shape can be found in the results of the analyses with weights couples g (Table 2) as well as in the results of prediction at moment t = 10s (shown in Fig. 7).

Another test is performed with a missing observation in the spatial domain (shown in Fig. 8a). Differing from previous assumptions, in this test, we assume that the missing observation is in the right position, and shape of the first guess is correctly simulated but has the wrong position. The analyses obtained by both methods with weights couples j (Table 2) are compared in Fig. 8c. It can be concluded that using the proposed method optimally integrates the shape and position information and makes full use of them when the missing observation is available. The forecast at moment t = 20s (shown as a red solid line in Fig. 8e) is closer to that by using the true state, which can assist in optimizing rescue operations for real accidents to some extent.

Conclusion

In this paper, we proposed a novel data assimilation strategy, namely topological data assimilation, which integrates the level-set method with the Wasserstein metric in an optimal way. Moreover, several numerical tests are presented. First, a distance regularization-based approach is introduced. Then based on it, the Wasserstein distance is used to measure the discrepancy between the model variables, observation, and background, rather than using the Euclidean distance. Combined with the level-set method, the contour information can be assimilated, which is suitable for the observation provided with only the edges of the pollutants, such as an oil spill. Once the misplaced missing observation and background are available and there is uncertainty in the errors, the proposed approach considered as a type of shape optimization can overcome the "double penalty" effect in the local Euclidean metric method. An analysis after the optimal shape integration and modification in the position between the observation and background is obtained. The new method is particularly capable of dealing with obstacles (e.g., cloud cover) in a remote-sensing observer, which occur in a real pollutant forecast. However, the results show that the proposed method will smooth the analysis.

In future, we will consider some edge-preserving regularization methods, SDF-fitting methods and extend the proposed method to possible real applications such as oil spill, wildfire propagation, tumor growth, and/or typhoons.

Acknowledgement

The authors would like to thank Dr.

Appendix

In this section, we provide a brief proof of Theorem 3.1 and mathematical illustration of Remark 3.3.

Proof 7.1 Here, we only give the short proof for the case that is not included in [START_REF] Feyeux | Transport optimal pour l'assimilation de données images[END_REF].

Owing to the use of normalization operator N (•) and observation operators H φ→S (•), the following condition can be satisfied simultaneously

∀ small ∈ R, η ∈ T φ 0 P, M i (φ 0 ) ∈ P(Ω) and M i (φ 0 + η) ∈ P(Ω)
where

M i (•) = N (H φ→S [G i (•)]
). We have

J o,b W (φ 0 + η) = δo 2 N i=1 W(M i (φ 0 + η), N (C S (t i ))) 2 + δ b 2 W(T (φ 0 + η), N (B S )) 2 = δo 2 N i=1 W(M i (φ 0 ) + M i [φ 0 ]η + o( ), N (C S (t i ))) 2 + δ b 2 W(T (φ 0 ) + T[φ 0 ]η + o( ), N (B S )) 2
where M i and T are the Jacobian of the operators M i (•) and T (•), which are NH φ→S G i and NH φ→S . Based on the theory (Theorem 8.13 in [START_REF] Villani | Topics in optimal transportation[END_REF]), we can obtain

1 2 W(M i (φ 0 ) + M i [φ 0 ]η + o( ), N (C S (t i ))) 2 = 1 2 W(M(φ 0 ), N (C S (t i ))) 2 + η, Ψ i + o( )
as well as the term regarding the background, where Ψ i denotes the Kantorovich potential of the transport between N (H φ→S [G i (φ 0 )]) and N (C S (t i )). Hence,

J o,b W (φ 0 + η) = δo 2 N i=1 W(M i (φ 0 ), N (C S (t i ))) 2 + δ o N i=1 Ψ i , M i [φ 0 ]η + δ b 2 W(T (φ 0 ), N (B S )) 2 + δ b Ψ b , T[φ 0 ]η + o( )
Finally, according to the definition of the Gâteaux derivative, the gradient is obtained as follows

grad 2 J o,b W (φ 0 ) = δ o N i=1 (G i ) * (H φ→S ) * N * Ψ i + δ b (H φ→S ) * N * Ψ b
where * denotes the transpose of the linear operator.

For Remark 3.3, first, we give the derivative of the normalization operator.

Suppose f ∈ L 2 (Ω), then according to definition N (f ) = f Ω f , we have

∀ > 0, η ∈ L 2 (Ω), N (f + η) = f + η Ω f + η = η Ω f + f Ω f - f Ω η ( Ω f ) 2 + o( )
Based on the definition of the derivative of operator G(•) described by

G(f ), η = lim →0 G(f + η) -G(f )
the following can be obtained

N(f ), η = η Ω f -f Ω η ( Ω f ) 2
Next, we compute the continuous adjoint operator of N. Because

Nη, ϕ = Ω Nηϕ = Ω ηm-f Ω η m 2 ϕ = 1 m Ω ηϕ -Ω f ϕ( Ω η) m 2 = 1 m Ω ηϕ -Ω η( Ω f ϕ m 2 ) = η, ϕ m -Ω f ϕ m 2
where m = Ω f , we have,

N * ϕ = ϕ Ω f - 1 ( Ω f ) 2 Ω f ϕ Thus, if C belongs to constant space, N * C = 0 8. Reference
o = δ b N , δ o = 25δ b N , δ o = 2δ b N , δ o = 25δ b N , δ o = 2δ b N and δ o = 25δ b
N (ED method shown as a blue dashed line and named as Analysis-E; WD method shown as a red solid line and named as Analysis-W). The results of prediction using the corresponding analyses.

Figure 1 :

 1 Figure 1: Level-set method. (a): Concentration distribution of the pollutant spots at one moment; (b): Boundary of pollutant spots ∂Ω c ; (c): Level-set function φ (signed distance function based on ∂Ω c where ∂Ω c = {x ∈ Ω|φ(x ) = 0}).

  H φ→S and N are the Jacobian of operators G i (•), H φ→S [•], and N (•). The definitions of the Euclidean scalar product and l 2 -norm are described by ∀ξ, ζ ∈ T φ 0 P, ξ, ζ 2 = Ω ξζdx and ξ 2 L 2 (Ω) = Ω |ξ| 2 dx (3.21)The brief proof is given in the Appendix. The method of computing the Kantorovich potential numerically used in the gradient (3.20) is illustrated in Remark 3.2 and 3.3.

0

  to a signed distance function φ init k 0 and set φ init k+1 0 as φ k+1 0 . 5: k = k + 1. 6: end while Algorithm 3.2 : Original TDA Algorithm (OTDA) Input: Observations C obs (t i ), i = 1, • • • , N , background φ 0 0 = φ b 0 (a signed distance function), weight parameters δ o and δ b , stopping criterion τ . Output: Analysis φ a 0 . 1: while |J(φ k 0 ) -J(φ k-1 0 )| > τ do 2: Compute the gradient of J W with respect to φ k 0

Figure 2 :

 2 Figure 2: Flowchart of topological data assimilation (OTDA: Original TDA; STDA: Splitting method-based TDA).

  . In the Euclidean distance-based (ED) method and original topology data assimilation (OTDA), the number of observations is N = 5 and weight is λ = 0.3. The zero level sets of the corresponding backgrounds with position or shape errors are shown as black dashed lines and named as Background in Figs. 4a, 4d, and 4g. Time τ 0 in the re-initialization (3.18) is set to 9.5s. The nonlinear conjugate gradient algorithm is iterated to search for the optimal solution. Parameters in the second group of DA test: A simple 2D concentration transport model is used to test the efficiency of the proposed method. The model is described by  • ∇C = µ∆C, on Ω × [0, T ] C(0) = C 0 , ∂C ∂n = 0 on ∂Ω (4.22) where C 0 denotes the concentration of the initial pollutant spots. V = (u, v) represents the velocity, which is the same as that in the model (2.2), and µ is the diffusion coefficient. The homogenous Neumann boundary condition is used. In the test, µ = 0.01 and V = (u, v) is obtained by a 2D Euler-type fluid dynamics model. The general form currents of the 2D turbulence are displayed in Fig. 5a in which the vortex can be seen clearly. Numerical solutions C(t) after thresholding are taken as the observation images of the pollutants (shown partly in Fig. 5a), which can be expressed as Ω c (t, x ) = {x (t) ∈ Ω|C(t, x ) > δ}, (4.23)

  ) are defined as (1.3∆x) 2 . Parameters σ and a are unknown depending on the pollutants and need to be retrieved in the model (2.2), which may occur in a real event. In our examples, σ is known (σ = 0) for simplicity. The control variables are parameter a and initial value φ 0 . First, parameter a should be identified with the weight couple (δ o = 30, δ b = λ = 0) by using the variational model (2.8) under the assumption that a reliable background is provided. Number of concentration observations N = 9. For the initial value inversion, background with position or shape errors is given, whose zero level sets are shown as black dashed lines and named as Background in Figs. 6a and 8b. Here, in the ED method, number of concentration observations N = 5 and weight λ = 0.3.

  ), two groups of DA test cases are formed. 4.2.1. The first group of DA test These are three tests (Test 1-3) in the first group of DA test. In test 1, observation with one simple spot and background with multiple interests are provided. The boundaries of the initial observation and zero level sets of the background are shown in Fig. 4a (shown as pink solid lines and black dashed lines).

4. 2

 2 .2. The second group of DA test For a more realistic experiment, the full and partially missing concentration observation from the 2D transport model and currents from the Euler equation are provided to test the efficiency of the proposed method. Here the parameters in the second group of DA test described in Section 4.1 are used. The concentration images and currents at four times are presented in Fig. 5a. Before the correction of the initial contours, constant diffusion parameter a in the level-set model is retrieved by using N = 9 concentration observation images C obs , with σ set to 0 for simplicity. In reality, the diffusion parameters depend on the pollutant. The change in a with the iteration number is shown in Fig. 5b, with first guess a = 0. It converges to a = 0.1507 after 10 iterations in this test.

[ 1 ]

 1 Sasaki Y 1970 Some Basic formalisms in numerical variational analysis Monthly Weather Review 98 875-883 [2] Le Dimet F X and Talagrand O 1986 Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects Tellus A: Dynamic Meteorology and Oceanography 38 97-110 [3] Ghil M, Cohn S, Tavantzis J, Bube K and Isaacson E 1981 Applications of estimation theory to numerical weather prediction In Dynamic Meteorology: Data Assimilation Methods Springer 139-224 [4] Evensen G 1994 Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte

Figure 3 :

 3 Figure 3: The results of the proposed WD method (STDA) compared with ED method, with observations and background containing position and shape errors. (a): Boundary of the initial observation and zero level set of the background (shown as a pink solid line and black dashed line); (b): Zero level sets of the corresponding analyses (ED method shown as a blue dashed line and named as Analysis-E; WD method shown as a red solid line and named as Analysis-W); (c): Cost function J W ; (d): l 2 -norm of the gradient.

Figure 4 :

 4 Figure 4: The results of the proposed WD method (OTDA) compared with ED method, with observations and background of multiple structures containing position and shape errors. (a), (d) and (g): Boundary of the initial observation and zero level set of the background (shown as a pink solid line and black dashed line); (b), (c), (e), (f), (h) and (i): Zero level sets of the corresponding analyses with weights δo = δ b N , δ o = 25δ b N , δ o = 2δ b N , δ o = 25δ b N , δ o = 2δ b N and δ o = 25δ bN (ED method shown as a blue dashed line and named as Analysis-E; WD method shown as a red solid line and named as Analysis-W).

Figure 5 :

 5 Figure 5: Identification of parameter a in the level-set model using the contour information in concentration observations. (a): Partial observation sequences C obs (t) from the 2D concentration model, currents and boundary of the pollutant (Pollutant spots transportation over time); (b): Retrieved parameter a in the level-set model using such data; (c): l 2 -norm of the gradient.

Figure 6 :

 6 Figure 6: The results of the proposed WD method (STDA) compared with ED method, with observations from the 2D concentration model containing position errors, under the assumption that the background is in the right position but has inaccurate shape and the observations are with accurate shape but in the wrong position. (a): Boundary of the initial observation and zero level set of the background (shown as a pink solid line and black dashed line); (b)-(d): Zero level sets of the corresponding analyses with weights δ o = 2δ b N , δ o = 15δ b N and δ o = 0.05δ b N(ED method shown as a blue dashed line and named as Analysis-E; WD method shown as a red solid line and named as Analysis-W).

Figure 7 :

 7 Figure 7: The results of prediction. (a): The results of prediction using the boundary of the initial observation and background (shown as a pink solid line and black dashed (b)-(d): The results of prediction using the corresponding analyses with weights δ o = 2δ b N , δ o = 15δ b N and δ o = 0.05δ b N(ED method shown as a blue dashed line and named as Analysis-E; WD method shown as a red solid line and named as Analysis-W).

Figure 8 :

 8 Figure 8: The results of the proposed WD method (STDA) compared with ED method, with partially missing observations (in spatial domain) and background containing position errors, under the assumption that the background is in the wrong position but has accurate shape. (a): Observation sequence covered by cloud; (b): Boundary of the true state (reference), boundary of the initial observation and zero level set of the background (shown as a pink solid line, black dashed line, and orange solid line); (c): Zero level sets of the corresponding analyses with weights δ o = 3δ b N (ED method shown as a blue dashed line and named as Analysis-E; WD method shown as a red solid line and named as Analysis-W); (d): The results of prediction using the boundary of the true state, boundary of the initial observation and background; (e):The results of prediction using the corresponding analyses.

Table 1 :

 1 Weights in the first group of DA test

	DA test case: No. of weight	Background	Background weight δ b Observation
	group 1	couples	weight δ b (ED)	(WD-OTDA)	weight δ o
	Test 1	a	1e+4	3e+2

Table 2 :

 2 Weights in the second group of DA test

	DA test case: No. of weight	Background	Background weight δ b Observation
	group 2	couples	weight δ b (ED)	(WD-STDA)	weight δ o
		g	1e+4	1e+3	2δ b N
	Test 4	h	2e+3	2e+2	15δ b N
		i	2.5e+4	2.5e+3	0.05δ b N
	Test 5	j	1e+4	1e+3	3δ b N

  Nelson Feyeux for his active discussion and helpful suggestions. The work is supported in part by National Key Research and Development Program of China under Grant 2017YFB0202902, and NSFC under Grant 41625017 and Grant 91730306, and China Scholarship Council.