Florian Faucher
email: ffaucher@math.purdue.edu

Maarten V De Hoop
email: mdehoop@purdue.edu

Henri Calandra
email: henri.calandra@total.com

Lingyun Qiu ¶

Christian Rivera
email: christian.rivera@total.com

MULTI-LEVEL, MULTI-FREQUENCY, COMPRESSED FULL WAVEFORM INVERSION *

We study the inverse boundary value problem for the wave equation and recovery of the wavespeed upon taking a time-Fourier transform of the data. We design a hierarchical compressed reconstruction in a multilevel scheme for the inverse boundary value problem associated with the Helmholtz equation using the Dirichletto-Neumann map, or the single-layer potential operator, at selected frequencies as the data. The compression is based on a domain partitioning of the subsurface, while the hierarchy is straightforwardedly established through refinement. The coefficients are assumed to be piecewise constant functions following the domain partitioning, thus allowing the presence of conormal singularities (causing reflections and diffractions); these guarantee a Lipschitz stability estimate for the inverse problem which gives a radius of convergence, related to the stability constant, of the scheme even in the case of partial data. The stability constant grows exponentially with the number of subdomains in the domain partitioning, whence the compression. The frequencies are selected to control the approximation errors in the data due to the compression at the different levels in the scheme. The stability constants (number of subdomains) and approximations errors (frequencies) are coupled through a condition between levels for convergence of the multi-level scheme. We use Haar wavelets for the compression in alternative strategies providing a gradual increase in the number of subdomains in the domain parititioning via adaptive, local multi-scale refinement. We include visco-acoustic behavior and carry out numerical experiments.

1. Introduction. We study the inverse boundary value problem for the wave equation and recovery of the wavespeed, including attenuation, upon taking a time-Fourier transform of the data. We follow a hierarchical compressed reconstruction in a multi-level scheme for the inverse boundary value problem associated with the Helmholtz equation using the Dirichlet-to-Neumann map as the data [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 effectiveness and convergence of this scheme was analysed in [START_REF] Beretta | Lipschitz stability of an inverse boundary value problem for a Schr\"{o}dinger type equation[END_REF][START_REF] Beretta | Inverse boundary value problem for the helmholtz equation with multi-frequency data[END_REF]. The wavespeed can be complex and frequency dependent and, hence, visco-acoustic behaviour is included. Density can, in principle, be recovered via the reconstruction of a potential (using its frequency dependency) in the corresponding Schrödinger equation. We do not address the inverse source (signature) problem here.

Iterative methods for the above mentioned inverse boundary value problem have been collectively referred to, in reflection seismology, as full waveform inversion (FWI). (The term 'full waveform inversion' was supposedly introduced in [START_REF] Pan | Full-waveform inversion of plane-wave seismograms in stratified acoustic media: Theory and feasibility[END_REF] with reference to the use of full seismograms information.) 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][START_REF]Inverse Problem Theory: methods for data fitting and and model parameter estimation[END_REF] introduced the formulation of the seismic inverse problem as a local optimization problem encompassing a least-squares (L 2 ) minimization of a misfit functional. We also mention the original work of Bamberger, Chavent & Lailly [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] in the onedimensional case. Tarantola noticed that the evaluation of the gradient of the mentioned functional corresponds with an imaging procedure with a cross correlation imaging condition. Initial numerical experiments in the two-dimensional case were carried out by Gauthier [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: Numerical results[END_REF]. Since then, a range of alternative misfit functionals have been considered; we mention, here, the criterion derived from the instantaneous phase [START_REF] Bozdag | Misfit functions for full waveform inversions based on instantaneous phase and envelope measurements[END_REF].

The time-harmonic formulation was initially promoted by Pratt and other collaborators in [START_REF] Pratt | Inverse theory applied to multi-source cross-hole tomography. part 1: Acoustic wave-equation method[END_REF][START_REF] Pratt | Combining wave-equation imaging with traveltime tomography to form highresolution images from crosshole data[END_REF] and specifically for wide-angle reflection data in [START_REF] Pratt | Two-dimensional velocity models from wide-angle seismic data by wavefield inversion[END_REF]. The Newton method using the discretized Helmholtz equation as a point of departure was developed by Pratt, Shin & Hicks [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF]. Later the use of complex frequencies was studied in [START_REF] Shin | Waveform inversion in the laplace-fourier domain[END_REF][START_REF] Ha | Acoustic full waveform inversion of synthetic land and marine data in the laplace domain[END_REF]; misfit functions using a logarithmic norm were explored in [START_REF]Inverse Problem Theory: methods for data fitting and and model parameter estimation[END_REF][START_REF] Shin | Waveform inversion using a logarithmic wavefield[END_REF][START_REF] Shin | Comparison of waveform inversion, part 1: conventional wavefield vs logarithmic wavefiled[END_REF].

The nonlinearity of FWI has motivated studies to develop hierarchical multiscale strategies to mitigate this nonlinearity. In the time domain, Bunks et al. [START_REF] Bunks | Multiscale seismic waveform inversion[END_REF] proposed successive inversion of data subsets of increasing frequency contents following the intuition that low frequencies are less sensitive to cycle skipping. This multiscale approach can, for example, be related to the subspace search method advocated by Kennett, Sambridge & Williamson [START_REF] Kennett | Subspace methods for large inverse problems with multiple parameter classes[END_REF]. The time-harmonic formulation, however, has been recognized as the more natural one for this purpose. The idea of frequency progression has been considered, amongst others, by Sirgue & Pratt [START_REF] Sirgue | Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies[END_REF] and extensively studied in the electromagnetic waves case, for example, by Bao and Li [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]. The application of wavelet bases in compressing the successive models in the iteration has been considered by [START_REF] Loris | Tomographic inversion using l1-norm regularization of wavelet coefficients[END_REF][START_REF] Loris | Nonlinear regularization techniques for seismic tomography[END_REF] in wave-equation tomography (sparsity promoting optimization) and in FWI by Lin, Abubakar & Habashy [START_REF] Lin | Seismic full-waveform inversion using truncated wavelet representations[END_REF] (for the purpose of reducing the size of the Jacobian). In our approach the compression, in particular, the compression rate, plays a fundamental role in guaranteeing convergence of our multi-level scheme.

Malinowski, Operto & Ribodetti [START_REF] Malinowski | High-resolution seismic attenuation imaging from widaperture onshore data by visco-acoustc frequency-domain full-waveform inversion[END_REF] demonstrate high-resolution imaging of attenuation and phase velocity in the visco-acoustic case. In this context, we also mention the work of Askan, Akcelik, Bielak & Ghattas [START_REF] Askan | Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures[END_REF]. Our convergence analysis applies to complex wavespeeds.

In FWI one commonly applies a 'nonlinear' conjugate gradient method, a Gauss-Newton method, or a quasi-Newton method (L-BFGS; for a review, see Brossier [START_REF] Brossier | Imagerie sismique à deux dimensions des milieux viso-élastiques par inversion des formes d'ondes : développements méthodologiques et applications[END_REF]). For the application of multiscale Newton methods, see Akcelik [START_REF] Akcelik | Parallel multiscale gauss-newton-krylov methods for inverse wave propagation[END_REF]. In this paper, we carry out numerical experiments with the nonlinear conjugate gradient method. In the gradient method, the step length is typically estimated by a simple line search for which a linearization of the direct problem is used (Gauthier, Virieux & Tarantola [START_REF] Gauthier | Two-dimensional nonlinear inversion of seismic waveforms: Numerical results[END_REF]). This estimation is challenging in practice and may lead to a failure of convergence. For this purpose, research based on trust region has been studied by Eisenstat & Walker in [START_REF] Eisenstat | Globally convergent inexact newton methods[END_REF], for FWI see Métivier and others in [START_REF] Métivier | Full waveform inversion and the truncated newton method[END_REF], the method is detailed in [START_REF] Conn | Trust Region Methods[END_REF]. In various approaches one accounts for the Hessian, at least its diagonal, in the Gauss-Newton scheme [START_REF] Shin | Improved amplitude preservation for prestack depth migration by inverse scattering theory[END_REF]. The second term in the Hessian (Newton method) can be interpreted in terms of multiple scattering (see, for example, Pratt, Shin & Hicks [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF]). In certain earthquake seismology applications, one builds the Fréchet derivative or Jacobian (sensitivity, for example, [START_REF] Chen | Full three-dimensional tomography: a comparison between the scatteringintegral and adjoint-wavefield methods[END_REF]) explicitly and then applies LSQR.

1.1. Time-harmonic waves. We formulate FWI in the frequency domain as an inverse boundary value problem for the Helmholtz equation. We consider a bounded domain Ω of R 3 representing the region of interest in the subsurface; the data are acquired on an open portion of the boundary, ∂Ω. Time-harmonic waves are solutions, u, of (1.1)

-(∆ + q(x)) u = 0 on Ω, u = g on ∂Ω.
Here, q = q(x) = k(x) 2 = ω 2 c(x) -2 with ω representing the (possibly complex) angular frequency and c = c(x) the wavespeed; g = g(x, ω) is a boundary source. For the discretization of the Helmholtz operator we employ a finite difference method with a compact stencil, and a massively parallel structured multifrontal solver, Hsolver [START_REF] Wang | Acoustic inverse scattering via helmholtz operator factorization and optimization[END_REF][START_REF]On 3d modeling of seismic wave propagation via a structured parallel multifrontal direct helmholtz solver[END_REF][START_REF] Wang | Massively parallel structured multifrontal solver for timeharmonic elastic waves in 3-d anisotropic media[END_REF]. This solver is naturally designed for solving the Helmholtz equation for many right-hand sides.

Attenuation. To incorporate visco-acoustic behaviour of subsurface materials, we introduce complex wavespeeds leading to attenuation. For a comprehensive review of models of attenuation, we refer to [START_REF] Ursin | Comparison of seismic dispersion and attenuation models[END_REF]. Suppressing the x-dependencies in the notation, the wave number k becomes:

k(ω) = ω c(ω) = ω c p (ω) + iα(ω),
where c p represents the phase velocity. We consider the Kolsky-Futterman model ( [START_REF] Kolsky | The propagation of stress pulses in viscoelastic solids[END_REF], [START_REF] Futterman | Dispersive body waves[END_REF]) which has been widely used in applications and is given by

           1 c p (ω) = 1 c r + log ω r ω πc r Q r , α(ω) = |ω| 2c r Q r ,
where c r and Q r represent the values of the the phase velocity c p and the quality factor Q at a reference angular frequency, ω r . If αc p 2|ω|,

Q(ω) = Q r + 1 π log ω r ω .
If Q and c p were both frequency independent, then c p = c r , Q r = Q, and the equation for the complex wavespeed may be simplified to

c(ω) = 4c p Q 2 -2i sgn(ω) c p Q 4Q 2 + 1 .
In the case of weak attenuation, 4Q 2 1, the expression on the right-hand side can be approximated as

(1.2) c(ω) = c p 1 -i sgn(ω) 1 2Q .
This simplified representation was used in [START_REF] Malinowski | High-resolution seismic attenuation imaging from widaperture onshore data by visco-acoustc frequency-domain full-waveform inversion[END_REF] while studying full waveform inversion with attenuation.

1.2. Modelling the data and mismatch functional. The data generated by the boundary sources, g, can be represented by the so-called single layer potential operator, S q . That is,

(1.3) S q : g → ∂Ω G q (x, y)g(y) dσ(y), x ∈ ∂Ω,
where G q is the fundamental solution of -(∆ + q(x)) in R 3 . This operator is discussed and analyzed in [START_REF] Park | On time-harmonic seismic data and blending in full waveform inversion[END_REF]. We note that the single-layer potential operator and the so-called Dirichlet-to-Neumann map represent equivalent data [START_REF] Park | On time-harmonic seismic data and blending in full waveform inversion[END_REF]. The action of S q has been referred to as source blending [START_REF] Berkhout | Changing the mindset in seismic data acquisition[END_REF]. The forward problem is then described by

F : ω 2 c(x) 2 = q(x) → S q
In [START_REF] Park | On time-harmonic seismic data and blending in full waveform inversion[END_REF] it was shown that S q1 -S q2 is a Hilbert-Schmidt (HS) operator if c 2 (x) = c 1 (x) in a small neighborhood of ∂Ω. Thus, if q † represents the 'true' model and, hence, F (q † ) = S q † the data, we can introduce the mismatch functional (1.4) J (q) = 1 2 S q -S q † 2 using the HS norm

(1.5) S q -S q † = ∞ j=1 (S q -S q † )ψ j 2 1/2
, where {ψ j } is a countable basis in H 1/2 (∂Ω), and could be chosen to be right singular functions of the residual operator. Note that the norm of the residual operator, S q -S q † , does not depend on the choice of basis {ψ j }. As an alternative, the log criterion replacing (S q -S q † )ψ j 2 by log(S q ψ j ) -log(S q † ψ j ) 2 , has been proposed [START_REF]Inverse Problem Theory: methods for data fitting and and model parameter estimation[END_REF][START_REF] Shin | Waveform inversion using a logarithmic wavefield[END_REF]. This is well-defined considering S q -S q † as a densely defined operator and using a basis {ψ j } of H 3/2 (Ω). Intuitively, the log criterion may capture better the components of the residual which belong to the small singular values; however, it does not give a proper norm for the residual operator. We will consider both functionals in our examples, where we will truncate the sum at N , a fixed number of terms.

2. Iterative scheme.

2.1. Gradient descent. We consider time-harmonic waves, described by solutions, u say, of the Helmholtz equation on R n , n ≥ 3. We write q(x) = -ω 2 c -2 (x, ω); ω ∈ R. We introduce a bounded open (computational) domain Θ ⊂ R n with boundary in C (1,1) . We assume that Θ = R n \Θ is connected. We consider the problem, (2.1)

         (-∆ -k 2 )u = 0, x ∈ Θ , (-∆ + q(x))u = f, x ∈ Θ , lim r→∞ r (n-1)/2 ∂u ∂r -iku = 0 . Here, q ∈ L ∞ (Θ) while -k 2 = -ω 2 c -2 ∞
, where c ∞ is constant. We replace (2.1) by the equivalent problem on Θ,

(2.2)    (-∆ + q(x))u = f, x ∈ Θ , ∂u ∂ν = -Λ e u, x ∈ ∂Θ ,
where Λ e is the exterior Dirichlet-to-Neumann map for the Helmholtz equation on Θ . We introduce an open bounded domain, D, with D ⊂ Θ and a boundary in C (1,1) , and D = R n \D being connected. Let Σ be an open portion of ∂D. In the above, furthermore, f represents volume sources supported in ∂D; this property will be implicit through the introduction of an appropriate Dirac measure and factoring out ψ (defined on ∂D) in f . The solution, u, of (2.2), in H 1 (Θ), also solves the variational equation

(2.3) b(q; u, v) = s(v) for all v ∈ H 1 (Θ),
where

(2.4) b(q; u, v) = Θ (∇u • ∇v -q uv) dx, and
(2.5) s(v) = Θ f v dx - ∂Θ (Λ e u) v dσ.
We consider the one-parameter family of models, q t , with q 0 = q, and a set of source functions,

(f k (x)) N k=1 ; we collect (u k,t (x)) N k=1 satisfying (2.6) b(q t ; u k,t , v) = s k (v) for all v ∈ H 1 (Θ), with (2.7) s k (v) = Θ f k v dx - ∂Θ (Λ e u) v dσ.
(The first integral on the right-hand side is effectively an integral over ∂D.) We write the solution, u k,t , restricted to ∂D as an operator, S(q t ), acting on f k , and introduce the functional, (2.8)

J (q t ) = 1 2 N k=1 P (S(q t )f k -d f,k ) 2 L 2 (∂D) ;
here, (d f,k (x)) N k=1 represents the data, and P is an elliptic operator such that

P d L 2 (∂D) = d H 3/2 (∂D) .
Taking the derivative with respect to t yields

(2.9) d dt J (q t ) t=0 = N k=1 Re Σ [P * P(S(q t )f k -d f,k )] d dt S(q t )f k t=0 dσ. We introduce v = (v k (x)) N k=1 , Lagrange multipliers (w k (x)) N k=1
and a family of functionals, (2.10)

L(q t ; v, w) = J (q t ) + N k=1 Re (b(q t ; v k , w k ) -s k (w k )),
where

S(q t )f k in J (q t ) is replaced by v k . If v k = u k,t is the solution of the direct problem (2.6), L(q t ; u t , w) = J (q t ) for all w k ∈ H 1 (Θ), k = 1, . . . , N, then d dt L(q t ; u t , w) t=0 = N k=1 Re Σ [P * P(S(q t )f k -d f,k )] d dt S(q t )f k t=0 dσ +b q; d dt S(q t )f k t=0 , w k + d dt b(q t ; u k , w k ) t=0 for all w k ∈ H 1 (Θ). If w k solves b(q; v, w k ) = - Σ [P * P(S(q t )f k -d f,k )] v dσ for all v ∈ H 1 (Θ), k = 1, . . . , N, (2.11) then (2.12) d dt L(q t , u t , w) t=0 = N k=1 Re d dt b(q t ; u k , w k ) t=0 .
We can write (2.11) in the form of the direct problems:

b(q; wk , v) = - Σ [P * P(S(q t )f k -d f,k )] v dσ for all v ∈ H 1 (Θ), k = 1, . . . , N ; (2.13)
here, ( wk ) N k=1 is the so-called adjoint state, with wk ∈ H 1 (Θ) being the solution to

(2.14)          (-∆ -k 2 ) wk = 0, x ∈ Θ , (-∆ + q(x)) wk = [P * P(S(q t )f k -d f,k )] δ Σ , x ∈ Θ , lim r→∞ r (n-1)/2 ∂ wk ∂r -ik wk = 0 .
Again, we have the equivalent problems (2.15)

(-∆ + q(x)) wk = [P * P(S(q t )f k -d f,k )] δ Σ , x ∈ Θ , ∂ wk ∂ν = -Λ e wk , x ∈ ∂Θ .
We have made use of the fact that Λ * e ( wk ) = Λ e w k . We then evaluate

(2.16) d dt b(q t ; u k , w k ) t=0 = - Θ u k wk d dt q t t=0 dx.
We identify d dt q t t=0 with γ.

Thus

(2.17)

dJ (q, γ) = - N k=1 Re Θ u k wk γ dx =: (ρ, γ) L 2 (Θ) .
In the case of attenuation in d dt q t above, we can distinguish its real and imaginary parts. Instead, we distinguish between ∂q ∂cp d dt c p;t and ∂q ∂Q d dt Q t , where

                   ∂q ∂c p = - 2ω 2 c 3 p 1 - i sgn(ω) 2Q 2 , ∂q ∂Q = - i sgn(ω)ω 2 c 2 p Q 2 1 - i sgn(ω) 2Q 3 .
A comprehensive review for gradient computations may be found in [START_REF] Pratt | Gauss-newton and full newton methods in frequency-space seismic waveform inversion[END_REF][START_REF] Plessix | A review of the adjoint-state method for computing the gradient of a functional with geophysical applications[END_REF].

Convergence.

The convergence rate and convergence radius of our iterative scheme are derived based on a conditional Lipschitz-type stability estimate for the inverse problem. In [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation with multi-frequency data[END_REF] it was proven that (2.18) q -q † ≤ C F (q) -F (q † ) , in particular, for piecewise constant functions,

q(x) = J j=1 q j χ Dj (x),
incorporating interfaces, pinchouts etc., where the q j are constant coefficients and the χ characteristic functions; here, we pre-assume a domain partitioning {D j } J j=1 of Ω. Moreover, in [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation with multi-frequency data[END_REF] it was shown that C increases exponentially with increasing J. Increasing J implies an increasing spatial resolution.

The convergence radius is also determined by the lower bound of DF (q † ) , where DF signifies the Fréchet derivative of F . Indeed, such a bound implies the above mentioned Lipschitz stability estimate. In the case of weak attenuation, it is clear that the lower bound of the Fréchet derivative with respect to Q decreases significantly, which implies a rapidly increasing stability constant, with decreasing frequency.

We mitigate the above mentioned growth of the stability constant by introducing a multi-level approach based on a hierarchical compression of the model during the reconstruction. The compression at any level leads to an approximation error: Let q † J denote the compressed approximation of q † using J domains, then the error, F (q † J ) -F (q † ) , decreases as ω 2 as ω ↓ 0. A condition on the coupling between approximation errors and stability constants from level to level guarantees convergence of our scheme, see [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]. This condition is satisfied if q † J -q † decreases at least exponentially with J. Given this compression rate, we can then choose an appropriate set of frequencies which control the approximation errors at the different levels. A straightforward choice is facilitated by a gradual increase in J. For the analysis, precise conditions and expressions for the convergence radii, see [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation with multi-frequency data[END_REF].

In our implementation we use the nonlinear conjugate gradient method including projections. The step size is computed using a backtracking line search. (By successively increasing the step size one finds the smallest one that reduces the functional, which is then used to update the model.)

Adaptive domain partitioning and model compression.

In our applications, we will employ domain partitionings based on boxes and piecewise constant functions derived from the Haar basis. The boxes naturally fit a grid. Haar wavelets have been used for compressed model representations, for example, in tomography [START_REF] Loris | Tomographic inversion using l1-norm regularization of wavelet coefficients[END_REF][START_REF] Loris | Nonlinear regularization techniques for seismic tomography[END_REF]. In subsection 3.1 we briefly review the basis of Haar wavelets and the notion of scale refinement. The increase of the number of domains by refining the scale one level is large and hence an obstruction to a proper use of our multi-level approach. To arrive at a more gradual increase of the number of domains, still using Haar wavelets, we introduce in subsection 3.2 two adaptive strategies based on local refinement. We give numerical examples for a test model. We will apply the adaptive strategies to the gradient, periodically, in the iterative scheme; that is, we use the current gradient to determine, periodically, the next domain partitioning in our multi-level approach; see subsection 3.3. [START_REF] Haar | Zur theorie der orthogonalen funktionensysteme[END_REF] is defined using a scaling function φ and a mother wavelet ψ. In the one-dimensional case,

Haar wavelets. The Haar basis

ψ(x) =            1 if 0 ≤ x < 1 2 -1 if 1 2 ≤ x < 1 0 elsewhere.
and the basis is formed by functions

ψ j,k (x) = 2 -j/2 ψ(2 -j x -k), j, k ∈ Z. This basis is orthogonal in L 2 .
The scaling function is simply the characteristic function on the interval [0; 1], while

φ j,k (x) = |I j,k | -1/2 if x ∈ I j,k 0 elsewhere,
where

I j,k = [2 j (k -1), 2 j k].
The Haar coefficients of a function f are then given by d j k = f, ψ j,k , where ., . denotes the inner product in L 2 , and the 'averages' of f on I j,k by

s j k = f, φ j,k .
These satisfy the relationships,

d j+1 k = 1 √ 2 (s j 2k-1 -s j 2k )
and

s j+1 k = 1 √ 2 (s j 2k-1 + s j 2k ),
which are used in the design of a fast transform [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms i[END_REF]; [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. We follow an extension to higher dimensions based on taking tensor products.

In this section, we use a model example on a regular grid which we refer to as louro3d. The model contains irregularly shaped salt bodies, is of size 124 × 79 × 31, and represents a wavespeed. We will illustrate multi-scale approximations of the model by injecting it into the original grid, which corresponds with the 0 level compression. In Figure 1 we illustrate the decomposition of louro3d into Haar wavelets over 4 levels of compression. The number of coefficients for each level is given in Table 3.1. It is immediate from the structure of the Haar basis that refining the scale increases J in the domain partitioning by a factor 8. For the multi-level scheme, which requires a gradual increase of J, we now introduce an adaptive local refinement.

3.2. Local refinement. We discuss two methods, which assume a given coefficient model on a (computational) grid.

Method 1. In this method, we simply refine only where the deviation from the local average of the wavespeed evaluated through the scaling functions appears to be above a certain threshold.

In case a refinement is applied, the box (domain) is split into 8 new boxes. Furthermore, we can distinguish the directions and subject them to refinement independently.

We apply this method to the louro3d model. We illustrate it for compression level 2 highlighting the implied domain partitionings. In Figure 2 we show a horizontal slice of the results; the left figure corresponds with Figure 1 (level 2) while the middle figure represents the method discussed here. The number of domains or coefficients has been moderately increased, from 4123 (1.4% of the original) to 7455 (2.5% of the original), still far below the number for Haar level 1 (12% of the original). In Figure 4 we show the logarithmic approximation error in L 2 as a function of the number of domains, J, in the domain partitioning, and illustrate the effect of a more gradual increase (the red dots).

Method 2. The strategy of our second method is based on gluing boxes together forming larger domains. The criterion for gluing follows from comparing adjacent averages at a current level. If the averages are sufficiently close we glue to underlying domains together. The threshold in the criterion essentially controls the increase in the number of domains in the domain partitioning from level to level. We show an example using the the louro3d model in Figure 2 (right). Here, the level 1 Haar compression yields the initial paritioning. We apply the gluing in the horizontal directions only. In Figure 4 we show the logarithmic approximation error in L 2 as a function of the number of domains, J, in the domain partitioning, and illustrate the effect of a more gradual increase (the green dots). While the first method keeps the structure of the Haar basis, the second method, of course, provides more flexibility, while both methods can be applied directionally. We illustrate this in Figure 5 for a 2-dimensional grid. Method 1 starts from a coarser scale representation (left subfigure); a local refinement to a finer scale is illustrated in the two top subfigures. Method 2 starts from a finer scale scale representation (right subfigure); the gluing to domains corresponding with a coarser scale scale is illustrated in the two bottom subfigures.

Gradient adaptive compression.

In our multi-level iterative scheme, we refine the domain partitioning at each new frequency in the (discrete) frequency progression. In the refinement, we use the current gradient as a point of departure and project it onto two neighboring scales; the finest of these two scales corresponds with the next level. Method 1 uses the coarsest scale to construct the new domain partitioning while Method 2 uses the finest scale for this purpose.

4. Numerical experiments. We carry out numerical experiments illustrating the performance of our multi-level, multi-frequency iterative scheme using modelled data. We consider the data to be noise-free, though, in view of the mentioned stability estimates, additive noise through the approximation errors will not influence the results in a significant way. As mentioned before, the underlying theory including a condition for convergence can be found in [START_REF] Beretta | Inverse boundary value problem for the helmholtz equation with multi-frequency data[END_REF]. We invoke the gradient adaptive refinement of the domain partitioning. We start with the lowest available frequency and a model representation on the coarsest scale using the computational grid as a reference for the finest scale. After a prescribed number of iterations, the frequency is increased; we carry out a domain re-partitioning as outlined in the previous section, and then restart the iterations. This process progresses through a set of frequencies determined by the condition on the approximation errors which involves the upper bound of the stability constants, and the stability constants themselves. These implicitly depend on the 'true' solution which we do not know in practice. However, through some trial iterations, a reasonable choice of frequencies can be obtained. In subsections 4.1 and 4.2 the wavespeed is real-valued and the stability constants and approximation errors control the convergence of our multi-level, multi-frequency iterative scheme; in subsection 4.3 we consider complex wavespeeds and attenuation where the (frequency dependent) lower bound of the Fréchet derivative also strongly influences the convergence in particular for the recovery of the quality factor.

4.1. Salt bodies: louro3d. The louro3d model was illustrated in Figure 1 (bottom); the step size is 20m and the wavespeed varies between 1517 m/s and 4527 m/s. We generate data between 1 Hz and 15 Hz with our Hsolver.

We employ, here, the non-adaptive approach using the Haar basis for the domain partitioning and refinement. The initial model is shown Figures 6 (right) and 7 (right), and does not contain any indication of the salt bodies. We begin with illustrating the role of compression in the convergence of our scheme. This role becomes pronounced at relatively higher frequencies, as expected from the explicit expressions for the convergence rates [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]. In Figure 8, the mismatch functionals are plotted as a function of the iteration, at 2 Hz (left) and 3 Hz (right), with compression (level 2) and without compression (level 0). We proceed with applying our multi-level iterative scheme progressively for frequencies in the set {2, 3, 4, 5, 6, 7, 8, 9, 10} (in Hz). Here, we simply use the Haar basis for compression and domain partitioning. At 2 Hz we start with level 2; at 4 Hz we refine the domain partitioning to level 1 and at 8 Hz to level 0. Figure 9 shows a horizontal slice of the (partial) recovery at level 2 and Figure 10 shows the same horizontal slice of the (partial) recovery at level 1. In Figures 11 and12 we show the (partial) recovery obtained at level 0. The model consists of multiple disconnected salt bodies of different sizes. Even after the first iterations at the coarsest level, in Figure 9, the largest bodies appear with a wavespeed close to the true wavespeed at their centers. Following the progression in our multi-level, multi-frequency scheme we begin to resolve the smaller bodies, as well as the boundaries of the bodies and the wavespeeds in their interiors; see Figures 10,[START_REF] Bozdag | Misfit functions for full waveform inversions based on instantaneous phase and envelope measurements[END_REF] Figure 13 presents the residuals corresponding with a centrally located shot (real part) during the 5 Hz reconstructive iterations, Figure 14 during the 10 Hz iterations. The left figure shows the starting residual, in the middle after 5 iterations at the considered frequency, on right after 20 iterations. The residuals are decreasing while the salt bodies are partly being resolved. The residual does not decrease near the boundary of the computational domain; indeed, we cannot resolve the wavespeed near this boundary. The set of frequencies used in our scheme is {3, [START_REF]About the stability of the inverse problem in the 1-d wave equation[END_REF][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][START_REF] Beretta | Lipschitz stability of an inverse boundary value problem for a Schr\"{o}dinger type equation[END_REF][START_REF] Beretta | Inverse boundary value problem for the helmholtz equation with multi-frequency data[END_REF][START_REF] Berkhout | Changing the mindset in seismic data acquisition[END_REF][START_REF] Beylkin | Fast wavelet transforms and numerical algorithms i[END_REF][START_REF] Brossier | Imagerie sismique à deux dimensions des milieux viso-élastiques par inversion des formes d'ondes : développements méthodologiques et applications[END_REF][START_REF] Conn | Trust Region Methods[END_REF][START_REF] Futterman | Dispersive body waves[END_REF][START_REF] Ha | Acoustic full waveform inversion of synthetic land and marine data in the laplace domain[END_REF]25} (in Hz). We carry out a maximum of 20 iterations per frequency. The initial model is built from a one-dimensional wavespeed profile in depth. Nonlinear compressed reconstruction using the Haar basis. The Haar level 3 is the initial compression. It is switched to level 2 at 4 Hz, to level 1 at 6 Hz and to level 0 at 10 Hz. Figures 17 and18 show the finite reconstruction at 25 Hz. Figure 19 gives a three-dimensional view and also shows wavespeed profiles in depth. We note the wavespeed in the high wavespeed layer has not completely been resolved. We expect that an iterative scheme based on high-frequency data will further resolve such a detail. Nonlinear compressed reconstruction using method 1. Again, the Haar level 3 is the initial compression. We apply method 1 for the domain re-partitioning. Table 4.1 shows the implied gradual increase in the number of domains in the multi-level approach while increasing the frequency. Figures 20,22 (left) Nonlinear compressed reconstruction using method 2. Similarly we now apply method 2 for the domain re-partitioning. The number of domains in the multi-level approach is represented by the Table 4 The differences between the results obtained with the Haar basis and method 1 or 2 for domain re-partitioning seem to be small. However, differences show up in the mismatch functional; see Figure 24. The differences become less in the progression to higher frequencies for a given domain partitioning; we illustrate this for the final iterations at 20 -25 Hz. Method 1 gives the best result. The blue curve corresponds with the refinement using the Haar basis while the red curve corresponds with the refinement using method 1. The green curve corresponds with the refinement using method 2. The domain partitioning is kept fixed in both cases.

4.3. Attenuation. Here, we consider complex wavespeeds and attenuation. We study the partial recovery of the phase velocity and quality factor (cf. (1.2)), or the real and imaginary parts of the wavespeed. For our numerical experiments, we use a model containing an artificial geobody with the shape of a cross; see Figures 25 and26. Outside the body, the phase velocity increases linearly in depth; also, the phase velocity is significantly higher inside the body than outside it. The quality factor is chosen to increase with increasing phase velocity. The model is of size 124 × 79 × 31 with a step size of 20 m.

There is no fundamental obstruction to recover a complex wavespeed, and, indeed, our stability estimates hold for this case. However, the key challenge, here, is that the Fréchet derivative with respect to the quality factor becomes significant in magnitude only for relatively high frequencies. We encounter the effect of the (lower) bound of the Fréchet derivative on the radius of convergence and convergence rates. In our model, a sufficiently high frequency is reached at 15 Hz. At this frequency, partial recovery of the quality factor becomes possible; see Figure 28. In Figure 27 we show the partial recovery of the phase velocity. An appropriate set of frequencies follows to be {3, 4, 5, 6, 7, 8, 10, 12, 15} Hz. We begin with the Haar basis at level 2; at 4 Hz we refine the domain partitioning to level 1 and then to level 0 at 8 Hz. 5. Discussion. We have implemented and applied the analysis and construction of a multilevel, multi-frequency iterative scheme including conditions for convergence developed 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 analysis provided an understanding and a precise characterization of, for example, frequency progression in full waveform inversion; it also showed that this progression should be coupled to hierarchical model compression in a particular way to guarantee convergence. The analysis originated from a Lipschitz stability estimate for the inverse problem [START_REF] Beretta | Lipschitz stability of an inverse boundary value problem for a Schr\"{o}dinger type equation[END_REF] using the single layer potential operator as the data. This estimate was obtained for piecewise constant (blocky) models incorporating interfaces, pinchouts, etc., corresponding with a given domain partitioning. We introduced two adaptive methods of (directional) compression and compared their performance in computational experiments using two geologically different models. The compression overcomes in part the need for very-low frequency data while allowing initial models not containing prior information about the structure of the 'true' model. Our analysis and application are based on the behaviors of the stability constants, and approximation errors at low frequencies. The stability analysis, which will be local, and, hence, convergence, at high frequencies requires the use of different techniques, which is a current subject of research.

Fig. 1 .

 1 Fig. 1. Decomposition into Haar wavelets, and implied domain partitionings, from level 4 (upper left figure) to level 0 (bottom figure). The level 0 corresponds with the original louro3d model.

Fig. 2 .

 2 Fig. 2. A horizontal slice of the compressed louro3d model at level 2; the color bar indicates wavespeed in m/s. Left: Haar basis (4123 domains); middle: method 1 (7 455 domains); right: method 2 (9704 domains).

Fig. 3 .

 3 Fig. 3. Representation of the local refinement methods for compression and domain partitioning. Left: method 1 (7455 domains); right: method 2 (9704 domains).

Fig. 4 .

 4 Fig.4. The error, in L 2 norm, as a function of the number of domains in the domain partitioning. The blue dots correspond with the different levels following the compression with Haar wavelets. The red dots correspond with method 1, and the green dots correspond with method 2; these methods generate intermediate domain partitionings as compared with the one using the Haar wavelets.

Fig. 5 .

 5 Fig. 5. Example of the local refinement methods in 2D. Method 1 is illustrated in the upper part, while method 2 is illustrated in the bottom part. Left: the Haar level 1 represention; right: the Haar level 0 representation. Method 1 takes the coarser level as a starting point, selects which part should be refined and then refines as if it were Haar with one level down on the selected part (upper right subfigure). Method 2 starts from the finer level, selects which domains should be glued together (two in this example) and merges domains accordingly (lower left subfigure).

Fig. 6 .

 6 Fig. 6. Left: True wavespeed model; right: initial wavespeed model.

Fig. 7 .

 7 Fig. 7. Vertical slices of the true wavespeed model (left) and initial wavespeed model (right).

Fig. 8 .

 8 Fig. 8. The mismatch functionals (logarithmic scale) as a function of iteration at 2 Hz (left) and 3 Hz (right), with compression (level 2, the red curve) and without (level 0, the blue curve).

Fig. 9 .

 9 Fig. 9. A horizontal slice of the true wavespeed model (left) and the (partial) recovery at level 2 (right).

Fig. 10 .

 10 Fig. 10. A horizontal slice of the true wavespeed model (left) and the (partial) recovery at level 1 (right).

Fig. 11 .

 11 Fig. 11. The true wavespeed model (left) and the (partial) recovery at level 0 (right).

Fig. 12 .

 12 Fig. 12. A vertical slice of the true wavespeed model (left) and the (partial) recovery at level 0 (right).

  and 12.

Fig. 13 .

 13 Fig. 13. The residual for a central shot at different stages of the 5Hz frequency iterative reconstruction of the louro3d model. The color bar indicates the difference, all Figures have the same scale. Left: at the first iteration; middle: after the fifth iteration; right: after the twentieth iteration.

Fig. 14 .

 14 Fig. 14. The residual for a central shot at different stages of the 10Hz frequency iterative reconstruction of the louro3d model. The color bar indicates the difference, all Figures have the same scale. Left: at the first iteration; middle: after the fifth iteration; right: after the twentieth iteration.

4. 2 .

 2 Statoil model. We use a second test model (courtesy Statoil) representing a different geological environment with high wavespeed contrasts to illustrate our adaptive compression methods. The model is of size 256 × 146 × 123 and the stepsize of the grid is 10 m; see Figures 15 and 16

  .

Fig. 15 .

 15 Fig. 15. True wavespeed model (left) and initial model (right).

Fig. 16 .

 16 Fig. 16. A vertical slice of the true wavespeed model (left) and the initial model (right).

Fig. 17 .

 17 Fig. 17. Horizontal slice of the true wavespeed model (left) and the reconstruction at 25 Hz (right).

Fig. 18 .

 18 Fig. 18. Vertical slice of the true wavespeed model (left) and the reconstruction at 25 Hz (right).

Fig. 19 .

 19 Fig. 19. The partial reconstruction at 25 Hz (left) and wavespeed profiles in depth at x = 1.3km; y = 400m: the true wavespeed in blue, the initial model in green, and the reconstruction using Haar basis in red circles.

Fig. 20 .

 20 Fig.[START_REF] Ha | Acoustic full waveform inversion of synthetic land and marine data in the laplace domain[END_REF]. The partial reconstruction at 25 Hz using method 1 (left) and wavespeed profiles in depth at x = 1.3km; y = 400m: the true wavespeed in blue, the initial model in green, and the reconstruction in red circles using method 1.

Fig. 21 .

 21 Fig.[START_REF] Haar | Zur theorie der orthogonalen funktionensysteme[END_REF]. The partial reconstruction at 25 Hz using method 1 (left) and wavespeed profiles in depth at x = 1.3km; y = 400m: the true wavespeed in blue, the initial model in green, and the reconstruction in red circles using method 2.

Fig. 22 .

 22 Fig. 22. Vertical slice of the true wavespeed model (middle) and the reconstruction at 25Hz using method 1 (left) and method 2 (right).

Fig. 23 .

 23 Fig. 23. Horizontal slice of the true wavespeed model (middle) and the reconstruction at 25Hz using method 1 (left) and method 2 (right).

Fig. 24 .

 24 Fig.[START_REF] Lailly | The seismic inverse problem as a sequence of before stack migrations[END_REF]. Convergence rates: The mismatch functional (logarithmic scale) as a function of iterate. Left: at 20 Hz; right: at 25 Hz. The blue curve corresponds with the refinement using the Haar basis while the red curve corresponds with the refinement using method 1. The green curve corresponds with the refinement using method 2. The domain partitioning is kept fixed in both cases.

Fig. 25 .

 25 Fig. 25. Horizontal slice of the true phase velocity (left) and quality factor (middle right) and the initial model (middle left and right).

Fig. 26 .

 26 Fig. 26. Vertical slice of the true phase velocity (left) and quality factor (right) and the initial model (middle).

Fig. 27 .

 27 Fig. 27. Partial reconstruction of the phase velocity at 15 Hz.

Fig. 28 .

 28 Fig. 28. Partial reconstruction of the quality factor at 15 Hz.

Table 3 . 1

 31 Number of coefficients or domains in the domain partitioning depending on the level of compression for the louro3d model shown in Figure1.

  and 23 (left) show the finite reconstruction obtained with method 1.

	frequency level of compression number of domains percentage
	3Hz	Haar level 3	8 640	0.19%
	4Hz	Haar level 2	69 120	1.5 %
	5Hz	method 1	348 847	7.6 %
	6Hz	Haar level 1	569 984	12.4 %
	7Hz	method 1	1 242 950	27 %
	8Hz	method 1	2 090 685	45.5 %
	9Hz	method 1	3 052 765	66.4 %
	10Hz-25Hz	level 0	4 597 248	100 %

Table 4 . 1

 41 Number of domains versus frequency in the multi-level scheme for the model depicted in Figure 15 of original size 256 × 146 × 123.

Table 4 . 2

 42 .2. The finite reconstruction is shown in Figures 21, 22 (right) and 23 (right). Number of domains versus frequency in the multi-level scheme for the model depicted in Figure 15 of original size 256 × 146 × 123.

	frequency level of compression number of domains percentage
	3Hz	Haar level 3	8 640	0.19%
	4Hz	Haar level 2	69 120	1.5 %
	5Hz	method 2	469 224	10.2 %
	6Hz	Haar level 1	569 984	12.4 %
	7Hz	method 2	1 232 733	26.8 %
	8Hz	method 2	2 176 304	47 %
	9Hz	method 2	3 080 955	67 %
	10Hz-25Hz	level 0	4 597 248	100 %

* This research was supported by Total.