Compound Regularization of Full-waveform Inversion for Imaging Piecewise Media

TL;DR

This paper introduces compound regularization techniques combining Tikhonov and total variation functionals within wavefield reconstruction inversion to improve seismic imaging of piecewise media, demonstrating superior reconstruction accuracy.

Contribution

It develops and implements hybrid regularizers using infimal convolution and convex combination in WRI, enhancing model reconstruction for piecewise media in seismic imaging.

Findings

IC-based regularization yields lowest reconstruction error.

Compound regularization outperforms individual Tikhonov or TV methods.

WRI with compound regularization improves imaging accuracy.

Abstract

The nonlinear and ill-posed nature of full waveform inversion (FWI) requires us to use sophisticated regularization techniques to solve it. In most applications, the model parameters may be described by physical properties (e.g., wave speeds, density, attenuation, anisotropic parameters) which are piecewise functions of space. Compound regularizations are thus necessary to reconstruct properly such parameters 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 prevent the so-called 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 of parameters is 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 types of compound regularizations in the WRI optimization problem using the alternating direction method of multipliers (ADMM). Then, we assess our regularized WRI in the framework of seismic imaging applications. Using a wide range of subsurface models, we conclude that 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.

Tables

Table 1. TABLE I : Different regularization function.

Nick name	Expression of $\Phi \left(𝕞\right)$
DMP	${\Vert 𝕞\Vert }_2^2$
Tikhonov	${\Vert {\nabla }^2𝕞\Vert }_2^2$
TV	${\Vert \nabla 𝕞\Vert }_1$
JTT	$\left(1-\alpha \right){\Vert {\nabla }^2𝕞\Vert }_2^2+\alpha {\Vert \nabla 𝕞\Vert }_1$
TT	$\underset{𝕞={𝕞}_1+{𝕞}_2}{\min }\{\left(1-\alpha \right)\Vert {{\nabla }^2{𝕞}_2\parallel }_2^2+\alpha {\Vert \nabla {𝕞}_1\Vert }_1\}$
TGV	$\underset{𝕞={𝕞}_1+{𝕞}_2}{\min }\{\left(1-\alpha \right)\Vert {{\nabla }^2{𝕞}_2\parallel }_1+\alpha {\Vert \nabla {𝕞}_1\Vert }_1\}$

Table 2. TABLE 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

Equations

$$\begin{array}{rl}& \underset{u,m\in \mathcal{C}}{\min }\Phi (m)\\ & \text{subject to}\mathbb{A}(m)u=b,\mathbb{P}u=d,\end{array}$$

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

$$\underset{m\in \mathcal{C}}{\min }\Phi (m)+\frac{\lambda }{2}||\mathbb{P}{\mathbb{A}}^{-1}(m)b-d|{|}_2^2,$$

$$\underset{u,m\in \mathcal{C}}{\min }\Phi (m)+\frac{\lambda }{2}\Vert \mathbb{P}u-d{\Vert }_2^2+\frac{\gamma }{2}\Vert \mathbb{A}(m)u-b{\Vert }_2^2,$$

$$\begin{array}{rl}\underset{u,m\in \mathcal{C}}{\min }\underset{v_1,v_2}{\max }\Phi (m)& +\frac{\lambda }{2}\Vert \mathbb{P}u-d{\Vert }_2^2+\frac{\gamma }{2}\Vert \mathbb{A}(m)u-b{\Vert }_2^2\\ & +v_1^T[\mathbb{P}u-d]+v_2^T[\mathbb{A}(m)u-b],\end{array}$$

$$\begin{array}{rl}\underset{u,m\in \mathcal{C}}{\min }\underset{v_1,v_2}{\max }\Phi (m)& +\frac{\lambda }{2}\Vert \mathbb{P}u-d+\frac{1}{\lambda }{v}_1{\Vert }_2^2-\frac{\lambda }{2}\Vert v_1{\Vert }_2^2\\ & +\frac{\gamma }{2}\Vert \mathbb{A}(m)u-b+\frac{1}{\gamma }{v}_2{\Vert }_2^2-\frac{\gamma }{2}\Vert v_2{\Vert }_2^2,\end{array}$$

$$\begin{array}{rl}& \underset{u,m\in \mathcal{C}}{\min }\Phi (m)+\frac{\lambda }{2}\Vert \mathbb{P}u-d-d^k{\Vert }_2^2\\ & +\frac{\gamma }{2}\Vert \mathbb{A}(m)u-b-b^k{\Vert }_2^2,\\ & d^{k+1}=d^k+d-\mathbb{P}u,\\ & b^{k+1}=b^k+b-\mathbb{A}(m)u,\end{array}$$

$$u^{k+1}$$

$$m^{k+1}$$

$$d^{k+1}$$

$$b^{k+1}$$

$$\Vert x{\Vert }_2^2=\sum _{i=1}^n{x}_i^2,$$

$$\Vert x{\Vert }_1=\sum _{i=1}^n|x_i|,$$

$$\{\begin{array}{l}({\nabla }_{x}f{)}_{i,j}=f_{i,j}-f_{i,j-1},\\ ({\nabla }_{z}f{)}_{i,j}=f_{i,j}-f_{i-1,j},\end{array}$$

$$(|\nabla f|{)}_{i,j}=\sqrt{({\nabla }_{x}f{)}_{i,j}^2+({\nabla }_{z}f{)}_{i,j}^2}.$$

$$\{\begin{array}{l}({\nabla }_{xx}f{)}_{i,j}=f_{i,j-1}-2f_{i,j}+f_{i,j+1},\\ ({\nabla }_{zz}f{)}_{i,j}=f_{i-1,j}-2f_{i,j}+f_{i+1,j},\end{array}$$

$$(|{\nabla }^{2}f|{)}_{i,j}=\sqrt{({\nabla }_{xx}f{)}_{i,j}^2+({\nabla }_{zz}f{)}_{i,j}^2}.$$

$$({\nabla }_{xz}f{)}_{i,j}=f_{i,j}-f_{i,j-1}-f_{i-1,j}+f_{i-1,j-1}.$$

$$(|{\nabla }^{2}f|{)}_{i,j}=\sqrt{({\nabla }_{xx}f{)}_{i,j}^2+2({\nabla }_{xy}f{)}_{i,j}^2+({\nabla }_{yy}f{)}_{i,j}^2},$$

$${\Phi }_{\alpha }(x)={\alpha }_1{\Phi }_1(x)+...+{\alpha }_r{\Phi }_r(x),$$

$${\alpha }_1+{\alpha }_2,...,+{\alpha }_r=1.$$

$${\Phi }_{\alpha }(x)=(1-\alpha )\Vert x{\Vert }_2^2+\alpha \Vert x{\Vert }_1,$$

$${\Phi }_{\alpha }(x)=\underset{x=x_1+...+x_r}{\min }\{{\alpha }_1{\Phi }_1(x_1)+...+{\alpha }_r{\Phi }_r(x_r)\}.$$

$${\Phi }_{\alpha }(x)=\underset{z}{\min }\{(1-\alpha ){\Phi }_1(x-z)+\alpha {\Phi }_2(z)\},$$

$${\Phi }_{\alpha }(x)=\underset{z}{\min }\{(1-\alpha )\Vert x-z{\Vert }_2^2+\alpha \Vert z{\Vert }_1\},$$

$$z=\max \left(1-\frac{\alpha }{2(1-\alpha )|x|},0\right)\circ x,$$

$${\Phi }_{\alpha }(x)=\{\begin{array}{l}(1-\alpha )|x{|}^2\text{if}|x|\le \frac{\alpha }{2(1-\alpha )}\\ \alpha |x|-\frac{{\alpha }^2}{4(1-\alpha )}\text{if}|x|>\frac{\alpha }{2(1-\alpha )}\end{array},$$

$${\Phi }_{\alpha }^{\text{TT}}(x)=\underset{x=x_1+x_2}{\min }(1-\alpha )\Vert {\nabla }^2{x}_2{\Vert }_2^2+\alpha \Vert \nabla x_1{\Vert }_1,$$

$${\Phi }_{\alpha }^{\text{TGV}}(x)=\underset{x=x_1+x_2}{\min }(1-\alpha )\Vert {\nabla }^2{x}_2{\Vert }_1+\alpha \Vert \nabla x_1{\Vert }_1,$$

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

$$\mathbb{L}=\frac{\partial \mathbb{A}(m)}{\partial m}u={\omega }^2\mathbb{C}\text{diag}(\mathbb{B}u),$$

$$\underset{\begin{array}{c}m=m_1+m_2\\ m\in \mathcal{C}\end{array}}{\mathrm{argmin}}{\Phi }_{\alpha }^{\text{TT}}(m_1,m_2)+\frac{\gamma }{2}\Vert \mathbb{L}[m_1+m_2]-y{\Vert }_2^2,$$

Full text

I Introduction

Full waveform inversion (FWI) seeks to estimate constitutive parameters of a medium with a wavelength-scale resolution by minimization of a distance between the recorded and simulated data subject to the wave equation is satisfied. This technology was originally developed in geophysical imaging [1], and has spread more recently in other fields of imaging sciences such as medical imaging [2] or oceanography [3]. FWI is a partial-differential equation (PDE) - constrained nonlinear optimization problem, which is classically solved with local reduced-space optimization methods [4]. In this linearized framework, a challenging source of non linearity is the so-called cycle skipping pathology. When FWI is solved with a classical least-squares difference-based distance [5, 6], cycle skipping drives FWI to spurious local minima as soon as the initial parameter model does not allow to match the data with a kinematic error smaller than half a period. 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 [7] rely on the smoothness assumption and hence fail to recover sharp interfaces of such media. Conversely, total variation (TV) regularizations [8] are based on blockiness assumption and hence are more suitable to image large contrasts. However, they generate undesirable staircase imprint in smooth regions. Regions characterized by smoothly-varying properties and those containing sharp contrasts have indeed different statistical properties. The former are characterized by the normal prior, while the latter by a heavy tailed prior. 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 [9, 10]. A naive approach consists of the simple additive coupling or convex combinations (CC) of regularizations. Alternatively, [11] proposed to decompose explicitly 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 such a compound regularizations based upon infimal convolution (IC) outperforms the one based upon additive coupling [12].

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) [13]. IR-WRI extends the search space of FWI and decreases cycle skipping accordingly by relaxation of the wave-equation constraint [14, 15, 16]. 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. The linearity of the parameter-estimation subproblem provides a suitable framework to implement sophisticated nonsmooth regularizations.

In this study, in the following of [11] and [13], 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 [13], we jointly update the blocky and the smooth components through a variable projection process rather than in 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 $i$th component of the column vector $\mathbb{x}$ is shown by $\mathbb{x}_{i}$ and its absolute value is returned by $|\mathbb{x}_{i}|$. For the real-valued $n$-length column vectors $\mathbb{x}$ and $\mathbb{y}$ the dot product is defined by $\langle\mathbb{x},\mathbb{y}\rangle=\mathbb{x}^{T}\mathbb{y}=\sum_{i=1}^{n}\mathbb{x}_{i}\mathbb{y}_{i}$ and their Hadamard product, denoted by $\mathbb{x}\circ\mathbb{y}$, is another vector made up of their component-wise products, i.e. $(\mathbb{x}\circ\mathbb{y})_{i}=\mathbb{x}_{i}\mathbb{y}_{i}$. The $\ell_{1}$- and $\ell_{2}$-norms of $\mathbb{x}$ are, respectively, defined by $\|\mathbb{x}\|_{2}=\sqrt{\langle\mathbb{x},\mathbb{x}\rangle}=\sqrt{\sum_{i=1}^{n}\mathbb{x}_{i}^{2}}$ and $\|\mathbb{x}\|_{1}=\sum_{i=1}^{n}|\mathbb{x}_{i}|$. For a real-valued $n\times m$ matrix $\mathbb{A}$, we use the (entrywise) $\ell_{1}$- and $\ell_{2}$-norms similar to vector norms as $\|\mathbb{A}\|_{2}=\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{m}\mathbb{A}_{ij}^{2}}$ and $\|\mathbb{A}\|_{1}={\sum_{i=1}^{n}\sum_{j=1}^{m}|\mathbb{A}_{ij}|}$.

III Method

In this section, we briefly review the frequency-domain FWI as a bi-convex feasibility problem and describe the reduced and extended forms of FWI. We show how the latter can be solved with the alternating direction method of multipliers (ADMM) [17] for a general regularization functional.

III-A Full-waveform inversion

Frequency-domain FWI with a general regularization term and bounding constraints can be formulated as [15, 18]

[TABLE]

where $\mathbb{m}\in\mathbb{R}^{n\times 1}$ gathers unknown squared slowness, $\Phi(\mathbb{m})$ is an appropriate regularization term which we assume to be convex, $\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}}$.

The first constraint in (1), $\mathbb{A(m)u}=\mathbb{b}$, is a partial-differential equation (PDE) wherein $\mathbb{u}\in\mathbb{C}^{n\times 1}$ is the wavefield and $\mathbb{b}\in\mathbb{C}^{n\times 1}$ is the source term. In this study, $\mathbb{A(m)}\in\mathbb{C}^{n\times n}$ is the discretized PDE Helmholtz operator [19, 20] given by

[TABLE]

with $\omega$ the angular frequency and $\mathbb{\Delta}$ the discretized Laplace operator. $\mathbb{C}$ embeds boundary conditions and can be dependent or independent on $\mathbb{m}$ depending on the kinds of absorbing boundary conditions (radiation versus sponge) [15]. Also, $\mathbb{B}$ is 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 [21, 22, 23].

The second constraint in (1), $\mathbb{Pu}=\mathbb{d}$, is the observation equation, in which $\mathbb{d}\in\mathbb{C}^{m\times 1}$ is the recorded seismic data and $\mathbb{P}\in\mathbb{R}^{m\times n}$ is a linear operator that samples the wavefield at the receiver positions.

III-A1 Reduced approach to solving (1)

The reduced approach, which is more commonly use 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 [19, 24, 5]

[TABLE]

where $\lambda>0$ is the penalty parameter. A number of methods have been proposed to solve the optimization problems of the form (3), either for unregularized form $\Phi(\mathbb{m})=0$ [19], or the regularized form [25, 26, 27]. Although this reduced approach is more computationally tractable than the full-space approach, the highly-oscillating nature of the inverse PDE operator $\mathbb{A}^{-1}$ makes the inverse problem highly nonlinear, and hence prone to convergence to a spurious local minima when the initial $\mathbb{m}$ is not accurate enough [28, 5]. The extended approach described below is an alternative way which is more immune to local minima.

III-A2 WRI approach to solving (1)

The extended approach, known as wavefield reconstruction inversion (WRI) [14], recasts the constrained optimization problem, equation (1), as an unconstrained problem where both constraints are implemented with quadratic penalty functions.

[TABLE]

where $\lambda,\gamma>0$ are the penalty parameters. For the unregularized case where $\Phi(\mathbb{m})=0$ and also without the bounding constraint, [14] solved this biconvex minimization problem with an alternating-direction algorithm, whereby the joint minimization over $\mathbb{u}$ and $\mathbb{m}$ is replaced by an alternating minimization over each variable separately. The main property of the penalty formulation given by equation (4) 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 [14]. Its main drawback, however, is the difficulty related to the adaptive tuning of the penalty parameter, which is common to all penalty methods [29].

III-A3 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 [30, 29].

[TABLE]

where $\mathbb{v}_{1}$ and $\mathbb{v}_{2}$ are the dual variables (the Lagrangian multipliers). The min max problem (5) can also be written in a more compact form (the scaled form AL) as

[TABLE]

Applying a gradient ascent to (6) with respect to the duals, after a simple change of variables $\mathbb{d}^{k}=-\mathbb{v}^{k}_{1}/\lambda$ and $\mathbb{b}^{k}=-\mathbb{v}^{k}_{2}/\gamma$, gives the following iteration:

[TABLE]

beginning with $\mathbb{d}^{0}=0$ and $\mathbb{b}^{0}=0$. Capitalizing on the bilinearity of the wave equation in $\mathbb{m}$ and $\mathbb{u}$, ADMM [17] is an efficient method to solve this kind of multivariate optimization problem. ADMM updates $\mathbb{m}$ and $\mathbb{u}$ separately through a Gauss-Seidel like iteration, i.e., fixing $\mathbb{m}$ and solving for $\mathbb{u}$ and vice versa. Accordingly, beginning with an initial guess $\mathbb{m}^{0}$, we end up with the following iteration to solve (7) [15, 31, 18]:

[TABLE]

The subproblem (8a) associated with the wavefield reconstruction is quadratic and admits a closed-form solution. It relaxes the requirement to satisfy exactly the wave equation ($\mathbb{A}(\mathbb{m}^{k})\mathbb{u}=\mathbb{b}$) for the benefit of an improved data fitting ($\mathbb{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 exactly satisfy the wave equation, $\mathbb{u}_{r}=\mathbb{A}(\mathbb{m}^{k})^{-1}\mathbb{b}$, which makes the classical FWI highly non-convex [15].

Equation (7) 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 [29].

In the next section we focus on the solution of the model subproblem (8b) when compound regularizations are used as the regularization term.

IV The model Subproblem

This section presents different forms of the regularization functional including simple and compound functionals and the details of our approach to solve the model subproblem (8b) with these functionals.

IV-A Simple regularizers

The two most widely used regularizers rely on the (squared) $\ell_{2}$ and $\ell_{1}$-norms. The squared $\ell_{2}$-norm, defined as

[TABLE]

promotes smooth reconstruction, since the minimization of the squared value of components will penalize large components more severely than small ones.

In contrast, the $\ell_{1}$-norm, defined as

[TABLE]

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 $\ell_{1}$- or $\ell_{2}$-norms of the first and/or second order differences of the model. The first order forward differences in $x$- and $z$-direction are denotes by $\nabla_{x}f$ and $\nabla_{z}f$, with

[TABLE]

with appropriate boundary conditions, where $i$ and $j$ run over the domain of the model parameters. Accordingly, the discrete first order operator is defined as $\nabla=\begin{bmatrix}\nabla_{x}^{T}&\nabla_{z}^{T}\end{bmatrix}^{T}$ with

[TABLE]

The $\ell_{2}$-norm of $|\nabla f|$ gives the first order Tikhonov regularization [32], which returns a flat regularized model (with a small gradient), while its $\ell_{1}$-norm gives the total variation regularization [8], which returns a piecewise constant model (with a sparse gradient).

Analogously, the second order forward differences are denoted by $\nabla_{xx}f$ and $\nabla_{zz}f$, with

[TABLE]

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 $\nabla^{2}=\begin{bmatrix}\nabla_{xx}^{T}&\nabla_{zz}^{T}\end{bmatrix}^{T}$ with

[TABLE]

The $\ell_{2}$-norm of $|\nabla^{2}f|$ gives the second-order Tikhonov regularization, which returns a smooth regularized model (with a small Laplacian), while its $\ell_{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 $\nabla_{xz}f\equiv\nabla_{z}\nabla_{x}f$ with

[TABLE]

A discrete second-order operator, which includes mixed differences is defined as $\nabla^{2}=\begin{bmatrix}\nabla_{xx}^{T}&\sqrt{2}\nabla_{xz}^{T}&\nabla_{zz}^{T}\end{bmatrix}^{T}$ with

[TABLE]

which equals the Frobenius norm of the Hessian matrix [33, 34, 35].

Simple regularizations are effective for recovering models which can be characterize by a single prior. In the next section we propose compound regularizers which are more effective for recovering complicated models.

IV-B Compound regularizers

Compound regularizers 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).

IV-B1 Convex combination of simple regularizers

A CC of $r$ simple regularizer functionals $\Phi_{1},...,\Phi_{r}$ is a compound regularizer functional of the form

[TABLE]

where weights $\alpha_{i}$ satisfy $\alpha_{i}\geq 0$ and

[TABLE]

Definitely, if all of the functions $\Phi_{1},...,\Phi_{r}$ are convex then $\Phi$ 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 $\ell_{1}$- and squared $\ell_{2}$-norms ($\ell_{1}+\ell_{2}$), which known as elastic net [36, 37], is

[TABLE]

with $0\leq\alpha\leq 1$. The convexity of $\ell_{1}$- and $\ell_{2}$-norms implies that $\Phi_{\alpha}(\mathbb{x})$ in (19) is convex. One may also construct a compound regularizer functional by a CC of two $\ell_{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 [9].

IV-B2 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 $\Phi_{1},...,\Phi_{r}$ is a compound functional of the form

[TABLE]

In the case of two functionals, $\Phi_{\alpha}$ in (20) takes the form

[TABLE]

which is similar to the classical formula of convolution, and hence the term infimal convolution.

The IC of $\ell_{1}$- and $\ell_{2}$-norms ($\ell_{1}\oplus\ell_{2}$) is

[TABLE]

which is a denoising problem and the solution of which is unique and is obtained simply by the well-known soft-threshold function [38]:

[TABLE]

Putting $\mathbb{z}$ from (23) into (22) gives that

[TABLE]

which is nothing other than the Huber function [39]. As seen, this function has a hybrid behavior: it has a quadratic behavior for small values of $|\mathbb{x}|$ and linear behavior for large values. The parameter $\frac{\alpha}{2(1-\alpha)}$ determines where the transition from quadratic to linear behavior takes place.

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

In this paper, we consider (21) in the following settings, though other configurations are possible:

[TABLE]

and

[TABLE]

where in both (25) and (26), the norms are applied on the absolute valued components of $\nabla^{2}\mathbb{x}_{2}$ ((14) and (16)) and $\nabla\mathbb{x}_{1}$ (12). The compound regularizer $\Phi^{\text{TT}}_{\alpha}$ is a combination of the second order Tikhonov and TV (TT) regularizations [11] and $\Phi^{\text{TGV}}_{\alpha}$ is a combination of the first and second order TV regularizations, called total generalized variation (TGV) [40, 41]. The former is suitable for recovering piecewise-smooth models, while the latter is better suited for piecewise linear models. Next section gives a solution procedure to solve (8b) with these regularizers.

IV-C Solving the subproblem (8b)

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

[TABLE]

where

[TABLE]

where we assume that $\mathbb{C}$ does not depend on $\mathbb{m}$ (this is the case for perfectly-matched absorbing boundary conditions [15]). From the explicit decomposition $\mathbb{m}=\mathbb{m}_{1}+\mathbb{m}_{2}$ and (27), the solution of the optimization problem (8b) can be expressed as

[TABLE]

where $\mathbb{y}=\mathbb{b+b}^{k}-\Delta\mathbb{u}^{k+1}$. Applying ADMM to (29) leads to the following iteration

[TABLE]

where

[TABLE]

and $\mathbb{p}=\nabla\mathbb{m}_{1}$ and $\mathbb{q}$ are auxiliary primal variables that are introduced to decouple the $\ell_{1}$ and the $\ell_{2}$ minimization problems and solve the former ones with proximal algorithms following the split Bregman method [42]. The dual variables $\tilde{\mathbb{g}}$, $\tilde{\mathbb{m}}$ are updated through a gradient ascent step according to the method of multipliers [29]

[TABLE]

We now discuss how to solve the subsubproblems given in (30).

IV-C1 The subsubproblem (30a)

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

[TABLE]

with

[TABLE]

and

[TABLE]

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

[13] broke down the $2n\times 2n$ problem (33) into two smaller $n\times n$ systems and updates $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$ in alternating mode at the expense of convergence speed [13, 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 (33), we find that

[TABLE]

and plugging this into the second equation of (33) we get the following:

[TABLE]

Interestingly, $\mathbb{L}$ is diagonal, implying that $\mathbb{G}_{12}$ is also diagonal. Thus we only need to solve an $n\times n$ system to estimate $\mathbb{m}_{1}$, from which $\mathbb{m}_{2}$ easily follows.

IV-C2 The subsubproblem (30b)

The sub-problem for $\mathbb{p}$, equation (30b), is a denoising problem and is straightforward to solve. Note that $\mathbb{p}$ has two components associated with the gradient in each direction:

[TABLE]

Equation (30b) is solved with a generalized proximity operator [43] leading to

[TABLE]

where

[TABLE]

and

[TABLE]

IV-C3 The subsubproblem (30c)

The optimization problem (30c) has also a entrywise solution given by

[TABLE]

where the projection operator projects its argument onto the desired box $[\mathbb{m}_{l},\mathbb{m}_{u}]$ according to $\text{proj}_{\mathcal{C}}(\bullet)=\min(\max(\bullet,\mathbb{m}_{l}),\mathbb{m}_{u})$.

It should be noted that the total algorithm consists of two levels of iterations: an outer iteration given in (8) and an inner iteration given in (30) corresponding to the model subproblem (8b). Numerical results, however, show that only one inner iteration suffices for convergence of the algorithm, hence significantly reducing the total computational cost.

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 [44] when a crude initial model and realistic frequencies are used as starting points.

V-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 [44], Marmousi II [45], SEG/EAGE overthrust [46], SEG/EAGE salt [46] and synthetic Valhall [47] 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 a single frequency 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 3-point 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 proposed by [18] for bound-constrained TV regularized IR-WRI. Moreover, for a fair comparison of the compound regularizers (JTT, TT and TGV), we select for each of them the optimum value of $\alpha$ among a range of preset values that minimizes the error in the models estimated by the IR-WRI. Also, we set the parameter bounds $\mathbb{m}_{l}$ and $\mathbb{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 using a maximum number of iteration equal to 100 as stopping criterion. The average mean error of the estimated velocity profiles for the five benchmark models and the different regularizations are plotted in Fig. 2. 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-up of the reconstructed profiles in Fig. 3 as well as the profiles reconstructed the DMP regularization in the first column. 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.

V-B 2004 BP salt model

We now consider a more realistic application with a target of the challenging 2004 BP salt model [44]. 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 with PML boundary conditions using a 10 Hz Ricker wavelet as source signature. A line of receivers with a 25 m spacing are deployed at the surface leading to a 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. 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, meaning that we initiate the inversion from scratch (Fig. 4b). We start with inverting the first batch of frequencies ({3, 3.5} Hz) with noiseless data using a maximum number of iteration equal to 45 as stopping criterion of iteration. 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 and 5c, 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 are 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. 5(d-f) and Fig. 6b. Clearly, the TT regularizer better captures the long wavelengths of the salt body and the smooth subsalt background model. The joint evolution in iterations of the observation-equation ($\|\mathbb{Pu-d}\|_{2}$) and wave-equation errors ($\|\mathbb{A(m)u-b}\|_{2}$) (Fig. 7a) and the relative model errors in iterations (Fig. 7b) further confirm the relative performances 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. 7(a-b). As already highlighted by [15], this results from the dynamic balancing in iterations of the observation-equation 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

[TABLE]

where $k_{max}$ denotes the maximum iteration count, $\varepsilon_{b}$=1e-3, and $\varepsilon_{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], [6, 13] 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. The number of iterations that have been performed with DMP regularization (Fig. 8a) is 426. Also, the number of iterations with Tikhonov (Fig. 8b) and TV regularizations (Fig. 8c) are 448 and 399, respectively, while it is 415 for the JTT regularization (Fig. 8d). Also, the number of iterations that have been performed with TT (Fig. 8e) and TGV (Fig. 8f) are 361 and 394, respectively. 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 large-contrast 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. Moreover, the TT has a better convergence rate compare to the TGV. It is also worth noting the significant differences between the JTT and TT IR-WRI models in particular in the deep part of the model.

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 and 11). The real part of wave-equation error (Fig 10) and data residuals (Fig 11) are plotted at the first and final iterations of the inversion. Both the final data and source residuals suggest that the TT regularizer slightly outperforms the TGV counterpart at low frequencies.

VI Conclusions

In this study, we first show how to implement efficiently 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 reconstruct accurately large-contrast subsurface media when starting from very crude initial models. This compound regularization is suitable for seismic imaging of the subsurface as this later can be often represented by piecewise smooth media. We show that the infimal convolution (IC) of the Tikhonov and TV regularizers captures much more efficiently the blocky and smooth components 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 infimal-convolution 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 search-space 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

This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by AKERBP, CGG, CHEVRON, EQUINOR, EXXON-MOBIL, JGI, PETROBRAS, 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.”