L Métivier

F Bretaudeau

R Brossier

S Operto

J Virieux

Full Waveform Inversion and the truncated Newton method: quantitative imaging of complex subsurface structures

Keywords: Full Waveform, Theory, Computing aspects, Numerical study, Imaging

Full Waveform Inversion (FWI) is a powerful tool for quantitative seismic imaging from wide-azimuth seismic data. The method is based on the minimization of the misfit between observed and simulated data. This amounts to the resolution of a large-scale nonlinear minimization problem. The inverse Hessian operator plays a crucial role in this reconstruction process. Accounting accurately for the e↵ect of this operator within the minimization scheme should correct for illumination deficits, restore the amplitude of the subsurface parameters, and help to remove artifacts generated by energetic multiple reflections. Conventional preconditioned gradient-based minimization methods only roughly approximate the e↵ect of this operator. We are interested in this study to another class of minimization methods, named as truncated Newton methods. These methods are based on the computation of the model update through a matrix-free conjugate gradient resolution of the Newton linear system. The aim of this study is to present a feasible implementation of this method for the FWI problem, based on a second-order adjoint state 2 L. Métivier, F. Bretaudeau, R. Brossier, S. Operto, J. Virieux formulation for the computation of Hessian-vector products. We compare this method with the nonlinear conjugate gradient and the l-BFGS method within the context of 2D acoustic frequency FWI for the reconstruction of P-wave velocity models. Two test cases are investigated. The first is the synthetic BP 2004 model, representative of the Gulf Of Mexico geology with high velocity contrasts associated with the presence of salt structures. The second is a 2D real data-set from the Valhall oil field in North sea. These tests emphasize the interesting properties of the truncated Newton method regarding conventional optimization methods within the context of FWI.

INTRODUCTION

Full Waveform Inversion (FWI) is a powerful seismic imaging tool, dedicated to quantitative estimations of subsurface parameters such as P-wave and S-wave velocities, density, impedance, or anisotropy parameters. The method is based on the minimization of a misfit function that measures the distance between recorded seismic data and predicted data computed through the numerical simulation of wave propagation. An initial subsurface model is iteratively updated to produce the final estimation.

The formalism of the FWI method has been introduced by [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF] and [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF], based on a time domain discretization of the wave equation. Its first application to 2D synthetic data in the acoustic approximation was performed by [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: numerical results[END_REF].

Later on, a hierarchical frequency domain approach has been introduced by Pratt for crosshole tomography [START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. Part I: acoustic wave-equation method[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part II : elastic wave-equation method[END_REF]). During the past ten years, the simultaneous advances in acquisition systems (development of wide-azimuth seismic surveys for instance) and high performance computing facilities have made possible the successful application of FWI to surface data, both in the 2D acoustic or 2D elastic approximation to reconstruct one or several parameters [START_REF] Operto | Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions: a case study[END_REF][START_REF] Operto | High-resolution crustal seismic imaging from OBS data by full-waveform inversion: application to the eastern-Nankai trough[END_REF][START_REF] Ravaut | Multiscale imaging of complex structures from multi-fold wide-aperture seismic data by frequencydomain full-wavefield inversions: application to a thrust belt[END_REF][START_REF] Gao | Waveform tomography at a groundwater contamination site: surface reflection data[END_REF][START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF][START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[END_REF][START_REF] Plessix | Full waveform inversion and distance separated simultaneous sweeping: a study with a land seismic data set[END_REF]Prieux et al. 2013a,b).

Applications of FWI to real surface data in the 3D acoustic approximation have also been FWI and the truncated Newton method: quantitative imaging of complex subsurface structures 3 performed [START_REF] Sirgue | 3D Frequency Domain Waveform Inversion using Time Domain Finite Di↵erence Methods[END_REF][START_REF] Vigh | 3d full waveform inversion on a gulf of mexico waz data set[END_REF][START_REF] Plessix | Full waveform inversion of a deep water ocean bottom seismometer dataset[END_REF][START_REF] Plessix | Full waveform inversion and distance separated simultaneous sweeping: a study with a land seismic data set[END_REF].

For an overview on the FWI methodology and its applications to synthetic and real case studies, the reader is referred to the survey proposed by [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF].

The simplest optimization methods used in the context of FWI are gradient-based algorithms, such as the steepest descent or the nonlinear conjugate gradient algorithms. From a given initial model, the sequence of updates yielding the final model is defined by the gradient of the misfit function. However, it is well known that these methods have poor convergence properties.

Conversely, Newton-based methods possess better convergence properties (superlinear to quadratic convergence rate). These methods are based on a model update given by the multiplication of the gradient by the inverse Hessian operator ? . The importance of this operator in the context of FWI has been emphasized by [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF]. The inverse Hessian operator acts as a deconvolution operator that accounts for the limited bandwidth of the seismic data and corrects for the loss of amplitude of poorly illuminated subsurface parameters. In addition, it helps to remove artifacts that the second order reflected waves may generate on the model update.

However, because of the large-scale aspect of the FWI problem, which easily involves millions of discrete unknowns in 2D up to billions discrete unknowns in 3D, explicit computation of the inverse Hessian operator is beyond current computational capabilities. As a consequence, research e↵orts have been mainly directed toward direct approximation of this operator.

A first possibility consists in approximating the diagonal of the Hessian. For instance, [START_REF] Operto | Crustal imaging from multi-fold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern Nankai trough[END_REF] compute the diagonal terms of the Gauss-Newton approximation of the Hessian, which requires some extra-computation. A cheaper strategy based on the so-called pseudo-Hessian operator is also proposed by [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF].

A second possibility consists in approximating the inverse Hessian operator using previous values of the gradient of the misfit function. Among this class of methods, known as quasi-Newton methods, the l-BFGS method is quite popular [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. Instead of approximating only the diagonal elements, a positive definite approximation of the full inverse Hessian is computed.

The approximation of the diagonal elements of the Hessian and the l-BFGS strategy can be combined to produce a more accurate approximation of the inverse Hessian operator.

Indeed, the accuracy of the l-BFGS approximation of the inverse Hessian operator can be improved when based on a first estimation of the inverse Hessian operator. This method has been applied in the framework of 2D elastic FWI: [START_REF] Brossier | Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion[END_REF] implemented a l-BFGS optimization using a diagonal pseudo-Hessian as initial guess. This method shows good convergence properties compared to preconditioned non-linear conjugate-gradient. The diagonal estimation can even be updated along the iterations [START_REF] Nocedal | Numerical Optimization[END_REF].

The truncated Newton represents an alternative to these already described optimization methods. At each iteration, the model update is computed as an approximate solution of the Newton equations through a linear iterative solver (namely a conjugate gradient solver) [START_REF] Nash | A survey of truncated Newton methods[END_REF]. Implemented in a "matrix-free" fashion, this iterative solver only requires to compute Hessian-vector products. It is not necessary to form the Hessian operator explicitly.

Although this class of methods is well known in the numerical optimization community, the application of truncated Newton method in the FWI context has still not been fully investigated. Given the importance of the inverse Hessian operator in the FWI reconstruction scheme, we believe that this method could benefit from a better approximation of the inverse Hessian e↵ect, and provides more accurate subsurface parameter estimations than standard optimization schemes. Therefore, the ambition of this study is to focus on the two following points:

• Designing a feasible implementation of the truncated Newton method for FWI in terms of computational time.

• Compare the performance of this method with conventional methods (nonlinear conjugate gradient, l-BFGS) on two realistic test cases.

In Section 2, we describe in more details the principle of the truncated Newton method FWI and the truncated Newton method: quantitative imaging of complex subsurface structures 5 compared to preconditioned gradient-based methods. In Section 3, we present how the method can be implemented e ciently in the FWI context. In Section 4, we present two application cases. The first case is based on the synthetic BP 2004 model, partly inspired from the deep water Gulf of Mexico geology. The presence of salt structures in a marine environment is responsible for high velocity contrasts which make the seismic imaging task di cult. The second test case concerns a 2D line of a real Ocean Bottom Cable (OBC) data-set, acquired in a shallow-water environment at the Valhall oil field, in the North Sea.

This test case is investigated to emphasize how the truncated Newton method behaves when the data is noise-contaminated. Conclusion and perspectives are given in Section 5.

THE TRUNCATED NEWTON SCHEME

For the sake of clarity, the mathematical results presented in this section are formulated for a number of sources equal to 1. The extension to a multi-source context is straightforward, as only explicit summations over the sources are required.

Problem settings

We consider the frequency-domain forward problem

A(m)u = s (1)
where

• m 2 M denotes the subsurface model;

• u 2 W is the complex-valued seismic wavefield;

• s is a source term;

• A(m) is a discretized partial di↵erential operator related to the wave equation (from the acoustic dynamics to the visco-elastic anisotropic dynamics).

The FWI problem is defined as the minimization over the parameter space of a distance between the data predicted by the forward problem and the recorded data.

min m2M f (m) = 1 2 kRu(m) dk 2 , ( 2
)
where

• u(m) is the solution of the forward problem (1) for the source term and the subsurface parameter s and m;

• R is a mapping of the wavefield to the receivers locations;

• d is the data set associated to the source s;

• k.k is a norm in the data space D.

For practical reasons, the use of the L 2 norm is common. However, more general L p norm could be also selected [START_REF] Tarantola | Inverse Problem theory and methods for model parameter estimation[END_REF]. The L 1 norm is, for instance, a good choice when high-amplitude noise (outliers) corrupts the data [START_REF] Brossier | Which data residual norm for robust elastic frequency-domain full waveform inversion?[END_REF]. More complex measurements of the distance between data sets can also be proposed to mitigate the sensitivity of FWI to the initial model. This is, however, beyond the scope of the work presented here.

Preconditioned gradient based-methods

From a numerical point of view, FWI is a large-scale nonlinear minimization problem. The high number of discrete parameters prevents from using global or semi-global optimization techniques to solve this problem. Therefore, we focus on local optimization methods, which are based on the following recurrence: from an initial guess m 0 , a sequence m k is computed such that

m k+1 = m k + k m k , (3)
where m k is the model update and k is a scalar parameter computed through a linesearch or a trust-region procedure [START_REF] Bonnans | Numerical Optimization, Theoretical and Practical Aspects[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF].

Within the framework of Newton algorithms, the increment m k is defined as

H(m k ) m k = rf (m k ). ( 4
)
where H(m) denotes the Hessian operator and rf (m) is the gradient of the misfit function with respect to the model parameter m.

The exact resolution of this linear system at each iteration is beyond current computational capacities in the context of FWI. This is why many large-scale optimization schemes rely on an approximation of the inverse Hessian operator. These schemes are based on the computation of m k as

m k = P k rf (m k ), (5)
where P k is an approximation of the inverse Hessian operator H(m k ) 1 . We shall refer to these schemes as preconditioned gradient-based methods in what follows.

Preconditioned steepest-descent, nonlinear conjugate gradient and l-BFGS methods fall into the scope of these minimization algorithms. These methods only di↵er in the way the matrix P k is computed. The preconditioned steepest-descent uses a "direct" approximation of the inverse Hessian (for instance an approximation of its diagonal [START_REF] Operto | Crustal imaging from multi-fold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern Nankai trough[END_REF][START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF]). The nonlinear conjugate gradient method and the l-BFGS method [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] use former values of the gradient to estimate the inverse Hessian operator † . As mentioned in the introduction, nonlinear conjugate gradient and l-BFGS techniques can be used in combination of a prior estimation of the diagonal to compute the matrix P k [START_REF] Nocedal | Numerical Optimization[END_REF].

Truncated Newton algorithm

Instead of building an approximation P k , the linear system (4) can be solved using a matrixfree version of the conjugate gradient algorithm [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. This only requires the capability of computing Hessian-vector products H(m k )v where v is an arbitrary vector in the model space M. The truncated Newton method thus results in a two nested loops algorithm: † Only the value of the gradient at the previous iteration is used for the nonlinear conjugate gradient, while l previous values are used for the l-BFGS algorithm.

• The external loop consists in the iterative update of the current subsurface parameter estimation, following equation (3).

• The internal loop consists in the iterative resolution of the linear system (4), in order to compute the model update m k . An approximate solution of this system is computed, as described in the following

The advantage of using a truncated Newton method is two-fold:

• The approximation of the inverse Hessian operator is local, while the approximation computed in l-BFGS or nonlinear conjugate gradient method depends on previous iterations.

For minimizing strongly nonlinear misfit functions, this can be advantageous, since the information on the local curvature carried out by previous iterate can be erroneous. Setting the memory parameter l is thus known to be a di cult issue as it is strongly problem dependent [START_REF] Nocedal | Numerical Optimization[END_REF].

• The truncation strategy in the inner loop consists in accounting only for the higher eigenvalues of the inverse Hessian operator. This has an intrinsic regularization e↵ect on the computation of the model update.

In addition, the truncated Newton method o↵ers the possibility of using the approximations of the inverse Hessian operator which have been developed in the context of preconditioned gradient-based methods as preconditioners of the inner linear systems, yielding the resolution of the preconditioned linear system

P k H(m k ) m k = P k rf (m k ), (6)
at each outer iteration. Hence, this method could be seen as a complementary evolution of previously developed minimization strategies.

However, compared to the preconditioned gradient-based methods, the truncated Newton method involves an additional computational expense associated with the iterative resolution of the linear system (4). This should be balanced by an improvement of the convergence speed of the external loop. Therefore, as mentioned by [START_REF] Nash | A survey of truncated Newton methods[END_REF], an e cient implementation of the truncated Newton method relies on the reduction of this additional cost. This can be achieved in the context of FWI by:

• Defining second-order adjoint formulas for the e cient computation of Hessian vector products H(m k )v for any v 2 M.

• Defining an adapted stopping criterion for the approximate resolution of the linear system (4) to limit as much as possible the number of iterations of the conjugate gradient algorithm at each step of the external loop.

• Using an appropriate preconditioner to accelerate the convergence of the resolution of the linear system (4).

The truncated Newton strategy can also be implemented using the Gauss-Newton approximation B(m) of the Hessian operator H(m), which consists in neglecting the secondorder terms of the Hessian operator. The Gauss-Newton approximation is

B(m) = J(m) † R † RJ(m), (7)
where J(m) is the Jacobian matrix:

J(m) = @ m u(m), (8)
and † denotes the transpose conjugate operator.

The Gauss-Newton method is appealing for the following reasons:

• B(m) is positive definite by construction, therefore the conjugate gradient algorithm is well adapted for the resolution of the linear system (4);

• Close to the solution, the residuals are small, therefore the matrix B(m) should be a good approximation of the full Hessian matrix H(m), since the second-order part which is neglected in the Gauss-Newton approximation is proportional to the residuals.

From the implementation point of view, the truncated Gauss-Newton procedure di↵ers only in the computation of matrix-vector products B(m)v instead of H(m)v from the truncated Newton method.

In the next section, we investigate how the truncated Newton and Gauss-Newton methods can be e ciently implemented in the context of FWI.

IMPLEMENTATION OF THE TRUNCATED NEWTON METHOD

Computation of the model update

The resolution of the linear systems (4) first requires to compute the right-hand side rf (m k ).

The computation of rf (m k ) is e ciently achieved through the first-order adjoint-state method, introduced by Lions (1968). For the sake of generality, we will use here the notations introduced by [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF] in his review of the adjoint-state technique for seismic imaging. In this framework, the forward problem should be rewritten as

F (m, u) = 0, (9)
where

F (m, u) = A(m)u s.
(10)

Computation of the gradient and first-order adjoint method

In the context of FWI, the first-order adjoint method amounts to compute rf (m k ) as the zero-lag cross-correlation of the incident wavefield and the adjoint wavefield [START_REF] Chavent | Identification of parameter distributed systems, in Identification of function parameters in partial di↵erential equations[END_REF].

The adjoint wavefield is defined as the solution of

@F (m, u) † @u = R † (d Ru(m)), (11)
where u(m) is the solution of (1). Based on the definition of the adjoint state , the gradient can be expressed as

rf (m) = Re ✓ @F (m, u) † @m (m) ◆ , ( 12
)
where Re denotes the real part operator.

Using this method, the computation cost of the gradient amounts to the resolution of two wave propagation problems per shot: one forward problem (1) and one adjoint problem

Computation of H(m)v through second-order adjoint method

The computation of Hessian-vector products through second-order adjoint methods is a topic that has already been investigated in the field of data assimilation and weather forecasting [START_REF] Wang | The second order adjoint analysis: Theory and applications[END_REF]). However, the control variable in data assimilation is an initial condition for the system, whereas in seismic imaging, the control variable is a coe cient of the partial di↵erential equation that describes the system. A formula for the computation of Hessian-vector products have been given by [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF] in the seismic imaging context, for the Gauss-Newton approximation in the discrete frequency domain. [START_REF] Fichtner | Hessian kernels of seismic data functionals based upon adjoint techniques[END_REF] also propose more general formulas for the computation of Hessian kernels. [START_REF] Epanomeritakis | A Newton-CG method for large-scale three-dimensional elastic full waveform seismic inversion[END_REF] give the formulas corresponding to the elastic case.

We propose here a general framework in the frequency domain for deriving these formulas, with no assumption on the discretization and the kind of partial di↵erential equations that are used for the wave propagation description. The method can be straightforwardly adapted to the time-domain formulation by adding proper initial and final conditions, and boundary conditions. In addition, no prior assumption on the linearity of the forward problem is required, as it was the case in [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF] and a previous work [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF].

We first define the functional h v (m) as

h v (m) = (rf (m), v) M , (13)
where (., .) M denotes the scalar product on the parameter space M. By definition, the

gradient of the functional h v (m) is rh v (m) = H(m)v. ( 14
)
We use the Lagrangian formalism to compute rh v (m). We introduce the Lagrangian oper-ator L(m, u, , g, µ 1 , µ 2 , µ 3 ) such that

L v (m, u, , g, µ 1 , µ 2 , µ 3 ) W = (g, v) M + Re ⇣ g @F (m,u) † @m , µ 1 ⌘ W + Re ⇣ @F (m,u) † @u R † (d Ru) , µ 2 ⌘ W + Re (F (m, u), µ 3 ) W . (15)
In this expression, (., .) W denotes the scalar product in the wavefield space W. In addition, u and play the role of the incident and adjoint wavefields respectively, while g plays the role of the gradient. Let u(m), (m), g(m) satisfying the constraints

F (m, u) = 0 @F (m, u) † @u = R † (d Ru) g = @F (m, u) † @m . ( 16
)
We have

L v (m, u, , g, µ 1 , µ 2 , µ 3 ) = h v (m), (17)
and

@L v @m (m, u, , g, µ 1 , µ 2 , µ 3 ) = rh v (m). ( 18
)
In addition

@L v @m (m, u, , g, µ 1 , µ 2 , µ 3 ) = Re ✓ @ 2 F (m, u) † @m 2 ◆ † µ 3 ! + Re ✓ @ 2 F (m, u) † @m@u ◆ † µ 2 ! + Re ✓ @F (m, u) † @m µ 1 ◆ + @L v @g (m, u, g, µ 1 , µ 2 , µ 3 ) @g(m) @m + @L v @ (m, u, g, µ 1 , µ 2 , µ 3 ) @ (m) @m + @L v @u (m, u, g, µ 1 , µ 2 , µ 3 ) @u(m) @m (19)
We can choose µ i , i = 1, 3 such that for any perturbations dg, d , du we have

8 > > > > > < > > > > > : @L v @g (m, u, g, µ 1 , µ 2 , µ 3 ).dg = 0 @L v @ (m, u, g, µ 1 , µ 2 , µ 3 ).d = 0 @L v @u (m, u, g, µ 1 , µ 2 , µ 3 ).du = 0. ( 20
)
This yields

8 > > > > > > > > > > > > < > > > > > > > > > > > > : µ 1 = v @F (m, u) @u µ 2 = @F (m, u) @ m µ 1 @F (m, u) @u † µ 3 = @ 2 F (m, u) @u 2 † ! † µ 2 @ 2 F (m, u) @u@m † ! † µ 1 R † Rµ 2 . ( 21
)
Using the relation (10) which defines our forward problem, we see that

@F (m, u) @u = A(m), @ 2 F (m, u) @u 2 = 0. ( 22
)
This yields 8 > > > < > > > :

A(m)µ 2 = @F (m, u) @m v A(m) † µ 3 = @ 2 F (m, u) @u@m † ! † µ 1 R † Rµ 2 . ( 23
)
We thus obtain the three-terms Hessian-vector product formula

H(m)v = Re ✓ @ 2 F (m, u) † @m 2 ◆ † µ 1 ! + Re ✓ @ 2 F (m, u) † @m@u ◆ † µ 2 ! + Re ✓ @F (m, u) † @m µ 3 ◆ . ( 24
)
where

@ 2 F (m, u) @m@u (25)
is the radiation matrix.

The adjoint wavefield µ 1 is equal to v. The computation of one Hessian-vector products thus requires to compute the the radiation wavefield u(m) the adjoint wavefield (m), as well as two additional wavefield µ 2 (m) and µ 3 (m). The two latter are computed through the resolution of respectively one forward and one adjoint problem with new source terms.

Gauss-Newton approximation

From the definition of the Gauss-Newton operator B(m) (7), we see that only first-order derivatives of the wavefield u(m) with respect to the model parameter m (also called Jacobian matrix or Fréchet derivatives) are taken into account. The computation of B(m)v can thus be derived from equations ( 23) and ( 24) by neglecting the contribution of all the second-order terms. This yields the following simplifications:

8 > < > : A(m)µ 2 = @F (m, u) @m v A(m) † µ 3 = R † Rµ 2 , (26)
and the following formula for the Gauss-Newton approximation

B(m)v = Re ✓ @F (m, u) † @m µ 3 ◆ . ( 27
)
This formula is consistent with the one derived in [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF] or [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF].

Note however that in these two articles, the method used to derive this formula relies on the assumption of a linear forward problem, which is not the case here.

Computation cost

Consider the resolution of the linear system (4) for the computation of the model update m k . The right hand side in (4) is the opposite gradient rf (m). Using the first-order adjoint state method, it can be computed at the expense of one forward problem for u(m)

and one adjoint problem for (m).

Provided these wavefields can be stored, the computation of H(m)v of B(m)v only requires to solve two additional problems: one forward problem for the computation of µ 2 (m), and one adjoint problem for the computation of µ 3 (m).

Note that in practice, the computation of B(m)v does not make use of the adjoint wavefield (m). However, its computation is imposed by the computation of the gradient, which is the right-hand side of the linear system (4). The computation cost of the action of the Hessian operator or its Gauss-Newton approximation on an arbitrary vector is thus the same in terms of number of wave equations to be solved. Nonetheless, only the incident wavefield u(m) have to be stored in the Gauss-Newton approximation.

The overall computation cost of the truncated Newton method in terms of wave propagation simulation is thus given by

C = N ext (2 + 2 ⇥ N int,k ), ( 28
)
where N ext is the total number of iterations of the external loop and N int,k is the number of conjugate gradient iterations performed at the k th iteration of the external loop. The choice of an appropriate stopping criterion for the conjugate gradient helps to reduce the quantities N int,k . Note that this computational cost depends on the ability of storing u(m) and (m). This is a reasonable assumption for 2D applications. This is a more complex issue in 3D: more sophisticated memory and I/O management methods should be required. In particular, the increase of the number of sources is critical, as one incident and one adjoint wavefield u(m) and (m) have to be stored per sources. The use of source encoding techniques should therefore be crucial in this context.

Definition of an adapted stopping criterion

The Newton method is an iterative minimization of local quadratic expansions of the misfit function. Indeed, the resolution of the system (4) amounts to the minimization of the quadratic form

q k ( m) = f (m k ) + (rf (m k ), m) + 1 2 (H(m k ) m, m). ( 29
)
The definition of the stopping criterion for the truncated Newton method is related to the accuracy of these quadratic expansions. The idea exploited by [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] is the following. Consider a stopping criterion for the CG iterations of the form

kH(m k ) m k + rf (m k )k  ⌘ k krf (m k )k. ( 30
)
where ⌘ k is the forcing term. The role devoted to this forcing term is to account for the accuracy of the local quadratic approximation. When this accuracy increases, ⌘ k should decrease, so as to require to solve the linear system (4) more accurately. Conversely, when the accuracy decreases, ⌘ k should increase, so as to allow a less precise resolution of the linear system (4). This is achieved by defining ⌘ k as the measure of the distance between the first order Taylor expansion of the gradient at the iteration k 1 and the gradient at iteration k:

⌘ k = krf (m k ) rf (m k 1 ) k 1 H(m k 1 ) m k 1 k krf (m k 1 )k . ( 31
)
The definition of the stopping criterion is complemented with an appropriate strategy to deal with the detection of negative eigenvalues of the Hessian operator. The conjugate gradient algorithm is designed for the resolution of symmetric definite positive systems.

However, far from the solution, the full Hessian operator H(m k ) may be indefinite. Therefore, the iterative construction of the solution of (4) in the Krylov space

r 0 , H(m k )r 0 , H(m k ) 2 r 0 , . . . , (32)
where r 0 is the initial residual

H(m k ) m 0 k + rf (m k ), (33)
may use a ascent direction associated with a negative eigenvalue of the operator H(m k ).

In this case, the linear iterations are stopped and the last value of the model update m k which is computed is returned. If this ascent direction is met at the very first linear iteration, the steepest-descent direction is returned. This strategy, proposed by [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF], ensures superlinear convergence properties far from the solution, and quadratic convergence when entering the attraction basin of the minimum.

Preconditioning

In order to speed-up the convergence of the resolution of the linear system (4), it is natural to introduce a preconditioning matrix. In this study we focus on the special preconditioner related to the FWI problem proposed by [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF]. The diagonal elements of the Gauss-FWI and the truncated Newton method: quantitative imaging of complex subsurface structures 17

Newton part of the Hessian B(m) are approximated using the pseudo-Hessian approach, which appears to be relevant for surface seismic survey. Let us denote ↵ j (m) = @ m j u(m) the column j of the Jacobian matrix. Deriving the forward problem with respect to the jth component of the model parameter m j yields

A(m) @u @ m j + A(m) @A @ m j u = 0. ( 34
)
Thus, ↵ j (m) is the solution of the forward problem

A(m)↵ j = @ m j A(m)u. ( 35
)
Using this formula, an exact computation of the entire Jacobian matrix J(m) would thus require to solve m forward problems, which is intractable from a computational cost point of view. Nonetheless, a cheap approximation can be built by approximating the forward problem operator A(m) as the identity matrix I in the left-hand side of equation ( 35). This leads to the definition of the pseudo-Hessian matrix entries:

e H ij (m) = ⇣ [@ m i A(m)u(m)] T ⇥ @ m j A(m)u(m) ⇤ ⌘ , i, j = 1, . . . , M. (36)
The preconditioner used by [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF] is defined by

P k = diag 1 e H ii (m k ) ! i = 1, . . . , M. (37)
However, because of the fast decrease of the wavefield with depth, very small values appear on the diagonal entries of e H(m) corresponding to deep subsurface parameters. Therefore, using directly P k as a preconditioner may yield numerical instabilities. We thus introduce a threshold parameter ✓ 2 R, the constant C k 2 R such that

C k = max j e H jj (m k ), ( 38
)
and we define the matrix P ✓ k such that

P ✓ k = diag 1 e H ii (m k ) + ✓C k ! , i = 1, . . . , M. (39)
Finally the norm of the misfit gradient rf (m) should be preserved by the preconditioner ? .

? Since the stopping criterion for the linear system (4) is based on the reduction of the linear residuals with respect to the norm of the gradient, it appears natural that the preconditioner conserves the gradient norm.

Therefore, we introduce

P ⌫,✓ k = ⌫P ✓ k , (40)
such that

⌫ = krf (m k )k kP ✓ k rf (m k )k . ( 41
)
We use P ⌫,✓ k as a preconditioner of the linear system (4) at each iteration of the external loop. The preconditioner P ⌫,✓ k can also be used as a preconditioner of the nonlinear conjugate gradient method and the l-BFGS method.

CASE STUDIES

4.1 Numerical framework

Forward problem

The numerical tests we present are performed in the 2D frequency domain. In the first case study, we use an isotropic acoustic approximation of the wave propagation. In the second case study, we consider the propagation of acoustic anisotropic waves in a VTI media . In both cases an optimized second-order finite di↵erences scheme with a compact stencil is used [START_REF] Hustedt | Mixed-grid and staggered-grid finite di↵erence methods for frequency domain acoustic wave modelling[END_REF][START_REF] Operto | Finite-di↵erence frequencydomain modeling of visco-acoustic wave propagation in two-dimensional TTI media[END_REF]. Perfectly Matched Layers (PML) [START_REF] Berenger | A perfectly matched layer for absorption of electromagnetic waves[END_REF][START_REF] Métivier | Interlocked optimization and fast gradient algorithm for a seismic inverse problem[END_REF]) are introduced to avoid fictitious reflections on the boundaries of the computation domain, excepted on top, where a free surface condition is implemented.

The numerical resolution of the forward problem amounts to the resolution of a sparse linear system. This is performed through a parallel LU factorization using the MUMPS algorithm [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF]. The LU factorization of the impedance matrix associated with the discretization of the forward problem is reused to solve the adjoint problems, as MUMPS o↵ers the possibility of solving the adjoint system once the factorization has been performed. This is especially important when the number of sources is large: the same LU factorization is used to solve the forward and adjoint problems associated with each source.

This interesting feature is one of the reason for working in the frequency domain: provided the LU factorization of the impedance matrix can be stored, this approach largely reduces the computational costs, compared to the time domain approach. In the particular case of the truncated Newton method, this also makes possible to use the factorization of the impedance matrix for the computation of the Hessian-vector products. The resolution of the extra forward and adjoint problems induced by the truncated Newton method can thus be performed without a new factorization of the impedance matrix.

Minimization scheme settings

In the following two tests an estimation of the P-wave velocity model is computed using a FWI scheme. For each of these two tests, we compare the performances of four minimization strategies:

• nonlinear conjugate gradient method;

• l-BFGS method;

• truncated Newton method using the full Hessian operator;

• truncated Newton method using the Gauss-Newton approximation;

In our implementation of the truncated Newton algorithm, we complement the Eisenstat stopping criterion by setting to 30 the maximum number of inner iterations that can be performed in the first case study, and only 3 in the second case study. These rather small values is chosen to enhance the smoothing e↵ect related to the truncation strategy, which is appropriate for the inversion of noisy data.

The stopping criterion we use is the following: the iterations end as soon as

f (m k )/f (p 0 ) < ✏. ( 42
)
The quantity ✏ is set to 10 2 for the first experiment. Such an accuracy can be reached since in this case we work with synthetic data without noise. In the second experiment, we use real noisy data, and ✏ is set to 6 ⇥ 10 1 , which is adapted to the expected decrease of the misfit function.

This stopping criterion is complemented with a maximum number of nonlinear iteration, which is set to 100 for the BP 2004 test case and 20 for the Valhall test case. In addition, if no acceptable step length is found after 20 linesearch iterations, the minimization is stopped and an error flag is returned.

Finally, note that the memory parameter l for the two l-BFGS methods, which corresponds to the number of gradient stored to compute the approximation of the inverse Hessian, is set to l = 20 in the first case and l = 5 in the second case.

Linesearch algorithm and nonlinear conjugate gradient implementation

The linesearch algorithm we use is designed so that the step k computed at each nonlinear iteration satisfies the Wolfe criterion (see equation ( 4)). Because we want the four minimization methods to use the same linesearch algorithm, we implement a particular form of the nonlinear conjugate gradient method. Indeed, as mentioned by [START_REF] Nocedal | Numerical Optimization[END_REF], standard implementations of this algorithm (such as Fletcher-Reeves or Polak-Ribière implementations) require to use a linesearch algorithm satisfying the strong Wolfe conditions to guarantee global convergence toward local minima. We select instead the nonlinear conjugate gradient algorithm proposed by [START_REF] Dai | A nonlinear conjugate gradient method with a strong global convergence property[END_REF], which is compatible with a linesearch process that only enforces the standard Wolfe conditions.

Preconditioner

The four minimization methods use the same preconditioner P ⌫,✓ k . For the two case studies, a trial-and-error approach has been used to determine a common usable values for ✓. Pragmatical values range from 10 1 to 10 5 . In practice, the value ✓ = 10 2 has been used along all the experiments, for all the methods.

The 2004 BP model

Presentation

The 2004 BP model has been originally designed as a benchmark model for testing sub-salt seismic imaging methods [START_REF] Billette | The 2004 BP velocity benchmark[END_REF]. We perform a decimation of the original model taking one parameter value each ten grid points, and we choose a 25 m

Estimated models

The models estimated by the four minimization methods are presented in figure 2. From km 1 to km 6 in the horizontal direction, the top salt-structure is correctly delineated in the four estimations. The reconstruction of the basin between x = 0 km and x = 1 km seems more di cult. This basin is responsible for high amplitude multi-scattered waves di cult to interpret, and is located at one extremity of the acquisition. The nonlinear conjugate gradient method and the l-BFGS method seem to be the most a↵ected by this particular configuration. The geometry of the basin is not recovered and is filled with high amplitude velocities. These perturbations are also responsible for obscuring the sub-salt targets, and creating erroneous slow velocity anomalies.

Conversely, the results provided by the truncated Gauss-Newton or the truncated Newton method seem more reliable. The best estimation is provided by the truncated Newton method. The geometry of the basin is better recovered, and the sub-salt slow velocity anomalies better reconstructed. Note also that only these two methods are able to provide details on the salt structure itself, namely the very small heterogeneities located at x = 2km, z = 2km.

The possible enhancement of the inverse Hessian approximation yielded by the truncated Newton method may explain this improvement in the resolution and the stability of the inversion.

Convergence profiles

The convergence profiles of the four methods are presented for the frequency group 1, 4 and 8 in figure 3. Excepted for the frequency group 1 for which the l-BFGS method converges slightly faster, the nonlinear conjugate gradient and the l-BFGS method converge slower than the two truncated Newton methods in terms of nonlinear iterations. This is satisfactory since each nonlinear iteration is more expensive for the truncated Newton method. In addition, only the two truncated Newton method satisfies the convergence criterion before reaching the maximum number of nonlinear iteration, set to 100.

In terms of number of forward problems resolution, the slowest method is the truncated Newton method. This indicates that the speed-up in terms of nonlinear iteration does not compensate the extra computation cost required. However, as suggested by the profile aspect, reducing the total number of authorized internal iterations could possibly improve the convergence in terms on forward problem resolutions. Indeed, we can detect plates which correspond to non-converging inner iterations, which reach the maximum allowed number of iterations, set to 30 for this experiment. In each case, the fastest method in terms of forward problem resolution is the l-BFGS method. The truncated Gauss-Newton method and the nonlinear conjugate gradient are in between.

The speed-up obtained in terms of convergence with respect to the number of nonlinear iterations does not compensate for the additional cost related to the resolution of the inner linear systems. However, compared to preconditioned gradient-based method, the truncated Newton method appears to provide better subsurface estimations.

2D Valhall case study

Description

The Valhall oil field is located in the North Sea, in production since the beginning of the 80s. In this shallow water environment field, the water depth reaches only 100 m. This particular configuration is adapted for the use of ocean bottom cables (OBC) seismic recorders.

Four components sensors have thus been disposed at the sea bottom level. Several threedimensional seismic surveys have been performed to follow the time-evolution of the oil field. The Valhall seismic data have been investigated by several authors in the context of 2D multi-parameter FWI [START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[END_REF](Prieux et al. , 2013a,b) ,b) and 3D mono-parameter FWI [START_REF] Sirgue | Full waveform inversion: the next leap forward in imaging at Valhall[END_REF][START_REF] Etienne | Threedimensional acoustic full waveform inversion: algorithm and application to Valhall[END_REF].

In this study, we are interested in the mono-parameter inversion of a 2D line from the Valhall OBC data-set. This data-set involves 320 sources located at 5 meters depth in the water layer, and 210 receivers located on the sea bottom, between 68 and 72 meter depth.

Sources and receivers are equally spaced each 50 m. Only the hydrophone component is used. The data is interpreted down to 3.5 Hz. The signal-over-noise ratio (SNR) for lower frequencies is too weak for the data to be exploited in this frequency band.

Below the sea level, the presence of soft shale sediments is responsible for strong attenu- The VTI forward modeling requires to define initial models not only for P-wave velocity, density and attenuation, but also for the Thomsen anisotropy parameters ✏ and . The initial models for P-wave velocity and Thomsen parameters are determined by reflection travel-time tomography. The density initial model is derived from the P-wave velocity initial model through the Gardner law. The quality factor model is taken constant equal to 200.

The corresponding models are presented in figure 4. In our experiments, the initial models for density, attenuation and Thomsen parameters remain the same, and we only focus on the reconstruction of the P-wave velocity model.

Regularization strategy

As the real seismic data is noise contaminated, an appropriate regularization strategy has to be designed. We choose to add a standard Tikhonov regularization term T (m) to the misfit function f (m) to define the regularized misfit function f T (m)

f T (m) = f (m) + ↵ 2 T (m) ( 4 3 )
where

T (m) = ↵ x k@ x mk 2 + ↵ z k@ z mk 2 (44)
The parameter ↵ accounts for the influence of the regularization term. The parameters ↵ x and ↵ z can be used to enforce a stronger regularization in the direction x or z. In the sequel, we use the settings

↵ x = 1, ↵ z = 0.5 ( 4 5 )
to enforce smoother variations in the horizontal direction.

In addition, the truncation strategy for the computation of the descent direction within the truncated Newton/Gauss-Newton methods is used as an additional regularization (see the work of [START_REF] Kaltenbacher | Iterative Regularization Methods for Nonlinear Problems[END_REF] on this particular topic). The maximum number of inner conjugate gradient iteration is set to 3.

In order to investigate the sensitivity of the methods to the regularization parameter ↵, we have performed experiments for 2 di↵erent values of ↵, reflecting what we may call a "strong" ↵ 1 and a "weak" ↵ 2 regularization.

• For nonlinear conjugate gradient and l-BFGS, we use ↵ 1 = 10 4 and ↵ 2 = 5 ⇥ 10 4 .

• For the truncated Newton/Gauss-Newton method, we use ↵ 1 = 10 6 and ↵ 2 = 5⇥10 6 .

The combination of two types of regularization for the truncated Newton/Gauss-Newton strategies makes possible to decrease the value of the parameter ↵ by two orders of magnitude for these methods. This reduces the smoothing e↵ect yielded by the Tikhonov regularization term, and possibly the resolution loss associated with this smoothing.

FWI results

We presents the results obtained by our four optimization schemes in figures 5 and 6. Below shallow layers of shale sediments (from z = 0 to z = 1.5 km), we can see two low velocity zones, corresponding to the presence of gas layers (between z = 1.5 km and z 2.5 km,

x = 8 km and x = 12 km). Below these gas layers (between z = 2.5 km and z = 4 km), we can locate the cap rock of the reservoir and stronger reflectors.

The results we obtain demonstrate the di culty of finding a suitable trade-o↵ between regularization and resolution power. For nonlinear conjugate gradient and the l-BFGS method, changing from weak to strong regularization yields significant di↵erences in the results. For weak regularization, the final results is contaminated by noise. For strong regular- ization, the results appears significantly smoother, to the cost of a potential loss of resolution in depth.

Conversely, the truncated Newton method seems to be more stable with respect to a modification of the regularization parameter. The two results obtained for weak and strong regularization do not show strong di↵erences. In addition, even if the regularization parameter used is two order of magnitude smaller than the one used for the two latter methods, the final results is less a↵ected by noise. This tends to confirm the regularization e↵ect associated with the truncation strategy.

This analysis is confirmed by the well-logs presented in figures 7 and 8. The results from the truncated Newton methods provide a globally better fit to the well-log data with strong or weak regularization. We are particularly interested in the jump in depth at z = 3.5 km (green circles), which can be quite clearly identified in the results provided by these two methods. This jump correspond to the presence of an actual reflector (see the migrated section using the initial velocity model 9). For the preconditioned gradient-based methods, the results obtained with a weak regularization present strong oscillations in depth which make di cult to identify this reflector. The results obtained with a strong regularization smooth out the results, and make equivalently di cult the identification of the reflector.

Resolution analysis

The analysis we have led suggests that a better trade-o↵ between regularization and resolution power can be obtained using the truncated Newton method. In order to investigate this particular point we compare the sensitivity of the preconditioned gradient-based methods and the truncated Newton method to a perturbation in the final estimations.

We add a positive perturbation of 200 m.s 1 in the gas cloud localized around x = 11 km, z = 2 km. We compute synthetic data using the same surface acquisition geometry in this perturbed model. We invert this synthetic data starting from the unperturbed model, and we perform only one nonlinear iteration. In this case, the nonlinear conjugate gradient and the l-BFGS method are equivalent: at the first iteration, the l-BFGS approximation is just given by the preconditioning matrix.

We compare the reconstruction of the perturbation for the nonlinear conjugate gradient model and for the truncated Newton model. For the truncated Newton experiment, we compute the perturbation using

• 3 inner linear iterations

• 7 inner linear iterations

The results are presented in figure 10. Not surprisingly, the focusing achieved by the nonlinear conjugate gradient method is poorer than the one obtained using the truncated Newton method. The amplitude of the perturbation reconstructed by the nonlinear conjugate gradient only reaches 10 m.s 1 while it reaches more than 30 m.s 1 for the truncated Newton method. More interestingly, we can see the refocusing e↵ect associated with the inverse Hessian operator on the model update obtained after 7 inner linear iterations. The amplitude of the reconstructed perturbation is similar, but the artifacts around the perturbation tends to vanish.

Convergence profiles

The convergence profiles for the Valhall case study are displayed in figure 11. They correspond to the inversion of the frequency groups 1, 3, and 5 in the strong regularization case.

Among the four methods, the l-BFGS method seems to be the less robust, as it fails to converge for the three frequency groups displayed. It means that the method terminates on a linesearch failure: the minimization of the misfit function does not reach the bound fixed at 0.6 or the maximum number of nonlinear iteration fixed at 20.

Compared to l-BFGS, the nonlinear conjugate gradient method seems to be more robust.

For the frequency group 1 and 3, the maximum number of nonlinear iterations is reached and the method provides a better decrease of the misfit function.

The truncated Newton/Gauss-Newton methods systematically reach the maximum number of authorized nonlinear iterations and provide a better decrease of the misfit function at each frequency group. However, the speed-up in terms of nonlinear iteration compared to the gradient-based methods is limited, excepted for the first frequency group.

In terms of number of resolution of wave propagation problems, the truncated Newton/Gauss-Newton methods appear not surprisingly as more expensive as their gradient-based counterpart. However, this extra computation cost is controlled by the maximum number of inner iteration, which is limited to 3. Note that for this case study, the results provided by the truncated Newton and the truncated Gauss-Newton are very similar (for the first frequency group, the convergence curves of these two methods are even superposed).

CONCLUSION AND PERSPECTIVES

Conclusion

In this study, we are interested in the implementation of the truncated Newton minimization scheme for FWI. The algorithm is based on a two nested loops architecture: the external loop consists in updating an initial model up to the final estimation. The inner loop consists in the computation of the model updates from an inexact resolution of the linear system associated with the computation of the Newton descent direction. This inexact resolution is performed with a matrix-free conjugate gradient solver. Only the action of the Hessian operator on an arbitrary vector has to be computed. This can be e ciently achieved using second-order adjoint state formulas: only two additional wave propagation problems have to be solved to compute this quantity. The use of an adaptive stopping criterion for the inner loop and a suitable preconditioner enhance the e ciency of the algorithm.

Compared to conventional preconditioned gradient-based methods operator such as the nonlinear conjugate gradient or the l-BFGS method, it seems that using the truncated Newton method could be advantageous. The synthetic BP 2004 case study emphasizes the di culty of salt and sub-salt imaging. The high velocity contrasts caused by the presence of large salt structures is responsible for the presence of energetic multiple reflections in the data. This cause di culties to standard minimization methods to provide correct images of the salt and sub-salt targets. In this case, the truncated Newton method seems to be more robust than the nonlinear conjugate gradient and the l-BFGS method. In particular, the truncated Newton method provide better results than the truncated Gauss-Newton method. A possible explanation is that the part of the Hessian operator which is neglected in the Gauss-Newton approximation is related to double scattered waves. In this context of multiple-scattering between the salt and the free surface, it could be crucial to account for the whole Hessian operator.

The Valhall real data case study also demonstrates the potential interest of using the truncated Newton method compared to the conventional methods. The noise-contamination of the data requires to use a suitable regularization. We use a standard Tikhonov regularization of the first order derivatives of the model. We identify the trade-o↵ between regularization and the expected resolution of the final estimation. The nonlinear conjugate gradient and the l-BFGS method are highly sensitive to the regularization hyper-parameter which controls the amount of regularization. Conversely, the truncated Newton method possesses an inherent smoothing related to the truncation strategy for computing the model update [START_REF] Kaltenbacher | Iterative Regularization Methods for Nonlinear Problems[END_REF]. This yields the possibility of choosing a regularization hyper-parameter 2 orders of magnitude lower than for the preconditioned gradient-based methods. It also reduces the sensitivity to this hyper-parameter. A better trade-o↵ between resolution and regularization is thus achieved. A point-spread function test in the final estimated models clearly indicates that the information on the local curvature carried out by the inverse Hessian operator yields a better resolution.

Perspectives

The numerical performances of the di↵erent minimization schemes we have compared strongly depend on an accurate estimation of the inverse Hessian operator as a preconditioner. For now, we have only used the diagonal preconditioner from the pseudo-Hessian approach developed by [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF]. The accuracy of this approximation is however limited and the use of this preconditioner is restricted to surface acquisition configurations. Therefore we are interested in the development of better approximation of the inverse Hessian. Several techniques could be used.

For instance, [START_REF] Chiu | Matrix probing and its conditioning[END_REF] have proposed approximation of the inverse Hessian operator through matrix probing. However, this method only approximates the Gauss-Newton part of the operator, and relies on the acoustic approximation. [START_REF] Bekas | An estimator for the diagonal of a matrix[END_REF] propose a general algebraic method to compute approximation of the diagonal of a given matrix in a matrix-free fashion. The explicit construction of the entire matrix is not required. This makes this method adapted to FWI, as opposed with incomplete LU or Cholesky factorization [START_REF] Benzi | Preconditioning techniques for large linear systems: A survey[END_REF].

Another way of improving the resolution of the inner linear systems is related to deflation strategies [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. From the Krylov subspace created during the resolution of the inner linear system at iteration k 1, an estimation of the smallest eigenvectors of the Hessian operator can be computed. These eigenvectors can be inserted in the starting Krylov suspace at iteration k in order to improve the convergence. Related methods developed within the data assimilation community could also be investigated. [START_REF] Gratton | A reduced and limited-memory preconditioned approach for the 4d-var data-assimilation problem[END_REF] propose for instance to compute an approximate initial guess for the inner linear systems in a reduced space given by a spectral decomposition of the initial Hessian operator.

The importance of an accurate estimation of the inverse Hessian operator should be even more important in the context of multi-parameter FWI. The simultaneous reconstruction of parameters such as the P-wave velocity and the density for instance is a↵ected by strong trade-o↵s. The perturbation in the P-wave models can be erroneously interpreted as perturbations in the density model, and vice-versa. A more accurate estimation of the inverse Hessian operator should help to mitigate these trade-o↵ e↵ects. In this context, the truncated Newton method could thus produce more reliable results. This topic will be investigated in future studies.

Figure 1 .

 1 Figure 1. BP 2004 synthetic exact model (a), initial model (b).

Figure 2 .

 2 Figure 2. Estimated models for the BP case study. Nonlinear conjugate gradient (a), l-BFGS method (b), truncated Gauss-Newton method (c), truncated Newton method (d).

Figure 3 .

 3 Figure 3. Misfit function decrease for the BP 2004 case study. Decrease with respect to the the number of nonlinear iterations: first frequency group (a), second frequency group (b), third frequency group (c). Decrease with respect to the number of wave equation problem solved: first frequency group (d), second frequency group (e), third frequency group (f). PCG: nonlinear conjugate gradient, PLB: l-BFGS, PTN: truncated Newton, PTGN: truncated Gauss-Newton method.

  ation and anisotropy of the wave propagation. As demonstrated in the work of[START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[END_REF], it is necessary to account for this anisotropy in the modeling to correctly invert the Valhall data-set. For this case study, we thus change the forward problem from standard isotropic acoustic modeling to the VTI acoustic modeling proposed by[START_REF] Operto | Finite-di↵erence frequencydomain modeling of visco-acoustic wave propagation in two-dimensional TTI media[END_REF].

From 3. 5

 5 Hz to 7.45 Hz, we define 6 overlapping frequency groups to be inverted sequen-

Figure 4 .

 4 Figure 4. Initial models for P-wave velocity (a), density (b). Dimensionless Thomsen parameters ✏ (c), (d)

Figure 5 .

 5 Figure 5. Weak regularization: computed P-wave velocity models using the nonlinear conjugate gradient (a), l-BFGS (b), truncated Gauss-Newton (c), truncated Newton (d).

Figure 6 .

 6 Figure 6. Strong regularization: computed P-wave velocity models using the nonlinear conjugate gradient (e), l-BFGS (f), truncated Gauss-Newton (g), truncated Newton (h).

Figure 7 .

 7 Figure 7. P-wave velocity logs compared to a reference well-log at x = 9.5 km, ), with weak regularization. Results obtained using the nonlinear conjugate gradient (a), l-BFGS (b), truncated Gauss-Newton (c), truncated Newton (d)

Figure 8 .

 8 Figure 8. P-wave velocity logs compared to a reference well-log at x = 9.5 km, ), with strong regularization. Results obtained using the nonlinear conjugate gradient (a), l-BFGS (b), truncated Gauss-Newton (c), truncated Newton (d)

Figure 9 .

 9 Figure 9. Reverse Time Migration result using the starting velocity model.

Figure 10 .

 10 Figure 10. Reconstruction of a small amplitude perturbation in the final estimation. Zoom on the perturbation. Nonlinear conjugate gradient (a), truncated Newton method with 3 inner linear iterations (b), 7 inner linear iterations (c)

Figure 11 .

 11 Figure 11. Misfit function decrease for the Valhall case study for the strong regularization case. Decrease with respect to the number of nonlinear iterations: frequency group 1 (a), frequency group 3 (b), and frequency group 5 (c). Decrease with respect to the number of wave equation problems solved: frequency group 1 (d), frequency group 3 (e), and frequency group 5 (f). PCG: nonlinear conjugate gradient, PLB: l-BFGS, PTGN: truncated Gauss-Newton method, PTN: truncated Newton.

Table 1 .

 1 Frequency group strategy for the BP 2004 case study 2012), we use here a smoother initial model. The characteristic length used for this smoothing is set to 500 m instead of 375 m for this former study. The resulting initial model is

presented in figure

1

. The use of this smoother initial model requires a careful hierarchical frequency strategy. We generate 27 data sets, from 2 Hz frequency to 19.5 Hz gathered into 8 overlapping subgroups, as presented in table 1.

Table 2 .

 2 Frequency group strategy for the Valhall case study

? The Hessian is the matrix of the second-order derivatives of the misfit function.

ACKNOWLEDGMENTS

This study was funded by the SEISCOPE II consortium (http : //seiscope2.osug.f r), sponsored by BP, CGG-VERITAS, CHEVRON, EXXON-MOBIL, JGI, SAUDI ARAMCO, SHELL, STATOIL and TOTAL. This study was granted access to the HPC facilities of CIMENT (Université Joseph Fourier Grenoble), and of GENCI-CINES under Grant 2011-046091 of GENCI (Grand Equipement National de Calcul Intensif). The authors would like to thank V. Tcheverda for his invitation to publish this study. They also thank BP Norge AS and their Valhall partner Hess Norge AS, for allowing access to the Valhall dataset as well as the well-log velocities.


Najo A Jomaa

Marwan Haddad

Grace Y Adwane

Amira J Zaylaa

M.D Abdallah Rahbani

Dedicated Magnetic Resonance Imaging Sequences: Contribution in the Diagnosis of Focal Epilepsy in the Lebanese Population

Keywords: Epilepsy, magnetic resonance imaging, dedicated sequences, focal, epileptogenic lesions, hippocampal sclerosis, cortical dysplasia

There have been many advances in the diagnosis and management of focal epilepsies particularly with neuroimaging techniques and magnetic resonance imaging (MRI). Special MRI sequences can be employed to localize and resect the epileptogenic lesions responsible for focal epilepsy. This study aims to show the benefit of dedicated epilepsy MRI sequences in the diagnosis of focal epilepsies. A general electric 1.5 Tesla MRI machine was used with standard and special sequences. Also, a Nihon Kohden electroencephalography (EEG) machine was used. This is a prospective observational study that included 51 patients with focal epilepsies who had an initial negative brain imaging. They underwent epilepsy MRI sequences along with a prolonged video EEG monitoring to localize the lesion, and then results were analyzed statistically using SPSS 22 program. The majority of patients were males (62.75%) with a mean age of 30 years. The grand majority of patients (74.5%, p value of 0.001) had their lesion localized by the epilepsy MRI. The most commonly found pathology was mesial temporal sclerosis. A significant number of patients (23.5%) were sent for an epilepsy surgery (p value 0.002). This study shows the significant impact of dedicated epilepsy MRI sequences on the diagnosis and management of focal epilepsy in the Lebanese population.

Introduction

An epileptic seizure is defined as a transient alteration in mentation associated with abnormal movements as a result of excessive bursts of brain electrical activity [START_REF] Fisher | The New Classification of Seizures by the International League Against Epilepsy 2017[END_REF]. However, epilepsy itself is a chronic condition in which seizures recur, and this can During the last decade, there have been many advances in our understanding and management of epilepsies. With the emerging neuroimaging modalities including specialized magnetic resonance imaging (MRI) sequences, the total number of presumed non-idiopathic cryptogenic partial epilepsies is drastically decreasing [START_REF] Degnan | Neuroimaging of Epilepsy: A Review of MRI Findings in Uncommon Etiologies and Atypical Presentations of Seizures[END_REF][START_REF] Ramli | Neuroimaging in Refractory Epilepsy: Current Practice and Evolving Trends[END_REF]. be complicated by neurobiological, cognitive, psychosocial, and occupational consequences [START_REF] Athira | Role of Magnetic Resonance Imaging Brain in Evaluation of Seizures[END_REF][START_REF] Hartman | Epilepsy Imaging: Approaches and Protocols[END_REF].

The 2017 International League Against Epilepsy (ILAE) classification divides epilepsy into 3 types: those with focal onset, generalized onset, and unknown onset. Focal seizures were further divided into those with and without an impaired awareness [START_REF] Fisher | The New Classification of Seizures by the International League Against Epilepsy 2017[END_REF]. Focal epilepsy is caused by an abnormal electrical discharge from an epileptogenic zone, which is the cortical area indispensable for the seizure generation [START_REF] Friedman | Epilepsy Imaging in Adults: Getting It Right[END_REF].

The majority of epileptic patients are satisfactorily treated medically with antiepileptic drugs, but around 15-30% continue to have seizures despite therapy [START_REF] Ramli | Neuroimaging in Refractory Epilepsy: Current Practice and Evolving Trends[END_REF][START_REF] Athira | Role of Magnetic Resonance Imaging Brain in Evaluation of Seizures[END_REF][START_REF] Hartman | Epilepsy Imaging: Approaches and Protocols[END_REF].

Epilepsy is usually considered to be refractory if it cannot be controlled after the use of multiple drugs for more than 2 years [START_REF] Saini | Susceptibility Weighted Imaging in the Diagnostic Evaluation of Patients with Intractable Epilepsy[END_REF][START_REF] French | Refractory Epilepsy: Clinical Overview[END_REF]. In this case, surgery is a major therapeutic alternative that can render 70% of patients seizure free [START_REF] Athira | Role of Magnetic Resonance Imaging Brain in Evaluation of Seizures[END_REF][START_REF] Saini | Susceptibility Weighted Imaging in the Diagnostic Evaluation of Patients with Intractable Epilepsy[END_REF][START_REF] Vonoertzen | Standard Magnetic Resonance Imaging Is Inadequate for Patients with Refractory Focal Epilepsy[END_REF]. This is much easier to obtain when the lesion is well identified on brain imaging prior to the surgery [START_REF] Friedman | Epilepsy Imaging in Adults: Getting It Right[END_REF]. Thus, neuroimaging is an essential step in the presurgical workup of localization-related epilepsies in order to well delineate the epileptogenic zone [START_REF] Ramli | Neuroimaging in Refractory Epilepsy: Current Practice and Evolving Trends[END_REF][START_REF] Saini | Susceptibility Weighted Imaging in the Diagnostic Evaluation of Patients with Intractable Epilepsy[END_REF].

MRI is considered to be the imaging modality of choice in defining the structural lesions that are responsible for seizures [START_REF] Hartman | Epilepsy Imaging: Approaches and Protocols[END_REF]. This is due to its better soft tissue contrast, multiplanar imaging and safety profile [START_REF] Ramli | Neuroimaging in Refractory Epilepsy: Current Practice and Evolving Trends[END_REF][START_REF] Chang | Call for Quality Epilepsy Neuroimaging-An Illustrative Case of Missed Epileptic Substrate[END_REF].

Earlier studies have shown that around 20% of epilepsies have no lesion that is identified on MRI [START_REF] Degnan | Neuroimaging of Epilepsy: A Review of MRI Findings in Uncommon Etiologies and Atypical Presentations of Seizures[END_REF][START_REF] Tellez-Zenteno | Surgical Outcomes in Lesions and Nonlesional Epilepsy: A Systemic Review and Meta-Analysis[END_REF]. However, a study done on children having refractory epilepsy found that a high-resolution MRI was able to detect lesions not previously detected on a standard MRI in more than half of the cases [START_REF] Ramli | Neuroimaging in Refractory Epilepsy: Current Practice and Evolving Trends[END_REF].

Important pathologies that must be addressed in focal epilepsy include Mesial Temporal Sclerosis (MTS), Focal Cortical Dysplasias (FCD), developmental cortical malformations, arteriovenous malformation, tumors and infections [START_REF] Hartman | Epilepsy Imaging: Approaches and Protocols[END_REF][START_REF] Friedman | Epilepsy Imaging in Adults: Getting It Right[END_REF][START_REF] Colombo | Focal Cortical Dysplasias: MR Imaging, Histopathologic, and Clinical Correlation in Surgically Treated Patients with Epilepsy[END_REF][START_REF] Colombo | Focal Cortical Dysplasia Type IIa and IIb: MRI Aspects in 118 Cases Proven by Histopathology[END_REF][START_REF] Turjman | Epilepsy Associated with Cerebral Arteriovenous Malformations: A Multivariate Analysis of Angioarchitectural Characteristics[END_REF]. Of these underlying pathologies, MTS is the most common entity in the drug-resistant epilepsy [START_REF] Zentner | Clinical Relevance of Amygdala Sclerosis in Temporal Lobe Epilepsy[END_REF].

In this regard, dedicated MRI sequences have been developed. These include high-resolution thin cuts using coronal FLAIR, T2, T2* Gradient Echo and 3D sequences [START_REF] Athira | Role of Magnetic Resonance Imaging Brain in Evaluation of Seizures[END_REF][START_REF] Hasboun | MR Determination of Hippocampal Volume: Comparison of Three Methods[END_REF][START_REF] Urbach | MR Imaging in the Presurgical Workup of Patients with Drug-Resistant Epilepsy[END_REF][START_REF] Urbach | Imaging of the Epilepsies[END_REF].

The role of Electroencephalography (EEG) is also a major one as the presence of a lesion on MRI does not necessarily indicate its epileptogenicity, unless the result is congruent with the EEG findings [START_REF] Ramli | Neuroimaging in Refractory Epilepsy: Current Practice and Evolving Trends[END_REF].

It is important to note that many patients with epilepsy frequently undergo standard MRI or CT in their initial workup in non-specialized centers and many of them have a missed diagnosis [START_REF] Vonoertzen | Standard Magnetic Resonance Imaging Is Inadequate for Patients with Refractory Focal Epilepsy[END_REF]. In Lebanon and other developing countries, this scenario is of a higher occurrence because of the limited resources in the management of epilepsy [START_REF] Chang | Call for Quality Epilepsy Neuroimaging-An Illustrative Case of Missed Epileptic Substrate[END_REF].

The aim of our study is to show the additional value of dedicated epilepsy MRI sequences in the detection of brain lesions in patients having cryptogenic focal epilepsy in the Lebanese population. Our study will hopefully guide the medical professionals to do a more extended work-up of the patients presenting with focal seizures.

Materials

The materials we used in this study are mentioned below. These included also the different types of MRI sequences we used for imaging.

(1) A 1. In case a neoplasm or neurocutaneous syndrome was suspected a contrast agent was injected intravenously: Gadodiamide (Omniscan) at a dose of 0.1 mmol/kg of body weight.

(2) A 128-channel video-EEG monitoring machine in an epilepsy monitoring unit. It is used as a part of a comprehensive epilepsy program at the Lebanese Geitaoui University Hospital.

Methods

According to our exploratory study, we took into consideration the following patients demographics, ethical considerations and data collection methods.

Patients Demographics

This is a prospective study that included patients of all ages whether males or females, having focal epilepsy with a negative initial brain imaging.

Patients were enrolled in the study according to the following inclusion and exclusion criteria.

Inclusion Criteria:  Partial epilepsy;  Non-idiopathic;  Negative baseline neuroimaging (CT scan or routine MRI).

Exclusion Criteria:  Generalized epilepsy;  Idiopathic;  Positive baseline imaging for a pathologic lesion.

Ethical Consideration

Patients included in this study were not subjected to any new procedure or medication. They underwent the basic exams that comprehensive epilepsy programs would include in order to diagnose and treat epilepsy. Thus subjects were not at risk of harm by any step of the study. Besides, confidential and private information concerning the patients was not exposed to anyone else but the treating physician.

Collection of Data

This study is a prospective observational study in which we included 51 patients presenting with focal epilepsy and a negative baseline brain imaging (CT scan or MRI). These patients were enrolled over a period of around 5 months extending from April to August 2017.

All the patients underwent an initial prolonged video-EEG monitoring in the setting of a university hospital and were included in a comprehensive epilepsy program. This was done with the aim to detect a clinical seizure or electrical epileptiform discharges. Then, patients were referred to a specialized radiology department to do dedicated MRI epilepsy sequences. MRI sequences were analyzed by the radiologists to detect a focal pathologic lesion. These results were then compared to the initial imaging results in order to analyze statistically the significance of these specialized sequences in focal epilepsy. Finally a PET scan was ordered for those patients who had a negative dedicated MRI with a high suspicion for a focal lesion (See Fig. 1).

Statistical Analysis Method

In our study, all the statistical analysis was carried out with SPSS Statistics version 22 (IBM, Inc, Armonk, NY). We calculated the mean only for the age. For the other categorical data, we calculated frequencies and percentages. For the p-value, we used a one sample t-test where a p-value of less than 0.05 was considered to be significant.

Results

In this section below, the different results regarding the distribution of the types of seizures, pathologies, and their correlation with MRI results in addition to the different patients' characteristics will be presented.

Patients' Characteristics

A total of 51 patients with focal epilepsy were included in this study. Our patients were divided into 32 males (62.7%) and 19 females (37.2%). Ages ranged between 2 years and 85 years with a mean age of 30 years.

Types of Seizures

Upon presentation, the majority of patients had temporal electrical discharges on the EEG (37.25%), followed in frequency by frontal discharges (21.56%). Moreover, 7 patients out of 51 had multifocal discharges (13.73%) and 3 patients had epileptic encephalopathy on the EEG (5.88%) (Table 1).

Epilepsy MRI Results

Positive MRI Results

Out of the 51 patients included in the study, 38 patients (74.5%) had a confirmed diagnosis by an epilepsy MRI which is statistically significant (p-value = 0.001). In these patients, specialized MRI sequences were able to detect the underlying epileptogenic lesion responsible for their epilepsy.

Types of Pathologies

The majority of lesions that were detected on the specialized MRI consisted mainly of hippocampal sclerosis and cortical dysplasias. A wide variety of other lesions were also detected ranging from congenital abnormalities such as heterotopias and polymicrogyrias to low grade tumors such as DNET and gliomas (Table 2).

One of the patients included in the study was found to have a dual pathology on imaging. These consisted of hippocampal sclerosis and focal cortical dysplasia (Figs. 2 and3).

The majority of the detected pathologies, as shown in Table 3, are non-neoplastic.

Table 3 demonstrates the difference in number between congenital and acquired causes. As shown below, the majority of pathologies were consisting of acquired causes.

Age Distribution versus Pathology

The great majority of patients who were found to have a focal cortical dysplasia in our study were in the first or second decade, while almost all of the patients with mesial temporal sclerosis were above 20 years of age.

4.5Positron Emission Tomography (PET) Scan

Four patients out of the 12 cases with a negative MRI underwent PET scan and were shown to have an epileptogenic lesion. One of these underwent surgical resection of the lesion and another patient underwent a VNS (Vagal Nerve Stimulation) surgery.

Disease Management

Thirteen patients were sent for an epileptic surgery based on their dedicated MRI results (12 patients, 23.5%), and 3 patients were sent for VNS (5.8%). Almost all of these patients are still seizure free postoperatively. Another larger group of patients (22 patients, 43.1%) with positive MRI were managed with an appropriate medical treatment that was suitable for their underlying lesion.

Discussion

Epilepsy is still one of the major challenging Central Nervous System (CNS) disorders. Reaching a definite diagnosis has a key value for the prognosis and therapeutic decisions. This is largely dependent on the clinical history and findings, along with a dedicated and specialized neuroimaging and electrophysiologic procedures.

Our results showed that when the clinical picture was highly in favor of a non-idiopathic focal epilepsy, dedicated MRI sequences for epilepsy were able to well localize epileptogenic lesions in 74.51% of the studied patient population in whom the basic neuroimaging studies were not conclusive. This number is considered significant (with a p-value of 0.001) and it is associated with a large impact on the diagnostic and therapeutic approach of epilepsy patients, especially in the Lebanese population.

This fact was highlighted by the significant number of patients who underwent a major change in their therapeutic strategy. This mainly consisted of epilepsy surgery and VNS.

It is important to note that patients were selected for surgery based on many factors including the operability of the lesion (whether it is present in an eloquent cortex or not) in addition to the validity and correlation of the clinical, radiological, and electrophysiological findings, i.e. patients lacking this correlation were considered for much further work up (for example PET scan) and were treated medically.

Our results were congruent with the previous studies showing that mesial temporal sclerosis is still the most frequent pathology detected on brain MRI in focal epilepsy patients. Also similar to previous studies, the majority of the detected lesions were acquired (60.53%) and the majority of these too were non-neoplastic (81.55%).

Our results showed the interest of asking for an epilepsy protocol MRI in patients presenting with seizures (especially in those presenting with focal seizures) by their primary care physicians, neurologists and epileptologists. This is important considering the fact that many patients in our country Lebanon and other developing countries are not undergoing this basic step in managing their epilepsy. This is mainly due to lack of insight on the importance of these very specific techniques, or due to the lack of financial medical coverage of patients. Many patients for example seek medical care and imaging in nonspecialized centers for financial reasons.

On the other hand, many centers in Lebanon are already equipped with high resolution MRI machines that are able to perform all the previously mentioned specialized epilepsy sequences. This is in addition to the fact that also a large number of radiologists are very well trained in this regard.

Conclusion

Refractory seizures are a major cause of morbidity in epileptic patients causing devastating consequences on several aspects of their lives. Knowing that a significant number of epileptic patients have focal epilepsy makes it very important to well localize the causative epileptogenic lesion. This would help in referring patients for a more effective surgical therapy.

This study showed the large influence a dedicated epilepsy MRI can have in the diagnosis and management of epileptic patients with focal epilepsy.

This study showed that a large number of patients in the Lebanese population can be falsely diagnosed as having cryptogenic epilepsy with an unknown causative lesion.

In this regard, it is important to send these patients to do a specific epilepsy protocol MRI instead of nonspecific brain imaging in their initial work up.

Fig. 1

 1 Fig.1Steps followed in the methodology of this study.

Fig. 2 AFig. 3 A

 23 Fig. 2 A coronal FLAIR sequence with thin cuts showing bilateral hippocampal sclerosis (white arrows).

Table 1 Distribution of the types of seizures or discharges present on EEG.

 1

	Type of seizure	Number of studied EEGs (n = 51)	Percentage %
	Temporal	19	37.25
	Frontal	11	21.56
	Multifocal	7	13.73
	Parietal	6	11.76
	Epileptic encephalopathy 3	5.88
	Occipital	2	3.92
	Gelastic	1	1.96
	Psychomotor	1	1.96
	Insular	1	1.96

Table 2 Different types of epileptogenic lesions detected on specialized MRI.

 2

		Number of studied
	Type of lesion	specialized MRIs	Percentage %
		(n = 38)
	Hippocampal sclerosis	10	26.3
	Cortical dysplasias	9	23.7
	Infarct with gliosis	4	10.5
	DNET	3	7.9
	Polymicrogyria	2	5.3
	Tuberous sclerosis	2	5.3
	Low grade glioma	2	5.3
	Meningioma	1	2.6
	Hypothalamic hamartoma 1	2.6
	Sturge weber	1	2.6
	Posterior reversible
	encephalopathy syndrome	1	2.6
	(PRES)
	Heterotopia	1	2.6
	Porencephalic cavity	1	2.6

Table 3 The distribution of the detected lesions between neoplastic, non-neoplastic etiologies, congenital and acquired causes.

 3

	Type of etiology	of studied etiologies (n = 38)	Percentage %
	Neoplastic	7	18.42
	Non-neoplastic 31	81.58
	Congenital	15	39.47
	Acquired	23	60.53