Full Reciprocity-Gap Waveform Inversion in the frequency domain, enabling sparse-source acquisition

TL;DR

This paper introduces Full Reciprocity-gap Waveform Inversion (FRgWI), a new seismic imaging method that improves subsurface property reconstruction by allowing flexible, sparse-source data acquisition and demonstrating superior robustness and accuracy over traditional methods.

Contribution

The paper presents FRgWI, a reciprocity-based waveform inversion method that enables flexible source and measurement configurations, reducing computational costs and improving robustness in seismic imaging.

Findings

FRgWI outperforms traditional FWI in 3D experiments.

It allows for sparse computational or observational acquisitions.

FRgWI is more robust against cross-talk than least-squares methods.

Abstract

The quantitative reconstruction of sub-surface Earth properties from the propagation of waves follows an iterative minimization of a misfit functional. In marine seismic exploration, the observed data usually consist of measurements of the pressure field but dual-sensor devices also provide the normal velocity. Consequently, a reciprocity-based misfit functional is specifically designed, and defines the Full Reciprocity-gap Waveform Inversion (FRgWI ) method. This misfit functional provides additional features compared to the more traditional least-squares approaches with, in particular, that the observational and computational acquisitions can be different. Therefore, the positions and wavelets of the sources from which the measurements are acquired are not needed in the reconstruction procedure and, in fact, the numerical acquisition (for the simulations) can be arbitrarily chosen. Based on three-dimensional experiments, FRgWI is shown to behave better than Full Waveform Inversion (FWI) in the same context. Then, it allows for arbitrary numerical acquisitions in two ways: when few measurements are given, a dense numerical acquisition (compared to the observational one) can be used to compensate. On the other hand, with a dense observational acquisition, a sparse computational one is shown to be sufficient, for instance with multiple-point sources, hence reducing the numerical cost. FRgWI displays accurate reconstructions in both situations and appears more robust with respect to cross-talk than the least-squares shot-stacking.

Equations

$$-\mathrm{i}\omega \rho (\mathbf{x})\mathbf{v}(\mathbf{x})$$

$$-\mathrm{i}\omega \kappa (\mathbf{x}{)}^{-1}p(\mathbf{x})$$

$$p(\mathbf{x})$$

$${\partial }_{\mathbf{\nu }}p(\mathbf{x})-\frac{\mathrm{i}\omega }{c(\mathbf{x})}p(\mathbf{x})$$

$$c(\mathbf{x})=\sqrt{\frac{\kappa (\mathbf{x})}{\rho (\mathbf{x})}}.$$

$$\mathcal{F}_{\omega}^{(f)}(m)=\Big{\{}p^{(f)}_{\omega}\mid_{\Sigma}\,,\,\,v_{\boldsymbol{\nu},\omega}^{(f)}\mid_{\Sigma}\Big{\}}.$$

$$p\mid_{\Sigma}=\big{\{}\,p(\boldsymbol{x}_{1}),\,\ldots,\,p(\boldsymbol{x}_{n_{\text{rcv}}})\,\big{\}},$$

$$\mathcal{R}\big{(}p\big{)}=p\mid_{\Sigma}~{},\quad\quad\mathcal{R}\big{(}v_{\boldsymbol{\nu}}\big{)}=v_{\boldsymbol{\nu}}\mid_{\Sigma}.$$

$$d_{\omega}^{(f)}=\Big{\{}\,d_{\omega,p}^{(f)}\,,\,\,\,d_{\omega,v}^{(f)}\,\Big{\}},$$

$$\mathcal{J}_{L2}(m)=\dfrac{1}{2}\sum_{i=1}^{n_{\text{src}}}\Big{\|}\mathcal{R}\big{(}p^{(f_{i})}\big{)}-d_{p}^{(f_{i})}\Big{\|}^{2}_{2}+\dfrac{\eta}{2}\Big{\|}\mathcal{R}\big{(}v_{\boldsymbol{\nu}}^{(f_{i})}\big{)}-d_{v}^{(f_{i})}\Big{\|}^{2}_{2},$$

$${\mathcal{J}}_{\mathfrak{r}}(m)$$

$$\displaystyle=\dfrac{1}{2}\sum_{i=1}^{n_{\text{src}}^{\text{obs}}}\sum_{j=1}^{n_{\text{src}}^{\text{sim}}}\Bigg{\|}\sum_{k=1}^{n_{\text{rcv}}}\bigg{(}d_{v}^{(f_{i})}(\boldsymbol{x}_{k})\,p^{(g_{j})}(\boldsymbol{x}_{k})-d_{p}^{(f_{i})}(\boldsymbol{x}_{k})\,v_{\boldsymbol{\nu}}^{(g_{j})}(\boldsymbol{x}_{k})\bigg{)}\Bigg{\|}^{2}_{2},$$

$$f_k=\delta (\mathbf{x}-{\mathbf{x}}_k),\text{point source in xk.}$$

$$f_j=\sum _{k=1}^{n_{\text{pt}}}\delta (\mathbf{x}-{\mathbf{x}}_{j,k}^{(f)}),\text{multiple-point source.}$$

$$\mathcal{J}_{L2}^{\text{stack}}(m)=\sum_{j=1}^{n_{\text{stack}}}\dfrac{1}{2}\Big{\|}\sum_{k=1}^{n_{\text{pt}}}\mathcal{R}\Big{(}p^{(f_{j,k})}\Big{)}-\sum_{k=1}^{n_{\text{pt}}}d_{p}^{(f_{j,k})}\Big{\|}^{2}+\dfrac{\eta}{2}\Big{\|}\sum_{k=1}^{n_{\text{pt}}}\mathcal{R}\Big{(}v_{\boldsymbol{\nu}}^{(f_{j,k})}\Big{)}-\sum_{k=1}^{n_{\text{pt}}}d_{v}^{(f_{j,k})}\Big{\|}^{2},$$

$$\boldsymbol{x}_{j,\bullet}^{(g)}=\big{\{}\boldsymbol{x}_{i}^{(f)}\big{\}}_{i=i_{1}}^{i_{2}}\,,\qquad\text{ with }i_{1}=(j-1)n_{\text{pt}}+1\quad\text{ and }\quad i_{2}=i_{1}+n_{\text{pt}}-1,$$

$$p{\mid }_{{\Omega }^{+}}=p^{+},p{\mid }_{{\Omega }^{-}}=p^{-},$$

$$-\mathrm{i}\omega \rho (\mathbf{x})\mathbf{v}(\mathbf{x})$$

$$-\mathrm{i}\omega \kappa (\mathbf{x}{)}^{-1}p(\mathbf{x})$$

$$p(\mathbf{x})$$

$${\partial }_{\mathbf{\nu }}p(\mathbf{x})-\frac{\mathrm{i}\omega }{c(\mathbf{x})}p(\mathbf{x})$$

$$p^{+}$$

$${\mathbf{v}}^{+}$$

$$\begin{array}{rl}& {\int }_{{\Omega }^{-}}-\mathrm{i}\omega {\rho }_1{\mathbf{v}}_1\cdot {\mathbf{v}}_2+\nabla p_1\cdot {\mathbf{v}}_2\mathrm{d}{\Omega }^{-}=0,\\ & {\int }_{{\Omega }^{-}}-\mathrm{i}\omega {\kappa }_1^{-1}{p}_1{p}_2+(\nabla \cdot {\mathbf{v}}_1)p_2\mathrm{d}{\Omega }^{-}=0,\end{array}$$

$$\begin{array}{rl}{\int }_{{\Omega }^{-}}\mathrm{i}\omega {\rho }_1{\mathbf{v}}_1\cdot {\mathbf{v}}_2+p_1\nabla \cdot ({\mathbf{v}}_2)& -\mathrm{i}\omega {\kappa }_1^{-1}{p}_1{p}_2-{\mathbf{v}}_1\cdot \nabla p_2\mathrm{d}{\Omega }^{-}\\ & ={\int }_{\partial {\Omega }^{-}}{p}_1({\mathbf{v}}_2\cdot \mathbf{\nu })-({\mathbf{v}}_1\cdot \mathbf{\nu })p_2\mathrm{d}\partial {\Omega }^{-}.\end{array}$$

$$\displaystyle\begin{split}\int_{\Omega^{-}}\mathrm{i}\omega\rho_{1}\boldsymbol{v}_{1}\cdot\boldsymbol{v}_{2}&+p_{1}\cdot\nabla\cdot(\boldsymbol{v}_{2})\,-\mathrm{i}\omega\kappa_{1}^{-1}p_{1}\,p_{2}-\boldsymbol{v}_{1}\cdot\nabla p_{2}\,\,\mathrm{d}{\Omega^{-}}\\ &=\int_{\Omega^{-}}\mathrm{i}\omega\rho_{1}\boldsymbol{v}_{1}\cdot\boldsymbol{v}_{2}+p_{1}\,(\mathrm{i}\omega\kappa_{2}^{-1}p_{2})\,-\mathrm{i}\omega\kappa_{1}^{-1}p_{1}\,p_{2}-\boldsymbol{v}_{1}\cdot(\mathrm{i}\omega\rho_{2}\boldsymbol{v}_{2})\,\,\mathrm{d}{\Omega^{-}}\\ &=\int_{\Omega^{-}}\mathrm{i}\omega(\rho_{1}-\rho_{2})\boldsymbol{v}_{1}\cdot\boldsymbol{v}_{2}+\mathrm{i}\omega\big{(}\kappa_{2}^{-1}-\kappa_{1}^{-1}\big{)}\,p_{1}\,p_{2}\,\,\mathrm{d}{\Omega^{-}}.\end{split}$$

$$\begin{array}{rl}{\int }_{\partial {\Omega }^{-}}{p}_1({\mathbf{v}}_2\cdot \mathbf{\nu })& -({\mathbf{v}}_1\cdot \mathbf{\nu })p_2\mathrm{d}\partial {\Omega }^{-}\\ =& {\int }_{{\Gamma }_2^{-}}{p}_1({\mathbf{v}}_2\cdot \mathbf{\nu })-({\mathbf{v}}_1\cdot \mathbf{\nu })p_2\mathrm{d}{\Gamma }_2^{-}+{\int }_{\Sigma }{p}_1({\mathbf{v}}_2\cdot \mathbf{\nu })-({\mathbf{v}}_1\cdot \mathbf{\nu })p_2\mathrm{d}\Sigma .\end{array}$$

$$\displaystyle\begin{split}&\int_{\Gamma_{2}^{-}}p_{1}\,(\boldsymbol{v}_{2}\cdot\boldsymbol{\nu})-(\boldsymbol{v}_{1}\cdot\boldsymbol{\nu})\,p_{2}\,\,\mathrm{d}\Gamma_{2}^{-}\\ &=\int_{\Gamma_{2}^{-}}p_{1}\,\big{(}(\mathrm{i}\omega\rho_{2})^{-1}\partial_{\boldsymbol{\nu}}p_{2}\big{)}-\big{(}(\mathrm{i}\omega\rho_{1})^{-1}\partial_{\boldsymbol{\nu}}p_{1}\big{)}\,p_{2}\,\,\mathrm{d}\Gamma_{2}^{-}\\ &=\int_{\Gamma_{2}^{-}}p_{1}\,\big{(}(c_{2}\rho_{2})^{-1}p_{2}\big{)}-\big{(}(c_{1}\rho_{1})^{-1}p_{1}\big{)}\,p_{2}\,\,\mathrm{d}\Gamma_{2}^{-}\\ &=\int_{\Gamma_{2}^{-}}\big{(}(c_{2}\rho_{2})^{-1}-(c_{1}\rho_{1})^{-1}\big{)}\,p_{1}\,p_{2}\,\,\mathrm{d}\Gamma_{2}^{-}.\end{split}$$

$$\displaystyle\begin{split}&\int_{\Omega^{-}}\mathrm{i}\omega(\rho_{1}-\rho_{2})\boldsymbol{v}_{1}\,\boldsymbol{v}_{2}+\mathrm{i}\omega\big{(}\kappa_{2}^{-1}-\kappa_{1}^{-1}\big{)}\,p_{1}\,p_{2}\,\,\mathrm{d}{\Omega^{-}}\\ &-\int_{\Gamma_{2}^{-}}\big{(}(c_{2}\rho_{2})^{-1}-(c_{1}\rho_{1})^{-1}\big{)}\,p_{1}\,p_{2}\,\,\mathrm{d}\Gamma_{2}^{-}=\int_{\Sigma}p_{1}\,(\boldsymbol{v}_{2}\cdot\boldsymbol{\nu})-(\boldsymbol{v}_{1}\cdot\boldsymbol{\nu})\,p_{2}\,\,\mathrm{d}\Sigma\,.\end{split}$$

$$\underset{m}{\min }{\mathcal{J}}_{\bullet }(m)\text{ subject to }\sum _i\mathcal{A}(U^{(f_i)})=F_i,\text{for i∈{1,…,nsrc}},$$

$$\mathcal{L}(m,\tilde{U},\tilde{\gamma})=\mathcal{J}_{\bullet}\,+\,\big{<}\mathcal{A}\tilde{U}-F\,,\,\tilde{\gamma}\big{>}_{\Omega}\,.$$

$${\nabla }_m\mathcal{L}(m,\tilde{U}=U,\tilde{\gamma })={\nabla }_m{\mathcal{J}}_{\bullet }.$$

$${\mathcal{A}}^{\ast }\gamma =-{\partial }_U{\mathcal{J}}_{\bullet },$$

$$\nabla_{m}\mathcal{J}_{\bullet}=\nabla_{m}\mathcal{L}(m,\tilde{U}=U,\tilde{\gamma}=\gamma)=\mathrm{Re}\bigg{(}\big{<}\partial_{m}\mathcal{A}\,U\,,\,\gamma\big{>}_{\Omega}\bigg{)}.$$

$$\mathcal{A}^{*}\begin{pmatrix}\gamma_{v_{x}}^{(f_{i})}\\ \gamma_{v_{y}}^{(f_{i})}\\ \gamma_{v_{z}}^{(f_{i})}\\ \gamma_{p}^{(f_{i})}\\ \end{pmatrix}=-\mathcal{R}^{*}\begin{pmatrix}\eta\big{(}\mathcal{R}(v_{\boldsymbol{\nu}}^{(f_{i})})-d_{v}^{(f_{i})}\big{)}\boldsymbol{\nu}_{x}\\ \eta\big{(}\mathcal{R}(v_{\boldsymbol{\nu}}^{(f_{i})})-d_{v}^{(f_{i})}\big{)}\boldsymbol{\nu}_{y}\\ \eta\big{(}\mathcal{R}(v_{\boldsymbol{\nu}}^{(f_{i})})-d_{v}^{(f_{i})}\big{)}\boldsymbol{\nu}_{z}\\ \mathcal{R}(p^{(f_{i})})-d_{p}^{(f_{i})}\end{pmatrix}\qquad\text{adjoint-state problem for $\mathcal{J}_{L2}$. }$$

Full text

1 Introduction

The Full Waveform Inversion (FWI) method has been extensively developed in the last decades for quantitative recovery of sub-surface Earth media in seismic exploration. The concept of FWI is to minimize with respect to the Earth parameters a misfit criterion defined from the simulations of wave propagation and the measured seismograms (i.e., the ‘full waveform’). The method was originally introduced in the work of Bamberger et al., (1977, 1979) for one-dimensional wave equation, and was extended in the work of Lailly, (1983) and Tarantola, (1984); Tarantola, 1987b . The method was first used with time-domain wave propagation, and the frequency-domain formulation of FWI, which requires a Fourier transform of the original time-dependent seismic traces, was established by Pratt et al., (1996, 1998).

In marine seismic, the data usually consist of measurements of the pressure field from hydrophones, but new devices, the dual-sensors, have been recently deployed and also give access to the vertical velocity, see Carlson et al., (2007) and Tenghamn et al., (2007). This additional information has shown advantages to reduce noise for image processing, cf. Whitmore et al., (2010) and Rønholt et al., (2015), and has motivated new analysis, such as in the work of Alessandrini et al., (2018, 2019), and specific numerical methodology by Faucher, (2017); Zhong and Liu, (2018) for the seismic inverse problem. In our work, we implement a misfit functional dedicated to the dual-sensors data and demonstrate its efficiency for the recovery of sub-surface physical parameters.

The FWI relies on an iterative minimization of the misfit functional which, in the traditional least-squares approach, is simply the $L^{2}$ difference between the observations and the simulations. In addition to the numerical challenges induced by the large scale domain, one main difficulty of FWI (which is a nonlinear and ill-posed inverse problem) is that the misfit functional suffers from local minima, in particular when the low frequency data are missing and/or in absence of a priori information, as it has been observed by, e.g., Gauthier et al., (1986); Luo and Schuster, (1991); Bunks et al., (1995) and Faucher et al., 2020a . In order to mitigate the ‘cycle-skipping’ effect, selection of increasing frequency content in the data is commonly employed, cf. Bunks et al., (1995); Sirgue and Pratt, (2004). In our work, we follow the frequency domain formulation, where this approach is natural, and further employ increasing sequential frequency, as advocated by Faucher et al., 2020a .

Several alternatives to the least-squares have been studied to enhance the convexity of the misfit functional, such as logarithmic function, mentioned by Tarantola, 1987a , that is particularly appropriate in complex-frequency FWI, see Shin and Min, (2006); Shin and Cha, (2008) and Faucher, (2017). Comparisons of misfit using the phase and amplitude of signals are carried out in the work of Shin et al., (2007); Bednar et al., (2007) and Pyun et al., (2007). Fichtner et al., (2008) studies the use of a misfit based upon the phase, correlation, or the envelope; the later is also advocated by Bozdağ et al., (2011), and is shown to be efficient for the reconstruction of attenuation properties in global Earth seismology by Karaoğlu and Romanowicz, (2017). Brossier et al., (2010) study the use of $L^{1}$ criterion. The optimal transport distance is considered in the work of Métivier et al., (2016) and Yang et al., (2018) and is shown to improve the misfit functional convexity. In order to avoid the local minima, one can also rely on a specific model parametrization such as the migration-based travel time (MBTT) method which decomposes the velocity model into a (smooth) background profile and a reflectivity in the data-space, see Clément et al., (2001) and Barucq et al., (2019). This approach is shown to increase the size of the attraction basins for background velocity reconstruction by Barucq et al., (2019).

In the case where Cauchy data are available on the boundary, that is, both a field and its normal derivative (e.g., the pressure and the normal derivative of the pressure, which actually relates to the normal velocity, see the appendix) the reconstruction can use a combination of the Dirichlet and Neumann traces, as in the work of Kohn and Vogelius, (1985) and Colton and Haddar, (2005) in inverse scattering, or Alessandrini et al., (2018) for seismic. This approach is labeled as “reciprocity-gap” in the literature, see, e.g., de Hoop and de Hoop, (2000); Colton and Haddar, (2005). Lipschitz-stability can be obtained for the partial Cauchy data inverse problem, in particular in the seismic context with surface measurements, as proved by Alessandrini et al., (2018, 2019). This result further holds for piecewise-linear parameters, employed for the numerical experiments of Alessandrini et al., (2019) and Faucher, (2017). In the context of acoustic waves governed by the Euler’s equations, we shall see that Cauchy data are in fact equivalent to the pressure field and the normal velocity, i.e., dual-sensors data. In the first appendix, we further connect the boundary measurements to the volume properties of the media. Hence, the minimization of the misfit functional defined from surface data is connected to the recovery of the volume medium parameters. In our work, we employ the reciprocity-gap formulation and further combine all sources, defining the Full Reciprocity-gap Waveform Inversion (FRgWI) framework for seismic imaging using dual-sensors data.

In the time-harmonic formulation, the reciprocity-gap results in a misfit functional where observations and simulations are multiplied. While our work is conducted in the frequency domain, it can similarly be carried out in the time domain where it relates to the family of correlation misfits. Correlation-based misfit functionals have been studied in seismic tomography, for instance by Van Leeuwen and Mulder, (2010) and in the work of Choi and Alkhalifah, (2011) and Zhang et al., (2016), however using a single type of measurements (i.e., the acoustic pressure fields). The use of the velocity fields is studied in the time-domain by Zhong and Liu, (2019), assuming that all directional velocities are available and convolving the same fields. A fundamental difference with these earlier studies, and in particular with Zhong and Liu, (2019), is that in our work, we correlate different fields (i.e. pressure with velocity). This is the essence of reciprocity-gap and it is crucial in order to relate surface measurements to global model reconstruction using Green’s identity (see the first appendix). In addition, our framework is naturally designed for the normal velocity, in accordance with the dual-sensors devices. In the work of Menke and Levin, (2003) and Bodin et al., (2014) in seismology, observations and simulations are also combined, but using a cross-convolution formula made of the vertical and horizontal components of the waves. More generally, approaches based upon a specific filtering of the data, such as in the work of Warner and Guasch, (2016) and Guasch et al., (2019), also rely on a minimization where the fields are not directly the quantities compared.

The main feature of FRgWI is to allow different observational and numerical acquisitions, cf. Faucher, (2017); Alessandrini et al., (2019); Zhong and Liu, (2019). Minimal information regarding the observational sources is required: the source function and the source positions are not needed to conduct the reconstruction. This is due to the definition of the misfit functional which does not compare an observation with a simulation directly, but products of an observation with a simulation. This flexibility in the choice of numerical probing sources opens up many perspectives and, in particular, to use denser or sparser computational acquisitions than the given observational one. The use of sparse acquisition relates to the shot-stacking approach for data decimation, which sums several single point-source data using the linearity of the wave equation. It has early been employed in seismic to reduce the computational cost by Mora, (1987). It is based upon the redundancy of information in the data, but has a major drawback as it is difficult to avoid the cross-talk between the encoded sources, cf. Krebs et al., (2009); Zhang et al., (2018). It has motivated several works to efficiently assemble the encoded sources in seismic, i.e., source blending (Berkhout,, 2008). We mention, for instance, the random combination of sources changing with iteration by Krebs et al., (2009), the approach based upon compress sensing by Li et al., (2012), while wavelet encoding is used by Zhang et al., (2018); see also the references therein. In our applications, we will see that FRgWI, by using arbitrary numerical sources, can work with multiple-point sources (for computational or observational acquisition) while being naturally robust to the cross-talk. In fact, it is not exactly cross-talk in this context because the key of FRgWI is that the measurements are not modified, i.e., they are still tested independently, one by one, with respect to the numerical simulations.

In this paper, we study the seismic inverse problem associated with time-harmonic waves, using dual-sensors data. We first state the mathematical problem and define the two misfit functionals that we analyze for the iterative reconstruction: the traditional least-squares difference and the reciprocity-based functional. Next, we detail the features provided by the reciprocity-gap formulation. We carry out three-dimensional reconstruction experiments: first with a layered medium, where we compare the performance of reciprocity-gap and full waveform inversion. We further demonstrate the efficiency of FRgWI with simultaneous point-sources to reduce the computational cost, and its robustness compared to traditional shot-stacking. Finally, we carry out a larger-scale experiment including salt domes using the SEAM benchmark. In particular, we highlight that FRgWI performs equally well in the case of acquisitions that are sparse for observations/dense for simulations and when acquisitions are dense for observations/sparse for simulations. Our experiments are carried out using data acquired in the time-domain and, even though our reconstruction algorithm is conducted in the frequency-domain, which might present some scale limitations for field configurations, we believe that our method can be implemented with the time-domain wave equation. This can be more appropriate for larger scales and it only requires a slight modification to our model. A discussion about it is carried out in the conclusions.

2 Time-harmonic seismic inverse problem with dual-sensors data

We work with time-harmonic wave propagation for the identification of the physical parameters in a seismic context. The quantitative reconstruction is conducted using an iterative minimization of a misfit functional. We first give the misfit as the traditional $L^{2}$ difference and further design the reciprocity-gap version of the functional, which combines pressure and normal velocity data.

2.1 Acoustic wave equation, forward problem from dual-sensors

We consider a three-dimensional domain $\Omega\subset\mathbb{R}^{3}$ with boundary $\Gamma$. The propagation of waves in acoustic media is represented with the scalar pressure and vectorial velocity fields, that satisfy the Euler’s equations (Colton and Kress,, 1998; Kirsch,, 1996; de Hoop and de Hoop,, 2000),

[TABLE]

The frequency is denoted by $\omega$, $f$ is the (scalar) interior (harmonic) source term and $\partial_{\boldsymbol{\nu}}$ denotes the normal derivative. The propagation is characterized by the two physical parameters of the medium: the density $\rho$ and the bulk modulus $\kappa$. In addition, the velocity $c$ is defined such that

[TABLE]

For boundary conditions, we follow a geophysical context where the boundary is separated into two: $\Gamma=\Gamma_{1}\cup\Gamma_{2}$. We denote by $\Gamma_{1}$ the interface between the air and the acoustic medium, where a free surface boundary condition holds, Eq. 1c. On the other part of the boundary, $\Gamma_{2}$, we implement an Absorbing Boundary Conditions (ABC), Engquist and Majda, (1977), Eq. 1d, to ensure that waves reaching the boundary are not reflected back into the domain. It corresponds to the fact that the domain $\Omega$ is a numerical restriction of the Earth.

The dual-sensor devices have been recently introduced in marine seismic exploration, and allow the recording of both the pressure field and the vertical velocity, see Carlson et al., (2007) and Tenghamn et al., (2007). Consequently we define the forward problem (which maps the parameters to the data) $\mathcal{F}_{\omega}^{(f)}$ at frequency $\omega$ for a source $f$ such that

[TABLE]

The model parameters are referred to by $m=(\kappa,\,\rho)$, the normal velocity by $v_{\boldsymbol{\nu}}$ and $\Sigma$ corresponds with the (discrete) set of receivers location:

[TABLE]

where $\boldsymbol{x}_{i}$ is the position of the $i^{\text{th}}$ receiver for a total of $n_{\text{rcv}}$ receivers. For notation, we introduce the restriction operator $\mathcal{R}$, which reduces the fields to the set $\Sigma$, such that

[TABLE]

2.2 Misfit functionals

The inverse problem aims the quantitative reconstruction of the subsurface medium parameters ($\kappa$ and $\rho$) from data measured at the receivers location. Using pressure and vertical velocity measurements, we denote the data at frequency $\omega$ for a source $f$ by

[TABLE]

where $d_{\omega,p}^{(f)}$ and $d_{\omega,v}^{(f)}$ are vectors of $\mathbb{C}^{n_{\text{rcv}}}$ and respectively refer to the pressure and normal velocity records, following Eq. 3. We further denote by $d_{\omega}^{(f)}(\boldsymbol{x}_{i})$ the data recorded for the source $f$ at the $i^{\text{th}}$ receiver.

In the following, we omit the frequency index and space dependency for the sake of clarity (in the experiments, we use increasing sequential frequency, as suggested by Faucher et al., 2020a ). We introduce two misfit functionals which evaluate the difference between the observations and simulations.

Firstly, the functional $\mathcal{J}_{L2}$, which follows the traditional least-squares, is

[TABLE]

where $\eta$ is a scaling factor to adjust between the amplitudes of the pressure and velocity. In our applications, it is taken such that $\|d_{p}\|=\eta\|d_{v}\|$.

Secondly, we define an alternative misfit functional based upon the reciprocity-gap, motivated by Green’s identity:

[TABLE]

where T denotes the transpose. The misfit functional is motivated by the Green’s identity and has been introduced in the context of inverse scattering, from Cauchy data, as mentioned in the introduction, cf. Kohn and Vogelius, (1985); Colton and Haddar, (2005), and used with partial seismic data by Alessandrini et al., (2019). In the first appendix, we justify the formulation of the misfit functional $\mathcal{J}_{\mathfrak{r}}$ using variational formulation of Problem 1, and note that the mathematical foundation of the reciprocity-gap formulation relies naturally on the normal velocity. Therefore, it perfectly matches the dual-sensors data, and does not necessitate the specific directional components ($v_{x}$, $v_{y}$ or $v_{z}$).

2.3 Iterative reconstruction procedure

The reconstruction procedure follows an iterative minimization of the selected misfit functional. For least-squares formulation such as Eq. 7, it is traditionally referred to as the Full Waveform Inversion (FWI) method in seismic (as one makes use of the full data seismograms) see the review of Virieux and Operto, (2009) and the references therein. Consequently, we shall refer to the minimization of Eq. 8 as Full Reciprocity-gap Waveform Inversion: FRgWI. In any case, the iterative minimization follows successive updates of the physical models, such that, at iteration $l$, the new model is given by $m_{l+1}=m_{l}-\alpha_{l}s_{l}$, where $\alpha$ is the scalar step size typically computed via a line search method (Nocedal and Wright,, 2006), and $s$ is the search direction.

The search direction depends on the gradient of the cost function, which is computed using the adjoint state method. The method has its foundation in the body of work of Lions, (1971), and is reviewed for seismic application by Plessix, (2006). Application with complex fields is further described by Barucq et al., (2018, 2019). The adjoint-state method for reciprocity-gap functional is briefly reviewed in the second appendix, see also Alessandrini et al., (2019). In our implementation, the search direction only depends on the gradient of the misfit functional, and we employ the Limited-BFGS algorithm, see Nocedal, (1980) and Nocedal and Wright, (2006). We review the steps for the iterative minimization in Algorithm 1.

2.4 Discretization with Hybridizable Discontinuous Galerkin method

In the numerical implementation, both the pressure and velocity fields must be computed to feed the misfit functional. We discretize Problem 1 using the Hybridizable Discontinuous Galerkin (HDG) method, which is specifically designed for first-order problems, and avoid oversize linear systems. Indeed, the global matrix using HDG is only composed of the degrees of freedom associated with the trace of the pressure field, that is, only the degrees of freedom on the faces of the elements of the discretization mesh, see Cockburn et al., (2009); Griesmaier and Monk, (2011) and Bonnasse-Gahot et al., (2017). Then, local (small) systems are solved to calculate the volume solutions of both the pressure and the velocity fields, computed with similar accuracy. We refer to Faucher and Scherzer, (2020) for the precise implementation.

With more traditional discretization methods such as Continuous Galerkin, Finite Differences of Internal Penalty Discontinuous Galerkin methods, one has to create a linear system which size is the total number of degrees of freedom for all unknowns (i.e., the pressure and the three components of the velocity), possibly leading to cumbersome systems. Alternatively, one can solve only for the scalar pressure field, and post-process the solution to obtain the velocity but the computed velocity looses one order of accuracy compared to the discretized pressure field due to the derivative in Eq. 1a. In the HDG method, we obtain both the pressure and velocity fields with the same accuracy, while the global linear system only contains the degrees of freedom of a scalar unknown (Faucher and Scherzer,, 2020). In addition, only the degrees of freedom on the faces of the elements are taken into account, hence removing all interior ones. Consequently, the HDG method has been shown to be more efficient (i.e., less memory consumption for the matrix factorization due to the smaller global matrix) compared to other approaches by Kirby et al., (2012) and Bonnasse-Gahot et al., (2017), for high order discretization.

In the context of large scale seismic applications, the use of HDG is crucial (particularly in the frequency domain) to efficiently compute the pressure and the velocity, because it eventually leads to linear systems which sizes are not larger than for second-order wave problems. Then, For the resolution of the subsequent linear system, we use a direct solver (Amestoy et al.,, 2001; Liu et al.,, 2020) for the multiple right-hand sides feature, in particular, the solver Mumps, designed for sparse matrices (Amestoy et al.,, 2006). Therefore, the matrix factorization of the forward problem is reused for the backward one that serves to compute the gradient of the misfit functional depicted in the second appendix.

3 Features of FRgWI and data acquisition

Following Faucher, (2017) and Alessandrini et al., (2019), the fundamental feature of the reciprocity-gap misfit functional Eq. 8 compared to Eq. 7 is that the set of computational sources is separated from the observational ones (respectively $g_{j}$ and $f_{i}$ in Eq. 8). It implies that

we do not have to know the observational source positions (the location of the $f_{i}$ in Eq. 8) for the reconstruction algorithm; 2. 2.

we do not have to know the observational source signatures (wavelet);

Consequently, and contrary to least-squares misfit such as Eq. 7, minimal information is required regarding the observational acquisition (only the positions of the receivers). The recovery of the source wavelet in parallel to the iterative minimization procedure usually performs well in seismic exploration, using the method prescribed by Pratt, (1999), see Remark 1. However, incorrect knowledge of the position of the sources is a strong difficulty which can lead to the failure of the reconstruction procedure. Namely, FRgWI increases robustness by being free of such considerations. In addition,

the set of computational sources can be arbitrarily taken, hence can differ from the observational one (respectively $n_{\text{src}}^{\text{obs}}$ and $n_{\text{src}}^{\text{sim}}$ in Eq. 8).

Therefore, it provides high flexibility regarding the choice of computational sources. It is important to notice that whatever computational sources are selected, the observed data are still independently tested one by one against the simulations, i.e., the measurements are not modified in any way.

3.1 Data acquisition

In our experiments, we consider the source term for the wave propagation ($f$ or $g$) to be a Ricker function in time and a delta-Dirac function in space. In the frequency domain, it translates to a delta-Dirac function $\delta$ (with value equal to the discrete Fourier transform of the time-domain signal). We refer to a point-source when the source is localized at the position $\boldsymbol{x}_{k}$:

[TABLE]

We further refer to multiple-point source when it is composed of $n_{\text{pt}}$ point-sources, such that

[TABLE]

The observational setup uses a fixed lattice of receivers, and sources located slightly above, as illustrated in Figure 1a. Furthermore, the position of the receivers remains the same for all sources. This configuration is consistent with the TopSeis acquisition system for marine seismic developed by CGG (Compagnie Générale de Géophysique) and Lundin Norwary AS. This principle has been recently deployed on field, showing benefits for near-offset coverage, and it has also won the ‘Exploration Innovation Prize 2019’111We refer to https://www.cgg.com/en/What-We-Do/Offshore/Products-and-Solutions/TopSeis and https://expronews.com/2019/05/23/topseis-a-worthy-winner/. . Namely, it consists in having two boats: one that moves and carries the sources, and one that carries the receivers and remains fixed.

For the reconstruction using FRgWI, we can employ arbitrary numerical sources, which do not have to coincide with the observational acquisition. Here, we investigate the two following configurations.

3.2 Dense point-source computational acquisition for sparse multiple-point measurements

In the case where the observational acquisition is composed of few sources (e.g., to reduce the cost of field experiments), FRgWI can use a dense coverage of computational point-sources to enhance the sensitivity. The sources for the measurements can be multiple-points, e.g., if several air-guns are excited at the same time. In our experiments, the choice of positions for the multiple-points follows a structured decomposition: adjacent sources are grouped, with a fixed distance in $x$ and $y$ between each of them, cf. Figure 1b. It appears more natural to employ such structured decision, e.g. the boat will carry the air-gun sources all together. In addition, note that we have observed in our experiments that this structured partition gives a better performance than using a random combination of sources.

3.3 Sparse multiple-point computational acquisition for dense point-source measurements

In the case where the observational acquisition is composed of many (point) sources (as it can happen in exploration seismic), FRgWI can instead use multiple-point sources in the computational acquisition to reduce the computational cost. It consists of a sparsification of the observational acquisition.

In order to compare with the traditional FWI approach, one can use the linearity of the wave equation and the multiple-point source can be assimilated to the well-known shot-stacking approach, which rewrites the misfit functional Eq. 7 such that

[TABLE]

where we have a total of $n_{\text{stack}}$ multiple-point sources, each composed of $n_{\text{pt}}$ points.

Therefore, we will investigate in the following numerical sections different situations where one of the acquisitions (observational or computational) is composed of point-sources while the other is made of multiple-point sources. The main advantage of FRgWI is that it does not require any modification of the data, which are tested one by one against the computational acquisition in Eq. 8. However, the shot-stacking version of FWI requires the summation of data in Eq. 11, which results in the possible loss of information, i.e., cross-talk (Zhang et al.,, 2018).

4 Numerical experiment 1: layered medium

In this section and the next, we carry out three-dimensional experiments of geophysical reconstruction to study the performance of FRgWI. In this first test, we compare with the traditional least-squares functional, and we employ multiple-point sources to probe the robustness of arbitrary source positions; we also test the combination of sparse and dense acquisitions.

4.1 Velocity model

We consider a three-dimensional acoustic medium (provided by Statoil) of size $2.55\times 1.45\times 1.22$ ${\mathrm{km}}^{3}$, in the $x$, $y$, and $z$ axes respectively (i.e., $z$ is the medium depth). We assume a constant density $\rho=1000~{}$\mathrm{kg}\text{,}{\mathrm{m}}^{-3}$$, and the subsurface velocity is pictured in Figure 2. It consists of different geophysical layers, with non-monotone variations, from $1500$ to $5200$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. The first 500 $\mathrm{m}$ are mostly constant with a velocity of about $1600$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$.

For the iterative reconstruction, we start with a one-dimensional velocity profile that is pictured in Figure 3. This initial guess does not encode any a priori information, and only has an increasing velocity with the depth. Moreover, the range of values is lower compared to the true medium of Figure 2.

4.2 Time-domain data with noise

We work with time-harmonic wave propagation but follow a seismic context where measurements are obtained in the time-domain. Furthermore, we incorporate noise in the synthetic seismograms to have a more realistic setup. The observational acquisition consists of 160 point-sources, excited one by one. They are positioned at the depth $z=10$ $\mathrm{m}$, and forms regular lattice such that sources are every $160$ $\mathrm{m}$ along the $x$ axis, and every $150$ $\mathrm{m}$ along the $y$ axis. There is a total of $1376$ receivers, which are located at a fixed depth of $100$ $\mathrm{m}$, every $60$ $\mathrm{m}$ along the $x$ axis, and every $50$ $\mathrm{m}$ along the $y$ axis. Note that the receivers remain at the same position for all sources, cf. Figure 1.

We generate the time-domain seismograms222We have used the parallel time-domain code Hou10ni, see https://team.inria.fr/magique3d/software/hou10ni/; it relies on Internal Penalty Discontinuous Galerkin discretization (while we use HDG). Also, the meshes are different between the time-domain modeling and the harmonic inversion. and incorporate noise in the resulting traces, with a signal-to-noise ratio of $10$ $\mathrm{dB}$. Then, we proceed with the discrete Fourier transform to feed Algorithm 1. In Figure 4, we picture the noisy time-domain pressure trace for a single-point source, and the corresponding Fourier-transform that we employ for the reconstruction. We respect the seismic constraint that the low-frequencies are not available from the time-domain data (because of noise and listening time) and we only work with data from $5$ to $15$ $\mathrm{Hz}$ frequency. In Figure 5, we picture two-dimensional time-space sections of the traces for the pressure and normal velocity, and illustrate the effect of the added noise (10 $\mathrm{dB}$ signal-to-noise ratio in our experiments).

4.3 Performance comparison of the misfit functionals

We first compare the performance of both misfit functionals, $\mathcal{J}_{L2}$ and $\mathcal{J}_{\mathfrak{r}}$, in the same context: the numerical acquisition is taken to be the same as the observational one (which is anyway mandatory for $\mathcal{J}_{L2}$). Therefore, we take $n_{\text{src}}^{\text{sim}}=n_{\text{src}}^{\text{obs}}$ and the same set for $f$ and $g$ in Eq. 8. We use Algorithm 1, using sequential frequency progression from $5$ to $15$ $\mathrm{Hz}$, more precisely, we use $\{5,\,6,\,7,\,8,\,9,\,10,\,12,\,15\}$ $\mathrm{Hz}$ data, and $30$ minimization iterations per frequency. In Figures 6 and 7, we show the final reconstruction, after $15$ $\mathrm{Hz}$ iterations when minimizing $\mathcal{J}_{L2}$ and $\mathcal{J}_{\mathfrak{r}}$ respectively. For the discretization, we employ a mesh of about $75$ thousands tetrahedra, and polynomials of order three to five (depending on the frequency) for accuracy (note the mesh differs from the one employed to generate the time-domain data).

Remark 1**.**

For the minimization of $\mathcal{J}_{L2}$ Eq. 7, one has to use the same source wavelet for the observations and simulations. Because the observational source wavelet is not precisely known, we need to reconstruct the source during the iterative process as well. We employ the update formula given by Pratt, (1999) for the iterative point-source reconstruction. This can however induce additional difficulty in the case where the source characteristic is not well recovered or, more important, when the source positions are not precisely known.

In this experiment, we use the same acquisition context for the two choices of misfit functional, to observe intrinsic differences between the two methods. We observe the following:

• –

both approaches provide a good reconstruction, where the layers of velocity appear (see the vertical section), and the correct speed values are retrieved.

• –

The full reciprocity-gap waveform inversion provides slightly better results, in particular for the parts that are near the boundaries (see the vertical and horizontal sections of Figures 6 and 7). This can be explained as FRgWI formulation Eq. 8 tests every simulation source with each observational one (the two sums in Eq. 8) and it somehow compensates for the limited illumination of the boundary zones.

• –

Regarding the computational time, both approaches are similar (the only difference is the two sums in the misfit functional Eq. 8, which is a computationally cheap operation compared to the matrix factorization and linear system resolution).

Therefore, in this first test-case, we have demonstrated the efficiency of FRgWI, which produces better reconstruction compared to the traditional approach, without incurring any increase of computational time. But more important, the reciprocity-gap offers flexibility for the choice of computational acquisition, which we now study.

4.4 Comparison of multiple-point sources FRgWI with shot-stacking

The main feature of FRgWI is to enable the use of arbitrary probing sources for the numerical simulations ($g$ in Eq. 8). Here we investigate the use of multiple-point sources in order to reduce the computational cost. The observational acquisition is composed of $160$ point-sources, i.e., each source function $f_{i}(\boldsymbol{x})$ in Eq. 8 and Eq. 7 corresponds with a delta-Dirac function in $\boldsymbol{x}_{i}^{(f)}$, according to Eq. 9. For the computational acquisition ($g$ in Eq. 8), we now consider $n_{\text{stack}}$ multiple-point sources, each composed of $n_{\text{pt}}$ points, cf. Eq. 10.

We assume that the multiple-point sources all have the same number of points. The positions of the multiple-points are taken to coincide with the observational acquisition, and we consider a group of sources in a structured partition of the original configuration, as illustrated in Figure 1. We have

[TABLE]

where $\boldsymbol{x}_{j}$ corresponds with the multiple-point source, composed of the point-sources located in $\boldsymbol{x}_{i}$. The FWI counterpart is obtained using shot-stacking, following Eq. 11, i.e., it requires the summation of the observed data.

We perform the iterative reconstruction using $n_{\text{stack}}=5$ multiple-point sources, each of them composed of $n_{\text{pt}}=32$ points. The reconstruction using the shot-stacking version of FWI Eq. 11 is shown Figure 8 and the result using FRgWI (i.e., where only the computational acquisition is changed) is shown Figure 9.

We observe that FRgWI performs better than the shot-stacking in FWI, in this experiment of relatively small scale:

• –

the vertical and horizontal sections shown in Figure 9 are actually very close from the reconstruction obtained using the same computational and observational acquisition, see Figure 7. The pattern of increasing and decreasing velocity is well captured, with the appropriate values.

• –

On the other hand, FWI with shot-stacking looses accuracy, in particular the vertical section in Figure 8 where the non-monotone variation is not even recovered.

FRgWI, by testing each of the observational sources independently with the arbitrary computational ones, is robust with multiple-point sources, and allows a sharp decrease in the number of numerical sources for the simulations.

Remark 2**.**

The results using the shot-stacking with FWI could be improved by using more advanced strategy of shot blending as mentioned in the introduction (e.g., Krebs et al., (2009); Li et al., (2012)). We illustrate here that even a naive approach (i.e. simple to implement numerically) of shot-stacking is sufficient for the full reciprocity-gap waveform inversion to produce satisfactory results. Also, note that we have tried to use a random selection of positions for the multiple-point but it was not performing as well as the structured criterion. This would however require more testing to confirm.

We can further reduce the number of computational sources. The reconstruction using FRgWI with $n_{\text{stack}}=2$, $n_{\text{pt}}=80$ is shown in Figure 10. In Figure 11, we show the reconstruction using only one source: $n_{\text{stack}}=1$, $n_{\text{pt}}=160$. We see that FRgWI still behaves quite well, in particular for two sources the reconstruction carries similar accuracy than before: the vertical structures are accurately discovered, with pertinent variations of velocity. For the reconstruction using only one source, Figure 11, we see some deterioration in the recovered layers but it remains more precise than the shot-stacking FWI using $5$ sources (Figure 8).

The essence of FRgWI is that it does not modify the observational acquisition and instead it probes every measured seismogram independently. It increases the robustness of the reconstruction procedure, and allows for the design of arbitrary computational acquisitions to reduce the computational cost. In this experiment, we have used a naive approach of shot summation, which requires no effort for its implementation, and it already shows accurate reconstructions.

4.5 Sparse observational acquisition

We now investigate a different situation: when the measurements are obtained from only a few set of multiple-point sources. This happens, for example, in marine seismic when several air guns are excited at the same time, i.e. the source boat carries an array of air guns. It consequently reduces the amount of data obtained, and the cost of the field acquisition.

Hence, we now consider that the measurements are obtained from $5$ multiple-point sources. we follow the configuration of the previous subsection but inverting the computational and observational acquisitions: it is now $f$ in Eq. 8 that is represented with a multiple-point sources Eq. 12 while the computational sources $g$ consist of $160$ single point sources. It means that the observed measurements are acquired from the physical phenomenon corresponding with multiple-point sources. Note that from this type of measurements, the FWI reconstruction coincides with the shot-stacking result of Figure 8. In Figure 12, we show the reconstruction where the observations are obtained from $n_{\text{stack}}=5$ multiple-point sources and the computational acquisition consists of $160$ single-point sources. In Figure 13, we show the reconstruction where the measured data result from one multiple-point source, composed of $160$ points (i.e., all sources are excited at the same time) while the numerical acquisition remains with $160$ point-sources.

We observe that the reconstructions obtained from a drastically reduced set of measurements, resulting from multiple-point sources, retrieve the appropriate velocity variation. It displays similar accuracy compared to the use of single-point source measurements with multiple-point source simulations (respectively Figures 9 and 11 for $5$ and $1$ computational sources). Namely, it seems that FRgWI is insensitive to which acquisition is made sparse (i.e., with multiple-point sources).

5 Numerical experiment 2: salt-body SEAM model

In this experiment, we consider a three-dimensional velocity model encompassing salt-domes. The velocity is extracted from the seismic SEAM333SEG Advanced Modeling Corporation, see https://seg.org/News-Resources/Research-and-Data/SEAM benchmark, and is of size $7\times 6.5\times 2.1$ ${\mathrm{km}}^{3}$. The velocity model is depicted in Figure 14, and varies from $1500$ to $4800$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. In this experiment, the density is heterogeneous, and pictured in Figure 15. Compared to previous test-case, it is of different nature (salt-domes), it is larger and it has a heterogeneous density. Therefore, we expect this numerical experiment to be more challenging.

5.1 Time-domain data

Similarly to the previous experiment, we work with noisy time-domain data and compute their Fourier transform to generate the harmonic data used in the iterative reconstruction. There is a total of $272$ point-sources in the observational acquisition, which are located at $10$ $\mathrm{m}$ depth. The sources are placed on a two-dimensional plane with $400$ $\mathrm{m}$ between each source along the $x$ and $y$-axes. For multiple-point acquisition, we follow the structured combination illustrated in Figure 1b. In this experiment, we generate the time-domain data and we use $20$ $\mathrm{dB}$ signal-to-noise ratio. The resulting seismograms are illustrated in Figure 16 for a single source, where we show the pressure and vertical velocity measurements. The data are acquired by a total of $2805$ receivers. For the reconstruction, we perform a Fourier transform of the time-domain noisy data, and use frequency content from $2$ to $7$ $\mathrm{Hz}$.

5.2 Reconstruction using FRgWI

For the reconstruction, we start with initial guesses that correspond with smooth version of the true models, they are shown in Figures 17 and 18 for the starting velocity and density respectively. We actually only focuses on the reconstruction of the velocity model, and keep the density as its initial representation of Figure 18. The density is known to be more complicated to recover than the velocity because of lack of sensitivity in the data, Jeong et al., (2012), but it should not prevent us from recovering the velocity, see, e.g., Faucher, (2017).

We use the FRgWI algorithm with the minimization of Eq. 8 using sequential frequency content from $2$ to $7$ $\mathrm{Hz}$. We perform $30$ iterations per frequency which leads to a total of $180$ iterations. We analyze the two following configurations.

• –

We first use the dense observational acquisition composed of $272$ point-sources, and multiple-point sources for the computational acquisition ($g$ in Eq. 8). The multiple-point sources use $n_{\text{pt}}=8$ points, resulting in $n_{\text{stack}}=34$ sources. The reconstructed velocity after the iterations at $7$ $\mathrm{Hz}$ is shown in Figure 19.

• –

Next, we consider the opposite situation: we assume a sparse observational acquisition of $n_{\text{stack}}=34$ sources each composed of $n_{\text{pt}}=8$ points. Here, the computational acquisition is made dense, with $272$ point-sources. The recovered velocity is shown in Figure 20.

The minimization performs well in both situations: the salt dome is appearing with appropriate values (about $4500$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$). The upper boundary of the salt is accurately recovered (cf. the right of Figures 19 and 20) while the deepest parts are more difficult to obtain. It is possible that the high velocity contrast prevents us from imaging below the salt. The horizontal section (bottom of Figures 19 and 20) also shows the appropriate pattern, in particular with the decrease of velocity in the middle, and only the parts near the boundary are slightly less accurate. It confirms the results of our first experiment that FRgWI is relatively insensitive to whether the observational or computational acquisitions is taken sparse or dense: the reconstructions display similar resolution for the two cases. In terms of numerical cost, sparse computational acquisition reduces the computational time but, in the frequency domain with the use of direct solver, this gain remains marginal. In order to further improve the reconstruction, we shall include a regularization parameter in the minimization, e.g. with Total Variation approach, which can indeed be supported by our misfit functional; see also Faucher et al., 2020b and the references therein.

6 Analysis of results and perspectives

We have implemented the iterative minimization method based upon the reciprocity-gap functional in the frequency domain, with the hybridizable discontinuous Galerkin (HDG) discretization to efficiently address the first-order system. To probe the method, we have carried out two three-dimensional experiments, one with a velocity made of layers, and another with a medium containing a salt dome, extracted from the SEAM benchmark. Firstly, our experiments have shown that the Full Reciprocity-gap Waveform Inversion (FRgWI) method performs better than FWI in the same configuration, that is, when we keep the observational acquisition for the numerical one, in particular by improving the illumination on the sides.

These promising results motivate the implementation of the method for larger, field benchmarks, in order to probe the scalability of FRgWI. In this case, the time-domain formulation should be considered to continue the performance analysis we have started.

6.1 Flexibility in the computational acquisition

The main feature of the FRgWI is its robustness regarding the absence of information on the observational acquisition, and the consequent freedom in the computational one. This is due to the misfit functional Eq. 8 that works with a product of a datum with a simulation, and that relates with the correlation-based family of methods. For instance, the wavelet used for the computational source does not have to match the observational one, as long as there is an overlap in their frequency contents, and we have not observed any other requirement on the choice of the source function. The use of less traditional sources such as plane waves is theoretically possible, and should be investigated in the future.

Our main result comes from the positions of the sources used for the simulations that can be arbitrarily chosen. This flexibility is in the core of the misfit functional which does not compare directly an observation with a simulation, but tests the product of any observation with any simulation, hence enforcing the domain illumination. It is crucial to note that the measurements are not modified and are independently tested one by one against each of the (chosen) computations. In our work, we remain with illumination from one side (i.e., simulations and observations are generated near the surface), motivated by Alessandrini et al., (2019) which prove the stability of the method. This choice is also motivated from the derivation of the variational formulation we give in appendix, to avoid the appearance of additional integrals (see the discussion at the end of the first appendix). The use of probing sources from other part of the domain (e.g., deeper in the domain) has to be carefully investigated, and is part of our ongoing research.

6.2 Using sparse acquisitions

The use of multiple-point sources allows for the reduction of the computational time or of the field acquisition process. Here, the key is that only one of the two acquisitions is reduced, while the other remains dense. Once again, as every datum in the dense acquisition is tested with respect to each of the sparse one, it appears to overall produce the proper illumination of the domain. This is possible due to the misfit formula which offers high flexibility, contrary to more traditional approaches where both acquisitions must be the same.

We have experimented with the two configurations (sparse computational acquisition with dense observational one or the opposite) and it results in similar accuracy in the reconstruction. Because FRgWI tests independently all combinations of an observation with a simulation, it further appears robust with respect to shot-stacking. With FRgWI, we do not exactly observe the cross-talk effect, as the measurements are not modified before entering the misfit, instead, we observe some inherent difficulties when the acquisition is too sparse. This is unavoidable, for instance, if the original field data have been obtained from a few sources only: we cannot artificially enhance the information contained in the measurements. On the other hand, the difficulties that come from the summation of data (i.e., having multiple-point sources in a sparse acquisition) and that results in cross-talk are diminished by the misfit functional which takes the product with each datum in the dense acquisition.

In the context of shot-stacking, there is a “reference” solution which is the use of the non-summed fields. Here, we raise the following question for future applications of the FRgWI method: how to select the numerical acquisition to make the best use of the observed measurements?

6.3 Data coverage

It is important to remind that the reciprocity-gap misfit functional relates, in terms of receivers, to the discretization of a surface integral, as depicted in the first appendix. In our experiments, we have a dense array of receivers, according to the standards in exploration geophysics. In the case where the surface is not properly covered, we certainly need to appropriately weight the records, using, for instance, existing techniques from Earth seismology, such as the one prescribed by Montagner et al., (2012). Similarly, the use of specific weight (e.g., from quadrature rules) could give a first indication on how to select the computational acquisition in order to ensure the equal illumination of the domain with minimal points.

6.4 Towards practical applications

Our reconstruction algorithm has been carried out in the frequency domain, in which it is difficult to address the scales encountered in field applications, due to the operations of linear algebra (mainly the memory required for the matrix factorization). While the efficiency of the frequency domain is improved by the use of HDG discretization (to reduce the size of linear systems) and by new techniques developed for the solver (Liu et al.,, 2020), time or hybrid domain inversion should be implemented to further probe the method in field studies. In this case, the multiplication of fields in the misfit functional Eq. 8 turns into convolution products, maintaining the same features regarding the flexibility of the numerical acquisition. We do not expect a difference in the performance of the method when implemented in the time-domain for the situation we have studied here, especially if one follows the (usual) progression of increasing frequency-content in the data, as introduced by Bunks et al., (1995). On the other hand, it is crucial to experiment on larger, field configurations, to further evaluate the method compared to the existing techniques. Here, the impact of the (limited) data coverage has to be carefully investigated.

7 Conclusion

In this paper, we have performed waveform inversion using a reciprocity-gap misfit functional which is adapted to the dual-sensors devices which capture both the pressure field and the normal velocity of the waves. In this case, the measurements can be seen as Cauchy data, motivating our choice of misfit from the point of view of inverse scattering theory. Our method enables the use of different acquisitions for the observations and the simulations, by relying on a product of data. We have investigated its performance with sparse acquisitions made of multiple-point sources and have demonstrated its efficiency compared to the traditional FWI. We have used a sparse acquisition either for the observational one (to reduce the cost of the field acquisition) or for the computational one (to reduce the numerical cost of the iterative minimization), while not altering the resolution of the reconstruction.

We have implemented the method in the frequency domain and carried out two experiments of different natures. These represent the first steps of validation, while time and hybrid-domain inversion should now be considered. The implementation of the misfit with time signals does not result in any technicality and should not lead to a difference in performance for the scales we have studied with our frequency-domain setup. However, it is necessary to further probe the behaviour of FRgWI in larger, practical test-cases, and to analyze how it competes with state-of-the art methods.

The FRgWI method allows for arbitrary computational source positions that now needs to be further analyzed. For instance, the definition of criteria for the design of the computational acquisition, and the study of the source positions (e.g., not only close to the near surface area, with quadrature rules to ensure equal illumination) is the subject our future research. Note also that the concept of full reciprocity-gap waveform inversion extends readily to elasticity.

Acknowledgments

The authors thank Jean Virieux and Andreas Fichtner for thoughtful and encouraging discussions. FF is funded by the Austrian Science Fund (FWF) under the Lise Meitner fellowship M 2791-N. The work of GA and ES was performed under the PRIN grant No. 201758MTR2-007. The research of HB has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 777778. The research of HB and FF is also supported by the Inria–TOTAL strategic action DIP. MVdH was supported by the Simons Foundation under the MATH + X program, the National Science Foundation under grant DMS-1815143, and by the members of the Geo-Mathematical Imaging Group at Rice University. The work of RG was supported by the Higher Education Authority (HEA) Government of Ireland International Academic Mobility Program 2019. ES has also been supported by the Individual Funding for Basic Research (FFABR) granted by MIUR. The numerical experiments have been performed on the CEA’s cluster Irene, as part of the GENCI allocation AP010411013; and on the cluster PlaFRIM (Plateforme Fédérative pour la Recherche en Informatique et Mathématiques, https://www.plafrim.fr/fr). The code used for the experiments, hawen, is developed by the first author and available at https://gitlab.inria.fr/ffaucher/hawen.

Appendix A Reciprocity-gap formula with Euler’s equations in seismic

In this appendix, we motivate the reciprocity-gap misfit, constructed from the variational formulation of Problem 1. In the context of seismic, we assume that the data (i.e. pressure and normal velocity) are acquired on a line $\Sigma$ which is slightly underneath the surface $\Gamma_{1}$. The domain is decomposed into above and below area from the line of receivers, such that $\Omega=\Omega^{+}\cup\Omega^{-}$, we illustrate in Figure 21.

Let us first rewrite Problem 1 introducing continuity conditions at $\Sigma$, using exponent + and - for the fields that are in $\Omega^{+}$ and $\Omega^{-}$ respectively, such that

[TABLE]

and similarly for the velocity. Problem 1 is equivalent to,

[TABLE]

We consider two couples $(p_{1},\,\boldsymbol{v}_{1})$ and $(p_{2},\,\boldsymbol{v}_{2})$ solutions to Problem Eq. 14. They are respectively associated to two sources, $f$ and $g$, and two different sets of physical parameters: $(\kappa_{1},\,\rho_{1})$ and $(\kappa_{2},\,\rho_{2})$ respectively. We take $p_{k}\in H^{1}(\Omega)$ ad $\boldsymbol{v}_{k}\in(H^{1}(\Omega))^{3}$, for $k=1$, $2$, where $H$ refers to the Hilbert space. For the derivation of the reciprocity-gap formula, we write the variational formulation on $\Omega^{-}$ only. The variational formulation of Eq. 14 for $(p_{1},\,\boldsymbol{v}_{1})$, using for the test functions $(p_{2},\,\boldsymbol{v}_{2})$ gives

[TABLE]

where we omit the - exponent in the fields $p$ and $\boldsymbol{v}$ for the sake of clarity. Furthermore we assume that the source $f$ is supported outside $\Omega^{-}$, e.g., by using delta-Dirac functions in some position $\boldsymbol{x}_{f}\in\Omega^{+}$, according to Figure 21.

We use an integration by part for both equations, and subtract the first one from the second one to get

[TABLE]

In the volume integral, we replace $(\nabla\cdot\boldsymbol{v}_{2})$ and $(\nabla p_{2})$ using the fact that $(p_{2},\,\boldsymbol{v}_{2})$ solves Eq. 14:

[TABLE]

Next, the integral on the boundary in Eq. 16 is decomposed between $\Gamma_{2}^{-}=\Gamma_{2}\cap\partial\Omega^{-}$ and $\Sigma$ such that

[TABLE]

On $\Gamma_{2}^{-}$, we first use Eq. 1a to replace $(\boldsymbol{v}_{1}\cdot\boldsymbol{\nu})$ and $(\boldsymbol{v}_{2}\cdot\boldsymbol{\nu})$ and then the absorbing boundary condition Eq. 1d:

[TABLE]

Eventually, injecting the new formulas for the volume and integral equations in Eq. 16, we obtain

[TABLE]

The right-hand side coincides with our choice misfit functional for full reciprocity-gap waveform inversion cf. Eq. 8 (where we take the squared of the expression to have a positive functional), and it equates to zero if and only if $\rho_{1}=\rho_{2}$ and $\kappa_{1}=\kappa_{2}$ (meaning that $c_{1}=c_{2}$) over the whole domain $\Omega^{-}$. We further refer to Alessandrini et al., (2018, 2019) for stability results. We see that the formula Eq. 20 does not involve the sources, because we have positioned them above of the receivers line (see Figure 21). In the case where sources are below, or if one wants to create a numerical acquisition with sources inside the domain, it results in an additional integral in Eq. 15, which then has to be included in the misfit.

Appendix B Gradient formulation using adjoint-state method

The gradient of the misfit functional is computed using the adjoint-state method, which comes from the work of Lions in optimal control, cf. Lions, (1971), with early applications in Chavent, (1974). The method is popular in seismic as it avoids the computation of the Jacobian matrix, thus requiring limited numerical effort; it is reviewed in Plessix, (2006).

For generality, we consider the misfit functional $\mathcal{J}_{\bullet}$, which is either $\mathcal{J}_{L2}$ of $\mathcal{J}_{\mathfrak{r}}$, and define the constrained minimization problem,

[TABLE]

where $\mathcal{A}$ is the linear wave operator corresponding with the Euler’s equations of Problem 1. We denote by $f_{i}$ the sources and use the notation $U^{(f_{i})}=\{\boldsymbol{v}^{(f_{i})},\,p^{(f_{i})}\}=\{v_{x}^{(f_{i})},\,v_{y}^{(f_{i})},\,v_{z}^{(f_{i})},\,p^{(f_{i})}\}$ for the solution associated with the volume source $F_{i}=\{0,\,0,\,0,\,f_{i}\}$, in accordance with Problem 1. For the sake of notation, we first consider a single source and drop the index $f_{i}$, the formulation with Lagrangian gives

[TABLE]

Here, $<.\,,\,.>_{\Omega}$ denotes the complex inner product in $L^{2}(\Omega)$ such that $<a,\,b>_{\Omega}=\int_{\Omega}\overline{a}\,b\,\mathrm{d}\Omega$, with $\overline{\phantom{a}}$ the complex conjugate.

The first step of the adjoint-state method is to take $\tilde{U}=U$ such that

[TABLE]

Then, the adjoint-state $\gamma$ is selected such that the derivative of the Lagrangian with respect to $U$ equates zero, i.e., $\gamma$ is solution to

[TABLE]

where ∗ denotes the adjoint (transposed of the complex conjugate). With this choice of adjoint-state, the gradient of the misfit functional (which coincides with the one of the Lagrangian from Eq. 23) is

[TABLE]

We further refer to the Appendix A of Barucq et al., (2019), Faucher, (2017) and Barucq et al., (2018) for more details on the complex-variable adjoint-state, and to Faucher and Scherzer, (2020) for the specificity with HDG discretization.

When we incorporate back the different sources, the adjoint problem Eq. 24 for the misfit functional $\mathcal{J}_{L2}$ is, for each source $f_{i}$ in the acquisition,

[TABLE]

Here we use $\boldsymbol{\nu}$ to indicate the normal direction (in the case of a flat surface, i.e. $x-y$ plane, $\boldsymbol{\nu}_{x}=\boldsymbol{\nu}_{y}=0$ and $\boldsymbol{\nu}_{z}=1$). The gradient using all sources is

[TABLE]

For the misfit functional $\mathcal{J}_{\mathfrak{r}}$ of Eq. 8, the adjoint-state solves, for each computational source $g_{j}$,

[TABLE]

with the scalar $\xi_{i,j}$ given by

[TABLE]

The gradient is

[TABLE]

We refer to Alessandrini et al., (2019) for more details on the adjoint-state method using Cauchy data. Therefore, we see that the adjoint-state for a single source in the computational acquisition for $\mathcal{J}_{\mathfrak{r}}$ encodes information from all measurements because of the reciprocity-gap (the sum over $i$ in Eq. 28) while it only takes the current source with $\mathcal{J}_{L2}$, see Eq. 26.