Target-oriented least-squares reverse-time migration using Marchenko double-
focusing: reducing the artifacts caused by overburden multiples–LABEL:lastpage

# Target-oriented least-squares reverse-time migration using Marchenko double-
focusing: reducing the artifacts caused by overburden multiples

Aydin Shoja1 Joost van der Neut1 Kees Wapenaar1
1 Delft University of Technology  Department of Geoscience and Engineering
Delft  The Netherlands.
Email<EMAIL_ADDRESS>
###### keywords:

Inverse theory; Numerical modelling; Acoustic properties; Controlled source
seismology; Wave scattering and diffraction

Geophysicists have widely used Least-squares reverse-time migration (LSRTM) to
obtain high-resolution images of the subsurface. However, LSRTM is
computationally expensive and it can suffer from multiple reflections.
Recently, a target-oriented approach to LSRTM has been proposed, which focuses
the wavefield above the target of interest. Remarkably, this approach can be
helpful for imaging below complex overburdens and subsalt domains. Moreover,
this approach can significantly reduce the computational burden of the problem
by limiting the computational domain to a smaller area. Nevertheless, target-
oriented LSRTM still needs an accurate velocity model of the overburden to
focus the wavefield accurately and predict internal multiple reflections
correctly. A viable alternative to an accurate velocity model for internal
multiple prediction is Marchenko redatuming. This method is a novel data-
driven method that can predict Green’s functions at any arbitrary depth,
including all orders of multiples. The only requirement for this method is a
smooth background velocity model of the overburden. Moreover, with Marchenko
double-focusing, one can make virtual sources and receivers at a boundary
above the target and bypass the overburden.

This paper proposes a new algorithm for target-oriented LSRTM, which fits the
Marchenko double-focused data with predicted data. The predicted data of the
proposed method is modeled by a virtual source term created by Marchenko
double-focusing on a boundary above the target of interest. This virtual
source term includes all the interactions between the target and the
overburden. Moreover, the Marchenko double-focused data and the virtual source
term are free of multiples generated in the overburden. Consequently, our
target-oriented LSRTM algorithm suppresses the multiples purely generated
inside the overburden. It correctly accounts for all orders of multiples
caused by the interactions between the target and the overburden, resulting in
a significant reduction of the artifacts caused by the overburden internal
multiple reflections and improves amplitude recovery in the target image
compared to conventional LSRTM.

##  1 Introduction

In seismic exploration, geophysicists estimate the Earth’s wave propagation
parameters, such as the subsurface velocity and density. To evaluate these
parameters, an array of sources and receivers is deployed at the surface of
the Earth to generate and record seismic waves. A common approach to retrieve
the subsurface parameters is to consider the model of the subsurface as a
superposition of a background model ($m_{0}$), which corresponds to the long
wavelengths, and a short-wavelength reflectivity model ($\delta m$). This
scale separation is based on a weak-scattering assumption [Schuster 2017]. In
seismic migration, we aim to obtain a structural image of $\delta m$.

Migration is an imaging process based on a linear relation between the
recorded data and the reflectivity model ($\delta m$). One of the common
migration techniques is called Reverse Time Migration (RTM) [Baysal et al.
1983, Zhou et al. 2018, Zhang et al. 2018]. In RTM, the underlying linear
relation is the Born integral, so the RTM image is produced by applying the
adjoint of the Born operator to the observed data. However, this approach may
suffer from many errors such as noise [Dutta et al. 2017], temporal band
limitations, limited recording aperture [Tang 2009], spatial aliasing, and
multiple reflections [Liu et al. 2016].

One way to reduce some of these errors is to resolve $\delta m$ by least-
squares inversion. This approach, known as Least-Squares reverse-time
migration (LSRTM), is the most common algorithm for producing a high-
resolution reflectivity model of the subsurface. Nevertheless, LSRTM is
computationally expensive [Dai et al. 2012] as it requires solving the wave
equation and its adjoint for a large model in many iterations and storing
large data volumes [Tang 2009, Herrmann & Li 2012, Farshad & Chauris 2021].
Reducing the computational burden of LSRTM is possible by reducing the
computational domain. In other words, by limiting ourselves to a relatively
small target region of interest. This so-called ”target-oriented” approach
[Valenciano et al. 2006, Haffinger et al. 2013, Willemsen et al. 2016, Yuan et
al. 2017, Zhao & Sen 2018, Guo & Alkhalifah 2020] divides the medium into two
parts: an overburden and a target. Target-oriented LSRTM (TOLSRTM) achieves
this goal by bypassing the overburden and focusing the wavefield above the
target of interest. This process, in which a wavefield is extrapolated from
its current acquisition surface to another surface, is called redatuming
[Berryhill 1984, Schuster & Zhou 2006, Luiken & van Leeuwen 2020, Barrera et
al. 2021].

Target-oriented LSRTM, in addition to its lower computational costs, can also
be helpful in cases the overburden produces strong internal multiple
reflections that mask the reflections of the target region. Nonetheless, for
this method to work correctly, a model of the overburden, which can produce
all orders of reflections, is needed. Researchers have shown a strong interest
in solving the problems caused by inaccurate overburden models and the
artifacts that internal multiple reflections can create in LSRTM images. For
example, Guo and Alkhalifah Guo2020,Guo2019 introduced a least-squares
waveform redatuming method to redatum the data to the boundary of the target
and combined it with full waveform inversion.

Among different redatuming methods, Marchenko redatuming showed itself as a
promising data-driven tool for predicting the Green’s functions at the
boundary of the target region with all orders of internal multiples [Ravasi et
al. 2016, Vargas et al. 2021, Dukalski & de Vos 2018, Wapenaar et al. 2014].
The Marchenko redatuming can predict the internal multiples by only having the
direct arrival between the surface and the target boundary in a smooth
background velocity model of the overburden. Once we solve the Marchenko
equations, we get the response of the medium stimulated by virtual sources at
the boundary of the target and recorded by receivers at the surface, including
all orders of multiple reflections. Hence, the first step of Marchenko
redatuming can be interpreted as source-side redatuming in the physical medium
(including reflections in the overburden that are not encoded in the
background model). Diekmann et al. Diek used Marchenko equations to predict
the in-volume Green’s function of a target of interest and used this wavefield
as the incident wavefield inside the scattering integral. Cui et al. Cui used
the reciprocity theorem to derive a forward modeling equation based on
Marchenko wavefields for Full waveform inversion. Moreover, one of the
byproducts of the Marchenko integrals is the downgoing part of the focusing
function, which is the inverse of the transmission response of the overburden.
This focusing function can also be applied at the receiver side to achieve
double-sided (meaning source and receiver) redatuming [Staring et al. 2018].
The process of double-sided redatuming with the help of the Marchenko focusing
function is called Marchenko double-focusing.

In this paper we consider the upgoing Marchenko double-focused data as the
”observed data” that we feed to an LSRTM algorithm. To complete our TOLSRTM
algorithm, we formulate a forward modeling to match the Marchenko double-
focused data with the predicted data of the target region. This forward
modeling is based on the downgoing wavefields from the Marchenko equation and
a Born model of the target. This Marchenko downgoing wavefield acts as virtual
source for the modeling and the inversion.

The organization of this paper is as follows: First, we briefly review least-
squares reverse-time migration and target-oriented LSRTM. Second, Marchenko
redatuming and double-focusing are explained. Then, our method of combining
target-oriented LSRTM and Marchenko double-focusing is introduced. In the
results section, we show numerical examples of our method and we exhibit the
consequences of using an inaccurate velocity model for redatuming and TOLSRTM.
Finally, results and potential future directions are discussed.

##  2 Theory

###  2.1 Least-squares reverse-time migration

We start to formulate our target-oriented algorithm by briefly reviewing
classical reverse-time migration with sources and receivers at the acquisition
surface. Classical RTM is based on the Born approximation where the incident
wavefield ($P^{inc}$) is approximated by the background Green’s function

$P^{inc}(\textbf{x},\textbf{x}_{s},\omega)\approx
G_{0}(\textbf{x},\textbf{x}_{s},\omega)W(\omega).$ (1)

Here $W(\omega)$ is the wavelet signature, and $\omega$ is the angular
frequency. Moreover, $\textbf{x}=(x_{1},x_{2},x_{3})$ is a location inside the
medium, $\textbf{x}_{s}$ is the source location, and
$G_{0}(\textbf{x},\textbf{x}_{s},\omega)$ is the Fourier transform of the
causal solution of the scalar wave equation in the background velocity model,
which in the frequency domain and with a constant density ($\rho_{0}$) is
equal to the Helmholtz equation:

$\nabla^{2}G_{0}+k_{0}^{2}G_{0}=-i\omega\rho_{0}\delta(\textbf{x}-\textbf{x}_{s}),$
(2)

where $k_{0}(\textbf{x})=\frac{\omega}{c_{0}(\textbf{x})}$ is the wavenumber,
$c_{0}$ is background velocity, and $i$ is the imaginary unit.

The reflectivity model is defined as $\delta
m=(\frac{1}{c^{2}}-\frac{1}{c_{0}^{2}})$, where $c$ is the perturbed velocity.
Here, $\delta m$ is connected to the scattered data ($P^{scat}$) through the
following linear equation [Born & Wolf 1999]:

$P^{scat}_{pred}(\textbf{x}_{r},\textbf{x}_{s},\delta
m,\omega)=\frac{\omega^{2}}{\rho_{0}}\int_{V}G_{0}(\textbf{x}_{r},\textbf{x},\omega)\delta
m(\textbf{x})G_{0}(\textbf{x},\textbf{x}_{s},\omega)W(\omega)\,d\textbf{x}.$
(3)

This integral is taken over the volume of the model ($V$). Subscripts ”$r$”
and ”$s$” are used to specify receiver and source, respectively. Equation 3
can be written in the operator form as:

$P^{scat}_{pred}(\textbf{x}_{r},\textbf{x}_{s},\delta
m,\omega)=\mathcal{L}\delta m.$ (4)

Here $\mathcal{L}$ is the forward Born operator.

In the standard reverse-time migration method, the reflectivity model is
approximately retrieved by applying the adjoint of $\mathcal{L}$ to the
observed scattered data:

$\delta m^{mig}(\textbf{x})=\mathcal{L}^{\dagger}P^{scat}_{obs}.$ (5)

Green’s functions ($G_{0}$) which constitute the kernel $\mathcal{L}$, only
consider wave propagation in the smooth background and do not contain
reflections from the perturbations. As a result, this kernel is unable to
interpret the multiple-scattered waves correctly. Moreover, since the adjoint
of this kernel is only an approximation of its inverse, the resolution of the
retrieved reflectivity model is low.

To overcome the resolution issue, researchers have utilized a least-squares
approach by replacing the adjoint ($\mathcal{L}^{\dagger}$) with a damped
least-squares solution [Marquardt 1963, Dutta et al. 2017]:

$\delta
m^{mig}=[\mathcal{L}^{\dagger}\mathcal{L}+\epsilon]^{-1}\mathcal{L}^{\dagger}P^{scat}_{obs}.$
(6)

However, computing the Hessian matrix ($\mathcal{L}^{\dagger}\mathcal{L}$) and
its inverse is prohibitive. Alternatively, updating the reflectivity model is
preferred by an iterative algorithm that minimizes the L2-norm of the
difference between observed and predicted data

$C(\delta m)=\frac{1}{2}\big{\|}P^{scat}_{pred}(\delta
m)-P^{scat}_{obs}\big{\|}_{2}^{2}.$ (7)

We can solve this optimization problem for instance by a conjugate gradient
algorithm. Since in least-squares reverse-time migration, the background
velocity model ($c_{0}$) is kept unchanged, and only the reflectivity model
($\delta m$) is updated, the Green’s functions of Equation 3 are only computed
once. For an overview of least-squares reverse-time migration, we refer to
Schuster Schus.

###  2.2 Target-oriented least-squares reverse time migration

Researchers have proposed to implement the imaging and inversion in a target-
oriented manner [Valenciano et al. 2006, Haffinger et al. 2013, Willemsen et
al. 2016, Yuan et al. 2017, Zhao & Sen 2018, Guo & Alkhalifah 2020]. In this
way, the wavefield can be focused below a complex overburden to produce
virtual sources and receivers right above the target [Berryhill 1984, Schuster
& Zhou 2006, Luiken & van Leeuwen 2020, Barrera et al. 2021]. These redatumed
wavefields can be used for imaging and inversion. The advantages of this
approach are as follows: 1) The wavefield is focused below specific
structures, such as salt domes, and 2) the computational costs are lower
because the computational domain can be reduced significantly. However,
current methods need a reasonable estimation of the velocity model of the
overburden (i.e., the part of the model above the target of interest), which
can predict internal multiple reflections. Otherwise, these overburden-
generated multiple reflections produce artifacts inside the retrieved image of
the target.

Marchenko redatuming enables us to take these overburden-generated multiples
into account and improve the quality of the target-oriented LSRTM images. The
following section gives a brief review of Marchenko redatuming by double-
focusing.

###  2.3 Marchenko redatuming by double-focusing

Marchenko redatuming is a novel data-driven method that can retrieve the
Green’s function at a surface above the target area with all orders of its
multiple-scattered events. The only requirements for this method are the
reflection response at the surface and a smooth background velocity model of
the overburden that can predict the direct arrival from the surface to the
redatuming level.

In order to retrieve Green’s functions at the redatuming level, the coupled
Marchenko-type representations are solved iteratively. These equations are
[Wapenaar et al. 2014]

$\begin{split}G^{-}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)&=\int_{\mathcal{D}_{acq}}R(\textbf{x}_{r},\textbf{x}_{s},\omega)f_{1}^{+}(\textbf{x}_{s},\textbf{x}_{v},\omega)\,d\textbf{x}_{s}\\\
&-f_{1}^{-}(\textbf{x}_{r},\textbf{x}_{v},\omega),\end{split}$ (8)

and

$\begin{split}G^{+}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)&=-\int_{\mathcal{D}_{acq}}R(\textbf{x}_{r},\textbf{x}_{s},\omega)f_{1}^{-}(\textbf{x}_{s},\textbf{x}_{v},\omega)^{*}\,d\textbf{x}_{s}\\\
&+f_{1}^{+}(\textbf{x}_{r},\textbf{x}_{v},\omega)^{*}.\end{split}$ (9)

In these equations $\mathcal{D}_{acq}$ is the acquisition surface where
$\textbf{x}_{s}$ and $\textbf{x}_{r}$ are located, $G^{-}_{Mar}$ and
$G^{+}_{Mar}$ are up-going and down-going parts of the Marchenko redatumed
Green’s function (Fig. 1(a) and 1(b)), respectively. Further, $f_{1}^{-}$ is
the up-going part, and $f_{1}^{+}$ is the down-going part of the focusing
function, and the subscript ”$v$” stands for a virtual point, which is located
on the redatuming level. We denote this redatuming level by
$\mathcal{D}_{tar}$. Moreover, $R(\textbf{x}_{r},\textbf{x}_{s},\omega)$ is
the dipole response of the medium at the acquisition surface. $R$ is related
to the upgoing Green’s function ($G^{-}$) with the following relation:

$R(\textbf{x}_{r},\textbf{x}_{s},\omega)=\frac{\partial_{3,s}G^{-}(\textbf{x}_{r},\textbf{x}_{s},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{s})},$
(10)

where $\partial_{3,s}$ is the partial derivative in the downward direction
taken at $\textbf{x}_{s}$. Note that before inserting
$R(\textbf{x}_{r},\textbf{x}_{s},\omega)$ into Equations 8 and 9, we need to
remove horizontally propagating waves and surface-related multiples. For more
information on the derivation of these integrals and their solution for
computing the focusing functions and Green’s functions, please refer to
Wapenaar et al. Kees and Thorbecke et al. Jan.

The above-mentioned equations can be interpreted as source-side redatuming, as
the receivers are still located at the acquisition surface. To relocate
receivers to the target level, Staring et al. Myr suggested convolving the up-
going and down-going parts of the Marchenko redatumed Green’s function with
the (filtered) downgoing focusing function in a multi-dimensional sense:

$G^{-,+}_{df}(\textbf{x}_{v},\textbf{x}^{\prime}_{v},\omega)=\int_{\mathcal{D}_{acq}}G^{-}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)\mathcal{F}_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)\,d\textbf{x}_{r},$
(11)

and

$G^{+,+}_{df}(\textbf{x}_{v},\textbf{x}^{\prime}_{v},\omega)=\int_{\mathcal{D}_{acq}}G^{+}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)\mathcal{F}_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)\,d\textbf{x}_{r},$
(12)

where

$\mathcal{F}_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)=\frac{\partial_{3,r}f_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{r})}.$
(13)

Here the vertical derivative is taken with respect to $\textbf{x}_{r}$. In
Equations 11 and 12, the first superscript on the left-hand side expresses the
direction of the propagation at the receiver location, the second one shows
the same at the source location, and ”$df$” means ”double-focused.” The above-
mentioned procedure is called ”Marchenko double-focusing.”

Through Marchenko double-focusing, we retrieve a down-going Green’s function
($G^{+,+}_{df}$) which contains a band-limited delta function and interactions
between target and overburden, and an up-going Green’s function
($G^{-,+}_{df}$), which is the up-going response to $G^{+,+}_{df}$ from the
target at the redatuming level, still including interactions between target
and overburden on the source side (Fig. 1(c) and 1(d)). In comparison,
”conventional double-focusing” is the act of using the inverse of the direct
arrival of the transmission response of the overburden instead of the
downgoing Marchenko focusing function for double-focusing. Considering that
the inverse of the direct arrival of the transmission response of the
overburden does not contain multiple reflections, it is unable to predict and
remove the multiples generated by the overburden. In the next section,
”double-focusing” is a general term referring to both methods, and it is
mentioned wherever a distinction is necessary.

(a)

(b)

(c)

(d)

Figure 1: The Green’s functions resulting from Marchenko redatuming and
double-focusing. a) Downgoing part of the Marchenko Green’s function, b)
upgoing part of the Marchenko Green’s function, c) downgoing Marchenko double-
focused Green’s function and d) upgoing Marchenko double-focused Green’s
function

###  2.4 Target-oriented LSRTM by Marchenko double-focusing

The up-going Green’s function is related to the down-going Green’s function at
the redatuming level with this equation [Wapenaar et al. 2000]:

$G^{-,+}(\textbf{x}_{vr},\textbf{x}^{\prime}_{vs},\omega)=\int_{\mathcal{D}_{tar}}R^{tar}(\textbf{x}_{vr},\textbf{x}_{vs},\omega)G^{+,+}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)\,d\textbf{x}_{vs}.$
(14)

Here $\textbf{x}_{vs}$ and $\textbf{x}_{vr}$ are located at
$\mathcal{D}_{tar}$, and ”$vr$” and ”$vs$” stand for virtual receiver/source.
In addition, $R^{tar}(\textbf{x}_{vr},\textbf{x}_{vs},\omega)$ is the dipole
reflection response of the target at the redatuming level with a reflection-
free half-space as the overburden. Substitution of the double-focused Green’s
functions into Equation 14 lays the foundation of our target-oriented scheme.
However, double-focused data suffers from limited aperture, meaning we can not
retrieve the far-offset section of the data at the depth ($R^{tar}$) simply
because the large propagation angles do not reach the surface within the
limits of the acquisition aperture.

To address this issue, consider we would put actual sources and receivers at
the target boundary. The Green’s functions of this configuration are related
to the double-focused Green’s functions via

$G^{-,+}_{df}(\textbf{x}^{\prime}_{vr},\textbf{x}_{vs},\omega)\approx\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)G^{-,+}(\textbf{x}_{vr},\textbf{x}_{vs},\omega)\,d\textbf{x}_{vr}.$
(15)

Here

$\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)=\int_{\mathcal{D}_{acq}}G_{d}^{+}(\textbf{x}^{\prime}_{vr},\textbf{x}_{s},\omega)^{-1}G_{d}^{+}(\textbf{x}_{vr},\textbf{x}_{s},\omega)\,d\textbf{x}_{s}$
(16)

is a point spread function (PSF) with characteristics of a band-limited delta
function that imposes temporal and spatial band limitations to the upgoing
Green’s function. In this equation, $G_{d}^{+}$ is the first arrival of the
Green’s function between the target boundary and the surface. Gamma is a
filter that removes high propagation angles from the data. In other words,
Gamma takes the full band reflection response of the target ($R^{tar}$) once
to the surface by applying $G_{d}^{+}$ to it, and then takes it back to the
focusing level by applying the inverse of the $G_{d}^{+}$ to the result of the
previous step. In this way, $R^{tar}$ goes under the same process as the other
double-focused wavefields.

To make Equation 14 compatible with the double-focused Green’s functions, we
assume

$G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)\approx
G^{+,+}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega),$ (17)

and then we convolve both sides of it with this PSF and use Equation 15 to
reach

$\begin{split}&G^{-,+}_{df}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\omega)=\\\
&\int_{\mathcal{D}_{tar}}\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)R^{tar}(\textbf{x}_{vr},\textbf{x}_{vs},\omega)G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)\,d\textbf{x}_{vr}\,d\textbf{x}_{vs}.\end{split}$
(18)

In other words, $G^{-,+}_{df}$ is related to the $G^{+,+}_{df}$ with a
temporal and spatial band-limited version of
$R^{tar}(\textbf{x}_{vr},\textbf{x}_{vs},\omega)$, which we can write as:

$\hat{R}^{tar}(\textbf{x}^{\prime}_{vr},\textbf{x}_{vs},\omega)=\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)R^{tar}(\textbf{x}_{vr},\textbf{x}_{vs},\omega)\,d\textbf{x}_{vr},$
(19)

and subsequently rewrite equation 18 as:

$G^{-,+}_{df}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\omega)=\int_{\mathcal{D}_{tar}}\hat{R}^{tar}(\textbf{x}^{\prime}_{vr},\textbf{x}_{vs},\omega)G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)\,d\textbf{x}_{vs}.$
(20)

By taking $\hat{R}^{tar}$ as a function of $\delta m$ of the target, this
equation will become the forward modeling equation for our target-oriented
least-squares reverse-time migration. We take the observed scattered wavefield
as

$\hat{P}^{scat}_{obs}=G^{-,+}_{df}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\omega)W(\omega),$
(21)

and the predicted scattered wavefield as

$\begin{split}&\hat{P}^{scat}_{pred}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\delta
m,\omega)=\\\
&\int_{\mathcal{D}_{tar}}\hat{R}^{tar}_{pred}(\textbf{x}^{\prime}_{vr},\textbf{x}_{vs},\delta
m,\omega)G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)W(\omega)\,d\textbf{x}_{vs}.\end{split}$
(22)

Additionally, it is possible to take the $G^{+,+}_{df}$ as a function of
$\delta m$ and solve a non-linear problem for applications where the target is
subject to changes [Shoja et al. 2020, Wapenaar & staring 2018].

To demonstrate the relation between Equation 22 and Equation 3, first,
consider Equation 19 and insert it into Equation 22

$\begin{split}&\hat{P}^{scat}_{pred}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\delta
m,\omega)=\\\
&\int_{\mathcal{D}_{tar}}\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)R^{tar}_{pred}(\textbf{x}_{vr},\textbf{x}_{vs},\delta
m,\omega)G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)W(\omega)\,d\textbf{x}_{vs}\,d\textbf{x}_{vr}.\end{split}$
(23)

Next, consider

$\begin{split}R^{tar}_{pred}(\textbf{x}_{vr},\textbf{x}_{vs},\delta
m,\omega)=\frac{\omega^{2}}{\rho_{0}}\int_{V}G_{0}(\textbf{x}_{vr},\textbf{x},\omega)\delta
m(\textbf{x})\frac{\partial_{3,vs}G_{0}(\textbf{x},\textbf{x}_{vs},\omega)}{\frac{1}{2}i\omega\rho_{0}}\,d\textbf{x},\end{split}$
(24)

where we used Equation 10 to rewrite Equation 3 for reflection response
(dipole response) instead of scattered pressure wavefield (monopole response).
Then, we insert Equation 24 into Equation 23, and then we rearrange the
integrals to reach

$\begin{split}&\hat{P}^{scat}_{pred}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\delta
m,\omega)=\\\
&\frac{\omega^{2}}{\rho_{0}}\int_{\nu}\left[\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)G_{0}(\textbf{x}_{vr},\textbf{x},\omega)\,d\textbf{x}_{vr}\right]\\\
&\times\delta
m(\textbf{x})\left[\int_{\mathcal{D}_{tar}}\frac{\partial_{3,vs}G_{0}(\textbf{x},\textbf{x}_{vs},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{vs})}G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)W(\omega)\,d\textbf{x}_{vs}\right]\,d\textbf{x}.\end{split}$
(25)

We redefine the term in the first square bracket as

$\hat{G}_{0}(\textbf{x}^{\prime}_{vr},\textbf{x},\omega)=\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)G_{0}(\textbf{x}_{vr},\textbf{x},\omega)\,d\textbf{x}_{vr},$
(26)

which is the spatial and temporal band limited Green’s function of the target
area, and the term in the second square bracket as

$\begin{split}P^{inc}_{df}(\textbf{x},&\textbf{x}^{\prime}_{vs},\omega)=\\\
&\int_{\mathcal{D}_{tar}}\frac{\partial_{3,vs}G_{0}(\textbf{x},\textbf{x}_{vs},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{vs})}G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)W(\omega)\,d\textbf{x}_{vs},\end{split}$
(27)

which is the incident wavefield computed by a Marchenko double-focused
downgoing wavefield, so that it includes the interactions between the target
and the overburden. Then we reach the following equation, which has a similar
structure as Equation 3:

$\hat{P}^{scat}_{pred}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\delta
m,\omega)=\frac{\omega^{2}}{\rho_{0}}\int_{\nu}\hat{G}_{0}(\textbf{x}^{\prime}_{vr},\textbf{x},\omega)\delta
m(\textbf{x})P^{inc}_{df}(\textbf{x},\textbf{x}^{\prime}_{vs},\omega)d\textbf{x},$
(28)

where $\nu$ is the target volume. Finally, the new cost function is

$C(\delta m)=\frac{1}{2}\big{\|}\hat{P}^{scat}_{pred}(\delta
m)-\hat{P}^{scat}_{obs}\big{\|}_{2}^{2},$ (29)

which we solve with a conjugate gradient algorithm.

It is possible to compute Equations 27 and 28 in two ways. The first way is to
compute the background Green’s function ($G_{0}$) with an algorithm of choice,
e.g. finite-difference, and then insert it in Equations 27 and 28. The second
way is to inject $G^{+,+}_{df}$ as a downward radiating dipole line source for
each $\textbf{x}^{\prime}_{vs}$ location in the target background medium with
a numerical PDE solver to compute $P^{inc}_{df}$. Then inject
$\frac{\omega^{2}}{\rho_{0}}\delta m(\textbf{x})P^{inc}_{df}$ as a contrast
source to compute $\hat{P}^{scat}_{pred}$. The decision of choosing the
appropriate computation method is based on the dimensions of the target and
the number of virtual sources.

##  3 Numerical results

###  3.1 Marchenko double-focusing vs. Conventional double-focusing

To validate our theory, we use two different synthetic models. We compute the
reflection data sets at the acquisition surface using a finite-difference
algorithm [Thorbecke & Draganov 2011] and a dipole (vertical force) source
with a flat spectrum wavelet (0 to 100 $Hz$). The injection rate is
represented as an $i\omega$ factor in Equation 2. Further, we use a spatial
sampling of 2.5 $m$ in both directions and a time sampling of $0.4$ $ms$. To
reduce the data size, we set the receiver temporal sampling to $4$ $ms$. For
the sake of better visualization, we convolve all time-space records with a 30
$Hz$ Ricker wavelet.

For computational simplicity, we define the target area of both models as a
homogeneous background and constant density with an embedded velocity
perturbation. Consequently, we use the analytical solution of the Helmholtz
equation for a homogeneous 2D medium to compute the Green’s functions for
TOLSRTM [Morse & Feshbach 1953, van den Berg 2021]:

$G_{0}(\textbf{x},\textbf{x}_{s},\omega)=\frac{1}{2\pi}K_{0}(\frac{-i\omega}{c_{0}}|\textbf{x}-\textbf{x}_{s}|).$
(30)

Here $K_{0}$ is the modified Bessel function of the second kind and zeroth
order. In the TOLSRTM algorithm, the spatial sampling is set to 5 $m$ and the
time sampling is 4 $ms$.

Note that we do not commit an inverse crime in these examples. As we mentioned
earlier, to compute $\hat{P}^{scat}_{obs}$, we compute the reflection response
at the acquisition surface by a finite-difference algorithm. Then we apply
Marchenko double-focusing to this data set. On the other hand, to compute
$\hat{P}^{scat}_{pred}$ we use Equations 30, 28, and 27. In all of the
mentioned steps, we deal with convolutions that lead to a loss of resolution.
To avoid it, we compute everything with a flat spectrum wavelet (0 to 100
$Hz$), and then we convolve both $\hat{P}^{scat}_{obs}$ and
$\hat{P}^{scat}_{pred}$ with a 30 $Hz$ Ricker wavelet.

For comparison, we compute the redatuming operator in two different ways for
both models. First, we isolate the direct arrival of the transmission response
of the overburden with a time window and use it as an estimation of the
redatuming operator and redatum the surface data with it. This approach is
equivalent to conventional redatuming of primary reflections [Berryhill 1984].
Second, we solve the coupled Marchenko equations with help of this direct
arrival to compute the redatumed Green’s and focusing functions [Thorbecke et
al. 2017]. We calculate the direct arrival in this way to avoid kinematic
errors. We discuss the issue of kinematic errors later in section 3.3.

####  3.1.1 Simple model

As a first example, we define a model with a single-layered overburden and a
single small diffractor with dimensions of 5 meters by 5 meters inside a
homogeneous background as the target (Fig. 2). The background velocity is set
to 1500 $m/s$, and the velocity of the overburden layer is 3000 $m/s$. The
density is constant and equal to 2500 $kg/m^{3}$ everywhere. For this model,
101 sources and receivers with a spacing of 10 meters are used. The number and
spacing of real and virtual sources and receivers are the same. In Figure 2
blue crosses are the actual source and receiver locations, and the red dots
are the virtual source and receivers locations. Moreover, the overburden is
designed such that a higher-order multiple reflection coincides with the
primary event of the diffractor. We design it in this way to show that our
method can handle the overburden-related multiple reflections correctly.

Figure 2: Velocity model with a single-layered overburden and a small
diffractor in the target.

In Figure 3, the common-source record at the acquisition surface is shown for
this model. The diffractor response of the simple model is visible around 0.5
seconds in this figure and it partly coincides with a higher-order multiple
reflection.

Figure 3: Common-source record at the surface of a source located at lateral
distance = 0 in the simple model. A linear time-varying ($\frac{t}{dt}$) gain
is applied to this record to amplify weaker events.

In Figure 4, we show the double-focused data of this model. From here onward,
we use the subscript ’Cdf’ for the data produced by ’conventional double-
focusing’, and ’Mdf’ for ’Marchenko double focused’ data. In Figure 4(b), it
is visible that the overburden-related multiple reflections have been
suppressed by applying Marchenko double focusing, yet the interactions between
the overburden and the target are still present. Moreover, in Figure 5, we
show the double-focused down-going Green’s function (Eq. 12). The downgoing
wavefield of the Marchenko double-focusing (Fig. 5(b)) contains the higher-
order interactions between the target and the overburden, whereas the
downgoing wavefield of the conventional double-focusing does not include these
interactions. Injecting the downgoing Marchenko double-focused wavefield into
the target from the focusing depth reproduces these high-order interactions in
the predicted data of the TOLSRTM algorithm.

(a)

(b)

Figure 4: Double-focused (upgoing) data of model with a single-layered
overburden corresponding to a virtual source at the middle. a) Conventional
double-focused data ($G^{-,+}_{Cdf}$) by applying conventional redatuming
operators, and b) Marchenko double-focused data ($G^{-,+}_{Mdf}$). The maximum
and the minimum value of the grey-level scale of both figures are the same and
a linear time-varying ($\frac{t}{dt}$) gain is applied to both records to
amplify weaker events.

(a)

(b)

Figure 5: Virtual source wavefield of the model with a single-layered
overburden at the middle virtual source location. a) the conventional double-
focused source ($G^{+,+}_{Cdf}$) by applying conventional redatuming
operators, and b) Marchenko double-focused source ($G^{+,+}_{Mdf}$). The
maximum and the minimum value of the grey-level scale of both figures are the
same and a linear time-varying ($\frac{t}{dt}$) gain is applied to both
records to amplify weaker events.

In Figure 6 the predicted data of both approaches after 20 iterations are
shown. The predicted data of the conventional double focusing is not a good
fit for the conventional double-focused data. In contrast, the predicted data
of the Marchenko double-focusing approach is a better fit for Marchenko
double-focused data. In Figure 6(b), the weak multiple reflections generated
by interactions between target and overburden can be seen inside the red box.
This confirms that our method is not only able to suppress multiple
reflections generated inside the overburden significantly but also able to
predict the interactions between the target and the overburden.

(a)

(b)

Figure 6: Predicted data with a virtual source located at the middle of the
model with a single-layered overburden after 20 iterations. a) Conventional
double-focusing, and b) Marchenko double-focusing. Weak multiple reflections
that are predicted by our method can be seen inside the red box. The maximum
and the minimum value of the grey-level scale of both figures are the same and
a linear time-varying ($\frac{t}{dt}$) gain is applied to both records to
amplify weaker events.

Figures 7 and 8 show the imaging result of these two approaches for the first
iteration and after 20 iterations for the model with single-layered
overburden. Note the horizontal event in Figure 6, coming from the overburden
multiple, which overlaps the image of the target. In figure 7 this horizontal
event is strongly suppressed. To demonstrate the resolution improvement of our
method, we compare the horizontal and vertical sections of the image of
different approaches in Figures 9(a) and 9(b). As these figures are showing,
our approach has a superior resolution. A comparison of the convergence rate
of the cost functions of these approaches is also shown in Figure 10.

(a) RTM

(b) LSRTM

Figure 7: Retrieved target image of the model with a single-layered overburden
by conventional double-focusing operators. a) RTM, and b) LSRTM after 20
iterations.

(a) RTM

(b) LSRTM

Figure 8: Retrieved target image of the model with a single-layered overburden
by Marchenko double-focusing. a) RTM, and b) LSRTM after 20 iterations.

(a) Horizontal cross-section

(b) Vertical cross-section

Figure 9: A Comparison of the Horizontal and vertical cross-section of the
retrieved image of the simple model with different approaches. a) A horizontal
cross-section at depth of 450 $m$, and b) A vertical cross section at lateral
distance zero. Figure 10: A comparison between convergence rate of cost
functions of different approaches to TOLSRM for the simple model.

####  3.1.2 Syncline model

As a second example, we design a model with more layers and a syncline-shaped
layer in its overburden (Fig. 11(a)) and a line reflector with a width of 5
meters in the target region. In this model, the background velocity of the
target area is 2400 $m/s$, and the velocity of the reflector is 2700 $m/s$.
Moreover, density variations are introduced in the overburden (Fig. 11(b)) to
generate substantial internal multiples, but the density of the target region
is constant and equal to 1000 $kg/m^{3}$. For this model, 301 sources and
receivers with the same spacing as in the first model are deployed. The number
and spacing of real and virtual sources and receivers are the same. In Figures
11(a), and 11(b) blue crosses are the actual source and receiver locations,
and the red dots are the virtual source and receivers locations.

(a) Velocity model

(b) Density model

Figure 11: Velocity and density of the model with the syncline overburden and
a line reflector in the target.

Figure 12 shows the middle common-source record at the acquisition surface for
the model with the syncline overburden.

Figure 12: Common-source record at the surface of a source located at lateral
distance = 0 in the model with the syncline overburden. A linear time-varying
($\frac{t}{dt}$) gain is applied to this record to amplify weaker events.

In Figure 13, we show the conventional double-focused data (Fig. 13(a)), the
predicted data in case of conventional double-focusing (Fig. 13(c)), the
Marchenko double-focused data (Fig. 13(b)) and the predicted data in case of
Marchenko double-focusing (Fig. 13(d)) for the model with the syncline
overburden. In Figure 13(b), it can be seen, compared to conventional double-
focused data (Fig. 13(a)), that the multiples purely generated in overburden
are suppressed by applying Marchenko double focusing. In addition, by
comparing Figure 13(c) and Figure 13(d) we see that the predicted data of
Marchenko double-focusing not only predicts the primary reflection, but also
predicts the multiples that are generated by interactions between the target
and the overburden.

(a)

(b)

(c)

(d)

Figure 13: A comparison between Double-focused data and predicted data of the
model with the syncline overburden corresponding to the middle virtual source.
a) Conventional double-focused data, b) Marchenko double-focused data, c)
predicted data of conventional double-focusing, and d) predicted data of
Marchenko double-focusing. The first few time samples are set to zero to mute
focusing artifacts. The maximum and the minimum value of the grey-level scale
of all figures are the same.

Moreover, in Figure 14, we show the double-focused down-going Green’s function
(Eq. 12). We see that the virtual source wavefield for the conventional
double-focused predicted data contains only a bandlimited delta function,
whereas the virtual source wavefield for the Marchenko double-focused
predicted data contains additional events to predict the target and overburden
interactions in the Marchenko double-focused data in Figure 13(b).

(a)

(b)

Figure 14: Virtual source wavefield of the model with the syncline overburden
for the middle virtual source. a) The conventional double-focused source
($G^{+,+}_{Cdf}$) by applying conventional redatuming operators, and b)
Marchenko double-focused source ($G^{+,+}_{Mdf}$). The maximum and the minimum
value of the grey-level scale of both figures are the same.

Consequently, the conventional double-focusing approach is unable to fit the
predicted data, whereas the Marchenko double-focusing approach can produce
predicted data with an acceptable fit to the Marchenko double-focused data.
Note that the predicted data of our proposed method can predict the
interactions between the target and the overburden, without the knowledge of
the overburden reflectivity model, by taking $G^{+,+}_{Mdf}$ as the virtual
source wavefield. This source contains these interactions, as can be seen in
Figure 14(b).

Figures 15 and 16 show the imaging result of these two approaches for the
first iteration and after 20 iterations for the model with the syncline
overburden. To demonstrate the improvement of vertical resolution and lateral
continuity, we compare the horizontal and vertical sections of the image of
different approaches in Figures 17(a) and 17(b). As these figures are showing,
our approach has a superior vertical resolution and lateral continuity. A
comparison of the cost functions of these different approaches to TOLSRTM is
also shown in Figure 18. TOLSRTM converges faster when Mdf (instead of Cdf) is
used for redatuming, since the multiple reflections are correctly taken into
account.

(a) RTM

(b) LSRTM

Figure 15: Retrieved target image of the model with the syncline overburden by
conventional double-focusing operators. a) RTM, and b) LSRTM after 20
iterations

(a) RTM

(b) LSRTM

Figure 16: Retrieved target image of the model with the syncline overburden by
Marchenko double-focusing. a) RTM, and b) LSRTM after 20 iterations.

(a) Horizontal cross-section

(b) Vertical cross-section

Figure 17: A Comparison of the Horizontal and vertical cross-section of the
retrieved image of the model with the syncline overburden with different
approaches. a) A horizontal cross-section at depth of 1200 $m$, and b) A
vertical cross-section at lateral distance zero. Figure 18: A comparison
between convergence rate of cost functions of different approaches to TOLSRM
for the model with the syncline overburden.

###  3.2 Investigating the effect of the point-spread function

In this section, we investigate the effect of the PSF and the necessity of
applying it to the reflection response of the target area. In Figure 19, we
show the effect of the PSF (Eq. 16) on the true reflection response (as
modeled during each iteration of TOLSRTM) of the target area of the simple
model and compare it with Marchenko double-focused data. Figure 19(a) shows
the modeled data (Eq. 28) in the true perturbation model without PSF. Figure
19(b) shows the same after the PSF has been applied. Figure 19(c) shows the
Marchenko double-focused data. We can see the effect of the PSF by analyzing
Figures 19(b) and 19(c). The PSF not only imposes a band limitation on the
data but also dampens the far-offsets such that the predicted data fits the
observed data better. Thus, applying the PSF to the reflection response is
necessary for better convergence. Moreover, in Figure 20 the cost function of
TOLSRTM with and without PSF is shown. From this Figure, we conclude that
including the PSF in our formulation successfully improves the convergence.

(a)

(b)

(c)

Figure 19: Comparison of the reflection response of target area of the model
with a single-layered overburden. a) Without PSF, and b) With PSF, and c)
double-focused data. The maximum and the minimum value of the grey-level scale
of figures are set to the maximum and the minimum of the input. Figure 20:
Comparison of the cost function of the TOLSRTM with and without PSF.

###  3.3 Inaccurate velocity model

To determine the sensitivity of our algorithm to kinematic errors, we define a
model similar to the model with a syncline in the overburden but with a high-
velocity inclusion in the overburden and a more complex target (Fig. 21(a)).
We use a smooth version of this model (Fig. 21(b)) to estimate the direct
arrival between the surface and the focusing level. In addition, we use the
target part of this smooth velocity model to compute the background Green’s
functions with a finite difference algorithm for the forward and adjoint
kernels. The density model is shown in Figure 21(c). In Figure 22 we show
images as obtained by applying our methodology with the inaccurate velocity
model. The yellow arrows and ellipses in Figure 22 indicate discontinuities
and depth estimation errors caused by velocity inaccuracies. Moreover, the red
arrows in this figure designate a couple of artifacts that are removed by our
method, and the white arrows and ellipses show the amplitude improvements
resulting from our method. Despite the fact that we use a wrong velocity
model, our algorithm which is based on Marchenko double-focusing is still able
to correctly predict and reduce the artifacts caused by the overburden-related
multiple reflections and improves the quality and amplitude recovery of the
LSRTM image. The observation that Marchenko imaging is not sensitive to errors
in the velocity model is not new, see for instance Broggini et al. Bro.

(a)

(b)

(c)

Figure 21: Velocity model: a) True model, b) smooth model, and c) density
model

(a)

(b)

Figure 22: Target-oriented LSRTM with the Erroneous velocity model. a)
Conventional double-focusing, and b) Marchenko double-focusing. The red arrows
show some overburden artifacts that are removed, the white arrows and ellipses
show amplitude improvement and the yellow ones show some of the effects of
velocity errors. The value of the grey-level scale is set to 25% of the
maximum value.

##  4 Discussion

Our numerical results suggest that our Marchenko target-oriented LSRTM can
recover a clearer image of the target of interest than conventional TOLSRTM.
The retrieved image of the model with single-layered overburden using
Marchenko double-focusing (Fig. 8(b)) has no (or at least a reduced) imprint
of the overburden-related spurious reflector. In contrast, the recovered image
using conventional redatuming operators (Fig. 7(b)) could not handle this
spurious reflector even after 20 iterations. A further investigation into the
resolution improvement of our LSRTM algorithm compared to the conventional
double-focused LSRTM, shows that our algorithm is able to improve the
resolution in both horizontal and vertical directions (Fig. 9). The model with
the syncline overburden demonstrates the power of using Marchenko double-
focusing for target-oriented LSRTM in removing the strong overburden-related
multiple reflection artifacts. The image that is obtained with the
conventional redatuming operator (Fig. 15(b)) is overwhelmed by artifacts from
the overburden. On the other hand, our algorithm is able to remove these
artifacts and retrieves a clear image of the target area (Fig. 16(b)). In
addition, our TOLSRTM algorithm can model data that can predict the
interactions between the target area and the overburden without the knowledge
of the overburden reflectivity model. Moreover, the Marchenko equation has
been derived for elastic media too [Reinicke et al. 2020, Wapenaar & Slob
2014, da Costa Filho et al. 2014]. Hence, we expect that the current
methodology can also be extended to elastic LSRTM.

Furthermore, a comparison between Figure 4(b) and Figure 19 shows the
necessity of applying the PSF to the modeled reflection data of the target.
Without imposing this PSF on the modeled data, the predicted data (Fig. 19(a))
do not match the observed data (Fig. 19(c)). Hence, the TOLSRTM without PSF
converges slowly compared to the case where the PSF is applied to it (Fig.
20).

In our last example, we addressed the issue of velocity inaccuracies in the
overburden. As Figure 22 shows our method is able to improve the quality of
the image by suppressing the overburden-generated multiples and predicting the
interactions between target and overburden, even if the overburden velocity
model is inaccurate.

The Marchenko redatuming algorithm used in this paper requires surface-related
multiples to be removed from the acquired reflection data. However,
researchers exploited the possibility of including the surface-related
multiple reflections within the Marchenko framework [Dukalski & de Vos 2018,
Ravasi 2017, Singh et al. 2017]. Consequently, it is quite straightforward to
include surface-related multiples in our algorithm.

##  5 Conclusion

We have introduced a target-oriented LSRTM algorithm that can correctly handle
the internal multiple reflections generated in the overburden and interactions
between target and overburden. To achieve this, we solved the Marchenko
equations to obtain the Marchenko double-focused data and the downgoing part
of the focusing function. By having these Green’s functions and focusing
operators, we can formulate a Born integral that can be used as the forward
modeling operator and an adjoint modeling operator can also be constructed for
target-oriented LSRTM. Importantly, we have to apply a point-spread function
to the reflection response of the target to account for the finite spatial
bandwidth caused by the overburden and the finite recording aperture at the
acquisition surface.

Since LSRTM and full waveform inversion are closely related, the question may
arise: ”Can we modify our algorithm for target-oriented full waveform
inversion?”. To address this question, we have to mention that the double-
focused data does not include the diving waves, so it is dominated by the
medium-to-short wavelengths and does not include long wavelengths. These long
wavelengths are crucial for velocity recovery in the full waveform inversion
process. Moreover, the evanescent waves, such as refracted waves, are not
handled by the current redatuming methodologies. The aforementioned issues put
a limit on modifying our algorithm to be applied for target-oriented full
waveform inversion, which requires long-to-medium wavelengths. However,
recently some advancement has been made towards including evanescent waves
within the Marchenko framework [Wapenaar et al. 2021, Wapenaar 2020, Diekmann
& Vasconcelos 2021]. Nevertheless, the double-focused data contains all target
reflections correctly (but not the other waveforms). Therefore, another
possible future direction of this research is to investigate the potential of
our method in improving inversion methods that rely on reflections only, such
as joint migration inversion [Verschuur et al. 2016] and reflection waveform
inversion [Brossier et al. 2015].

In recent years, the focus of the seismic exploration community has shifted
toward fast and high-resolution imaging and inversion methods. Our proposed
target-oriented LSRTM is able to produce such images in a relatively small
target region and has a relatively faster convergence rate while correctly
accounting for multiple reflections caused by the overburden.

###### Acknowledgements.

We thank the reviewers, Katrin Löer and Ole Edvard Aaker, for their
constructive comments, which helped us to improve the paper. This work has
received funding from the European Union’s research and innovation programme:
European Research Council (grant agreement 742703).

This study uses numerical data only. The codes used will be shared on
reasonable request to the corresponding author.

## References

  * [Baysal et al. 1983] Baysal, E., Kosloff, D.D., and Sherwood, J.W.C. 1983. reverse-time migration, GEOPHYSICS, 48: 1514-1524.
  * [Berryhill 1984] Berryhill, J.R. 1984. Wave‐equation datuming before stack, GEOPHYSICS, 49: 2064-2066.
  * [Born & Wolf 1999] Born, M., and Wolf, E. 1999. Principles of Optics: Electromagnetic Theory of Propagation, Interference, and Diffraction of Light, Cambridge University Press, London.
  * [Broggini et al. 2014] Broggini, F., Snieder, R., and Wapenaar, K. 2014. Data-driven wavefield focusing and imaging with multidimensional deconvolution: Numerical examples for reflection data with internal multiples, GEOPHYSICS 79: WA107-WA115.
  * [Brossier et al. 2015] Brossier, R., Operto, S., and Virieux, J. 2015. Velocity model building from seismic reflection data by full-waveform inversion, Geophysical Prospecting 63: 354-367.
  * [Barrera et al. 2021] Barrera, D., Schleicher, J., and Brakenhoff, J. 2021. Interferometric redatuming by deconvolution and correlation-based focusing, GEOPHYSICS 86: Q1-Q13.
  * [Cui et al. 2020] Cui, T., James, R., Vasconcelos, I, and Veitch, B. 2020. Target-oriented full-waveform inversion using Marchenko redatumed wavefields, Geophys. J. Int., Volume 223: 792–810.
  * [da Costa Filho et al. 2014] da Costa Filho, C.A., Ravasi, M., Curtis, M., and Meles, G.A. 2014. Elastodynamic Green’s function retrieval through single-sided Marchenko inverse scattering, Physical Review E., 90: 063201.
  * [Dai et al. 2012] Dai, W., Fowler, P., and Schuster G.T. 2012. Multi-source least-squares reverse-time migration Geophysical Prospecting 60: 681-695.
  * [Diekmann et al. 2021] Diekmann, L., Vasconcelos, I., Cummings, D., and Curtis, A. 2021. Towards exact linearized full-waveform inversion via Marchenko redatuming, SEG Technical Program Expanded Abstracts: 3380-3384.
  * [Diekmann & Vasconcelos 2021] Diekmann, L., and Vasconcelos, I. 2021. Focusing and Green’s function retrieval in three-dimensional inverse scattering revisited: A single-sided Marchenko integral for the full wave field, Phys. Rev. Research, 3: 013206.
  * [Dukalski & de Vos 2018] Dukalski, M., and de Vos, K. 2018. Marchenko inversion in a strong scattering regime including surface-related multiples, Geophys. J. Int., volume 212: 760–776.
  * [Dutta et al. 2017] Dutta, G., Giboli, M., Agut, C., Williamson, P., and Schuster G.T. 2017. Least-squares reverse-time migration with local Radon-based preconditioning, GEOPHYSICS 82: S75-S84.
  * [Farshad & Chauris 2021] Farshad, M., and Chauris, H. 2021. Sparsity-promoting multiparameter pseudoinverse Born inversion inacoustic media, GEOPHYSICS 86: S205-S220.
  * [Guo & Alkhalifah 2019] Guo, Q., and Alkhalifah, T. 2019. Datum-based waveform inversion using a subsurface-scattering imaging condition, GEOPHYSICS 84: S251-S266.
  * [Guo & Alkhalifah 2020] Guo, Q., and Alkhalifah, T. 2020. Target-oriented waveform redatuming and high-resolution inversion: Role of the overburden, GEOPHYSICS 85: R525-R536.
  * [Haffinger et al. 2013] Haffinger, P., Gisolf, A., and van den Berg, P. 2013. Towards high resolution quantitative subsurface models by full waveform inversion, Geophys. J. Int., Volume 193: 788–797.
  * [Herrmann & Li 2012] Herrmann, F.J., and Li, X. 2012. Efficient least-squares imaging with sparsity promotion and compressive sensing, Geophysical Prospecting 60:696-712.
  * [Liu et al. 2016] Liu, Y., Liu, X., Osen, A., Shao, Y., Hu, H., and Zheng, Y. 2016. Least-squares reverse-time migration using controlled-order multiple reflections, GEOPHYSICS 83: S347-S357.
  * [Luiken & van Leeuwen 2020] Luiken, N., and van Leeuwen, T. 2020. Seismic wavefield redatuming with regularized multi-dimensional deconvolution, Inverse Problems 36: 095010.
  * [Marquardt 1963] Marquardt, D.W. 1963. An Algorithm for Least-Squares Estimation of Nonlinear Parameters, J. Soc. Indust. Appl. Math. vol.11, No.2.
  * [Morse & Feshbach 1953] Morse, P.M., Feshbach, H. 1953. Methods of theoretical physics, McGraw-Hill Book Company, INC.
  * [Ravasi 2017] Ravasi, M. 2017. Rayleigh-Marchenko redatuming for target-oriented, true-amplitude imaging, GEOPHYSICS 82: S439-S452.
  * [Ravasi et al. 2016] Ravasi, M., Vasconcelos, I., Kritski, A., Curtis, A., da Costa Filho, C.A., and Meles, G.A. 2016. Target-oriented Marchenko imaging of a North Sea field, Geophys. J. Int., Volume 205: 99–104.
  * [Reinicke et al. 2020] Reinicke, C., Dukalski, M., and Wapenaar, K. 2020. Comparison of monotonicity challenges encountered by the inverse scattering series and the Marchenko demultiple method for elastic waves, GEOPHYSICS 85: Q11-Q26.
  * [Schuster 2017] Schuster, G.T. 2017. Seismic inversion, Society of Exploration Geophysicists.
  * [Schuster & Zhou 2006] Schuster, G., and Zhou, M. 2006. A theoretical overview of model-based and correlation-based redatuming methods, GEOPHYSICS 71: SI103-SI110.
  * [Shoja et al. 2020] Shoja, S.M.A., Meles, G.A. and Wapenaar, K. 2020. A Proposal for Marchenko-Based Target-Oriented Full Waveform Inversion, Conference Proceedings, EAGE 2020 Annual Conference & Exhibition: 1-5.
  * [Singh et al. 2017] Singh, S., Snieder, R., van der Neut, J., Thorbecke, J., Slob, E., and Wapenaar, K. 2017. Accounting for free-surface multiples in Marchenko imaging, GEOPHYSICS 82: R19-R30.
  * [Staring et al. 2018] Staring, M., Pereira, R., Douma, H., van der Neut, J., and Wapenaar, K. 2018. Source-receiver Marchenko redatuming on field data using an adaptive double-focusing method, GEOPHYSICS 83: S579-S590.
  * [Tang 2009] Tang, Y. 2009 Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian, GEOPHYSICS 74: WCA95-WCA107.
  * [Thorbecke & Draganov 2011] Thorbecke, J.W., and Draganov, D. 2011. Finite-difference modeling experiments for seismic interferometry, GEOPHYSICS 76: H1-H18.
  * [Thorbecke et al. 2017] Thorbecke, J.W., Slob, E., Brackenhoff, J., van der Neut, J., and Wapenaar, K. 2017. Implementation of the Marchenko method, GEOPHYSICS 82: WB29-WB45.
  * [Valenciano et al. 2006] Valenciano, A.A., Biondi, B., and Guitton, A. 2006. Target-oriented wave-equation inversion, GEOPHYSICS 71: A35-A38.
  * [van den Berg 2021] van den Berg, P.M. 2021. Forward and Inverse Scattering Algorithms based on Contrast Source Integral Equations, John Wiley & Sons, Inc.
  * [Vargas et al. 2021] Vargas, D., Vasconcelos, I., Sripanich, Y., and Ravasi, M. 2021. Scattering-based focusing for imaging in highly complex media from band-limited, multicomponent data, GEOPHYSICS 86: WC141-WC157.
  * [Verschuur et al. 2016] Verschuur, E., Staal, X.R., and Berkhout, A.J. 2016. Joint migration inversion: Simultaneous determination of velocity fields and depth images using all orders of scattering, The Leading Edge 35: 1037-1046.
  * [Wapenaar et al. 2000] Wapenaar, K., Fokkema, J., Dillen, M., Scherpenhuijsen, P. 2000. One‐way acoustic reciprocity and its applications in multiple elimination and time‐lapse seismics, SEG Technical Program Expanded Abstracts: 2377-2380.
  * [Wapenaar 2020] Wapenaar, K. 2020. The Marchenko method for evanescent waves, Geophys. J. int., Volume 223: 1412–1417.
  * [Wapenaar & Slob 2014] Wapenaar, K., Slob, E. 2014. On the Marchenko equation for multicomponent single-sided reflection data, Geophys. J. Int., 199: 1367-1371.
  * [Wapenaar et al. 2021] Wapenaar, K., Snieder, R., de Ridder, S., and Slob, E., 2021. Green’s function representations for Marchenko imaging without up/down decomposition, Geophys. J. Int., 227: 184–203.
  * [Wapenaar & staring 2018] Wapenaar, K., and Staring, M. 2018. Marchenko-based target replacement, accounting for all orders of multiple reflections, Journal of Geophysical Research: Solid Earth, 123: 4942– 4964.
  * [Wapenaar et al. 2014] Wapenaar, K., Thorbecke, J., van der Neut, J., Broggini, F., Slob, E., and Snieder, R. 2014. Marchenko imaging, GEOPHYSICS 79: WA39-WA57.
  * [Willemsen et al. 2016] Willemsen, B., Malcolm, A., and Lewis, W. 2016. A numerically exact local solver applied to salt boundary inversion in seismic full-waveform inversion, Geophys. J. Int., Volume 204: 1703–1720.
  * [Yuan et al. 2017] Yuan, S., Fuji, N., Singh, S., and Borisov, D. 2017. Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies, Geophys. J. Int., Volume 209: 1699–1717.
  * [Zhao & Sen 2018] Zhao, Z., and Sen, M.K. 2018. Fast image-domain target-oriented least-squares reverse-time migration, GEOPHYSICS 83: A81-A86.
  * [Zhang et al. 2018] Zhang, L., Slob, E., van der Neut, J., and Wapenaar, K. 2018. Artifact-free reverse time migration, GEOPHYSICS 83: A65-A68.
  * [Zhou et al. 2018] Zhou, H.W., Hu, H., Zou, Z., Wo, Y., and Young, O. 2018. Reverse time migration: A prospect of seismic imaging methodology, Earth-science Reviews, 179: 207-227.