Hossein S Aghamiry
email: h.aghamiry@ut.ac.ir

Ali Gholami
email: agholami@ut.ac.ir.

Stéphane Operto
email: op-erto@geoazur.unice.fr.

Compound Regularization of Full-waveform

Inversion for Imaging Piecewise Media Hossein S. Aghamiry, Ali Gholami, Stéphane Operto, Abstract-Full waveform inversion (FWI) is an iterative nonlinear waveform matching procedure, which seeks to reconstruct unknown model parameters from partial waveform measurements. The nonlinear and ill-posed nature of FWI requires sophisticated regularization techniques to solve it. In most applications, the model parameters may be described by physical properties (e.g., wave speeds, density, attenuation, anisotropy) which are piecewise functions of space. Compound regularizations are thus beneficial to capture these different functions by FWI. We consider different implementations of compound regularizations in the wavefield reconstruction inversion (WRI) method, a formulation of FWI that extends its search space and prevents the socalled cycle skipping pathology. Our hybrid regularizations rely on Tikhonov and total variation (TV) functionals, from which we build two classes of hybrid regularizers: the first class is simply obtained by a convex combination (CC) of the two functionals, while the second relies on their infimal convolution (IC). In the former class, the model parameters are required to simultaneously satisfy different priors, while in the latter the model is broken into its basic components, each satisfying a distinct prior (e.g. smooth, piecewise constant, piecewise linear). We implement these compound regularizations in WRI using the alternating direction method of multipliers (ADMM). Then, we assess our regularized WRI for seismic imaging applications. Using a wide range of subsurface models, we conclude that the compound regularizer based on IC leads to the lowest error in the parameter reconstruction compared to that obtained with the CC counterpart and the Tikhonov and TV regularizers when used independently.

I. INTRODUCTION

F ULL waveform inversion (FWI) seeks to estimate consti- tutive parameters by nonlinear minimization of a distance between recorded and simulated wavefield measurements. This technology was originally developed in geophysical imaging [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF], and has spread more recently into other fields of imaging sciences such as medical imaging [START_REF] Sandhu | Frequency domain ultrasound waveform tomography: breast imaging using a ring transducer[END_REF] and oceanography [START_REF] Wood | Full waveform inversion of reflection seismic data for ocean temperature profiles[END_REF]. This partial-differential equation (PDE)-constrained nonlinear inverse problem is classically solved with local reduced-space optimization methods [START_REF] Haber | On optimization techniques for solving nonlinear inverse problems[END_REF]. In this linearized framework, a challenging source of non linearity is the so-called cycle skipping pathology which occurs when the initial model does not allow to match the data with a kinematic error smaller than half a period [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF], [START_REF] Virieux | An introduction to Full Waveform Inversion[END_REF]. Other sources of error are noise, approximate wave physics and ill-posedness resulting from parameter cross-talk, coarse acquisition sampling and uneven illumination of the targeted structure. Designing regularization techniques that mitigate these source of errors is therefore a key challenge for the success of FWI applications.

A proper regularization should be driven by the shape and statistical characteristics of the medium to be imaged. For example, in geophysical imaging, the subsurface can be represented by a piece-wise smooth medium, that is a model which contains smoothly varying and blocky components. The widespread Tikhonov regularizations [START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF] rely on the smoothness assumption and hence fail to recover sharp interfaces of such media. Conversely, first order total variation (TV) regularizations [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] are based on a blockiness assumption and hence are more suitable to image large contrasts. However, they generate undesirable staircase imprints in smooth regions [START_REF] Gholami | A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals[END_REF]. Hereafter we refer TV to as first order TV regularization, unless we explicitly mention the order of the TV regularization. Regions characterized by smoothly-varying properties and those containing sharp contrasts have different statistical properties. The former are characterized by the normal prior, while the latter by a heavy tailed prior [START_REF] Polson | Bayesian regularization: From Tikhonov to horseshoe[END_REF]. Consequently, simultaneous recovery of both properties is difficult when one type of regularization is used (Tikhonov, TV, etc). To overcome this issue, a combination of different regularizations can be used [START_REF] Gholami | Regularization of linear and nonlinear geophysical ill-posed problems with joint sparsity constraints[END_REF], [START_REF] Benning | Modern regularization methods for inverse problems[END_REF], [START_REF] Bioucas-Dias | An iterative algorithm for linear inverse problems with compound regularizers[END_REF]. A naive approach consists of the simple additive coupling or convex combinations (CC) of regularizations. Alternatively, [START_REF] Gholami | A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals[END_REF] proposed to explicitly decompose the model into several components of different statistical properties and use an appropriate regularization to reconstruct each component. Using this strategy, they combined Tikhonov and TV regularizations (referred to as TT regularization) to reconstruct piece-wise smooth media. The smooth components are captured by the Tikhonov regularization, while the blocky ones are determined by the TV counterpart. In many applications, it has been shown that a compound regularization based upon infimal convolution (IC) outperforms the one based upon additive coupling [START_REF] Bergmann | Priors with coupled first and second order differences for manifold-valued image processing[END_REF].

TT regularization based upon IC has been successfully applied to FWI for seismic subsurface imaging in the framework of iteratively-refined wavefield reconstruction inversion (IR-WRI) [START_REF] Aghamiry | Hybrid Tikhonov + total-variation regularization for imaging large-contrast media by fullwaveform inversion[END_REF]. IR-WRI extends the search space of FWI and decreases cycle skipping through a relaxation of the waveequation constraint [START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the search space[END_REF], [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF], [START_REF]Improving full-waveform inversion based on wavefield reconstruction via Bregman iterations[END_REF]. Taking advantage of the bilinearity of the wave equation, IR-WRI breaks down FWI into two linear subproblems which are solved in an alternating mode: wavefield reconstruction driven by the observables and model-parameter estimation by minimization of the source residuals the relaxation generated. Furthermore, [START_REF]ADMM-based multi-parameter wavefield reconstruction inversion in VTI acoustic media with TV regularization[END_REF] extended the method to acoustic multiparameter inversion. The linearity of the parameter-estimation subproblem provides a suitable framework to implement sophisticated nonsmooth regularizations.

In this study, following [START_REF] Gholami | A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals[END_REF] and [START_REF] Aghamiry | Hybrid Tikhonov + total-variation regularization for imaging large-contrast media by fullwaveform inversion[END_REF], we develop a general framework to combine a couple of regularization terms in IR-WRI through IC. Then, we specifically develop this framework for Tikhonov and TV regularizations, which are suitable for seismic subsurface imaging applications. Compared to [START_REF] Aghamiry | Hybrid Tikhonov + total-variation regularization for imaging large-contrast media by fullwaveform inversion[END_REF], we jointly update the blocky and the smooth components through a variable projection process rather than in an alternating mode. We first show that our new IC-based TT regularization outperforms the CC-based counterpart with several well-documented numerical benchmarks in the field of seismic imaging. We also compare the results obtained with these two TT regularizations with those obtained with total generalized variation (TGV) regularization, a combination of first and second order TV, and those obtained with Tikhonov and TV regularizations when used independently.

II. NOTATION

The mathematical symbols adopted in this paper are as follows. We use italics for scalar quantities, boldface lowercase letters for vectors, and boldface capital letters for matrices and tensors. We use the superscript T to denote the adjoint of an operator. The ith component of the column vector x is shown by x i and its absolute value is returned by |x i |. For the realvalued n-length column vectors x and y the dot product is defined by x, y = x T y and their Hadamard product, denoted by x • y, is another vector made up of their component-wise products, i.e. (x • y) i = x i y i . The 2 -and 1 -norms of x are, respectively, defined by

x 2 = x, x = n i=1 |x i | 2 and x 1 = n i=1 |x i |.

III. METHOD

In this section, we briefly review the frequency-domain FWI as a bi-convex feasibility problem and describe the extended forms of FWI. We show how the problem can be solved with the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction of multipliers[END_REF] for a general regularization function.

A. Full-waveform inversion

Frequency-domain FWI with a general regularization term and bounding constraints can be formulated as [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF], [START_REF] Aghamiry | Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media[END_REF]

min u,m∈C Φ(m) subject to A(m)u = b, Pu = d, (1)
where m ∈ R n×1 gathers unknown squared slowness, n is the number of discrete grid points, Φ(m) is an appropriate regularization term which we assume to be convex,

C = {x ∈ R n×1 | m l ≤ x ≤ m u }
is the set of all feasible models bounded by the lower bound m l and the upper bound m u . The first constraint in (1), A(m)u = b, is a partialdifferential equation (PDE) wherein u ∈ C n×1 is the wavefield and b ∈ C n×1 is the source term. In this study, A(m) ∈ C n×n is the discretized PDE Helmholtz operator [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF], [START_REF] Chen | An optimal 9-point finite difference scheme for the Helmholtz equation with PML[END_REF] given by

A(m) = ∆ + ω 2 C(m)diag(m)B, (2)
with ω the angular frequency and ∆ the discretized Laplace operator. The diagonal matrix C embeds boundary conditions and can be dependent or independent on m depending on the kinds of absorbing boundary conditions (radiation versus sponge) [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF]. Also, B is used to spread the mass term ω 2 C(m)diag(m) over all the coefficients of the stencil to improve its accuracy following an anti-lumped mass strategy [START_REF] Marfurt | Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[END_REF], [START_REF] Jo | An optimal 9-point, finite-difference, frequency-space 2D scalar extrapolator[END_REF], [START_REF] Hustedt | Mixed-grid and staggeredgrid finite difference methods for frequency domain acoustic wave modelling[END_REF].

The second constraint in (1), Pu = d, is the observation equation, in which d ∈ C m×1 is the recorded seismic data, m is the number of recorded data and P ∈ R m×n is a linear operator that samples the wavefield at the receiver positions.

1) WRI approach to solving (1): The extended approach, known as wavefield reconstruction inversion (WRI) [START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the search space[END_REF], recasts the constrained optimization problem, equation [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF], as an unconstrained problem where both constraints are implemented with quadratic penalty functions.

min u,m∈C Φ(m) + λ 0 2 Pu -d 2 2 + λ 1 2 A(m)u -b 2 2 , (3)
where λ 0 , λ 1 > 0 are the penalty parameters. For the unregularized case where Φ(m) = 0 and also without the bounding constraint, [START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the search space[END_REF] solved this biconvex minimization problem with an alternating-direction algorithm, whereby the joint minimization over u and m is replaced by an alternating minimization over each variable separately. The main property of the penalty formulation given by equation ( 3) is that the PDE constraint in the original problem is replaced by a quadratic penalty term, which enlarges the search space and mitigates the inversion nonlinearity accordingly [START_REF] Van Leeuwen | Mitigating local minima in full-waveform inversion by expanding the search space[END_REF]. Its main drawback, however, is the difficulty related to the adaptive tuning of the penalty parameter, which is common to all penalty methods [START_REF] Nocedal | Numerical Optimization[END_REF].

2) IR-WRI approach to solving (1): To overcome the above limitation, the iteratively-refined WRI (IR-WRI) implements the original constrained problem (1) with the augmented Lagrangian (AL) method [START_REF] Nocedal | Numerical Optimization[END_REF], [START_REF] Hestenes | Multiplier and gradient methods[END_REF].

min u,m∈C max v0,v1 Φ(m) + λ 0 2 Pu -d 2 2 + λ 1 2 A(m)u -b 2 2 + v T 0 [Pu -d] + v T 1 [A(m)u -b], (4)
where v 0 ∈ C m×1 and v 1 ∈ C n×1 are the dual variables (the Lagrangian multipliers). The min-max problem (4) can also be written in a more compact form (the scaled form AL) as

min u,m∈C max v0,v1 Φ(m) + λ 0 2 Pu -d + 1 λ 0 v 0 2 2 - λ 0 2 v 0 λ 0 2 2 + λ 1 2 A(m)u -b + 1 λ 1 v 1 2 2 - λ 1 2 v 1 λ 1 2 2 . (5)
Applying a gradient ascent to [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF] with respect to the duals, after a simple change of variables

d k = -v k 0 /λ 0 and b k = -v k 1 /λ 1 ,
gives the following iteration:

min u,m∈C Φ(m) + λ 0 2 Pu -d -d k 2 2 + λ 1 2 A(m)u -b -b k 2 2 , d k+1 = d k + d -Pu, b k+1 = b k + b -A(m)u, (6)
beginning with d 0 = 0 and b 0 = 0. Capitalizing on the bilinearity of the wave equation in m and u, ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction of multipliers[END_REF] is a powerful method to solve this kind of multivariate optimization problem. ADMM updates m and u separately through a Gauss-Seidel like iteration, i.e., fixing m and solving for u and vice versa. Accordingly, beginning with an initial guess m 0 , we end up with the following iteration to solve (6) [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF], [START_REF] Aghamiry | Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media[END_REF]:

u k+1 = arg min u λ0 λ1 P A(m k ) u - λ0 λ1 (d + d k ) b + b k 2 2 (7a) m k+1 = arg min m∈C Φ(m) + λ 1 2 A(m)u k+1 -b -b k 2 2 , (7b)
d k+1 = d k + d -Pu k+1 , (7c)
b k+1 = b k + b -A(m k+1 )u k+1 . (7d)
The subproblem (7a) associated with the wavefield reconstruction is quadratic and admits a closed-form solution. It relaxes the requirement to satisfy exactly the wave equation (A(m k )u = b) for the benefit of improved data fitting (Pu = d). This is achieved by reconstructing the wavefields that best jointly fit the observations and satisfy the wave equation in a least-squares sense. While wavefields generated by the reduced approach (Appendix A) satisfy exactly the wave equation, u r = A(m k ) -1 b, this makes classical FWI highly non-convex. Equation [START_REF] Virieux | An introduction to Full Waveform Inversion[END_REF] shows that the duals are updated with the running sum of the data and source residuals in iterations and are used to update the right-hand sides in the penalty functions of the scaled AL. These error correction terms in the AL method are the key ingredients that allow for a constant penalty parameter to be used in iterations, while guaranteeing convergence to accurate minimizer [START_REF] Nocedal | Numerical Optimization[END_REF]. In the next section we focus on the solution of the model subproblem (7b) when compound regularizations are used as the regularization term.

IV. THE MODEL SUBPROBLEM

This section presents compound regularization functionals and the details of our approach to solve the model subproblem (7b) with these functionals.

A. Compound regularizers

Simple regularizers, Appendix B, are effective for recovering models which can be characterized by a single prior and structure e.g. smooth, blocky, piecewise linear, etc. The compound regularizers are more effective for recovering complicated models that are represented by more than one prior. They are constructed by combining two or more separate simple regularizers. This can be done by either a convex combination (CC) or an infimal convolution (IC).

1) Convex combination of simple regularizers: A CC of r simple regularizer functionals Φ 1 , ..., Φ r is a compound regularizer functional of the form

Φ α (x) = α 1 Φ 1 (x) + ... + α r Φ r (x), (8)
where weights α i satisfy α i ≥ 0 and

α 1 + α 2 , ..., +α r = 1. (9)
Definitely, if all of the functions Φ 1 , ..., Φ r are convex then Φ is so. In CC models, the regularized solution is forced to satisfy the individual priors simultaneously. As an example, a compound regularizer functional constructed by a CC of 1and squared 2 -norms ( 1 + 2 ), known as an elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], [START_REF] Gholami | Sparse time-frequency decomposition and some applications[END_REF], is

Φ α (x) = (1 -α) x 2 2 + α x 1 , (10)
with 0 ≤ α ≤ 1. The convexity of 1 -and 2 -norms implies that Φ α (x) in ( 10) is convex. One may also construct a compound regularizer functional by a CC of two 1 -norms which are applied in different domains, such as those spanned by two different wavelet transforms, or those spanned by a wavelet transform and the gradient operator [START_REF] Gholami | Regularization of linear and nonlinear geophysical ill-posed problems with joint sparsity constraints[END_REF].

2) Infimal convolution of simple regularizers: In IC models, the solution is decomposed into simple components and then each component is regularized by an appropriate prior. Accordingly, the IC of r simple regularizer functionals Φ 1 , ..., Φ r is a compound functional of the form

Φ α (x) = min x=x1+...+xr {α 1 Φ 1 (x 1 ) + ... + α r Φ r (x r )}. ( 11
)
In the case of two functionals, Φ α in [START_REF] Benning | Modern regularization methods for inverse problems[END_REF] takes the form

Φ α (x) = min z {(1 -α)Φ 1 (x -z) + αΦ 2 (z)}, (12)
which is similar to the classical formula of convolution, and hence the term infimal convolution.

The IC of 1 -and (squared)

2 -norms ( 1 2 ) is Φ α (x) = min z {(1 -α) x -z 2 2 + α z 1 }, (13)
which is a denoising problem whose solution is unique and is given by the well-known soft-threshold function [START_REF] Donoho | De-noising by soft-thresholding[END_REF]:

z = max 1 - α 2(1 -α)|x| , 0 • x, (14)
Putting z from ( 14) into [START_REF] Gholami | A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals[END_REF] gives that

Φ α (x) = (1 -α)|x| 2 if |x| ≤ α 2(1-α) α|x| -α 2 4(1-α) if |x| > α 2(1-α) , (15)
which is nothing other than the Huber function [START_REF] Huber | Robust regression: Asymptotics, conjectures, and Monte Carlo[END_REF]. As seen, this function has a hybrid behavior: it has a quadratic behavior for small values of |x| and linear behavior for large values. The parameter α 2(1-α) determines where the transition from quadratic to linear behavior takes place.

Geometrical illustrations of the 1 -norm, 2 -norm, ( 1 + 2 )norm, and ( 1 2 )-norm for α = 0.7 are shown in Fig. 1. This figure shows that the 1 -and 2 -norms have a uniform behavior for all values, while the CC norm (the ( 1 + 2 )-norm) has a hybrid behavior: it approaches the 1 -norm near zero, where it behaves as a linear function, but approaches the 2 -norm for large values, where it behaves as a quadratic function. Unlike 1 + 2 , the IC function 1 2 approaches the 2 -norm near zero but is linear and approaches the 1 -norm for large values.

In this paper, we consider [START_REF] Bioucas-Dias | An iterative algorithm for linear inverse problems with compound regularizers[END_REF] in the following settings, though other configurations are possible:

Φ TT α (x) = min x=x1+x2 (1 -α) ∇ 2 x 2 2 2 + α ∇x 1 1 , (16)
and

Φ TGV α (x) = min x=x1+x2 (1 -α) ∇ 2 x 2 1 + α ∇x 1 1 , (17)
where, in both ( 16) and ( 17), the norms are applied on the absolute valued components of ∇ 2 x 2 (( 39) and ( 41)) and ∇x 1 [START_REF] Billette | The 2004 BP velocity benchmark[END_REF]. The compound regularizer Φ TT α is a combination of the second order Tikhonov and TV (TT) regularizations [START_REF] Gholami | A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals[END_REF] and Φ TGV α is a combination of the first and second order TV regularizations, called total generalized variation (TGV) [START_REF] Bredies | Total generalized variation[END_REF], [START_REF] Setzer | Infimal convolution regularizations with discrete 1-type functionals[END_REF]. The former is suitable for recovering piecewise-smooth models, while the latter is better suited for piecewise linear models. The next section gives a solution procedure to solve (7b) with these regularizers.

B. Solving the subproblem (7b)

In this section we present how to solve the subproblem (7b) with TT regularization. The solution procedure for the TGV regularizer follows easily. From the definition of A in (2), we get that

A(m)u = ∆u + Lm, (18)
where

L = ∂A(m) ∂m u = ω 2 Cdiag(Bu), (19)
and we assume that C does not depend on m (this is the case for perfectly-matched absorbing boundary conditions [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF]).

From the explicit decomposition m = m 1 + m 2 and ( 18), the solution of the optimization problem (7b) can be expressed as

arg min m=m1+m2 m∈C Φ TT α (m 1 , m 2 ) + λ 1 2 L[m 1 + m 2 ] -y 2 2 , (20)
where y = b + b k -∆u k+1 . Defining auxiliary variables p = ∇m 1 ∈ R 2n×1 and q = m 1 + m 2 ∈ R n×1 and recasting [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction of multipliers[END_REF] as a constrained problem and then applying ADMM leads to the following iteration [21, Section 2.2.2]

m k+1 1 m k+1 2 = arg min m1,m2 C(m 1 , m 2 , q k , qk , p k , pk ), (21a)
p k+1 = arg min p α p 1 + γ 1 2 ∇m k+1 1 -p -pk 2 2 , ( 21b
)
q k+1 = arg min q∈C γ 0 2 m k+1 1 + m k+1 2 -q -qk 2 2 , ( 21c
)
where γ 0 , γ 1 > 0 and

C(m 1 , m 2 , q k , qk , p k , pk ) = λ 1 2 L[m 1 + m 2 ] -y 2 2 + (1 -α) ∇ 2 m 2 2 2 + γ 1 2 ∇m 1 -p k -pk 2 2 + γ 0 2 m 1 + m 2 -q k -qk 2 2 .
The auxiliary primal variables p and q are introduced to decouple the 1 and the 2 minimization problems and solve the former ones with proximal algorithms following operator splitting methods. The dual variables p and q are updated through a gradient ascent step according to the method of multipliers [START_REF] Nocedal | Numerical Optimization[END_REF]

pk+1 = pk + p k+1 -∇m k+1 1 , (22a) qk+1 = qk + q k+1 -[m k+1 1 + m k+1 2 ], (22b)
We now discuss how to solve the subsubproblems given in [START_REF] Aghamiry | Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media[END_REF].

1) The subsubproblem (21a): A solution of subsubproblem (21a) occurs at the point where the derivatives of the objective function C with respect to m 1 and m 2 vanish simultaneously. Accordingly, we end up with the following linear system of equations:

G 11 G 12 G 21 G 22 m 1 m 2 = h 1 h 2 , ( 23
) with      G 11 = λ 1 L T L + γ 1 ∇ T ∇ + γ 0 I, G 12 = G 21 = λ 1 L T L + γ 0 I, G 22 = λ 1 L T L + (1 -α)(∇ 2 ) T ∇ 2 + γ 0 I, and
h 1 = λ 1 L T y + γ 1 ∇ T [p k + pk ] + γ 0 [q k + qk ], h 2 = λ 1 L T y + γ 0 [q k + qk ],
where I is the identity matrix. [START_REF] Aghamiry | Hybrid Tikhonov + total-variation regularization for imaging large-contrast media by fullwaveform inversion[END_REF] broke down the 2n × 2n problem (23) into two smaller n × n systems and updates m 1 and m 2 in alternating mode at the expense of convergence speed [15, their eqs. 10 and 11]. Instead, we solve here the original system exactly using a variable projection scheme, thus leading to faster convergence and more accurate results. From the first equation of ( 23), we find that

m 2 = G -1 12 [h 1 -G 11 m 1 ] (24)
and plugging this into the second equation of ( 23) we get the following:

m 1 = [G 11 -G 22 G -1 12 G 11 ] -1 [h 2 -G 22 G -1 12 h 1 ]. (25
) Interestingly, L is diagonal, implying that G 12 is also diagonal. Thus we only need to solve an n × n system to estimate m 1 , from which m 2 easily follows.

2) The subsubproblem (21b): The sub-problem for p, equation (21b), is a denoising problem and is straightforward to solve. Note that p has two components associated with the gradient in each direction:

p = p x p z . (26)
Equation ( 21b) is solved with a generalized proximity operator [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] leading to

p k+1 = prox α/γ1 (z) = ξ • z x ξ • z z , (27)
where

z = ∇m k+1 1 -pk = z x z z , (28)
and

ξ = max(1 - α γ 1 z 2 x + z 2 z , 0). ( 29
)
3) The subsubproblem (21c): The optimization problem (21c) also has an entrywise solution given by

q k+1 = proj C (m k+1 1 + m k+1 2 -qk ), (30)
where the projection operator projects its argument onto the desired box [m l , m u ] according to proj C (•) = min(max(•, m l ), m u ). Based on the above, the proposed ADMM based TT regularized IR-WRI algorithm is summarized in Algorithm 1. It should be noted that the total algorithm consists of two levels of iterations: an outer iteration (lines 4-13) given in [START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF] and an inner iteration (lines [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF][START_REF] Virieux | An introduction to Full Waveform Inversion[END_REF][START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Polson | Bayesian regularization: From Tikhonov to horseshoe[END_REF][START_REF] Gholami | Regularization of linear and nonlinear geophysical ill-posed problems with joint sparsity constraints[END_REF][START_REF] Benning | Modern regularization methods for inverse problems[END_REF] given in [START_REF] Aghamiry | Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media[END_REF] corresponding to the model subproblem (7b). Numerical results, however, show that only one inner iteration suffices for convergence of the algorithm, hence significantly reducing the total computational cost [START_REF] Goldstein | The split Bregman method for L1regularized problems[END_REF], [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF]. The main computational cost of the algorithm is on lines 4 and 5 where we need to solve the augmented PDE and the n × n sparse system given in equation 25. The other steps don't introduce significant computational overheads. The computational overhead introduced by compound regularizer compared to single regularizer results from the larger numerical bandwidth of the matrix

[G 11 -G 22 G -1 12 G 11 ]
in equation 25 compared to that associated with a single regularizer, either G 11 or G 22 .

Algorithm 1 ADMM based TT regularized IR-WRI algorithm 1: Initialize: set the dual variables b 0 , d 0 , p0 and q0 equal to 0 2: Input: m 0 (initial model parameters) 3: while convergence criteria not satisfied do 4: qk+1 ← qk + q k+1 -m k+1 12:

u k+1 ← λ 0 λ 1 P T P + A(m k ) T A(m k ) -1 λ 0 λ 1 P T [d + d k ] + A(m k ) T [b + b k ] 5: m k+1
p k+1 ← prox α/γ 1 (∇m k+1 1 -pk ) 9: q k+1 ← proj C (m k+1 -qk ) 10: pk+1 ← pk + p k+1 -∇m k+1
d k+1 ← d k + d -Pu k+1 13: b k+1 ← b k + b -A(m k+1 )u k+1 14: end while V. NUMERICAL EXAMPLES
We assess the performance of our algorithm against 1D and 2D mono-parameter synthetic examples. In Table I we give different regularization functions which are applied for stabilizing the FWI solution. We start with zero-offset Vertical-Seismic-Profiling (VSP) examples (1D IR-WRI) where the targeted wave speed profiles are selected from well-documented 2D benchmark subsurface velocity models in exploration seismic. To tackle more realistic applications, we proceed with a target of the 2D challenging 2004 BP salt model [START_REF] Billette | The 2004 BP velocity benchmark[END_REF] with noiseless and noisy data when a crude initial model and realistic frequencies are used as starting points.

A. Performance comparison using 1D test on benchmark models

First, we assess the performance of our regularized IR-WRI against 1D mono-parameter synthetic examples when the true models are 100 vertical profiles selected from the 2004 BP salt [START_REF] Billette | The 2004 BP velocity benchmark[END_REF], Marmousi II [START_REF] Martin | Marmousi2: An elastic upgrade for Marmousi[END_REF], SEG/EAGE overthrust [START_REF] Aminzadeh | 3-D Salt and Overthrust models[END_REF], SEG/EAGE salt [START_REF] Aminzadeh | 3-D Salt and Overthrust models[END_REF] and synthetic Valhall [START_REF] Prieux | On the footprint of anisotropy on isotropic full waveform inversion: the Valhall case study[END_REF] benchmark velocity models (we extracted 20 profiles from each benchmark model). For all of the experiments, a single source is used at the surface and the receivers are evenly deployed along the entire profile. A single frequency, whose value is set so that the reduced-space inversion is prone to cycle skipping, is considered for inversion. The model dimension, the inverted frequency and the receiver spacing are outlined for each model in Table II. We perform forward modeling with a 3point O(∆x 2 ) staggered-grid finite-difference stencil and PML absorbing boundary conditions at the two ends of the model. The starting model for IR-WRI is a homogeneous velocity model in which the velocity is the mean value of each profile. We set the penalty parameters according to the guideline given in Appendix C. Moreover, for a fair comparison of the compound regularizers (JTT, TT and TGV), we select for each of them the optimum value of α among a range of preset values that minimizes the error in the models estimated by the IR-WRI. Also, we set the parameter bounds m l and m u equal to 50% and 150% of the minimum and maximum velocities of the true model, respectively. The monochromatic inversion is performed with noiseless data when a maximum number of iterations, equal to 100, is used as a stopping criterion. The average error of the estimated velocity profiles for the five benchmark models and the different regularizations are plotted in Fig. 2. In this paper, the model error is defined as the energy of the difference between the true model and the estimated one compared to the energy of the true model. The errors in each model for different regularizations are normalized to 1 for sake of clarity (the error of DMP regularizer is not shown because of its worse performance). Fig. 2 clearly shows that the compound regularizations based upon infimal convolution (TT and TGV) always behave better than the CC regularization and the single regularization functionals (TV and Tikhonov). To emphasize the effects of the different regularization functions, we plot some close-ups of the reconstructed profiles in Fig. 3. These results show that TT provides the most accurate reconstruction for the 2004 BP salt (Fig. 3a) and Overthrust (Fig. 3c) models. This is consistent with the fact that the velocity trends of these two models match well the piecewise smooth prior. In contrast, TGV behaves slightly better than TT for the Valhall model, whose velocity trend is the closest one to the piecewise linear prior (Fig. 3d). For Marmousi II (Fig. 3b), TT and TGV give similar results.

DMP m 2 2 Tikhonov ∇ 2 m 2 2 TV ∇m 1 JTT (1 -α) ∇ 2 m 2 2 + α ∇m 1 TT min m=m 1 +m 2 {(1 -α) ∇ 2 m 2 2 2 + α ∇m 1 1 } TGV min m=m 1 +m 2 {(1 -α) ∇ 2 m 2 1 + α ∇m 1 1 }

B. 2004 BP salt model

We now consider a more realistic application with a target of the challenging 2004 BP salt model [START_REF] Billette | The 2004 BP velocity benchmark[END_REF]. The 2004 BP salt model is representative of the geology of the deep offshore Gulf of Mexico and mainly consists of a simple background with a complex rugose multi-valued salt body, sub-salt slow velocity anomalies related to over-pressure zones and a fast velocity anomaly to the right of the salt body. The selected subsurface model is 16250 m wide and 5825 m deep, and is discretized with a 25 m grid interval (Fig. 4a). We used 108 sources spaced 150 m apart on the top side of the model. We perform forward modeling with a staggered-grid 9-point finite-difference method [START_REF] Chen | An optimal 9-point finite difference scheme for the Helmholtz equation with PML[END_REF] with PML boundary conditions along the four edges of the model. The source signature is a 10 Hz Ricker wavelet. A line of receivers with a 25 m spacing are deployed all along the surface leading to a longoffset stationary-receiver acquisition. We used small batches of two frequencies with one frequency overlap between two consecutive batches, moving from the low frequencies to the higher ones according to a classical frequency continuation strategy. We use α = 0.7 and set the rest of hyperparameters according to the guidelines reviewed in appendix C. We set the parameter bounds m l and m u equal to the minimum and maximum velocity of the true model, respectively. The starting and final frequencies are 3 Hz and 13 Hz and the sampling interval in one batch is 0.5 Hz. The initial velocity model is a crude laterally-homogeneous velocity-gradient model with velocities ranging between 1.5 to 4.5 km/s (Fig. 4b). We start with inverting the first batch of frequencies ({3, 3.5} Hz) with noiseless data using a maximum number of iterations equal to 45 as a stopping criterion. To highlight the specific role of bound constraints, we activate them after 20 iterations. To emphasize the effect of regularization, the result of bound constrained IR-WRI with a simple DMP regularization is shown in Fig. 5a, while the bound-constrained IR-WRI results with Tikhonov and TV regularizations are shown in Figs. 5b and5c, respectively. Although the TV reconstruction is better than the Tikhonov one, it provides a velocity model which is far from the optimal one. A direct comparison between the true model, the starting model and the estimated models is shown in Fig. 6a along three vertical logs at 2.5, 9.0 and 15.0 km distance (as depicted with dashed white lines in Fig. 4a). We continue with compound regularization results which are shown in Fig. 5d-f and Fig. 6b. Clearly, the TT regularizer better captures the long wavelengths of the salt body and the smooth subsalt 7a-b, and the relative model errors in iterations, Fig. 7c-d, further confirm the relative performance of each regularizer during the inversion of the first frequency batch. Note the complex zigzag path followed by the inversion to jointly minimize the data residuals and the wave equation error in Fig. 7a-b. As already highlighted by [START_REF] Aghamiry | Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers[END_REF], this results from the dynamic balancing in iterations of the observationequation and wave-equation constraints performed by the dual updates with the data and source residuals.

We continue the inversion at higher frequencies using the final models of the {3, 3.5} Hz inversion, Fig. 5a-f, as initial models when the stopping criteria is either k max = 15 or

A(m k+1 )u k+1 -b 2 ≤ ε b , Pu k+1 -d 2 ≤ ε d , (31)
where k max denotes the maximum iteration count, ε b =1e-3, and ε d =1e-5. We perform three paths through the frequency batches to improve the IR-WRI results, using the final model of one path as the initial model of the next one (these cycles can be viewed as outer iterations of IR-WRI). The starting and finishing frequencies of the paths are [3.5, 6], [4, 8.5], [START_REF] Virieux | An introduction to Full Waveform Inversion[END_REF][START_REF] Gholami | A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals[END_REF] Hz respectively, where the first element of each pair shows the starting frequency and the second one is the finishing frequency. The bound-constrained IR-WRI models obtained from noiseless data are shown in Fig. 8. As for the inversion of the first batch, direct comparison between the true model, the starting model and the estimated models are shown in Fig. 9 along three vertical logs at 2.5 km, 9.0 km and 15 km distance (vertical dashed lines in Fig. 4a). The TT and TGV regularizers lead to high-quality velocity models, that capture both the fine-scale structure of the rugose largecontrast salt body and the high-velocity shallow anomaly on the right, as well as the smoother sub-salt background model including the low-velocity over-pressure structure. It is also Fig. 6: Noiseless 2004 BP salt case study. Direct comparison along the logs shown in Fig. 4a between the true velocity model (black), the initial model (dashed line) and the estimated models obtained after the {3, 3.5} Hz inversion (Fig. 5) with (a) simple regularizations (DMP in olive-green, Tikhonov in blue and TV in red) and (b) compound regularizations (JTT in orange, TT in green and TGV in pink).

worth noting the significant differences between the JTT and TT IR-WRI models in particular in the deep part of the model. Moreover, the number of iterations performed by IR-WRI for each regularization shows that TT has the best convergence speed (Table III).

As a final quality control of the different IR-WRI models, it is instructive to check the wave-equation and data residuals left by the different regularization methods for the starting 3-Hz frequency (Fig. 10 and11). The real part of waveequation error (Fig. 10) and data residuals (Fig. 11) are plotted at the first and final iterations of the inversion. Both of the final data and source residuals suggest that the TT regularizer slightly outperforms the TGV counterpart at low frequencies. To further illustrate the ability of compound regularizations to manage the blocky and smooth components of the subsurface, we show separately the two model components of IC based regularizations (m 1 and m 2 ) estimated with TT (Fig. 12a-b) and TGV regularizations (Fig. 12c-d) as well as histograms of ∇m 1 and ∇ 2 m 2 (Fig. 13). The trend of the reconstructed m 1 and m 2 shows that the IC based compound regularizations mostly succeeded in de-coupling the reconstruction of the blocky structure from that of the smooth background for TT regularization (Fig. 12a-b) and de-coupling the reconstruction of the blocky structure from that of the piecewise linear back- ground for TGV regularization (Fig. 12c-d). This statement is further supported by the long-tail shape and the Gaussian shape of the histograms of the ∇m 1 and ∇ 2 m 2 components of TT regularization (Fig. 13a-b) and the long-tail shapes of the ∇m 1 and ∇ 2 m 2 components of TGV regularization (Fig. 13c-d).

We continue by assessing the resilience of the different regularization strategies to noise when data are contaminated with a Gaussian random noise with a SNR=10 db, where SNR is defined as

SNR = 20 log A signal A noise , (32)
in which A denotes root mean square (RMS) amplitude. We use the same setup and the same initial velocity model (Fig. 4b) as those used for the noiseless case. The stopping criterion is defined by [START_REF] Donoho | De-noising by soft-thresholding[END_REF], where ε d is now set to the noise level.

The final models of bound-constrained IR-WRI obtained from noisy data are shown in Fig. 14. The number of iterations performed by IR-WRI with the different regularizations are outlined in Table III. As for the inversion of the first batch, results further confirm that TT regularization provides the most reliable results and illustrate the resilience of this compound regularization to noise.

VI. CONCLUSIONS

In this study, we first show how to efficiently implement different kinds of regularization and bound constraints in the wavefield reconstruction inversion method with the alternating direction method of multipliers (ADMM). Then, we show the capability of IR-WRI when equipped with compound Tikhonov and TV regularizations to accurately reconstruct large-contrast subsurface media when starting from a very crude initial model. This compound regularization is suitable for seismic imaging of the subsurface as it can often be represented by piecewise smooth media. We show that the infimal convolution (IC) of the Tikhonov and TV regularizers captures much more accurately the blocky and smooth com- ponents of the subsurface than the convex combination of the two regularizers. It also outperforms the Tikhonov and TV regularizers when used alone. We also show how the infimalconvolution regularizer can be efficiently implemented by jointly updating the smooth and blocky subsurface components through variable projection. Alternatively, TGV regularized IR-WRI can be a suitable tool to reconstruct piecewise linear media and provides similar results than TT IR-WRI. We conclude that such hybrid regularizations in the extended searchspace IR-WRI potentially provide a suitable framework to reconstruct, without cycle skipping, large-contrast subsurface media from ultra-long offset seismic data. It should also find applications in other fields of imaging sciences such as medical imaging.

VII. ACKNOWLEDGMENTS

We would like to thank the Editor Mauricio D. Sacchi, Gian Matharu and two anonymous reviewers for their comments which help improving the manuscript. This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by AKERBP, CGG, CHEVRON, EQUINOR, EXXON-MOBIL, JGI, PETRO-BRAS, SCHLUMBERGER, SHELL, SINOPEC and TOTAL. This study was granted access to the HPC resources of SIGAMM infrastructure (http://crimson.oca.eu), hosted by Observatoire de la Côte d'Azur and which is supported by the Provence-Alpes Côte d'Azur region, and the HPC resources of CINES/IDRIS/TGCC under the allocation A0050410596 made by GENCI." The reduced approach, which is more commonly used for sake of computational efficiency, strictly enforces the PDE constraint at each iteration by projection of the full space onto the parameter search space, leading to the following unconstrained optimization problem [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF], [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]

min m∈C Φ(m) + λ 0 2 ||PA -1 (m)b -d|| 2 2 , (33)
where λ 0 > 0 is the penalty parameter. A number of methods have been proposed to solve the optimization problems of the form [START_REF] Bredies | Total generalized variation[END_REF], either for unregularized form Φ(m) = 0 [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF], or the regularized form [START_REF] Asnaashari | Regularized seismic full waveform inversion with prior model information[END_REF], [START_REF] Esser | Total variation regularization strategies in Full-Waveform Inversion[END_REF]. Although this reduced approach is more computationally tractable than the full-space approach, the highly-oscillating nature of the inverse PDE operator A -1 makes the inverse problem highly nonlinear, and hence prone to convergence to a spurious local minima when the initial m is not accurate enough [START_REF] Virieux | An overview of full waveform inversion in exploration geophysics[END_REF], [START_REF] Symes | Migration velocity analysis and waveform inversion[END_REF]. The extended approach described in this paper (section III) is an alternative way which is more immune to local minima.

B. Simple regularizers

The two most widely used regularizers rely on the (squared)

2 and 1 -norms. The squared 2 -norm, defined as

x 2 2 = n i=1 |x i | 2 , (34)
promotes smooth reconstruction, since the minimization of the squared value of components will penalize large components more severely than small ones.

In contrast, the 1 -norm, defined as

x 1 = n i=1 |x i |, (35)
promotes sparse reconstruction (with many zero components), since the minimization of the absolute value of components will penalize small components more severely than the large counterparts.

The priors can be defined under a suitable transformation. For example, one may minimize the 1 -or 2 -norms of the first and/or second order differences of the model. The first order forward differences for discrete scalar field f in xand z-direction are denoted by ∇ x f and ∇ z f , with

(∇ x f ) i,j = f i,j -f i,j-1 , (∇ z f ) i,j = f i,j -f i-1,j , (36)
with appropriate boundary conditions, where i and j run over the domain of the model parameters. Accordingly, the discrete first order operator in 2D is defined as

∇ = ∇ T x ∇ T z T with (|∇f |) i,j = (∇ x f ) 2 i,j + (∇ z f ) 2 i,j . (37)
The squared 2 -norm of |∇f | gives the first order Tikhonov regularization [START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF], which returns a flat regularized model (with a small gradient), while its 1 -norm gives the total variation regularization [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], which returns a piecewise constant model (with a sparse gradient). Analogously, the second order forward differences are denoted by ∇ xx f and ∇ zz f , with

(∇ xx f ) i,j = f i,j-1 -2f i,j + f i,j+1 , (∇ zz f ) i,j = f i-1,j -2f i,j + f i+1,j , (38)
with appropriate boundary conditions, where again i and j run over the domain of the model parameters. Accordingly, the discrete second order operator is defined as

∇ 2 = ∇ T xx ∇ T zz T with (|∇ 2 f |) i,j = (∇ xx f ) 2 i,j + (∇ zz f ) 2 i,j . (39)
The squared 2 -norm of |∇ 2 f | gives the second-order Tikhonov regularization, which returns a smooth regularized model (with a small Laplacian), while its 1 -norm gives the second order TV regularization, which returns a piecewise linear model (with a sparse Laplacian).

Mixed second-order differences can also be constructed as

∇ xz f ≡ ∇ z ∇ x f with (∇ xz f ) i,j = f i,j -f i,j-1 -f i-1,j + f i-1,j-1 .
(40) ferences is defined as

∇ 2 = ∇ T xx √ 2∇ T xz ∇ T zz T with
(|∇ 2 f |) i,j = (∇ xx f ) 2 i,j + 2(∇ xy f ) 2 i,j + (∇ yy f ) 2 i,j , (41) which equals the Frobenius norm of the Hessian matrix [START_REF] Lefkimmiatis | Hessian-based norm regularization for image restoration with biomedical applications[END_REF], [START_REF] Gholami | Fast 3D blind seismic deconvolution via constrained total variation and gcv[END_REF], [START_REF] Gholami | 3D Dix inversion using boundconstrained TV regularization[END_REF].

C. Parameter tuning

Here, we provide some guidelines to tune the different hyperparameters in the regularized IR-WRI method. The reader is also referred to [START_REF] Aghamiry | Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media[END_REF] for more details. We start with γ 1 (step 8 of Algorithm 1), which controls the soft thresholding performed by the TV regularization, and set it equal to 0.02 × max |∇m k 1 -pk |. This tuning can be refined according to prior knowledge of the geological structure, coming from well logs for example. Also, we use the same weight for the bound constraints and the TV regularization: γ 0 = γ 1 . Once we set γ 1 , we define λ 1 such that γ 1 /λ 1 is a percentage of mean absolute value of the diagonal coefficients of L T L. Parameter λ 1 may be increased during iterations to reduce the weight of TV regularization and bound constraints near the convergence point. Finally, we set λ 0 such that λ = λ 1 /λ 0 is a small fraction of the highest eigenvalue ξ of the normal operator A(m) -T P T PA(m) -1 during the wavefield reconstruction subproblem according to the criterion proposed by [START_REF] Van Leeuwen | A penalty method for PDEconstrained optimization in inverse problems[END_REF]. In all the numerical tests, we use λ = 1e-5ξ and λ = 1e-3ξ for noiseless and noisy data, respectively.

Fig. 1 :

 1 Fig. 1: Geometrical illustration of different regularizers. (a) the 1 -norm, (b) the 2 -norm, (c) the ( 1 + 2 )-norm, and (d) the ( 1 2 )-norm.

1 ←

 1 update according to eq.[START_REF] Jo | An optimal 9-point, finite-difference, frequency-space 2D scalar extrapolator[END_REF]

Fig. 2 :

 2 Fig. 2: Zero offset VSP test. Average model error in estimated 1D profiles for different velocity models and different regularization functions.

Fig. 3 :

 3 Fig. 3: Zero offset VSP test. Some part of (a) 2004 BP salt, (b) Marmousi II, (c) Overthrust, and (d) synthetic Valhall models estimated with different regularizations.

Fig. 4 :Fig. 5 :

 45 Fig. 4: 2004 BP salt case study. (a) True velocity model. The vertical dashed lines indicate the location of vertical logs of Figs. 6, 9 and 15. (b) The velocity-gradient initial model.

Fig. 7 :

 7 Fig. 7: Noiseless 2004 BP salt case study. {3, 3.5} Hz inversion with the velocity-gradient initial model. (a-b) convergence history of the algorithm in the ( Pu kd 2 -A(m k )u kb 2 ) plane for (a) simple regularizations and (b) compound regularizations. The black arrow points the starting point. (c-d) evaluation of m km * 2 / m * 2 during the iteration where m * is the true model. The panels (a) and (b) as well as (c) and (d) are plotted with the same horizontal and vertical scale.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Noiseless 2004 BP salt case study. Final inversion results with the velocity models of Fig. 5a-f as initial models. The display of the panels is the same as that of Fig. 5.

Fig. 10 :

 10 Fig. 10: Noiseless 2004 BP salt case study. For a source at x=8.12 km and the 3-Hz frequency: (a) Real part of waveequation residual at first iteration of wavefield reconstruction, namely (A(m 0 )u 1 -b). (b-g) Real part of wave-equation residual at the final iteration achieved respectively with DMP, Tikhonov, TV, JTT, TT, and TGV regularization.

Fig. 11 :

 11 Fig. 11: Noiseless 2004 BP salt case study. For a source at x=8.12 km and the 3-Hz frequency: (a) Real part of data residual (Pu -d) at first iteration. (b-g) Real part of data residual at the final iteration achieved respectively with DMP, Tikhonov, TV, JTT, TT, and TGV regularization.

Fig. 12 :Fig. 13 :

 1213 Fig. 12: Noiseless 2004 BP salt case study. Final blocky and smooth components, m 1 and m 2 , reconstructed by IR-WRI with (a-b) TT and (c-d) TGV regularizations. The corresponding velocity models are shown in Figs. 8e and 8f. (a,c) m 1 . (b,d) m 2 . Note that m 1 and m 2 are parametrized with squared slownesses.

Fig. 14 :

 14 Fig. 14: Noisy 2004 BP salt case study. Final inversion results with Fig. 4b as initial model. (a-f) bound constrained IR-WRI with (a) DMP, (b) Tikhonov, (c) TV, (d) JTT, (e) TT and (f) TGV regularization.

Fig. 15 :

 15 Fig. 15: Noisy 2004 BP salt case study. Direct comparison (along the logs shown in Fig. 4a) between the true velocity model (black), the initial model (dashed line) and the estimated models with (a) simple regularizations (DMP in olive-green, Tikhonov in blue and TV in red) and (b) compound regularizations (JTT in orange, TT in green and TGV in pink).

TABLE I :

 I Different regularization functions.

	Abbreviation	Expression of Φ(m)

TABLE II :

 II Experimental setup of 1D model tests

		Length (km)	Inverted frequency (Hz)	Grid interval (m)	Receiver interval (m)
	2004 BP salt	11.46	5	6	180
	Marmousi II	3.75	12	5	85
	Overthrust	4.6	12	20	120
	SEG/EAGE salt model	4.2	10	20	120
	Synthetic valhall	5.22	5	25	175

TABLE III :

 III Total number of IR-WRI iterations for each regularization.

		DMP	Tikhonov	TV	JTT	TT	TGV
	Noiseless data	426	448	399	415	361	394
	Noisy data	285	293	264	271	270	274