Multi-parameter inverse scattering

# Efficient scattering approach to seismic full-waveform inversion in
anisotropic elastic media with variable density

Kui Xiang111footnotemark: 122footnotemark: 2 Morten Jakobsen1,2 Ujjwal
Shekhar2 Kjersti Solberg Eikrem3 and Geir Nævdal3

###### Abstract

This paper introduces a novel matrix-free approach for full waveform inversion
in anisotropic elastic media, incorporating density variation through the
utilization of the distorted Born iterative method. This study aims to
overcome the computational and storage challenges associated with the
conventional matrix-based distorted Born iterative inversion method while
accurately capturing the subsurface’s anisotropic properties and density
variations. An elastic integral equation is utilized to account for the
anisotropic nature of elastic wave propagation, enabling more precise modeling
of subsurface complexities. This integral equation is efficiently solved by a
fast Fourier transform accelerated Krylov subspace method. Leveraging the
integral equation with the distorted Born approximation, a linear relationship
between the scattered wavefield and the model parameter perturbation is
formulated for an integrated inversion scheme. To address the inherent ill-
posedness of each linear inversion step, we formulate the normal equation with
a regularization term. This is achieved by minimizing an objective function
using the generalized Tikhonov method. Therefore, we can find an adequate
solution for the inverse scattering problem by solving the normal equation.
Following the physical interpretation of Green’s function, the Fréchet and
adjoint operators within the normal equation can be employed in a matrix-free
manner, allowing for significant improvement of the computational efficiency
and memory demand without compromising accuracy. The proposed matrix-free full
waveform inversion framework is thoroughly validated through extensive
numerical experiments on synthetic datasets, showcasing its ability to
reconstruct complex anisotropic structures and accurately recover stiffness
parameters and density. This study contributes to advancing seismic imaging
techniques and holds promise for applications in various domains, including
hydrocarbon exploration, geological hazard assessment, and geotechnical
engineering. Integrating stiffness parameters and density variation within the
matrix-free full waveform inversion framework paves the way for a more
comprehensive understanding of subsurface properties.

##  1 Introduction

Seismic full waveform inversion (FWI) is one of the most powerful tools for
inferring the properties and structure of the subsurface by iteratively
minimizing the residuals between calculated and observed data (Tarantola,,
1984, 2005; Virieux and Operto,, 2009). It can provide a high-resolution
reconstructed model by using the full data content, including diving waves,
reflections, and multiples. Depending on the chosen forward-modeling approach,
FWI can be carried out either in the time domain (Tarantola,, 1984; Mora,,
1988; Bunks et al.,, 1995) or in the frequency domain (Pratt,, 1990; Zhou et
al.,, 1995; Liao and McMechan,, 1996). From the mathematical physics
perspective, FWI can be categorized as an inverse scattering problem. Weglein
et al., (2003) provides a comprehensive overview of the essential logical
steps and the foundational mathematical-physics principles behind the inverse
scattering series used in seismic exploration. The inverse scattering problem
has been extensively explored in physics, engineering, and applied sciences
(Pike and Sabatier,, 2001). There exist several highly developed methods for
solving linear or nonlinear inverse scattering problems (Weglein et al.,,
1981; Stolt and Jacobs,, 1980; Van Den Berg and Kleinman,, 1997; Abubakar and
van den Berg,, 2004; Innanen and Lira,, 2010; Jakobsen and Ursin,, 2012;
Haffinger et al.,, 2013; Jakobsen and Wu,, 2016; Osnabrugge et al.,, 2016;
Malovichko et al.,, 2018; Song et al.,, 2020; van den Berg,, 2021). Therefore,
adapting these methods for application in seismic FWI is promising.

One of these methods is the so-called distorted Born iterative (DBI) method.
It is based on the integral equation method and was initially proposed for
electromagnetic inverse scattering (Wang and Chew,, 1989; Chew and Wang,,
1990). Although the implementation and performance of these methods are
different, Remis and Van den Berg, (2000) and Oristaglio and Blok, (2012)
reported that the DBI method is equivalent to the Gauss–Newton method in
optimization. The concept underlying the DBI method involves approximating the
nonlinear inverse scattering problem with a sequence of linear problems. Then,
each linear problem is solved to update the model incrementally. The DBI
method has several notable characteristics: (1) it requires discretization
solely for the anomalous region, enhancing efficiency for time-lapse inversion
(Malovichko et al.,, 2017; Huang et al.,, 2019); (2) it can provide the
Fréchet derivative (sensitivity kernel) directly via the use of Green’s
function, making it suitable for computing the Hessian information, which is
crucial for reducing crosstalk effects in multi-parameter inversion (Jakobsen
et al.,, 2020); and (3) when Green’s function pertains to homogeneous
background media, its computation can be expedited using the fast Fourier
transform (FFT) technique (Beylkin et al.,, 2009). In the geophysical
community, Jakobsen and Ursin, (2015) introduced a T-matrix variant of the DBI
method for solving the nonlinear problems in seismic FWI. Huang et al., (2019)
verified that the T-matrix based DBI method is naturally target-oriented and
successfully implemented it to the time-lapse seismic waveform inversion of a
target region. Eikrem et al., (2019) proposed a Bayesian version of the DBI
method for FWI, in which they also estimate the uncertainty surrounding the
maximum a posteriori solution using the iterated extended Kalman filter. The
main disadvantages of the DBI method are its large computational cost and
memory requirements. However, these disadvantages can be circumvented by
implementing the Fréchet derivative and its adjoint in a fast and efficient
manner (Xiang et al.,, 2023), in which the matrix-vector multiplications are
achieved by operation of the Green’s function with FFT.

The Earth’s subsurface is defined by an intricate combination of parameters,
including seismic velocities (P-wave and S-wave), density, attenuation, and
anisotropy. Each of these parameters affects seismic wave propagation in
unique ways. With the development of geophysical exploration technology,
simple reservoirs with straightforward geological features become increasingly
rare. Instead, we face intricate oil and gas reservoirs characterized by
complex subsurface structures, including intricate tectonics, stratigraphy,
and lithology. The imaging requirements of these complex reservoirs are such
that we cannot ignore the impact of multiple parameters on the inversion
results. By updating multiple parameters concurrently, multiparameter FWI can
provide more accurate and high-resolution subsurface models (Operto et al.,,
2013). In light of significant development in computer technology over recent
decades, multiparameter FWI has been widely studied in numerous publications
(Brossier et al.,, 2009; Köhn et al.,, 2012; Prieux et al.,, 2013; Lin and
Huang,, 2014; Wang and Cheng,, 2017). Some publications have shown the results
of multiparameter FWI in the anisotropic situation (Burridge et al.,, 1998;
Lee et al.,, 2010; Warner et al.,, 2013; Kamath and Tsvankin,, 2016; Oh and
Alkhalifah,, 2016; Pan et al.,, 2016; Alkhalifah et al.,, 2016; He and
Plessix,, 2017; Rusmanugroho et al.,, 2017; Yang et al.,, 2019).

In the seismic full-waveform inversion algorithm utilizing the DBI method,
Jakobsen et al., (2020) generalized the DBI T-matrix method so that it can be
used in conjunction with the elastodynamic equation for the anisotropic
elastic medium. Huang et al., (2020) extended the Bayesian framework of Eikrem
et al., (2019), which incorporated the DBI method and the iterated extended
Kalman filter, to anisotropic elastic media. Although these implementations of
the DBI method in seismic multi-parameter inversion obtained good inversion
results, the computational cost and memory requirements remain substantial due
to the matrix formulation employed. Jakobsen et al., (2023) generalized the
matrix-free DBI method from mono-parameter (velocity) to multi-parameter
(velocity and density) FWI and extended the implementation of the matrix-free
DBI method to seismic inversion and medical imaging. In this paper, we further
extend the matrix-free DBI method to FWI in elastic anisotropic media.
Compared to previous studies, the approach proposed in this paper has several
innovative aspects: (1) we only use a single displacement-related integral
equation instead of two coupled integral equations for the particle
displacement and the strain fields; (2) this is matrix-free for FWI in elastic
anisotropic media; (3) the reconstruction of density is included in both
theory and numerical examples; (4) we employ finite differences to compute the
derivative of the Green’s function rather than the previous analytical
formulation for the derivative of the Green’s function, which offers a reduced
computational cost compared to the analytical formulation; and (5) we extend
the application of the DBI method for FWI in elastic anisotropic media from 2D
to 3D.

This paper is organized as follows. We first provide the integral equation of
the forward problem for elastic anisotropic media, which considers the 21
elastic parameters and the mass density. We then describe a fast forward
solver, in conjunction with the Krylov subspace method and the FFT, for
solving the integral equation, which will be used for simulating the synthetic
data. Next, we derive the mathematical formulation of the elastic DBI method
and provide the matrix-free expressions of the Fréchet derivative and its
adjoint operators. The normal equation required for solving the inverse
problem is also described. Finally, we use 2D and 3D numerical examples to
test the performance of our method for elastic FWI in the transversely
isotropic (VTI) media with a vertical symmetry axis.

## Methodology

###  1.1 Forward modeling

The propagation of seismic waves within a heterogeneous anisotropic elastic
medium occupying a bounded domain $D\subset\mathbf{R}^{d}$, $d=2$, 3, is
governed by the elastodynamic wave equation (Červenỳ,, 2001; Jakobsen et al.,,
2020):

$\omega^{2}\rho(\mathbf{x})\mathbf{u}(\mathbf{x},\omega)+\nabla\cdot\left[\mathbf{C}(\mathbf{x}):\nabla\mathbf{u}(\mathbf{x},\omega)\right]=-\mathbf{S}(\mathbf{x},\omega),$
(1)

where $\rho\ $ is the mass density, $\mathbf{C}$ is the elastic stiffness
tensor, $\omega$ is the angular frequency, $\mathbf{u}\subset\mathbf{R}^{d}$,
is the displacement, and $\mathbf{S}\subset\mathbf{R}^{d}$ is the source. The
double dot product of the fourth rank tensor $\mathbf{C}$ and the second rank
tensor $\nabla\mathbf{u}$ is defined by the summation over the two repeated
indices (Auld,, 1973). For clarity in notation, we will omit the
$\omega$-dependency within the displacement $\mathbf{u}$, the source
$\mathbf{S}$ and the following Green’s function. The displacement field
$\mathbf{u}$ can be expressed as an integral over all space (Morse and
Feshbach,, 1954; Gubernatis et al.,, 1977; Jakobsen et al.,, 2020), often
referred to as the source representation, given by:

$\mathbf{u}(\mathbf{x})=\int
d\mathbf{x}^{\prime}\mathbf{G}\left(\mathbf{x},\mathbf{x}^{\prime}\right)\mathbf{S}\left(\mathbf{x}^{\prime}\right),$
(2)

where $\mathbf{G}\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ is the Green’s
function, defined by

$\omega^{2}\rho(\mathbf{x})\mathbf{G}(\mathbf{x},\mathbf{x}^{\prime})+\nabla\cdot\left[\mathbf{C}(\mathbf{x}):\nabla\mathbf{G}(\mathbf{x},\mathbf{x}^{\prime})\right]=-\delta(\mathbf{x}-\mathbf{x}^{\prime}),$
(3)

where the delta function $\delta(\mathbf{x}-\mathbf{x}^{\prime})$ represents a
unit point source.

Consider a homogeneous anisotropic elastic model described by stiffness
parameters $\mathbf{C}^{(0)}$ and density $\rho^{(0)}$ and introduce the
perturbations of elastic parameters and density, $\delta\mathbf{C}^{(0)}$ and
$\delta\rho^{(0)}$, as follows:

$\mathbf{C}(\mathbf{x})=\mathbf{C}^{(0)}+\delta\mathbf{C}^{(0)}(\mathbf{x}),\qquad\rho(\mathbf{x})=\rho^{(0)}+\delta\rho^{(0)}(\mathbf{x}).$
(4)

Inserting (4) into (1), we obtain

$\omega^{2}\rho^{(0)}\mathbf{u}(\mathbf{x})+\nabla\cdot\left[\mathbf{C}^{(0)}:\nabla\mathbf{u}(\mathbf{x})\right]=-\mathbf{S}(\mathbf{x})-\omega^{2}\delta\rho^{(0)}\left(\mathbf{x}\right)\mathbf{u}(\mathbf{x})-\nabla\cdot\left[\delta\mathbf{C}^{(0)}\left(\mathbf{x}\right):\nabla\mathbf{u}(\mathbf{x})\right].$
(5)

If we treat the right-hand side of (5) as a new source term, then
$\mathbf{u}(\mathbf{x})$ can be interpreted as the displacement in the
homogeneous background model due to the new source. Using the volume integral
(2), we can represent equation (5) in the integral form:

$\displaystyle\mathbf{u}\left(\mathbf{x}\right)=$
$\displaystyle\mathbf{u}^{(0)}\left(\mathbf{x}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\rho^{(0)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}\left(\mathbf{x}^{\prime}\right)$
(6) $\displaystyle+\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\nabla_{\mathbf{x}^{\prime}}\cdot\left[\delta\mathbf{C}^{(0)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right],$

where

$\mathbf{u}^{(0)}(\mathbf{x})=\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\mathbf{S}\left(\mathbf{x}^{\prime}\right)$
(7)

is the displacement in the homogeneous background medium caused by
$\mathbf{S}\left(\mathbf{x}^{\prime}\right)$ and
$\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)$ is the Green’s
function for the homogeneous medium. The Green’s function for a homogeneous
medium can be derived from equation (3) by substituting the actual medium with
the homogeneous medium. The analytical expressions of 2D and 3D elastodynamic
Green’s function for homogeneous media are given in Appendix LABEL:op_Green.
Since Green’s function for the homogeneous medium is translation invariant, we
present it as a function of the vector difference
$\mathbf{x}-\mathbf{x}^{\prime}$. More details about the Green’s function can
be found in Arfken and Weber, (1999) and Červenỳ, (2001).

Equation (6) involves the spatial derivative on perturbations
$\delta\mathbf{C}^{(0)}$, which is not convenient for subsequent inversion
tasks. Therefore, we need to rewrite equation (6) in a more suitable form
(Červenỳ,, 2001). Using

$\displaystyle\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\nabla_{\mathbf{x}^{\prime}}\cdot\left[\delta\mathbf{C}^{(0)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right]=$
$\displaystyle\nabla_{\mathbf{x}^{\prime}}\cdot\left[\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\mathbf{C}^{(0)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right]$
(8)
$\displaystyle-\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\cdot\left[\delta\mathbf{C}^{(0)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right],$

inserting this into equation (6), we obtain

$\displaystyle\mathbf{u}\left(\mathbf{x}\right)=$
$\displaystyle\mathbf{u}^{(0)}\left(\mathbf{x}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\rho^{(0)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}\left(\mathbf{x}^{\prime}\right)$
(9) $\displaystyle-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\cdot\left[\delta\mathbf{C}^{(0)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right].$

Note that the volume integral over the first term on the right hand of
equation (8) has been transformed into a surface integral which is identical
to zero under the assumption that the displacement field $\mathbf{u}$
approaches zero at infinity.

Equation (9) is the desired integral equation for our computational approach.
Shekhar et al., (2023) proposed utilizing this single integral equation (9) to
solve for the displacement field instead of relying on two coupled integral
equations as presented in Jakobsen et al., (2020). In the work presented by
Shekhar et al., (2023), the derivative of the Green’s function
$\nabla\mathbf{G}^{(0)}$ is obtained analytically, while the derivative of the
displacement $\nabla\mathbf{u}$ is computed using the finite difference
method. In this paper, the derivative of Green’s function and the derivative
of the displacement are computed using the finite difference method. Once the
perturbations $\delta\mathbf{C}^{(0)}$ and $\delta\rho^{(0)}$, the incident
wavefield $\mathbf{u}^{(0)}$, and the Green’s function $\mathbf{G}^{(0)}$ are
known, equation (9) enables us to calculate the displacement $\mathbf{u}$. To
solve $\mathbf{u}$ from (9), we first define an operator as

$\displaystyle\mathcal{A}\left(\cdot\right)=$
$\displaystyle\mathcal{I}\left(\cdot\right)-\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\rho^{(0)}\left(\mathbf{x}^{\prime}\right)\left(\cdot\right)$
(10) $\displaystyle+\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\cdot\left[\delta\mathbf{C}^{(0)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\left(\cdot\right)\right],$

where $\mathcal{I}$ is the identity operator. Substituting (10) into (9), we
have

$\mathcal{A}\left\\{\mathbf{u}\left(\mathbf{x}\right)\right\\}=\mathbf{u}^{(0)}\left(\mathbf{x}\right).$
(11)

Equation (11) is a linear equation in the form of
$\mathbf{A}\mathbf{x}=\mathbf{b}$, which can be iteratively solved using
Krylov subspace methods, such as the conjugate gradient method, BiCGSTAB
(biconjugate gradient stabilized) method, GMRES (generalized minimum residual)
method, among others. In operator (10), the Green’s function and its
derivative, when multiplied by their subsequent vectors, can be interpreted as
circular convolutions. These circular convolutions can be computed efficiently
by using the Fast Fourier Transform. The computations are written symbolically
as

$\mathbf{G}^{(0)}\left\\{\mathbf{v}\right\\}=\mathrm{FFT}^{-1}\\{\mathrm{FFT}\\{\mathbf{G}^{(0)}\\}\odot\mathrm{FFT}\\{\mathbf{v}\\}\\},$
(12)

where $\mathbf{v}$ can be any vector and $\odot$ denotes element-wise
multiplication.

###  1.2 The distorted Born iterative method

Let us consider a known heterogeneous reference model described by elastic
parameters $\mathbf{C}^{(b)}$ and density $\rho^{(b)}$. The differences in the
elastic parameters and density between the true and the reference model are
defined as

$\delta\rho\left(\mathbf{x}\right)=\rho\left(\mathbf{x}\right)-\rho^{(b)}\left(\mathbf{x}\right),\qquad\delta\mathbf{C}\left(\mathbf{x}\right)=\mathbf{C}\left(\mathbf{x}\right)-\mathbf{C}^{(b)}\left(\mathbf{x}\right).$
(13)

The discrepancy between the true displacement field $\mathbf{u}$ in the true
model and the reference displacement field $\mathbf{u}^{(b)}$ in the reference
model is defined as

$\delta\mathbf{u}\left(\mathbf{x}\right)=\mathbf{u}\left(\mathbf{x}\right)-\mathbf{u}^{(b)}\left(\mathbf{x}\right).$
(14)

By replacing the homogeneous background model in equation (9) with the
heterogeneous reference model, we have

$\displaystyle\mathbf{u}\left(\mathbf{x}\right)=$
$\displaystyle\mathbf{u}^{(b)}\left(\mathbf{x}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{x}^{\prime}\right)\delta\rho\left(\mathbf{x}^{\prime}\right)\mathbf{u}\left(\mathbf{x}^{\prime}\right)$
(15) $\displaystyle-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{x}^{\prime}\right)\cdot\left[\delta\mathbf{C}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right],$

where $\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ is the
Green’s function for a heterogeneous background medium that can be numerically
solved from equation (3). The heterogeneous Green’s function
$\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ can also be
solved from the Dyson equations in Jakobsen and Ursin, (2015). However, in
this paper, we will employ operations only involving the homogeneous Green’s
function as an alternative to the direct use of the heterogeneous Green’s
function, thereby circumventing the need for storage and updates of the
Green’s function associated with the model update. When performing inversion,
the data residual of the displacement field can only be observed at the
receiver position $\mathbf{r}$. By combining this information with equation
(15) and (14), we obtain

$\delta\mathbf{u}\left(\mathbf{r}\right)=\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\delta\rho\left(\mathbf{x}^{\prime}\right)\mathbf{u}\left(\mathbf{x}^{\prime}\right)-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\cdot\left[\delta\mathbf{C}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}\left(\mathbf{x}^{\prime}\right)\right].$
(16)

In equation (16), it is evident that to solve for $\delta\rho$ and
$\delta\mathbf{C}$, we require knowledge of the displacement
$\mathbf{u}\left(\mathbf{x}^{\prime}\right)$ within the unknown true model,
which is unavailable to us. However, we know the heterogeneous background
model and its corresponding displacement field. Therefore, we assume that the
known heterogeneous background model is close to the unknown true model, which
results in the reference displacement being close to the true displacement.
This assumption allows equation (16) to be approximated as

$\delta\mathbf{u}\left(\mathbf{r}\right)\approx\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\delta\rho\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(b)}\left(\mathbf{x}^{\prime}\right)-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\cdot\left[\delta\mathbf{C}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(b)}\left(\mathbf{x}^{\prime}\right)\right].$
(17)

In the distorted Born iterative method method, the displacement residual
$\delta\mathbf{u}$, the background Green’s function $\mathbf{G}^{(b)}$, and
the background displacement $\mathbf{u}^{(b)}$ are computed first based on a
known background model. After that, we can solve (17) for $\delta\mathbf{C}$
and $\delta\rho$, and then iteratively update the known background model.

###  1.3 Fréchet derivative and adjoint operators

To further progress, we must decompose the stiffness tensor perturbations into
a spatially invariant tensor structure part and a spatially variable scalar
function for each independent model parameter. Thus, we decompose
$\delta\mathbf{C}\left(\mathbf{x}\right)$ as

$\delta\mathbf{C}\left(\mathbf{x}\right)=\sum_{p=1}^{21}\mathbf{B}^{(p)}\delta
m^{(p)}\left(\mathbf{x}\right),$ (18)

where $\delta m^{(p)},\quad p=1,2,...,21$ are the model perturbations of the
21 independent elastic parameters, and $\mathbf{B}$ is the tensor field
structure related to the position of the elastic parameter in the stiffness
tensor (Jakobsen et al.,, 2020). To ensure consistency with the representation
of $\delta m^{(p)}$ in terms of elastic parameters, we introduce the notation
$\delta m^{(0)}=\delta\rho$ for the density. Inserting (18) into (17), we
obtain

$\displaystyle\delta\mathbf{u}\left(\mathbf{r}\right)\approx$
$\displaystyle\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\delta
m^{(0)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(b)}\left(\mathbf{x}^{\prime}\right)$
(19) $\displaystyle-\sum_{p=1}^{21}\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\cdot\left[\mathbf{B}^{(p)}\delta
m^{(p)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(b)}\left(\mathbf{x}^{\prime}\right)\right].$

From equation (19), the Fréchet derivative operators corresponding to
different parameters can be defined as

$\displaystyle\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)=\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\mathbf{u}^{(b)}\left(\mathbf{x}^{\prime}\right)\delta
m^{(0)}\left(\mathbf{x}^{\prime}\right),$ (20)
$\displaystyle\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right)=-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\cdot\left[\mathbf{B}^{(p)}:\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(b)}\left(\mathbf{x}^{\prime}\right)\delta
m^{(p)}\left(\mathbf{x}^{\prime}\right)\right].$

Here, $\mathcal{F}^{(0)}$ is related to density perturbation and
$\mathcal{F}^{(p)}$ is related to elastic parameters perturbation. Inserting
(20) into (19) yields

$\delta\mathbf{u}\left(\mathbf{r}\right)\approx\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)+\sum_{p=1}^{21}\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right).$
(21)

Equation (19) can be further rewritten as

$\delta\mathbf{u}\approx\mathcal{F}\delta\mathbf{m},$ (22)

where
$\mathcal{F}=[\mathcal{F}^{(0)},\mathcal{F}^{(1)},...,\mathcal{F}^{(21)}]$ and
$\delta\mathbf{m}=[\delta\mathbf{m}^{(0)},\delta\mathbf{m}^{(1)},...,\delta\mathbf{m}^{(21)}]^{T}$.

In most seismic applications, equation (22) is ill-posed due to the Fréchet
derivative operator being many-to-one rather than one-to-one. As a result, the
solution of equation (22) is not unique. To obtain an appropriate solution to
(22), the generalized Tikhonov method is employed to ensure the stability of
the solution. In the generalized Tikhonov method, a regularized solution to
(22) is found by minimizing the objective function, as described in (Menke,,
2012):

$\mathcal{E}\left(\delta\mathbf{m}\right)=\left\|\delta\mathbf{u}-\mathcal{F}\delta\mathbf{m}\right\|_{2}^{2}+\lambda\left\|\delta\mathbf{m}\right\|_{2}^{2},$
(23)

in which $\left\|\cdot\right\|_{2}^{2}$ is the $L_{2}$ norm and $\lambda$ is
the regularization parameter. This paper determines the regularization
parameter using a self-adaptive cooling scheme (Jakobsen and Ursin,, 2015). An
initial value for the regularization parameter is set at the beginning of the
iteration process and is progressively reduced as the iterations proceed. To
find the minimizer of $\delta\mathbf{m}$ of (23), the following normal
equation is utilized:

$(\mathcal{H}+\lambda\mathcal{I})\delta\mathbf{m}=\mathcal{F}^{\dagger}\delta\mathbf{u},$
(24)

where $\mathcal{H}=\mathcal{F}^{\dagger}\mathcal{F}$ is the corresponding
approximate Hessian operator, $\mathcal{I}$ is the identity operator, and
$\mathcal{F}^{\dagger}$ is the adjoint of the Fréchet operator. By utilizing
the definition of adjoint (Tarantola,, 2005; Claerbout and Fomel,, 2012), we
have

$\displaystyle\left[(\mathcal{F}^{(0)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[\omega^{2}\int
d\mathbf{r}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{r}\right)\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*},$
(25)
$\displaystyle\left[(\mathcal{F}^{(p)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[-\nabla_{\mathbf{x}}\int
d\mathbf{r}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{r}\right)\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\cdot\mathbf{B}^{(p)}:\nabla_{\mathbf{x}}\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*},$

where $\left[\cdot\right]^{*}$ denotes complex conjugation. More details on
the derivation of the adjoint operator can be found in Appendix LABEL:op_ad.

###  1.4 Matrix-free formulations of the Fréchet and adjoint operator

In the normal equation (24), the model perturbation $\delta\mathbf{m}$ can be
solved by inverting $\mathcal{H}+\lambda\mathcal{I}$. However, direct
construct and inversion of $\mathcal{H}+\lambda\mathcal{I}$ become
computationally expensive and memory-intensive when dealing with realistic-
scale problems. Instead, we solve the normal equation (24) via a Krylov
subspace method in conjunction with the matrix-free implementation of the
Fréchet and adjoint operator. By considering $\mathbf{u}^{(b)}\delta m^{(0)}$
and $\mathbf{B}^{(p)}:\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(b)}\delta
m^{(p)}$ in equation (20) as two virtual sources, denoted as $\mathbf{Y}$ and
$\mathbf{M}$, we obtain

$\displaystyle\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)$
$\displaystyle=\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\mathbf{Y}\left(\mathbf{x}^{\prime}\right),$
(26)
$\displaystyle\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right)$
$\displaystyle=-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}^{\prime}\right)\cdot\mathbf{M}\left(\mathbf{x}^{\prime}\right).$

Here, $\mathbf{Y}$ is a vector source, while $\mathbf{M}$ is a second rank
moment tensor source. Drawing inspiration from the works of Hesford and Chew,
(2006, 2010) and Jakobsen et al., (2023), the physical interpretations of the
two sub-equations in (26) can be interpreted as the observed displacement
fields from the background media due to sources $\mathbf{Y}$ and $\mathbf{M}$.
Consequently, equation (26) can be expressed as:

$\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)=\omega^{2}\mathbf{u}^{(Y)}\left(\mathbf{r}\right),\qquad\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right)=-\mathbf{u}^{(M)}\left(\mathbf{r}\right),$
(27)

where the observed displacement fields
$\mathbf{u}^{(Y)}\left(\mathbf{r}\right)$ and
$\mathbf{u}^{(M)}\left(\mathbf{r}\right)$ can be solved from

$\displaystyle\begin{aligned}
\mathbf{u}^{(Y)}\left(\mathbf{r}\right)=&\mathbf{u}^{(Y,0)}\left(\mathbf{r}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{r}-\mathbf{x}^{\prime}\right)\delta\rho^{(b)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(Y)}\left(\mathbf{x}^{\prime}\right)\\\
&-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{r}-\mathbf{x}^{\prime}\right)\cdot\delta\mathbf{C}^{(b)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(Y)}\left(\mathbf{x}^{\prime}\right),\end{aligned}$
(28) $\displaystyle\begin{aligned}
\mathbf{u}^{(M)}\left(\mathbf{r}\right)=&\mathbf{u}^{(M,0)}\left(\mathbf{r}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{r}-\mathbf{x}^{\prime}\right)\delta\rho^{(b)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(M)}\left(\mathbf{x}^{\prime}\right)\\\
&-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{r}-\mathbf{x}^{\prime}\right)\cdot\delta\mathbf{C}^{(b)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(M)}\left(\mathbf{x}^{\prime}\right),\end{aligned}$

where

$\delta\rho^{(b)}=\rho^{(b)}-\rho^{(0)},\qquad\delta\mathbf{C}^{(b)}=\mathbf{C}^{(b)}-\mathbf{C}^{(0)}.$
(29)

To solve equation (28) involves knowledge of the displacement fields
$\mathbf{u}^{(Y)}\left(\mathbf{x}^{\prime}\right)$ and
$\mathbf{u}^{(M)}\left(\mathbf{x}^{\prime}\right)$ within the heterogeneous
background media, which we need to first compute by solving the forward-
scattering equation

$\displaystyle\begin{aligned}
\mathbf{u}^{(Y)}\left(\mathbf{x}\right)=&\mathbf{u}^{(Y,0)}\left(\mathbf{x}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\rho^{(b)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(Y)}\left(\mathbf{x}^{\prime}\right)\\\
&-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\cdot\delta\mathbf{C}^{(b)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(Y)}\left(\mathbf{x}^{\prime}\right),\end{aligned}$
(30) $\displaystyle\begin{aligned}
\mathbf{u}^{(M)}\left(\mathbf{x}\right)=&\mathbf{u}^{(M,0)}\left(\mathbf{x}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\rho^{(b)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(M)}\left(\mathbf{x}^{\prime}\right)\\\
&-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\cdot\delta\mathbf{C}^{(b)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(M)}\left(\mathbf{x}^{\prime}\right),\end{aligned}$

where

$\displaystyle\mathbf{u}^{(Y,0)}\left(\mathbf{x}\right)=\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\mathbf{Y}\left(\mathbf{x}^{\prime}\right),$
(31) $\displaystyle\mathbf{u}^{(M,0)}\left(\mathbf{x}\right)=\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\mathbf{M}\left(\mathbf{x}^{\prime}\right)$

are the reference displacement fields.

In the same way, following the physical interpretations of the sub-equations
in (25), the adjoint operators can be rewritten as

$\displaystyle\left[(\mathcal{F}^{(0)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[\omega^{2}\mathbf{u}^{(a)}\left(\mathbf{x}\right)\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*},$
(32)
$\displaystyle\left[(\mathcal{F}^{(p)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[-\nabla_{\mathbf{x}}\mathbf{u}^{(a)}\left(\mathbf{x}\right)\cdot\mathbf{B}^{(p)}:\nabla_{\mathbf{x}}\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*},$

where

$\mathbf{u}^{(a)}\left(\mathbf{x}\right)=\int
d\mathbf{r}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{r}\right)\delta\mathbf{u}^{*}\left(\mathbf{r}\right).$
(33)

Here, $\mathbf{u}^{(a)}$ is the back-propagating displacement field within the
heterogeneous background media from source $\delta\mathbf{u}^{*}$, which can
be solved from

$\displaystyle\begin{aligned}
\mathbf{u}^{(a)}\left(\mathbf{x}\right)=&\mathbf{u}^{(a,0)}\left(\mathbf{x}\right)+\omega^{2}\int
d\mathbf{x}^{\prime}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\delta\rho^{(b)}\left(\mathbf{x}^{\prime}\right)\mathbf{u}^{(a)}\left(\mathbf{x}^{\prime}\right)\\\
&-\int
d\mathbf{x}^{\prime}\nabla_{\mathbf{x}^{\prime}}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{x}^{\prime}\right)\cdot\delta\mathbf{C}^{(b)}\left(\mathbf{x}^{\prime}\right):\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(a)}\left(\mathbf{x}^{\prime}\right),\end{aligned}$
(34)

where

$\displaystyle\mathbf{u}^{(a,0)}\left(\mathbf{x}\right)=\int
d\mathbf{r}\mathbf{G}^{(0)}\left(\mathbf{x}-\mathbf{r}\right)\delta\mathbf{u}^{*}\left(\mathbf{r}\right).$
(35)

Solving equations (30) and (34) are forward modeling problems that can be
effectively solved utilizing the methods mentioned in the preceding section.
The integrals involving Green’s function and its derivative in all the above
equations can be treated as convolutions and can be efficiently and accurately
computed using the fast Fourier transform. This eliminates the need to store
the matrix of the homogeneous Green’s function and compute the Green’s
function for the heterogeneous background medium. The current formulations of
the Fréchet derivative operator (27) and its adjoint operator (32) only
involve vectors: $\mathbf{u}^{(Y)}$, $\mathbf{u}^{(M)}$, $\mathbf{u}^{(a)}$,
$\mathbf{u}^{(b)}$, $\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(a)}$, and
$\nabla_{\mathbf{x}^{\prime}}\mathbf{u}^{(b)}$, making our method matrix-free.

###  1.5 Abbreviated subscript notation for implementation

The components of the elastic stiffness tensor exhibit symmetries
($C_{ijkl}=C_{jikl}=C_{ijlk}=C_{jilk}$); thus, employing the abbreviated
subscript notation (Auld,, 1973), the tensor can be compactly represented as a
6$\times$6 matrix:

$\mathbf{C}=\left[\begin{array}[]{cccccc}C_{11}&C_{12}&C_{13}&C_{14}&C_{15}&C_{16}\\\
C_{12}&C_{22}&C_{23}&C_{24}&C_{25}&C_{26}\\\
C_{13}&C_{23}&C_{33}&C_{34}&C_{35}&C_{36}\\\
C_{14}&C_{24}&C_{34}&C_{44}&C_{45}&C_{46}\\\
C_{15}&C_{25}&C_{35}&C_{45}&C_{55}&C_{56}\\\
C_{16}&C_{26}&C_{36}&C_{46}&C_{56}&C_{66}\end{array}\right].$ (36)

In the same way, by using the abbreviated subscript notation, the tensors
representing the derivatives of the Green’s function $\nabla\mathbf{G}^{(0)}$
and the displacement gradient $\nabla\mathbf{u}$ can be represented as:

$\nabla_{\mathbf{x}}\mathbf{G}^{(0)}=\left[\begin{array}[]{cccccc}\frac{\partial
G^{(0)}_{11}}{\partial x_{1}}&\frac{\partial G^{(0)}_{12}}{\partial
x_{2}}&\frac{\partial G^{(0)}_{13}}{\partial x_{3}}&\frac{\partial
G^{(0)}_{12}}{\partial x_{3}}+\frac{\partial G^{(0)}_{13}}{\partial
x_{2}}&\frac{\partial G^{(0)}_{11}}{\partial x_{3}}+\frac{\partial
G^{(0)}_{13}}{\partial x_{1}}&\frac{\partial G^{(0)}_{11}}{\partial
x_{2}}+\frac{\partial G^{(0)}_{12}}{\partial x_{1}}\\\ \frac{\partial
G^{(0)}_{21}}{\partial x_{1}}&\frac{\partial G^{(0)}_{22}}{\partial
x_{2}}&\frac{\partial G^{(0)}_{23}}{\partial x_{3}}&\frac{\partial
G^{(0)}_{22}}{\partial x_{3}}+\frac{\partial G^{(0)}_{23}}{\partial
x_{2}}&\frac{\partial G^{(0)}_{21}}{\partial x_{3}}+\frac{\partial
G^{(0)}_{23}}{\partial x_{1}}&\frac{\partial G^{(0)}_{21}}{\partial
x_{2}}+\frac{\partial G^{(0)}_{22}}{\partial x_{1}}\\\ \frac{\partial
G^{(0)}_{31}}{\partial x_{1}}&\frac{\partial G^{(0)}_{32}}{\partial
x_{2}}&\frac{\partial G^{(0)}_{33}}{\partial x_{3}}&\frac{\partial
G^{(0)}_{32}}{\partial x_{3}}+\frac{\partial G^{(0)}_{33}}{\partial
x_{2}}&\frac{\partial G^{(0)}_{31}}{\partial x_{3}}+\frac{\partial
G^{(0)}_{33}}{\partial x_{1}}&\frac{\partial G^{(0)}_{31}}{\partial
x_{2}}+\frac{\partial G^{(0)}_{32}}{\partial x_{1}}\\\ \end{array}\right],$
(37)

and

$\nabla_{\mathbf{x}}\mathbf{u}=\left[\begin{array}[]{cccccc}\frac{\partial
u_{1}}{\partial x_{1}}&\frac{\partial u_{2}}{\partial x_{2}}&\frac{\partial
u_{3}}{\partial x_{3}}&\frac{\partial u_{2}}{\partial x_{3}}+\frac{\partial
u_{3}}{\partial x_{2}}&\frac{\partial u_{1}}{\partial x_{3}}+\frac{\partial
u_{3}}{\partial x_{1}}&\frac{\partial u_{1}}{\partial x_{2}}+\frac{\partial
u_{2}}{\partial x_{1}}\end{array}\right]^{T},$ (38)

where $T$ denotes the transpose. The matrix-free distorted Born iterative
method proposed in this study has the theoretical capability to invert all 21
elastic parameters and density. However, in practical applications, inverting
such a vast number of parameters simultaneously is not feasible due to the
significant computation and storage resources required. Therefore, we restrict
our numerical tests to the transversely isotropic (VTI) media with a vertical
symmetry axis with variable density. The VTI media are generally described by
five independent elastic parameters: $C_{11}$, $C_{33}$, $C_{55}$, $C_{66}$,
and $C_{13}$ (Carcione,, 2014). The simplified elasticity tensor corresponding
to VTI media is presented as follows:

$\mathbf{C}^{VTI}=\left[\begin{array}[]{cccccc}C_{11}&C_{12}&C_{13}&0&0&0\\\
C_{12}&C_{11}&C_{13}&0&0&0\\\ C_{13}&C_{13}&C_{33}&0&0&0\\\
0&0&0&C_{55}&0&0\\\ 0&0&0&0&C_{55}&0\\\
0&0&0&0&0&C_{66}\end{array}\right],\quad C_{12}=C_{11}-2C_{66}.$ (39)

Despite our focus on inverting only five independent elastic parameters and
density, as opposed to the full set of 22 parameters, this remains an
immensely challenging task due to the substantial computational cost and
crosstalk issues inherent to multi-parameter inversion. In the abbreviated
subscript notation, the tensor $\mathbf{B}$ in equation (18) has also been
reformulated as a $6\times 6$ constant matrix. The $\mathbf{B}$ matrices for
different elastic parameters of the VTI media are shown in Appendix LABEL:B.

## Numerical examples

###  1.6 2D VTI Reservoir model

We first tested the proposed method on a simple 2D VTI Reservoir model. This
model is described by four elastic parameters ($C_{11}$, $C_{33}$, $C_{55}$,
$C_{13}$) and density $\rho$, as shown in Figure LABEL:fig:true_res_2D4. The
size of this model is 3600 m in the horizontal dimension and 900 m in the
depth dimension. It has been discretized into 180 $\times$ 45 grid blocks for
numerical computation. There are 90 receivers and 45 sources uniformly
distributed at the top of this model. A Ricker wavelet with a central
frequency of 10 Hz has been used to generate the incident wave. In this
example, we employ the sequential frequency inversion scheme, in which
frequencies are inverted individually from the lowest to the highest. The nine
frequencies used for this example are 3 Hz, 5 Hz, 7 Hz, 9 Hz, 11 Hz, 13 Hz, 15
Hz, 17 Hz, and 19 Hz. We use the GMRES (generalized minimum residual) method
with the fast Fourier transform acceleration to solve the integral equation
(9) for the calculation of displacement fields. To accurately quantify the
disparity between the calculated data $\mathbf{d}_{cal}$ and the observed data
$\mathbf{d}_{obs}$ as well as the difference between the true model
$\mathbf{m}_{true}$ and the updated model $\mathbf{m}$ during the inversion
process, we define the normalized data error $\epsilon_{d}$ and normalized
model difference $\epsilon_{m}$ as

$\epsilon_{d}=\frac{\left\|\mathbf{d}_{cal}-\mathbf{d}_{obs}\right\|}{\left\|\mathbf{d}_{obs}\right\|},$
(40)
$\epsilon_{m}=\frac{\left\|\mathbf{m}-\mathbf{m}_{true}\right\|}{\left\|\mathbf{m}_{true}\right\|}.$
(41)

The stopping criteria of model update for each frequency are either (1) data
error $\epsilon_{d}$ less than 0.001 or (2) the number of iterations reaches
10. We created the initial models, as shown in Figure LABEL:fig:ini_res_2D4,
by applying a Gaussian smoothing filter with a standard deviation of 15 to the
true models.

The inverted results are shown in Figure LABEL:fig:inv_res_2D4. All of the
model parameters are well recovered. However, we can see that the inverted
resolution of $C_{13}$ and $\rho$ in Figure LABEL:fig:inv_res_2D4 is not as
good as the other parameters. This is because the physical parameters are
coupled with each other, and the perturbation of the elastic parameters is
much larger than the perturbation of the density. In this case, the strong
mapping from other parameters greatly influences $C_{13}$ and $\rho$ during
the inversion process. This can be reduced through different parameterization
techniques (Operto et al.,, 2013; Prieux et al.,, 2013), which we will try in
our further work. Figure LABEL:fig:conv_res_2D3 shows the convergence
performance of our method. The upper and middle plots in Figure
LABEL:fig:conv_res_2D3 show the normalized data error and model difference
change with iteration. The lower plot in Figure LABEL:fig:conv_res_2D3 gives
the related frequency for each iteration. From Figure LABEL:fig:conv_res_2D3,
we can see that the data error and model difference decrease with iteration
during the inversion process at each frequency. All these results illustrate
that our method has been successfully implemented in the full waveform
inversion for the elastic anisotropic media.

true_res_2D4width=1 2D reservoir model described by four elastic parameters
and density.

ini_res_2D4width=1 Initial model with four elastic parameters and density
obtained by smoothing the true model in Figure LABEL:fig:true_res_2D4.

inv_res_2D4width=1 Inverted results of the 2D reservoir model.

conv_res_2D3width=0.9 Convergence diagrams of the 2D reservoir model test: (a)
normalized data difference $\epsilon_{d}$ versus the number of iterations, (b)
normalized model error $\epsilon_{m}$ versus the number of iterations, and (c)
the frequency corresponding to each iteration.

###  1.7 2D VTI Hess model

We have used a resampled version of the 2D VTI Hess model to test our method
on a more practical and complicated model (Figure LABEL:fig:true_Hess_2D3).
The size of this model is 3220 m in the horizontal dimension and 860 m in the
depth dimension. It has been discretized into 161 $\times$ 43 grid blocks for
numerical computation. There are 80 receivers and 40 sources uniformly
distributed at the top of this model. A Ricker wavelet with a peak frequency
of 10 Hz has been used to generate the incident wave. In this example, we
employ eight frequency components for inversion: 3 Hz, 5 Hz, 7 Hz, 9 Hz, 11
Hz, 13 Hz, 15 Hz, and 17 Hz. All the displacement fields are computed by
equation (9) with the fast-Fourier-transform-accelerated GMRES (generalized
minimum residual) method. The stopping criteria of model update for each
frequency are either (1) data residual $\epsilon_{d}$ less than 0.001 or (2)
the number of iterations reaches 10. We generate the initial models (Figure
LABEL:fig:ini_Hess_2D3) by filtering the true models with a Gaussian smoothing
kernel with a standard deviation of 14.

The inverted stiffness parameters and density are shown in Figure
LABEL:fig:inv_Hess_2D3. All parameters have been accurately estimated, and the
structure of this model has been well constructed. The fault is evident on the
right side of the inverted model. Figure LABEL:fig:conv_Hess_2D shows the
normalized data error and model difference change with iterations for
different frequencies. At each frequency, the data error and model difference
decrease with iteration during the inversion process. These inversion results
illustrate that our method can deal with a complicated model.

true_Hess_2D3width=1 Resampled 2D VTI Hess model.

ini_Hess_2D3width=1 Initial model obtained by smoothing the true model in
Figure LABEL:fig:true_Hess_2D3.

inv_Hess_2D3width=1 Inverted results of the resampled 2D VTI Hess model.

conv_Hess_2Dwidth=0.9 Convergence diagrams of the resampled 2D VTI Hess model
test: (a) normalized data difference $\epsilon_{d}$ versus the number of
iterations, (b) normalized model error $\epsilon_{m}$ versus the number of
iterations, and (c) the frequency corresponding to each iteration.

###  1.8 Modified 3D VTI Hess model

Next, we extended the test of our method from 2D to 3D. In this test, a
modified 3D Hess model is utilized, which expands upon the 2D Hess model in
the y-direction. The 3D VTI model are described by six physical parameters
($C_{11}$, $C_{33}$, $C_{55}$, $C_{66}$, $C_{13}$, and $\rho$), as shown in
Figure
LABEL:fig:true_C11_3Dvis,true_C33_3Dvis,true_C55_3Dvis,true_C66_3Dvis,true_C13_3Dvis,true_rho_3Dvis.
Here we set $C_{66}=C_{55}$, because the 2D model does not include $C_{66}$.
This modified 3D VTI Hess model’s dimensions are 2000 m $\times$1000 m
$\times$ 800 m. We discretize this model into 80$\times$40$\times$32 grid
blocks. Each grid block has a size 25m$\times$25m$\times$25m. There are 100
sources and 400 receivers uniformly distributed on the top of this model, as
shown in Figure LABEL:fig:SR7. A Ricker wavelet with a dominant frequency of
10 Hz is used to generate the incident wavefield. The initial model (Figure
LABEL:fig:ini_C11_3Dvis,ini_C33_3Dvis,ini_C55_3Dvis,ini_C66_3Dvis,ini_C13_3Dvis,ini_rho_3Dvis)
is a smooth version of the true model with a 3D Gaussian smoothing kernel with
standard deviation 11. Seven frequencies, 3 Hz, 5 Hz, 7 Hz, 9 Hz, 11 Hz, 13
Hz, and 15 Hz, are used for this inversion.

Figure
LABEL:fig:inv_C11_3Dvis,inv_C33_3Dvis,inv_C55_3Dvis,inv_C66_3Dvis,inv_C13_3Dvis,inv_rho_3Dvis
shows the 3D inverted results. To show the inner structure of the 3D inverted
results, we also show three slices in the y-direction for each inverted
parameter in Figure
LABEL:fig:inv_C11_3D7,inv_C33_3D7,inv_C55_3D7,inv_C66_3D7,inv_C13_3D7,inv_rho_3D7.
The parameters are accurately estimated, and the model’s structure is well-
defined. The inverted model clearly delineates the salt dome and the fault
structure. Figure LABEL:fig:conv_3D7 shows the normalized data error and model
difference change during the inversion process at different frequencies. At
each frequency, there is a consistent decrease in both data error and model
difference throughout the inversion. All these figures show that the proposed
method can also be used for the inversion of a 3D VTI model.

2true_C11_3Dvis,true_C33_3Dvis,true_C55_3Dvis,true_C66_3Dvis,true_C13_3Dvis,true_rho_3Dviswidth=0.45
Modified 3D VTI Hess model:(a) $C_{11}$, (b) $C_{33}$, (c) $C_{55}$, (d)
$C_{66}$, (e) $C_{13}$, and (f) $\rho$. SR7width=0.9 Distribution of sources
and receivers on the top of the 3D model.
2ini_C11_3Dvis,ini_C33_3Dvis,ini_C55_3Dvis,ini_C66_3Dvis,ini_C13_3Dvis,ini_rho_3Dviswidth=0.45
3D initial model:(a) $C_{11}$, (b) $C_{33}$, (c) $C_{55}$, (d) $C_{66}$, (e)
$C_{13}$, and (f) $\rho$.
2inv_C11_3Dvis,inv_C33_3Dvis,inv_C55_3Dvis,inv_C66_3Dvis,inv_C13_3Dvis,inv_rho_3Dviswidth=0.45
3D inverted model:(a) $C_{11}$, (b) $C_{33}$, (c) $C_{55}$, (d) $C_{66}$, (e)
$C_{13}$, and (f) $\rho$.
2inv_C11_3D7,inv_C33_3D7,inv_C55_3D7,inv_C66_3D7,inv_C13_3D7,inv_rho_3D7width=0.45
Slice visualization of the 3D inverted model:(a) $C_{11}$, (b) $C_{33}$, (c)
$C_{55}$, (d) $C_{66}$, (e) $C_{13}$, and (f) $\rho$.

conv_3D7width=0.9 Convergence diagrams of the modified 3D VTI Hess model test:
(a) normalized data difference $\epsilon_{d}$ versus the number of iterations,
(b) normalized model error $\epsilon_{m}$ versus the number of iterations, and
(c) the frequency corresponding to each iteration.

##  2 Concluding remarks

We have extended the application of the matrix-free distorted Born iterative
method to multi-parameter full waveform inversion that can reconstruct elastic
parameters and density simultaneously from frequency-domain waveform data. The
key idea of our method is the matrix-free implementations of the Fréchet
derivatives of different parameters and their adjoint operators, which can
significantly reduce the computational cost and memory demand. In multi-
parameter inversion, the Hessian information is crucial for reducing the
crosstalk effects among different parameters. The newly proposed algorithm
effectively utilizes the Hessian information by incorporating Fréchet
derivatives and their adjoint operators, eliminating the need for forming and
inverting the entire Hessian matrix. The Fréchet derivatives and their adjoint
operators are formulated as vector operations according to the physical
meaning of Green’s function. In our formulations, applying the Fréchet
derivative and its adjoint does not require knowledge of Green’s function in
the heterogeneous background media. Instead, a set of virtual displacement
fields is computed using an integral equation solver with fast Fourier
transform acceleration. All these improvements make the distorted Born
iterative method more practical for seismic applications in multi-parameter
full waveform inversion. Numerical results show that our inversion algorithm
provides promising results for seismic inversion applicable to realistic model
sizes. The 3D numerical test shows that the proposed algorithm can deal with
the 3D model (more than 0.6 million unknowns) in an ordinary computer (with an
Intel i7-7700 CPU and 64 GB of RAM). This is impossible for the conventional
distorted Born iterative method.

Although we have some favorable results, certain aspects still require
attention in future work. Firstly, it is essential to incorporate an absorbing
boundary in the numerical test to reduce boundary reflections. However,
implementing an absorbing boundary within the integral equation method remains
challenging, so it has yet to be included in this study. Secondly, it would be
interesting to explore various parametrizations (for example, Thomsen
parameters) to assess whether they can help further mitigate the crosstalk
issue.

##  3 ACKNOWLEDGEMENTS

[op_Green]2D and 3D elastic Green’s function The analytical Green’s tensor in
a 3D homogeneous isotropic medium is given by Aki and Richards, (1980) as

$\displaystyle G_{ij}(r,\omega)=$ $\displaystyle\frac{e^{i\omega
r/\alpha}}{4\pi\rho\alpha^{2}r}\left[\gamma_{i}\gamma_{j}+\left(3\gamma_{i}\gamma_{j}-\delta_{ij}\right)\left(\frac{-\alpha}{i\omega
r}\right)\right.\left.+\left(3\gamma_{i}\gamma_{j}-\delta_{ij}\right)\left(\frac{-\alpha}{i\omega
r}\right)^{2}\right]$ (42) $\displaystyle-\frac{e^{i\omega
r/\beta}}{4\pi\rho\beta^{2}r}\left[\left(\gamma_{i}\gamma_{j}-\delta_{ij}\right)+\left(3\gamma_{i}\gamma_{j}-\delta_{ij}\right)\left(\frac{-\beta}{i\omega
r}\right)+\left(3\gamma_{i}\gamma_{j}-\delta_{ij}\right)\left(\frac{-\beta}{i\omega
r}\right)^{2}\right],$

where

$r=\left|\mathbf{x}-\mathbf{x}^{\prime}\right|,\qquad\gamma_{i}=\frac{\mathbf{x}-\mathbf{x}^{\prime}}{r}\cdot\mathbf{e}_{i},\qquad\delta_{ij}=\mathbf{e}_{i}\cdot\mathbf{e}_{j}.$
(43)

In a 2D homogeneous isotropic medium, the analytical formula of Green’s tensor
is given by Sánchez-Sesma and Campillo, (1991) as

$G_{ij}(r,\omega)=\frac{1}{i8\rho}\left\\{A\delta_{ij}-B\left(2\gamma_{i}\gamma_{j}-\delta_{ij}\right)\right\\},$
(44)

where

$A=\frac{H_{0}^{(2)}(\omega r/\alpha)}{\alpha^{2}}+\frac{H_{0}^{(2)}(\omega
r/\beta)}{\beta^{2}}\text{ and }B=\frac{H_{2}^{(2)}(\omega
r/\alpha)}{\alpha^{2}}-\frac{H_{2}^{(2)}(\omega r/\beta)}{\beta^{2}},$ (45)

and $H_{0}^{(2)}$ and $H_{2}^{(2)}$ are the second kind Hankel functions of 0
order and 2 order, respectively. In the above equations, $\alpha$, $\beta$,
and $\rho$ are the homogeneous isotropic background medium’s P wave velocity,
S wave velocity, and density, respectively.

[op_ad]Derivation of the adjoint operators The Fréchet operator and its
adjoint operator satisfy the inner product rule:

$\left\langle\delta\mathbf{u},\mathcal{F}\delta\mathbf{m}\right\rangle_{D}=\left\langle\mathcal{F}^{\dagger}\delta\mathbf{u},\delta\mathbf{m}\right\rangle_{\Omega},$
(46)

where $\left\langle\cdot,\cdot\right\rangle$ is the inner product on Hilbert
space. According to the definition of inner product, we rewrite equation (46)
as

$\displaystyle\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)\equiv\int
d\mathbf{x}\left[(\mathcal{F}^{(0)})^{\dagger}\delta\mathbf{u}\right]^{*}\left(\mathbf{x}\right)\delta
m^{(0)}\left(\mathbf{x}\right),$ (47) $\displaystyle\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right)\equiv\int
d\mathbf{x}\left[(\mathcal{F}^{(p)})^{\dagger}\delta\mathbf{u}\right]^{*}\left(\mathbf{x}\right)\delta
m^{(p)}\left(\mathbf{x}\right),$

where $\left[\cdot\right]^{*}$ denotes complex conjugation. Inserting (20)
into the left hand of (47) yields

$\displaystyle\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)=\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\omega^{2}\int
d\mathbf{x}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}\right)\mathbf{u}^{(b)}\left(\mathbf{x}\right)\delta
m^{(0)}\left(\mathbf{x}\right),$ (48) $\displaystyle\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right)=-\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\int
d\mathbf{x}\nabla_{\mathbf{x}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}\right)\cdot\mathbf{B}^{(p)}:\nabla_{\mathbf{x}}\mathbf{u}^{(b)}\left(\mathbf{x}\right)\delta
m^{(p)}\left(\mathbf{x}\right).$

By reordering the integral, we obtain

$\displaystyle\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\left[\mathcal{F}^{(0)}\delta\mathbf{m}^{(0)}\right]\left(\mathbf{r}\right)=\int
d\mathbf{x}\omega^{2}\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}\right)\mathbf{u}^{(b)}\left(\mathbf{x}\right)\delta
m^{(0)}\left(\mathbf{x}\right),$ (49) $\displaystyle\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\left[\mathcal{F}^{(p)}\delta\mathbf{m}^{(p)}\right]\left(\mathbf{r}\right)=-\int
d\mathbf{x}\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\nabla_{\mathbf{x}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}\right)\cdot\mathbf{B}^{(p)}:\nabla_{\mathbf{x}}\mathbf{u}^{(b)}\left(\mathbf{x}\right)\delta
m^{(p)}\left(\mathbf{x}\right).$

Comparing the right hand term of (49) and the right hand term of (47), we have

$\displaystyle\left[(\mathcal{F}^{(0)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[\omega^{2}\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}\right)\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*},$
(50)
$\displaystyle\left[(\mathcal{F}^{(p)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[-\int
d\mathbf{r}\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\nabla_{\mathbf{x}}\mathbf{G}^{(b)}\left(\mathbf{r},\mathbf{x}\right)\cdot\mathbf{B}^{(p)}:\nabla_{\mathbf{x}}\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*}.$

Using the reciprocity of the Green’s function and moving the differentiation
operator on the Green’s function out of the integral yields

$\displaystyle\left[(\mathcal{F}^{(0)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[\omega^{2}\int
d\mathbf{r}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{r}\right)\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*},$
(51)
$\displaystyle\left[(\mathcal{F}^{(p)})^{\dagger}\delta\mathbf{u}\right]\left(\mathbf{x}\right)=\left[-\nabla_{\mathbf{x}}\int
d\mathbf{r}\mathbf{G}^{(b)}\left(\mathbf{x},\mathbf{r}\right)\delta\mathbf{u}^{*}\left(\mathbf{r}\right)\cdot\mathbf{B}^{(p)}:\nabla_{\mathbf{x}}\mathbf{u}^{(b)}\left(\mathbf{x}\right)\right]^{*}.$

[B]$\mathbf{B}$ matrices for a VTI medium The VTI medium is characterized by
five independent elastic parameters: $C_{11}$, $C_{33}$, $C_{55}$, $C_{66}$,
and $C_{13}$. Their corresponding $\mathbf{B}$ matrices are shown as follows:

$\displaystyle\mathbf{B}^{11}=\left[\begin{array}[]{cccccc}1&1&0&0&0&0\\\
1&1&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&0&0\end{array}\right],\quad\mathbf{B}^{33}=\left[\begin{array}[]{cccccc}0&0&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&0&0\end{array}\right],\qquad$ (52)
$\displaystyle\mathbf{B}^{55}=\left[\begin{array}[]{cccccc}0&0&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&0\end{array}\right],\qquad\mathbf{B}^{66}=\left[\begin{array}[]{cccccc}0&-2&0&0&0&0\\\
-2&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&0&1\end{array}\right],$
$\displaystyle\mathbf{B}^{13}=\left[\begin{array}[]{cccccc}0&0&1&0&0&0\\\
0&0&1&0&0&0\\\ 1&1&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&0&0\end{array}\right].$

## References

  * Abubakar and van den Berg, (2004) Abubakar, A., and P. M. van den Berg, 2004, Iterative forward and inverse algorithms based on domain integral equations for three-dimensional electric and magnetic objects: Journal of computational physics, 195, 236–262.
  * Aki and Richards, (1980) Aki, K., and P. G. Richards, 1980, Quantitative seismology: Theory and methods: W. R. Freeman and Company.
  * Alkhalifah et al., (2016) Alkhalifah, T., N. Masmoudi, and J.-W. Oh, 2016, A recipe for practical full-waveform inversion in orthorhombic anisotropy: The Leading Edge, 35, 1076–1083.
  * Arfken and Weber, (1999) Arfken, G. B., and H. J. Weber, 1999, Mathematical methods for physicists: American Association of Physics Teachers.
  * Auld, (1973) Auld, B. A., 1973, Acoustic fields and waves in solids: Wiley Interscience Pub.
  * Beylkin et al., (2009) Beylkin, G., C. Kurcz, and L. Monzón, 2009, Fast convolution with the free space Helmholtz Green’s function: Journal of Computational Physics, 228, 2770–2791.
  * Brossier et al., (2009) Brossier, R., S. Operto, and J. Virieux, 2009, Seismic imaging of complex onshore structures by 2d elastic frequency-domain full-waveform inversion: Geophysics, 74, WCC105–WCC118.
  * Bunks et al., (1995) Bunks, C., F. M. Saleck, S. Zaleski, and G. Chavent, 1995, Multiscale seismic waveform inversion: Geophysics, 60, 1457–1473.
  * Burridge et al., (1998) Burridge, R., M. V. de Hoop, D. Miller, and C. Spencer, 1998, Multiparameter inversion in anisotropic elastic media: Geophysical Journal International, 134, 757–777.
  * Carcione, (2014) Carcione, J. M., 2014, Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media, 3rd edition: Elsevier Science.
  * Červenỳ, (2001) Červenỳ, V., 2001, Seismic ray theory: Cambridge University Press.
  * Chew and Wang, (1990) Chew, W. C., and Y.-M. Wang, 1990, Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method: IEEE transactions on medical imaging, 9, 218–225.
  * Claerbout and Fomel, (2012) Claerbout, J. F., and S. Fomel, 2012, Image estimation by example: geophysical soundings image construction: multidimensional autoregression: Citeseer.
  * Eikrem et al., (2019) Eikrem, K. S., G. Nævdal, and M. Jakobsen, 2019, Iterated extended Kalman filter method for time-lapse seismic full-waveform inversion: Geophysical Prospecting, 67, 379–394.
  * Gubernatis et al., (1977) Gubernatis, J., E. Domany, and J. Krumhansl, 1977, Formal aspects of the theory of the scattering of ultrasound by flaws in elastic materials: Journal of Applied Physics, 48, 2804–2811.
  * Haffinger et al., (2013) Haffinger, P., A. Gisolf, and P. v. d. Berg, 2013, Towards high resolution quantitative subsurface models by full waveform inversion: Geophysical Journal International, 193, 788–797.
  * He and Plessix, (2017) He, W., and R.-É. Plessix, 2017, Analysis of different parameterisations of waveform inversion of compressional body waves in an elastic transverse isotropic earth with a vertical axis of symmetry: Geophysical Prospecting, 65, 1004–1024.
  * Hesford and Chew, (2006) Hesford, A. J., and W. C. Chew, 2006, A frequency-domain formulation of the Fréchet derivative to exploit the inherent parallelism of the distorted Born iterative method: Waves in Random and Complex Media, 16, 495–508.
  * Hesford and Chew, (2010) ——–, 2010, Fast inverse scattering solutions using the distorted Born iterative method and the multilevel fast multipole algorithm: The Journal of the Acoustical Society of America, 128, 679–690.
  * Huang et al., (2020) Huang, X., K. S. Eikrem, M. Jakobsen, and G. Nævdal, 2020, Bayesian full-waveform inversion in anisotropic elastic media using the iterated extended Kalman filter: Geophysics, 85, C125–C139.
  * Huang et al., (2019) Huang, X., M. Jakobsen, G. Nævdal, and K. S. Eikrem, 2019, Target-oriented inversion of time-lapse seismic waveform data: Communications in Computational Physics, 27, 1–27.
  * Innanen and Lira, (2010) Innanen, K. A., and J. E. Lira, 2010, Direct nonlinear Q-compensation of seismic primaries reflecting from a stratified, two-parameter absorptive medium: Geophysics, 75, V13–V23.
  * Jakobsen et al., (2020) Jakobsen, M., E. Ivan, I. Psencik, and B. Ursin, 2020, Transition operator approach to seismic full-waveform inversion in arbitrary anisotropic elastic media: Communications in Computational Physics, 27, 1–31.
  * Jakobsen and Ursin, (2012) Jakobsen, M., and B. Ursin, 2012, Nonlinear seismic waveform inversion using a Born iterative T-matrix method: SEG Technical Program Expanded Abstracts 2012, 1–5.
  * Jakobsen and Ursin, (2015) ——–, 2015, Full waveform inversion in the frequency domain using direct iterative T-matrix methods: Journal of Geophysics and Engineering, 12, 400–418.
  * Jakobsen and Wu, (2016) Jakobsen, M., and R.-S. Wu, 2016, Renormalized scattering series for frequency-domain waveform modelling of strong velocity contrasts: Geophysical Journal International, 206, 880–899.
  * Jakobsen et al., (2023) Jakobsen, M., K. Xiang, and K. van Dongen, 2023, Seismic and medical ultrasound imaging of velocity and density variations by nonlinear vectorial inverse scattering: The Journal of the Acoustical Society of America, accepted.
  * Kamath and Tsvankin, (2016) Kamath, N., and I. Tsvankin, 2016, Elastic full-waveform inversion for vti media: Methodology and sensitivity analysis: Geophysics, 81, C53–C68.
  * Köhn et al., (2012) Köhn, D., D. De Nil, A. Kurzmann, A. Przebindowska, and T. Bohlen, 2012, On the influence of model parametrization in elastic full waveform tomography: Geophysical Journal International, 191, 325–345.
  * Lee et al., (2010) Lee, H.-Y., J. M. Koo, D.-J. Min, B.-D. Kwon, and H. S. Yoo, 2010, Frequency-domain elastic full waveform inversion for vti media: Geophysical Journal International, 183, 884–904.
  * Liao and McMechan, (1996) Liao, Q., and G. A. McMechan, 1996, Multifrequency viscoacoustic modeling and inversion: Geophysics, 61, 1371–1378.
  * Lin and Huang, (2014) Lin, Y., and L. Huang, 2014, Acoustic-and elastic-waveform inversion using a modified total-variation regularization scheme: Geophysical Journal International, 200, 489–502.
  * Malovichko et al., (2017) Malovichko, M., N. Khokhlov, N. Yavich, and M. Zhdanov, 2017, Approximate solutions of acoustic 3D integral equation and their application to seismic modeling and full-waveform inversion: Journal of Computational Physics, 346, 318–339.
  * Malovichko et al., (2018) ——–, 2018, Acoustic 3D modeling by the method of integral equations: Computers & Geosciences, 111, 223–234.
  * Menke, (2012) Menke, W., 2012, Geophysical data analysis: discrete inverse theory: Academic press, 45.
  * Mora, (1988) Mora, P., 1988, Elastic wave-field inversion of reflection and transmission data: Geophysics, 53, 750–759.
  * Morse and Feshbach, (1954) Morse, P. M., and H. Feshbach, 1954, Methods of theoretical physics: McCraw-Hill Book Company Inc.
  * Oh and Alkhalifah, (2016) Oh, J.-W., and T. Alkhalifah, 2016, Elastic orthorhombic anisotropic parameter inversion: An analysis of parameterization: Geophysics, 81, C279–C293.
  * Operto et al., (2013) Operto, S., Y. Gholami, V. Prieux, A. Ribodetti, R. Brossier, L. Metivier, and J. Virieux, 2013, A guided tour of multiparameter full-waveform inversion with multicomponent data: From theory to practice: The leading edge, 32, 1040–1054.
  * Oristaglio and Blok, (2012) Oristaglio, M., and H. Blok, 2012, Wavefield imaging and inversion in electromagnetics and acoustics: Cambridge University Press.
  * Osnabrugge et al., (2016) Osnabrugge, G., S. Leedumrongwatthanakun, and I. M. Vellekoop, 2016, A convergent Born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media: Journal of computational physics, 322, 113–124.
  * Pan et al., (2016) Pan, W., K. A. Innanen, G. F. Margrave, M. C. Fehler, X. Fang, and J. Li, 2016, Estimation of elastic constants for hti media using gauss-newton and full-newton multiparameter full-waveform inversion: Geophysics, 81, R275–R291.
  * Pike and Sabatier, (2001) Pike, E. R., and P. C. Sabatier, 2001, Scattering, two-volume set: Scattering and inverse scattering in pure and applied science: Elsevier.
  * Pratt, (1990) Pratt, R. G., 1990, Inverse theory applied to multi-source cross-hole tomography.: Part 2: Elastic wave-equation method: Geophysical Prospecting, 38, 311–329.
  * Prieux et al., (2013) Prieux, V., R. Brossier, S. Operto, and J. Virieux, 2013, Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. part 1: Imaging compressional wave speed, density and attenuation: Geophysical Journal International, 194, 1640–1664.
  * Remis and Van den Berg, (2000) Remis, R. F., and P. Van den Berg, 2000, On the equivalence of the Newton-Kantorovich and distorted Born methods: Inverse Problems, 16, L1.
  * Rusmanugroho et al., (2017) Rusmanugroho, H., R. Modrak, and J. Tromp, 2017, Anisotropic full-waveform inversion with tilt-angle recovery: Geophysics, 82, R135–R151.
  * Sánchez-Sesma and Campillo, (1991) Sánchez-Sesma, F. J., and M. Campillo, 1991, Diffraction of P, SV, and Rayleigh waves by topographic features: A boundary integral formulation: Bulletin of the seismological Society of America, 81, 2234–2253.
  * Shekhar et al., (2023) Shekhar, U., M. Jakobsen, E. Iversen, I. Berre, and F. A. Radu, 2023, Integral equation method for microseismic wavefield modelling in anisotropic elastic media: arXiv preprint arXiv:2301.12836.
  * Song et al., (2020) Song, X., M. Li, F. Yang, S. Xu, and A. Abubakar, 2020, Study on 3-D acoustic imaging for human thorax based on contrast source inversion: IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 67, 1533–1543.
  * Stolt and Jacobs, (1980) Stolt, R. H., and B. Jacobs, 1980, Inversion of seismic data in a laterally heterogeneous medium: SEP Rep, 24, 135–52.
  * Tarantola, (1984) Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49, 1259–1266.
  * Tarantola, (2005) ——–, 2005, Inverse problem theory and methods for model parameter estimation: SIAM.
  * van den Berg, (2021) van den Berg, P. M., 2021, Forward and inverse scattering algorithms based on contrast source integral equations: John Wiley & Sons.
  * Van Den Berg and Kleinman, (1997) Van Den Berg, P. M., and R. E. Kleinman, 1997, A contrast source inversion method: Inverse problems, 13, 1607.
  * Virieux and Operto, (2009) Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74, WCC1–WCC26.
  * Wang and Cheng, (2017) Wang, T., and J. Cheng, 2017, Elastic full waveform inversion based on mode decomposition: The approach and mechanism: Geophysical Journal International, 209, 606–622.
  * Wang and Chew, (1989) Wang, Y., and W. C. Chew, 1989, An iterative solution of the two-dimensional electromagnetic inverse scattering problem: International Journal of Imaging Systems and Technology, 1, 100–108.
  * Warner et al., (2013) Warner, M., A. Ratcliffe, T. Nangoo, J. Morgan, A. Umpleby, N. Shah, V. Vinje, I. Stekl, L. Guasch, C. Win, et al., 2013, Anisotropic 3d full-waveform inversion, geophysics, 78, r59–r80: LIST OF FIGURES.
  * Weglein et al., (2003) Weglein, A. B., F. V. Araújo, P. M. Carvalho, R. H. Stolt, K. H. Matson, R. T. Coates, D. Corrigan, D. J. Foster, S. A. Shaw, and H. Zhang, 2003, Inverse scattering series and seismic exploration: Inverse problems, 19, R27.
  * Weglein et al., (1981) Weglein, A. B., W. Boyse, and J. Anderson, 1981, Obtaining three-dimensional velocity information directly from reflection seismic data: An inverse scattering formalism: Geophysics, 46, 1116–1120.
  * Xiang et al., (2023) Xiang, K., M. Jakobsen, K. S. Eikrem, and G. Nævdal, 2023, A matrix-free variant of the distorted Born iterative method for seismic full-waveform inversion: Geophysical Prospecting, 71, 431–442.
  * Yang et al., (2019) Yang, Q., A. Malcolm, H. Rusmanugroho, and W. Mao, 2019, Analysis of radiation patterns for optimized full waveform inversion in fluid-saturated porous media: Geophysical Journal International, 216, 1919–1937.
  * Zhou et al., (1995) Zhou, C., W. Cai, Y. Luo, G. T. Schuster, and S. Hassanzadeh, 1995, Acoustic wave-equation traveltime and waveform inversion of crosshole seismic data: Geophysics, 60, 765–773.