Target-Enclosed Least-Squares Seismic Imaging

TL;DR

This paper introduces a target-enclosed LSRTM method that leverages acoustic reciprocity and Marchenko redatuming to efficiently produce high-resolution seismic images of specific subsurface regions.

Contribution

It develops a novel target-enclosed LSRTM algorithm incorporating wavefield representations and explores the use of virtual receivers for improved imaging efficiency.

Findings

Effective high-resolution imaging of targeted regions.

Reduced computational costs compared to full-volume LSRTM.

Potential of Marchenko redatuming for virtual receiver implementation.

Abstract

Least-Squares Reverse-Time Migration (LSRTM) is a method that seismologists utilize to compute a high-resolution subsurface image. Nevertheless, LSRTM is a computationally demanding problem. One way to reduce the computational costs of the LSRTM is to choose a small region of interest and compute the image of that region. However, finding representations that account for the wavefields entering the target region from the surrounding boundaries is necessary. This paper confines the region of interest between two boundaries above and below this region. The acoustic reciprocity theorem is employed to derive representations for the wavefields at the upper and lower boundaries of the target region. With the help of these representations, a target-enclosed LSRTM algorithm is developed to compute a high-resolution image of the region of interest. Moreover, the possibility of using virtual receivers created by Marchenko redatuming is investigated.

Equations

$$P^{scat}({\text{x}}^{\prime })={\int }_{\mathcal{V}}{\gamma }_0^2(\text{x})G_0({\text{x}}^{\prime },\text{x}){\chi }^c(\text{x})P^{inc}(\text{x})d\text{x}.$$

$${\text{P}}_{pred}^{scat}(\delta \text{m})=\text{L}\delta \text{m}.$$

$$\delta {\text{m}}^{img}={\text{L}}^{†}{\text{P}}_{obs}^{scat},$$

$$J(\delta \text{m})=\frac{1}{2}\Vert {\text{P}}_{pred}^{scat}(\delta \text{m})-{\text{P}}_{obs}^{scat}{\Vert }_2^2.$$

$$\begin{gathered} {\int }_{\mathfrak{d}{\mathcal{V}}_u}{\rho }^{-1}(p_A^{+}({\partial }_3{p}_B^{-})+p_A^{-}({\partial }_3{p}_B^{+}))d\text{x}= \\ {\int }_{\mathfrak{d}{\mathcal{V}}_l}{\rho }^{-1}(p_A^{+}({\partial }_3{p}_B^{-})+p_A^{-}({\partial }_3{p}_B^{+}))d\text{x}. \end{gathered}$$

$${\nabla }^{2}G+k(\text{x}{)}^{2}G=-\rho \delta (\text{x}-{\text{x}}_s),$$

$$\underset{x_{u,3}\to x_{u,3^{+}}}{\lim }{\partial }_3{G}_{tar}^{+}(\text{x},{\text{x}}_u)=-\frac{\rho }{2}\delta ({\text{x}}_{\mathrm{H}}-{\text{x}}_{u,\mathrm{H}}).$$

$$\begin{gathered} G^{-}({\text{x}}_u,{\text{x}}_s)={\int }_{\mathfrak{d}{\mathcal{V}}_u}{G}^{+}(\text{x},{\text{x}}_s)\frac{2{\partial }_3}{\rho (\text{x})}{G}_{tar}^{-}(\text{x},{\text{x}}_u)d\text{x} \\ +{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}^{-}(\text{x},{\text{x}}_s)\frac{2{\partial }_3}{-\rho (\text{x})}{G}_{tar}^{+}(\text{x},{\text{x}}_u)d\text{x}. \end{gathered}$$

$$\begin{array}{r}{\int }_{\mathfrak{d}\mathcal{V}}{p}_A^{\pm }(\text{x}){\partial }_3{p}_B^{\mp }(\text{x})d\text{x}=-{\int }_{\mathfrak{d}\mathcal{V}}{p}_B^{\mp }(\text{x}){\partial }_3{p}_A^{\pm }(\text{x})d\text{x}\end{array}$$

$$G_{tar}^{-}({\text{x}}_u^{\prime },{\text{x}}_u)=G_{tar}^{-}({\text{x}}_u,{\text{x}}_u^{\prime })$$

$$G_{tar}^{+}({\text{x}}_l,{\text{x}}_u)=G_{tar}^{-}({\text{x}}_u,{\text{x}}_l),$$

$$\begin{gathered} G^{-}({\text{x}}_u,{\text{x}}_s)W(\omega )={\int }_{\mathfrak{d}{\mathcal{V}}_u}{G}_{tar}^{-}({\text{x}}_u,{\text{x}}_u^{\prime })S_u({\text{x}}_u^{\prime },{\text{x}}_s)d{\text{x}}_u^{\prime } \\ +{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{tar}^{-}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l, \end{gathered}$$

$$S_u({\text{x}}_u^{\prime },{\text{x}}_s)=\frac{2{\partial }_{3,u}^{\prime }}{-\rho ({\text{x}}_u^{\prime })}{G}^{+}({\text{x}}_u^{\prime },{\text{x}}_s)W(\omega )$$

$$S_l({\text{x}}_l,{\text{x}}_s)=\frac{2{\partial }_{3,l}}{\rho ({\text{x}}_l)}{G}^{-}({\text{x}}_l,{\text{x}}_s)W(\omega )$$

$$\begin{gathered} {\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{tar}^{-}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l= \\ {\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{0,tar}^{-}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l \\ +{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l. \end{gathered}$$

$$\begin{gathered} G^{-}({\text{x}}_u,{\text{x}}_s)W(\omega )-{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{0,tar}^{-}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l= \\ {\int }_{\mathfrak{d}{\mathcal{V}}_u}{G}_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_u^{\prime })S_u({\text{x}}_u^{\prime },{\text{x}}_s)d{\text{x}}_u^{\prime } \\ +{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l. \end{gathered}$$

$$\begin{gathered} P_{obs}^{scat,TE}({\text{x}}_u,{\text{x}}_s)= \\ {G}^{-}({\text{x}}_u,{\text{x}}_s)W(\omega )-{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{0,tar}^{-}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l, \end{gathered}$$

$$\begin{gathered} G_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_u^{\prime })= \\ {\int }_{\mathcal{V}}{\gamma }_0^2(\text{x})G_{0,tar}({\text{x}}_u,\text{x}){\chi }^c(\text{x})G_{0,tar}(\text{x},{\text{x}}_u^{\prime })d\text{x}, \end{gathered}$$

$$\begin{gathered} G_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_l)= \\ {\int }_{\mathcal{V}}{\gamma }_0^2(\text{x})G_{0,tar}({\text{x}}_u,\text{x}){\chi }^c(\text{x})G_{0,tar}(\text{x},{\text{x}}_l)d\text{x}, \end{gathered}$$

$$\begin{gathered} {\int }_{\mathfrak{d}{\mathcal{V}}_u}{G}_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_u^{\prime })S_u({\text{x}}_u^{\prime },{\text{x}}_s)d{\text{x}}_u^{\prime } \\ +{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{tar}^{-,scat}({\text{x}}_u,{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l= \\ {\int }_{\mathcal{V}}{\gamma }_0^2(\text{x})G_{0,tar}({\text{x}}_u,\text{x}){\chi }^c(\text{x})P_{inc,TE}(\text{x},{\text{x}}_s)d\text{x}, \end{gathered}$$

$$\begin{gathered} P^{inc,TE}(\text{x},{\text{x}}_s)= \\ {\int }_{\mathfrak{d}{\mathcal{V}}_u}{G}_{0,tar}(\text{x},{\text{x}}_u^{\prime })S_u({\text{x}}_u^{\prime },{\text{x}}_s)d{\text{x}}_u^{\prime } \\ +{\int }_{\mathfrak{d}{\mathcal{V}}_l}{G}_{0,tar}(\text{x},{\text{x}}_l)S_l({\text{x}}_l,{\text{x}}_s)d{\text{x}}_l. \end{gathered}$$

$$\begin{gathered} P_{pred}^{scat,TE}({\text{x}}_u,{\text{x}}_s)= \\ {\int }_{\mathcal{V}}{\gamma }_0^2(\text{x})G_{0,tar}({\text{x}}_u,\text{x}){\chi }^c(\text{x})P^{inc,TE}(\text{x},{\text{x}}_s)d\text{x}. \end{gathered}$$

$$J(\delta \text{m})=\frac{1}{2}\Vert {\text{P}}_{pred}^{scat,TE}(\delta \text{m})-{\text{P}}_{obs}^{scat,TE}{\Vert }_2^2.$$

$$\begin{gathered} G_{Mar}^{-}({\text{x}}_v,{\text{x}}_s)={\int }_{\mathfrak{d}{\mathcal{V}}_0}R({\text{x}}_s,{\text{x}}_s^{\prime })f_1^{+}({\text{x}}_s^{\prime },{\text{x}}_v)d{\text{x}}_s^{\prime } \\ -f_1^{-}({\text{x}}_s,{\text{x}}_v), \end{gathered}$$

$$\begin{gathered} G_{Mar}^{+}({\text{x}}_v,{\text{x}}_s)=-{\int }_{\mathfrak{d}{\mathcal{V}}_0}R({\text{x}}_s,{\text{x}}_s^{\prime })f_1^{-}({\text{x}}_s^{\prime },{\text{x}}_v{)}^{\ast }d{\text{x}}_s^{\prime } \\ +f_1^{+}({\text{x}}_s,{\text{x}}_v{)}^{\ast }. \end{gathered}$$

$$R({\text{x}}_s,{\text{x}}_s^{\prime })=\frac{2{\partial }_{3,s}^{\prime }}{-\rho ({\text{x}}_s^{\prime })}{G}^{-}({\text{x}}_s,{\text{x}}_s^{\prime }).$$

$$G^{+}({\text{x}}_u,{\text{x}}_s)\approx G_{Mar}^{+}({\text{x}}_u,{\text{x}}_s),$$

$$G^{-}({\text{x}}_u,{\text{x}}_s)\approx G_{Mar}^{-}({\text{x}}_u,{\text{x}}_s),$$

$$G^{-}({\text{x}}_l,{\text{x}}_s)\approx G_{Mar}^{-}({\text{x}}_l,{\text{x}}_s).$$

Full text

††publicationid: pubid: 1558-0644 © 2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html

I Introduction

Wavefield migration is the art of computing the medium reflectivity from the recorded wavefield passing through the medium. Many migration algorithms are available, such as Kirchhoff migration [1], one-way wave equation migration [2], and Reverse-Time Migration (RTM) [3, 4, 5]. RTM is one of the common migration algorithms. RTM is commonly implemented by applying the adjoint of the Born operator to the recorded data [6]. However, the inverse of the operator is needed for the true image. Consequently, the migration result suffers from amplitude and resolution issues [7, 8, 9]. One way to address this problem is to solve the migration problem with a least-squares solution called Least-Squares Reverse-Time Migration (LSRTM). The least-squares solution is usually applied as an iterative optimization problem.

However, LSRTM is computationally expensive [10, 8, 11, 12]. To reduce the computation cost of the LSRTM, one can reduce the computation domain by confining the model to a Region of Interest (ROI) by recording or computing the wavefields at the boundary of this region. The process of migrating for a ROI is called target-oriented migration [13, 14, 15, 16, 17, 18, 19]. However, deploying receivers on the boundaries of ROI usually is not possible due to physical obstacles.

Due to the aforementioned accessibility issue, the typical approach for target-oriented migration is to opt for redatuming algorithms and only consider the upper horizontal boundary of the ROI or target [20, 18, 21, 22, 23, 24, 25, 26]. In cases where the ROI is enclosed between two boundaries, i.e. when wavefields are entering the ROI from the underburden through the lower boundary [27, 28], the shortcoming of only considering the upper boundary is that any wavefield entering the ROI from the medium below the lower boundary of the ROI is unaccounted for, hindering the convergence of the inversion process. Moreover, including the lower boundary in the algorithm can add transmission information to the inversion. However, including the lower boundary in the inversion process has rarely been studied directly. For instance, Cui et al. [29] derive a representation with a reciprocity theorem and Marchenko redatuming to include surrounding boundaries in target-oriented full waveform inversion (FWI), Diekmann et al. [30] use a Marchenko retrieved Green’s function of the ROI and insert it in the Lippmann-Schwinger integral to create a linear inversion process, and van der Neut et al. [31] design a target-enclosed imaging algorithm with the help of a reciprocity theorem. Of the above-mentioned papers, only [31] directly studies the consequences of including the lower boundary in the imaging process, and the others only implicitly imply the effects of it.

This paper, which is an extension of [31], studies the contribution of the lower boundary by introducing a target-enclosed LSRTM algorithm. To derive this algorithm, we start by explaining the LSRTM briefly. Next, we derive a target-enclosed representation for Green’s functions on the upper and lower boundaries of the ROI by using the reciprocity theorem. Then, we connect this target-enclosed representation to LSRTM to complete our algorithm. After deriving the required equations, we test our algorithm with numerical examples. First, we use physical receivers at the boundaries of the ROI to check the results in the ideal situation. Then, we briefly introduce Marchenko redatuming to explore the possibility of using virtual receivers in our algorithm. For both physical and virtual receivers cases we use both a homogeneous background and a smooth background velocity as the velocity model for migration. Finally, we finish the paper by discussing the results and providing a conclusion.

II Theory

To develop the theory, we start with a brief discussion of LSRTM. Next, a representation of the target-enclosed Green’s function is given with the help of the reciprocity theorem. Finally, to derive our target-enclosed LSRTM formulation, we combine LSRTM with the target-enclosed representations.

In the entire theory section, we are in the frequency-space domain, and for simplicity, we drop the dependency on angular frequency ($\omega$).

II-A Least-Squares Reverse-Time Migration

We start the explanation of LSRTM by investigating the Born integral for the scattered wavefield by a scattering potential [32, 6, 33]. Here we follow the convention of [33]:

[TABLE]

In this equation, $\textbf{x}^{\prime}$ is the observation location, x is a location inside the computation volume ($\mathcal{V}$), $P^{scat}(\textbf{x}^{\prime})$ is the scattered pressure field at the observation point, $P^{inc}(\textbf{x})$ is the incident pressure field at the computation point and $G_{0}(\textbf{x}^{\prime},\textbf{x})$ is the background Green’s function between x and $\textbf{x}^{\prime}$. Moreover, $\gamma_{0}(\textbf{x})=\frac{-i\omega}{c_{0}(\textbf{x})}$ and $\chi^{c}(\textbf{x})=1-\frac{c^{2}_{0}(\textbf{x})}{(c^{scat}(\textbf{x}))^{2}}$ is the propagation velocity perturbation, where $c_{0}(\textbf{x})$ and $c^{scat}(\textbf{x})$ are the background and scatterer’s velocity, respectively.

It is possible to rewrite this equation in matrix form:

[TABLE]

Here, L is matrix form of the integral operator of equation 1 and $\delta\textbf{m}$ is a vector, containing the perturbation $\chi^{c}(\textbf{x})$.

To obtain an estimation of $\delta\textbf{m}$, we can apply the adjoint of L to the observed data:

[TABLE]

where $\dagger$ denotes complex conjugate transposition.

We can go one step further to obtain a high-resolution estimation of the image by minimizing the following objective function [34]:

[TABLE]

Different optimization algorithms, such as conjugate gradient, can minimize this function. This optimization problem is known as Least-Squares Reverse-Time Migration (LSRTM).

II-B Target-enclosed representations

LSRTM is a computationally expensive algorithm. In order to reduce its computational burden, seismologists usually opt for a target-oriented algorithm, limiting the medium to a smaller region. As we mentioned in section I, target-oriented algorithms usually redatum the data to the upper boundary of the target and ignore any information coming from the lower boundary.

In this section, the idea is to find a representation that can account for a heterogeneous medium above the upper and below the lower boundary of the target area. To derive this representation, we follow [31]. The starting point of the derivation is the acoustic reciprocity theorem of the convolution type [21], which connects the wavefields of two different states via

[TABLE]

Here, we consider a volume $\mathcal{V}$, which is limited by two infinite horizontal surfaces. States A and B (Fig. 1) are defined in two different media which are identical inside the volume $\mathcal{V}$ with boundaries denoted by $\mathfrak{d}\mathcal{V}_{u}$ (upper) and $\mathfrak{d}\mathcal{V}_{l}$ (lower), and arbitrary outside of this volume. In addition, $p^{+}$ and $p^{-}$ are decomposed wavefields on the boundaries where $+$ means downgoing and $-$ means upgoing, and $\partial_{3}$ is the partial derivative in direction $x_{3}$ (downward).

To continue the derivation we define Green’s function $G(\textbf{x},\textbf{x}_{s})$ as the solution of the following Helmholtz equation:

[TABLE]

where $k(\textbf{x})=\frac{\omega}{c(\textbf{x})}$ is the wavenumber, and $c$ is the propagation velocity.

A representation for Green’s function at the upper boundary of the target ($\mathfrak{d}\mathcal{V}_{u}$) can be derived by defining state A in the actual medium and state B in a medium identical to medium A inside the volume $\mathcal{V}$ (target) and homogeneous outside of it. We denote the Green’s functions of state A with $G(\textbf{x},\textbf{x}_{s})$ and Green’s function of state B with $G_{tar}(\textbf{x},\textbf{x}_{u})$. In state A we define $p_{A}^{\pm}=G^{\pm}(\textbf{x},\textbf{x}_{s})$, where $\textbf{x}_{s}$ is a location at Earth’s surface $\mathfrak{d}\mathcal{V}_{0}$. For state B we have $p_{B}^{\pm}=G_{tar}^{\pm}(\textbf{x},\textbf{x}_{u})$. Here, $\textbf{x}_{u}$ is located at the upper boundary $\mathfrak{d}\mathcal{V}_{u}$ and $G_{tar}^{\pm}(\textbf{x},\textbf{x}_{u})$ is stimulated by an impulsive point source at $\textbf{x}_{u}$. Since the half-space above the target is homogeneous in state B, the vertical derivative of the $G_{tar}^{+}(\textbf{x},\textbf{x}_{u})$ at $\mathfrak{d}\mathcal{V}_{u}$ is

[TABLE]

Here, ${x_{u,3}\to x_{u,3^{+}}}$ means the limit from below the boundary, where $\textbf{x}_{H}$ and $\textbf{x}_{u,H}$ denote the horizontal coordinates of x and $\textbf{x}_{u}$, respectively. In addition, since in state B the medium below the target is reflection-free, $G^{-}_{tar}(\textbf{x},\textbf{x}_{u})$ and its derivative disappear at the lower boundary $\mathfrak{d}\mathcal{V}_{l}$. By substituting all of the ingredients into equation 5, the following can be reached:

[TABLE]

This equation is the base for our target-enclosed LSRTM derivation. The first integral on the right-hand side of Equation 8 accounts for anything entering the medium from the upper boundary, and the second integral accounts for anything that comes from the lower boundary.

II-C Target-enclosed LSRTM

To merge this representation with LSRTM, we use

[TABLE]

[35], and we use Green’s function reciprocity:

[TABLE]

[TABLE]

where $\textbf{x}^{\prime}_{u}$ is an element of $\mathfrak{d}\mathcal{V}_{u}$ and $\textbf{x}_{l}$ is an element of $\mathfrak{d}\mathcal{V}_{l}$. We rewrite Equation 8 as follows:

[TABLE]

where $W(\omega)$ is the source signature,

[TABLE]

is the dipole source from the upper boundary of the target, which accounts for reflections from above the upper boundary of the target. Further

[TABLE]

is the dipole source term from the lower boundary of the target, which accounts for the reflections generated below the target. Fig. 2 represents the right-hand side of Equation 12.

As said above, the second integral on the right-hand side of equation 12 is the contribution of the medium below the target to the data. This integral can be split into two terms: 1) the arrival from the lower boundary to the upper one in the background, and 2) the forward scatterings inside the target region. In mathematical terms:

[TABLE]

Here, $G^{-}_{0,tar}(\textbf{x}_{u},\textbf{x}_{l})$ is the Green’s function in the background model of the target, and $G^{-,scat}_{tar}(\textbf{x}_{u},\textbf{x}_{l})$ is the Green’s function that contains the scattered events. By substituting equation 15 in equation 12 and taking the background contribution to the left-hand side we end up with:

[TABLE]

Importantly, the Green’s function $G^{-}_{tar}(\textbf{x}_{u},\textbf{x}^{\prime}_{u})$ in the first integral on the right-hand side of equation 12 is the scattered Green’s function inside the target, so we rename it to $G^{-,scat}_{tar}(\textbf{x}_{u},\textbf{x}^{\prime}_{u})$.

To obtain our target-enclosed LSRTM algorithm, we take three steps: First, we assign

[TABLE]

where ”TE” stands for ”Target-Enclosed”. Second we compute the scattered Green’s functions ($G^{-,scat}_{tar}(\textbf{x}_{u},\textbf{x}^{\prime}_{u})$ and $G^{-,scat}_{tar}(\textbf{x}_{u},\textbf{x}_{l})$) with equation 1

[TABLE]

and

[TABLE]

where x is a location inside the target volume. Using this in the right-hand side of Equation 16, we obtain

[TABLE]

where

[TABLE]

Finally, we complete our derivation by assigning

[TABLE]

Our approach computes an incident wavefield, which contains all of the information from the surrounding medium of the target of interest. This means our approach is not limited to one kind of parametrization, and it can be implemented for any parametrization choice. Further, It is also possible to inject $S_{u}$ and $S_{l}$ as dipole sources using a finite-difference algorithm instead of solving Equations 21 and 22. Finally, we can solve the following least-squares problem, i.e., minimizing the objective function:

[TABLE]

II-D Marchenko Green’s function retrieval

We obtained a target-enclosed LSRTM algorithm in the previous section. Nevertheless, in most real-world situations, one does not have physical access to the boundaries of the region of interest. The alternative to physical receivers inside the medium is to create virtual receivers with redatuming. Marchenko redatuming is a state-of-the-art data-driven approach that can compute Green’s functions at any depth level with all orders of multiple reflections from the reflection response at the surface and a smooth background model of the medium.

To summarise, these redatumed Green’s functions are retrieved by iteratively solving the Marchenko-type representations. These representations are [21]

[TABLE]

and

[TABLE]

Here, $\mathfrak{d}\mathcal{V}_{0}$ is the surface, $\textbf{x}_{s}$ and $\textbf{x}^{\prime}_{s}$ are locations at the surface, and $\textbf{x}_{v}$ is a virtual location on an arbitrary depth. Moreover, $f^{\pm}_{1}$ are upgoing (-) and downgoing (+) parts of the focusing function. In addition, $R(\textbf{x}_{s},\textbf{x}^{\prime}_{s})$ is the reflection response at the surface which is related to the upgoing Green’s function of the medium via

[TABLE]

Here, $\partial^{\prime}_{3,s}$ is the vertical partial derivative at $\textbf{x}^{\prime}_{s}$. We refer to [36] for a comprehensive explanation of the derivation and numerical algorithms for solving these equations.

Consequently, we can substitute Marchenko Green’s functions with the target boundaries Green’s functions as follows:

[TABLE]

[TABLE]

and

[TABLE]

Equations 27 and 29 can be used in equations 13 and 14 to obtain the source terms $S_{u}$ and $S_{l}$, whereas equation 28 can be used in equation 17 to obtain the observed target-oriented scattered response.

III Numerical results

III-A * single-sided algorithm vs. double-sided algorithm*

In this section, we aim to visualize the performance of target-enclosed LSRTM. Here we use direct modeling of the Green’s functions with receivers inside the medium; in section B we use the Marchenko method to retrieve these Green’s functions from the reflection response at the surface. A model with dimensions of 1000 $m$ by 650 $m$ is designed as shown in Fig. 3. The spatial grid sampling is 5 $m$ in both directions. The target region consists of a rectangular velocity anomaly embedded in a homogeneous background and a constant density. A total of 201 sources are placed at the surface, and 402 receivers are on the target’s upper (250 $m$) and lower (550 $m$) boundaries. The required Green’s functions and wavefields are computed by a finite-difference algorithm [37] and a Ricker wavelet with a dominant frequency of 30 $Hz$, where the recording time sampling of the receivers is set to 4 $ms$. According to the theory section, the wavefields at the receiver positions are decomposed into upgoing and downgoing components. Fig. 4 shows the upgoing component of the data at the upper boundary with a source located at $\textbf{x}_{s}=(0,0)$, corresponding to the first term on the right-hand side of Equation 17. We study the effects of including the lower boundary with two different cases: 1) a homogeneous background velocity and 2) a smooth background velocity for migration.

III-A1 Homogeneous background velocity

In this section, we show the results using a homogeneous background velocity in the target area. The true perturbation model of this case is shown in Fig. 5. We design two scenarios to demonstrate the performance and consequences of the target-enclosed algorithm. For the first scenario, we only consider the upper boundary of the target in the inversion process and completely ignore the lower boundary contribution. We call this scenario ”single-sided algorithm”. For this scenario, the observed data is the same as Fig. 4. For the second one, we include the lower boundary contribution, which is the ”double-sided algorithm” explained in the previous section. The observed data for this scenario is shown in Fig. 6(a), which corresponds to the left-hand side of Equation 17, whereas Fig. 4 shows $G^{-}(\textbf{x}_{u},\textbf{x}_{s})W(\omega)$. The black arrows in Figures 4 and 6(a) indicate the full reflection from the reflector below the lower boundary and the forward scattered part of it, respectively. Figures 6(b) and 6(c) show the predicted data of single-sided and double-sided algorithms, respectively, after 30 iterations of LSRTM. A detailed investigation of Fig. 6 proves that the double-sided algorithm can predict the forward scattered event that is passing through the perturbation which is indicated by a black arrow in Fig. 6(c).

To move our investigation further, we show the imaging results of both algorithms in Fig. 7. Figures 7(a) and 7(b) show the RTM images of both approaches. As we can see, the RTM result of the double-sided algorithm faintly reveals the long wavelength part of the model. Moving to the LSRTM results in Figures 7(c) and 7(d), we observe an interesting outcome. The LSRTM result of the single-sided algorithm shows that it cannot recover the long wavelength part of the model (Fig. 7(c)). In contrast, the double-sided algorithm can incorporate the information embedded inside the forward scattered field (Fig. 7(d)), and it recovers the long wavelength parts of the volume perturbation. Moreover, Fig. 8 shows the horizontal cross-section of the retrieved perturbation. In this figure, we can see the double-sided approach can recover the vertical boundaries of the perturbation. Nevertheless, since the background velocity for the migration is not updated during LSRTM, the fit of the reflected event from below the target is not accurate.

To conclude this section, we compare the cost functions of both approaches in Fig. 9. The cost function of target-oriented LSRTM shows a slow convergence rate. In comparison, the target-enclosed approach includes the extra information coming from the lower boundary, so its cost function converges faster and to a lower minimum.

III-A2 Smooth background velocity

This section uses a smooth background velocity for migration. The setup is exactly the same as before, except for the background velocity. The perturbation model for this case is shown in Fig. 10. Again, we do the same two scenarios as before, i.e. ”single-sided” and ”double-sided.” For the single-sided scenario, the observed data is the same as before (Fig. 4). However, since the right-hand side of the equation 17 is computed in a different background velocity, the last primary in the observed data for the double-sided approach (Fig. 11(a)) is slightly different.

Comparing the results of both approaches shows (Fig. 11) that the double-sided approach (Fig. 11(c)) can predict the reflected event coming from the lower boundary. Since we use a smooth background velocity here, this prediction is more accurate than the previous section’s results (Fig. 6(c)). In the image domain (Fig. 12), the double-sided approach recovers a faint image of the vertical sides of the rectangular perturbation (Fig. 12(d)) whereas the single-sided results in a more standard image (Fig. 12(c)). Moreover, Fig. 13 shows the horizontal cross-section of the retrieved perturbation. In this figure, we can see the double-sided approach can recover the vertical boundaries of the perturbation. The double-sided image is more comparable to the true perturbation in Fig. 10. Finally, investigating the cost functions (Fig. 14) of these approaches shows that the double-sided approach converges better since it can predict the event coming from below the target.

III-B * Virtual receivers*

Here, we use the same setup as before but replace the Green’s functions at the boundaries with their Marchenko counterparts. In other words, we create virtual receivers with the help of Marchenko redatuming. For a study about the benefits of using Marchenko redatuming instead of a more conventional redatuming algorithm for target-oriented LSRTM we refer to [26]. In this section, we only focus on including the lower boundary by utilizing virtual receivers created by Marchenko redatuming, and instead of using ”double-sided”, we use the ”target-enclosed” term. Moreover, we show the results for both homogeneous and smooth background velocities in this section.

III-B1 Homogeneous background velocity

Fig. 15 shows the data obtained by Marchenko redatuming. In Fig. 16 the observed data calculated by Equation 17 (Fig. 16(a)) and the predicted data after 30 iterations (Fig. 16(b)) are shown. A comparison between Fig. 16(b) and Fig. 6(c) shows that even in the presence of redatuming error, such as limited aperture and lack of certain parts of the wavelength spectrum, our algorithm can predict acceptable data.

In the image domain, Fig. 17 shows the RTM (Fig. 17(a)) and LSRTM (Fig. 17(b)) image resulting from Marchenko redatumed data. Comparing Fig. 17(b) with Fig. 7(d) reveals our target-enclosed algorithm with redatumed data as input can not recover the long wavelength part of the model. This is due to the fact that the direct arrival of the Marchenko-based Green’s function is incorrect since it is computed in the background model. Consequently, the forward-scattered waveforms responsible for the long wavelength updates in Figure 7d can no longer be utilized.

III-B2 Smooth background velocity

In Fig. 18 the observed data calculated by Equation 17 (Fig. 18(a)) and the predicted data after 30 iterations (Fig. 18(b)) are shown. Similar to the homogeneous case, a comparison between Fig. 18 and Fig. 11 shows that our algorithm successfully predicts the data.

In the image domain, Fig. 19 shows the RTM (Fig. 19(a)) and LSRTM (Fig. 19(b)) image resulting from Marchenko redatumed data. Comparing Fig. 19(b) with Fig. 12(d) shows that redatumed data reveals an acceptable perturbation model. However, the faint recovered vertical interfaces are not presented in Fig. 19(b) since forward-scattered waveforms are not processed correctly in the retrieved Marchenko Green’s function at the lower boundary.

III-B3 LSRTM results for the entire medium

To make a fair comparison, we include the results of standard LSRTM for the entire medium with the smooth background model in this section. Figure 20 shows the LSRTM image after 30 iterations. Figure 20(a) shows the image of the entire medium and Figure 20(b) singled out the target area of it.

IV Discussion

In Section II, we develop a theory for target-enclosed LSRTM that can limit the computation domain by confining the target between two boundaries. Equations 17 and 22 enable us to account for any wavefield entering the target region by including the upper and lower boundaries in the inversion process.

Further, in section III-A, several numerical tests are designed to demonstrate the advantage of incorporating the lower boundary in the conventional target-oriented LSRTM. From the data point of view, our double-sided target-enclosed LSRTM, compared to a conventional single-sided target-oriented LSRTM, not only removes the background arrival from the lower boundary to the upper boundary but also can predict the forward scattered field inside the target (Fig. 6, Fig. 11, and Fig. 16). Additionally, a comparison between the resulting images of both algorithms with a homogeneous migration velocity (Fig. 7), shows that single-sided target-oriented LSRTM can not update the long wavelength part of the model. In contrast, the double-sided target-enclosed LSRTM updates the image according to the Born inversion criteria by integrating the forward scattered wavefield information. In the case of a smooth migration velocity, our algorithm recovers a higher-resolution image and a faint estimation of the vertical sides of the rectangular perturbation (Fig. 12).

Moreover, in section III-B, we investigate the possibility of using virtual receivers created by Marchenko redatuming in our target-enclosed algorithm. Fig. 17 shows that in the case of a homogeneous background migration velocity, it is hardly possible to update the long wavelengths with virtual receiver data, and only short wavelength parts of the perturbation are recovered. Moreover, for the smooth background case, our algorithm increases the image’s resolution by updating the short wavelengths in this case.

Ultimately, we show the standard LSRTM image for the entire medium. Comparing Figure 20 with previous cases shows that our method is superior in imaging the target area in any case. Our method adds valuable information by explicitly incorporating transmitted wavefields. Further, with our setup and hardware, the computational time of a single iteration of LSRTM for the entire medium, which has 201 by 131 grid points, is about 45 seconds, and a single iteration of target-enclosed LSRTM, which has 201 by 61 grid points, is about 25 seconds. We do not consider the computational cost of the Marchenko redatuming method here since we only do it once and compared to the total time of LSRTM it is negligible. One disadvantage we can mention is saving the extra redatumed wavefields and focusing functions on the disk. However, relative to the reduction of the memory need by the reduction of the target dimensions, this can be neglected.

V Conclusion

This paper proposes a target-enclosed seismic imaging algorithm that can account for the wavefields entering the target region from the upper and lower boundaries of the region. The three main advantages of this paper’s algorithm are 1) It significantly reduces the computational domain by limiting it to a smaller domain, 2) it removes interactions with the part of the medium above the upper boundary, and 3) it can incorporate the transmission information from the lower boundary to the upper one.

Nevertheless, our algorithm has also some limitations. First, we need access to the lower boundary of the target to deploy receivers at the boundaries of the target. Second, we need a background model of the target that can predict the arrival time from the lower boundary to the upper. It is possible to overcome the first limitation by using virtual seismology methods such as Marchenko redatuming to create virtual receivers around the target region [38, 39, 40, 41, 24, 42, 21, 26] as we showed in a numerical example in section III-B. To address the second limitation, a reformulation of the target-enclosed LSRTM is possible to make it compatible with full waveform inversion to update the background velocity model [43].

The need for high-resolution images is increasing daily, which demands more computational power. This paper’s proposed target-enclosed LSRTM can produce less computationally demanding high-resolution images by focusing on a relatively small target of interest and including all the interactions between this region and the outside environment.

conflict of interest

All authors declare that they have no conflicts of interest.

VI Biography Section