Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media

TL;DR

This paper introduces a robust regularized full waveform inversion method using bound constraints and total variation within an ADMM framework, effectively handling cycle skipping and noise in high-resolution subsurface imaging.

Contribution

It develops a novel implementation of bound constraints and TV regularization in frequency-domain WRI using ADMM, enhancing robustness and resolution in complex media.

Findings

Successfully reconstructs salt and subsalt structures with high fidelity.

Resilient to cycle skipping and noise in large-contrast models.

Achieves high-resolution imaging starting from crude initial models.

Abstract

Full waveform inversion (FWI) is a waveform matching procedure, which can provide a subsurface model with a wavelength-scale resolution. However, this high resolution makes FWI prone to cycle skipping, which drives the inversion to a local minimum when the initial model is not accurate enough. Other sources of nonlinearities and ill-posedness are noise, uneven illumination, approximate wave physics and parameter cross-talks. All these sources of error require robust and versatile regularized optimization approaches to mitigate their imprint on FWI while preserving its intrinsic resolution power. To achieve this goal, we implement bound constraints and total variation (TV) regularization in the so-called frequency-domain wavefield-reconstruction inversion (WRI) with the alternating direction method of multipliers (ADMM). In the ADMM framework, WRI relies on an augmented Lagrangian function, a combination of penalty and Lagrangian functions, to extend the FWI search space by relaxing the wave-equation constraint during early iterations. Moreover, ADMM breaks down the joint wavefield reconstruction plus parameter estimation problem into a sequence of two linear subproblems, whose solutions are coordinated to provide the solution of the global problem. The decomposability of ADMM is further exploited to interface in a straightforward way bound constraints and TV regularization with WRI via variable splitting and proximal operators. The resilience of our regularized WRI formulation to cycle skipping and noise as well as its resolution power are illustrated with two targets of the large-contrast BP salt model. Starting from a 3Hz frequency and a crude initial model, the extended search space allows for the reconstruction of the salt and subsalt structures with a high fidelity.

Figures

Equations

$$\text{Find}m\text{and}u$$

$$\text{subject to}\mathbb{F}(m)u=s$$

$$s=\left[\begin{array}{c}d\\ b\end{array}\right],\mathbb{F}(m)=\left[\begin{array}{c}\mathbb{P}\\ \mathbb{A}(m)\end{array}\right],$$

$$\mathbb{A}(m)$$

$$\underset{u,m}{\min }0$$

$$\text{subject to}\mathbb{F}(m)u=s.$$

$$\underset{u,m}{\min }0+\Vert \mathbb{F}(m)u-s{\Vert }_{\Lambda }^2,$$

$$\mathbb{\Lambda}=\left[\begin{array}[]{r@{}c|c@{}l}&\begin{smallmatrix}\lambda_{0}&&0\\ &\ddots&\\ 0&&\lambda_{0}\rule[-4.30554pt]{0.0pt}{8.61108pt}\end{smallmatrix}&\mbox{0}&\hbox to0.0pt{\kern 14.22636pt$M$\hss}\\ \hline\cr&\mbox{0}&\begin{smallmatrix}\rule{0.0pt}{8.61108pt}\lambda_{1}&&0\\ &\ddots&\\ 0&&\lambda_{1}\end{smallmatrix}&\hbox to0.0pt{\kern 14.22636pt$N$\hss}\end{array}\right].$$

$$\underset{u,m}{\min }\Vert \mathbb{P}u-d{\Vert }_2^2+\lambda \Vert \mathbb{A}(m)u-b{\Vert }_2^2,$$

$$\underset{u,m}{\min }\underset{v}{\max }0+v^T[\mathbb{F}(m)u-s]+\frac{1}{2}\Vert \mathbb{F}(m)u-s{\Vert }_{\Lambda }^2,$$

$$\underset{u,m\in \mathcal{C}}{\min }\Vert m{\Vert }_{\text{TV}}$$

$$\text{subject to}\mathbb{F}(m)u=s,$$

$${\mathcal{L}}_A(m,u,v)=\Vert m{\Vert }_{\text{TV}}+v^T[\mathbb{F}(m)u-s]+\frac{1}{2}\Vert \mathbb{F}(m)u-s{\Vert }_{\Lambda }^2,$$

$${\mathcal{L}}_A(m,u,\bar{s})=\Vert m{\Vert }_{\text{TV}}+\frac{1}{2}\Vert \mathbb{F}(m)u-s-\bar{s}{\Vert }_{\Lambda }^2-\frac{1}{2}\Vert \bar{s}{\Vert }_{\Lambda }^2,$$

$$m^{k+1},u^{k+1}$$

$${\bar{s}}^{k+1}$$

$$u^{k+1}$$

$$m^{k+1}$$

$$s^{k+1}$$

$$u^{k+1}$$

$$s^{k+\frac{1}{2}}$$

$$m^{k+1}$$

$$s^{k+1}$$

$$\left[\begin{array}{c}{\lambda }_0^{\frac{1}{2}}\mathbb{P}\\ {\lambda }_1^{\frac{1}{2}}\mathbb{A}(m^k)\end{array}\right]u^{k+1}=\left[\begin{array}{c}{\lambda }_0^{\frac{1}{2}}[d+d^k]\\ {\lambda }_1^{\frac{1}{2}}[{b+b}^k]\end{array}\right],$$

$$\mathbb{u}^{k+1}=\big{[}\lambda_{0}\mathbb{P}^{T}\mathbb{P}+\lambda_{1}\mathbb{A}(\mathbb{m}^{k})^{T}\mathbb{A}(\mathbb{m}^{k})\big{]}^{-1}\big{[}\lambda_{0}\mathbb{P}^{T}[\mathbb{d}+\mathbb{d}^{k}]+\lambda_{1}\mathbb{A}(\mathbb{m}^{k})^{T}[\mathbb{b+b}^{k}]\big{]}.$$

$$\mathbb{A}(m)u^{k+1}$$

$$m^{k+1}=\underset{m\in \mathcal{C}}{\mathrm{argmin}}\Vert m{\Vert }_{\text{TV}}+\frac{{\lambda }_1}{2}\Vert \mathbb{L}(u^{k+1})m-y^k{\Vert }_2^2.$$

$$y^k=b+b^k-{\Delta u}^{k+1}.$$

$$J(m,p)=\sum \Vert (p_1{p}_2)\Vert +\frac{{\lambda }_1}{2}\Vert \mathbb{L}(u^{k+1})m-y^k{\Vert }_2^2,$$

$$\underset{m,p,p_0\in \mathcal{C}}{\mathrm{argmin}}J(m,p)$$

$$\text{subject to}p=\nabla m,$$

$$\nabla =\left[\begin{array}{c}\mathbb{I}\\ {\nabla }_1\\ {\nabla }_2\end{array}\right]\in {\mathbb{R}}^{3N\times N},p=\left[\begin{array}{c}{p}_0\\ p_1\\ p_2\end{array}\right]\in {\mathbb{R}}^{3N\times 1},$$

$${\mathcal{L}}_A(m,p,\bar{q})=J(m,p)+\frac{1}{2}\Vert \nabla m-p-{\bar{q}}^k{\Vert }_{\Gamma }^2-\frac{1}{2}\Vert {\bar{q}}^k{\Vert }_{\Gamma }^2,$$

$$\mathbb{\Gamma}=\left[\begin{array}[]{r@{}c|c|c@{}l}&\begin{smallmatrix}\gamma_{0}&&0\\ &\ddots&\\ 0&&\gamma_{0}\rule[-4.30554pt]{0.0pt}{8.61108pt}\end{smallmatrix}&\mbox{0}&\mbox{0}&\hbox to0.0pt{\kern 14.22636pt$N$\hss}\\ \hline\cr&\mbox{0}&\begin{smallmatrix}\rule{0.0pt}{8.61108pt}\gamma_{1}&&0\\ &\ddots&\\ 0&&\gamma_{1}\rule[-4.30554pt]{0.0pt}{8.61108pt}\end{smallmatrix}&\mbox{0}&\hbox to0.0pt{\kern 14.22636pt$N,$\hss}\\ \hline\cr&\mbox{0}&\mbox{0}&\begin{smallmatrix}\rule{0.0pt}{8.61108pt}\gamma_{2}&&0\\ &\ddots&\\ 0&&\gamma_{2}\end{smallmatrix}&\hbox to0.0pt{\kern 14.22636pt$N$\hss}\end{array}\right]$$

$$m^{k+1},p^{k+1}$$

Full text

K****eywords FWI, WRI, TV regularization, Bound constraints, ADMM, Split Bregman.

1 Introduction

During the last decade, full waveform inversion (FWI) has been used to estimate subsurface parameters (P and S wavespeeds, density, attenuation, anisotropic parameters) with a resolution close to the seismic wavelength by matching recorded and synthetic seismograms (Tarantola, 1984; Pratt et al., 1998; Virieux & Operto, 2009). From the numerical optimization viewpoint, the data-fitting/parameter-estimation problem underlying FWI is a nonlinear partial differential equation (PDE)-constrained optimization problem, where the equality constraint is the wave equation and the optimization parameters are embedded in the coefficients of the PDE. Due to the computational burden of multiple source modelling and the size of the data and parameter spaces, this PDE-constrained optimization problem is solved with iterative local (linearized) optimization techniques, namely gradient-based methods (Nocedal & Wright, 2006). Moreover, it is often solved with a reduced-space formulation, which means that the full search space that encompasses the unknown wavefield and the subsurface parameters is first projected onto the parameter space by computing exactly the incident wavefields in the current subsurface model before updating this later (Haber et al., 2000; Askan et al., 2007; Epanomeritakis et al., 2008). It is well acknowledged that the oscillating nature of seismic signals makes the reduced-space formulation highly nonlinear as the modelled seismograms computed in the current subsurface model may be too far away from the recorded ones to satisfy the cycle-skipping criterion, that is the modelled seismograms should predict the recorded traveltimes with an error lower than half a period (e.g. Virieux & Operto, 2009).

Beyond cycle skipping, other sources of nonlinearity and ill-posedness such as noise, uneven subsurface illumination, approximate wave physics and parameter cross-talks in multi-parameter reconstruction require the use of stabilizing or regularization techniques that drive the inversion towards subsurface models that satisfy some a priori assumptions. Among the penalization techniques, Tikhonov regularization is probably the most popular one and seeks to penalize the roughness of the subsurface model to force smooth reconstruction (see Benning & Burger (2018) for a review). As the subsurface may be better represented by piecewise smooth media with potentially sharp contrasts as in presence of salt, edge-preserving techniques such as total-variation (TV) regularization have been proposed to steer the inversion to the space of blocky structured models. TV regularization has been applied on several geophysical applications such as FWI (Askan et al., 2007; Anagaw & Sacchi, 2011; Guitton, 2012; Maharramov & Biondi, 2015; Peters & Herrmann, 2017; Brandsberg-Dahl et al., 2017; Esser et al., 2018; M. Kalita & Alkhalifah, 2018), seismic tomography (Gholami & Siahkoohi, 2010; Loris & Verhoeven, 2012), impedance inversion (Gholami, 2015, 2016), amplitude versus offset (AVO) inversion (Gholami et al., 2018), and seismic deconvolution (Gholami & Sacchi, 2013). For FWI applications, TV regularization can be implemented as a penalty function (Askan et al., 2007; Anagaw & Sacchi, 2011; Brandsberg-Dahl et al., 2017; Kazei et al., 2017) or as a constraint (Peters & Herrmann, 2017; Esser et al., 2018). Choosing the most suitable implementation strategy may depend on the prior information on the TV norm of the model and on the optimization method that is used to minimize the objective function (Alkhalifah et al., 2018). If the information about the value of TV norm of the model is available, TV regularization can be implemented as a constraint in the FWI objective function. Otherwise, one may resort to a penalty method with the difficulty to design an adaptive penalty parameter, which optimally balances over iterations the relative weight of the data misfit and the total variation of the model in the objective.

In this context, the objective of this study is to present a novel implementation of TV regularization and bound-constraints in frequency-domain FWI based upon wavefield reconstruction inversion (WRI). WRI has been originally proposed by van Leeuwen & Herrmann (2013) to extend the search space and mitigate the risk of cycle skipping accordingly. WRI recasts the PDE-constrained optimization problem underlying FWI into an unconstrained quadratic penalty method, where the penalty term is the $\ell_{2}$ norm of the source residuals (namely, the PDE-constraint violation) that is weighted by a positive penalty parameter $\lambda$. The penalty method relaxes the wave-equation constraint at the benefit of the data fitting during early iterations, hence mitigating the risk of cycle skipping. To make WRI computationally tractable, van Leeuwen & Herrmann (2013) perform the wavefield reconstruction and the subsurface parameter estimation in an alternating way: first, keeping the current subsurface model fixed, the wavefields, which best jointly fit the observations and satisfy the wave equation in a least-squares sense, are reconstructed for each source; Second, keeping the previously-reconstructed wavefields fixed, the subsurface parameters are estimated by least-squares minimization of the source residuals the wave-equation relaxation generated. This cycle being iterated until convergence. A nice property of the alternating-direction strategy is to linearise the parameter-estimation subproblem around the reconstructed wavefield because the wave equation constraint is bilinear. However, a significant pitfall of WRI resides in the tuning of the penalty parameter $\lambda$. Ideally, increasing values should be used during iterations to progressively enforce the wave-equation constraint and, hence satisfy the first-order optimality conditions of the original constrained problem with acceptable precision at the minimizer. A significant issue is that this continuation approach is tedious to implement and the Hessian is ill conditioned for large $\lambda$. Therefore, van Leeuwen & Herrmann (2013) implement WRI with a small preset value of $\lambda$, which leads to slow convergence and a subsurface model of limited accuracy.

Later, van Leeuwen & Herrmann (2016) reformulated WRI as a reduced penalty method implemented with a variable projection approach: the closed-form expression of the extended-domain reconstructed wavefield is injected as a function of the subsurface parameters in the penalty function instead of using this wavefield as a passive variable (i.e., independent to the subsurface parameters). Unlike the alternating-direction approach, this variable projection leaves the parameter-estimation subproblem non linear. van Leeuwen & Herrmann (2016) assess their method with a Gauss-Newton method (by opposition to the full Newton counterpart) to mitigate the computational burden. Moreover, using a sparse approximation of the Gauss-Newton Hessian makes the descent direction of the reduced approach identical to that of the alternating-direction WRI of van Leeuwen & Herrmann (2013). Also, Aravkin et al. (2017) analysed the convergence properties of the reduced penalty method when the full Hessian is taken into account and concluded that the variable projection penalty method is insensitive to the penalty parameter. However, this convergence property still needs to be verified against realistic numerical experiments.

To make the alternating-direction WRI of van Leeuwen & Herrmann (2013) (referred to as WRI in the following for sake of brevity) more independent to the penalty parameter, Aghamiry et al. (2019) have replaced the penalty method by an augmented Lagrangian method (Nocedal & Wright, 2006), leading to the so-called iteratively-refined (IR)-WRI method. As in WRI, IR-WRI performs the primal wavefield and parameter updates in an alternating mode, while the Lagrange multipliers (i.e., the dual variables) are updated with a gradient ascent method. As above mentioned, this alternating direction strategy makes the parameter estimation subproblem linear due to the bilinearity of the wave equation constraint. It follows from this linearization that the alternating direction strategy combined with the augmented Lagrangian method is equivalent to an extension of the alternating direction method of multiplier (ADMM) to biconvex problem (Boyd et al., 2010). Also, using a scaled form of the augmented Lagrangian, we recast IR-WRI as a penalty method where the right-hand sides (the data and the sources) in the objectives of the penalty function are iteratively updated with the running sum of the data and source residuals (the dual gradient steps). This reformulation of IR-WRI as a penalty method with right-hand side updating clearly draws some similarities and differences with WRI. The right-hand side updating makes IR-WRI largely insensitive to the penalty parameter for a wide range of preset values (Aghamiry et al., 2019, Their Figures 2 and 3). Using a moderate value of the penalty parameter allows for significant wave equation error and improved data fitting during early iterations for search space extension, without preventing the fulfilment of the wave equation constraint with small error at the minimizer (Nocedal & Wright, 2006, Chapter 17, Theorem 17.6). This adaptivity makes IR-WRI resilient to cycle skipping as WRI with however a much faster convergence toward a more accurate minimizer. The reader is referred to Aghamiry et al. (2019) for a thorough comparative convergence and accuracy analysis of WRI and IR-WRI based upon toy and complex large-contrast synthetic examples.

The objective of this study is to show how to interface TV regularization and bound constraints (hereafter, we refer to it as BTV regularization) with IR-WRI by taking advantage of the alternating-direction strategy implemented in ADMM and the split-Bregman variable splitting scheme developed by Goldstein & Osher (2009). More precisely, we recast the BTV regularized IR-WRI as a TV minimization problem subject to constraints, that are the modelled wavefield fit the observables and satisfy the wave equation with prescribed errors, and the model parameters preset bounds. As in WRI and IR-WRI, we solve the wavefield and subsurface parameter subproblems in an alternating mode. However, the later one involves now a combination of $\ell_{1}$ and $\ell_{2}$ norms related to the TV minimization and wave-equation error minimization, respectively, with additional bound constraints. This is managed by the split-Bregman variable splitting scheme, which de-couples the $\ell_{1}$ and $\ell_{2}$ components and bound constraints of the functional through the introduction of auxiliary variables and solves each related subproblem in sequence (Goldstein & Osher, 2009).

We first apply our method on a toy example corresponding to a high-velocity box-shape anomaly embedded in a background model where the velocity increases with depth. Then, we consider two more realistic examples corresponding to the left and central parts of the large-contrast 2004 BP salt model (Billette & Brandsberg-Dahl, 2004). We show that the BTV regularized IR-WRI converges to accurate minimizers when we start from a crude initial model and a realistic 3 Hz frequency. We also compare the results of WRI and IR-WRI without any priors, and when performed only with bound constraints and with BTV regularization to highlight the impact of each ingredient upon the quality of the results and the computational burden. We also assess the resilience of the method to noise by comparing the results that are obtained with noiseless and noisy data.

This paper is organized as follows. In the first section, we review the principles of WRI and IR-WRI. We first recast FWI as a feasibility problem and review different approaches that are suitable to solve PDE-constrained optimization problems such as penalty and augmented Lagrangian methods. Then, we review how we can easily interface some stabilizing terms with the feasibility problem through variable splitting and ADMM. In the second part, we present the results of the synthetic examples involving the inclusion model with two different starting models and the two targets of the BP2004 model with noiseless and noisy data. The results confirm that the combined use of TV regularization and bound constraint in the ADMM-based IR-WRI method defines a suitable framework to make high-resolution FWI immune to cycle skipping in large-contrast media.

2 Method

In the following, we first recast frequency-domain FWI as a bi-convex feasibility problem, which can be formulated as a constrained optimization problem with identically-zero objective function (Aghamiry et al., 2019, their appendix A). Then, we review penalty and augmented Lagrangian methods as optimization techniques to solve this constrained optimization problem with extended search-space, leading to WRI and IR-WRI, respectively. Finally, we interface bound constraints and isotropic TV regularization with IR-WRI by replacing the identically-zero objective function by the TV norm of the subsurface model, and we show how to solve efficiently the regularized IR-WRI with ADMM and split Bregman iterations.

2.1 WRI and IR-WRI principles

FWI can be formulated in the frequency domain as the following bi-convex feasibility problem (Aghamiry et al., 2019):

[TABLE]

with

[TABLE]

where $\mathbb{m}\in\mathbb{R}^{N\times 1}$ denotes the vector of discrete model parameters (here, the squared slowness), $\mathbb{u}\in\mathbb{C}^{N\times 1}$ the wavefield, $\mathbb{b}\in\mathbb{C}^{N\times 1}$ the source term, $\mathbb{d}\in\mathbb{C}^{M\times 1}$ the recorded wavefield (data) at receiver locations, and $\mathbb{P}\in\mathbb{R}^{M\times N}$ is a linear observation operator that samples the modelled wavefield at the receiver positions. The matrix $\mathbb{A(m)}\in\mathbb{C}^{N\times N}$ represents the discretized PDE Helmholtz operator (Pratt et al., 1998; Plessix, 2007; Chen et al., 2013).

[TABLE]

where $\omega$ is the angular frequency and $\mathbb{\Delta}$ is the discretized Laplace operator. The operator $\mathbb{C}$ encloses boundary conditions, which can be a function of $\mathbb{m}$ (e.g., Robin paraxial conditions, Engquist & Majda, 1977) or independent from $\mathbb{m}$ (e.g., sponge-like absorbing boundary conditions such as perfectly-matched layers, Bérenger, 1994). Also, the linear operator $\mathbb{B}$ can be used to spread the "mass" term $\omega^{2}\mathbb{C}(\mathbb{m})\text{diag}(\mathbb{m})$ over all the coefficients of the stencil to improve its accuracy following an anti-lumped mass strategy (Marfurt, 1984; Jo et al., 1996; Hustedt et al., 2004).

In the feasibility problem 1, we just want to find $\mathbb{m}$ and $\mathbb{u}$ that satisfy the constraint (we assume that the feasible set is non-empty, namely the constraint is consistent). The feasibility problem can be formulated as the following constrained optimization problem with identically-zero objective function (Boyd & Vandenberghe, 2004, Page 128)

[TABLE]

2.1.1 WRI penalty method

WRI implements this constrained optimization problem with a penalty method.

[TABLE]

where $\|\mathbb{x}\|_{\mathbb{Q}}^{2}:=\mathbb{x}^{T}\mathbb{Qx}$ for a vector $\mathbb{x}$ and a square matrix $\mathbb{Q}$ with superscript $T$ denoting matrix transposition. In equation 5, $\mathbb{\Lambda}$ is a diagonal matrix which includes the penalty parameters $\lambda_{0}$, $\lambda_{1}>0$ on its main diagonal,

[TABLE]

The objective function in equation 5 is a compact writing of the penalty function of van Leeuwen & Herrmann (2013)

[TABLE]

where $\lambda=\lambda_{1}/\lambda_{0}$. They solve this biconvex minimization problem with an alternating-direction approach to break down the full problem into a sequence of two linear sub-problems: A cycle of the algorithm first reconstructs, for each source, the wavefield $\mathbb{u}$ that best fits the data and satisfies the wave equation in a least-square senses for the current subsurface model. Then, the subsurface model is updated by minimization of the source residuals with a Gauss-Newton algorithm keeping the reconstructed wavefields fixed. A difficulty with the penalty method given by equation 7 resides in the tuning of the penalty parameter $\lambda$ during iterations. An increasing values of $\lambda$ should be used during iterations, known as the penalty algorithm, to progressively enforce the wave-equation constraint in iterations and hence satisfy the Karush-Kuhn-Tucker (KKT) optimality conditions (Nocedal & Wright, 2006) associated with the original constrained problem with acceptable precision. The main problem is that this continuation strategy is tedious to implement and a large $\lambda$ makes the problem severely ill-conditioned.

2.1.2 IR-WRI augmented Lagrangian method

To bypass this difficulty, IR-WRI implements the original constrained problem, equation 4, with the augmented Lagrangian (AL) method (Hestenes, 1969; Nocedal & Wright, 2006; Boyd et al., 2010; Bertsekas, 2016).

[TABLE]

where $\mathbb{\Lambda}$ is defined as in equation 6 and $\mathbb{v}\in\mathbb{C}^{(M+N)\times 1}$ is the Lagrangian multiplier (known as dual variable). Comparing the penalty function, equation 5, and the augmented Lagrangian function, equation 8, clearly shows that the augmented Lagrangian method combines a penalty method with a Lagrangian method. A first advantage of the augmented Lagrangian method relative to the penalty method is to prevent ill-conditioning by introducing explicit estimate of the Lagrange multiplier in the optimization (Nocedal & Wright, 2006, chapter 17). Moreover, the Lagrange multiplier gives the augmented Lagrangian method one more way of improving the accuracy of the minimizer in addition to the penalty parameter, hence allowing for a fixed value to be used for this latter (Nocedal & Wright, 2006, Theorem 17.6). Applying the alternating direction strategy of WRI on the augmented Lagrangian method leads to an adaptation of the alternating-direction method of multiplier (ADMM) to biconvex problem (Boyd et al., 2010, Section 9.2). One ADMM iteration first minimizes the augmented Lagrangian function with respect to the primal variables $\mathbb{u}$ and $\mathbb{m}$ via a single Gauss-Seidel like iteration (namely, fix one variable and solve for the other) and then update the Lagrangian multiplier via a gradient ascent method.

In the following section, we review each step of the ADMM-based IR-WRI algorithm when equipped with TV regularization and bound constraints. The reader is also referred to Aghamiry et al. (2019) for the detailed IR-WRI algorithm when no regularization is used.

2.2 BTV-regularized ADMM-based IR-WRI

To implement TV regularization and bound constraints in IR-WRI, we first recast FWI as a constrained TV minimization problem given by

[TABLE]

where $\mathbb{F}$ and $\mathbb{s}$ are defined as in equation 2, $\|\mathbb{m}\|_{\text{TV}}=\sum\sqrt{|{\mathbb{\nabla}}_{1}\mathbb{m}|^{2}+|{\mathbb{\nabla}}_{2}\mathbb{m}|^{2}}$ is the blockiness-promoting isotropic TV norm (Rudin et al., 1992), and $\mathbb{\nabla}_{1}$ and $\mathbb{\nabla}_{2}$ are first-order finite-difference operators in the horizontal and vertical directions, respectively. With notation abuse, the absolute sign, square power, and the square root operations are done component-wise, and the sum runs over all elements (the domain of parameters). Also $\mathcal{C}=\{\mathbb{x}\in\mathbb{R}^{N\times 1}~{}|~{}\mathbb{m}_{{l}}\leq\mathbb{x}\leq\mathbb{m}_{{u}}\}$ is the set of all feasible models bounded by the lower bound $\mathbb{m}_{{l}}$ and the upper bound $\mathbb{m}_{{u}}$.

Compared to the FWI definition given in the previous section, equation 4, we have replaced the identically-zero objective by the TV norm of the model and restricted the space of feasible models to $\mathcal{C}$. Accordingly, the augmented Lagrangian function for the problem defined by equation 9 is

[TABLE]

where the same notations as those of equation 8 are used.

Equation 10 can also be written in a more compact form as

[TABLE]

where $\mathbb{\bar{s}}=-\mathbb{\Lambda}^{-1}\mathbb{v}$ is the scaled dual variable and equation 11 is the scaled form of the augmented Lagrangian ((Boyd et al., 2010, Page 15) and Appendix A). The scaled form recasts the augmented Lagrangian method as a quadratic penalty method where the right-hand sides are updated with the scaled dual variables. This highlights similarities and differences between WRI and IR-WRI, since this right-hand side updating is lacking in the former.

The method of multipliers seeks to find the saddle point of the scaled augmented Lagrangian 11 through a primal descent - dual ascent updating resulting in the following iteration:

[TABLE]

for $k=0,1,...$ beginning with a prior estimate $\mathbb{\bar{s}}^{0}=\mathbb{0}$. The iteration 12 can be viewed as follows: we begin with a prior estimate of the dual $\mathbb{\bar{s}}^{0}$, and minimize the objective function with respect to the primal variables $\mathbb{m}$ and $\mathbb{u}$, equation 12a. Subsequently, we maximize the objective function with respect to the dual variable $\mathbb{\bar{s}}$ with a gradient ascent method when $\mathbb{m}$ and $\mathbb{u}$ are kept fixed, equation 12b. The steepest-ascent step, equation 12b, shows that the scaled dual variable $\mathbb{\bar{s}}$ is updated with the residual constraint violation of the current iteration. Remembering that the constraint gathers the observation equation $\mathbb{Pu}=\mathbb{d}$ and the wave equation $A(\mathbb{m})\mathbb{u}=\mathbb{b}$, equation 2, the scaled dual variable $\mathbb{\bar{s}}$ updates the right-hand sides of the quadratic penalty function, equation 12a, with the running sum of the data and source residuals in iterations. This right-hand side updating describes the well-known iterative solution refinement procedure for ill-posed linear inverse problems as reviewed by Aghamiry et al. (2019, Their appendix B). This process is iterated until convergence, i.e., when $\mathbb{F}(\mathbb{m}^{k+1})\mathbb{u}^{k+1}=\mathbb{s}$. In the following, we remove the bar of $\mathbb{\bar{s}}^{k}$ for the sake of simplicity.

Solving the subproblem 12a jointly for the primal variables ($\mathbb{m}$, $\mathbb{u}$) is computationally too intensive. A splitting method is useful here to break down the joint optimization over $\mathbb{m}$ and $\mathbb{u}$ into two subproblems (the readers can refer to Glowinski et al. (2017) for an overview of splitting methods). The alternating-direction method of multipliers (ADMM) (Boyd et al., 2010) provides a simple framework to achieve this goal via a single Gauss-Seidel like iteration leading to the following iteration:

[TABLE]

The ADMM iteration has decomposed the full problem into two subproblems associated with primal variables $\mathbb{u}$ and $\mathbb{m}$ by passing the primal update of one subproblem as a passive variable for the next subproblem. We show below that, taking advantage of the bi-convexity of the problem, this alternating-direction approach has linearized the primal subproblem for $\mathbb{m}$ around the reconstructed wavefield $\mathbb{u}$. The fact that the primal update of one subproblem is passed to the next subproblem implies obviously that the two subproblems are solved in sequence rather than in parallel as in ADMM for linear separable problems.

A second modification related to ADMM resides in the updating of the dual variable. In ADMM, the dual variables are updated only once per iteration after the primal-variable updates, equation 13c. A variant of ADMM, referred to as the Peaceman-Rachford splitting method (PRSM) (Peaceman & Rachford, 1955), consists of updating the Lagrange multipliers several times, once after the update of each primal variable (see He et al. (2014, Compare equations 1.3 and 1.4). One issue with PRSM relative to ADMM is that PRSM requires more restrictive assumptions to ensure its convergence, while it is always faster than ADMM whenever it is convergent (He et al., 2014) . This issue prompted He et al. (2014) to implement a relaxation factor (or, step length) $\alpha\in(0,1)$ to guarantee the strict contraction of the PRSM iterative sequence. Applying the strictly contractive PRSM algorithm to our minimization problem gives

[TABLE]

The reader is referred to Aghamiry et al. (2019, Their Figure 4) for a comparative numerical analysis of the convergence speed of ADMM and PRSM in IR-WRI. In our numerical tests we found that $\alpha=0.5$ can serve as a suitable value.

We now provide the closed-form solution of the two primal subproblems associated with $\mathbb{u}$ and $\mathbb{m}$.

2.2.1 Solving for $\mathbb{u}$

The primal subproblem associated with $\mathbb{u}$, equation 14a, is a linear optimization problem whose solution satisfies in a least-squares sense the following system of linear equations:

[TABLE]

where $\mathbb{d}^{k}$ and $\mathbb{b}^{k}$ are the components of the dual variable $\mathbb{s}^{k}$ associated with the observation-equation and wave-equation constraints, and are formed by the running sum of the data and source residuals in iteration (see equation 12). The closed-form expression of the reconstructed wavefield is given by

[TABLE]

The reconstructed wavefield can be computed numerically with linear algebra methods (direct or iterative methods) suitable for sparse matrices as reviewed by van Leeuwen & Herrmann (2016).

2.2.2 Solving for $\mathbb{m}$

In order to solve the bound-constrained TV regularized nonlinear problem described by equation 14c, we first tackle the nonlinearity issue by considering the special structure of the Helmholtz operator $\mathbb{A}$ given in equation 3. Using the following approximation

[TABLE]

we linearise the operator $\mathbb{A}$ with respect to $\mathbb{m}$ by building the matrix $\mathbb{C}$ from $\mathbb{m}^{k}$ to manage potential nonlinear boundary conditions. Note that this linearization step is not necessary when the PML absorbing conditions are used. We also exploit the bilinearity of the wave equation to permute $\mathbb{B}\mathbb{u}$ and $\mathbb{m}$ in the operator $\mathbb{L}(\mathbb{u}^{k+1})$. With these manipulations, the sub-problem associated with $\mathbb{m}$ is recast as

[TABLE]

where

[TABLE]

We solve the constrained optimization problem described by equation 18 with the split Bregman method (Goldstein & Osher, 2009) to decouple the TV minimization subproblem (first term in equation 18) from the $\ell_{2}$ subproblem (second term in equation 18) and force the box constraint. We introduce the auxiliary variables $\mathbb{p}_{0}$, $\mathbb{p}_{1}$ and $\mathbb{p}_{2}$ to perform this splitting, in which these variables being related to $\mathbb{m}$ by means of a simple equality constraint. The auxiliary variable $\mathbb{p}_{0}$ is used to enforce the box constraint, while $\mathbb{p}_{1}$ and $\mathbb{p}_{2}$ are the variables of the TV minimization problem. Let us define $(\mathbb{p}_{1}\;\mathbb{p}_{2})$ as a two-column matrix, $\|(\mathbb{p}_{1}\;\mathbb{p}_{2})\|=\sqrt{\mathbb{p}_{1}^{2}+\mathbb{p}_{2}^{2}}$ be a vector which contains $\ell_{2}$ norm of each row of $(\mathbb{p}_{1}\;\mathbb{p}_{2})$, and $\sum\|(\mathbb{p}_{1}\;\mathbb{p}_{2})\|$ be the mixed $\ell_{2,1}$ norm ($\ell_{1}$ norm of $\|(\mathbb{p}_{1}\;\mathbb{p}_{2})\|$), which promotes sparsity. Then by defining the objective function $J$ as

[TABLE]

the unconstrained optimization problem described by equation 18 can be written in a split and constrained form as

[TABLE]

where

[TABLE]

and $\mathbb{I}$ is the identity matrix.

The scaled augmented Lagrangian function (Appendix A) for the problem defined by equation 21 is

[TABLE]

where $\mathbb{\bar{q}}$ is the scaled Lagrangian multipliers and

[TABLE]

with the penalty parameters $\gamma_{0}$, $\gamma_{1}$, $\gamma_{2}>0$. Again, applying the method of multipliers to find the saddle point of the problem 22 gives

[TABLE]

where the bar of $\mathbb{\bar{q}}^{k}$ is removed for simplicity.

Substituting the explicit expression of $J$, equation 20, into equation 24a leads to the following PRSM iteration:

[TABLE]

Note that the weight $\frac{1}{2}$ in equations 25b and 25d has a similar role as $\alpha$ in equations 14b and 14d. Now we come up with a linear inverse subproblem for $\mathbb{m}$, equation 25a. Accordingly, the update $\mathbb{m}^{k+1}$ is obtained by solving the following system of linear equations in a least-squares sense:

[TABLE]

where we have substituted $\mathbb{y}^{k}$ by its explicit expression, equation 19. In equation 26, the first line describes the information carried out by the reconstructed wavefield to update $\mathbb{m}$ via the wave-equation rewriting, while the second line describes the action of the TV regularization and bound constraints on $\mathbb{m}$ via its linear relation with the auxiliary variable $\mathbb{p}$.

The closed-form expression of $\mathbb{m}$ is given by

[TABLE]

As for the linear system 16, $\mathbb{m}$ can be computed numerically with any suitable sparse linear algebra method.

The sub-problem for $\mathbb{p}$, equation 25c, is straightforward to solve. The objective function is separable with respect to the variable $\mathbb{p}_{0}$ and the variables $\mathbb{p}_{1}$ and $\mathbb{p}_{2}$ (i.e., the optimization can be performed for $\mathbb{p}_{0}$ and $\mathbb{p}_{1}$, $\mathbb{p}_{2}$ separately). The variable $\mathbb{p}_{0}$ is solution of the following linear inverse problem

[TABLE]

The $i$th element of the solution, $\mathbb{p}^{k+1}_{0}(i)$, is the closest element of $\mathbb{m}^{k+1}(i)-\mathbb{q}^{k}_{0}(i)$ to the desired set $[\mathbb{m}_{l}(i),\mathbb{m}_{u}(i)]$. Therefore,

[TABLE]

where the projection operator is $\text{proj}_{\mathcal{C}}(\bullet)=\min(\max(\bullet,\mathbb{m}_{l}),\mathbb{m}_{u})$.

The variables $\mathbb{p}_{1}$ and $\mathbb{p}_{2}$ are updated via the following proximity operator:

[TABLE]

where $\gamma=\gamma_{1}=\gamma_{2}$, and $\mathbb{z}_{i}=\nabla_{i}\mathbb{m}^{k+1}-\mathbb{q}_{i}^{k}~{}i=1,2$. Proximity operators are generalization of projection operators (Combettes & Pesquet, 2011). Equation 30 describes a separable optimization problem with respect to $\mathbb{p}_{1}$ and $\mathbb{p}_{2}$. Furthermore, $\mathbb{p}_{1}$ and $\mathbb{p}_{2}$ have closed-form expressions (Goldstein & Osher, 2009)

[TABLE]

where

[TABLE]

It can be seen that, for a single vector, the proximity operator defined in equation 32 reduces to soft thresholding.

Considering all the above-mentioned processes, a pseudocode for the BTV-regularized IR-WRI algorithm is summarized in Algorithm 1. Note that the lines 6-11 of the algorithm correspond to one ADMM iteration of the model-parameter updating. These operations could be iterated in an inner loop to update the BTV-regularized model several times after the wavefield reconstruction at each outer iteration. However, we observed numerically that only a single iteration of the inner loop guarantees the most efficient convergence of the full algorithm. This property has been noticed by Goldstein & Osher (2009) and is discussed more extensively in the framework of IR-WRI by Aghamiry et al. (2019). The inefficiency of the inner iterations can be understood by the fact that the original nonlinear problem is solved with an alternating-direction strategy (managed by the outer loop). This implies that each subproblem is solved from potentially inaccurate passive variables, this inaccuracy preventing an efficient minimization of the objective during inner iterations. Furthermore, in order to drive the algorithm, we assumed that the constraint is feasible. However, it has been shown that in the case of infeasible linear constraints the ADMM iteration can still produce approximate solutions that are stable (Frick et al., 2011; Jiao et al., 2016).

3 Numerical examples

3.1 Experimental setup and parameter tuning

We assess the performance of our BTV regularized IR-WRI against 2D mono-parameter synthetic examples. We start with a toy example built with a high-velocity inclusion model that is embedded in a background medium where velocity linearly increases with depth . To tackle more realistic applications, we proceed with two targets of the challenging 2004 BP salt model (Billette & Brandsberg-Dahl, 2004). With the BP salt case study, we seek to illustrate the potential of IR-WRI equipped with our BTV regularization to image salt bodies and sub-salt structures starting from crude initial models and realistic frequencies.

For all the numerical examples, forward modelling is performed with a 9-point stencil implemented with anti-lumped mass and PML absorbing boundary conditions (Chen et al., 2013). In this setting the diagonal matrix $\mathbb{C}$ contains the damping PML coefficients and does not depend on $\mathbb{m}$. With this setting, equation 17 does not require any approximation for linearization.

We will compare the results of WRI and IR-WRI to highlight the improved convergence history of IR-WRI resulting from the iterative updating of the right-hand sides in the penalty function associated with the scaled-form augmented Lagrangian, equations 11-14. We assume that our IR-WRI algorithm, when this right-hand side updating is not activated, is representative of the WRI penalty method (van Leeuwen & Herrmann, 2013). For a fair comparison, we will use the same experimental setup (penalty parameters and stopping criterion of iteration) for the two methods. We also compare the WRI and IR-WRI results when they are obtained without any priors ($\gamma_{0}=\gamma=0$), with only bound constraints ($\gamma=0$) and with BTV regularization ($\gamma_{0}\neq 0$ and $\gamma\neq 0$), where it is reminded that $\gamma_{0}$ and $\gamma$ are the penalty parameters that control the weight of the bound constraints and TV regularization, respectively, in the objective function (see equation 23).

We tune the different penalty parameters according to the following guideline. We start from the last subproblem of the splitting procedure and set the parameter $\gamma$, which controls the soft thresholding performed by the TV regularization, equation 31-32. In this study, we find that $\gamma=2\%\max{\|(\mathbb{z}_{1}\;\mathbb{z}_{2})\|}$ was a good pragmatical value. This tuning can be refined according to prior knowledge of the geological structure, coming from well log for example. In this study, we use the same weight for the bound constraints and the TV regularization: $\gamma_{0}=\gamma$. Once we set $\gamma$, we define $\lambda_{1}$ such that $\gamma/\lambda_{1}$ is a percentage of mean absolute value of the diagonal coefficients of $\mathbb{L}^{T}\mathbb{L}$ during the parameter estimation subproblem, equation 27. This percentage is set according to the weight that we want to assign to the TV regularization and bound constraints relatively to the wave equation constraint during the parameter estimation. Parameter $\lambda_{1}$ may be increased during iterations to reduce the weight of TV regularization and bound constraints near the convergence point. We found this adaptation useful when we start from very crude initial models. Finally, we set $\lambda_{0}$ such that $\lambda=\lambda_{1}/\lambda_{0}$ is a small fraction of the highest eigenvalue $\xi$ of the normal operator $\mathbb{A(m})^{-T}\mathbb{P}^{T}\mathbb{P}\mathbb{A(m})^{-1}$ during the wavefield reconstruction subproblem, equation 16, according to the criterion proposed by van Leeuwen & Herrmann (2016). In all the numerical tests, we use $\lambda=1\text{e-5}\xi$ and $\lambda=1\text{e-3}\xi$ for noiseless and noisy data, respectively. This tuning of $\lambda$ is indeed important because it controls the extension of the search space. A too high value of $\lambda$ reduces the weight of $\|\mathbb{Pu-d}\|_{2}^{2}$ during the wavefield reconstruction and makes IR-WRI to behave like a reduced approach. Conversely, using a small value for $\lambda$ fosters data fitting and expends the search space accordingly. However, a too small value can lead to a prohibitively high number of iterations of the augmented Lagrangian method before the wave equation constraint is fulfilled with a sufficient accuracy. Moreover, when data are contaminated by noise, a too small value for $\lambda$ will make the wavefield reconstruction to over-fit the data and drive WRI toward poor minimizer. We always use $\lambda$ as a fixed percentage of $\xi$ in iterations for both WRI and IR-WRI. This does not prevent IR-WRI to converge towards accurate minimizers thanks to the iterative error correction performed by the Lagrange multiplier updating. The reader is referred to Aghamiry et al. (2019) for a more thorough sensitivity analysis of IR-WRI to the penalty parameter $\lambda$.

3.2 Inclusion model

The subsurface model contains a sharp box-shape anomaly of sides 0.2 $\times$ 0.3 km with a velocity ($V_{P}$) of 5 km/s. It is embedded in a smooth background model where $V_{P}$ increases linearly with depth from 1.5 to 3.5 km/s (Figure 1). The model is 1.5 km long and 1 km deep, and is discretized with a 10 m grid interval. The regular surface acquisition consists of five sources (as depicted with yellow stars in Figure 1) and 65 receivers deployed on the surface. The source signature is a Ricker wavelet with a 5 Hz dominant frequency. We start the inversion from the true background model and invert simultaneously three frequency components (2.5, 5 and 7 Hz) with noiseless data. This frequency bandwidth has been selected to cover a significant band of vertical wavenumbers in the waveform-inversion sensitivity kernels, considering the limited aperture illumination provided by the surface acquisition. Moreover, a realistic starting frequency of 2.5Hz allows us to assess the resilience to cycle skipping of IR-WRI. A maximum number of iterations set to 70 is used as stopping criterion for all of the tests shown in Figure 2. When bound constraints are used, the bounds $\mathbb{m}_{l}$ and $\mathbb{m}_{u}$, equation 29, is set to the true minimum and maximum square slownesses, respectively.

We first compare WRI and IR-WRI results when bound constraints and TV regularization are not activated, i.e. $\gamma_{0}=\gamma=0$ (Figure 2(a-b)). WRI and IR-WRI reconstruct only the top of the anomaly with strongly overestimated velocities. Then, we add bound constraints in WRI and IR-WRI, using $\gamma_{0}/\lambda_{1}=0.01\zeta$, where $\zeta$ is the mean absolute value of the diagonal coefficients of $\mathbb{L}^{T}\mathbb{L}$. The bound-constrained WRI and IR-WRI only reconstruct the top of the anomaly as in Figure 2(a-b). However, the inclusion velocities are now well controlled by the bound constraints (Figure 2(c-d)). These first two tests show that IR-WRI reconstructs better the shape of the anomaly than WRI with however more significant artifacts on both sides of the anomaly. These artifacts may result from the deficit of horizontal-wavenumber illumination provided by the sparse limited-offset surface acquisition and by multi-scattering pollutions.

Then, we apply BTV regularization with $\gamma/\lambda_{1}=\gamma_{0}/\lambda_{1}=0.01\zeta$ (Figure 2(e-f)). Since the initial model matches the true velocity-gradient background model, we use a small value of $\gamma/\lambda_{1}$ (i.e., a high value of $\lambda_{1}$) and keep it constant in iterations to preserve the smooth components of the subsurface model. We show that the BTV regularized WRI still fails to reconstruct the full anomaly (Figure 2e). In contrast, BTV regularized IR-WRI keeps on improving the reconstruction of the anomaly in depth, while efficiently mitigating the oscillating artifacts (Figure 2f).

A vertical profile across the reconstructed anomaly also highlights some limitations of the BTV regularization (Figure 2f): below the anomaly, the BTV regularization superimposes staircase artifacts on the velocity gradient, consistently with the piecewise constant assumption underlying TV regularization.

To emphasize the resilience to cycle skipping of the BTV regularized WRI and IR-WRI, we repeat this toy example using a 2.2 km/s homogeneous velocity model as initial model (Figure 3). We perform a first test without any prior. Compared to the previous test, we just stabilize the inversion by adding a small damping term (=$0.01\zeta$) to $\mathbb{L}^{T}\mathbb{L}$. Compared to Figure 2(a-b), the artifacts have a much stronger imprint due to the inaccuracy of the starting model (Figure 3(a-b)). Then, we move to bound-constrained and TV regularized tests (Figure 3(c-f)). To decrease the above-mentioned artifacts, we assign a high initial weight to the BTV regularization ($\gamma/\lambda_{1}=\gamma_{0}/\lambda_{1}=\zeta$) and decrease it by a factor 2 every 10 iterations until it reach a minimal value set to $0.01\zeta$ (i.e., the constant value previously used ). In accordance with the former test, bound constraints alone are not sufficient to reconstruct the bottom part of the anomaly and cancel out the oscillating artifacts (Figure 3(c-d)). In contrast, BTV IR-WRI achieves these two goals, although it leaves a significant staircase footprint below the anomaly (Figure 3f). Compared to Figure  2f, the edges of the anomaly are better reconstructed at the expense of the background velocity-gradient model. This results because more aggressive TV regularization was used during the early iterations of this test, allowing for a better reconstruction of the blocky components of the medium, while injecting undesired staircase footprint on its smooth components. As for the former test, IR-WRI clearly outperforms WRI due to the more efficient solution refinement procedure resulting from the right-hand side updating.

We also show the joint evolution in iterations of the observation-equation and wave-equation errors (Figure 4) and the wavefield and subsurface model errors (Figure 5), when the initial model is the true velocity-gradient background model and the homogeneous model. IR-WRI fits the data and wave equation better than WRI after 70 iterations with both initial models because the right-hand side updating embedded in IR-WRI cancels out more efficiently the data and source residuals in iterations and, hence better refines the solution accordingly. Moreover, BTV regularization in IR-WRI further improves the data and wave equation fit for both initial models because it reduces more efficiently the oscillating artifacts in the reconstructed velocity model and better reconstructs the edges of the anomaly (Figure 2f). As above mentioned, the oscillating artifacts may result from the deficit of horizontal wavenumber illumination generated by the limited-offset surface acquisition. In this framework, the prior contained in the BTV regularization efficiently narrows the null space of the inversion. The more complex zigzag path followed by IR-WRI relative to WRI in the ($\|\mathbb{Pu}^{k}-\mathbb{d}\|_{2}-\|\mathbb{A(m}^{k})\mathbb{u}^{k}-\mathbb{b}\|_{2}$) plane highlights how the joint updating of the data and source by their associated residuals dynamically balances the weight of the two objective functions in iterations. This more complex convergence history of IR-WRI, which has been already noticed in Aghamiry et al. (2019), suggests that the right-hand side updating perform a self-adaptive weighting of the two competing objective functions driven by the relative reduction of the data and source residuals in iterations. This zigzag convergence trend translates also into more complex path in the $(\|\mathbb{u}^{k}-\mathbb{u}^{*}\|_{2}/\|\mathbb{u}^{*}\|_{2}-\|\mathbb{m}^{k}-\mathbb{m}^{*}\|_{2}/\|\mathbb{m}^{*}\|_{2})$ plane ($\mathbb{m}^{*}$ and $\mathbb{u}^{*}$ denote the true model and wavefield, respectively), which illustrates how the solution refinement is pushed toward the wavefield reconstruction or the velocity model estimation according to the self-adaptive weighting of the observation-equation and wave-equation objective functions (Figure 5). Note that the relative model errors increase in the case of the initial velocity-gradient model (Figure 5a,b). This results because the smooth background model is degraded by the oscillating artifacts and the staircase footprint during the sharp inclusion reconstruction. Indeed, this degradation of the smooth components has a much higher weight in the $\ell_{2}$ misfit function than the more accurate reconstruction of the blocky components. However, this increase is much more moderate in IR-WRI (Figure 5b) than in WRI (Figure 5a). Finally, we plot the TV norm of the reconstructed models in iterations for the different tests (Figure 6). As expected, the models reconstructed with bound constraints and TV regularization match better the TV of the true model. The BTV WRI model matches slightly better the TV of the true model than the BTV IR-WRI one, in particular when the initial model is the true background model. Indeed, this does not reflect that WRI better reconstructs the anomaly than IR-WRI. Instead, it reflects the slower convergence of WRI relatively to IR-WRI which contributes to keep the background model smooth (namely, which a TV close to 0 and equal to that of the true model).

3.3 2004 BP salt model - central target

We now consider a more realistic application with a first target of the challenging 2004 BP salt model. 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 (Billette & Brandsberg-Dahl, 2004). The first selected target corresponds to the central part of the 2004 BP salt model characterized by a deeply rooted salt body (Figure 7a).

The subsurface model is 8.8 km wide and 2.9 km deep, and is discretized with a 25 m grid interval. We used 50 sources spaced 175 m apart on the top side of the model. The source signature is a Ricker wavelet with a 10 Hz dominant frequency. A line of receivers with a 50 m spacing are deployed at the surface leading to a stationary-receiver acquisition.

We start the inversion from a smoothed version of true velocity model where the imprint of the salt body and other structures were cancelled out (Figure 7b) and invert the 3-Hz frequency with noiseless data.

We compare the results of WRI and IR-WRI with bound constraints and BTV regularization. To highlight the specific role of bound constraints, we activate them after 21 iterations. Since we start from a rough initial model, we set $\gamma_{0}/\lambda_{1}=\gamma/\lambda_{1}=\zeta$ and decrease them during iterations in a manner similar to the box-shape anomaly test with homogeneous starting model. Also, we add a damping term to $\mathbb{L}^{T}\mathbb{L}$ with a weight equal to $0.01\zeta$ to further stabilize the inversion. We stop inversion after 70 iterations. The estimated models are shown in Figure 7 together with horizontal profiles at 1.65 km depth and vertical profiles at 4.15 km distance extracted from the true, initial and reconstructed models. As for the inclusion test, WRI fails to reconstruct the salt body and the subsalt structure because the data and source residuals are not re-injected in the right-hand sides of the penalty function at each iteration as in equation 13, leading to a stagnant convergence of the inversion (Figure 7c,e). When IR-WRI is applied with bound constraints alone, the reconstructed model is affected by noise with a periodic horizontal pattern. This noise likely results from the monochromatic nature of the inversion, multi-scattering within the salt body and limited illumination of the horizontal wavenumbers of the salt body leading to wraparound (Figure 7d). The BTV regularized IR-WRI mitigates efficiently this noise without degrading the resolution of the salt body and the sub-salt structures (Figure 7f).

Figure 8(a) shows the joint evolution in iterations of the data misfit and the wave-equation error. As for the inclusion test, note the zigzag path followed by the IR-WRI objective functions over iterations. Also, the joint evolution of wavefield and subsurface model errors and the evolution of TV norm over iteration are shown in Figure 8(b,c). The TV norm evolution emphasizes how the bound constraints fasten the convergence of TV-regularized IR-WRI after iteration 20.

We continue the inversion at higher frequencies using the final models of the 3 Hz inversion as initial model (Figure 7c-f). We used small batches of two frequencies with one frequency overlapping between two consecutive batches, moving from the low frequencies to the higher ones according to a classical frequency continuation strategy. We set $\gamma_{0}/\lambda_{1}=\gamma/\lambda_{1}=0.01\zeta$ and remove the damping of $\mathbb{L}^{T}\mathbb{L}$. The starting and final frequencies are 3.5 Hz and 12 Hz and the sampling interval in one batch is 0.5 Hz. The algorithm performs at most 15 iterations per frequency batch and the number of iterations that have been performed is 170. The inversion results are shown in Figure 9. WRI with bound constraints and BTV (Figure 9(a,c)) fails to converge toward satisfactory results, while BTV-regularized IR-WRI converges to accurate velocity model, although a significant imprint of the TV regularization is shown (Figure 9d). When IR-WRI is performed with only bound constraints, the oscillating artifacts are not cancelled out (Figure 9b). This highlights the role of TV regularization in reconstructing blocky structures and removing wraparound artifacts.

3.4 2004 BP salt model - Left target

We now consider a second target of the 2004 BP salt model located on the left side of the model (Figure 10a). This target was previously used by among others Métivier et al. (2016), Brandsberg-Dahl et al. (2017), Esser et al. (2018) and M. Kalita & Alkhalifah (2018) for FWI applications. The subsurface model is 16250 m wide and 5825 m deep, and is discretized with a 25 m grid interval. We used 108 sources spaced 150 m apart on the top side of the model. The source signature is a Ricker wavelet with a 10 Hz dominant frequency. A line of receivers with a 25 m spacing are deployed at the surface leading to a stationary-receiver acquisition. We perform IR-WRI with bound constraints alone and with BTV regularization, for noiseless and noisy data.

We used a crude laterally-homogeneous velocity-gradient model as initial model (Figure 10b). 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. The starting and final frequencies are 3 Hz and 13 Hz and the sampling interval in one batch is 0.5 Hz. The stopping criterion of iteration for each batch is given by

[TABLE]

where $k_{max}$ denotes the maximum iteration count, $\varepsilon_{b}$=1e-3, and $\varepsilon_{d}$=1e-5 for noiseless data and $\varepsilon_{b}$=1e-3 , $\varepsilon_{d}$= noise level of batch for noisy data. 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 frequency of the second and third path is 6 Hz and 8.5 Hz, respectively. The IR-WRI models inferred from noiseless data with bound constraints and with BTV regularization are shown in Figure 10(c-d). The number of iterations that have been performed with bound constraints and with BTV regularization are 441 and 340, respectively. Direct comparison between the true model, the starting model and the IR-WRI models along three vertical logs cross-cutting the salt body at 5 km, 7.5 km and 10 km distance (vertical dashed lines in Figure 10a) are shown in Figure 11a. The results show the resilience of IR-WRI to cycle skipping with a pretty accurate reconstruction of the salt body and sub-salt structures (Figure 10(c-d)). However, the model obtained with bound constraints alone shows high-frequency noise in the salt body and below (Figure 10c). This noise can result from Gibbs phenomenon resulting from the frequency decimation and artifacts resulting from multi scattering. The BTV regularization efficiently removes these artifacts except those resulting from truncation of the acquisition near the left end of the model at 5-km depth (Figure 10d). The inversion captures reasonably well the subsalt structures, including the low-velocity over-pressure zone at (x,z)=(7.5km,4km) as well as the smooth velocity variations.

When noisy data are used (Figure 10(e-f) and 11b), the number of iterations that have been performed with bound constraints alone and with BTV regularization is 263 and 254, respectively. As for the noiseless case, a direct comparison between the true model, the starting model and the two IR-WRI models along three vertical logs cross-cutting the salt body at 5 km, 7.5 km and 10 km distance is shown in Figure 11b. The artifacts have now a more significant imprint when only bound constraints are used (Figure 10e). Accordingly, we apply a more aggressive TV regularization to obtain the results shown in Figure 10f. The artifacts have been efficiently removed with however a more obvious imprint of the piecewise-constant approximation underlying BTV regularization. This blocky pattern is clearly visible in deep part of the vertical profiles of Figure 11b where velocity gradients have been replaced by stack of constant-velocity layers.

4 Discussion

We have implemented bound constraints and TV regularization in the wavefield reconstruction inversion (WRI) method (van Leeuwen & Herrmann, 2013) which has been recently improved by Aghamiry et al. (2019) in the framework of the alternating-direction method of multiplier (ADMM), leading to the iteratively-refined WRI (IR-WRI). To do this, we formulate IR-WRI as a TV minimization problem subject to observation-equation, wave-equation and bound constraints. We use a scaled form of the method of multiplier to recast the augmented Lagrangian function as a quadratic penalty function in which the scaled Lagrange multipliers, here the running sum of the data and source residuals, are added to the right-hand sides (data and source) at each iteration. At each iteration of the workflow, we perform a first ADMM step to break down the wavefield reconstruction and the parameter estimation into a sequence of two subproblems. The wavefield reconstruction subproblem is a linear inverse problem which can be solved with direct or iterative methods for sparse linear systems (Duff et al., 1986; Saad, 2003). Once one iteration of the wavefield reconstruction has been performed, we tackle the parameter estimation subproblem which involves the mixed $\ell_{1,2}$ TV norm of the model, the $\ell_{2}$ wave-equation objective and the bound constraints. We apply a second ADMM step to decompose this multi-objective optimization problem into a sequence of simpler subproblems. A first interesting property is the bilinearity of the wave equation constraint with respect to the wavefield and the parameter, which makes the $\ell_{2}$ wave-equation objective quadratic. A second key point is to introduce auxiliary variables, which are linearly related to the subsurface parameters and describes the TV of the model and the bound constraints. They allows for the de-coupling between the $\ell_{1}$ and $\ell_{2}$ components of the penalty function following the split Bregman method proposed by Goldstein & Osher (2009). After this de-coupling, the subsurface parameter are first updated by minimizing the source residuals through the resolution of a sparse BTV-regularized Gauss-Newton system, before updating the auxiliary variables with proximal operators. This cycle is iterated until convergence.

There are some key differences between our implementation of BTV-regularized WRI and previous ones based upon projected-gradient method (Peters & Herrmann, 2017) and primal dual hybrid gradient (PDHG) methods (Esser et al., 2018; Yong et al., 2018) that we review below. In the projected-gradient method, Peters & Herrmann (2017) first update model parameters with FWI or WRI and then project the updated model into the intersection of the TV and bound constraint with Dykstra projection algorithm (Boyle & Dykstra, 1986), an alternating projection onto some constraints until satisfying all of them. Therefore, their workflow is subdivided in two different parts, that are the model update followed by the projection onto intersection of all constraints (Peters & Herrmann, 2017, their equation 7). Unlike Peters & Herrmann (2017) method which rely on independent update and projection steps, all the ingredients of BTV-regularized IR-WRI (i.e., wavefield reconstruction, parameter estimation, TV regularization and bound constraints) are consistently integrated in the theoretical framework of ADMM optimization (the readers can also refer to Maharramov & Levin, 2015).

Esser et al. (2018); Yong et al. (2018) implement TV regularization and bound constraints in the reduced variable projection WRI of van Leeuwen & Herrmann (2016) with PDHG. PDHG is a method to solve constrained optimization problems by alternating gradient descent (for primal variable) and gradient ascent (for dual variable), which can be interpreted as linearized ADMM (Goldstein et al., 2015). PDHG can be helpful if the least squares minimizations embedded in ADMM are difficult to solve efficiently (Goldstein et al., 2015). This is not really the case in WRI, which mainly requires to solve two sparse linear systems for $\mathbb{u}$ and $\mathbb{m}$. Meanwhile, the selection of step size in PDHG that guarantees fast convergence, or even convergence at all, is not intuitive at all and can make PDHG impractical. This issue prompted Goldstein et al. (2015) to develop step size tuning rules which contribute to make PDHG self-adaptive as illustrated recently by Yong et al. (2018) in the frame of WRI. Beyond PDHG method, Esser et al. (2018); Yong et al. (2018) minimize data and source residuals with a penalty method, which lacks the convergence property of the augmented Lagrangian method promoted in IR-WRI. Moreover, they implement TV regularization as a hard constraint in the parameter-estimation subproblem through a Lagrangian function, whose saddle point is estimated with PDHG. In contrast, we implement the TV regularization as a soft constraint after introducing the auxiliary variable $\mathbb{p}$ and solve the regularized sub-problem for $\mathbb{m}$ with the split Bregman method (or equivalently ADMM), equation 22. This gives us the necessary flexibility to implement aggressive regularization during early iterations and relax it progressively, when very crude initial models are used. With more accurate initial models, the augmented Lagrangian embedded in the split Bregman method, equation 22, allows for constant penalty parameter to be used. These penalty parameters are used as step lengths in the augmented Lagrangian method, hence leading to a self-adaptive TV regularization implementation rid of tedious TV-norm ball continuation strategies and/or adaptive step lengths (Esser et al., 2018; Yong et al., 2018). Third, while Esser et al. (2018, equation 18) and Yong et al. (2018, equation 41) implement bound constraints as hard constraints by projection at each iteration of the TV-regularized perturbation model onto the feasible set defined by the bound constraints, we implement bound constraints consistently with TV regularization in the framework of the method of multiplier, equations 26-30. Finally, our approach does not require prior guess of TV-norm ball because we implement TV regularization as a minimization problem rather than as a constraint. That being said, Esser et al. (2018) show promising results with the $\tau$ continuation approach and asymmetric TV norm on the BP salt model starting from a very crude initial model. It will be interesting in future work to assess the benefit of asymmetric TV norm in BTV IR-WRI.

Indeed, implementation of BTV-regularized IR-WRI requires to setup different parameters, which have been discussed at the beginning of the section Experimental setup and parameter tuning. The penalty parameter which controls the relative weight between the wave equation objective and the data fitting objective during the wavefield-reconstruction subproblem has been discussed in length in Aghamiry et al. (2019). Aghamiry et al. (2019) have concluded that a fixed penalty parameter can be used during iterations because the accuracy of the minimizer in the method of multiplier is controlled both by the penalty parameter and the accuracy of the updated multiplier (Nocedal & Wright, 2006, Theorem 17.6). In this study, we suggest good pragmatical values of this penalty parameter as a percentage of the maximum eigenvalue of the augmented wave-equation normal operator for noiseless and noisy data. Optimal values may be refined by trial and error to prevent noisy data over-fitting and keep the iteration count within reasonable limits. The TV and bound parameters should be easily determined from well logs or a priori geological knowledge. Finally, the relative weight between the TV regularization and the wave-equation constraint, which is controlled by $\lambda_{1}$, needs also to be estimated during the parameter-estimation subproblem. As above mentioned, this penalty parameter is kept fixed during iterations or is progressively decreased to relax the TV regularization and bound constraints near the convergence point, depending of the accuracy of the initial model. The reader is also referred to Goldstein & Osher (2009, section 2.2) who discuss the sensitivity of the split Bregman method (an optimization method similar to ADMM (Esser, 2009)) to the penalty parameter for $\ell_{1}$-regularized problems.

One drawback of the BTV regularization is related to the piecewise constant approximation underlying TV regularization, which tends to superimpose some blocky patterns on the smooth part of the subsurface. To overcome this issue, ongoing work seeks to optimally combine Tikhonov and TV regularization in IR-WRI (Gholami & Hosseini, 2013; Aghamiry et al., 2018).

Other perspective works involve extension to multiparameter reconstruction, 3D geometries and application to real data to further assess the potential and limits of IR-WRI. For 3D applications, the wavefield reconstruction in the frequency domain requires to solve a large-scale linear algebra problem, equation 16. Operto et al. (2015), Amestoy et al. (2016), and Operto & Miniussi (2018) have shown the computational efficiency of 3D frequency-domain FWI based on sparse direct solvers in the 3.5-10 Hz frequency band for dense stationary-recording ocean-bottom cable acquisitions. Mary (2017) showed that block low-rank multifrontal solver allows one to tackle numerical problems involving up to 100 million unknowns. The symmetry of the normal operator, equation 16, should balance the computational overhead resulting from the higher number of non-zero coefficients relative to the impedance matrix ${\bf{A}}$. Note that ${\bf{A}}$ was processed as an unsymmetric matrix in the above-mentioned references due to discretization issues, although more suitable discretizations which preserve the symmetry of the impedance matrix may be considered in the future (Pratt & Smithyman, 2018). Alternatively, domain decomposition methods suitable for Helmholtz problems can be interfaced with hybrid direct-iterative solvers to perform wavefield reconstruction in bigger computational domains (Dolean et al., 2015). These approaches may be suitable for stationary-receiver acquisition involving a more limited number of reciprocal sources sich as ocean bottom seismometer acquisitions.

Let’s add that BTV IR-WRI should be easily implemented in the time-domain formulation of WRI recast as an extended-source waveform inversion (Wang et al., 2016; Huang et al., 2018).

5 Conclusion

We have presented a new method to implement TV regularization and bound constraints in frequency-domain FWI based on wavefield reconstruction (WRI). In a previous study, we have reformulated WRI in the framework of the alternating-direction method of multiplier (ADMM), leading to the iteratively-refined wavefield-reconstruction inversion method (IR-WRI). We have shown how the augmented Lagrangian embedded in ADMM makes IR-WRI weakly sensitive to a wide range of penalty parameter thanks to the Lagrange multiplier updating. Using a small value of this penalty parameter efficiently extends the search space during early iterations to foster data fitting without preventing the wave-equation constraint to be fulfilled at the convergence point with a preset prescribed error. IR-WRI performs a first ADMM step to alternate wavefield reconstruction and subsurface parameter estimation as in the original WRI method. When BTV regularization is used, we perform a second ADMM step to decompose the BTV-regularized parameter estimation sub-problem into a sequence of two simpler subproblems through the introduction of auxiliary variables. This variable splitting allows for the de-coupling between the $\ell_{1}$ and the $\ell_{2}$ components of the penalty function according to the so-called split Bregman method. An interesting property exploited by WRI methods is the bilinearity of the wave equation constraint with respect to the wavefield and the parameter, which makes the $\ell_{2}$ wave-equation objective of the second sub-problem quadratic. Our implementation of BTV-regularization in IR-WRI with ADMM (or, equivalently split Bregman) provides a versatile framework to cascade constraints and regularization of different nature and is reasonably easy-to-tune due to the limited sensitivity of the augmented Lagrangian method to the choice of the penalty parameters. For challenging subsurface targets with contrasted structures such as salt bodies, we have shown that our BTV-regularized WRI shows a high resilience to cycle skipping and noise and efficiently mitigates high-frequency artifacts associated with incomplete illumination and multi scattering without detriment to the resolution of the imaging.

6 Acknowledgements

This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by AkerBP, CGG, CHEVRON, EXXON-MOBIL, JGI, PETROBRAS, SCHLUMBERGER, SHELL, SINOPEC, EQUINOR and TOTAL. This study was granted access to the HPC resources of SIGAMM infrastructure (http://crimson.oca.eu), hosted by Observatoire de la Cote d’Azur and which is supported by the Provence-Alpes Côte d’Azur region, and the HPC resources of CINES/IDRIS under the allocation 046091 made by GENCI.

Appendix A Scaled form of augmented Lagrangian

Let’s start with the following constrained problem

[TABLE]

The augmented Lagrangian function for the problem is (Nocedal & Wright, 2006)

[TABLE]

The problem 35 can be written in a more compact form by introducing the scaled dual variable $\bar{\mathbb{q}}=-\frac{\mathbb{v}}{\lambda}$ and adding and subtracting the term $\frac{\lambda}{2}\|\bar{\mathbb{q}}\|_{2}^{2}$ to the augmented Lagrangian 35. In this case, we arrive at the following scaled-form of the method of multipliers:

[TABLE]