Alexandre Janon

Maëlle Nodet

Clémentine Prieur

Goal-oriented error estimation for the reduced basis method, with application to sensitivity analysis

Keywords: reduced basis method, surrogate model, reduced order modelling, response surface method, scientific computation, sensitivity analysis, Sobol index computation, Monte-Carlo method AMS MSC: 65M15

The reduced basis method is a powerful model reduction technique designed to speed up the computation of multiple numerical solutions of parametrized partial differential equations. We consider a quantity of interest, which is a linear functional of the PDE solution. A new probabilistic error bound for the reduced model is proposed. It is efficiently and explicitly computable, and we show on different examples that this error bound is sharper than existing ones. We include application of our work to sensitivity analysis studies.

Introduction

A large number of mathematical models are based on partial differential equations (PDEs). These models require input data (e.g., the physical features of the considered system, the geometry of the domain, the external forces...) which enter in the PDE as parameters. In many applications (for instance, design optimization, data assimilation, or uncertainty quantification), one has to numerically compute the solution of a parametrized partial differential equation for a large number of values of the parameters. In such a case, it is generally interesting, in terms of computation time, to perform all possible parameterindependent computations in an offline phase, which is done only once, and to call an online phase for each required value of the parameter, during which the information gathered in the offline phase can be used to speed up the computation of an approximate solution of the PDE, and, hence, to reduce the marginal (ie., per parameter) computation cost. The reduced basis method [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF] is a way of specifying such offline and online phases, which has been successfully applied to various well-known PDEs [START_REF] Grepl | A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations[END_REF][START_REF] Janon | Certified reduced-basis solutions of viscous Burgers equations parametrized by initial and boundary values[END_REF][START_REF] Knezevic | Reduced basis approximation and a posteriori error estimation for the parametrized unsteady boussinesq equations[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: Rigorous reduced-basis a posteriori error bounds[END_REF]. One should note that, in the reduced basis (RB) method, the online phase does not compute a solution which is strictly identical to the numerical PDE solution, but an approximation of it, obtained by projecting the original discretized equations onto a well-chosen basis. In the application cases given above, however, one is not interested in the solution by itself, but rather in a quantity of interest, or model output, which is a functional of this solution. Taking this functional into account when performing the model reduction leads to a so-called goal-oriented method. For instance, goal-oriented basis choice procedures have been tried with success in the context of dynamical systems in [START_REF] Ilak | Modeling of transitional channel flow using balanced proper orthogonal decomposition[END_REF][START_REF] Willcox | Balanced model reduction via the proper orthogonal decomposition[END_REF], where the basis is chosen so as to contain the modes that are relevant to accurately represent the output of interest, and in a general context in [START_REF] Bui-Thanh | Goal-oriented, model-constrained optimization for reduction of large-scale systems[END_REF], where the basis is chosen so as to minimize the overall output error. All those papers showed that using an adapted basis could lead to a great improvement of reduction error. This paper is about goal-oriented error estimation, that is, the description of a rigorous and computable error bound between the model output and the reduced one. Two different reduced model outputs can be considered: the first (which we call uncorrected reduced output) is simply the output functional evaluated at the reduced output. The second (called corrected reduced output), described in [START_REF] Machiels | A general framework for finite element a posteriori error control: application to linear and nonlinear convection-dominated problems[END_REF][START_REF] Maday | A blackbox reduced-basis output bound method for noncoercive linear problems[END_REF][START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF], is the same, up to a correction term obtained from the solution of an auxiliary (dual) problem. The rate of convergence of the corrected output is better than the uncorrected one but the computation of the correction involves the application of the RB method to the dual problem, and this has generally the drawback of doubling offline and online computational times. Regarding output error estimation, an error bound for the difference between the corrected reduced output and the original output is provided in the papers cited above. In this paper, we propose two new goal-oriented error bounds: one for the uncorrected reduced output, and one for the corrected reduced output. We also show, in numerical examples, that our bound is more precise than the existing bound. This paper is organized as follows: in the first part, we describe our output error bounds and explain how to compute them; in the second part, we see how to apply our error bound to sensitivity analysis studies; finally, the third and fourth parts present numerical applications.

Methodology

Preliminaries

Reference problem We begin by setting up the context of the reduced basis method for affine-parametrized linear partial differential equations presented in [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF]. Our reference problem is the following: given a parameter tuple µ ∈ P ⊂ R p , and the vector space X = R N (for N ∈ N), find u(µ), the solution of:

A(µ)u(µ) = f (µ), (1)
where A(µ) is an invertible square matrix of dimension N , and f (µ) ∈ X, then compute the output :

s(µ) = S(u(µ)) (2)
where S : X → R is a linear form on X.

Choice of the inner product: We suppose that X = R N (with the standard basis), is endowed with the standard Euclidean inner product:

u, v = u t v, with associated norm u = u, u
Affine decomposition hypothesis: We suppose that A(µ) and f (µ) admit the following so-called affine decomposition [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF]:

∀µ ∈ P, A(µ) = Q q=1 Θ q (µ)A q , f (µ) = Q ′ q ′ =1 γ q ′ (µ)f q ′ (3)
where Q, Q ′ ∈ N * , Θ q : P → R and γ q ′ : P → R (for q = 1, . . . , Q, q ′ = 1, . . . , Q ′ ) are given functions, A q are square matrices of dimension dim X and f q ′ ∈ X. This hypothesis is required by the reduced basis method.

Reduced basis method

The dimension of the finite element subspace dim X is generally fairly large, so that the numerical computation of u(µ) from the inversion of A(µ) is expensive. The reduced basis aims at speeding up "many queries", that is, the computation of u(µ) for all parameters µ ∈ P 0 where P 0 is a finite but "large" subset of the parameter set P.

Reduced problem:

We consider a subspace X of X, and a matrix Z whose columns are the components of a basis of X in a basis of X. This basis of X is called the reduced basis in the sequel. We denote by u(µ) the components, in the reduced basis, of the solution of the projection of (1) onto X, that is, the solution of:

Z t AZ u(µ) = Z t f (µ) (4)
(where, for any matrix M , M t is the transpose of M ).

Choice of the reduced subspace: There are different techniques for choosing the reduced basis (the Z matrix). This paper does not focus on this topic, but we cite the POD method (Proper orthogonal decomposition) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part i-ii[END_REF], and the Greedy method [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF].

Offline-online decomposition:

The many-query computation can then be split into two parts: the first part (usually called the "offline phase"), which is done only once, begins by finding a reduced subspace, then the Q parameterindependent matrices:

A q = Z t A q Z, q = 1, . . . , Q and the Q ′ vectors: f q ′ = Z t f q ′ , q ′ = 1, . . . , Q ′
are computed and stored. In the second part (the "online phase"), we compute, for each value of the parameter µ:

A(µ) = Q q=1 Θ q (µ) A q , f (µ) = Q ′ q ′ =1 γ q (µ) f q ′ (5)
and solve for u(µ) satisfying:

A(µ) u(µ) = f (µ). (6)
The key point is that the operations in ( 5) and ( 6) are performed on vectors and matrices of size dim X, and that the complexity of these operations is totally independent from the dimension of the underlying "truth" subspace X. In many cases, the smoothness of the map µ → u(µ) allows to find (in a constructive way, ie., compute) X so that dim X ≪ dim X while keeping u(µ) -Z u(µ) small, hence enabling significant computational savings.

Output approximation: The output s(µ) can also be approximated from u(µ) using an efficient offline-online procedure: let l ∈ X be so that:

s(u(µ)) = l, u ∀u ∈ X;
in the offline phase we compute and store:

l = Z t l
and in the online phase we take:

s(µ) = l, u(µ)
as an approximation for s(µ).

Reduced-basis error bounds

Bound on u: Under additional coercivity hypothesis on A, the reduced basis method [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF] also provides an efficient offline-online procedure for computing ǫ u (µ) so that the approximation can be certified. This bound is based on the dual norm of the residual :

ρ(µ) = r(µ) RB = A(µ)Z u(µ) -f (µ) RB
where • RB is a suitably chosen norm on X (not necessarily • ), and a stability constant bound, which can be written as:

0 < α(µ) ≤ inf v∈X, v RB =1 v t A(µ)v (7)
when A is symmetric. The inequality sign is due to the fact that the exact infimum can be costly to evaluate in the online stage; usually a procedure such as the successive constraints method [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF] is used in order to quickly find a lower bound.

The bound ǫ u (µ) reads:

∀µ ∈ P u(µ) -Z u(µ) RB ≤ ρ(µ) α(µ) := ǫ u (µ),
The online procedure for the computation of ǫ(µ) is also of complexity independent of dim X.

Lipschitz bound on s: This online error bound can in turn be used to provide a certification on the output:

∀µ ∈ P |s(µ) -s(µ)| ≤ l RB ǫ u (µ) =:ǫ L (µ) (8)
We call this bound the "Lipschitz" bound, and denote it by ǫ L (µ). It is wellknown that this bound is very pessimistic.

The aim of Section 1.2 is to bound |s(µ)s(µ)| by a quantity which is smaller than ǫ L (µ) of ( 8) and can be computed using an efficient offline-online procedure which does not require computation of ǫ u (µ), described in Section 1.3. In Section 1.4, we consider a better approximation of s(µ) (denoted by s c (µ)) which also depends on the solution of the adjoint equation of (1) projected on a suitably selected dual reduced basis, and we see how the proposed bound for |s(µ)s(µ)| can be modified in order to bound |s(µ)s c (µ)|.

Probabilistic error bound

In this section, we give the expression of our output error bound. We recall the notation for the residual r(µ):

r(µ) = A(µ)Z u(µ) -f (µ) ∈ X,
and the adjoint problem solution (which will naturally appear in the proof of Theorem 1.1) by w(µ):

w(µ) = A(µ) -t l.
Let, for any orthonormal basis

Φ = {φ 1 , . . . , φ N } of X, any N ∈ N * , and i = 1, . . . , N , D i (µ, Φ) = w(µ), φ i .
We take a partition {P 1 , . . . , P K } of the parameter space P, that is:

P = ∪ K k=1 P k and k = k ′ ⇒ P k ∩ P k ′ = ∅.
We set, for i = 1, . . . , N and k = 1, . . . , K:

β min i,k (Φ) = min µ∈P k D i (µ, Φ), β max i,k (Φ) = max µ∈P k D i (µ, Φ),
and:

β up i (µ, Φ) = β max i,k(µ) (Φ) if r(µ), φ i > 0 β min i,k(µ) (Φ) else,
β low i (µ, Φ) = β min i,k(µ) (Φ) if r(µ), φ i > 0 β max i,k(µ) (Φ) else,
where k(µ) is the only k in {1, . . . , K} so that µ ∈ P k . We also set:

T low 1 (µ, N, Φ) = N i=1 r(µ), φ i β low i (µ, Φ), T up 1 (µ, N, Φ) = N i=1 r(µ), φ i β up i (µ, Φ), T 1 (µ, N, Φ) = max T low 1 (µ, N, Φ) , |T up 1 (µ, N, Φ)| .
Finally, we suppose that µ is a random variable on P and set:

T 2 (N, Φ) = E µ N i=N +1 w(µ), φ i r(µ), φ i .
We have the following theorem:

Theorem 1.1. For any α ∈]0; 1[ and for any N ∈ N * , we have:

P |s(µ) -s(µ)| > T 1 (µ, N, Φ) + T 2 (N, Φ) α ≤ α.
Proof: We begin by noticing that:

A(µ) -1 r(µ) = Z u(µ) -u(µ) so that: s(µ) -s(µ) = l, Z u(µ) -u(µ) = l, A(µ) -1 r(µ) = w(µ), r(µ) .
We expand the residual in the Φ basis:

r(µ) = i≥1 r(µ), φ i φ i .
Hence:

s(µ) -s(µ) = i≥1 l, A(µ) -1 φ i r(µ), φ i = i≥1 w(µ), φ i r(µ), φ i . (9)
We clearly have that for any N ∈ N * :

N i=1 r(µ), φ i β low i (µ, Φ) ≤ N i=1 r(µ), φ i w(µ), φ i ≤ N i=1 r(µ), φ i β up i (µ, Φ)
and this implies:

N i=1 r(µ), φ i w(µ), φ i ≤ T 1 (µ, N, Φ). (10)
So we have:

P |s(µ) -s(µ)| > T 1 (µ, N, Φ) + T 2 (N, Φ) α ≤ P |s(µ) -s(µ)| > N i=1 r(µ), φ i w(µ), φ i + T 2 (N, Φ) α by (10) = P |s(µ) -s(µ)| - N i=1 r(µ), φ i w(µ), φ i > T 2 (N, Φ) α ≤ P N i=N +1 r(µ), φ i w(µ), φ i > T 2 (N, Φ) α by (9)
≤ α thanks to Markov's inequality.

Choice of Φ The error bound given in Theorem 1.1 above is valid for any orthonormal basis Φ. For efficiency reasons, we would like to choose Φ so that the parameter-independent part T 2 (N, Φ) is the smallest possible, for a fixed truncation index N ∈ N * . To our knowledge, minimizing T 2 (N, Φ) over orthonormal bases of X is an optimization problem for which no efficient algorithm exists. However, we can minimize an upper bound of T 2 (N, Φ).

We define an auto-adjoint, positive operator G : X → X by:

∀φ ∈ X, Gφ = 1 2 E µ ( r(µ), φ r(µ) + w(µ), φ w(µ)) . (11)
Let

λ 1 ≥ λ 2 ≥ . . . λ N ≥ 0 be the eigenvalues of G. Let, for i ∈ {1, 2 . . . , N }, φ G
i be an unit eigenvector of G associated with the i th eigenvalue, and

Φ G = {φ G 1 , . . . , φ G N }. We can state that: Theorem 1.2. T 2 (N, Φ G ) ≤ N i=N +1 λ 2 i .
Proof. We have:

T 2 (N, Φ) ≤ 1 2 E µ N i=N +1 w(µ), φ i 2 + N i=N +1 r(µ), φ i 2 =: T sup 2 (N, Φ) = N i=N +1 Gφ i , φ i
Using Theorem 1.1 of [START_REF] Volkwein | Proper orthogonal decomposition and singular value decomposition. Spezialforschungsbereich F003 Optimierung und Kontrolle, Projektbereich Kontinuierliche Optimierung und Kontrolle[END_REF], we get that the minimum of T sup 2 (N, Φ) is attained for Φ = Φ G , and that minimum is

N i=N +1 λ 2 i .
This theorem suggests to use Φ = Φ G , so as to control T 2 (N, Φ).

Monte-Carlo approximation of the error bound

In this Subsection, we present an implementable offline/online procedure for the estimation of the upper bound for | s(µ)s(µ)| presented in Theorem 1.1.

Estimation of φ G i

We fix a truncation index N ∈ N * , and we estimate {φ G i } i=1,...,N by using a modification of the method of snapshots used in Proper Orthogonal Decomposition [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part i-ii[END_REF]. This estimation is performed during the offline phase. We begin by estimating the G operator by G, then we approximate φ G i by the appropriate eigenvectors of G.

Estimation of G:

We take a finite (large), subset of parameters Ξ ⊂ P, randomly sampled from the distribution of the parameter, and we approximate the G operator by:

Gφ = 1 2#Ξ µ∈Ξ ( r(µ), φ r(µ) + w(µ), φ w(µ))
In other words, G is a Monte-Carlo estimator of G. We take { φ G i } i=1,...,N as the unit eigenvectors associated with the N largest eigenvalues of G.

Computation of the eigenvalues of G:

The operator G admits the following matrix representation:

G = 1 2#Ξ W W t + RR t ,
where W (resp. R) is the matrix whose columns are the components of w(µ) (resp. r(µ)) in a basis of X, for µ ∈ Ξ. These two matrices have #Ξ columns and dim X lines, which means that the matrix above is dim X × dim X.

In general, we take #Ξ ≪ dim X, and so it is computationally advantageous to notice that if φ is an eigenvector of G associated with a nonzero eigenvalue λ, then:

1 λ 1 2#Ξ (W W t φ + RR t φ) = φ, so that φ ∈ Im W + Im R =: V.
Hence, if V is the matrix of an orthonormal basis of V, then there exists ψ so that φ = V ψ and we have:

W W t φ + RR t φ = λφ =⇒ V t 1 2#Ξ (W W t + RR t ) V ψ = λψ.
As a consequence, it is sufficient to find the dominant eigenvectors

ψ G 1 , . . . , ψ G N of the matrix Σ = 1 2#Ξ V t (W W t + RR t )V (of size 2#Ξ), and to deduce φ G i from ψ G i by the relation φ G i = V ψ G i .
Besides, by writing Σ as:

Σ = 1 2#Ξ (V t W )(W t V ) + (V t R)(R t V ) ,
it is possible to compute and store Σ without storing nor computing any dense dim X × dim X matrix.

Computation of T 1 (µ, N, Φ) All the quantities intervening in T 1 (µ, N, Φ) can be straightforwardly deduced from β min,max i,k and r(µ), φ G i . Computation of β min i,k (Φ) and β max i,k (Φ): For i = 1, . . . , N and k = 1, . . . , K, the reals β min i,k (Φ) and β max i,k (Φ) can be computed during the offline phase, as they are parameter-independent. If the Θ q s functions are smooth enough, a quasi-Newton optimization such as L-BFGS [START_REF] Zhu | L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization[END_REF] can be used so as to compute these reals; they can also be approximated by a simple discrete minimization:

β min i,k (Φ) = min µ∈Ξ∩P k D i (µ, Φ), β max i,k (Φ) = max µ∈Ξ∩P k D i (µ, Φ). (12)
One should note here that one has to be careful, as these numerical optimisation procedures are not exact. Indeed, the discrete optimisation is only an approximation, and the Quasi-Newton optimisation may be handled carefully, so as to avoid being trapped in local extrema.

Computation of r(µ), φ G i :

We denote by {ζ 1 , . . . , ζ n } are the column vectors of Z, which form a basis of the reduced space X. We can, during the offline phase, compute the following parameter-independent quantities:

f q ′ , φ G i , A q ζ j , φ G i (i = 1, . . . , N, j = 1, . . . , n, q = 1, . . . , Q, q ′ = 1, . . . , Q ′ ).
Let a parameter µ ∈ P be given, and u 1 (µ), . . . , u n (µ) be the components of the reduced solution u(µ) in the reduced basis {ζ 1 , . . . , ζ n }. By using the relation:

r(µ), φ G i = Q q=1 Θ q (µ) n j=1 u j (µ) A q ζ j , φ G i - Q ′ q ′ =1 γ q ′ (µ) f q ′ , φ G i ,
the dot products between the residual and φ G i can be computed in the online phase, with a complexity of O(nQ+Q ′ ) arithmetic operations, O(Q) evaluations of Θ functions and O(Q ′ ) evaluations of γ functions, which is independent of dim X.

Estimation of T 2 (N, Φ) We approximate T 2 (N, Φ) by computing the following Monte-Carlo estimator:

T 2 (N, Φ) = 1 2#Ξ µ∈Ξ s(µ) -s(µ) - N i=1 w(µ), φ i r(µ), φ i .
As this quantity is µ-independent, it can be computed once and for all during the offline phase.

Final error bound By using Theorem 1.1, we get that for ǫ(µ, α, N, Φ) = T 1 (µ, N, Φ) + T 2 (N, Φ)/α, we have:

P (|s(µ) -s(µ)| ≥ ǫ(µ, α, N, Φ)) ≤ α,
so we may take, as estimated (computable) error bound with risk α,

ǫ(µ, α, N, Φ) = T 1 (µ, N, Φ) + T 2 (N, Φ) α . (13)
In the rest of the text, this computable error bound is designated as the error bound on the non-corrected output, by contrast to the bound described in Section 1.2. One may note that estimating T 2 (N, Φ) by T 2 (N, Φ) causes some error on the risk of the computable error bound ǫ(µ, α, N, Φ). This error is analyzed in Appendix A.

Bound on the corrected output

The reduced output s(µ) is a natural reduced output that approximates s(µ). It is possible to solve an auxiliary problem in order to compute an error correction that improves the order of convergence of the reduced output. As we will see, the error bound presented above can be easily modified so as to certify the corrected output.

Output correction

The idea of solving an adjoint problem in order to improve the order of convergence of a scalar output has originated in [START_REF] Machiels | A general framework for finite element a posteriori error control: application to linear and nonlinear convection-dominated problems[END_REF], with the first application in the reduced basis context in [START_REF] Maday | A blackbox reduced-basis output bound method for noncoercive linear problems[END_REF]. We introduced the so-called adjoint problem, whose solution w(µ) satisfies:

A(µ) t u d w(µ) = l,
and the solution ũd w(µ) of the reduced adjoint problem:

Z t d A(µ) t Z d ũd w(µ) = Z t d l,
where Z d is the selected matrix of the reduced basis for the adjoint basis.

The corrected output is:

sc (µ) = s(µ) -Z ũd w(µ), r(µ)
One shoud note that the corrected output can be computed using an efficient offline-online procedure which requires two reduced basis solutions, hence roughly doubling (when Z d has the same number of columns than Z) the offline and online computation times, except in the particular case where A is symmetric and l is proportional to f .

Existing error bound on the corrected output In [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF], it is shown that:

|s(µ) -sc (µ)| ≤ r(µ) ′ r d (µ) ′ α(µ) =: ǫ cc (µ) (14)
where α(µ) is the stability constant bound defined at (7) (in the symmetric case) and r d (µ) is the dual residual:

r d (µ) = A t (µ)Z d ũd w(µ) -l(µ).
Hereafter, the existing bound ǫ cc is called dual-based error bound.

Proposed error bound It is clear that the work performed in the above sections can be reused so as to provide a probabilistic bound on |s(µ)sc (µ)|, by simply replacing w(µ) by:

w c (µ) = w(µ) -Z ũd w(µ), (15)
and hence giving a competitor (called error bound on the corrected output) to the dual-based error bound.

Summary of the different bounds

To sum up, we have presented four output computable error bounds. Two of them are bounds for the uncorrected ouput error:

• the "Lipschitz" error bound ǫ L (8);

• the estimated error bound on the uncorrected output ǫ, that we propose in [START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF];

and two of them for the corrected output error:

• the existing dual-based error bound ǫ cc , defined at (14);

• the estimated error bound on the corrected output ǫ c , that is (13) amended with [START_REF] Rozza | Potential flow[END_REF].

Comparison of the various presented bounds. We can now discuss and compare the different considered error bounds (numerical comparisons are gathered in sections 3 and 4). Online costs: using ǫ cc or ǫ c will require a double computation time (except if we are in a "compliant" case, i.e. l = f ), when compared with using ǫ L or ǫ. However, this computational supplement enables error correction of reduced output and an improved error estimation. Except for the possible supplemental cost for the adjoint problem resolution, the online cost of ǫ c or ǫ should be comparable to the cost of ǫ cc , if not slightly better for the former, as no optimization problem has to be solved in the online phase (while ǫ cc , when using the Successive Constraints Method (SCM) [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF], requires a linear programming in order to compute an online α(µ)).

Offline costs: the SCM procedure used for ǫ cc and ǫ L requires the offline resolution of a number of large eigenproblems on the "finite element space" X.

Depending on the problem at hand, this may, or may not be more expensive than the optimization problems required to compute the β constants (see [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations[END_REF]) and the Monte-Carlo estimation of T 2 . One of the methods can also be feasible and the other not. Probabilistic vs. deterministic: one can also notice that the ǫ c and ǫ bounds are probabilistic in nature; they increase when the risk level decreases. In practice, as shown in the numerical experiments of Section 4, our probabilistic bound is much shaper, even while choosing a very small risk level. Accuracy: the classical bound introduces the dual norm of the residual, which causes a loss of accuracy. In this work, we avoid this step by using the majoration in [START_REF] Machiels | A general framework for finite element a posteriori error control: application to linear and nonlinear convection-dominated problems[END_REF] and a probabilistic argument.

To conclude this discussion, one can say that there may not be a definitive winner error bound, and that the best choice highly depends on the problem at hand (dimension of the X space, numbers of terms in the affine decompositions), the computational budget, the required precision, the number of online queries to perform, and the probability of failure of the error bound that one can afford.

Application to sensitivity analysis

Our error estimation method is applied in sensitivity analysis, so as to quantify the error caused by the replacement of the original model output by the reduced basis output during the Monte-Carlo estimation of the Sobol indices. For the sake of self-completeness, we briefly present the aim and the computation of these indices, and we refer to [START_REF] Saltelli | Wiley Series in Probability and Statistics[END_REF], [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF] and [START_REF] Janon | Uncertainties assessment in global sensitivity indices estimation from metamodels[END_REF] for details.

Definition of the Sobol indices

For i = 1, . . . , p, the i th Sobol index of a function of p variables s(µ 1 , . . . , µ p ) is defined by:

S i = Var (E(s(µ 1 , . . . , µ p )|µ i )) Var (s(µ 1 , . . . , µ p )) , (16)
the variances and conditional expectation being taken with respect to a postulated distribution of the (µ 1 , . . . , µ p ) input vector accounting for the uncertainty on the inputs' value. These indices are well defined as soon as s ∈ L 2 (P) and Var (s(µ 1 , . . . , µ p )) = 0. When µ 1 , . . . , µ p (stochastically) independent, the i th Sobol index can be interpreted as the fraction of the variance of the output that is caused by the uncertainty on the i th parameter µ i . All the Sobol indices lie in [0; 1]; the closer to zero (resp., one) S i is, the less (resp., the more) importance µ i 's uncertainty has on s's uncertainty.

Estimation of the Sobol indices

The conditional expectation and variances appearing in ( 16) are generally not amenable to analytic computations. In those cases, one can estimate S i by using a Monte-Carlo estimate: from two random, independent samples of size M of the inputs' distribution, we compute 2M appropriate evaluations {s j } and {s ′ j } of s, and estimate S i by:

S i = 1 M M j=1 s j s ′ j -1 M M j=1 s j 1 M M j=1 s ′ j 1 M M j=1 s 2 j -1 M M j=1 s j 2 . ( 17
)
When M and/or the required time for the evaluation of the model output are large, it is computationally advantageous to replace s by its surrogate model s. By using [START_REF] Saltelli | Wiley Series in Probability and Statistics[END_REF] on s (hence with reduced model outputs { s j } and { s ′ j }), one estimates the Sobol indices of the surrogate model rather than those of the true model. We presented in [START_REF] Janon | Uncertainties assessment in global sensitivity indices estimation from metamodels[END_REF], Sections 3.1 and 3.2, a method to quantify the error made in the Sobol index estimation when replacing the original model by the surrogate one. We defined two estimators S m i,αas/2 and S M i,1-αas/2 , relying on output error bound samples {ǫ j } and {ǫ ′ j }, and proved that: Theorem 2.1. If:

∀j = 1, . . . , M, |s j -s j | ≤ ǫ j and s ′ j -s ′ j ≤ ǫ ′ j ,
then we have:

P S i ∈ [ S m i,αas/2 ; S M i,1-αas/2 ] ≥ 1 -α as .
In our case, the output error bound ǫ(µ) of Theorem 1.1 does not satisfy the above hypothesis, but satisfies a weaker "probabilistic" statement. This is the object of the following Corollary:

Corollary 2.2. If: ∀j = 1, . . . , M, P (|s j -s j | ≥ ǫ j ) ≤ α and ∀j = 1, . . . , M, P s ′ j -s ′ j ≥ ǫ ′ j ≤ α,
then we have:

P S i ∈ [ S m i,αas/2 ; S M i,1-αas/2 ] ≥ (1 -α as ) × (1 -α) 2M .
Proof. We easily have that:

P S i ∈ [ S m i,αas/2 ; S M i,1-αas/2 ] ≥ P S i ∈ [ S m i,αas/2 ; S M i,1-αas/2 ] | ∀j, |s j -s j | < ǫ(µ) ×P (∀j, |s j -s j | < ǫ(µ)) ≥ (1 -α as ) × (1 -α) 2M .

Numerical results I: Diffusion

Important instances of problem (1) appear as discretizations of µ-parametrized linear partial differential equations (PDE); the X space is typically a finite element subspace (e.g., Lagrange P 1 finite elements), which we still see as identical to R N . A(µ) and f are given by Galerkin projection of the weak form of the PDE onto a suitable basis of X. The boundary conditions of the PDE are usually either encoded in X or in A(µ), and the inner product used to perform the Galerkin projection is typically the L 2 or H 1 inner product. The use of the standard Euclidean product is justified by the fact that the relevant functional inner product has already been used to write A(µ) and f (µ), and the discrete matricial problem can be considered using the Euclidean inner product.

Benchmark problem

Our benchmark problem [START_REF] Rozza | Potential flow[END_REF] is the following: given a parameter vector

µ = (µ 1 , µ 2 , µ 3 ) ∈ P = [0.25, 0.5] × [2, 4] × [0.1, 0.2],
we consider the domain Ω = Ω(µ) below:

1 µ 1 µ 2 4 µ 3 4 Γ N Γ D Ω Our continuous field variable u e = u e (µ) ∈ X e satisfies:        ∆u e = 0 in Ω u e = 0 on Γ D ∂ue ∂n = -1 on Γ N ∂ue ∂n = 0 on ∂Ω \ (Γ N ∪ Γ D ) (18)
where

X e = {v ∈ H 1 (Ω) s.t. v| ΓD = 0},
∆ denotes the Laplace operator, and ∂ ∂n is the normal derivative with respect to ∂Ω. This continuous variable denotes the potential of a steady, incompressible flow moving in a tube whose profile is given by Ω, with open ends on Γ N and Γ D . The Neumann boundary condition on Γ N states that the fluid enters by Γ N with velocity equal to one, the condition on ∂Ω \ (Γ N ∪ Γ D ) states that the velocity field is tangential to the boundary of the tube; finally the Dirichlet condition on Γ D guarantees well-posedness, as the potential field is determinated up to a constant. The variational formulation of (18) states as follows: find u e = u e (µ) ∈ X e so that:

Ω ∇u e • ∇v = - ΓN v, ∀v ∈ X e .
This variational problem is well-posed, as the bilinear form (u, v) → Ω ∇u • ∇v is coercive on X e (see, for instance, [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF], lemma A.14). The above variational problem is discretized using a finite triangulation T of Ω and the associated P 1 (T ) (see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] or [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]) finite element subspace: find u ∈ X so that

Ω ∇u • ∇v = - ΓN v ∀v ∈ X,
where X = {v ∈ P 1 (T ) s.t. v| ΓD = 0}. In our experiments, dim X = 525. The affine decomposition of the matrix of the bilinear form in the left-hand side of the above equation is obtained by using a piecewise affine mapping from Ω(µ) to a reference domain Ω as explained in [START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF], page 11. Our scalar output of interest is taken to be:

s(µ) = ΓN u(µ),
and µ has uniform distribution on P.

Results

We now present the numerical results obtained using the different error bounds on the output of the model described above. We report our bounds on the noncorrected and corrected outputs, as well as the dual-based output bound. Note that the stability constant α(µ) is taken as the exact inf; this clearly advantages the dual-based output bound. For the comparisons to be fair, one should compare the error bounds of same online cost. It is widely assumed that there exists a constant C so that this cost is C × 2(dim X) 3 for the dual-based method, and C(dim X) 3 for our method, since dual-based method involves online inversion of two linear systems of size dim X, and one system of the same size for our method. Hence, the reported reduced basis sizes for the dual method are multiplied by a factor 3 √ 2. In all cases, the reduced bases are computed using POD with snapshot size 80. To compute G, we use a snapshot of size 200. We also took K = 1 (ie., a trivial partition of P). The truncation index N is taken equal to 20. We used a discrete minimization procedure to estimate the β i constants.

In Figure 1, we compare the different error bounds on the non-corrected, and corrected output. For instance, for the error bound on the non-corrected output, we plot:

ǭ = 1 #S µ∈S ǫ(µ, α, N, Φ)
where S is a random subset of P with size 200 and ǫ(µ, α, N, Φ) is defined at [START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations and applications[END_REF]. Other error bound means are computed accordingly. We also computed the mean of the Lipschitz error bound ǫ L . It is not reported here as it was way higher than dual-based output error bound. We observe that our new output error bound outperforms the dual-based error bound, at least for finite reduced bases sizes. Two reasons can be seen to this superiority: the difference in nature (probabilistic vs. sure) between the two bounds, and the fact that we make a crucial use of expansion (9) instead of using a Cauchy-Schwarz (or duality norm) argument. The rate of convergence (slope of the curve) of the corrected output is better than the non-corrected one, and this fact is reported by the two error bounds. Also, the expectation of T 2 was estimated at 10 -12 , which allows to choose a low target risk and remain very competitive, as the α-dependency of the bound is in T 2 /α.

Application to sensitivity analysis

We estimate confidence intervals for the sensitivity indices of s(µ) by using the method described in [START_REF] Janon | Uncertainties assessment in global sensitivity indices estimation from metamodels[END_REF], together with Corollary 2.2, and the non-corrected output.

We take M = 1000 as sample size, B = 500 as number of bootstrap replications (this parameter is used in the procedure which provides the confidence intervals from B replications of S i , see [START_REF] Janon | Uncertainties assessment in global sensitivity indices estimation from metamodels[END_REF]), dim X = 10 as reduced basis size, α = 0.00001 as output error bound risk, and α as = 0.05 as Monte-Carlo risk. The level of the combined confidence interval S m i,αas/2 ; S M i,1-αas/2 is then (1α as )(1α) M > 0.93. The results are gathered in Table 1. The spread between S m i and S M i accounts for the metamodel-induced error in the estimation of the Sobol indices. The remaining spread between S m i,αas/2 and S M i,1-αas/2 is the impact of the sampling error (due to the replacement of the variances in the definition of the Sobol indices by their empirical estimators). We see that, in this case, the metamodelinduced error (certified by the use of our goal-oriented error bound) is very small with regard to the sampling error. We also notice that the estimate for the Sobol index for µ 3 is negative; this is not contradictory as it is the true value of the index that is in [0, 1]. For small indices, the estimate can be negative.

Numerical results II: transport equation

We now apply our error bound on a non-homogeneous linear transport equation. Compared to the previous example, the considered PDE is of a different kind (hyperbolic rather than elliptic).

Benchmark problem

In this problem, the continuous field u e = u e (x, t) is the solution of the linear transport equation: and boundary condition:

u e (x = 0, t) = 0 ∀t ∈ [0, 1].
The parameter µ is chosen in P = [0.5, 1] and P is endowed with the uniform measure.

We now choose a spatial discretization step ∆x > 0 and a time discretization step ∆t > 0, and we introduce our discrete unknown u = (u n i ) i=0,...,Nx;n=0,...,Nt where N x = 1 ∆x , and N t = 1 ∆t .

We note here that the considered PDE is hyperbolic and time-dependent, and that we perform the reduction on the space-time unknown u, of dimension

Conclusion

We have presented a new explicitly computable output error bound for the reduced-basis method. We have shown, on two different practical examples, that this bound is clearly better than the naive Lipschitz bound and that, at the expense of a slight, controllable risk, the performances of this new bound are better than the ones of the existing dual-based output error bound.

A Estimation of T 2 (N, Φ): error analysis

We now assess the error of estimation of T 2 (N, Φ), and its consequence on the risk of the error bound ǫ(µ, α, N, Φ). The error on the risk of ǫ(µ, α, N, Φ) is not easily attainable, because the same Ξ set of parameters is used to choose Φ and to compute T 2 (N, Φ), leading to probabilistic dependence in the family:

N i=N +1
w(µ), Φ i r(µ), Φ i ; µ ∈ Ξ .

However, one can take another random sample Ξ ′ ⊂ P of parameters, independent of Ξ, with size M = #Ξ ′ , and define T ′ 2 (N, Φ), the following estimator of T 2 (N, Φ):

T ′ 2 (N, Φ) = 1 M µ∈Ξ ′ N i=N +1
w(µ), Φ i r(µ), Φ i ,

  basis size Bound on non-corrected output Bound on corrected output Dual-based output bound True error on non-corrected output True error on corrected output

Figure 1 :

 1 Figure 1: Comparison of the mean error bound on the non-corrected output, the mean dual-based error bound (ǫcc) and the mean error bound on the corrected output (for risk α = 0.0001). The "equivalent" reduced basis sizes are in abscissae.

  ∂u e ∂t (x, t) + µ ∂u e ∂x (x, t) = sin(x) exp(-x)for all (x, t) ∈]0, 1[×]0, 1[, satisfying the initial condition:u e (x, t = 0) = x(1x) ∀x ∈ [0, 1],

  basis size Bound on non-corrected output Bound on corrected output Dual-based output bound True error on non-corrected output True error on corrected output

Figure 2 :

 2 Figure 2: Comparison between the mean error bound (corr. and non-corr. outputs), the mean dual-based error bound, and the true (corr. and non-corr.) errors, for different reduced basis sizes.

  output bound (with Greedy basis) Dual-based output bound (with POD basis)

Figure 3 :

 3 Figure 3: Mean of dual-based error bound, when using a POD or a Greedy basis selection procedure, for different reduced basis sizes.

  First, by repeating the proof of Theorem 1.1, and by replacing (for the application of Markov inequality) T 2 (N, Φ)/α by T 2 (N, Φ)/α, one getsP |s(µ)s(µ)| > T 1 (µ, N, Φ) + T 2 (N, Φ) α ≤ α T 2 (N, Φ) T 2 (N, Φ) .If T 2 (N, Φ) = 0, we have thatP |s(µ)s(µ)| > T 1 (µ, N, Φ) + T 2 (N, Φ) α = 0,hence the computable error bound has zero risk. So we can assume that T 2 (N, Φ) = 0.

Table 1 :

 1 Results of the application of Section 2 to the estimation of the Sobol indices of the output of our benchmark model.

	Input parameter	S m i ; S M i	S m i,αas/2 ; S M i,1-αas/2
	µ 1	[0.530352;0.530933]	[0.48132; 0.5791]
	µ 2	[0.451537;0.452099]	[0.397962;0.51139]
	µ 3	[0.00300247;0.0036825] [-0.0575764;0.0729923]

Acknowledgements. This work has been partially supported by the French National Research Agency (ANR) through COSINUS program (project COSTA-BRAVA nr. ANR-09-COSI-015). We thank Anthony Nouy (École Centrale de Nantes) and Yvon Maday (Université Paris 6) for fruitful discussions, and the two anonymous referees for their pertinent remarks, which have greatly improved the quality of the paper.

(N x + 1) • (N t + 1). This is different from reducing the space-discretized equation at each time step. The u vector satisfies the discretized initial-boundary conditions:

∀n, u n 0 = 0 [START_REF] Van Der | Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF] and the first-order upwind scheme implicit relation:

Let us denote by B = B(µ) (resp. y) the matrix (resp. the vector) so that ( 19),( 20) and ( 21) are equivalent to:

that is:

so that equation ( 23) is [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis[END_REF] with

The output of interest is: s(µ) = u Nt Nx . In the following, we take ∆t = 0.02 and ∆x = 0.05. As in the previous example, the true stability constants are computed for the dual-based error bound.

Results

Comparison of the bounds

We took a very low risk level α = 0.0001, a snapshot size of 200, N = 20 retained φ G i vectors and K = 1. The results (Figure 2) show that, once again, the error bounds we propose in this paper outperforms the dual-based error bound.

Choice of basis

For the comparison to be fair, we have also checked that using a POD basis does not penalize the dual-based error bound (ǫ cc ), by comparing the performance of this error bound when using a POD basis and a so-called "Greedy" [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis[END_REF] procedure. This procedure has a smaller offline cost than the POD, as it requires less resolutions of the reference problem [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis[END_REF]. The results, shown in Figure 3, show that the Greedy procedure yields to inferior performance.

Coverage of the error bound

We have also estimated the actual error risk by doing computations of the error bound (on the corrected output), together with the actual error on a random sample of parameters of size 20000, using 10-sized bases for primal and dual problem. On this random sample, the estimated error bound is always greater than the true error; hence the error risk estimation appears conservative. which gives in turn another error bound:

which is also computable in practice, but requires more solutions of the reference problem 1 during the offline phase. Notice that T ′ 2 (N, Φ) is a random variable with respect to the probability measure used to sample the Ξ ′ set. We denote by P ′ this probability measure. We denote by ρ the risk majorant of the computable error bound:

, which is also a random variable with respect to P ′ . We have the following theorem:

We have:

where -→ P ′ denotes convergence in P ′ -distribution when M → +∞ and N (0, σ) is the centered gaussian distribution with variance σ 2 .

Hence ρ converges (in a probabilistic sense) to α with rate 1/ √ M .

Proof. We have:

Hence, by the central limit theorem [START_REF] Van Der | Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF],

Now we define g(x) = αT 2 (N, Φ)/x. As T 2 (N, Φ) = 0, g is differentiable in T 2 (N, Φ) and one can use the Delta method [START_REF] Van Der | Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF] to write:

which proves the theorem.


Ratiba Derfoul

Sebastien Da Veiga

Christian Gout

Carole Le Guyader

Elodie Tillier

Image processing tools for better incorporation of 4D seismic data into reservoir models

Keywords: History matching, filtering, segmentation, image comparison

Reservoir engineers have to predict the behavior of a hydrocarbon reservoir by building a simulation model which can reproduce as best as possible the data collected on the field. These data fall into two types: static data, which are invariable in time, and dynamic data, which evolve according to fluid motions in the reservoir. In this paper, we focus on the integration of dynamic data related to 4D inverted seismic data. Such seismic data constitute an invaluable source of information on fluid displacement and geology over extensive areas of the reservoir. However, incorporating them in the reservoir model through a matching process is a challenging task. Classical formulations of the objective function, which computes the misfit between observed data and responses computed by the reservoir model, are not adapted to 4D inverted seismic data. For example, a least square based mismatch is not representative of the visual difference between two seismic images. In this paper, we define a new formulation of the objective function based on simplification of seismic data in order to extract relevant information. This simplification involves filtering and segmentation techniques, as well as image comparison methods rooted in image analysis, such as the local modified Haussdorff distance. We investigate the efficiency of such techniques in the context of seismic data, and illustrate their potential on a synthetic history matching reservoir example.

Introduction

History matching is an inverse problem where some parameters are modified in order to build a reliable reservoir model which can reproduce the field observed data. These data fall into two types: static data, which are invariable in time, and dynamic data, which evolve in time according to fluids motion in the reservoir. Traditionally, dynamic data consist of production data collected at wells. Nowadays, an other type of data called time-lapse seismic data can be acquired. These data are generated by using vibrations to capture a two-dimensional picture of rock layers beneath the surface at the same localisation but at different times. These data constitute a highly informative spatial information. The geologists interpret these data to make an estimated picture which can reveal the hydrocarbon motion under the surface without drilling. Incorporating such data in a reservoir model calls for a specific workflow, which consists of a chain of successive modelling steps. First, a geological model is built and populated with petrophysical properties such as porosity and permeability. Then, this model is upscaled to a reservoir grid used to perform a fluid flow simulation by resolving an hyperbolic equation system [START_REF] Aziz | Petroleum reservoir simulation[END_REF]. The result of these equations provides some simulated production data. Seismic attribute variations (velocity and acoustic impedances) can be derived from the computed pressure and saturation variations through a petro-elastic modeling. Typically, it involves linear equations called the Gassmann equations [START_REF] Gosselin | History Matching Using Time-lapse Seismic (HUTS)[END_REF] to model fluid effects and the Hertz-Mindlin equations [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF] to model pressure effects. Outputs of the petro-elastic model are then usually filtered for consistency with the typical seismic bandwidth. The resulting seismic attribute variations are called the simulated seismic attributes. Finally, comparison between simulated and observed data is performed through the evaluation of a so-called objective function. This objective function measures the misfit between available data and data simulated by a given reservoir model. Figure 1 summarizes the workflow used to build a reservoir model. The history matching process is thus an inverse problem where an optimisation method is used to adjust the unknown parameters in order to reduce the value of the objective function. If the formulation of the objective function is adapted, at the end of the process the simulated data produced by the matched model will be close to the available data. Consequently, the definition of the objective function is crucial for this workflow. Currently, the objective function is based on a least square formulation, denoted by OF LS and given by :

OF LS = 1 2 (d -d s ) t C -1 (d -d s )
where d is the reference data vector, d s is the simulated data vector and C is the covariance matrix, which represents the measurement uncertainties. This formulation encompasses both production and seismic data: as a result, all data are treated the same way. This formulation proved to be efficient for production data, in the sense that it is a good characterization of the error between simulated and real data. Moreover, during the history matching process, the objective function can be reduced significantly. On the contrary, this formulation is not suitable for seismic data. Indeed, the nature of inverted seismic data is totally different from production data. First, production data are few and are represented by curves, whereas seismic data are millions of points in a 3D cube with spatial information. Second, while the production data collected on the field are spoiled only by measurement errors, the available seismic attributes, such 2 as P-impedances, have two types of incertitudes. Besides measurement errors observed during the seismic survey, additional uncertainty comes from the inversion process, which transforms a seismic reflection data into a quantitative rock-property describing the reservoir. Therefore, seismic data are more uncertain than production data: it is pointless to try to match exactly each cell of a seimic impedance cube as is done with a least square comparison. Let us illustrate this fact with a simple example. In Figure 2, we display three seismic maps corresponding to a steam injection scenario in a heavy oil reservoir. The map on the left is the reference data to be matched and shows the growth of a steam chamber. The first simulated map displayed in Figure 2, middle, is the one obtained with an initial model before history matching. The last image on the left corresponds to a scenario without steam injection. Visually, the initial model is much closer to the reference data than the model without injection. Indeed, in this case the steam chamber is reproduced by simulation. However, the value of the objective function computed with the least square formulation between the reference data and the data produced by the initial model is higher (about 341) than the misfit calculated between the reference and the same model without injection (about 120). This example clearly illustrates that the formulation with the least square method doest not characterize the error between two seismic inverted data. Recently, Tillier et al. [START_REF] Tillier | Appropriate formulation of the objective function for the history matching of seismic attributes[END_REF] proposed to replace the least square formulation with a new measure based on image analysis tools. First, reference and simulated seismic data are filtered and clustered in order to extract the relevant information from the images. Then, these simplified images are compared with a modified local Hausdorff distance in order to capture spatial differences. Finally, the total error is computed as the euclidean norm of the Hausdorff distances. They show that this procedure yields a misfit measure which is able to discriminate seismic maps very efficiently on several reservoir cases [START_REF] Tillier | Appropriate formulation of the objective function for the history matching of seismic attributes[END_REF], [START_REF] Tillier | Simultaneous inversion of production data and seismic attributes: Application to a SAGD produced field[END_REF]. However, this new formulation consists of very elementary tools. The filtering step calls for a simple median filter, while the classification step is performed with the k-means algorithm, which does not account for spatial correlations in the image.

In this paper, we go one step further by proposing state-of-the-art methods for both the filtering and clustering steps. On the one hand, we investigate the potential of the NL means filter, which was shown to perform very well on several examples [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF]. Its efficiency still has to be studied in the context of seismic data. On the other hand, we focus on segmentation methods based on active contours without edges [START_REF] Chung | Energy Minimization Based Segmentation and Denoising Using a Multilayer Level Set Approach[END_REF] in order to replace the k-means algorithm. Since segmentation methods use spatial information to detect shapes inside an image, they are expected to behave better than simple clustering. In the context of 4D seismic maps, we will generally extract two or three regions of interest. As a result, we extend the modified local Hausdorff distance to deal with more than two levels, in order to compare images with three target areas. The paper is organized as follows. In section 2, we first describe the new objective function and study the impact of the new NL-means filtering and segmentation steps on the procedure. Section 3 is dedicated to a synthetic resevoir case. We illustrate the potential of this new formulation with respect to the least square approach.

2. Objective function adapted to 4D seismic data: simplification and dissimilarities of images

In this paper, according to the meshing of the reservoir, we choose to divide the three dimensional seismic map into several two dimensional seismic slices and consider each slice as an image. Since a seismic map contains a lot of information, we propose to first extract the most relevant information and focus on it. Then, the following step is to use an appropriate metric to measure the error between the relevant information derived from real seismic attributes and the one derived from the numerical outputs. The proposed formulation then consists of the three following steps:

1. Simplify images to extract the relevant information by filtering and segmentation. This process will be applied on both observed and simulated 4D seismic images. 2. Measure a local dissimilarity between the two simplified images.

Define a global dissimilarity to include it into the objective function.

In what follows, we will first study the NL-means filter [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF] and a variational segmentation technique [START_REF] Chan | Active Contours Without Edges[END_REF]. These tools will provide an efficient simplification step. Then, we adapt the local modified Haussdorff distance which was already proposed for step 2 by Tillier et al. [START_REF] Tillier | Appropriate formulation of the objective function for the history matching of seismic attributes[END_REF]. In order to account for more than two regions of interest, we use the distance transform for real values proposed by Molchanov and Terán [START_REF] Molchanov | Distance transforms for real-valued functions[END_REF]. Finally, step 3 will only consist in computing the L2 norm of the local Haussdorff distances and will not be detailed here.

Simplification of seismic attributes with filtering and image segmentation

Preliminary step: Smoothing images.

As mentioned in the introduction, inverted seismic attributes are noisy and are also tainted by uncertainty due to the inversion process. Moreover, simulated seismic maps can also exhibit small patterns which are somewhat irrelevant and difficult to interpret. In order to prepare the segmentation step, we thus propose to clean up the data by filtering them. The aim of filtering methods is to find the best estimation of an original image from its noisy version. A filtering method is based on the following decomposition:

u = D h (u) + n(D h , u) (1)
where u is the noisy image, D h (u) is the estimated image, n(D h , u) is the estimated noise affecting the data and h is a filtering parameter. Among the several available local methods, some are based on local smoothing filters (arithmetic average filter, Gaussian filter, median filter or filter based on calculus of variations [START_REF] Rudin | Total variation based image restoration with free local constraints[END_REF]). For example, the arithmetic average filter replaces the value of one pixel by an arithmetic average of its neighboring cells or patch, the size of the patch being the parameter of the filter. Although computationally efficient, this filter is known to have serious difficulties at preserving contours. The use of the median filter, which replaces the value of each pixel by the median average value of its patch, makes it possible to lessen this effect. However, a large amount of noise can not be removed by the median filter. Local filters generally share this drawback.

Other filters worth mentioning call for thresholding techniques (Wiener-Fourier filter, Yaroslavsky local adaptive DCT-based filters [START_REF] Yaroslavsky | Digital Picture Processing. An Introduction[END_REF] or wavelet thresholding method [START_REF] Donoho | Ideal spatial adaptation via wavelet shrinkage[END_REF]). Recently, Buades et al. [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF] introduced an efficient smoothing method known as the Non Local means filter (NL means filter). The main idea, first used by Efros and Leung [START_REF] Efros | Texture synthesis by non parametric sampling[END_REF] for texture synthesis, is to take advantage of self-similarities in images by comparing local neighbourhoods across the whole image.

In the case of 4D inverted seismic images, this approach seems to be promising: seismic images are not structured such as natural pictures and consist of geological patterns or fluid displacement over large areas of the image. Coming back to equation ( 1), the NL means estimate of the denoised image is given by

D h (u 0 ) = y∈I exp{-||U(x)-U(y)|| 2h }u(y) y∈I exp{-||U(x)-U(y)|| 2h } (2)
where U(x) and U(y) are the neighborhood windows, or patches, with respectively x and y as the central pixel. The distance between two patches is:

||U(x) -U(y)|| = 1 s | j|≤ s-1 2 |u(x + i) -u(y + i)| ( 3
)
where s is the size of patches. In practice, the filtering parameter h depends on the noise variance of the image. This variance is unknown for a real picture u, but a possible estimation is given by the empirical variance

σ 2 (u) = 1 | I | x∈I (u(x) -u) 2
where I is the set of all pixels of u and u = 1

|I| x∈I u(i).
In order to investigate the behavior of the NL means filter on seismic images, we propose to compare it with the arithmetic average and the median filter, which were already studied in this setting by Tillier et al. [START_REF] Tillier | Appropriate formulation of the objective function for the history matching of seismic attributes[END_REF]. In both examples we set the filtering parameter h to 0.03 and the size of the patch to 7.

The first illustrative example deals with the steam chamber reservoir model presented in the introduction. The noisy image is depicted in Figure 4, top left. Results obtained with the three filters are given in the top row. As expected, the moving average filter exhibits oversmoothing, hence producing a blurred image, while the median filter makes it possible to preserve the steam chamber contours. The NL means filter (Figure 4 -line 1 column 3) also keeps intact the contours. In addition, we provide in Figure 4, bottom row, the noise removed by each filter, i.e. the difference between the noisy image and its denoised version. The image structure can be observed in the noise removed by the arithmetic average filter, this explains why a lot of information is lost in this case. On the contrary, the median filter removes a very small amount of noise: structural information is preserved, but small irrelevant patterns are not cut off. Between these two unwanted behaviors, the NL means filter produces a very clean image, where almost all the noise has been removed while preserving the steam chamber.

The second example focuses on a reservoir model encompassing several channels, i.e. geological structures which are preferential ways for fluid-flow displacement. Here, objects are numerous and contours have to be identified correctly since connectivity between channels is crucial for fluid displacements. On the top left part of Figure 4, the reference image consists of a noisy version of a seismic image. The initial map, containing normalized values between 0 and 1, has been corrupted with a Gaussian white noise with standard deviation equal to 0.03. Red areas correspond to channels with gas displacement and blue areas to channels with water displacement. Once again, the arithmetic average filter over-smooths the channel objects, and a lot of structural and connexity information is lost, as can be seen in the removed noise. Channels are less blurred with the median filter, but many details are removed in the red diagonal part of the image. Finally, the noise removed by the NL means filter looks like a white noise and the channel structure is very well preserved.

These reservoir examples clearly illustrate the large improvement provided by the NL means filter in the seismic data framework. This is the reason why we propose to include it into our objective function formulation. From a computational perspective, let us mention that Buades et al. advise to compute the weighted sum in equation ( 2) only on a subset of the image called the search window, in order to reduce CPU time. However, since relevant information (geological structures, fluid displacements) are observed on extensive areas of the image, we suggest to make this computation on the full image.

Segmentation

Seismic maps usually contain a lot of information. As mentioned before, in order to match the simulated map with the real data, we want to simplify both of them to only keep the main features. Consequently, the number of data will be significantly reduced to compare only the relevant information. This step has to be done automatically without the operation of the user. Indeed, the history mathcing process usually involves hundreds of calls to the objective function during the iterations of the optimization algorithm. Previously in [START_REF] Tillier | Appropriate formulation of the objective function for the history matching of seismic attributes[END_REF], statistical clustering methods such as the k-means algorithm have been tested for the simplification of the data. However, the proposed methods only deal with the absolute value of seismic attributes and do not incorporate information on structures and objects. Here, we propose to improve this clustering step by applying a segmentation based on variational methods, which accounts for spatial correlations in the image. In the image processing community, segmentation is the decomposition of an image into different areas of interest. For example, watersheds transform [START_REF] Beucher | The watershed transformation applied to image segmentation[END_REF] is a basic tool of mathematical morphology for image segmentation. It is based on an analogy comparing images to topographic relief, where contours are defined as the water meeting points in a scenario of flood of the relief. This method uses the gradient of the image and can not be relevant when the edges in images are not well defined. Deformable models consist of algorithms which generate curves that move within images to find object boundaries by minimizing an energy function. Among the existing methods, active geodesic contours (Caselles et al. [START_REF] Caselles | Geodesic Active Contours[END_REF], [START_REF] Caselles | A geometric model for active contours in image processing[END_REF]) and level set methods developed by Osher and Sethian [START_REF] Osher | Fronts propagation with curvature dependent speed : Algorithms based on Hamilton-Jacobi formulations[END_REF] are the reference ones. Generally, they use an edge detector based on the gradient of the image. In our case, even after applying a denoising method, the edges of the objects contained in the seismic map are not accurate. As a result, gradient-based methods are not well suited for seismic maps. The active contours without edges proposed by Vese et al. ( [START_REF] Chan | Active Contours Without Edges[END_REF], [START_REF] Chan | A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model[END_REF], [START_REF] Chung | Energy Minimization Based Segmentation and Denoising Using a Multilayer Level Set Approach[END_REF]) makes it possible to circumvent this limitation. This model uses a stopping term based on Mumford-Shah segmentation technique and has several advantages. First, it detects edges with and without gradient and second, the initial curve does not necessarily have to start around the objects to be detected and instead can be placed anywhere in the image. Then, it gives a partition of the image into two regions, the first region formed by the set of the detected objects and the second one consisting of the background. In our case, the objects to be detected can be for instance channels, steam chambers or any fluid displacement identified in the reservoir. For illustration purposes, we first recall the reduced case introduced by Mumford and Shah [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associed variational problems[END_REF]. Given a function u 0 , which in the context of image segmentation is the image to be segmented, we are looking for its optimal piecewise-constant function approximation. This problem is called the "minimum partition problem". Mumford and Shah [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associed variational problems[END_REF] propose to solve it by minimizing the following functional:

F MS (K, u) = i Ω i | u 0 -c i | 2 +µ K dH 1 ( 4
)
where K = ∪ i ∂Ω i is a set of boundaries, u = c i inside each component Ω i of Ω\K and H 1 is the one-dimensional Haussdorff measure. In two dimensions H 1 is the length of the boundaries K. Generalization to three dimensions involves H 2 , the surface area of K. This energy is minimized with respect to the variable c i by setting it to the mean of the image over Ω i , i.e. c i = Ω i u 0 |Ω i | , where Ω i represents a region where u is constant, and ∂Ω i is the boundary of this region.

In order to minimize this energy over the set of piecewise-constant functions, Zhao et al. [START_REF] Zhao | Variational level set approach to multiphase motion[END_REF] propose to use a variational level-set approach. The unknown boundaries ∂Ω are represented by zero levels of a Lipschitz continuous function. Such a function φ : Ω → R is used to divide the domain Ω in two open regions, defined respectively by {x ∈ Ω, φ(x) < 0} and {x ∈ Ω, φ(x) > 0}. The edges of the objects to be identified in the image are then defined by {x ∈ Ω, φ(x) = 0}. Finally, the energy can be minimized with this following restriction:

{u(x) = c 1 H(φ(x)) + c 2 H(-φ(x))}
where H denotes the Heaviside function. Based on the Mumford-Shah equation, Chan and Vese [START_REF] Chan | Active Contours Without Edges[END_REF] proposed a new model for active contours involving minimization of the following energy:

F(c 1 , c 2 , C) = λ 1 inside(C) (u 0 (x, y) -c 1 ) 2 dx dy + λ 2 outside(C) (u 0 (x, y) -c 2 ) 2 dx dy + µ Length(C) + ν Area(C) ( 5
)
where C is an evolving curve in Ω, µ > 0 and ν > 0 are regularization parameters and λ 1 and λ 2 are the weights assigned to each region.

Let us illustrate this model on the steam chamber example, where the object to be identified is precisely the steam chamber appearing in the middle of the image. Note first that the pixel values have to be positive in order to properly compute c 1 and c 2 , the means inside our target regions (steam chamber and the background). Consequently, original seimsic attributes are linearly mapped to grey level pixels between 0 and 255. In Figure 5, left, we recall the image obtained after NL means filtering. The result obtained after segmentation is given in Figure 5, right. Here, the steam chamber is clearly identified and separated from the non-informative background. Note also that the injection point at the bottom of the image is extracted. In this case, segmentation performs very well and makes it possible to extract the relevant information, i.e. the steam chamber growth. However, in general, 4D seismic images exhibit three types of impedance variations:

• regions with negative variations, which correspond to gas saturation or pressure increase (e.g. steam chamber in previous example);

• regions with positive variations, involving pressure decrease or water saturation increase;

• regions with almost zero variation, where there is no fluid displacement.

Direct application of the previous segmentation methodology will then fail at retrieving these three types of relevant regions. In order to represent more that two regions, a multiphase level-set was introduced by Chan and Vese [START_REF] Chan | A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model[END_REF], where several functions φ i are used. Another generalization was proposed by Vese and Chung [START_REF] Chung | Energy Minimization Based Segmentation and Denoising Using a Multilayer Level Set Approach[END_REF], where they use the idea of multilayer techniques for modeling epitaxial growth from [START_REF] Gyure | Level-set methods for the simulation of epitaxial phenomena[END_REF] in order to represent several regions with fewer levelset functions. Here, we focus on this method and consider only one level-set function with an additional level. This model, which now divides the image into three regions, hinges upon the minimization of the following energy:

F(c 1 , c 2 , c 3 , φ) = Ω | u 0 (x) -c 1 | 2 H(-φ(x)) dx + Ω | u 0 (x) -c 2 | 2 H(φ(x))H(l -φ) dx + Ω | u 0 (x) -c 3 | 2 H(φ(x) -l) dx + µ Ω | ∇H(φ) | + Ω | ∇H(φ(x) -l) | (6)
where c 1 , c 2 and c 3 are the means over each region and µ > 0 is a regularization parameter.

Ω | ∇H(φ) | and Ω | ∇H(φ(x)l) | represent the length of the boundary between R 1 ,R 2 and R 2 , R 3 , respectively. The segmented image is finally given by:

u(x) = c 1 H(-φ(x)) + c 2 H(φ(x))H(l -φ(x)) + c 3 H(φ(x) -l).
In order to write the associate Euler Lagrange equation and to minimize this energy with respect to the level-set function φ, the Heaviside function is approximated with a regularized version H ∈ C 1 . According to [START_REF] Chan | Active Contours Without Edges[END_REF], we retain the following regularization of H:

H (x) = 1 2 1 + 1 π atan x
and its derivative is given by:

δ (x) = H = 1 π 2 + x 2 .
Minimizing the regularized energy with respect to φ, c 1 , c 2 and c 3 yields the associated Euler-Lagrange equations, where the descent direction is parameterized by an artificial time t ≥ 0:

φ(0, x) =φ 0 (x) c 1 (t) = Ω u 0 (x)H (-φ(t, x)) dx Ω H(-φ(t, x)) dx , c 2 (t) = Ω u 0 (x)H (φ(t, x))H(l -φ(t, x)) dx Ω H(φ(t, x))H(l -φ(t, x)) dx , c 3 (t) = Ω u 0 (x)H (φ(t, x) -l) dx Ω H (φ(t, x) -l) dx , ∂φ ∂t =δ (φ) | u 0 -c 1 | 2 -| u 0 -c 2 | 2 H (l -φ) + µ div ∇φ | ∇φ | +δ (φ -l) H (φ) | u 0 -c 2 | 2 -| u 0 -c 3 | 2 +µ div ∇φ | ∇φ | , (δ (φ) + δ (φ -l))∇φ | ∇φ | . n = 0 on δΩ, t > 0. (7)
Here, n is the exterior unit normal to δΩ. This energy is not convex; there is no guaranty to find the optimal minimizer. With this new formulation for segmentation on three regions, we can now investigate its efficiency on the channelized reservoir already presented. In this example, positive variations correspond to gas displacement and negative ones to water movement. As in the steam chamber example, some pre-treatment of the data is required for efficient computation of the means over regions. The linear mapping proposed before is no longer suited, as positive and negative variations do not share the same order of magnitude. Here, a fair balance between positive and negative regions is achieved through histogram equalization [START_REF] Acharya | Image Processing: Principles and Applications[END_REF]: original data are transformed so as to be distributed as close as possible to the uniform distribution on [0, 255]. The seismic image after NL means filtering and histogram equalization is given in Figure 6, left. The corresponding image after three-class segmentation is depicted in Figure 6, right. The regions identified by segmentation perfectly correspond to the expected areas, meaning that all fluid displacements are correctly extracted from the orginal seismic map.

Image comparison and computation of the objective function

In the context of history matching, the previous segmentation step yields two segmented images in two or three regions: one for the reference 4D inverted seismic image and one for the image simulated by the reservoir model. Now, the goal is to establish the differences between these two simplified images and compute an objective function characterizing their mismatch. In a previous work [START_REF] Tillier | Appropriate formulation of the objective function for the history matching of seismic attributes[END_REF], a computation of a Local Dissimilarity Map (LDM) has been efficiently used to compare two binary images. We extend here this formulation to three classes. In this section, we first recall the main ingredients to compute the LDM between two binary images. Then, we extend this method to compute the difference between images encompassing three classes. A global dissimilarities index is finally computed to obtain the value of the objective function between two images.

Local Dissimilarity Map

The Haussdorff distance is known to have an efficient application in image matching [START_REF] Dubuissonand | A Modiffed Hausdorff Distance for Object Matching[END_REF] or face recognition [START_REF] Jesorsky | Robust Face Detection Using the Hausdorff Distance[END_REF]. Given two non-empty sets of points A = (a 1 , ..., a n ) ∈ R 2 , B = (b 1 , ..., b n ) ∈ R 2 and an underlying distance d, the Hausdorff distance between A and B is given by

D H (A, B) = max(h(A, B), h(B, A))
where h(A, B) = max a∈A (min b∈B d(a, b)). In image processing applications, this distance can be viewed as a dissimilarity measure between two binary images A and B. The Hausdorff distance has good properties, but it is quite sensitive to noise since it accounts for the most distant points in A and B. Several modifications of the Hausdorff distance have been proposed to improve it. The interested reader is referred to Zhao [START_REF] Zhao | A new Hausdorff distance for image matching[END_REF] for some overview of these modified distances. However, all of them are global and do not take into account local dissimilarities. To consider such a local aspect, Baudrier et al. [START_REF] Baudrier | The Adaptative Local Hausdorff-Distance Map as a new Dissimilarity Measure[END_REF] introduced the local dissimilarity map (LDM), which compares the two images locally through a sliding window centered on each considered pixel. The authors also propose a criterion to determine, for each pixel, the optimal size of the window. Practically, the LDM has an efficient formulation through the so-called distance transform. The distance transform of a given set associates to every point its distance to the nearest point in this set. More precisely, for binary images, the distance transform at one pixel is the distance between this pixel and the nearest black pixel (if the white color is the background and the black color is the information of the binary image). Any metric could be used to measure the distance between two points, the most classical choice being the Euclidean distance. Using this distance transform, the LDM between two binary images A and B, at a pixel x, is given by

LDM A,B (x) =| A(x) -B(x) | max(dt A (x), dt B (x)) (8)
where dt A (x) and dt B (x) are the distance transform of image A and B at pixel x, respectively. For binary images, this expression rewrites as follows:

LDM A,B (x) = A(x)dt A (x) + B(x)dt B (x).
The definition of the LDM ( 8) can be directly applied to grey level images. However, while the definition of the distance transform is unique for binary images, this is no longer the case when dealing with grey level ones. Extensions of the distance transform for grey images are either achieved through an extended concept of distance between pixels, or with a new background definition. The grey-weighted distance transform and the weighted distance transform on curve space are two examples of such distance transform extending distances between pixels. Both distances are studied by Fouard and Gedda [START_REF] Fouard | An Objective Comparison between Gray Weighted Distance Transforms and Weighted Distance Transforms On Curved Spaces[END_REF]. Besides, an approach based on an extension of the concept of background is proposed by Arlandis [START_REF] Arlandis | Handwritten Character Recognition using the Continuous Distance Transformation[END_REF]. Here, we focus on the distance transform for real-valued functions proposed by Molchanov and Terán [START_REF] Molchanov | Distance transforms for real-valued functions[END_REF]. The main idea is to consider a grey level image I defined on R 2 as a binary volume defined on R 2 × R. Introduction of this additional dimension leads to the definition of the hypograph of image I given by: hypo I = {(x, t) ∈ R 2 , I(x) ≥ t}.

This set consists of binary images obtained with all possible thresholds applied to image I. The original distance transform is then computed for each binary image included in hypo I, and the new distance transform is finally given by their weighted sum over the hypograph. For instance, if I encompasses grey level pixels in [0, 1], the ν-weighted distance transform is defined as follows: for some weights w i . In both formulations above, d is a distance between a point and a set, classically given by d(x, F) = inf{ρ(x, y), y ∈ F} for some point x and a set F. Here, ρ is chosen as the Euclidean distance, but any distance could be used. This new distance transform is then used in equation [START_REF] Caselles | A geometric model for active contours in image processing[END_REF] to compute the LDM between the simulated and the reference segmented images.

Finally, the LDM is summarized in a scalar objective function through L2-norm averaging.

Summary

Finally, the new formulation adapted to 4D seismic data consists of the following steps:

1. Apply the NL means filter on real seismic attributes and the ones derived from the responses of the given reservoir model. 2. Extract relevant information from both maps using • Chan and Vese segmentation method [START_REF] Chan | Active Contours Without Edges[END_REF] for two regions • Vese and Chung segmentation method [START_REF] Chung | Energy Minimization Based Segmentation and Denoising Using a Multilayer Level Set Approach[END_REF] for three regions Note that preliminary transformation of the data may be needed before segmentation (linear mapping, histogram equalization). 3. Compute the LDM between the segmented images with the appropriate distance transform (depending on the number of classes). 4. The final scalar objective function is obtained by computing the L2 norm of the LDM.

Application on an history matching case 4. Conclusion

Figure 1 :

 1 Figure 1: Illustration of the history matching loop.

Figure 2 :

 2 Figure 2: Least square formulation for seismic data. Reference (left), simulated (middle) and non-informative map (right). The value of the objective function between the reference and the simulation is 341, while it is about 120 between the reference and the non-informative map.

Figure 3 :

 3 Figure 3: Comparison of three filters on a steam chamber reservoir case. First row: reference and filtered maps. Second row: noise removal. From left to right: reference, arithmetic average filter, median filter and NL means filter.

Figure 4 :

 4 Figure 4: Comparison of three filters on a channelized reservoir case. First row: reference and filtered maps. Second row: noise removal. From left to right: reference, arithmetic average filter, median filter and NL means filter.

Figure 5 :

 5 Figure 5: Chan and Vese segmentation on the steam chamber reservoir case. Left: Seismic map after NL means filtering. Right: Segmentation result (0: background, 1: detected object) with λ 1 = 1, λ 2 = 4, µ = 4500, ν = 0.

Figure 6 :

 6 Figure 6: Multilayer level set segmentation on the channelized reservoir case. Left: Seismic map after NL means filtering and histogram equalization. Right: Segmentation result (0: first detected region, 1: second detected region, 2: background).

  , I t )ν(dt)where ν is a finite measure on [0, 1] and I t = {x ∈ R 2 , I ≥ t}. Similarly, if I is discretized on [0, 255] and has values between a and b, the ν-weighted distance transform becomes: , I i )w i

Preprint submitted to Journal of Computational and Applied Mathematics December 2, 2011