Hélène Barucq

Henri Calandra

Guy Chavent

Florian Faucher
email: florian.faucher@inria.fr

A priori estimates of attraction basins for velocity model reconstruction by time-harmonic Full Waveform Inversion and Data Space Reflectivity formulation

Keywords: Inverse problem, optimization, convergence analysis, geophysics, wave equations, FWI, MBTT, Helmholtz equation, reconstruction algorithm Problème inverse, optimisation, analyse de convergence, géophysique, équation des ondes, FWI, MBTT, algorithme de reconstruction

Determination of background velocity by Full Waveform Inversion is known to be hampered by the local minima of the data misfit caused by the phase shifts associated to background perturbations. We give in this paper an a priori study of the convergence basins around any nominal velocity model, where the data misfit is guaranteed to have only one localand hence global -minimum, provided the data error is below a given tolerable error level. The size of the attraction basin gives an upper bound to the error in the initial velocity model which ensures that local optimization algorithm will converge to the global mimimum. We use geometry in the data space in order to compute this size and the associated tolerable error without having to solve any FWI least squares minimization problem. Application is made to the inverse problem associated with wave equations in the time-frequency domain, using reflection data. We study first how, when FWI is used, frequency and search direction geometry influence the size of attraction basins, and how to foster it (e.g. complex frequency progression). Then, we compare the size of attraction basins for background velocities for the classical FWI and the Data Space Reflectivity MBTT formulation. We investigate the domain of validity of the MBTT reformulation, and show that it can increase the size of the background attraction basins by a factor of four to fifteen.

Estimations a priori de bassins d'attraction pour la reconstruction de modèles de vitesse pour l'équation de Hemholtz par inversion de forme d'onde (FWI) et par la formulation stockant la réflectivité dans l'espace des données (MBTT)

The Full Waveform Inversion (FWI) approach to seismic inversion consists in minimizing the misfit between the recorded and synthetic data with respect to the Earth parameters. It uses the full waveform of the recorded waves as data, which makes it free of horizon picking, and hence a good candidate for automatic inversion. It was first introduced by Bamberger et al. [START_REF] Bamberger | Une application de la théorie du contrôle à un problème inverse de sismique[END_REF][START_REF]About the stability of the inverse problem in the 1-d wave equation[END_REF] for one-dimensional model, and followed, in two dimensions, by Lailly [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF] and Tarantola [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF]. In particular, [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF] made the link between the gradient of the misfit functional, computed by the adjoint state method, and the well-known migration imaging operator used by geophysicists. Early numerical experimentations in [START_REF] Tarantola | Inversion of travel times and seismic waveforms[END_REF][START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms; numerical results[END_REF] already show the strength and the weakness of FWI: it can produce very clear and detailed images of the sub-subsurface (strength), under the imperative condition that a good background velocity model is available (weakness). In fact, the misfit functional with respect to background parameters shows local minima due to phase shift/cycle skipping, which hamper the reconstruction by local algorithms, unless the minimization starts inside the attraction basin. FWI is extended to elastic waves in [START_REF] Mora | Nonlinear two-dimensional elastic inversion of multioffset seismic data[END_REF], with implementation in [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF], to resolve the short wavelengths of P-and S-impedance. Initially carried out in the time-domain, FWI can also be used in the frequency domain, as promoted by Pratt and collaborators in [START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography[END_REF][START_REF] Pratt | Combining wave-equation imaging with traveltime tomography to form high-resolution images from crosshole data[END_REF][START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by wavefield inversion[END_REF]. The use of complex frequencies was later introduced by [START_REF] Shin | Waveform inversion in the laplace domain[END_REF][START_REF]Waveform inversion in the laplace fourier domain[END_REF], giving rise to the so-called Laplace and Laplace-Fourier domain FWI.

In the recent years, with the increase of the computational power, the solution of the FWI approach to seismic imaging by local optimization algorithms has been extensively developed. However, this approach does not solve the inherent problem of local minima: the determination of the full velocity model, including the low spatial frequency background, requires data with very low frequency (as unrealistic as 1 Hz), coupled with a resolution from low to high frequency to avoid 'cycle-skipping', as for example in [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF][START_REF] Sirgue | Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF] (similar behaviours are studied in electromagnetism in [START_REF] Bao | Inverse medium scattering for the helmholtz equation at fixed frequency[END_REF][START_REF]Numerical solution of an inverse medium scattering problem for Maxwell's equations at fixed frequency[END_REF]). One could also resort to random optimization techniques e.g., [START_REF] James Martin | A stochastic newton mcmc method for large-scale statistical inverse problems with application to seismic inversion[END_REF], which have the ability to find the global minimum even in presence of many parasitic local minima. But they require a very large number of misfit evaluations for the determination of a relatively small number of parameters: they are not adapted for large number of unknowns (in our experiments, we consider several hundred thousands of parameters).

Different misfit functionals have been considered to alleviate the phase shift/cycle skipping problem, such as the logarithmic norm, explored by [START_REF] Tarantola | Inversion of travel times and seismic waveforms[END_REF] and [START_REF] Shin | Waveform inversion using a logarithmic wavefield[END_REF][START_REF] Shin | Comparison of waveform inversion, part 1: conventional wavefield vs logarithmic wavefield[END_REF]. [START_REF] Bozdağ | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF] uses a criterion based on the envelope of the signal; [START_REF] Métivier | Measuring the misfit between seismograms using an optimal transport distance: application to full waveform inversion[END_REF][START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF][START_REF] Yang | Application of optimal transport and the quadratic wasserstein metric to full-waveform inversion[END_REF] are based upon an optimal transport distance. Then, Newton type algorithms provide a natural framework to solve the nonlinear minimization problem in FWI, cf. [START_REF] Pratt | Gauss-Newton and full Newton methods in frequencyspace seismic waveform inversion[END_REF]. However, they require the computation and inversion of the Hessian, which may be overwhelming due to the numerical cost. Thus, some algorithms simply account for the diagonal of the Gauss-Newton Hessian, e.g., [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF]; multi-scale Newton methods are studied in [START_REF] Akcelik | Parallel multiscale gauss-newton-krylov methods for inverse wave propagation[END_REF] for the same reason. Conjugate gradient method with Hessian-vector computation is implemented in [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF]. The estimation of the step length is another challenge for the iterative minimization problem and is often conducted using the line search method (e.g., [START_REF] Nocedal | Numerical Optimization[END_REF]). Alternatives for a more precise step exist, such as the trust region method [START_REF] Eisenstat | Globally convergent inexact newton methods[END_REF], the Maximum Projected Curvature method [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-bystep guide for applications[END_REF]. For Landweber iteration scheme, this step is further given in [START_REF] De Hoop | A convergence analysis of a multi-level projected steepest descent iteration for nonlinear inverse problems in banach spaces subject to stability constraints[END_REF]. The choice of search direction and line search method naturally influences the performance of FWI algorithm, comparisons are given in [START_REF] Barucq | Localization of small obstacles from back-scattered data at limited incident angles with full-waveform inversion[END_REF] in the context of inverse scattering. Yet, none of these local optimization methods can avoid local minima.

In order to overcome the presence of local minima when solving the FWI problem by local optimization algorithms, alternative parametrizations have been proposed early. They are based on a decomposition of the Earth parameter into separate background and reflectivity parameters, and result in an increase of the computational complexity. The Differential Semblance Optimization (DSO) of Symes [START_REF] Symes | Velocity inversion by differential semblance optimization[END_REF] extends the depth reflectivity model to account for the various illuminations in the data, defining a semblance criterion whose minimization produces the desired background model [START_REF] Symes | Velocity inversion by differential semblance optimization[END_REF]; The Migration Based Travel Time (MBTT) approach [START_REF] Chavent | Waveform inversion through MBTT formulation[END_REF][START_REF] Clément | Migration-based traveltime waveform inversion of 2-d simple structures: A synthetic example[END_REF][START_REF] Chavent | Reflection fwi in mbtt formulation[END_REF] parameterizes the depth reflectivity by a data-space reflectivity, which is used to store the reflectivity during the update of the background. Note that the so-called Reflection FWI [START_REF] Xu | Full waveform inversion for reflected seismic data[END_REF] retains from MBTT the migration-demigration process, but differs in that it does not use the data-space reflectivity concept, and uses as search direction the gradient of the classical FWI data misfit.

In this paper, we investigate the optimizability of the FWI problem associated with the time-frequency domain wave equation, and compare it with that of the MBTT reformulation (in terms of background and data-space reflectivity). The objective is to quantify the size of "attraction basins" and "tolerable level errors" which ensure that a local algorithm with an initial guess inside the basin will converge to the global minimum, provided the error on the data is less than the tolerable error. By construction, these quantities depend solely on the (parameter-to-synthetic) forward map, but not on the data, and can be computed without any minimization of the data misfit. Section 2 presents the time-harmonic wave equation to be inverted, together with the FWI minimization problem and its MBTT reformulation employed for the model parameters reconstruction. Section 3 presents the required geometrical tools of [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-bystep guide for applications[END_REF] for optimizability study. Sufficient and necessary optimizability conditions are given, and estimates of the size of the attraction basin in the parameter space, and of the tolerable error level in the data space are derived. Section 4 analyses the influence of (complex) frequency and search geometry on the size of attraction basins and tolerable error levels for the classical FWI formulation. In Section 5, we compare the optimizability of classical FWI and its MBTT reformulation. We expect that the migration-demigration included in MBTT eliminates phase shift/cycle skipping problems, leading to larger attraction basins for the minimum with respect to background, and investigate the effectiveness of this approach. In Section 6 we analyze the optimizability of least squares minimization algorithm for other context (elasticity, boundary conditions). In Appendix A, we review the gradient computation in the frequency domain, emphasizing the specificity of complex valued fields for the adjoint state method. The MBTT decomposition is specified in Appendix B. Appendix C provides numerical experiments of acoustic FWI to validate our estimates.

Time-harmonic inverse wave problem

We formulate the seismic inverse problem associated with the wave equation, for the recovery of acoustic or elastic materials from reflection data. We consider a bounded domain Ω of R 2 which represents the region of interest in the subsurface (note that we restrict ourselves to two dimensions for the numerical experiments but the analysis holds similarly in three dimensions). The boundaries are denoted ∂Ω = Γ 1 ∪ Γ 2 , where we distinguish the upper free surface (Γ 1 ) from the numerical boundary (Γ 2 ). We assume that partial data are acquired on a portion Σ of the domain, which possibly coincides with a part of the boundary. Figure 1

Time-harmonic wave equations

The propagation of waves is defined by the displacement field u, solution to the following partial differential equation,

-

ω 2 ρ(x)u(x) -∇ • σ(x) = g(x), in Ω. ( 1
)
We have introduced the (possibly complex) angular frequency ω, the medium density ρ and the (Cauchy) stress tensor of order two σ. The source is represented by g. We disregard boundary conditions for now. The stress tensor encompasses the different material properties and reduces, for example, when considering isotropy (see [START_REF] Thomsen | Weak elastic anisotropy[END_REF][START_REF] Tsvankin | Anisotropic parameters and p-wave velocity for orthorhombic media[END_REF]).

In the case of acoustic propagation (i.e. propagation in fluid), the wave equation is established for the scalar pressure field p instead. It is based upon the Euler equations (see, for example, [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF][START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF]). If we further assume a homogeneous density, then the wave equation coincides with the Helmholtz equation (see, e.g., [START_REF] Faucher | Contributions to seismic full waveform inversion for time harmonic wave equations: stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms[END_REF]), and the pressure field p is solution to

- ω 2 c 2 (x) p(x) -∆p(x) = h(x), in Ω, (2)
where c(x) is the velocity (wave speed) and h the scalar source. The upper part of the boundary, Γ 1 (see Figure 1), corresponds with a physical free surface, i.e. the interface between the air and the medium, the field satisfies p(x) = 0, on Γ 1 (free surface condition).

(3)

Due to the numerical truncation of the real domain (the Earth), appropriate conditions are imposed on Γ 2 (see Figure 1(b)). One possibility is the consideration of absorbing boundary condition (see [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF]) to ensure that waves that reach Γ 2 are not reflected back to the domain, it is given by

∂ ν p(x) -iωc -1 (x)p(x) = 0, on Γ 2 (absorbing boundary condition), (4)
where ν indicates the normal direction. As an alternative, one can consider instead Perfectly Matched Layers (PML, see [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]) on the sides and bottom boundary (see Figure 1(a)) where the derivatives become

∂ x → 1 + i σ(x) ω -1 ∂ x , in Ω Γx (Perfectly Matched Layer), (5)
and analogously for the other direction, on Ω Γz . In our implementation, the damping function σ is defined following the work of [START_REF] Turkel | Absorbing pml boundary layers for wave-like equations[END_REF][START_REF] Wang | On 3d modeling of seismic wave propagation via a structured parallel multifrontal direct helmholtz solver[END_REF].

The inverse problem aims the recovery of the physical parameter c(x) in (2), using observations of the wave phenomenon p(x) at some location. In the following, we mainly work with this acoustic inverse problem but we also offer some analysis for the isotropic elastic situation in Subsection 6.1 and for other boundary conditions in Subsection 6.2. Note that the methodology developed below is not restricted to inverse wave problem or geophysical setup, and can be applied in any context involving least-squares minimization schemes.

Model and data spaces

We denote by m ∈ M the subsurface model parameter to be recovered, here m = c -2 in (2), and by M ⊂ E the admissible parameter set, which is a subset of the parameter space E. Here, E = R n where n indicates the number of unknowns in the representation of m. We assume that M is a closed convex subset of E equipped with the L 2 norm . E :

m E = n i=1 m 2 i 1/2 , ( 6
)
where m i is the i th component of m.

Remark 1 (Model representation). For the computational resolution of the inverse problem in the discrete settings, it is natural to consider the parameter to be recovered with a piecewise constant representation, such that

m(x) = n i=1 m i χ(D i ). (7)
Here, m i , for i = 1 to n, are scalar coefficients (which represent the model) and χ(D i ) represents the characteristic function on the domain region D i . The subdomains D i (for i = 1 to n) are nonoverlapping and verify

Ω = ∪ n i=1 D i , D i ∩ D j = ∅, ∀i = j. ( 8
)
Inria Therefore the model is represented via n coefficients, and n is the number of unknowns in the FWI inverse problem. Practically, the representation usually corresponds with the space discretization employed for the resolution of the wave equation: one coefficient per cell for discretization based upon mesh (e.g., Finite Element, Discontinuous Galerkin) or one value per node for Finite Differences discretization. Moreover, the piecewise constant model representation [START_REF] Bao | Inverse medium scattering for the helmholtz equation at fixed frequency[END_REF] is also fundamental to show properties of stability, see [START_REF] Alessandrini | Lipschitz stability for the inverse conductivity problem[END_REF][START_REF] Beretta | Inverse boundary value problem for the helmholtz equation: quantitative conditional lipschitz stability estimates[END_REF] and Remark 4.

The data are recorded by the receivers on a portion Σ of the domain Ω, they contain both direct and reflected energy. Because one aims to image the deep structures of the earth, the direct arrivals need to be eliminated from the recorded data before they can be used in the inversion procedure. This is automatically done when a 'linearized' model is used, but in general, for the full Helmholtz problem (2), it requires subtracting the direct arrivals. So we shall denote by d the data deprived from their first arrivals, and by D the associated data space. In the time-harmonic case, d is a vector of q complex numbers and D = C q , with q = n f requency × n source × n receiver . The data space D is equipped with the norm . D ,

d D = q i=1 d i d i 1/2 , ( 9
)
where d i is the i th component of d and denotes the complex conjugate.

Inversion via classical FWI

The essence of FWI ( [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography[END_REF][START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]) is to try to reconstruct the subsurface image by minimizing a misfit functional defined as the difference between the observations and simulations, starting from an initial model. In order to focus on reconstruction of deeper structures, we have to remove the direct arrivals from the forward model. Let us denote by p

s,ω the solutions of ( 2) for m and m s , where m s is a smooth version of m. For the source h at frequency ω, we can define the forward operator at frequency ω for a source h by:

F (h) ω (m) = p (h) ω (x)| Σ -p (h) s,ω (x)| Σ = R p (h) ω (x) -R p (h) s,ω (x) , (10)
where R denotes the restriction operator to the discrete receivers location. When computing derivatives of F, it will be necessary to remember that m s depends also on m. FWI amounts to solve the non linear least squares minimization problem:

min m∈M J (m) = 1 2 ω h F (h) ω (m) -d (h) ω 2 D . (11)
where the data d

(h)
ω are made of the observations at frequency ω for a source h, deprived from first arrivals. For sake of conciseness, we use the vector notation

F(m) = F (h) ω (m), for all h, ω and d = d (h) ω , for all h, ω , (12)
and write the FWI minimization problem as:

min m∈M J (m) = 1 2 F(m) -d 2 D . (13)
As mentioned in the introduction, different norms can be used to build the objective function, and the minimization is usually performed by a Quasi-Newton algorithm, which requires only the gradient of the cost function:

∇J (m) = DF(m) * (F(m) -d), (14)
where DF stands for the Fréchet derivative of F and * is the adjoint. This gradient can be efficiently computed by the adjoint method, which does not require the formation of the Jacobian matrix. Appendix A describes a careful adaptation of this method to the case of complex variables (contrarily to the time-domain formulation, the data and the wavefields are, in the harmonic formulation, complex).

It is well known that the determination of the background (low spatial frequencies) of m by ( 13) is hampered by the many local minima of J , caused by phase shifts in the synthetics (see Figure 18). This can be overcome only if the data contain extremely low frequencies. Appendix C shows a few reconstructions of velocity models via FWI, which illustrate the difficulties inherent to this approach. These difficulties are a motivation for alternative techniques such as the MBTT reformulation of FWI below, and for the optimizability study developed in this paper.

Inversion via MBTT-FWI (background/data-space-reflectiv ity decomposition)

In the MBTT (Migration-Based Traveltime) approach, see [START_REF] Clément | Migration-based traveltime waveform inversion of 2-d simple structures: A synthetic example[END_REF], the Earth model m is parameterized by a smooth background p and a data-space reflectivity s using a migration operator:

m = m(p, s) = p + r = p + W DF(p) * s = p + ω W(ω)DF ω (p) * s(ω), ( 15
)
where r is the depth reflectivity associated to s and p; W is a scaling operator (which possibly depends on the frequency) and * denotes the adjoint. The weight W is meant to compensate for the lower amplitude of deep migrated events, see [START_REF] Plessix | A quantitative kirchhoff migration to esimate the 2d velocity distribution[END_REF]. In our experiments we use a simple scaling proportional to the square root of the depth. where (x, y, z) are the 3D space coordinates with z the depth, and β is a scalar coefficient. We refer to Appendix B for more details regarding the computational aspects of the decomposition. Note that the solution of the linearized version (Born approximation) of the FWI problem ( 13) is of the form [START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF], in which case parameterization [START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF] is not underparameterizing.

Hence, when the full model ( 2) is used, the parameterization by the data space reflectivity s will be able to generate all primary events of the data, but maybe not all multiple events (i.e., it will miss the events associated to multiple reflections involving at least one reflector which generates no primary reflection), cf. [START_REF] Chavent | Data Space Reflectivity and the Migration based Travel Time approach to FWI[END_REF].

When this decomposition is employed, the natural choice for the smooth version m s of m in (10) (to account for the direct arrivals in the forward map F) is simply m s = p. With this change of parameter, the forward map given by [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] [START_REF] Beretta | Lipschitz stability of an inverse boundary value problem for a schrödinger type equation[END_REF] rewrites:

F(p, s) def = F(m)
with F given by [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] and m by [START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF]. ( 16)

By construction, F contains no direct arrivals, which implies that, for a smooth enough background p, F satisfies:

F(p, 0) def = F(p) ≈ 0. (17)
The motivation for this parameterization is to eliminate phase shifts induced in the synthetics by changes in the background p: the events in the synthetic section F(p, s) are obtained from the dataspace reflectivity s by migration followed by simulation with the same kinematic, and hence are expected to have the same phase as those of s, as illustrated in Figure 18, Section 5.1. Besides controlling the phase, this migration-demigration process has the additional property that the stack involved in the migration turns, for a fixed data space reflectivity, the data misfit into a coherency measure for the current background [START_REF] Chavent | Data Space Reflectivity and the Migration based Travel Time approach to FWI[END_REF].

The MBTT-FWI minimization problem is then:

min p∈Ms,s∈D J(p, s) = 1 2 F(p, s) -d 2 D , (18)
where M s ⊂ M is the set of admissible smooth backgrounds and D is the data space. The iterative reconstruction follows alternative updates of the two components of the model:

1. Start the reconstruction with the recovery of the reflectivity part s, from the lowest available frequency.

2. Using the updated reflectivity, which we fix, we can iterate to reconstruct the background part of the model, p.

3. Continue to perform the alternative updates, in a succession of increasing frequency.

This approach has been shown successful in [START_REF] Chavent | Reflection fwi in mbtt formulation[END_REF], for the inversion of synthetic Marmousi data starting from a lower frequency of 5 Hz. Hence another motivation for the optimizability study of this paper is to quantify how far the MBTT reformulation of FWI succeeds in overcoming in step 2 the local minima problem inherent to classical FWI.

Optimizability and Attraction Basins

In this section, optimizability conditions are defined for general least-squares minimization problems, such as [START_REF] Bozdağ | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF] and [START_REF] Chavent | Identification of functional parameters in partial differential equations[END_REF]. By optimizability, we refer to the possibility for a local (deterministic) optimization algorithm to converge to a global minimum, i.e. there is no local minimum to hamper the search of the solution. This analysis follows the work of [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-bystep guide for applications[END_REF]. In addition to the notations of Subsection 2.2, we shall refer to F(M) as the attainable set of the least squares problem.

Definition 1 (Path). A curve P drawn on F(M) ⊂ D is a path of F(M) if it is of the form:

P : t ∈ [0, 1] → F (1 -t)m 1 + tm 2
where m 1 , m 2 are two models of M.

The following set of hypotheses is required to derive the optimizability conditions (cf. [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-bystep guide for applications[END_REF]):

-the forward map F : M → D is twice differentiable along segments of M,

-there exists C ≥ 0 such that ∀m 1 , m 2 ∈ M, ∀t ∈ [0, 1], D t F (1-t)m 1 +tm 2 D ≤ C m 2 -m 1 E
, where D t stands for the derivative with respect to t.

Note that F is indeed twice differentiable because our study is conducted in the finite (discrete) dimensional setting.

Definition 2 (Velocity and acceleration). P is twice differentiable and we denote by V and A (velocity and acceleration along P ) its two first derivatives:

V (t) = P (t), A(t) = P (t). (20)
For simplicity, we shall consider only paths for which V (t) = 0 for all t, so we can define the unit tangent velocity v, and the normal acceleration a by

v(t) = V (t) V (t) D , a(t) = A(t) -A(t), v(t) v(t) V (t) 2 D , (21)
where •, • is the real inner product in D = C q :

z 1 , z 2 = Re{z 1 z2 }. (22)
Let us point out that, due to the limited accuracy of the recording devices, model error and noise, the observed data d do not belong in general to the attainable set F(M). Therefore, it is important that the least squares misfit function does not have parasitic local minima for data d which are "not too far" from the attainable set. This property is made precise by the following definition.

Definition 3 (Optimizability/Attraction Basin). The least squares problem (13) is optimizable on M, or equivalently the parameter set M is an attraction basin for (13), if there exists a neighborhood V of F(M) such that -uniqueness: all data d ∈ V have a unique projection d † on F(M), -unimodality: for any d ∈ V, the distance to d has no parasitic stationary point over F(M), -convergence: if d ∈ V, any minimizing sequence d n ∈ F(M) of the distance to d is a Cauchy sequence for both the norm • D and the arc length distance (P ) along the path P defined by [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-bystep guide for applications[END_REF]. Hence d n converges in F to the unique projection d † of d onto F(M).

Therefore, an optimizable least squares problem (13) exhibits rightful dispositions for its resolution by a local gradient algorithm: uniqueness of the projection and absence of local minimum ensure that the algorithm will converge to a global minimizer whatever the initial guess in its basin of attraction M.

Global Radius of Curvature and Deflection

Following [19, pp. 167-172 and 300-308], we define the global radius of curvature and the deflection along a path P , and further give in Propositions 1 and 2 a characterization and a sufficient condition of optimizability.

Definition 4 (Radius of curvature). The (possibly infinite) radius of curvature R(t) of a path P at t is given by: 1

R(t) = a(t) D = A(t) D V (t) 2 D sin(A(t), V (t)). ( 23
)
The radius of curvature of the whole path P is then defined as:

1 R(P ) def = sup t∈[0,1] 1 R(t) . ( 24
)
It is straightforward to see that

1 R(t) ≤ A(t) D V (t) 2 D , for a.e. t ∈ [0, 1]. ( 25
)
Definition 5 (Global radius of curvature). The (possibly infinite) global radius of curvature R G (t, t ) of a path P at t seen from t , with t = t , is given by:

R G (t, t ) = N + /D if v(t), v(t ) ≥ 0 , N + if v(t), v(t ) ≤ 0 , (26)
where v(t), v(t ) are the normalized velocities along P defined by [START_REF] Chavent | Waveform inversion through MBTT formulation[END_REF], and

N + = max(N, 0) where N = sign(t -t) P (t ) -P (t), v(t ) , D = 1 -v(t), v(t ) 2 1/2 . ( 27
)
The global radius of curvature of the path P is then defined as:

R G (P ) def = inf t,t ∈[0,1] R G (t, t ) ≥ 0. (28)
The interest of global radius of curvature comes from the following proposition.

Proposition 1 (R G > 0 ⇐⇒ optimizability).

The least squares problem (13) is optimizable -or equivalently M is an attraction basin for (13)if and only if it exists R G > 0 such that R G (P ) ≥ R G > 0 for all path P of F(M). The associated neighborhood V is defined by:

V = d ∈ F | dist d, F(M) < R G . (29)
The proofs can be found in [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-bystep guide for applications[END_REF]. The global radius of curvature can be computed numerically using [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF] and ( 27), as will be done in following numerical sections. It can also be estimated via the usual radius of curvature depending on the value of the deflection, which we define now, and which is illustrated Figure 2(a). Definition 6 (Deflection). The deflection between two points t and t of the curve P is the angle between the two velocities V (t) and V (t ) (see Figure 2(a)). It is given by:

Θ(t, t ) = arccos V (t), V (t ) V (t) D V (t ) D ∈ [0, π[ . (30)
The deflection Θ(P ) of the curve P is defined as the largest angle Θ(t, t ) ∈ [0, π] between any two tangent vectors V (t) and V (t ) for any two points t and t of [0, 1]. An infinitesimal variation of the deflection dΘ satisfies

dΘ ≤ A(t) D V (t) D dt. ( 31
)
Denoting t 1 and t 2 the values of t for which the deflection is maximum, the deflection Θ(P ) along the curve P satisfies

Θ(P ) = t2 t1 dΘ ≤ 1 0 A(t) D V (t) D dt. ( 32
)
This majoration is sharp, but it is very conservative: equality holds only when P is an arc of circle with constant velocity V (t) : "when the path P turns always in the same direction with a constant radius". Therefore, the estimation of the deflection Θ(P ) of the path using its upper bound in [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms; numerical results[END_REF] is the worst case estimate.

The relation between global and local radii of curvature is then given by the following proposition.

From Definition 7, a FC/LD problem verifies that, using [START_REF] James Martin | A stochastic newton mcmc method for large-scale statistical inverse problems with application to seismic inversion[END_REF]:

R G (P ) = R(P ) ≥ R > 0 for all paths P , (36)
which shows that FC/LD problems (also referred to as weakly nonlinear inverse problem in [START_REF] Chavent | On weakly nonlinear inverse problems[END_REF]) are necessarily optimizable. Notice that Proposition 1 gives a characterization of optimizable problems, whereas Definition 7 provides only a sufficient condition.

Directional Attraction Basins

Numerical application of previous section to check whether or not a given least squares problem is optimizable becomes quickly intractable when the number of parameters n increases, as it is the case in seismic inversion. So we limit ourselves to directional (or one-dimensional ) parameter sets of the form:

M(m 0 , u, ∆) = [m 0 -∆u, m 0 + ∆u] with u E = 1. (37)
Here, m 0 is a nominal model, u a normalized perturbation direction, and ∆ gives the size of the domain of investigation. The associated attainable set is constructed from the path P defined by:

P : t ∈ [0, 1] ; F m 0 + (2t -1)∆u . ( 38
)
We refer to directional optimizability when the problem is optimizable for an interval such as [START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF], and this interval is a directional attraction basin. Directional optimizability is only a necessary condition for optimizability, but it will allow to analyze the behavior of seismic inverse problems and to compare formulations: the size ∆ of a directional attraction basin in a descent direction tells us how far one can move away in this direction without being stopped by parasitic local minima. Our objective now is to determine (see illustration Figure 2):

1. the size ∆ u m0 of the directional attraction basin centered at m 0 . The larger ∆ u m0 is, the better the least squares problem is amenable to minimization by local algorithm, because we allow a larger area for investigation. In our numerical experiments, we shall scale the estimate with the norm of the nominal model, m 0 E , to provide relative (to the model) quantity.

2. the associated tolerable error level R u G,m0 . It is the largest tolerable error on the data d which ensures the absence of parasitic local minima for the least squares objective function

t ∈ [0, 1] ; 1 2 F(m 0 + (2t -1)∆ u m0 u) -d 2 D ( 39
)
over [0, 1]. The larger R u G,m0 is, the better is the robustness of the minimization procedure to noise in the data. In our numerical experiments, we divide the estimates with the norm of the synthetic data d 0 = F(m 0 ) to provide relative (to the data) quantity.

t t V (t) t V (t ) Θ(t, t ) F(m0) F(m0 -∆u) F(m0 + ∆u) (a) F(m0) F(m0 -∆u) F(m0 + ∆u) d dist(d, F) (b)
Figure 2: A one-dimensional setup for least squares problems. The figure lives in the data space, the attainable set is the path P image of the interval M(m0, u, ∆) of the model space. (a) Illustration of the computation of the deflection Θ(t, t ) between two arbitrary points t and t . (b) The path has a finite curvature and the deflection is smaller than π/2, so the FC/LD Property 7 is satisfied, and RG = R > 0 by Proposition 2. Hence the "distance to d" function cannot have local minimum over P provided the data d is at a distance of the attainable set P = F(M(m0, u, ∆)) smaller than R.

We shall use two types of estimate:

• Θ-estimates of ∆ u m0 , where the optimizability over M is obtained by satisfying the sufficient condition Θ(P ) ≤ π/2 of Definition 7. In this case, R G (P ) = R(P ), so the tolerable error level is given by the minimum over the [0, 1] interval of R(t) given by [START_REF] Chavent | On weakly nonlinear inverse problems[END_REF].

• R G -estimates of ∆ u m0 , where optimizability is obtained by satisfying the R G > 0 characterization of optimizability of Proposition 1. In this case, the associated tolerable error level R G has to be computed by evaluating numerically the infimum in (28) using [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF] and [START_REF] De Hoop | A convergence analysis of a multi-level projected steepest descent iteration for nonlinear inverse problems in banach spaces subject to stability constraints[END_REF].

Remark 2. The size of attraction basins depends only on the choice of the forward map to be inverted, but not on the optimization algorithm used (e.g. Newton method, gradient descent, etc). The choice of the method naturally affects the rate of convergence and the speed at which the final solution is eventually reached, but it has no influence on the presence or absence of local minimum.

Local

Θ-estimate for ∆ u m 0 and associated tolerable error level R u m 0
We provide here a local Θ-estimate ∆ of the attraction basin, in the sense that it is based only on the velocity V and acceleration A at m 0 in the direction u. In order to ensure that the deflection of the path P defined by [START_REF] Li | Optimization in the complex domain for nonlinear adaptive filtering[END_REF] is smaller than π/2, we use the upper bound (32) of Θ(P ), which, according to the optimizability condition of definition 7, ensures in turn that M(m 0 , u, ∆) is an attraction basin.

With the notations of [START_REF] Cartan | Differential calculus[END_REF] for the directional derivative, the chain rule differentiation gives:

       V (t) = ∂F ∂m ∂m ∂t = 2∆DF(m)(u), A(t) = ∂ 2 F ∂m 2 ∂m ∂t 2 + ∂F ∂m ∂ 2 m ∂t 2 = 4∆ 2 D 2 F(m)(u, u), (40)
where u acts as the direction of derivation. Then we use a rectangle approximation in equation [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms; numerical results[END_REF], which gives the approximate upper bound Θ u m0 to the deflection Θ(P ):

Θ(P ) ≤ 1 0 A(t) D V (t) D dt ∼ A(1/2) D V (1/2) D = 2∆ D 2 F(m 0 )(u, u) D DF(m 0 )(u) D def = Θ u m0 . ( 41
)

Inria

This gives immediately a local Θ-estimate of the size ∆ of an attraction basin at m 0 in the direction u:

∆ u m0 = π 4 DF(m 0 )(u) D D 2 F(m 0 )(u, u) D , local Θ-estimate. ( 42
)
This estimate is an approximate (because of the rectangle approximation of the integral) lower bound (because it is based on the upper bound ( 32)) to the size of the largest attraction basins at m 0 in the direction u. It is computationally cheap, and we will use it in Section 4 to analyze the effect of frequency and perturbation direction on the optimizability of the FWI least squares problem.

The associated tolerable error level R u m0 is then the minimum of the radius of curvature along P = F(M(m 0 , u, ∆)), which is approximated simply by its value at m 0 , that is R(t = 1/2) given by ( 23):

R u m0 = V (t) 2 D A(t) D | sin A(t), V (t) | t= 1 2 , tolerable error level, (43)
where V (t) and A(t) have been defined in [START_REF] Mandache | Exponential instability in an inverse problem for the Schrödinger equation[END_REF].

Remark 3. One may use the directional attraction basins M(m 0 , u, ∆ u m0 ) to construct an approximate multidimensional attraction basin M(m 0 , ∆ m0 ) as follows: suppose that ∆ u m0 ≥ ∆ m0 ≥ 0 for all u for some ∆ m0 , and define:

M(m 0 , ∆ m0 ) = B(m 0 , ∆ m0 ), (44)
where B denotes a ball in the model space. The size ∆ m0 of an approximate multidimensional attraction basin M(m 0 , ∆ m0 ) follows immediately from the previous results by computing an upper bound Θ m0 to Θ u m0 using (41):

Θ u m0 = 2∆ D 2 F(m 0 )(u, u) DF(m 0 )(u) ≤ 2∆ D 2 F(m 0 ) u 2 λ DF min u = 2∆ D 2 F(m 0 ) λ DF min def = Θ m0 , (45)
where λ DF min denotes the lowest singular value of DF(m 0 ), and we use the Frobenius norm for the matrix norm of the numerator. Estimation of ∆ m0 follows:

∆ m0 = π 4 λ DF min D 2 F(m 0 ) . ( 46
)
Note that M(m 0 , ∆ m0 ) is only an approximate attraction basin, because condition R G > 0 (Proposition 1) or Θ ≤ π/2 (Definition 7) holds only for segments passing through its center m 0 . Then, (43) provides a lower bound R m0 independent of u to the radius of curvature of M(m 0 , ∆):

R m0 = λ DF (m0) min 2 D 2 F(m 0 ) . ( 47
)
3.4 Exact Θ-and R G -estimates of ∆ u m 0 and associated tolerable error levels

The determination of the exact Θand R G -estimates of the attraction basin centered at m 0 in a direction u inside an interval M(m 0 , u, ∆) of given size ∆ requires the numerical computation of the deflection Θ(t, t ) and the global radius of curvature R G (t, t ) between any two points

F(m 0 + tu) , -∆ ≤ t ≤ ∆ and F(m 0 + t u) , -∆ ≤ t ≤ ∆ (48)
of the path P , which is the image by F of the investigated interval M(m 0 , u, ∆).

For this purpose,deflection maps and global radius maps are computed, and they display the values of Θ(t, t ) (Definition 6) and of R G (t, t ) (Definition 5) between the points of M(m 0 , u, ∆). On the diagonal of the maps, where t = t , R G is not defined by [START_REF] Crase | Robust elastic nonlinear waveform inversion: Application to real data[END_REF] [START_REF] De Hoop | A convergence analysis of a multi-level projected steepest descent iteration for nonlinear inverse problems in banach spaces subject to stability constraints[END_REF], and we indicate instead the values of R(t) given by (21) [START_REF] Chavent | On weakly nonlinear inverse problems[END_REF], which represent the limits of R G (t, t ) when t → t. One can then read on these maps the exact Θ-estimate (resp. the exact R G -estimate) of ∆ u m0 by searching for the largest squares centered at (0,0) where

Θ(t, t ) ≤ π/2 (resp. where R G (t, t ) > 0).
During this process, when ∆ increases from 0 to its exact Θ-estimate, the associated exact tolerable error R u m0 = inf -∆≤t≤∆ R(t) decreases from its value R 0 at m 0 to the tolerable error R u m0 of the Θattraction basin. When ∆ increases further to its exact R G -estimate, the tolerable error is R G = inf -∆≤t,t ≤∆ R G (t, t ), which continues to decrease, until it reaches the value 0 of the R G -attraction basin.

These computationally intensive exact estimates will be used in Section 4 on the analysis of FWI, for comparison purpose with the local Θ-estimates. They will also be used in Section 5 to compare more precisely the optimizability properties of classical FWI and MBTT-FWI.

FWI: Numerical determination of convergence basins and tolerable error levels

Our objective in this section is to give a comprehensive understanding of the behavior of the seismic FWI problem [START_REF] Bozdağ | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF] with respect to parameters such as the frequency, the geometry of the target, the parametrization, etc. Therefore, we select a nominal model m 0 and compute, according to Subsection 3.2, the size ∆ u m0 of an attraction basin and the associated tolerable error level R u m0 for different (normalized) directions u. When m 0 is smooth, the direction u can be seen as the search direction, or the geometry of the unknown parameter to be reconstructed in the framework of inversion.

Influence of frequency ω and search direction u

We consider the two-dimensional geophysical setup of Figure 1 for the acoustic Helmholtz equation ( 2), and select the squared slowness m := c -2 as unknown parameter, which is expressed in s 2 m -2 . The forward map F ω at frequency ω for parameter m is a vector of C (nrcv×nsrc) with n rcv the number of receivers per source and n src the number of sources. In this experiment, we take 19 sources and 183 receivers associated with each source. Both are located near the surface (accordingly to Figure 1), note that the receivers remain in the same position for all sources. Thus, we work with reflection data obtained from a one side (the surface) illumination.

The nominal model m 0 is the smooth 9.2 × 3 km background velocity pictured in Figure 3, where the velocity varies linearly with depth, according to the fact that no information is known on the structures to reconstruct.

We investigate the behavior of attraction basins for different directions u, which are chosen to be based on straight reflectors, or extracted from the Marmousi model (probably the most popular velocity model in geophysics, synthetically designed by the Institut Français du Pétrole in 1988). They are shown in Figure 4: u 1 has a single, slim reflector (Figure 4(a)), u 2 has two slim reflectors (Figure 4(b)), u 3 has a single large reflector (Figure 4(c)) and u m is obtained by high-pass filtering from the Marmousi model of Figure 4(d). The amplitude for the directions are selected such that u k is normalized:

u k = 1, for all k.
Given a perturbation size ∆, one can think of m 1 = m 0 +∆u as being the target velocity, which one tries to retrieve starting from the smooth model m 0 by minimization of the FWI objective function in direction u by a local algorithm (it is analogous to the update formula in Newton-type algorithm). The question is: how large can ∆ be while still guaranteeing that the algorithm will not stop before the global minimum, provided that the error in the data is smaller than some R or R G ?

Local Θ-estimates of ∆ u m0 and R u m0
The determination of ∆ u m0 with the local Θ-estimate, [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF], and R u m0 with [START_REF] Mora | Nonlinear two-dimensional elastic inversion of multioffset seismic data[END_REF], requires the computation of the first and second order directional derivatives in the direction u, and we refer to Appendix A.3 for the computational details. Figure 5 shows the evolution of ∆ u m0 (top) and R u m0 (bottom) for frequencies between 0.5 to 15 Hz, for the four directions u of Figure 4.

We observe the following.

-In all cases, the size ∆ u m0 of the directional attraction basin decreases with increasing frequency, meaning that the largest basin of attraction is obtained for the lowest frequencies. This is the expected behaviour that indicates that lower frequencies give a better chance for the iterative minimization to converge, especially when no prior information is known for the initial model.

-Concerning the effect of the perturbation direction on attraction basins, Figure 5(a) shows that their size decreases in the order u 1 → u 2 → u m → u 3 . This can be understood as follows: Directions u 1 , u 2 and u m are essentially reflective, and contain no low frequency components. Hence increasing ∆ will not change the phase of backscattered data, but because of multiple reflections some nonlinearity will occur, which limits the size of the attraction basins. The directions u 1 and u 2 with the less multiple reflections have the largest basins of attraction. Then comes u m which has many reflectors and hence many multiples. On the contrary, direction u 3 has a strong low frequency component besides its reflectivity. Hence increasing ∆ will generate a phase shift in the backscattered data, which is likely to create local minima more quickly, and one expects the size ∆ u m0 of the attraction basin in the direction u 3 to be much smaller than in the three other directions.

-Concerning the largest tolerable error, Figure 5(b) shows that R u m0 increases with the frequency, meaning that the low frequencies are more likely to be affected in the presence of noise. We notice the exception of u 3 for which the allowed error remains relatively similar (and worse than the others) with frequency.

These results clearly show that it is harder to converge towards salt dome than it is to recover Frequency (Hz) Marmousi-like structures. Because of its low spatial frequency component, the former reduces the size of the attraction basin and is less tolerant to noise in the data. This situation is further illustrated in the context of model reconstruction in Appendix C. Note that the order of magnitude in Figure 5 are around 10 -1 of the norms of the model (top) and data (bottom).

∆ u m0 / m 0 u 1 u 2 u 3 u m (a) Evolution of ∆ u m 0 computed with local Θ-estimate.
Exact Θand R G -estimates of ∆ u m0 and R u m0

As an alternative to confirm these preliminary observations, we compute ∆ u m0 using the exact Θestimate and R G -estimate. It requires the computation of deflection and global radius maps on the investigated interval, from which the estimate is extracted, cf. Subsection 3.4. In Figures 6 and7, we show the resulting maps in the direction of the Marmousi reflectors u m of Figure 4(d) and of the single large reflector u 3 of Figure 4(c), using frequencies 4 and 7 Hz. In Table 1, we show the numerical values of the interval size centered around m 0 , depending on the estimation.

The maps confirm the results of the previous experiment: the exact Θ-estimates of ∆ u m0 (top of the figures) are of the same order of magnitude, but slightly larger than the local Θ-estimates of Figure 5, cf. Table 1. Similarly, the exact R G -estimates of ∆ u m0 ( bottom of the figures), are larger (Table 1), but at the price of a smaller tolerable error on the data. We also note a much more consistent pattern in the deflection for u m , compared to u 3 . Comparing the formulas, we observe that the inexpensive local Θ-estimates already provide accurate values for classical FWI problems, close to the more refined (and computationally more intensive) exact Θ-and R G -estimates. Therefore, we should only use the local Θ-estimate in the following of this section.

direction u m direction u 3 ∆ um m 0 m0 R um G,m 0 F (m0) ∆ u 3 m 0 m0 R u 3 G,m 0 F (m0)

Influence of the initial model

We now investigate the influence of the initial model m 0 on the size of attraction basins. We consider two initial models: the acoustic Marmousi medium m m and a model m c encompassing objects of high contrasts, given in Figure 8. The local Θ-estimate of ∆ u mj given by ( 42) is plotted in Figure 9, for the two directions u 1 (one reflector) and u 3 (one large reflector). We see that ∆ u mj is barely affected by the choice of initial model for the computation. We still observe the decrease in the size of the interval with increasing frequency, and the deterioration in the direction of the large reflector u 3 . Clearly, when comparing amplitudes of our estimates, the direction u has much more influence than

Approximate n-dimensional attraction basins

We compare the directional local Θ-estimates of Subsection 4.1, evaluated at the smooth model m 0 of Figure 3 in the Marmousi direction u m of Figure 4(d), with the approximate multidimensional attraction basins centered in m 0 introduced in Remark 3. The determination of the radius ∆ m0 of such a spherical attraction basin by [START_REF] Plessix | A quantitative kirchhoff migration to esimate the 2d velocity distribution[END_REF] and of the associated level of tolerable error R m0 by [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF] requires the computation of the lowest singular value of the Jacobian of F, and the computation of its full second order derivative matrix. In order to reduce the computational burden associated with this computation, we develop a simple model reduction process where the wave speed in adjacent cells is averaged to generate a coarse representation. Note that this approach is also employed to improve the stability of the inverse problem, cf. [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation: quantitative conditional lipschitz stability estimates[END_REF]. When applied to the smooth acoustic model of Figure 10(b) shows an important reduction of the size ∆ m0 of the n-dimensional basin compared to the directional estimates in the Marmousi direction u m using the local Θ-estimate (six orders of magnitude). Yet, it confirms the main property, which is the decrease of the attraction basin size with increasing frequencies. A similar observation holds for the tolerable data error level R m0 , see Figure 10(c), which suffers an important reduction compared to the directional estimates. In addition, we observe a different pattern between R m0 , which decreases with frequency, and the directional estimates R u m0 , which increases with frequency. As illustrated in Figure 5, R u m0 depends more on the direction than ∆ u m0 , and it is reflected on the n-dimensional estimates, which takes the worst case. ), the black squares ( ) show the directional estimates using um (reflectors from the Marmousi model, see Figure 4(d)).

From this first set of experiments, we clearly see that the frequency should evolve from low to high values in order to facilitate the convergence. When no prior information is known for the initial model, one should take advantage of the larger interval size ∆ u m0 given by low frequencies. Here we motivate the well-known frequency progression (e.g., [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF][START_REF] Sirgue | Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF]) from quantitative estimates of the size of the basin of attraction, which allows the quantification of the benefits. However, the low frequencies require more accuracy regarding the data, as illustrated with the estimates of R u m0 . In addition, we mention the fact that seismic data are particularly affected by noise at low frequencies, which are sometimes unusable. This highlights where resides the difficulty: the low frequencies that are required to converge may be unusable due to the noise, or the inaccuracy may violate the curvature condition.

Frequency bandwidth data

We have seen that the sequential progression in frequency must be conducted from low to high regime. We now investigate the case of frequency bandwidth: when the forward problem encompasses not only one, but several frequencies. This approach is particularly natural for the time domain inverse problem where several frequencies are simultaneously contained in the data. The forward problem is now a vector of C (nrcv×nsrc×nω) , where n ω is the number of frequencies taken in the subgroup.

We remain with the smooth background of Figure 3 for m 0 . We select groups of ten consecutive frequencies. The initial group encompasses the frequencies from 0.1 to 1 Hz using 0.1 Hz increment. Similarly, we design fifteen groups so that the last one is from 14.1 to 15 Hz, employing identical increment. The local Θ-estimate of the attraction basin size is shown Figure 11, together with the largest tolerable error R u m0 . We investigate two directions: the Marmousi one of Figure 4(d) and the large reflector of Figure 4(c). We picture curves and histograms, that correspond to sequential or group of frequencies respectively. Numerical values for the attraction basis are prescribed Table 2. 2 highlight that -the size of the attraction basin when using group of frequencies is dictated by the largest frequency in the group. Therefore, it is better to use the lowest frequency sequentially in order to have the largest size possible. For instance, using the single frequency 0.1 Hz provides a size of interval one order of magnitude larger than when including frequency up to 1 Hz.

-At a fixed frequency, taking into account lower frequencies does not affect the size of the interval.

This preference on sequential over multiple frequency is actually advocated in [START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF], motivated, as the frequency progression in [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF], by the intuition that it reduces the cycle-skipping effect. Here m 0 , computed with local Θ-estimate, with frequency depending on the use of sequential or group of frequencies. Here u is chosen to be the Marmousi direction, u = um, according to the Figure 4(d).

we justify this choice by a quantitative estimation of the size of the attraction basin. Moreover, we show that once the lowest frequencies are taken individually, incorporating them in higher frequency content does not modify the size of the attraction basin.

Regarding the largest tolerable error on data, Figure 11(b) shows that adding frequency increases the robustness in the direction u m , except for low frequency. Then, for the direction of the large contrasting object, u 3 , we do not observe any improvement, and the allowed distance remains relatively stable, as it was observed in Figure 5.

Selection of complex frequencies

The use of complex frequencies (also referred to as the "Laplace-Fourier" domain) has shown advantageous behavior when used during reconstruction algorithm, as mentioned in the introduction. In this case, the angular frequency ω in (2) writes as

-ω 2 = (σ + 2iπf ) 2 , iω = σ + 2iπf, ( 49
)
where f is the 'Fourier' component of the frequency (in Hz) and σ a damping factor, representative of the Laplace component. Here, we are interested in deciding how to select the progression of those coefficients to obtain the best convergence properties. Applications have been carried out with a progression that appears mostly intuitive in literature. In the case of zero-Fourier frequency (f = 0), it is proposed from low to high damping coefficients, [START_REF] Shin | Waveform inversion in the laplace domain[END_REF]. However, for the complex frequency progression (f = 0 and σ = 0), the evolution of the damping parameter at fixed f is proposed to be from high to low damping, see [START_REF] Shin | Sequentially ordered single-frequency 2-d acoustic waveform inversion in the laplace-fourier domain[END_REF][START_REF] Petrov | Three-dimensional inverse modelling of damped elastic wave propagation in the fourier domain[END_REF]. We compute the estimation of ∆ u m0 (using local Θ-estimate) and R u m0 with complex frequencies to see how it affects the basin of attraction and robustness to noise. We consider m 0 the smooth background velocity and the similar setup as employed in Subsection 4.1. The estimates are shown in Figures 12, 13 and14, where we compare the behaviours on a logarithmic scale for the twodimensional complex frequency plane. Here, both f and σ varies from 0 to 15 (excluding the case where both f and σ are zero).

It is clear that the incorporation of a damping coefficient in the frequency is beneficial for the convergence, because it increases the size of the attraction basin of several orders of magnitude. Then, the process of frequency selection, for convergence, is the following:

-for a fixed f , the interval size increases with increasing damping, with the exception of the case f = 0, cf. Figures 12(b) and 13(b).

-At fixed damping σ, the interval size decreases with increasing f .

-Following a large to small interval size gives, 1. initial iterations should use f = 0 (largest interval), and a damping coefficient progression from low to high (see Figures 12(b) and 13(b)).

2. Fix f = 0 as small as possible, the damping coefficient should evolve from high to low (see Figures 12(a) and 13(a)); then f increases progressively.

3. Finally, one can take σ = 0 and f evolves from low to high. Regarding the largest tolerable data error R u m0 , we see that -at fixed σ, increasing f allows an increase in the distance between the data and the attainable set.

-The pattern depends strongly on the type of the search direction geometry, as it has already been observed in the previous experiment (and contrary to ∆ u m0 ). The distance is allowed to increase with σ for the direction u m , but decreases in the direction of the strong reflection, u 3 . Also in the direction of the Marmousi reflectors, Figure 14(a), we observe a band where the allowed distance decreases, around σ = 1. On the other hand, for both directions, the case where σ = 0 appears the most robust (i.e. where the allowed distance between the data and the attainable set is maximal). -Regarding the magnitude, we clearly see that the high contrast object has globally much less allowance, except when the damping is close to zero.

It is important to notice that the use of complex frequency, especially when f = 0, improves the convergence by increasing notably the size of the attraction basin. Therefore, it provides a fundamental benefit for the reconstruction, in particular when no information is initially known. Note that it basically acts as a very low frequency (note that the magnitude of the estimates for complex frequencies eventually matches the very low frequency ones; we can see a continuation in the magnitudes). On the other hand, such frequencies are more sensitive to noise, and require accuracy in the data. We refer to [START_REF] Shin | Waveform inversion in the laplace domain[END_REF][START_REF] Shin | Sequentially ordered single-frequency 2-d acoustic waveform inversion in the laplace-fourier domain[END_REF][START_REF] Petrov | Three-dimensional inverse modelling of damped elastic wave propagation in the fourier domain[END_REF][START_REF] Faucher | Contributions to seismic full waveform inversion for time harmonic wave equations: stability estimates, convergence analysis, numerical experiments involving large scale optimization algorithms[END_REF] for illustrations of reconstruction.

Remark 4 (Stability and frequency progression).

We have shown that the reconstruction procedure requires low frequencies to ensure a larger size for the basin of attraction, hence allowing the better chance of convergence despite the possible lack of prior information. However, higher frequencies are still required for reasons of stability, and justify the traditional progression in frequency along with the minimization iterations. Stability of inverse problem has been studied in particular with the Calderón problem, where the data are the Dirichlet-to-Neumann map. Assumptions are commonly introduced to obtain a conditional Lipschitz-type stability results, see, e.g., [START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF][START_REF] Mandache | Exponential instability in an inverse problem for the Schrödinger equation[END_REF][START_REF] Alessandrini | Lipschitz stability for the inverse conductivity problem[END_REF][START_REF] Beretta | Lipschitz stability of an inverse boundary value problem for a schrödinger type equation[END_REF] and the references therein. For instance, the piecewise constant model representation introduced Remark 1 is a necessary ingredient, e.g., [START_REF] Alessandrini | Lipschitz stability for the inverse conductivity problem[END_REF][START_REF] Beretta | Inverse boundary value problem for the helmholtz equation: quantitative conditional lipschitz stability estimates[END_REF]. Lipschitz-type stability is obtained when, for two models m 0 and m † , it holds that

m 0 -m † ≤ C F(m 0 ) -F(m † ) , (50)
where C represents the stability constant. The iterative procedure for the reconstruction minimizes the difference between the forward map and the stability result ensures the accuracy of the reconstruction (depending on the values of C).

In [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation: quantitative conditional lipschitz stability estimates[END_REF], numerical estimates of the stability constant are obtained in the context of the geophysical inverse problem. It illustrates, for example, the influence of the number of coefficients n (in the model representation [START_REF] Bao | Inverse medium scattering for the helmholtz equation at fixed frequency[END_REF]) for stability. When the stability constant C is small, it indicates that the correspondence in the data guarantees the accuracy of the reconstruction. Here, we estimate the stability constant (following the methodology of [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation: quantitative conditional lipschitz stability estimates[END_REF]) to see its dependency with complex frequency, in order to compare with the convergence estimates of Figures 12 and13. We use m 0 the smooth background of Figure 3 and select m † to be the Marmousi model of Figure 8 -Incorporating damping with frequency deteriorates the stability constant (the numerical estimates increases). More precisely, low frequencies and damping (f and σ respectively) give the worst result in terms of stability. Therefore, the estimates for convergence (Figures 12 and13) and stability (Figure 15) behave oppositely. The convergence improves at low frequency and with damping meanwhile the stability deteriorates. Higher frequency reduces the convergence radius but ameliorates the stability (by reducing the wavelength, such frequencies allow to capture more details). Both aspects of the problem have to be taken into account during the reconstruction; it designs the well-known multi-frequency algorithm for the computation, starting from low frequency (to foster the convergence) towards high frequency (to foster stability) (see, e.g., Algorithm 1).

Optimizability: MBTT versus classical FWI

In this section, we apply the optimizability analysis tools of Section 3 to compare the MBTT formulation of FWI [START_REF] Chavent | Identification of functional parameters in partial differential equations[END_REF] (where the unknown model is parametrized by a smooth background p and a data space reflectivity s), with the original FWI problem [START_REF] Bozdağ | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF] (where the unknown model is the squared slowness m). The MBTT approach was originally designed to overcome the phase shift problem which hampers classical FWI, so the optimizability study of this section is expected to quantify the gainif any -of the MBTT formulation over the classical FWI in this respect.

Choice of the model

For a fair comparison, the two approaches have to be applied to the same model, so we construct a nominal model m 0 whose MBTT decomposition, p 0 , s 0 , is known exactly, i.e. which satisfies m 0 = m(p 0 , s 0 ) according to [START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF]. We have decided to construct a model m 0 inspired by the Marmousi model m m of Figure 8(a). This is achieved by defining p 0 and s 0 by:

p 0 = m 0 , ∀ω : s 0 (ω) = F ω (m m ) -F ω (m 0 ), (51)
Inria where m m is the Marmousi model of Figure 8(a) and m 0 is the smooth background of Figure 3. Then, m 0 is defined by:

m 0 def = p 0 + ω W(ω)DF * ω (p 0 )s 0 (ω) r0(ω) = p 0 + WDF * (p 0 )s 0 r0 = m(p 0 , s 0 ), (52)
where r 0 (ω) (respectively r 0 ) is the depth reflectivity associated to frequency ω1 (respectively to the sum of all frequencies). Definition [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF] guaranties that:

F(m 0 ) = F(p 0 , s 0 ). (53)
We choose the weight W proportional to the square root of depth, as proposed in Section 2, with a proportionality coefficient such that the model reflectivity level β at frequency ω defined by

β(ω) def = r 0 (ω) / p 0 ( 54
)
takes the value 1% at all frequencies. When β is small, the forward modeling part of the MBTT approach becomes close to the Born approximation.

In Figure 16, we illustrate the resulting models r 0 (ω) for three frequencies: 2, 4 and 7 Hz. We also show the model r 0 where the frequency sum contains frequencies between 0.5 to 15 Hz, with 0.5 Hz increment. We observe that the reflectivity, defined from the difference between observations and simulations using a smooth background, provides structures of size consistent with the selected frequency. For the global model, shown Figure 16(d), we see the contributions of all wavelengths, and we can distinguish some structures of the Marmousi medium given Figure 8(a). For simplicity, in the following, we restrict ourselves by studying only single frequency nominal models, which means that we only work with models resulting from a single, fixed ω:

m 0 (ω) = p 0 + W(ω)DF * ω (p 0 )s 0 (ω) = p 0 + r 0 (ω). (55)
It allows us to study the behaviour of both approaches (global velocity model reconstruction or MBTT decomposition) with individual frequency, preventing artifacts that may appear due to the frequency summation in the reflectivity.

Choice of perturbation directions

Background perturbation We consider a background perturbation in the direction of the onedimensional ramp u of Figure 17, and illustrate the effect of such a perturbation applied either to the global model m, as in the case of FWI, or to the background unknown p in the case of MBTT. Receiver index

F(m 0 ) F(m 0 + T u) F(p 0 + T u, s 0 )

Reflectivity perturbation

The FWI objective function is known to be nearly quadratical with respect to reflectivity, i.e. to the high spatial frequency part of m, and the same property holds by construction for the dependance of the MBTT objective function with respect to the data space reflectivity s. Hence one expects large basins of attraction with respect to s in the MBTT formulation. We choose for u s a random vector of the data space.

Inria

Comparison of Local Θ-estimates

The formulas for local Θ-estimates have already been derived in Section 3.3 for classical FWI. Application of these formulas to F(p, s) instead of F(m) gives immediately the local Θ-estimates of the sizes ∆ u p0 and ∆ us s0 of attraction basins for p and s at p 0 , s 0 in directions u, u s for the MBTT formulation of FWI:

       ∆ u p0 = π 4 D p F(p 0 , s 0 )(u) D 2 p F(p 0 , s 0 )(u, u) , ∆ us s0 = π 4 D s F(p 0 , s 0 )(u s ) D 2 s F(p 0 , s 0 )(u s , u s ) . ( 56
)
We have omitted the frequency index for clarity in the notation. Note that this is not misleading because we have chosen to work with single frequency content here, i.e., the reflectivity r only contains one selected frequency.

Figure 19 shows the evolution of the local Θ-estimates with frequency. The size of the corresponding attraction basins decreases with frequency, when the perturbation is applied on m, p and s. The left part of the figure shows a slightly larger attraction basin in the direction of the background perturbation u when it is applied to the propagator part p 0 of the MBTT parameterization rather than when applied directly to m 0 . But its not up to the expectations raised by the claim that the MBTT parameterization allows to overcome the phase shift problem [START_REF] Chavent | Data Space Reflectivity and the Migration based Travel Time approach to FWI[END_REF][START_REF] Chavent | Waveform inversion through MBTT formulation[END_REF]. Regarding s, the estimated size appears surprisingly small compared to the large attraction basin expected, see Section 5.2. But one has to remember that these are local Θ-estimate, which can be very pessimistic, as explained in Section 3.3, and postpone more definitive comments to the end of Subsection 5.4, where exact Θ-estimates are calculated. 42) and ( 56). Here p0 is the smooth velocity background of Figure 3, the direction u for p is given Figure 17, and the direction us for s is a random vector.

In the MBTT representation, the reflectivity uses only the selected frequency.

Comparison of exact Θ-and R G -estimates

We apply the method described in Section 3.4 for the case of classical FWI, which translates immediately to the case of MBTT by replacing the FWI forward map m ; F(m) by the MBTT forward map p, s ; F(p, s). This leads to the computation of deflection and global radius of curvature maps between the following points: FWI (attraction basin for m) : F(m 0 + tu) and F(m 0 + t u);

MBTT (attraction basin for p) : F(p 0 + tu, s 0 ) and F(p 0 + t u, s 0 ); MBTT (attraction basin for s) : F(p 0 , s 0 + tu s ) and F(p 0 , s 0 + t u s );

for all -∆ ≤ t, t ≤ ∆.

We first compute the deflection and global radius of curvature maps for m and p, using values of t and t in an interval [-∆, ∆] which is chosen to represent in each case about ±20% of the norm of m 0 or p 0 defined in [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF]. Figures 20 and21 show these maps at two selected frequencies: 4 and 7 Hz, and Table 3 summarizes the extracted exact estimates of ∆ u m0 , ∆ u p0 , together with the local estimates extracted from Figure 19(a). model m model p model s 3). Notice that the size of the Θ-attraction basin is divided by two for both FWI and MBTT when the frequency increases from 4 to 7 Hz.

∆ u m 0 m0 R u G,m 0 F (m0) ∆ u p 0 p0 R u G,p 0 F (m0) ∆ u s 0 s0 R u G,s 0 F (m0)
-The second observation concerns the strict positivity of the global radius of curvature R G , bottom of Figures 20 and21, which determines the R G -attraction basin characterized by a zero tolerable error (see Section 3.4). For the MBTT formulation, Figures 20(b) and 21(b), R G remains strictly positive all over the map, which shows that the R G -basin is larger than the investigated interval. Its size is of approximately 20% at both 4 and 7 Hz (Table 3). On the contrary, for the usual FWI formulation, Figures 20(a) and 21(a), R G decreases very rapidly to zero when one moves away of the diagonal, producing smaller R G -attraction basins, with size of 5% at 4 Hz and 2.5% at 7 Hz, smaller by a factor four to eight to the corresponding MBTT attraction basins.

-Concerning the magnitude of R G , whose minimum over the attraction basin gives the tolerable error level, one sees that it takes larger values for FWI near the main diagonal, i.e. for small attraction basins, than for MBTT over the whole map, i.e. for large attraction basins, which is confirmed by the values of R G in Table 3, which give the tolerable error level associated with the Θ-attraction basins (this level is zero by definition for the R G attraction basins).

To summarize, MBTT extends significantly the size of attraction basins with respect to background perturbations, at the price of a reduction in the admissible error level. This explains the success of MBTT's alternate minimization algorithm, as reported in [START_REF] Chavent | Reflection fwi in mbtt formulation[END_REF].

We compare now the above exact Θ-estimates with the local Θ-estimates of Subsection 5.3.

Barucq & Calandra & Chavent & Faucher

-For the FWI approach, Figures 19(a), 20(a) and 21(a) and Table 3, it shows that both local and exact Θ-estimates of ∆ u m are of the same size. In sight of the upper bound estimate (32) on which the local Θ-estimate is based, one can think that the FWI formulation corresponds to the worst situation, where the image of a segment in the background space is a curve close to an arc of circle in the data space.

-The situation is completely different for the MBTT formulation: Figures 19(a), 20(b) and 21(b) and Table 3, we see that the exact Θ-estimate ∆ u p is about ten times larger than its local Θ-estimate.

We can also determine the exact attraction basins in the MBTT formulation for the data space reflectivity s at s 0 in the direction u s , which we expect to be large because the forward map F is nearly linear with respect to s. Figure 22 shows the corresponding deflection and global radius of curvatures maps for values of t and t in an interval [-∆, ∆] which is chosen to represent in each case about ±35 times the norm of s 0 defined in [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF]. As expected, the exact Θ-attraction basin is large (23 to 54 times the norm of s 0 depending on frequency), and is 10 5 times larger than its local estimate, which, together with the previous results on the estimation of ∆ u p , confirms the necessity of exact estimates for accuracy.

Influence of reflectivity level

The previous comparisons at 4 Hz and 7 Hz have been performed on the nominal model which has a reflectivity level β of 1% (see [START_REF]Waveform inversion in the laplace fourier domain[END_REF]). We repeat the same computations, at 7 Hz, changing only the reflectivity of the model, which is fixed to a low level, 0.2% in Figure 23 (close to the Born approximation) and to a stronger value, 5% (close to Marmousi reflectivity) in Figure 24. Table 4 summarizes the sizes of and tolerable error for the corresponding attraction basins.

Increasing the reflectivity increases the relative importance of multiple reflections -whose phase are not controlled by the migration/demigration included in the MBTT approach -with respect to primary reflections, whose phase are controlled. It is then expected that the performance of MBTT, when applied to a full wave model as it is the case here, will deteriorate when the reflectivity level increases. This is confirmed by Table 4: the MBTT attraction basin for a reflectivity level of 5 % is three times smaller than at 1 %, which itself is two times smaller than at 0,2 %. On the other hand, the FWI attraction basin remains essentially unaffected by the reflectivity level, and the MBTT basins remain always larger, with a factor four to fifteen (with respect to the Θ-estimates, eight with the R G -estimates). 7 Hz model m model p

Inria

∆ u m 0 m0 R u G,m 0 F (m0) ∆ u p 0 p0 R u G,p 0 F (m0) low reflectivity β = 2 × 10 -3 Exact Θ-

Influence of background smoothness

Here we illustrate an important condition which has to be satisfied for the MBTT decomposition to work. When a background perturbation is applied to p 0 while s 0 is maintained fixed, the phase shifts of the reflections are controlled, provided these reflections are generated only by s 0 . Any backscattered energy generated by the updated background p 0 + tu -or worse, by p 0 itself -will lead to reflections whose phase is not controlled by the migration-demigration process included in MBTT, thus annihilating its benefits. This appears already in Figure 21(b), where one sees that the deflection reaches π/2 (resp. the global radius of curvature equates zero) much closer to the diagonal at the top right corner than at the bottom left one. Indeed, the large positive t, t correspond to strongly increasing backgrounds, which backscatter energy. At the opposite corner of the diagrams (bottom left), which corresponds to mildly increasing velocities, the amount of scattered energy is smaller, and the limiting values are further apart from the diagonal. We have considered for simplicity only attraction basins which are centered at the nominal model m 0 , but on most of the deflection and curvature maps, the attraction basin could, in fact, be extended on the side of negative t (where backgrounds are less steep). Then, if the unperturbed background itself scatters back some energy, as the one illustrated Figure 25, the situation is even worse, as illustrated in the corresponding deflection and curvature maps of Figure 26, at 4 Hz. These maps show no enlargement of the attraction basin in favor of MBTT, in fact, it is even smaller than for FWI. That is why it is fundamental to ensure the smoothness of the background in the MBTT decomposition by using an adapted parametrization, e.g., spline functions are employed in [START_REF] Chavent | Reflection fwi in mbtt formulation[END_REF] to represent p.

-0.1 0 0.1 -0.1 0 0.1 t/ m 0 t / m 0 0 π Θ(t, t ) -0.1 0 0.1 -0.1 0 0.1 t / m 0 0 0.3 R G (t, t )/ F (m 0 ) (a

Behavior of the data misfit

We illustrate in this section the behavior of the FWI and MBTT misfit functions in the model space direction u at the nominal model m 0 = m(p 0 , s 0 ) for some seismic data d. The cost functions are respectively J (m 0 + tu) and J(p 0 + tu, s 0 ) for FWI and MBTT, according to ( 13) and ( 18). This will allow to verify the theoretical results of Section 3: if the minimum of the misfit function is smaller than the tolerable error R u G associated to an attraction basin, the misfit should not have any local minimum over the basin. However, when the data exhibit an error larger than the tolerable error, local minima may appear, but they do not have to. In order to define a seismic data d to build the misfit, we first denote by d 0 = F(m 0 ) = F(p 0 , s 0 ) the synthetic data associated to the nominal model. Note that, due to the choices made in the construction of the nominal model in Subsection 5.1, d 0 is nothing but a migrated-demigrated version of the Marmousi data deprived of the direct arrivals. We define a scaled Marmousi data d m by scaling the synthetic Marmousi data d M = F(m m ) (where m m is given Figure 8(a)) in such a way that d m = d 0 . Then, the data that enter in the FWI and MBTT misfit functionals are defined by:

d = (1 -η)d 0 + ηd m , η ∈ [0, 1]. (57)
When η increases from 0 to 1, the distance of d to the attainable set increases from 0 to d md 0 , which is larger than the tolerable error associated to the attraction basin of interest (FWI or MBTT). The computations are done for the nominal model with 1% reflectivity level (the results for models with 0.2% and 5% reflectivity are similar). As in the previous sections, the nominal model follows Subsection 5.1, and the background perturbation direction u is that of Figure 17. Figure 27 shows the cost functions for FWI and MBTT at 7 Hz for different values of η. We also show by vertical lines the exact R G -attraction basin associated to zero tolerable error, and by horizontal lines the levels 1/2R 2 G corresponding to the tolerable error level associated to the exact Θ-attraction basin (cf. Table 4).

-The FWI and MBTT attraction basins are essentially independent of the error level η in the data. The size of the FWI R G -attraction basins is slightly larger than the one anticipated by ∆ u m0 . When η = 0, d = d 0 is attainable, with a zero error level, and one sees that the first local minimum is just after the right limit of the exact R G -attraction basin, hence consistent with the analysis of Section 3. The tolerable error level R u G,p0 on d for the Θ-attraction basin is attained at η = 0.1.

-The local minima appear much earlier with FWI than with MBTT and the sizes of intervals where there is no local minimum is well anticipated by the estimates ∆ u .

-In the numerical setting considered here, the MBTT approach appears to be robust with respect to data error levels much larger than the tolerable level, as we do not see local minimum appearing inside the attraction basin.

-0.20 -0.10 0.00 0.10 0.20 0 57) for a nominal model with a reflectivity level β = 1%. We also indicate by vertical dotted lines the size of the exact RG-attraction basins, and by horizontal dotted lines the maximal tolerable errors R u G,m 0 associated to the exact Θ-attraction basins, taken from Table 4.

2 4 +∆ u m 0 / m0 -∆ u m 0 / m0 +∆ u p 0 / p0 -∆ u p 0 / p0 1/2(R u G,m 0 ) 2 1/2(R u G,p 0 ) 2 t/ m 0 ; t/ p 0 J (m0 + tu) J(p0 + tu, s0)

Strategy for MBTT waveform inversion

The numerical experiments made in this section lead to the following recommendations, for the implementation of the MBTT reformulation of FWI with a full wave equation, in order to obtain large attraction basins for the background reconstruction:

-the background space for p has to be chosen smooth enough and the reflectivity level β = r / p large enough so that the energy backscattered by p is negligible compared to that backscattered by r;

-the reflectivity level β has to be chosen small enough for the importance of multiples to remain small compared to primary reflections in the synthetics.

This can be achieved by using an adapted parametrization for the background p, e.g. spline functions as employed in [START_REF] Chavent | Reflection fwi in mbtt formulation[END_REF], and by tuning the reflectivity level β by adjusting the migration weight W in [START_REF] Qiu | Full-waveform inversion with an exponentially encoded optimal-transport norm[END_REF].

We remind that the relevance of the chosen settings can be checked before performing any iterative optimization algorithm, by computing the size of the Θ and R G -attraction basins, and the associated tolerable error.

-The attraction basins tend to be larger in the direction of slowly increasing (with depth) backgrounds, which suggests it is better to underestimate the initial slope of the background velocity with respect to depth.

The dependency of the size of attraction basins with the reflectivity level suggests the following organization for MBTT waveform inversion:

1. scale down the given seismic data so that they can be resimulated by a reflectivity level β of circa 0.2-1%, which amounts to bring the forward model close to the Born approximation.

2. Perform alternative minimization with respect to s and p, starting with s first.

3. Rescale (up) the processed seismic data in order to get closer to the original given data, and return to step 2.

4. Continue until processed and given original data coincide.

Note that all the difficulties related to the choice of the background and the reflectivity level disappear if one uses a linearized model (Born approximation) as forward modeling.

Extension for other applications of least squares minimization

The quantities that we have defined above only depend on the least-squares minimization problem we are considering. In particular, we have used it to compare the performance of methods in the context of the acoustic inverse wave problem in seismic. In addition, it is straightforward to apply the computational estimates towards different situations, which do not involve the Helmholtz equation, or the geophysical discrete setup. In this section, other situations are experimented: the inverse problem associated with elastic medium reconstruction, and problems with different kind of domain boundary conditions.

Elastic medium reconstruction

In the case of isotropic elastic propagation, the medium is characterized by three physical models: the Lamé parameters λ and µ, and the density ρ. The propagation of waves in such media is expressed in terms of the vectorial displacement field u, see [START_REF] Akcelik | Parallel multiscale gauss-newton-krylov methods for inverse wave propagation[END_REF]. More precisely, for elastic isotropic media, the equation of wave propagation writes as (where we omit the space dependency for clarity)

-ρω 2 u -∇ λ∇ • u -∇ • µ ∇u + ∇u T = g, in Ω. (58)
In this context the unknown m is composed of the three parameters, m = {λ, µ, ρ} (which are space dependent). The forward problem is accordingly adapted and we write

F elastic (m) = u(x)| Σ , ( 59
)
where we omit the source and frequency dependencies for clarity. Hence, we have considered that the forward problem is given by the measurements of displacement in all directions2 . We consider the models m 0 to be smooth backgrounds, similarly to the previous acoustic case. They are illustrated in Figure 28 and represent a subsurface of size 17 × 3.5 km. The geophysical acquisition is composed of 19 sources and 168 receivers for each source, both being located at the surface. In this case, the derivatives have to be computed with respect to the three parameters. For the directions, we consider the elastic Marmousi2 medium and extract the different reflectors. The directions for the Lamé parameters and the density are presented in Figure 29 where the amplitude has been chosen such that u m = 1. Because of the difference of magnitude between the different quantities, we actually employ a scaling parameter to obtain dimensionless coefficients such that the parameters of interest are (λ/λ 0 , µ/µ 0 , ρ/ρ 0 ), where λ 0 , µ 0 and ρ 0 are scalar coefficients given by the maximal values of the respective quantities. In Figure 30, we show the evolution of the estimates ∆ u m0 and R u m0 with frequency. We see that the estimates for the size of the interval associated with the elastic wave equation, Figure 30(a) behave similarly to the acoustic case: the size of the attraction basin decreases with frequency. However, the maximal distance between the data and the attainable set, given Figure 30(b) appears now to be decreasing with frequency, contrary to the acoustic situation. It would be interesting to pursue the analysis in the elastic framework, in particular to identify how the different parameters are intertwined and how they should be taken.

Extension to different boundary conditions

The estimates we have derived can also be employed for other situations than the inverse wave problem affiliated with a geophysical setup. Let us briefly illustrate the consideration of alternative boundary conditions for the domain of interest. We acknowledge two situations, sketched in Figure 31: 1. completely absorbing medium where the upper free surface is replaced by a perfectly matched layer, m 0 and R u m 0 with frequency for elastic isotropic medium. The initial models m0 are represented with smooth elastic background of Figure 28, the direction is given by the Marmousi structures of Figure 29. The computation of the directional derivatives uses scaled quantities: (λ/λ0, µ/µ0, ρ/ρ0).

Dirichlet boundary condition set to zeros for all boundaries, generating reflection from all sides.

This situation is particularly appropriate in the context of medical imaging and Electrical Impedance Tomography (EIT).

We compute the estimates according to the new situations and show the evolution of ∆ u m0 (obtained with local Θ-estimate) and R u m0 in Figures 32 and33 respectively. The initial model is given by the smooth velocity background of Figure 3 and the directions are those introduced in Figure 4. Namely, we follow the setup of Subsection 4.1, and only the boundary conditions are changed.

The size of the attraction basin, Figure 32, decreases with increasing frequencies independently of the situation. However, in the case of Dirichlet boundaries on all sides, the behavior is much more distorted and reveals a much higher magnitude in the evolution. Namely it has four order of magnitude while only two for the other situation. Also in this case (Dirichlet boundary conditions), there does not seem to be any difference between the two directions, and the large obstacle direction, u 3 , behaves similarly as the Marmousi ones u m . Regarding the maximal distance to the attainable set (largest tolerable eroor on data), Figure 33, we observe a decrease of the distance for low frequencies and then stabilization for the absorbing medium. This medium (fully absorbing) is strongly influenced by the geometry of the direction and the medium with Dirichlet boundary condition renders a totally unsettled evolution. We can intuit that the patterns obtained with the Dirichlet boundary conditions are due to the multiple reflections introduced at the boundaries, which strongly affect the problem in terms of convergence and sensitivity to noise.

Conclusion

We have presented a theoretical and numerical toolbox for the apriori analysis of the optimizability of nonlinear least-squares minimization problems by local algorithms. Such problems arise in seismic Full Waveform Inversion (FWI), where the misfit function tends to exhibit local minima in the directions associated with low spatial frequencies perturbation of the background velocity. But the methodology is totally applicable for other least squares minimization problems.

We have defined attraction basins around a nominal model and the associated tolerable error level such that for any data below tolerable error, one is sure that the data misfit has a single local -and hence global -minimum over the basin. These quantities depend only on the forward map to be inverted, so one can obtain apriori information on the absence of local minima without having to experiment with the misfit function for different data. We have characterized attraction basins by the notion of global radius of curvature in the data space, and given formulas and algorithms for their numerical approximate and exact determination.

We have first applied this tool to the classical FWI problem for a wave equation in the timefrequency domain, and computed directional attraction basins which give quantitative insight into optimizability of FWI:

1. unsurprisingly, the size of attraction basins becomes smaller when frequency increases. In our tests, it is only of ±2% of the norm of the background at 4 Hz, and ±1% at 7 Hz, 2. the use of complex frequencies including damping enlarges the attraction basin, 3. the largest attraction basins obtained for low (or complex) frequency usually go with smaller tolerable error levels. Therefore, low and complex frequencies suffer more from noise.

4. solving the problem for a sequence of increasing single frequency data is better in terms of attraction basin than using an increasing sequence of frequency-bandwidth data, 5. salt domes produce smaller attraction basins, and are harder to recover.

These findings confirm the difficulty inherent to FWI, which requires unrealistically low frequency data for the determination of the full velocity model (Appendix C) In order to see if these limitations could be overcome, we have then applied our tool to the MBTT approach, which uses a background/data space reflectivity reparameterization of the velocity model. Our findings are: 6. in a direction of background perturbation, the largest attraction basin for MBTT shows an increase in size by a factor four to fifteen compared to FWI, and is only barely reduced when frequency increases.

7. Correlatively, the tolerable error for MBTT is divided approximately by 10 compared to FWI. Nonetheless, the comparison of misfit functions has shown that local minimum does not appear when the tolerable error is exceeded.

8. In the direction of data-space reflectivity perturbations, where the forward map is linear up to the multiple reflections, the attraction basin is very large, more than ±20 times the norm of the nominal data-space reflectivity. 9. Critical parameters for an efficient implementation of MBTT are the background smoothness and the reflectivity level.

This provides a strong incentive for the use of the MBTT decomposition to alleviate the low frequency requirement of FWI, despite a larger computational burden. For future work, we plan to implement the MBTT method with these settings to investigate further the efficiency in the context of iterative minimization with velocity model reconstruction.

Theorem 1. [14, Theorem 1] Let g : C × C → C be a function of a complex number z and its conjugate z and let g be analytic with respect to each variable (z and z) independently. Let h : R × R → C be the function of the real variables x and y such that g(z, z) = h(x, y) where z = x + iy. Then the partial derivative ∂ z g (treating z as a constant) gives the same result as

(∂ x h -i∂ y h)/2. Similarly, ∂ z g is equivalent to (∂ x h + i∂ y h)/2.
Corollary 1. Following the statement of Theorem 1, we have

∂g ∂z = ∂g ∂z . (63)
Proof. By direct application of Theorem 1,

∂g ∂z = 1 2 ∂f ∂x + i ∂f ∂y = 1 2 ∂f ∂x -i ∂f ∂y = ∂g ∂z . (64)
We straightforwardly apply the theorem to the misfit function where we consider p := z = x + iy.

J : (x, y) → 1 2 R(x + iy) -d 2 , (65)
where x, y and d can be assimilated with vectors in the discrete setting. Then by deriving independently with respect to x and y we obtain

       ∂J ∂x = 1 2 R * (R(p) -d ) T + 1 2 (R(p) -d ) * R, ∂J ∂y = - i 2 R * (R(p) -d ) T + i 2 (R(p) -d ) * R. (66)
We can further deduce the derivative of J with respect to p and p, where they are considered independent such that J = J(p, p),

       ∂J ∂p = 1 2 R * (R(p) -d ) T = 1 2 (R(p) -d ) T R, ∂J ∂p = 1 2 (R(p) -d ) * R = 1 2 (R(p) -d ) T R (67)
The following theorems give the framework of what can be seen as the chain rule for complex derivation. Theorem 2. Consider the complex-valued function f of a real parameter m and the real-valued functions g 1 and g 2 such that f (m) = g 1 (z(m), z(m)) + ig 2 (z(m), z(m)).

The alternative expression is obtained similarly but by replacing ∂ z g in (71), instead of ∂ z g.

Application of Theorem 3 gives the gradient of the cost function with respect to m,

∇ m J = ∂ ∂m J(m, p) T = 2 Re ∂J ∂p ∂p ∂m T = Re (R(p) -d ) * R ∂p ∂m T = Re ∂p ∂m * R * (R(p) -d ) , (75)
where T stands for the transposed.

A.2 Adjoint state method

In order to avoid the computation of the Jacobian, the gradient is computed with the first order adjoint state method. It has been introduced in the work of [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF], and implemented by [START_REF] Chavent | Identification of functional parameters in partial differential equations[END_REF] for the computation of a functional gradient. The formulation for the elastic wave problem has been carried out by [START_REF] Tarantola | Inversion of travel times and seismic waveforms[END_REF][START_REF] Tarantola | Theoretical background for the inversion of seismic waveforms including elasticity and attenuation[END_REF]. It is a relatively standard techniques nowadays, e.g. [START_REF] Kern | Numerical Methods for Inverse Problems[END_REF], see [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF] for a review in geophysical situations. Yet, the complex variable specification is less common in seismic literature.

In order to compute the derivative ∇J , we formulate the minimization problem (omitting the space dependency)

min m∈M J (m) = J(m, p) subject to Ap = h, (76)
where we introduce the wave operator A, which corresponds to the Helmholtz equation in the acoustic case,

A := -ω 2 c -2 -∆ = -ω 2 m -∆ . (77)
We consider one single source for now, to clarify the indexes, we shall later reintroduce the sum over the sources by linearity, cf. (85). The problem is recast into a formulation with Lagrangian to define

L(m, p, γ) = J(m, p) + A p -h, γ , (78)
where •, • stands for the complex inner product in L 2 such that v, w = v * w, with v * the adjoint. By taking p solution of Ap = h, we have that ∇ m L(m, p, γ) = ∇ m J (m). Furthermore by application of complex derivation Theorem

The adjoint state γ is now selected such that

Re ∂L ∂p + ∂L ∂p = 0, (81)
which gives,

Re ∂J ∂p + ∂J ∂p + A * γ = 0. (82)
We now incorporate [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF], and the adjoint state γ solves the problem

A * γ = -R * R(p) -d . (83)
Using this formulation for γ, the gradient reduces to

∇ m J = Re ∂ m Ap, γ T . (84)
We can reintroduce the sum over the sources, which gives

∇ m J = s Re (∂ m A)p (s) , γ (s) T , where γ (s) solves A * γ (s) = -R * R(p (s) ) -d (s) . ( 85
) Inria A.

Directional derivative computation

For the computation of the directional Fréchet derivative, required for our estimates, the technique is quite different. Note that from the forward operator definition, the Fréchet derivative consists in computing the derivative of the pressure field and to apply the linear restriction R. Thus we focus on the following in the computation of the directional derivative of the pressure field. Here we apply a model perturbation onto the wave equation, in a similar approach as the one prescribed in [START_REF] Demanet | Waves and imaging[END_REF].

Let us consider a model perturbation u for the Helmholtz equation ( 2), where the parameter m is identified with c

-2 , -ω 2 m(x) + u(x) p(x) -∆p(x) = h(x). (86)
The wave field solution is decomposed into

p := p(m + u) = p(m) + p u = p 0 + p u , (87)
where p 0 is solution to the original (non-perturbed) equation ( 2) and the space dependency is omitted for clarity. Developing (86) gives

-ω 2 m -∆ p -ω 2 up = h h + -ω 2 m -∆ p u -ω 2 up = h -ω 2 m -∆ p u = ω 2 up. (88)
Hence, p u is solution to the Helmholtz equation with some right-hand side depending on the total solution p. Assuming ω is not an eigenfrequency of the Helmholtz equation, we write

p u = -ω 2 m -∆ -1 ω 2 up . (89)
Consequently, using the definition p = p 0 + p u , we have

I --ω 2 m -∆ -1 uω 2 p = p 0 , I --ω 2 m -∆ -1 uω 2 -1 p 0 = p, (90)
where I denotes the identity. Direct application of Neumann series provides

I --ω 2 m -∆ -1 uω 2 -1 p 0 = p 0 + ∞ k=1 -ω 2 m -∆ -1 uω 2 k p 0 = p 0 + ∞ k=1 p k . (91)
Using the linearity of the wave equation, we have that p 1 is solution to

-ω 2 m -∆ p 1 = uω 2 p 0 , (92)
and repeatedly, ω 2 m -∆ p j = uω 2 p j-1 .

(93)

By assuming the perturbation u of the model to be small, the Taylor expansion allows us to express (with the notation of [START_REF] Cartan | Differential calculus[END_REF]

We can finally identify the directional derivative in terms of Neumann series and the directional derivative is expressed as solution to a forward problem, i.e. ∂ m p(u) = p 1 ; ∂ 2 m p(u, u) = 2p 2 ; etc. Note that obtaining the derivative from a forward problem can also be found in the context of elastic-fluid interaction in [START_REF] Azpiroz | Characterization of partial derivatives with respect to material parameters in a fluid-solid interaction problem[END_REF], where the derivation is conducted with respect to the Lamé parameters.

B MBTT (Data Space) model parameterization

In this appendix, we explicit how the computations are conducted in the case of MBTT model decomposition.

B.1 Expression of the reflectivity with MBTT

The reflectivity part in the MBTT model representation is given by, cf. [START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF],

r = WDF * 0 s. ( 95
)
The technique developed with the adjoint state method allows to compute of r without explicitly using DF 0 . Indeed, the adjoint state method provides (identification from (75) and ( 85))

where DF (s) stands for the forward operator associated with source s. The fields p (s) and γ (s) solve respectively the forward and ajoint problems, see ( 2) and (85). Proceeding by analogy with the reasoning of Appendix A, it is straightforward to see that r can be express as

r = W s (∂ p A 0 )p (s) 0 * γ (s) 0 T , (97)
where A 0 is the Helmholtz operator with zero reflectivity (i.e., m = p), The terms (∂ p p 0 )(u p ) and (∂ p γ 0 )(u p ) are computed using the same method as presented in Subsection A.3 (thus, each requires the resolution of the wave equation with specific right-hand side). After we have computed u m = (∂ p m)(u p ), we again use Subsection A.3 to compute the final ∂ m F(u m ).

A 0 := -ω 2 p -∆ ; (98)
Eventually, one can process similarly for s, adapting the chain rule. Regarding the second order derivatives, it is analogous with one degree more of derivation in the chain rule. Computationally speaking, it only requires the resolution of additional forward problems.

C Influence of the geometry and frequency on FWI

We illustrate the influence of the subsurface geometry with an acoustic FWI experiment, where we demonstrate the effect we have depicted in Section 4, Figure 5: 1. a high contrast object complicates the process of reconstruction compared to a Marmousi-like subsurface, by reducing the size of the basin of attraction, 2. the low frequencies, by increasing the size of the basin of attraction, provide benefits on the reconstruction.

We design an experiment for a two-dimensional area of size 12.3 × 3 km, considering an acoustic medium with constant density set to 1000 kg m -2 . The wave propagation is governed by the Helmholtz equation [START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF], one parameter has to be recovered: the velocity c. The forward operator, see [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF], provides pressure measurement at the receivers location. For the seismic acquisition, we take 135 sources and 306 receivers per source. We take synthetic data, avoiding noise, simply to illustrate our analysis.

C.1 Comparison between salt and Marmousi geometries

We compare the subsurface reconstruction using FWI depending on the geometry of the subsurface. We consider a medium involving a high contrasting object, and the Marmousi model. The iterative minimization is conducted following a nonlinear conjugate gradient algorithm (see [START_REF] Nocedal | Numerical Optimization[END_REF]) and using frequencies from 1 to 10 Hz, with 1 Hz increment. We perform 20 iterations per frequencies, which are taken sequentially. The corresponding FWI algorithm is shown in Algorithm 1.

For the experiment related to the Marmousi medium, the 'true' model has been given in Figure 8(a) and the initial guess in Figure 3. Here we artificially modify the dimension so that the models are now of size 12.3 × 3 km. For the model encompassing a salt dome, the initial and target models are given in Figures 34(a In Figures 34(c) and 35, we compare the reconstructions from FWI regarding the contrasting object or the Marmousi model. Both situations use media of same size, with the same frequency progression (sequential from 1 to 10 Hz), and start from initial models where no information is known (one dimensional wave speed variation, with depth only). However, only the Marmousi case is able to accurately capture the subsurface structures. It demonstrates the results of Section 4 that, namely, high contrast objects, by reducing the size of the basin of attraction, are more complicated to reconstruct.

Initial parameters: initial model c 1 and the frequency list of ω i , for i = 1, . . . , n ω . The frequencies are ordered such that ω 1 < ω 2 < . . . < ω nω , see Remark 4. Frequency loop for i ∈ {1, . . . , n ω } do Optimization loop for j ∈ {1, . . . , n iter } do set k := (i -1)n iter + j -solve the wave equation at frequency ω i with wave speed c k ; -compute the cost function by comparing the simulation and observation data; -compute the gradient of the cost function, using the adjoint-state method of Appendix A; -compute the search direction, s k with the nonlinear conjugate gradient method (see [START_REF] Nocedal | Numerical Optimization[END_REF]); -compute the step length α with line search method (see [START_REF] Nocedal | Numerical Optimization[END_REF]); -update the wave speed c k+1 = c kαs k ; end end Algorithm 1: FWI iterative minimization algorithm for the wave speed reconstruction.

C.2 Low frequencies compensation

It has been shown in Figure 5 that low frequencies enlarge the size of the radius. To illustrate this effect, we reproduce the previous experiment but incorporating nine frequencies from 0.1 to 0.9 Hz (sequentially using 0.1 Hz increment). Then, we employ the set of frequencies from 1 to 10 Hz (sequentially using 1 Hz increment). In Figure 36, we show the reconstruction starting from these very low (unrealistic) frequencies.

The FWI algorithm benefits from this incorporation of low frequencies, and is now able to accurately retrieve the contrasting objects, as expected. Here we have used synthetic experiment to demonstrate the basin of attraction dependency on the geometry and frequency. Note that we could as well have used complex frequencies to increase the basin of attraction in early stages, cf. Subsection 4.5. However, for realistic application, low frequencies are usually unavailable (due to the noise) and it is then impossible to reduce arbitrarily the initial frequency. As an alternative, techniques based on shape reconstruction could be very efficient in this context but, because the nature of the target is unknown at the beginning, it is hard to justify apriori.

  illustrates the configuration. Computational domain using Perfectly Matched Layers (PML) at the lateral and bottom boundaries.

  Computational domain using absorbing boundary condition on Γ 2 .

Figure 1 :

 1 Figure 1: Illustration of the two-dimensional area of interest. It characterizes a geophysical situation where sources and receivers are located in the upper part. The boundary is decomposed into a free surface at the top (Γ1); on the sides and bottom, conditions (Γ2) are implemented to avoid returning waves.

Figure 3 :

 3 Figure 3: Acoustic model c0 of size 9.2 × 3 km, the velocity is indicated in km s -1 . The estimates are calculated at the nominal model m0 = c -2 0 .

  u 2 : two reflectors.

  u 3 : single large reflector.

  um: reflectors following Marmousi structure.

Figure 4 :

 4 Figure 4: Four different directions for the estimation of the quantities of interest (R u m 0 and ∆ u m 0 ). The amplitude is determined such that u k = 1. The values of u k are given in (m s -1 ) -2 = s 2 m -2 .

  Evolution of the maximal distance between the data and the attainable set, R u m 0 .

Figure 5 :

 5 Figure 5: Evolution of ∆ u m 0 and R u m 0 with frequency using m0 as a smooth velocity background of Figure 3, and using different directions u. The blue squares ( ) employ one reflector (u1, Figure 4(a)), the green circles ( ) two reflectors (u2, Figure 4(b)) the red stars ( ) one large reflector (u3, Figure 4(c)) and the black squares ( ) the reflectors extracted from the Marmousi model (um, Figure 4(d)).

1 R 1 R

 11 G (t, t )/ F (m 0 ) (a) Using the Marmousi direction um and frequency 4Hz. G (t, t )/ F (m 0 ) (b) Using the Marmousi direction um and frequency 7Hz.

Figure 6 :

 6 Figure 6: Maps of the deflection (30) and global radius (26) between two models of M(m0, um, ∆), see (48). The selected direction um is Figure 4(d), m0 and Figure 3. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

1 R 1 R

 11 G (t, t )/ F (m 0 ) (a) Using the direction u 3 and frequency 4Hz. G (t, t )/ F (m 0 ) (b) Using the direction u 3 and frequency 7Hz.

Figure 7 :

 7 Figure 7: Maps of the deflection (30) and global radius (26) of M(m0, u3, ∆), see (48). The selected direction u3 is Figure 4(c), m0 and Figure 3. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

Figure 8 :

 8 Figure 8: Acoustic models of size 9.2 × 3 km used as initial models for the estimation of ∆ u m j , we take mj = c -2 j , for j = {m, c}. The velocity is indicated in km s -1 .

3 (

 3 a) Estimates using mm of Figure8(a).

3 (

 3 b) Estimates using mc of Figure 8(b).

Figure 9 :

 9 Figure 9: Evolution with frequency of ∆ u m from local Θ-estimate for two initial models mm (left) and mc (right), in the direction of u1 and u3 (Figure 4).

Figure 3 ,

 3 Figure3, this leads to the model of Figure10(a), where the number of coefficients has been reduced to 56. Using this reduced representation, we can compute all quantities to effectively obtain the lower bounds ∆ m0 and R m0 (i.e., we compute the matrix D 2 F and the singular values of DF).Figure10(b) shows an important reduction of the size ∆ m0 of the n-dimensional basin compared to the directional estimates in the Marmousi direction u m using the local Θ-estimate (six orders of magnitude). Yet, it confirms the main property, which is the decrease of the attraction basin size with increasing frequencies. A similar observation holds for the tolerable data error level R m0 , see Figure10(c), which suffers an important reduction compared to the directional estimates. In addition, we observe a different pattern between R m0 , which decreases with frequency, and the directional estimates R u m0 , which increases with frequency. As illustrated in Figure5, R u m0 depends more on the direction than ∆ u m0 , and it is reflected on the n-dimensional estimates, which takes the worst case.

  Compression of the acoustic starting model of size 9.2 × 3 km, where only 56 coefficients are kept to represent the velocity, indicated in km s -1 . Evolution of ∆ u m 0 computed with local Θ-estimate.

7 Frequency

 7 (Hz) R u m0 / d (c) Maximal distance between the data and the attainable set.

Figure 10 :

 10 Figure 10: Lower bounds and directional estimates of ∆ u m 0 and R u m 0 using a model with reduced dimensions. The respective lower bounds are given with the blue squares (), the black squares ( ) show the directional estimates using um (reflectors from the Marmousi model, see Figure4(d)).

  Evolution of ∆ u m 0 using local Θ-estimate. The large blue bars ( ) and the red dotted line ( ) employ the direction um. The smaller orange bars ( ) and the black dotted line ( ) employ direction u 3 . Evolution of R u m 0 , the large blue bars ( ) and the red dotted line ( ) employ the direction um The smaller orange bars ( ) and the black dotted line ( ) employ direction u 3 .

Figure 11 :

 11 Figure 11: Evolution of ∆ u m 0 and R u m 0 with sequential or group of frequency. We use m0 as a smooth velocity background of Figure 3. The histograms represent the estimations for groups of ten frequencies and cover the respective interval of frequency. The lines represent the estimations using sequential frequency. The different directions u are either extracted from the Marmousi model (um, Figure 4(d)) or from one large reflector (u3, Figure 4(c)).

  Estimates for the complex frequency plane. Estimates using fixed f .

Figure 12 :

 12 Figure 12: Map of the local Θ-estimate of ∆ u m 0 with complex frequencies using m0 as the smooth background velocity of Figure 3 and the direction um of Figure 4(d). The angular frequency is defined by -ω 2 = (σ + 2iπf ) 2 .

Figure 13 :

 13 Figure 13: Map of the local Θ-estimate of ∆ u m 0 with complex frequencies using m0 as the smooth background velocity of Figure 3 and the direction u3 of Figure 4(c). The angular frequency is defined by -ω 2 = (σ+2iπf ) 2 .

  d (b) Estimates using u 3 .

Figure 14 :

 14 Figure 14: Evolution of the tolerable error level R u m 0 with complex frequencies using m0 as the smooth background velocity of Figure 3. The angular frequency is defined by -ω 2 = (σ + 2iπf ) 2 .

  (a). The resulting stability constant estimates are pictured Figure 15. We observe that Barucq & Calandra & Chavent & Faucher -the stability constant estimates decrease (i.e. the expected accuracy of the reconstruction is improved) with increasing frequency f .

2 CFigure 15 :

 215 Figure 15: Numerical estimates for the stability constant C in (50) with complex frequencies using m0 as the smooth background velocity of Figure 3 and m † as the Marmousi model of Figure 8(a). The angular frequency is defined by -ω 2 = (σ + 2iπf ) 2 .

  Model r 0 at 2 Hz frequency. Model r 0 at 4 Hz frequency.

  Model r 0 at 7 Hz frequency.

  Model r 0 using frequencies from 0.5 to 15 Hz with 0.5 Hz increment.

Figure 16 :

 16 Figure 16: Reflectivity models r0 obtained from the MBTT representation defined by (52). The model s is defined from (51) as the difference between the data obtained from the Marmousi model Figure 8(a) and the smooth background Figure 3. The figures correspond with squared slowness and the values are given in (m s -1 ) -2 .

Figure 17 :

 17 Figure 17: Perturbation u used for the background model p. The amplitude is determined such that u = 1 and the values are given in (m s -1 ) -2 = s 2 m -2 .

Figure 18

 18 Figure18shows the corresponding unperturbed and perturbed synthetic data for a single, centrally located source at frequency 4 Hz. It corresponds to the solution of the Helmholtz equation (2) recorded at the receivers. Note that, from[START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] the direct arrivals are removed from the forward operator. We see that the modification of the background part p preserves the phase of the original signal, and only the amplitudes are modified. On the other hand, when the perturbation is applied to the global model m, both phase and amplitude of the signal are modified.

Figure 18 :

 18 Figure 18: Comparison of the synthetic data associated with a centrally located shot at 4 Hz using a model perturbed by the direction u of Figure 17 applied onto the global model m or on the part p using the MBTT model decomposition. The step for the perturbation is T = 5 × 10 -5 .

  Perturbation of the reflectivity s.

Figure 19 :

 19 Figure 19: Evolution with frequency of the local Θ-estimates of the size of the attraction basins in the context of FWI or MBTT model space parametrization, see (42) and[START_REF] Shin | Sequentially ordered single-frequency 2-d acoustic waveform inversion in the laplace-fourier domain[END_REF]. Here p0 is the smooth velocity background of Figure3, the direction u for p is given Figure17, and the direction us for s is a random vector. In the MBTT representation, the reflectivity uses only the selected frequency.

  FWI: perturbation of the global model m. MBTT: perturbation of the background p.

Figure 20 :

 20 Figure 20: Maps of the deflection (30) and global radius (26) between two perturbed velocity or background models at frequency 4 Hz. The perturbation direction is the ramp of Figure 17, it is either applied to the global model m or to the background parameter p. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

  , 21 and 22. -The first observation is that lower values of deflection are achieved when the background perturbation u is applied to p (MBTT) rather than to m (FWI): at 4 Hz, Figures 20(b), it never FWI: perturbation of the global model m. MBTT: perturbation of the background p.

Figure 21 :

 21 Figure 21: Maps of the deflection (30) and global radius (26) between two perturbed velocity or background models at frequency 7 Hz. The perturbation direction is the ramp of Figure 17, it is either applied to the global model m or to the background parameter p. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

Figure 22 :

 22 Figure 22: Maps of the deflection (30) (left) and global radius (26) (right) at frequency 7 Hz between two perturbed data-space reflectivities, for the perturbation direction us chosen in section 5.2 The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

  FWI: perturbation of the global model m. MBTT: perturbation of the background p.

Figure 23 :

 23 Figure 23: Maps of the deflection (30) and global radius (26) between two perturbed velocity or background models at frequency 7 Hz. Compared to Figure 21, the reflectivity β defined by (54) has been divided by five. The perturbation direction is the ramp of Figure 17, it is either applied to the global model m or to the background parameter p. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

  ) FWI: perturbation of the global model m. MBTT: perturbation of the background p.

Figure 24 :

 24 Figure 24: Maps of the deflection (30) and global radius (26) between two perturbed velocity or background models at frequency 7 Hz. Compared to Figure 21, the reflectivity β defined by (54) has been multiplied by five. The perturbation direction is the ramp of Figure 17, it is either applied to the global model m or to the background parameter p. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

Figure 25 :

 25 Figure 25: Acoustic smooth Marmousi model c0 where the velocity is indicated in km s -1 . The estimates are calculated using p0 = c -2 0 .

Figure 26 :

 26 Figure 26: Maps of the deflection (30) and global radius (26) between two models at frequency 4 Hz. The model p0 is shown Figure 25 and is the only difference in the computation compared to the maps shown Figure 20. The perturbation direction is the ramp of Figure 17, and is either applied to the global model m or to the background p. The black lines indicate when the deflection becomes higher than π/2, the white lines indicate when the global radius becomes 0.

  0 ; t/ p 0 (c) Using η = 1.

Figure 27 :

 27 Figure 27: Comparison of misfit functional J (m0 + tu) (blue) and J(p0 + tu, s0) (green) at 7 Hz, for different data d as defined in (57) for a nominal model with a reflectivity level β = 1%. We also indicate by vertical dotted lines the size of the exact RG-attraction basins, and by horizontal dotted lines the maximal tolerable errors R u G,m 0 associated to the exact Θ-attraction basins, taken from Table4.

  P-velocity (km s -1 ).

  S-velocity (km s -1 ).

  density (10 3 kg m -2 ).

Figure 28 :

 28 Figure 28: Elastic smooth model of size 17 × 3.5 km used for m0 in the estimates. The relations between the velocities and Lamé parameters are cp = (λ + 2µ)/ρ and cs = µ/ρ.

Figure 29 :

 29 Figure 29: The elastic directions extracted from the Marmousi model. The amplitude is accordingly selected to give u = 1 for all coefficients (λ, µ and ρ).

  Evolution of ∆ u m 0 computed with local Θ-estimate. Evolution of the maximal distance between the data and the attainable set, R u m 0 .

Figure 30 :

 30 Figure 30: Evolution of ∆ u

  Computational domain using Dirichlet boundary conditions only.

Figure 31 :

 31 Figure 31: Illustration of different boundary conditions for the computational domain used to estimate the quantity ∆ u m 0 .

Figure 32 :

 32 Figure 32: Evolution of ∆ u m 0 computed with the local Θ-estimate with frequency for different kinds of boundary condition problems. We use m0 as the smooth background model of Figure 3; the directions are given Figure 4.

  Medium where all boundaries have a Dirichlet condition set to 0.

Figure 33 :

 33 Figure 33: Evolution of the maximal distance between the data and the attainable set, quantity R u m 0 , with frequency for different kinds of boundary condition problems. We use m0 as the smooth background model of Figure 3; the directions are given Figure 4.

2 ∂ 2 2 ∂

 222 ), p(m + u) = p(m) + ∂p(m) ∂m (u) + 1 p(m) ∂m 2 (u, u) + . . . = p 0 + ∂p(m) ∂m (u) + 1 2 p(m) ∂m 2 (u, u) + . . .

s(

  ∂ m A)p (s) , γ (s) T = s ∂p (s) ∂m * R * (R(p (s) )d (s) ) ⇒ s ∂ m A)p (s) * γ (s) T = DF (s) (m) * (R(p (s) )d (s) ),

0s 2

 2 solves the forward problem with A 0 for the source s, and γ 0 solves for the source s, R * s (s) , (99) where s writes as s = {s(1) . . . s (nsrc) }. The model representation[START_REF] Brenders | Comparison of 3d time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy[END_REF] becomesm(p, s) = p + W s (∂ p A 0 )pProcess summaryIn order to express the model and the gradient of the cost function in this case, the workflow is as follows 1. start from prescribed representation p and s, 2. solve the Helmhotz equation with zero-reflectivity for all sources s, A 0 p R * s (s) , 4. compute the model m with (100), 5. solve the Helmholtz equation for p (s) : Ap (s) = h (s) , 6. solve the adjoint problem for γ (s) :A * γ (s) = -R * R(p (s) )d (s) .Finally, the gradient of the cost function, with respect to p or s is expressed from the adjoint state method,Re (∂ p A)p (s) , γ (s) T , ∇ s J = s Re (∂ s A)p (s) , γ (s) Directional derivative computationFor the estimation of the size of the basin of attraction and of the radius of curvature, we need the directional derivative of the forward operator. Using Neumann series and Taylor expansion, we have given in Subsection A.3 their expression using the traditional model formulation. When the model is decomposed using MBTT, one can process similarly, it simply requires a few more steps. For clarity we only focus on the parameter p, we start with the chain rule,

  ) and 34(b) respectively.

  Target wave speed.

  FWI wave speed reconstruction using frequencies from 1 to 10 Hz.

Figure 34 :

 34 Figure 34: 2D reconstruction of wave speed including salt domes with data covering frequencies from 1 to 10 Hz. The model is of size 12.3 × 3 km.

1 Figure 35 :

 135 Figure 35: 2D reconstruction of the Marmousi wave speed with synthetic data covering frequencies from 1 to 10 Hz. The dimensions are artificially increased to the size 12.3 × 3 km compared to the true and starting media, given Figures 8(a) and 3 respectively.

Figure 36 :

 36 Figure 36: 2D reconstruction of wave speed including salt domes with data covering frequencies from 0.1 to 10 Hz. The initial model used for the algorithm is given in Figure 34(a), and the true model Figure 34(b). The algorithm takes full benefits of the low frequency content in the data because we use synthetic measurements without noise.

Table 1 :

 1 Size ∆ of attraction basins centered at m0 using directions um and u3 (see Figure4) and corresponding maximal tolerable error RG, at 4 and 7 Hz. By construction, the RG-estimates correspond to the limit case of a zero tolerable error, the other values are extracted fromFigures 5, 6 and 7.

	4 Hz
	Local Θ-estimates	0.018	0.10	0.012	0.03
	Exact Θ-estimates	0.002	0.03	0.02	0.018
	Exact R G -estimates	0.003	0	0.05	0
	7 Hz
	Local Θ-estimates	0.01	0.15	0.008	0.03
	Exact Θ-estimates	0.01	0.05	0.01	0.016
	Exact R G -estimates	0.03	0	0.03	0

Table 2 :

 2 Evolution of ∆ u

Table 3 :

 3 Size

	4 Hz
	Local Θ-estimates Exact Θ-estimates	0.02 0.02	1.6 0.6	0.022 0.2	0.8 0.05	2 × 10 -3 54	6.5 6.5
	Exact R G -estimates	0.05	0.0	0.23	0	60	0
	7 Hz
	Local Θ-estimates Exact Θ-estimates	0.01 0.01	1.6 0.6	0.014 0.11	0.7 0.06	1 × 10 -3 23	4.9 4.9
	Exact R G -estimates	0.025	0.0	0.20	0	>35	0

∆ of attraction basins centered at m0 and corresponding maximal tolerable error RG for the different estimations, at 4 and 7 Hz. By construction, the RG-estimates correspond to the limit case of a zero tolerable error, the other values are extracted from

Figures 19,[START_REF] Chavent | Data Space Reflectivity and the Migration based Travel Time approach to FWI[END_REF]

Table 4 :

 4 Comparison of the size ∆ of attraction basins centered at m0 and corresponding maximal tolerable error RG for the different levels of reflectivity at 7 Hz. By construction, the RG-estimates correspond to the limit case of a zero tolerable error, the values are extracted fromFigures 21, 23 and 24.

	estimates	0.011	0.648	0.155	0.060
	Exact R G -estimates	0.025	0	0.214	0
	Exact Θ-estimates	0.011	medium reflectivity β = 1 × 10 -2 0.645 0.107	0.060
	Exact R G -estimates	0.025	0	0.202	0
	Exact Θ-estimates	0.008	strong reflectivity β = 5 × 10 -2 0.580 0.037	0.124
	Exact R G -estimates	0.028	0	0.042	0

RR n°9253

Barucq & Calandra & Chavent & Faucher

Note that with s 0 given by (51), the resulting r 0 (ω) is the gradient of the (L2 norm) cost function at frequency ω using the model p 0 , see Appendix A.RR n°9253

We have also tried experiments restricting the available data to the vertical direction (i.e., uz), which may be more relevant in seismic, and the results are the same. Hence, we only present the multi-direction data for simplicity.

avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

A Adjoint state using complex variables

The Full Waveform Inversion (FWI) method aims the reconstruction of the subsurface medium parameters by an iterative minimization of the cost function defined as the difference between simulation and observations. We follow the standard least squares formulation of [START_REF] Bozdağ | Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements[END_REF], and consider the Helmholtz equation ( 2) to write

where the forward problem is written with the restriction operator R, and we use the index (s) for the sources. For the minimization, one needs to obtain the gradient of the cost function, which is usually obtained using adjoint state method, see [START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF] for a review of the method in geophysical application.

We specify now the computations for complex-valued fields, omitting the source and frequency sums for clarity, and write

In the frequency domain, p is complex, which requires some precaution for the application of the adjoint state method. In particular, the functional

is not analytic (holomorphic) with respect to the field p. A workaround is relatively standard, see for example [START_REF] Brandwood | A complex gradient operator and its application in adaptive array theory[END_REF][START_REF] Li | Optimization in the complex domain for nonlinear adaptive filtering[END_REF][START_REF] Kreutz-Delgado | The complex gradient operator and the cr-calculus[END_REF], with elements of complex calculus based on Wirtinger calculus. We believe it is important to mention this aspect which is too often disregarded in seismic applications and hereby present the steps involved.

A.1 Complex derivation

The derivation of complex functional is conducted by taking independently the complex variable and its conjugate, respectively z and z, with z = x + iy.