Eigenvector Model Descriptors for Solving an Inverse Problem of Helmholtz Equation: Extended Materials

TL;DR

This paper introduces an eigenvector-based regularization method for seismic inverse problems, improving subsurface property recovery by reducing ill-posedness and effectively reconstructing complex media like salt-domes.

Contribution

The study proposes a novel eigenvector model descriptor approach for seismic inversion, demonstrating its effectiveness in 2D and 3D reconstructions, especially for challenging salt-dome media.

Findings

Eigenvector representation enhances seismic imaging accuracy.

The method compensates for low-frequency data limitations.

Effective in 2D and 3D seismic inverse problems.

Abstract

We study the seismic inverse problem for the recovery of subsurface properties in acoustic media. In order to reduce the ill-posedness of the problem, the heterogeneous wave speed parameter to be recovered is represented using a limited number of coefficients associated with a basis of eigenvectors of a diffusion equation, following the regularization by discretization approach. We compare several choices for the diffusion coefficient in the partial differential equations, which are extracted from the field of image processing. We first investigate their efficiency for image decomposition (accuracy of the representation with respect to the number of variables and denoising). Next, we implement the method in the quantitative reconstruction procedure for seismic imaging, following the Full Waveform Inversion method, where the difficulty resides in that the basis is defined from an initial model where none of the actual structures is known. In particular, we demonstrate that the method is efficient for the challenging reconstruction of media with salt-domes. We employ the method in two and three-dimensional experiments and show that the eigenvector representation compensates for the lack of low frequency information, it eventually serves us to extract guidelines for the implementation of the method.

Figures

Tables

Table 1. Table 1: List of formula for the coefficient of the diffusion operator which is used to decompose the image.

reference (name)	definition	$\beta \to 0$	$\beta \to \infty$
Perona & Malik (1988, 1990)	${\eta }_1(m,\beta )=\frac{\beta }{\beta +g_2\left(m\right)}$	0	1
Perona & Malik (1988, 1990)	${\eta }_2(m,\beta )=\exp \left(-\frac{g_2\left(m\right)}{\beta }\right)$	0	1
Geman & Reynolds (1992)	${\eta }_3(m,\beta )=\frac{2\beta }{{\left(\beta +g_2\left(m\right)\right)}^2}$	0	0
Green (1990)	${\eta }_4(m,\beta )=\tanh \left(\frac{g_1\left(m\right)}{\beta }\right)\left(\frac{1}{\beta g_1\left(m\right)}\right)$	$+\infty$	0
Charbonnier et al. (1994)	${\eta }_5(m,\beta )=\frac{1}{\beta }{\left(\frac{\beta +g_2\left(m\right)}{\beta }\right)}^{-1/2}$	$+\infty$	0
Grote & Nahum (2019)
(Lorentzian)	${\eta }_6(m,\beta )=\frac{\beta }{{\left(1+\beta g_2\left(m\right)\right)}^2}$	0	0
Grote & Nahum (2019)
(Gaussian)	${\eta }_7(m,\beta )={\left(\beta \exp \left(\frac{g_2\left(m\right)}{\beta }\right)\right)}^{-1}$	0	0
Rudin et al. (1992)
(Total Variation, TV)	${\eta }_8\left(m\right)=\frac{1}{g_1\left(m\right)}$	n/a	n/a
Tikhonov	${\eta }_9=1$ $\frac{1}{2}$	n/a	n/a

Table 2. Table 2: Minimal relative error obtained and associated scaling coefficient: ℰ ℰ \mathcal{E} ( β 𝛽 \beta ), for the decomposition of the Marmousi model Figure  3(a) . The definition of η 𝜂 \eta is given Table  1 .

Coeff.	$N=10$	$N=20$	$N=50$	$N=100$	$N=250$	$N=500$
${\eta }_1$	6% (${10}^{-6}$)	5% (${10}^{-6}$)	4% (${10}^{-6}$)	4% (${10}^{-6}$)	3% (${10}^{-6}$)	3% (${10}^{-7}$)
${\eta }_2$	14% (5.${10}^{-2}$)	13% (5.${10}^{-2}$)	12% (5.${10}^{-2}$)	9% (${10}^{-2}$)	7% (${10}^{-2}$)	5% (${10}^{-2}$)
${\eta }_3$	8% (${10}^{-4}$)	7% (${10}^{-3}$)	6% (${10}^{-3}$)	5% (${10}^{-3}$)	5% (${10}^{-3}$)	4% (${10}^{-2}$)
${\eta }_4$	14% (5.${10}^{-1}$)	14% (${10}^1$)	13% (${10}^1$)	13% (${10}^{-1}$)	12% (5.${10}^{-1}$)	10% (5.${10}^1$)
${\eta }_5$	13% (${10}^{-7}$)	12% (${10}^{-6}$)	12% (${10}^{-5}$)	10% (${10}^{-5}$)	10% (${10}^{-5}$)	6% (${10}^{-7}$)
${\eta }_6$	8% (5.${10}^3$)	7% (${10}^3$)	6% (5.${10}^2$)	5% (${10}^3$)	5% (${10}^3$)	4% (${10}^2$)
${\eta }_7$	14% (5.${10}^{-2}$)	13% (5.${10}^{-2}$)	12% (5.${10}^{-2}$)	11% (5.${10}^{-2}$)	9% (5.${10}^{-2}$)	7% (5.${10}^{-2}$)
${\eta }_8$	15% (n/a)	14% (n/a)	13% (n/a)	12% (n/a)	10% (n/a)	9% (n/a)
${\eta }_9$	14% (n/a)	14% (n/a)	14% (n/a)	13% (n/a)	12% (n/a)	11% (n/a)

Table 3. Table 3: Minimal relative error obtained and associated scaling coefficient: ℰ ℰ \mathcal{E} ( β 𝛽 \beta ), for the decomposition of the salt model Figure  3(b) . The definition of η 𝜂 \eta is given Table  1 .

Coeff.	$N=10$	$N=20$	$N=50$	$N=100$	$N=250$	$N=500$
${\eta }_1$	4% (${10}^{-3}$)	4% (${10}^{-3}$)	3% (${10}^{-2}$)	2% (${10}^{-2}$)	1% (${10}^{-2}$)	1% (${10}^{-2}$)
${\eta }_2$	9% (${10}^{-1}$)	7% (${10}^{-1}$)	5% (${10}^{-1}$)	5% (${10}^{-1}$)	4% (${10}^{-1}$)	3% (${10}^{-1}$)
${\eta }_3$	8% (5.${10}^{-2}$)	4%(5.${10}^{-2}$)	3% (5.${10}^{-2}$)	3% (5.${10}^{-2}$)	2% (5.${10}^{-2}$)	1% (5.${10}^{-2}$)
${\eta }_4$	14% (${10}^{-5}$)	12%(5.${10}^1$)	9% (5.${10}^1$)	8% (5.${10}^1$)	6% (${10}^1$)	3% (5.${10}^1$)
${\eta }_5$	3% (${10}^{-6}$)	3% (${10}^{-6}$)	2% (${10}^{-5}$)	1% (${10}^{-5}$)	1% (${10}^{-5}$)	1% (${10}^{-4}$)
${\eta }_6$	6% (50)	5% (${10}^1$)	4% (${10}^1$)	3% (${10}^1$)	2% (${10}^1$)	1% (${10}^1$)
${\eta }_7$	9% (${10}^{-1}$)	7% (${10}^{-1}$)	5% (${10}^{-1}$)	5% (${10}^{-1}$)	3% (${10}^{-1}$)	3% (${10}^{-1}$)
${\eta }_8$	59% (n/a)	22% (n/a)	16% (n/a)	13% (n/a)	11% (n/a)	7% (n/a)
${\eta }_9$	20% (n/a)	15% (n/a)	13% (n/a)	10% (n/a)	7% (n/a)	5% (n/a)

Table 4. Table 4: Minimal relative error obtained and associated scaling coefficient: ℰ ℰ \mathcal{E} ( β 𝛽 \beta ). The error is computed with respect to the noiseless model Figure  3(a) but the decomposition uses the noisy model of Figure  7(a) . The definition of η 𝜂 \eta is given Table  1 .

Coeff.	$N=10$	$N=20$	$N=50$	$N=100$	$N=250$
${\eta }_1$	16% (${10}^{-4}$)	15% (${10}^{-4}$)	14% (${10}^{-7}$)	13% (${10}^{-7}$)	11% (${10}^{-7}$)
${\eta }_2$	16% (5.${10}^{-2}$)	16% (5.${10}^{-2}$)	15% (5.${10}^{-2}$)	14% (5.${10}^{-2}$)	13% (5.${10}^{-2}$)
${\eta }_3$	16% (${10}^{-2}$)	14% (${10}^{-3}$)	13% (${10}^{-4}$)	12% (${10}^{-4}$)	10% (${10}^{-4}$)
${\eta }_4$	16% (${10}^1$)	16% (${10}^1$)	16% (5.${10}^1$)	15% (5.${10}^1$)	14% (${10}^2$)
${\eta }_5$	16% (${10}^{-6}$)	16% (${10}^{-7}$)	15% (${10}^{-7}$)	15% (${10}^{-7}$)	14% (${10}^{-6}$)
${\eta }_6$	16% (${10}^2$)	14% (${10}^3$)	13% (${10}^4$)	12% (${10}^4$)	10% (${10}^4$)
${\eta }_7$	16% (5.${10}^{-2}$)	16% (5.${10}^{-2}$)	15% (5.${10}^{-2}$)	14% (5.${10}^{-2}$)	13% (5.${10}^{-2}$)
${\eta }_8$	17% (n/a)	16% (n/a)	16% (n/a)	15% (n/a)	14% (n/a)
${\eta }_9$	17% (n/a)	17% (n/a)	16% (n/a)	16% (n/a)	14% (n/a)

Table 5. Table 5: Minimal relative error obtained and associated scaling coefficient: ℰ ℰ \mathcal{E} ( β 𝛽 \beta ). The error is computed with respect to the noiseless model Figure  3(b) but the decomposition uses the noisy model of Figure  7(b) . The definition of η 𝜂 \eta is given Table  1 .

Coeff.	$N=10$	$N=20$	$N=50$	$N=100$	$N=250$
${\eta }_1$	14% (${10}^{-6}$)	14% (${10}^{-6}$)	11% (${10}^{-5}$)	9% (${10}^{-5}$)	5% (${10}^{-3}$)
${\eta }_2$	17% (5.${10}^{-2}$)	15% (5.${10}^3$)	11% (5.${10}^{-2}$)	8% (5.${10}^{-2}$)	6% (5.${10}^{-2}$)
${\eta }_3$	10% (${10}^{-4}$)	10% (${10}^{-4}$)	8% (${10}^{-4}$)	6% (${10}^{-4}$)	5% (${10}^{-2}$)
${\eta }_4$	18% (${10}^{-3}$)	15% (${10}^0$)	12% (${10}^{-1}$)	10% (${10}^{-3}$)	6% (${10}^{-3}$)
${\eta }_5$	17% (${10}^{-7}$)	15% (5.${10}^3$)	12% (${10}^{-4}$)	9% (${10}^{-7}$)	6% (${10}^{-5}$)
${\eta }_6$	10% (${10}^4$)	9% (${10}^4$)	8% (${10}^4$)	6% (5.${10}^3$)	5% (${10}^2$)
${\eta }_7$	17% (5.${10}^{-2}$)	15% (${10}^5$)	11% (${10}^{-2}$)	8% (5.${10}^{-2}$)	6% (5.${10}^{-2}$)
${\eta }_8$	18% (n/a)	16% (n/a)	12% (n/a)	9% (n/a)	6% (n/a)
${\eta }_9$	21% (n/a)	15% (n/a)	13% (n/a)	11% (n/a)	8% (n/a)

Table 6. Table 6: Minimal relative error obtained and associated scaling coefficient: ℰ ℰ \mathcal{E} ( β 𝛽 \beta ), for the decomposition of the SEAM model Figure  10(a) . The definition of η 𝜂 \eta is given Table  1 .

Coeff.	$N=10$	$N=20$	$N=50$	$N=100$	$N=200$
${\eta }_1$	4% (${10}^{-6}$)	3% (${10}^{-6}$)	3% (${10}^{-6}$)	2% (${10}^{-6}$)	2% (${10}^{-7}$)
${\eta }_2$	8% (5.${10}^{-2}$)	6% (5.${10}^{-2}$)	5% ( ${10}^{-1}$)	3% ( ${10}^{-2}$)	3% ( ${10}^{-2}$)
${\eta }_3$	6% (${10}^{-3}$)	5% (${10}^{-2}$)	4% (${10}^{-3}$)	3% (${10}^{-2}$)	2% (${10}^{-3}$)
${\eta }_4$	10% ( 5)	9% ( 10)	7% (${10}^2$)	6% (${10}^2$)	5% ( 5)
${\eta }_5$	7% (${10}^{-7}$)	6% (${10}^{-6}$)	4% (${10}^{-6}$)	3% (${10}^{-5}$)	3% (${10}^{-5}$)
${\eta }_6$	6% (${10}^3$)	5% (${10}^2$)	4% (${10}^3$)	3% (${10}^2$)	2% (${10}^3$)
${\eta }_7$	8% (${5.10}^{-2}$)	6% (${5.10}^{-2}$)	5% (${10}^{-1}$)	4% (${5.10}^{-2}$)	3% (${5.10}^{-2}$)
${\eta }_8$	18% (n/a)	17% (n/a)	14% (n/a)	12% (n/a)	9% (n/a)
${\eta }_9$	12% (n/a)	10% (n/a)	7% (n/a)	7% (n/a)	6% (n/a)

Table 7. Table 8: List of experiments for the reconstruction of the two-dimensional salt dome model (Figure  13(b) ). For each combination of frequency and associated number of eigenvectors N 𝑁 N in the decomposition, 30 30 30 iterations are performed ( n iter subscript 𝑛 iter n_{\text{iter}} in Algorithm  2 ).

	frequency list	list of $N$	total iterations	Reference
Experiment 1	$2$ Hz	50	180	Figures 15, B1
Experiment 2	$2$ Hz	$\{50,60,70,80,90,100\}$	180	Figures 17, B2
Experiment 3	$\{2,3,4,5\}$Hz	50	120	Figures 17, B3
Experiment 4	$\{2,3,4,5\}$Hz	$\{50,60,70,80\}$	120	Figures 17, B4

Equations

$$\Big{(}-\Delta-\omega^{2}c^{-2}(\boldsymbol{x})\Big{)}p(\boldsymbol{x})=f(\boldsymbol{x}),\quad\quad\text{in $\Omega$},$$

$$\left\{\begin{aligned} \big{(}-\Delta-\omega^{2}c^{-2}(\boldsymbol{x})\big{)}p(\boldsymbol{x})&=f(\boldsymbol{x}),\quad\quad&\text{in $\Omega$},\\ p(\boldsymbol{x})&=0,\quad\quad&\text{on $\Gamma_{1}$},\\ \partial_{\nu}p(\boldsymbol{x})-\mathrm{i}\omega c^{-1}(\boldsymbol{x})p(\boldsymbol{x})&=0,\quad\quad&\text{on $\Gamma_{2}$},\\ \end{aligned}\right.$$

$${\mathcal{F}}_{\omega }:m\to p(\mathbf{x}){\mid }_{\Sigma }.$$

$${\mathbf{y}}_{\omega }={\mathcal{F}}_{\omega }(m_{†})+{\mathfrak{E}}_{\omega },$$

$$\underset{m}{\min }\mathcal{J}(m)=\frac{1}{2}\sum _{\omega }\Vert {\mathcal{F}}_{\omega }(m)-{\mathbf{y}}_{\omega }{\Vert }^2,$$

$$m^{(k+1)}=m^{(k)}+\mu s^{(k)},k>0.$$

$$\underset{m}{\min }{\mathcal{J}}_r(m)=\frac{1}{2}\sum _{\omega }\Vert {\mathcal{F}}_{\omega }(m)-{\mathbf{y}}_{\omega }{\Vert }^2+\mathcal{I}(m),$$

$$\mathcal{I}(m)={\int }_{\Omega }\varphi (|\nabla m|)\mathrm{d}\Omega ,$$

$$\nabla\cdot\bigg{(}\dfrac{\phi^{\prime}(|\nabla m|)}{|\nabla m|}\nabla m\bigg{)}=0\quad\quad\text{in $\Omega$.}$$

$$\mathcal{A}(\mathrm{y},\,\eta):=-\nabla\cdot\big{(}\eta(\mathrm{y})\nabla\big{)},\quad\quad\quad\text{with}\quad\quad\eta(\mathrm{y})=\dfrac{\phi^{\prime}(|\nabla\mathrm{y}|)}{|\nabla\mathrm{y}|}.$$

$$g_1(m,\mathbf{x})=\frac{|\nabla m(\mathbf{x})|}{{\gamma }_1},g_2(m,\mathbf{x})=\frac{|\nabla m(\mathbf{x}){|}^2}{{\gamma }_2},\text{with}$$

$$\displaystyle\gamma_{1}=\max\big{(}|\nabla m(\boldsymbol{x})|\big{)}\,,\quad\gamma_{2}=\max\big{(}|\nabla m(\boldsymbol{x})|^{2}\big{)}\,,\quad|\nabla m(\boldsymbol{x})|=\sqrt{\sum_{k=1}^{\mathfrak{d}}\bigg{(}\dfrac{\partial m}{\partial x_{k}}\bigg{)}^{2}},$$

$$\{\begin{array}{rlr}\mathcal{A}(m,\eta ){\mathfrak{m}}_0& =0,& \text{in Ω,}\\ {\mathfrak{m}}_0& =m,& \text{on Γ.}\end{array}$$

$$\{\begin{array}{rlr}\mathcal{A}(m,\eta ){\psi }_k& ={\lambda }_k{\psi }_k,& \text{in Ω,}\\ {\psi }_k& =0,& \text{on Γ.}\end{array}$$

$$\mathfrak{m}={\mathfrak{m}}_0+\sum _{k=1}^N{\alpha }_k{\psi }_k,$$

$$\mathbf{\psi }(m,\eta ,N)=\{{\psi }_k{\}}_{k=1,\dots ,N}$$

$$\eta$$

$${\mathfrak{m}}(m,N,{\boldsymbol{\psi}})\qquad\text{is the decomposition of the model $m$ using $N$ vectors from ${\boldsymbol{\psi}}$,~{}\eqref{eq:methodEV_decomposition}.}$$

$$\mathbf{\psi }=\mathbf{\psi }(m,\eta ,N_{\max }).$$

$${\mathfrak{m}}^{(0)}=\mathfrak{m}(m^{(0)},N_1,\mathbf{\psi })={\mathfrak{m}}_0+\sum _{l=1}^{N_1}{\alpha }_l^{(0)}{\psi }_l.$$

$${\mathfrak{m}}^{(k+1)}={\mathfrak{m}}_0+\sum _{l=1}^N{\alpha }_l^{(k+1)}{\psi }_l.$$

$$\frac{\partial \mathcal{J}}{\partial \alpha }=\frac{\partial \mathcal{J}}{\partial m}\frac{\partial m}{\partial \alpha }.$$

$$\mathcal{E}=100\frac{\Vert m-\mathfrak{m}\Vert }{\Vert m\Vert },\text{Relative Error (\%)},$$

Full text

1 Introduction

We consider the inverse problem associated with the propagation of time-harmonic waves which occurs, for example, in seismic applications, where the mechanical waves are used to probe the Earth. Following the non-intrusive geophysical setup for exploration, we work with measured seismograms that record the waves at the surface (i.e., partial boundary measurements) and one-side illumination (back-scattered/reflection data). In the last decades, this problem has encountered a growing interest with the increase in numerical capability and the use of supercomputers. However, the accurate recovery of the deep subsurface structures remains a challenge, due to the nonlinearity and ill-posedeness of the problem, the availability of partial reflection data only, and the large scale domains of investigation.

In the context of seismic, the quantitative reconstruction of physical properties using an iterative minimization of a cost function originally follows the work of Bamberger et al. (1979), Lailly (1983) and Tarantola (1984, 1987) in the time-domain, Pratt et al. (1990; 1996; 1998) for the frequency approach. The method is commonly referred to as Full Waveform Inversion (FWI), which takes the complete observed seismograms for data. One key of FWI is that the gradient of the misfit functional is computed using the adjoint-state method (Lions & Mitter, 1971; Chavent, 1974), to avoid the formation of the (large) Jacobian matrix; we refer to Plessix (2006) for a review in geophysical applications. Then, Newton-type algorithms represent the traditional framework to perform the iterative minimization. Due to the large computational scale of the domain investigated, seismic experiments may have difficulties to incorporate second order (Hessian) information in the algorithm, and alternative techniques have been proposed, for example, in the work of Pratt et al. (1998); Akcelik et al. (2002); Choi et al. (2008); Métivier et al. (2013); Jun et al. (2015). The quantitative (as opposed to qualitative) reconstruction methods based upon iterative minimization are naturally not restricted to seismic and we refer, among others, to Ammari et al. (2015); Barucq et al. (2018) and the references therein for additional applications using similar techniques.

The main difficulty of FWI (in both exploration and global seismology) lies in the high nonlinearity of the problem and the presence of local minima in the misfit functional, which are due to the time shifts and cycle-skipping effect (Bunks et al., 1995), in particular when the background velocity (the low-frequency profile) is not correctly anticipated (Gauthier et al., 1986; Luo & Schuster, 1991; Fichtner et al., 2008; Barucq et al., 2019b). For this reason, the phase information is included in the traveltime inversion by Luo & Schuster (1991) by using a cross-correlation function between measurements and simulations, where the relative phase shift is given by the maximum of the correlation. The method is further generalized by Gee & Jordan (1992), while Van Leeuwen & Mulder (2010) propose to select the phase shift using a weighted norm. The choice of misfit has further encountered a growing interest in the past decade: in application to global-scale seismology, Fichtner et al. (2008) compare the phase, correlation-based, and envelope misfit functionals, the latter being also studied by Bozdağ et al. (2011). In exploration seismic, comparisons of phase and amplitude inversion are performed by Bednar et al. (2007); Pyun et al. (2007). The $L1$ norm is studied by Brossier et al. (2010) while approaches based upon optimal transport are considered by Métivier et al. (2016); Yang et al. (2018). In the context where different fields are measured, Alessandrini et al. (2019); Faucher et al. (2019) advocate for a reciprocity-based functional, which further connects to the correlation-based formulas (Faucher et al., 2019). In the case of accurate knowledge of the background velocity, the inverse problem is close to linear or quasi-linear as the Born approximation holds and then, alternative methods of linear inverse problem can be applied, such as the Backus–Gilbert method (Backus & Gilbert, 1967, 1968). The difficulty to recover the background velocity variation has also motivated alternative parametrization of the inverse problem: for instance the MBTT (Background/Data-Space Reflectivity) reformulation of FWI (Clément et al., 2001; Barucq et al., 2019b).

In order to diminish the ill-posedness of the inverse problem, a regularization criterion can be incorporated. It introduces an additional constraint (in addition to the fidelity between observations and simulations), which, however, may be complicated to select a priori and problem dependent (with ‘tuning’ parameters). For instance, we refer to the body of work of Kirsch (1996); Isakov (2006); Kern (2016); Kaltenbacher (2018) and the references therein. In the regularization by discretization approach, the model representation plays the role of regularizing the functional, by controlling (and limiting) the number of unknowns in the problem, and possibly defining (i.e. constraining) the shape of the unknown (e.g., to force smoothness). Controlling the number of unknowns influences the resolution of the outcome, but also the stability and convergence of the procedure. The use of piecewise constant coefficients appears natural for numerical applications, and is also motivated by stability results (Alessandrini & Vessella, 2005; Beretta et al., 2016). However, such a decomposition can lead to an artificial ‘block’ representation (cf. Beretta et al. (2016); Faucher (2017)) which would not be appropriate in terms of resolution. For this reason, a piecewise linear model representation is explored by Alessandrini et al. (2018, 2019), still motivated by the stability properties. We also mention the wavelet-based model reductions, that offer a flexible framework and are used for the purpose of regularization in seismic tomography by Loris et al. (2007, 2010). In the work of Yuan et al. (2014; 2015), FWI is carried out in the time-domain with a model represented from a wavelet-based decomposition.

In our work, we will use a model decomposition based upon the eigenvectors of a chosen diffusion operator, as introduced by De Buhan & Kray (2013); Grote et al. (2017); Grote & Nahum (2019). Note that this decomposition is shown (with the right choice of operator) to be related with the more standard Total Variation (TV) or Tikhonov regularizations. The main difference in our work is that we study several alternatives for the choice of the operator following image processing techniques, which traditionally also relies on such diffusion PDEs (e.g. Weickert (1998)). We first investigate the performance of the decomposition depending on the choice of PDE, and, next, the performance of such a model decomposition as parametrization of the reconstruction procedure in seismic FWI. It shows that the efficient choice of PDE should change depending on the situation. In addition, we provide a series of experiment to extract the robust guidelines for the implementation of the method in seismic.

We specifically target the reconstruction of subsurface salt domes (i.e. media with high contrasting objects), which is particularly challenging, because (in addition to the usual restrictive data) of the change of the kinematics involved and the lack of low frequency data (Farmer et al., 1996; Chironi et al., 2006; Virieux & Operto, 2009; Barucq et al., 2019a). In such cases, the use of the Total-Variation regularization (Rudin et al., 1992) with FWI is becoming more and more prominent, and consists in incorporating an additional constraint on the model in the minimization problem. Its efficiency is shown in the context of acoustic media with salt-dome contrasts by, e.g., Brandsberg-Dahl et al. (2017); Esser et al. (2018); Kalita et al. (2019); Aghamiry et al. (2019). In our work, we study several alternatives and demonstrate that the model representation with the criterion extracted from Geman & Reynolds (1992) appears the most appropriate for the eigenvector decomposition method in the presence of salt-domes. We also show the limitations of the method, in particular, it appears that the decomposition fails to represent models which are composed of several thin structures.

In Section 2, we define the inverse problem associated with the Helmholtz equation and introduce the iterative method for the reconstruction of the wave speed. In Section 3, we review several possibilities for the model decomposition using the eigenvectors of the diffusion operators. The process of model (image) decomposition is illustrated in Section 4. Then, in Section 5, we carry out the iterative reconstruction with FWI experiments in two and three dimensions. Here, the model decomposition is based upon the initial model, which does not contain a priori information on the target, hence increasing the differences of performance depending on the selection of the basis. It allows us to identify the best candidate for the recovery of salt dome, and to extract some guidelines for applications in quantitative reconstruction.

2 Inverse time-harmonic wave problem

2.1 Forward problem

We consider a domain $\Omega$ in two or three dimensions, with $\Omega\subset\mathbb{R}^{2}$ or $\Omega\subset\mathbb{R}^{3}$. We focus on acoustic media where, for simplicity, the density is taken as a constant, leading us to the identification of a single heterogeneous parameter: the wave speed. The propagation of waves in an acoustic medium with constant density is given by the scalar pressure field $p$, solution to the Helmholtz equation

[TABLE]

where $c$ is the wave speed, $f$ the source, and $\omega$ the angular frequency. We now have to specify the boundary conditions to formulate the appropriate problem.

Following a seismic setup, the boundary of $\Omega$, $\Gamma$ is separated into two. The upper (top) boundary, $\Gamma_{1}$, represents the interface between the Earth and the air, and is subsequently represented via a free surface boundary condition, where the pressure is null. On the other hand, other part of the boundary corresponds to the numerical need for restricting the area of interest. Here, conditions must ensure that entering waves are not reflected back to the domain (i.e., because the area of interest is only a part of the Earth), see Figure 1. The two most popular formulations to handle such numerical boundary are either the Perfectly Matched Layers (PML, Bérenger (1994)), or outgoing artificial boundary conditions. In our case, we use Absorbing Boundary Conditions (ABC, Engquist & Majda (1977)) so that the complete problem writes as

[TABLE]

where $\partial_{\nu}$ is the normal derivative. We recall that $p$ denotes the pressure, $\omega$ denotes the frequency of the induced source and $c$ denotes the wave speed.

The inverse problem aims the recovery of the wave speed $c$ in (2), from a discrete set of measurements (i.e. partial data), which corresponds to observation of the wave propagations. More precisely, our data consist in measurements of the pressure field solution to (2), at the (discrete) device locations. We refer to $\Sigma$ for this set of positions, where the receivers are located, and define the forward map at frequency $\omega$, $\mathcal{F}_{\omega}$ (which links the model to the data), such that,

[TABLE]

We have introduced $m(\boldsymbol{x}):=c^{-2}(\boldsymbol{x})$, which is also our choice of parameter for the reconstruction, see Remark 2. In seismic, the data are further generated from several point sources (excited one by one) and all devices (sources and receivers) remain near the surface. All devices (sources and receivers) remain near the surface ($\Gamma_{1}$), as illustrated with $\Sigma$ in Figure 1.

2.2 Quantitative reconstruction using iterative minimization

The inverse problem aims the reconstruction of the unknown medium squared slowness $m$ (i.e. the wave speed) from data $\boldsymbol{y}_{\omega}$ that connects to the forward map $\mathcal{F}_{\omega}(m_{\dagger})$ for a reference (target) model $m_{\dagger}$ with

[TABLE]

where $\mathfrak{E}_{\omega}$ represents the noise in the measurements (from the inaccuracy of the devices, model error, etc), possibly frequency dependent.

For the reconstruction, we follow the Full Waveform Inversion (FWI) method (Tarantola, 1984; Pratt et al., 1998), which relies on an iterative minimization of a misfit functional defined as the difference between the observed data and the computational simulations:

[TABLE]

where we use the standard least-squares minimization but alternatives have also been studied, as indicated in the introduction.

Remark 1** (Multi-frequency algorithm).**

For the choice of frequency in Problem 5, applications commonly use a sequence of increasing frequencies during the iterative process (Bunks et al., 1995; Pratt & Worthington, 1990; Sirgue & Pratt, 2004; Brossier et al., 2009; Barucq et al., 2018; Faucher, 2017). Namely, one starts with a low frequency and minimize the functional for the fixed frequency content in the data. Once the misfit functional stagnates, or after a prescribed number of iterations, the frequency is updated (increased) and the iterations continue, cf. Algorithm 2. Moreover, the use of sequential frequency (instead of band of frequencies) is advocated by Barucq et al. (2019a), because it enlarges the size of the basins of attraction.

Remark 2** (Parametrization of the unknown).**

For the reconstruction, we invert the squared slowness $m=c^{-2}$ instead of the velocity. The choice of this parameter is first motivated by the Helmholtz equation (2). However, it (i.e., velocity, slowness or squared slowness inversion) can lead to an important difference in the efficiency of the reconstruction procedure. It is discussed, for example, by Tarantola (1986); Brossier (2011); Köhn et al. (2012); in particular, we motivate our choice from the comparison of reconstructions provided in the context of seismic by (Faucher, 2017, Section 5.4).

Then, an iterative minimization algorithm is used for the resolution of Problem 5 in the framework of the Newton methods. Starting with an initial guess $m^{(0)}$, the model is updated at each iteration $k$, using a search direction $s^{(k)}$, such that

[TABLE]

Several possibilities exist for the search direction (e.g., Newton, Gauss-Newton, BFGS, gradient descent, etc.) and we refer to Nocedal & Wright (2006) for an extensive review of the methods. The scalar coefficient $\mu$ is approximated using a line search algorithm (Nocedal & Wright, 2006). In our implementation, we rely on a gradient-based optimization, with the nonlinear conjugate gradient method, and a backtracking line search (Nocedal & Wright, 2006). A review of the performance of first order-based minimization algorithms and the influence of line search step selection is further investigated by Barucq et al. (2018) in the context of inverse scattering. The computation of the gradient of the misfit functional is carried out using the adjoint-state method (Chavent, 1974; Plessix, 2006), which specific steps for complex-valued fields can be found in (Barucq et al., 2019b, Appendix A)

3 Regularization by discretization: model decomposition

In this section, we introduce a representation of the unknown, i.e., the model decomposition, based upon the eigenvectors of a diffusion equation. The objective is to reduce the dimension of the unknown to mitigate the ill-posedness of the inverse problem. We provide several possibilities for the choice of the eigenvectors, following the literature in image processing.

3.1 Regularization and diffusion operators

The resolution of the inverse problem using a quantitative method introduces an optimization problem (5) where the misfit functional $\mathcal{J}$ accounts for a fidelity term. The resolution of such a problem is also common in the context of image processing (e.g., for denoising or edge enhancement) where the fidelity term corresponds to the matching between the original and processed images. It is relatively common (for both the quantitative reconstruction methods and in image processing) to incorporate an additional term in the minimization, for the purpose of regularization. The primary function of this additional term is to reduce the ill-posedness of the problem, by adding a constraint. It has been the topic of several studies, we refer to, e.g., Kirsch (1996); Isakov (2006); Charbonnier et al. (1994); Robert & Deriche (1996); Rudin et al. (1992); Vogel & Oman (1996); Lobel et al. (1997); Kern (2016); Qiu et al. (2016); Kaltenbacher (2018). The regularized minimization problem writes as

[TABLE]

where $\mathcal{I}$ stands for the regularization term.

In many applications such as image processing, $\mathcal{I}$ is usually defined to only depend on the gradient of the variable (image), such that

[TABLE]

where $\phi\in L^{2}(\Omega)$. In particular, the minimum of $\mathcal{I}$ with respect to $m$ verifies the Euler–Lagrange equations (Evans, 2010; Dubrovin et al., 1992). In one dimension, it is given by (Dubrovin et al., 1992, Theorem 31.1.2) and is extended for higher dimensions with (Dubrovin et al., 1992, Theorem 37.1.2) (further simplified in our case because $\mathcal{I}$ only depends on the gradient). It states that the minimizer of $\mathcal{I}$ is the solution of the diffusion equation:

[TABLE]

For the sake of clarity, we introduce the following notation:

[TABLE]

In the following, we present several choices for the diffusion PDE coefficient $\eta$, following image processing theory.

Remark 3**.**

The minimization of $\mathcal{J}_{r}$ in Problem (7) can be performed using traditional gradient descent or Newton type algorithms. Another alternative, in particular when rewriting with the Euler–Lagrange formulation in the context of image processing, is to recast the problem as a time dependent evolution one, see, e.g., the work of Weickert (1998); Catté et al. (1992); Rudin et al. (1992); Alvarez et al. (1992).

Remark 4**.**

The diffusion equation (9) is obtained using the fact that $\phi=\phi(|\nabla m|)$ only depends on $\nabla m$. In case of dependency of the function with $m$, or higher order derivatives, the Euler–Lagrange formulation must be adapted.

3.2 Diffusion coefficients from image processing

There exist several possibilities for the choice of diffusion coefficient $\eta$ (also referred to as the weighting function) in (10), inherited from image processing theory and applications. In the following, we investigate the most common formulations, see Table 1, for which we have mainly followed the ones that are reviewed by (Blanc-Féraud et al., 1995, Table 1) and Robert & Deriche (1996). Furthermore, we incorporate a scaling coefficient $\beta>0$ for the diffusion coefficient, which impacts on the magnitude. For consistency in the different models, the norms that we employed are scaled with the maximal values so that they remain between [math] and $1$. We define,

[TABLE]

where $\mathfrak{d}$ is the space dimension ($\mathfrak{d}=2$ or $\mathfrak{d}=3$ in our experiments) and $\boldsymbol{x}$ the space coordinates: $\boldsymbol{x}=\{x_{k}\}_{k=1}^{\mathfrak{d}}$. In order to simplify the formulas, we will omit the space dependency in the following. Note that in the numerical experiments, we calculate the eigenvalues and eigenvectors from the linear differential operator $\mathcal{A}(m,\eta)$ defined in (10), where the diffusion coefficient $\eta$ is taken from the nonlinear PDE model.

Remark 5**.**

We can make the following comments regarding the nine diffusion coefficients that are introduced in Table 1.

• –

The PDE (9) using the Tikhonov diffusion coefficient $\eta_{9}$ coincides with the Laplace equation.

• –

For the formulation of $\eta_{4}$ and $\eta_{8}$, we have to impose a threshold as the coefficient is not defined for the points where the gradient is zero. In the computations, we impose that $\eta_{4}=\eta_{8}=1$ for the points $\boldsymbol{x}_{i}$ where \nabla m(\boldsymbol{x}_{i})<\text{10^{-12}}.

• –

The first Perona–Malik formula $\eta_{1}$, is very similar to the Lorentzian approach, $\eta_{6}$: only the position of $\beta$ differs. Namely, the Perona–Malik formula would rather use small $\beta$ while the Lorentzian formula would use large $\beta$.

• –

The second Perona–Malik formula, $\eta_{2}$, is very similar to the Gaussian criterion $\eta_{7}$, which only includes an additional dependency on $\beta$.

• –

The formulation of $\eta_{8}$ corresponds to the Total Variation (TV) regularization (Vogel & Oman, 1996).

3.3 Eigenvector model decomposition in FWI

In our work, we employ the regularization by discretization approach: instead of adding the regularization term $\mathcal{I}$ in the minimization problem, we remain with Problem (5), and use a specific representation for the model (unknown). We follow the work of De Buhan & Osses (2010); De Buhan & Kray (2013); Grote et al. (2017); Grote & Nahum (2019) with the “Adaptive Inversion” or “Adaptive Eigenspace Inversion” method. Namely, the unknown is represented via a decomposition into the basis of eigenvectors computed from a diffusion PDE. The purpose is to control the number of unknowns in the representation, and consequently reduce the ill-posedness of the inverse problem. The decomposition uses the steps (given in Grote et al. (2017)) depicted in Algorithm 1.

Following Algorithm 1, we introduce the notation,

[TABLE]

and

[TABLE]

Therefore, the model is represented via $N$ coefficients $\alpha$ in (14) in the basis given by the diffusion operator. The reconstruction procedure follows an iterative minimization of (5), and performs successive update of the coefficients $\alpha$. The key is that $N$ is much smaller than the dimension of the original representation of $m$, but allows an accurate resolution, as we illustrate in Sections 4 and 5. Algorithm 2 details the procedure.

Remark 6** (Minimization algorithm).**

For the minimization procedure depicted in Algorithm 2, we use a non-linear conjugate gradient method for the search direction. This method has the advantage that it only necessitates the computation of the gradient of the cost function (Nocedal & Wright, 2006). Then, to control the update step $\mu$ in Algorithm 2, a line search algorithm is typically employed (Eisenstat & Walker, 1994; Nocedal & Wright, 2006; Chavent et al., 2015; Barucq et al., 2018). This operation is complex in practice because an accurate estimation would require intensive computational operations (with an additional minimization problem to solve). Here, we employ a simple backtracking algorithm (Nocedal & Wright, 2006).

Remark 7** (Gradient computation).**

The gradient of the cost function is computed using the first order adjoint-state method (Lions & Mitter, 1971; Chavent, 1974), which is standard in seismic application (Plessix, 2006). It avoids the formation of a dense Jacobian matrix and instead requires the resolution of an additional PDE, which is the adjoint of the forward PDE, with right-hand sides defined from the difference between the measurements and the simulations, see Plessix (2006); Faucher (2017); Barucq et al. (2018, 2019b) for more details.

In our implementation, the gradient is first computed with respect to the original (nodal) representation and we use the chain rule to retrieve the gradient with respect to the decomposition coefficients $\alpha$:

[TABLE]

It is straightforward, from (14), that the derivation for a chosen coefficient ${\alpha_{l}}$ gives $\partial_{\alpha_{l}}m=\psi_{l}$. Therefore, it is computationally easy to introduce the formulation with respect to the eigenvector decomposition from an already existing ‘classical’ (i.e. when the derivative with respect to the model is performed) formulation: it only necessitates one additional step with the eigenvectors.

Remark 8**.**

In Algorithm 2, the basis of eigenvectors remains the same for the complete set of iterations, and is extracted from the initial model. Only the number of vectors taken for the representation, $N_{i}$, changes. Namely, from $N_{1}$ to $N_{2}>N_{1}$, the decomposition using $N_{2}$ still has the same $N_{1}$ first eigenvectors in its representation (with different weights $\alpha$), and additional $(N_{2}-N_{1})$ eigenvectors. As an alternative, we investigate the performance of an algorithm where the basis changes at each frequency (i.e. it is recomputed from the current iteration model), see Appendix A.

3.4 Numerical implementation

Our code is developed in Fortran90, it uses both mpi and OpenMP parallelism and run on cluster111The experiments have been performed on the cluster PlaFRIM (Plateforme Fédérative pour la Recherche en Informatique et Mathématiques, https://www.plafrim.fr/fr) with the following node specification: 2 Dodeca-core Haswell Intel Xeon E5–2680 v3 (2.5GHz); 128Go RAM; Infiniband QDR TrueScale: 40Gb s*-1*, Omnipath 100Gb s*-1*. for efficiency. The forward wave operator is discretized using a Finite Differences scheme, e.g. Virieux (1984); Operto et al. (2009); Wang et al. (2011). The discretization of the Helmholtz operator generates a large sparse matrix, for which we use the direct solver Mumps (Amestoy et al. (2001, 2006)) for its factorization and the resolution of linear system. This solver is particularly optimized and designed for this type of linear algebra problems, i.e. large, sparse matrices. Our preference for direct solver instead of iterative ones is mainly motivated by two reasons:

seismic acquisition is composed of a large amount of sources, i.e. a large amount of right-hand sides (rhs) to be processed for the linear system. Using direct solver, the resolution time is very low once the factorization is performed, hence it is well adapted for the multi-rhs seismic configuration. 2. 2.

For the minimization algorithm, the gradient is computed via the adjoint-state method (see Remark 7). It means that an additional linear system has to be solved, which is actually the adjoint of the forward one. Here, the factors obtained from the factorization of the forward operator can be directly reused, and allow a reduced computational cost, see Barucq et al. (2018) and the references therein.

The next step is the computation of the eigenvectors associated with the smallest eigenvalues for the diffusion operator. We use the package Arpack222www.caam.rice.edu/software/ARPACK/, Arpack uses sequential computation, hence, contrary to the rest of our code, this part does not use parallelism. Future developments include the implementation of the parallel version of the package: Parpack., which is devoted to solve large sparse eigenvalue problems using iterative methods. More precisely, it uses implicitly restarted Lanczos or Arnoldi methods, respectively for symmetric and non-symmetric matrices, Lehoucq & Sorensen (1996). Several options are available in the package, including the maximum number of iterations allowed, or a tolerance parameter for the accuracy of acceptable solution333 We have observed important reduction of time cost when allowing some flexibility in the accuracy with this threshold criterion. However, in the computational experiments, we do not use this option, as the numerical efficiency is not the primary objective of our study. .

Remark 9** (Eigenvectors associated with the lowest eigenvalues).**

The Lanczos and Arnoldi methods are particularly efficient to compute the largest eigenvalues and associated eigenvectors of the matrix, and only require matrix vector multiplication. However, we are interested in the lowest eigenvalues for our decomposition. The idea is simply to use that the lowest eigenvalues of the discretized diffusion matrix, say $A$, are simply the largest eigenvalues of the matrix $A^{-1}$. Then, the matrix-vector multiplication, say $Av$ for a vector $v$, becomes a resolution of a linear system $A^{-1}v$. It may appear computationally expensive but it is not thanks to the use of the direct solver Mumps (see above), which, once the factorization is obtained, is very efficient for the resolution procedure. Hence, the computation follows the steps444Arpack has the possibility to compute the smallest eigenvalues using matrix-vector multiplication, however, we have observed a drastic increase of the computational time compared to using the inverse matrix and resolution of linear system.:

compute the (sparse) matrix discretization of the selected diffusion operator: $A$; 2. 2.

compute the factorization of the matrix $A$ using Mumps, 3. 3.

use the package Arpack to compute the largest eigenvalues of $A^{-1}$, by replacing the matrix-vector multiplication step in the iterations by the resolution of a linear system using Mumps.

Finally, the last step is to retrieve the appropriate coefficients $\alpha_{k}$ in (14) for the decomposition. It basically consists in the resolution of a dense linear system (from least squares method). We use Lapack, Anderson et al. (1999) (contrary to Mumps, Lapack is adapted to dense linear system). Note that, because we usually consider a few hundreds of coefficients for the decomposition, this operation remains relatively cheap compared to the eigenvectors computation. We compare the computational time for the eigenvectors computation and model decomposition in Figure 2 for different values of $N$ and $\eta$. We first note that the choice of $\eta$ does not really modify the computational time. Then, we see that the two operations are mostly linear in $N$, and that the time to solve the least squares problem with Lapack is much smaller than the time to compute the eigenvectors with Arpack, namely, hundred time smaller. For the largest case: $N=1000$ for a squared matrix of size 277221, it takes about 30min to retrieve the eigenvectors, and $10$ s to compute the $\alpha$ (in our applications we usually take $N<100$).

4 Illustration of model decomposition

First, we illustrate the eigenvector model decomposition with geophysical media in two dimensions. The original model is represented on a structured grid by $n_{x}\times n_{z}$ coefficients, and we have here $921\times 301=\text{277221}$ coefficients. We consider three media of different nature:

• –

the Marmousi velocity model, which consists in structures and faults, see Figure 3(a);

• –

a model encompassing salt domes: objects of high contrast velocity, see Figure 3(b);

• –

eventually, the SEAM Phase I velocity model which consists in both salt and layer structures, see Figure 10.

All three models uses the same number of coefficients for their representations, and the first two are actually of the same size ($9.2\times 3$ km).

We perform the decomposition of the models by application of Algorithm 1, and steps (12), (13) and (14) (and we recall that we use the linear PDE problem). We study the main parameters of the decomposition:

• –

the choice of $\eta$, with the possibilities given in Table 1,

• –

the choice of the scaling parameter $\beta$ in the formulation of $\eta$ (Table 1),

• –

the number of eigenvectors $N$ employed for the decomposition in (14).

The accuracy of the decomposition is estimated using the $L^{2}$ norm of the relative difference between the decomposition and the original representation such that

[TABLE]

where $m$ is the original model (Figure 3) and $\mathfrak{m}$ the decomposition using the basis of eigenvectors from Algorithm 1.

4.1 Decomposition of noise-free models

We decompose the salt and Marmousi models using the nine possibilities for $\eta$, that are given in Table 1. For the choice of scaling coefficient $\beta$ (which does not affect $\eta_{8}$ and $\eta_{9}$), we roughly cover an interval from $10^{-7}$ up to $10^{6}$, namely: \{$$10^{-7}, $10^{-6}$, $10^{-5}$, $10^{-4}$, $10^{-3}$, $10^{-2}$, 5$.10^{-2}$, $10^{-1}$, 5$.10^{-1}$, 1, 5, 10, $10^{2}$, $10^{3}$, $10^{4}$, $10^{5}$, 10^{6}$$\}. In Tables 2 and 3, we show the best relative error (i.e. minimal value) obtained for the Marmousi model of Figure 3(a) and the salt model of Figure 3(b). We test all choices of $\eta$ and values of $N$ between $10$ (coarse) and $500$ (refined). The corresponding values of the scaling parameter $\beta$ which gives the best (i.e. the minimal) error are also given in parenthesis.

As expected, we observe that the more eigenvectors are chosen (higher $N$), the better will be the decomposition. When using $500$ eigenvectors, which represents about 2% of the original number of coefficients ($921\times 301$ in Figure 3), the error is of a few percent only. This can be explained by the redundancy of information provided by the original fine grid where the model is represented (e.g. the upper part of Figure 3(b) and the three salt bodies are basically constant). Comparing the methods and models, we see that

• –

the Marmousi model (Table 2) is harder to retrieve than the salt model (Table 3) as it gives higher errors. In particular for low $N$, the salt model can be acutely decomposed (possibly $3$% error with $N=10$).

• –

For both models, it appears that four methods stand out: $\eta_{1}$ (Perona–Malik), $\eta_{3}$ (Geman & Reynolds (1992)), $\eta_{5}$ (Charbonnier et al. (1994)) and $\eta_{6}$ (Lorentzian), with a slight advantage towards $\eta_{1}$.

• –

The scaling coefficient that minimizes the error is consistent with respect to $N$, with similar amplitude. However, changing the model may require the modification of $\beta$: between the salt and Marmousi decomposition, the optimal $\beta$ is quite different for $\eta_{1}$, and also for $\eta_{6}$.

To investigate further the last point, we show the evolution of relative error $\mathcal{E}$ with respect to the scaling coefficient $\beta$ for the decomposition of the Marmousi and salt models in Figure 4, where we compare four selected formulations for $\eta$. We observe some flexibility in the choice of $\beta$ that gives an accurate decomposition for $\eta_{1}$ and $\eta_{5}$. On the other hand, $\eta_{3}$ and $\eta_{6}$ show sharp functions, which means that the selection of $\beta$ has more influence in these cases (and must be carefully taken). In addition, the range of efficient $\beta$ changes depending on the model decomposed, except for $\eta_{5}$. It demonstrates that the choice of $\beta$ for optimality is not trivial in general, and is model-dependent.

We then picture the resulting images obtained after the decomposition of both models, see Figures 5 and 6. The pictures illustrate correctly the observations of the tables and the differences between the formulation. The salt model is usually well recovered with all formulations, while the Marmousi model is more hardly discovered, except with $\eta_{1}$, $\eta_{3}$ and $\eta_{6}$. Those three formulations are the only ones able to capture the structures.

4.2 Decomposition of noisy models (denoising)

We incorporate noise in the representation for getting closer to the reality of applications where few information on the model is available. Hence, we reproduce the model decomposition, this time working with noisy pictures. For every nodal velocity (of Figure 3), we recast the values using an uniform distribution that covers $\pm 20\%$ of the noiseless value. The resulting media are illustrated in Figure 7.

We apply the model decomposition using the different formulations of $\eta$ and choice of scaling coefficient $\beta$, following the procedure employed for the noiseless model. In Tables 4 and 5, we show the evolution of best relative error with $N$, for the noisy Marmousi and salt models respectively. Here, the relative error is computed from the difference between the noiseless model and the decomposition of the noisy one. The objective of the regularization is to preserve the structures while smoothing out the noise effect.

The decomposition of noisy pictures requires more eignevectors for an accurate representation. Then, the salt model, with high contrast objects, still behaves better than the many structures of the Marmousi model. For the decomposition of the noisy Marmousi model, none of the formulations really stands out and the error never reaches below 10% using at most $N=250$.

In Figures 8 and 9, we picture the resulting decomposition for the two media. For the decomposition of the Marmousi model, we use $N=200$; and $N=50$ for the salt model. It corresponds to higher values compared to the pictures shown for the noiseless models (Figures 5 and 6).

The decomposition of the salt model remains acceptable, and we easily distinguish the main contrasting object. The smaller objects also appear, in a smooth representation. The formulations using $\eta_{1}$, $\eta_{3}$ and $\eta_{6}$ provide sharper boundary for the contrasting objects, in particular for the upper interface. Regarding the decomposition of the noisy Marmousi model, it illustrates the limitation of the method, where none of the formulations is really able to reproduce the structures, and most edges are lost. In particular, the central part of the model is mostly missing and the amplitude of the values has been reduced. It seems that $\eta_{1}$, $\eta_{3}$ and $\eta_{6}$ are slightly more robust and gives (relatively speaking) the best results. To conclude, these three formulations appear less sensitive (for both media) to noise than the other ones.

In the context of image decomposition, we have can draw the following conclusions for the the decomposition using the basis of eigenvectors.

• –

The method is efficient to represent media with contrasting shapes (e.g., salt domes), even when noise is contained in the images. In this case, the choice of $\eta$ does not really affect the representation of the objects, and all methods behave quite well, see Figure 9.

• –

The performance of the decomposition strongly depends on the media, and diminishes with thin structures as in the Marmousi model. In this case, an appropriate choice of formulation ($\eta_{1}$, $\eta_{3}$, $\eta_{6}$ from Table 2) can provide the accurate representation for noise-free picture but when incorporating noise, the performance deteriorates and the edge contrasts are lost.

4.3 Sub-surface salt with layers: SEAM Phase I model

We have seen that the decomposition behaves well when a contrasting object with sharp contrast belongs to the medium, while structures/layers are hardly represented. We pursue our investigation with a common geophysical configuration where both salt-domes and layers exist in the subsurface. We use a velocity model extracted from the SEAM (SEG Advanced Modeling Program ) Phase I benchmark555see https://wiki.seg.org/wiki/Open_data. a consider a medium of size $17.5\times 3.75$ km. Per consistency with the previous experiment, it is represented using a grid of $921\times 301$ points, see Figure 10.

Similarly to our previous experiment, we investigate the noiseless model and a noisy one, which incorporates $\pm 20$ % error. The relative error and corresponding scaling coefficients for both models are given in Tables 6 and 7. The relative error is of a few percent for high $N$, and we observe important differences between the formulation. Here, $\eta_{8}$, $\eta_{9}$ and $\eta_{4}$ give the worst results.

We further illustrate the decomposition in Figure 11, using $\eta_{1}$ only for the sake of clarity.

This model, which encompasses salt and layers, is well recovered with a decomposition using $N=20$ when there is no noise. In case of noise, it needs higher $N$ and the contrasting shapes are smoothed out. Nonetheless, compared to the Marmousi model, the decomposition is able to capture the main features.

5 Experiments of reconstruction with FWI

In this section, we perform seismic imaging following the FWI Algorithm 2 for the identification of the subsurface physical parameters. We focus on media encompassing salt domes, as we have shown that the decomposition is more appropriate, and it also gives a more challenging configuration in seismic applications. We follow a seismic context, where the forward problem is given by (2) and the data are restricted to be acquired near the surface (see Figure 1). The challenges of our experiments is threefold:

the recovery of salt domes is recognized to be difficult (Barucq et al., 2019a); 2. 2.

we consider an initial guess that has no information on the subsurface structures and where the background velocity amplitude is incorrect. 3. 3.

we avoid the use of the (unrealistically) low frequencies (below $2$ Hz in exploration seismic).

5.1 Reconstruction of two-dimensional salt model

We first consider a two-dimensional salt model of size $9.2\times 3$ km, which consists in three domes, see Figure 13(a). We generate the data using $91$ sources and $183$ receivers (i.e. data points) per source. Both devices are located near the surface: the sources are positioned on a line at $20$ m depth and the receivers at $80$ m depth. In order to establish a realistic situation despite having a synthetic experiment, the data are generated in the time-domain and we incorporate white noise in the measurements. The level of the signal to noise ratio in the seismic trace is of $10$ dB, the noise is generated independently for all receivers record associated with every source. Then we proceed to the discrete Fourier transform to obtain the signals to be used in the reconstruction algorithm. In Figure 12, we show the time-domain data with noise and the corresponding frequency data for one source, located at $20$ m depth, in the middle of the $x$-axis.

For the reconstruction of the salt dome model, the starting and true model are given in Figure 13. We do not assume any a priori knowledge of the contrasting object in the subsurface, and start with a one-dimensional variation, which has a drastically lower amplitude (i.e., the background velocity is unknown).

5.1.1 Fixed decomposition, single frequency reconstruction

We first only use $2$ Hz frequency data, and perform $180$ iterations for the minimization. In Figure 14, we show the reconstruction where the decomposition has not been employed, i.e. the model representation follows the original piecewise constant decomposition of the model (one value per node). In Figure 15, we compare the reconstruction using Algorithm 2 for the different formulations of $\eta$ given in Table 1, using a fixed $N=50$. We use $N=100$ for Figure 16. For the sake of clarity, we focus on $\eta_{3}$ (the most effective formulation), and $\eta_{8}$ (which relates to the Total Variation regularization) and move the complete pictures with comparison of all formulations in Appendix B, Figure B1 for $N=100$.

We observe that

• –

the traditional FWI algorithm (without decomposition), see Figure 14, fails to recover any dome. It only shows some thin layers of increasing velocity, with amplitudes much lower than the original ones.

• –

The decomposition using $N=50$ is able to discover the largest object with formulation $\eta_{1}$, $\eta_{3}$, $\eta_{5}$, $\eta_{6}$ and $\eta_{8}$, see Figure 15. The best result is given by $\eta_{3}$ which recovers the three domes; while $\eta_{1}$, $\eta_{5}$ and $\eta_{8}$ show artifacts in the lower right corner. The other decompositions fail. We note that, due to the lack of velocity background information, the positions of the domes are slightly above the correct ones to compensate for the low travel times.

• –

In this experiment, the method behaves much better with a restrictive number of eigenvectors. With $N=100$ (Figures 16 and B1), the iterative reconstruction only results in artifacts. The restrictive number of eigenvectors provides a regularization of the problem by reducing the number of parameters, which is crucial. For instance, the stability is known to deteriorate exponentially with the number of parameters in the piecewise constant case (Beretta et al., 2016).

Opposite to the decomposition of images (Section 4), the quantitative reconstruction using a model represented with a basis of eigenvectors from a diffusion PDE shows drastic differences between the formulations, where the procedure can fail depending on the choice of $\eta$. In addition, the number of eigenvectors for the representation has to be carefully selected, see Subsection 5.3.

5.1.2 Experiments with increasing $N$ and multiple frequencies

We investigate the performance of the eigenvector decomposition for multiple frequency data, and with progressive evolution of the number of eigenvectors in the representation $N$. We have a total of four different experiments, which are summarized in Table 8. The reconstructions, for $\eta_{3}$ and $\eta_{8}$, are shown Figure 17. The results for all $\eta$ of Table 1 are pictured in Appendix B, Figures B2, B3 and B4.

From these experiments using multiple $N$ and/or frequency contents, we observe the performance of the method. The best results are obtained using a single $2$ Hz frequency with either progression of increasing $N$ or constant $N$: Experiments 2 and 1. The progression of $N$, Experiment 2, appears the most robust.

In this experiment, using multiple frequencies does not improve the results. It is probably due to the lack of knowledge of the velocity background which prevents us from recovering finer scale (i.e., kinematic error, Bunks et al. (1995)). In particular, local minima in the misfit functional become more and more prominent in the high-frequency regime (Bunks et al., 1995; Barucq et al., 2019b). We illustrate the performance of the reconstruction with the results in the data-space: in Figure 18, we show the time-domain seismograms for the true, starting and reconstructed velocity models. We observe that when we filter out frequencies above 2 Hz (first line of Figure 18), the trace from the reconstructed model is indeed very similar to the measured one. However, when encompassing all frequency contents (bottom line of Figure 18), important differences arise, in particular, one can see the travel time of the first reflection which is earlier with the recovered model. This indicates that the location of the salt in the reconstructed velocity is above its ‘true’ position.

Then, while we incorporate the higher frequency in the minimization procedure, the FWI is not amenable to improve the results (see Figure 17) and it is most likely due to the missing velocity background which is not improved during the first iterations, and still missing. In Figure 19, we show the frequency-domain data at $2$ and $4$ Hz: the observed data at $2$ Hz are accurately obtained with the model reconstructed with the decomposition in eigenvectors, which confirms the pertinence of the method. Interestingly, at $4$ Hz, while the frequency is not even used in the inversion scheme (we only use $2$ Hz for Figure 17(a)), we already have a good correspondence near the source and only the parts further away show a shift. In order to overcome the issue of recovering the background velocity, one would need lower frequency content, or one could employ alternative strategies, such as the MBTT method, based upon the decomposition of the background velocity model and the reflectivity (Clément et al., 2001). Here, the decomposition in eigenvectors appropriately recover the reflectivity part (better than traditional FWI), but the background model remains missing.

5.1.3 On the choice of the number of eigenvectors

We have shown in Figures 15 and 16 that one should take an initial $N$ relatively low for the reconstruction algorithm to succeed. It remains to verify if the appropriate $N$ can be selected ‘a priori’, or based upon minimal experiments. In Figure 20(a), we show the evolution of the misfit functional with thirty iterations, for different values of $N$, from $10$ to $250$. We compare, in Figure 20(b), with the progression of $N$, which follows Experiment 2 of Table 8.

From Figure 20(a), we see that all of the choices of $N$ have the same pattern: first the decrease of the functional and then its stagnation. We notice that the good choice for $N$ is not reflected by the misfit functional. Indeed, it shows lower error for larger $N$, while they are shown to result in erroneous reconstructions (Figure 16 compared to Figure 15). It is most likely that using larger $N$ leads to local minima and/or deteriorates the stability (see, for the piecewise constant case, Beretta et al. (2016)). It results in the false impression (from the misfit functional point of view) that it would improve the efficiency of the method. Using a progression of $N$, Figure 20(b), eventually gives the same misfit functional value than the large $N$, but it needs more iterations. This increase of iterations and ‘slow’ convergence is actually required, because it leads to an appropriate reconstruction, see Figure 17.

Therefore, we cannot anticipate a good choice for $N$ a priori (with a few evaluation of the misfit functional). the guideline we propose, as a safe and robust approach, is the progression of increasing $N$, from low to high: it costs more in terms of iterations, but it converges properly.

5.2 Reconstruction of the SEAM Phase I model

We now consider the recovery of the SEAM Phase I model, which is expected to be more challenging as it contains both a salt-dome and sub-salt layers. The starting model for the reconstruction is shown in Figure 21, where, for the sake of clarity, we also picture the reference model which was previously used for the decomposition. This medium is of size $17.5\times 3.75$ km and the starting model we use is a smooth version of the reference one, where the contrasting objects and layers are missing.

We follow the same configuration as in the previous experiment: the data are generated in the time-domain, and noise is incorporated to the seismograms using a signal-to-noise ratio of 10 dB. Next, we proceed with the discrete Fourier-transform to employ the iterative procedure with time-harmonic waves. In this experiment, the smallest available frequency is $2$ Hz.

For the sake of conciseness, we only present the reconstruction results with the the representation using $\eta_{3}$ (which was the most efficient). We follow a slow increase of $N$ in a fixed basis (analogous to Experiment 2 of Table 8) which was the more stable approach. The reconstruction using $2$ Hz is shown in Figure 22.

We observe that the standard FWI algorithm gives artifacts over the medium, with oscillatory patterns, in particular on the sides. On the other hand, the reconstruction using the representation based upon the eigenvector decomposition is stable, and is able to capture accurately the upper boundary of the salt dome. Because of the frequency employed, only the long wavelengths are recovered at this stage. In Figure 23, we further illustrate the recovery by showing the frequency-domain data at $2$ Hz frequency using the different wave speed models. We can see that the data from the starting model are out of phase as soon as we move away from the source, with cycle-skipping effects. However, the reconstruction using the eigenvector decomposition is able to retrieve this information and accurately capture the oscillations of the signal, and only the amplitude is inaccurate.

We now continue the procedure using increasing frequencies. Because of the smoothing effect of the decomposition, we employ the algorithm without the eigenvector representation (alternatively, one could use the decomposition but with large $N$). Therefore, the velocity model obtained from the FWI with eigenvector decomposition of Figure 22(b) is used as an initial model for restarted FWI with multiple frequencies, from $2$ to $10$ Hz (we use the sequential progression advocated by Barucq et al. (2019a)). Eventually, the reconstruction after $10$ Hz iterations is pictured in Figure 24.

The reconstruction is able to capture the finer details of the velocity model: the salt dome is clearly defined, and the sub-salt layer starts to appear. Therefore, the eigenvector decomposition method can also serve to build initial models for the FWI algorithm, where it appears as an interesting alternative to overcome the lack of low frequency. To illustrate the different steps of the reconstruction, we show the time-domain seismograms at the different stages of the reconstruction in Figure 25 (while the $2$ Hz frequency-domain data are given in Figure 23). We compare the traces resulting from the initial and true wave speed models, from the partial reconstruction obtained with the eigenvector decomposition (Figure 22(b)), and from the final reconstruction (Figure 24).

The trace that uses the starting model mainly shows the first arrivals, with some minor reflections coming from the smoothing of the salt dome. The one using the partial reconstruction with eigenvenctor decomposition and only the $2$ Hz data accurately resolves the main multiple reflections between the salt upper boundary and the surface (the thick line in the center of the red ellipses in Figure 25) but it misses the events of smaller importance. Eventually, the final reconstruction, that uses up to $10$ Hz data contents is able to reproduce some of the events of smaller amplitudes. We also note that the amplitude of the trace for the final reconstruction is slightly high, which indicates that the contrasts are even too strong.

5.3 Implementation of the method

The eigenvector decomposition for the model representation depends on different choices of parameters, and it is not trivial to efficiently implement the method in the iterative reconstruction procedure. From the experiments we have carried out, we propose the following strategy.

1. (1).

Regarding the choice of weighting coefficient, $\eta_{3}$ from Geman & Reynolds (1992) is the most efficient in these applications, and supersedes the standard Total Variation approach (i.e., $\eta_{8}$).

Secondly, the difficulty resides in the number of eigenvectors to take for the decomposition: $N$. More important, it appears that the misfit functional does not provide us with a good indication (see Figure 20).

1. (2).

The number of eigenvector $N$ for the decomposition should initially takes a low value, and progressively evolves to higher values with iterations. It may not be the fastest convergence, but it is the most robust approach.

Finally, the reconstruction can serve as an initial model for multi-frequency data:

1. (3).

the (partial) reconstruction with eigenvector decomposition is used as a starting model for multi-frequency algorithm. It allows the recovery of the finer details, which depend on the smaller wavelengths and where the smoothing effect is misleading.

Namely, the decomposition is particularly efficient to overcome the lack of low-frequency in the data.

5.4 Three-dimensional experiment

The method extends readily for three-dimensional model reconstruction, simply incurring a larger computational cost (as larger matrices are involved for the eigenvector decomposition and the forward problem discretization). We proceed with a three-dimensional experiment, where we consider a subsurface medium of size $2.46\times 1.56\times 1.2$ km, encompassing several salt domes, illustrated in Figure 26. The seismic acquisition consists in 96 sources, positioned on a two-dimensional plane at 10 m depth; one thousand receivers are positioned at 100 m depth. Similar to the previous experiments, the data are first generated in the time-domain and we incorporate noise before we proceed to the Fourier transform. Figure 28 shows the time-domain data associated with a centrally located source, and the corresponding Fourier transform at $5$ Hz frequency. For the reconstruction, we start with a one-dimensional variation, in depth only, where none of the objects is intuited, see Figure 27, and the velocity background is incorrect.

In Figure 29, we show the reconstruction without employing the eigenvector decomposition, where the wave speed has a piecewise constant representation on a $124\times 79\times 61$ nodal grid. We only use $5$ Hz frequency data, and $30$ iterations. Next, we employ the eigenvector model representation with Algorithm 2. Following the discussion in Subsection 5.3, we select $\eta_{3}$, which is the most robust, and try two situations:

• –

single frequency ($5$ Hz), fixed $N$ reconstruction using $\eta_{3}$, $N=50$ and $30$ iterations, the final reconstruction is Figure 30;

• –

single frequency ($5$ Hz), multiple $N$ reconstruction using $\eta_{3}$, $N=\{20,30,50,75,100\}$ and $30$ iterations per $N$, i.e. $180$ iterations in total: the final reconstruction is Figure 31.

For visualization, we focus on the two-dimensional vertical and horizontal sections which illustrate the positions and shapes of the objects. This experiment is consistent with the two-dimensional results, and we observe the following.

• –

The classical FWI reconstruction fails to discover the subsurface objects, with only a narrow layer, misplaced, see Figure 29.

• –

The reconstruction with the representation with eigenvector decomposition is able to accurately capture the largest subsurface salt domes, see Figures 30 and 31.

• –

The best results are obtained when we use the progression of $N$, with a single frequency, see Figure 31. The main salt dome is captured and smaller ones start to appear, including near the boundary. Using a single $N$ is also a good candidate, Figure 30, as it necessitates much less iterations ($30$ instead of $180$) hence less computational time.

6 Conclusion

In this paper, we have investigated the use of a representation based upon a basis of eigenvectors for image decomposition and quantitative reconstruction in the seismic inverse problem, using two and three-dimensional experiments. We have implemented several diffusion coefficients, and compared their performance, depending on the target medium.

1. (1).

In the context of image decomposition, the case of contrasting objects (salt domes) is clearly more appropriate than the layered media (such as Marmousi). All of the diffusion coefficients behave well and provide satisfactory results for salt domes image decomposition, even in the presence of noise. 2. (2).

For the decomposition of images with layered patterns, only a few formulations perform well ($\eta_{1}$, $\eta_{3}$, $\eta_{6}$). It would be interesting to investigate further the performance of anisotropic or directional diffusion coefficients, mentioned by Weickert (1998); Grote & Nahum (2019).

Next, we have considered the quantitative reconstruction procedure in seismic, where only partial, backscattered (i.e. reflection) data are available, from one side illumination. We have probed the performance of the method by considering initial guesses with minimal a priori information, and by avoiding the low frequency data, which are not accessible in seismic applications. In this context, the FWI algorithm based upon the eigenvector model representation has shown promising results of subsurface $3$D salt dome media. The method only requires the preliminary computation of the basis of eigenvectors associated with the diffusion operator, and a trivial modification of the gradient computation. Namely, the overall additional cost of the method remains marginal compared to the cost of FWI. Our findings are the following.

1. (3).

For reconstruction, the result depends on the choice of diffusion coefficient. We recommend $\eta_{3}$, from Geman & Reynolds (1992), which was the most robust in our applications, even with a fixed $N$ and a few iterations. 2. (4).

We have shown that the choice of $N$ is not trivial, and one cannot rely on the misfit functional evaluation. Therefore, we have proposed a progression of increasing $N$, which appears to stabilize the reconstruction. 3. (5).

Because the method has a smoothing effect, it focuses on the long wavelength structures. Thus, the reconstruction using the decomposition can serve as an initial model to iterate with higher frequency contents, in order to improve the resolution.

Following these analyses, some difficulties remain regarding the optimal choice of parameters. For instance, it is possible that $\eta$ has to be selected differently depending on the model (as illustrated with the Marmousi decomposition). Similarly, the scaling coefficient $\beta$ would affect the performance but the acceptable range appears, hopefully, quite large.

For future work, it seems that the lack of background velocity information would not be overcome by the decomposition, possibly resulting in artifacts. Therefore, we envision the use of multiple basis to parametrize the velocity (e.g., using the background/reflectivity decomposition idea of Clément et al. (2001); Barucq et al. (2019b), with a dedicated smooth eigenvector basis to represent the background, and another to represent the reflectors).

Eventually, the method can readily extend to multi-parameter inversion (e.g. for elastic medium), upon taking a separate basis per parameter (e.g. one for each of the Lamé parameters in linear elasticity). However, a more appropriate approach would be the use of joint-basis by, e.g., considering a system of PDE instead of the scalar diffusion operator. We have in mind strategies such that joint-sparsity and tensor decomposition.

Acknowledgments

The authors would like to thank Prof. Marcus Grote for thoughtful comments and discussions. The research of FF is supported by the Inria–TOTAL strategic action DIP. FF acknowledges funding from the Austrian Science Fund (FWF) under the Lise Meitner fellowship M 2791-N. OS is supported by the FWF, with SFB F68, project F6807-N36 (Tomography with Uncertainties). OS also acknowledges support from the FWF via the project I3661-N27 (Novel Error Measures and Source Conditions of Regularization Methods for Inverse Problems). HB acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement Number 777778 (Rise action Mathrocks) and with the E2S–UPPA CHICkPEA project.

Appendix A Two-dimensional reconstruction with evolution of basis

In this appendix, we experiment the re-computation of the eigenvectors basis along with the iterations. We consider two additional experiments, which are derived from Experiments 1 and 2 of Table 8, with one major difference: the set of eigenvectors is recomputed from the current iteration model every $30$ iterations. This is advocated in De Buhan & Kray (2013); Grote & Nahum (2019) where, contrary to our experiments, the background velocity is mostly known.

Therefore, for the $180$ total iterations, the basis are

• –

computed from the initial model when starting the very first iteration,

• –

re-computed from the current reconstruction at the beginning of iterations $31$, $61$, $91$, $121$, $151$.

Compared to Algorithm 2, instead of using a fixed ${\boldsymbol{\psi}}_{\text{loc}}$ from the initial model, it is recomputed from the current ${\mathfrak{m}}^{(k)}$. The reconstructions using the update of basis are shown Figures A1 and A2, for $N=50$ and $N=\{50,60,70,80,90,100\}$ respectively.

We observe that changing the basis along with the iterations provides similar accuracy of the reconstruction compared to keeping the basis fixed. For some formulations, the salt dome appears slightly larger than with the fixed basis but the artifacts at the bottom are also stronger, with patches of high velocities.

We believe that the difficulty is due to the lack of prior information regarding the velocity background. Namely, the starting model has a background profile of low amplitude compared to the target model, see Figure 13(b). While our reconstruction captures the contrasting object, the background remains mostly erroneous (i.e. $\mathfrak{m}_{0}$ in (14)), impacting the decomposition. Then, when we update the decomposition, the new set of eigenvectors will try to encompass the background variation, creating the bottom artifacts.

Appendix B Extended set of figures

We provide here the reconstructions obtained for all $\eta$ for the two-dimensional salt dome FWI test case of Subsection 5.1.