Q-compensated least-squares reverse time migration using low-rank one-step wave extrapolation

Junzhe Sun; Sergey Fomel; Tieyuan Zhu; Jingwei Hu

doi:10.1190/geo2015-0520.1

ABSTRACT

Attenuation of seismic waves needs to be taken into account to improve the accuracy of seismic imaging. In viscoacoustic media, reverse time migration (RTM) can be performed with $Q$ -compensation, which is also known as $Q$ -RTM. Least-squares RTM (LSRTM) has also been shown to be able to compensate for attenuation through linearized inversion. However, seismic attenuation may significantly slow down the convergence rate of the least-squares iterative inversion process without proper preconditioning. We have found that incorporating attenuation compensation into LSRTM can improve the speed of convergence in attenuating media, obtaining high-quality images within the first few iterations. Based on the low-rank one-step seismic modeling operator in viscoacoustic media, we have derived its adjoint operator using nonstationary filtering theory. The proposed forward and adjoint operators can be efficiently applied to propagate viscoacoustic waves and to implement attenuation compensation. Recognizing that, in viscoacoustic media, the wave-equation Hessian may become ill-conditioned, we propose to precondition LSRTM with $Q$ -compensated RTM. Numerical examples showed that the preconditioned $Q$ -LSRTM method has a significantly faster convergence rate than LSRTM and thus is preferable for practical applications.

attenuation, least squares, Q, reverse time migration, wave equation

INTRODUCTION

Seismic attenuation is caused by the effective anelastic properties of the earth (Aki and Richards, 2002; Carcione, 2007) and may lead to poor illumination and misplacement of reflectors in a migration image. To directly compensate for seismic attenuation during reverse time migration (RTM) (Baysal et al., 1983; McMechan, 1983; Whitmore, 1983), Zhang et al. (2010) propose a viscoacoustic wave equation involving a pseudodifferential operator based on the constant- $Q$ model (Kjartansson, 1979) with decoupled effects of amplitude loss and velocity dispersion. Suh et al. (2012) extend the operator to vertically transversely isotropic (VTI) media. Bai et al. (2013) adopt a similar approach for attenuation compensation in RTM, but use a viscoacoustic wave equation without memory variables. Using fractional Laplacians, Zhu and Harris (2014) propose a constant- $Q$ viscoacoustic wave equation with separate terms accounting for amplitude loss and velocity dispersion, which is further applied for $Q$ -compensated RTM using synthetic and field data (Zhu et al., 2014; Zhu and Harris, 2015). Fletcher et al. (2012) and Sun and Zhu (2015) investigate stable approaches for Q-compensation in RTM.

The imaging problem can also be cast as an inverse problem, with the objective of minimizing the $L_{2}$ norm of the difference between recorded data and predicted data (Ronen and Liner, 2000). Such approaches are known as least-squares migration (Nemeth et al., 1999; Tang, 2009; Dai et al., 2011) and more specifically least-squares RTM (LSRTM) in the context of RTM (Wong et al., 2011; Dai et al., 2012; Dai and Schuster, 2013; Liu et al., 2013; Zhang et al., 2013; Sun et al., 2014a; Xue et al., 2014; Hou and Symes, 2015). LSRTM is capable of mitigating imaging artifacts and enhancing subsurface illumination and may have a correlative objective function to relax the amplitude matching requirement (Zhang et al., 2014). Pioneering works of linearized inversion in viscoacoustic and viscoelastic media are done by Ribodetti et al. (1995, 2000) using an asymptotic theory and by Blanch and Symes (1994, 1995) using the wave equation. Recently, Dutta and Schuster (2014) and Sun et al. (2014b) show that LSRTM can be applied for attenuation compensation in viscoacoustic media. Dutta and Schuster (2014) use the standard linear solid model (Robertsson et al., 1994; Blanch et al., 1995), with a simplified stress-strain relation and incorporated a single relaxation mechanism (Blanch and Symes, 1995). Sun et al. (2014b) use the low-rank one-step method to solve the constant- $Q$ wave equation, which allows for an efficient formulation involving fractional Laplacians (Carcione, 2010; Zhu and Harris, 2014). The computational cost of LSRTM depends on the number of iterations, which hinges on the conditioning of the wave-equation Hessian that it tries to invert. In acoustic media, RTM is an efficient approximation to the inverse of reverse time demigration (RTDM), the forward modeling operator, and it provides accurate kinematic information of subsurface structures (Symes, 2008). However, in viscoacoustic media, RTM is a poor approximation to the inverse of RTDM because the wave amplitude suffers from attenuation during forward and backward propagation (Zhu and Harris, 2014; Sun et al., 2015). As a result, the wave-equation Hessian becomes ill-conditioned and iterative LSRTM suffers from a slow convergence rate. To improve the convergence rate of LSRTM in viscoacoustic media, we propose to seek for a proper preconditioner that mitigates the effect of attenuation in the inversion.

To construct LSRTM in viscoacoustic media, we use the low-rank one-step wave extrapolation (Sun et al., 2016) and derive its adjoint operator based on nonstationary linear filtering theory (Margrave, 1998). Sun et al. (2015) successfully apply the low-rank one-step wave extrapolation operator to solve the constant- $Q$ wave equation with fractional Laplacians. To solve the problem of slow convergence of LSRTM in viscoacoustic media, we propose to construct a preconditioned formulation by replacing the viscoacoustic RTM operator, i.e., RTM based on the solution of the viscoacoustic wave equation forward and backward in time, with a better approximate inverse of the RTDM operator, i.e., the $Q$ -compensated RTM or $Q$ -RTM (Zhang et al., 2010; Suh et al., 2012; Bai et al., 2013; Zhu et al., 2014). $Q$ -RTM involves a modeling operator with separate control over amplitude and phase, and is designed to compensate for the amplitude loss along the attenuated wavepaths. As a result, the preconditioned wave-equation Hessian is well conditioned, helping the new framework to quickly converge to the true amplitude solution within only a few iterations. Because the inverted matrix is numerically non-Hermitian, we adopt the generalized minimum residual (GMRES) algorithm (Saad and Schultz, 1986), for iterative inversion. Using a synthetic model, we test the ability of the proposed $Q$ -LSRTM to dramatically enhance image quality at a reasonable cost.

THEORY

Wave extrapolation in viscoacoustic media

We first briefly review the basic theory of low-rank one-step wave extrapolation in viscoacoustic media.

A constant-

Q

model (Kjartansson, 1979) describes an attenuating medium whose quality factor

Q

is constant in frequency (but may vary in space), indicating that the attenuation coefficient is linear in frequency. Zhu and Harris (2014) derive the following approximate constant-

Q

wave equation with decoupled fractional Laplacians for modeling and imaging in viscoacoustic media:

\frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} = \nabla^{2} p + β_{1} {η (- \nabla^{2})^{γ + 1} - \nabla^{2}} p + β_{2} τ \frac{\partial}{\partial t} (- \nabla^{2})^{γ + 1 / 2} p,

(1)

where

η (x) = - c_{0}^{2 γ (x)} ω_{0}^{- 2 γ (x)} \cos (π γ (x)),

(2)

τ (x) = - c_{0}^{2 γ (x) - 1} ω_{0}^{- 2 γ (x)} \sin (π γ (x)),

(3)

c^{2} (x) = c_{0}^{2} (x) \cos^{2} (π γ (x) / 2),

(4)

γ (x) = \arctan (1 / Q (x)) / π .

(5)

Here,

γ

is a dimensionless parameter that ranges between 0 to 1/2;

p (x, t)

is the pressure wavefield, and

c_{0} (x)

is the acoustic velocity model defined at a reference frequency

ω_{0}

⁠. The

β_{1}

and

β_{2}

parameters act like on/off switches that control velocity dispersion and amplitude loss effects, respectively (Zhu and Harris, 2014). For simplicity of notation, in the rest of the paper, the fractional Laplacian operators are denoted as

L = {(- \nabla^{2})}^{γ + 1}

and

H = {(- \nabla^{2})}^{γ + 1 / 2}

⁠.

Setting

β_{1}

and

β_{2}

to one, equation leads to the viscoacoustic dispersion relation with fractional powers of the wavenumber:

\frac{ω^{2}}{c^{2}} = - η {| k |}^{2 γ + 2} - i ω τ {| k |}^{2 γ + 1},

(6)

Solving for

ω

in equation yields

ω = \frac{- i p_{1} + p_{2}}{2},

(7)

where

p_{1} = τ c^{2} {| k |}^{2 γ + 1},

(8)

p_{2} = \sqrt{- τ^{2} c^{4} {| k |}^{4 γ + 2} - 4 η c^{2} {| k |}^{2 γ + 2}} .

(9)

The phase function

ϕ (x, k, Δ t)

that determines the phase shift of the wavefield for propagation in time is then defined as

ϕ (x, k, Δ t) \approx \frac{- i p_{1} + p_{2}}{2} Δ t .

(10)

The one-step wave extrapolation provides an approximate solution to equation by incorporating the phase function defined in equation into the Fourier integral operator (FIO):

p (x, t + Δ t) = \int P (k, t) e^{i k \cdot x + i ϕ (x, k, Δ t)} d k,

(11)

where

P

is the spatial Fourier transform of

p

⁠, and

ċ

indicates the dot product between two vectors. The accuracy of the approximation increases with smaller

Δ t

(Fomel et al., 2013). The adjoint of operator in equation can be expressed as

P (k, t) = \int p (x, t + Δ t) e^{- i k \cdot x - i \bar{ϕ} (x, k, Δ t)} d x,

(12)

where

\overline{ϕ}

denotes the complex conjugate of

ϕ

⁠.

The FIOs introduced in equations 11 and 12 can be efficiently applied using the low-rank one-step wave extrapolation (Sun et al., 2015), which we also refer to as the low-rank phase shift plus interpolation (PSPI) operator because of its resemblance to the well-known PSPI method for solving the one-way wave equation (Gazdag and Squazzero, 1984; Kesinger, 1992; Margrave and Ferguson, 1999). The detailed formulation of low-rank PSPI operator, as well as the derivation of its adjoint operator, the low-rank NSPS operator, is shown in the Appendix A. LSRTM in viscoacoustic media can therefore be constructed using the forward and adjoint operators.

Viscoacoustic RTM and RTDM

To obtain a seismic image with an attenuated data from the $i$ th shot $d_{i} (x_{r}, t)$ ⁠, where $x_{r}$ denotes the receiver location, viscoacoustic RTM can be carried out in the following three steps:

Forward propagate the source wavefield $S_{i} (x, t)$ by solving
$\frac{1}{c^{2}} \frac{\partial^{2} S_{i} (x, t)}{\partial t^{2}} - η L S_{i} (x, t) - τ \frac{\partial}{\partial t} H S_{i} (x, t) = δ (x - x_{i}) f (t) .$
(13)
Backward propagate the receiver wavefield $R_{i} (x, t)$ by injecting the observed seismic record as the boundary condition $R_{i} (x_{r}, t) = d_{i} (x_{r}, t)$ and solving
$\frac{1}{c^{2}} \frac{\partial^{2} R_{i} (x, t)}{\partial t^{2}} - η L R_{i} (x, t) - τ \frac{\partial}{\partial t} H R_{i} (x, t) = 0 .$
(14)
Apply the crosscorrelation imaging condition (Claerbout, 1985):
$I (x) = \sum_{i} \sum_{t} S_{i} (x, t) {\bar{R}}_{i} (x, t),$
(15)
where $\overline{R}$ denotes the complex conjugate of $R$ (Sun and Fomel, 2013).

The RTDM in viscoacoustic media can be formulated as the adjoint of the RTM process:

Calculate the source wavefield $S_{i} (x, t)$ in the background velocity model in the same manner as RTM by solving equation 13.
Generate the receiver wavefield $R_{i} (x, t)$ by using the stacked image $I (x)$ as a secondary source and solving
$\frac{1}{c^{2}} \frac{\partial^{2} R_{i} (x, t)}{\partial t^{2}} - η L R_{i} (x, t) - τ \frac{\partial}{\partial t} H R_{i} (x, t) = I (x) S_{i} (x, t) .$
(16)
Extract the predicted seismic record (denoted by the hat) at receiver locations $x_{r}$ ⁠:
${\hat{d}}_{i} (x_{r}, t) = R_{i} (x_{r}, t) .$
(17)

Note that, to make the RTDM adjoint to RTM, the wave extrapolation operator used to solve equation 16 needs to be the adjoint of the operator used to solve equation 14. For example, if low-rank PSPI is used to solve equation 14, then low-rank NSPS (derived in Appendix A) needs to be used to solve equation 16. If we write the RTM process symbolically as $\hat{m} = A^{*} d$ ⁠, where $\hat{m}$ is the stacked image, $A^{*}$ is the viscoacoustic RTM operator, and $*$ denotes the adjoint, then the RTDM process corresponds to $\hat{d} = Am$ ⁠, where $\hat{d}$ represents the predicted data and $A$ is the viscoacoustic RTDM operator.

Least-squares RTM in viscoacoustic media

LSRTM aims to minimize the misfit between the observed data and predicted data measured by the quadratic function:

J (m) = \frac{1}{2} {‖ \hat{d} - d ‖}_{2}^{2} = \frac{1}{2} {‖ Am - d ‖}_{2}^{2} .

(18)

Because

A

is a linear operator, a gradient-based local optimization method, such as the conjugate gradient (CG) method, is usually applied to iteratively update the image (Dai and Schuster, 2013; Xue et al., 2014). The

J (m)

is minimized when

m

satisfies (Tarantola, 2005)

m = {(A^{*} A)}^{- 1} A^{*} d .

(19)

The square matrix

A^{*} A

is known as the wave-equation Hessian, and its condition number affects the convergence rate of LSRTM implemented as an iterative inversion (Plessix and Mulder, 2004). In acoustic media, RTM usually provides a good approximation to the inverse of RTDM, and the Hessian matrix is well-conditioned (Symes, 2008). However, in viscoacoustic media because both RTM and RTDM operators attenuate seismic waves, the image obtained by the aforementioned algorithm suffers from twice the amplitude loss accumulated along the reflection wavepath. Therefore, differently from the pure acoustic case, viscoacoustic RTM provides a poor approximation to the inverse of viscoacoustic RTDM, which makes the Hessian matrix

A^{*} A

ill-conditioned. In the presence of strong attenuation and without proper preconditioning, this slows down the convergence rate of an iterative method like CG and, in practice, may require a prohibitively large number of iterations to achieve a satisfactory result.

$Q$ -compensated LSRTM using the GMRES method

To compensate for attenuation in seismic images, Zhu et al. (2014) and Sun et al. (2015) propose the $Q$ -compensated RTM (⁠ $Q$ -RTM). $Q$ -RTM in general can be formulated as follows: Notice the sign reversal in front of $τ$ in equations 20 and 21 in comparison with equations 13 and 14. This reversal aims to recover the attenuated wavefield by undoing phase distortion and amplifying the amplitude. Practically, it may also exponentially increase noise through each time step. A low-pass filter can be applied to stabilize the extrapolation process. Another robust compensation strategy is developed by Sun and Zhu (2015) based on stable division between wavefields. The source and receiver wavefields need to be compensated to accumulate compensation along the entire reflection wavepath. Because $Q$ -RTM is capable of restoring the attenuated energy in the seismic image (Zhu et al., 2014; Sun et al., 2015), it is reasonable to expect that $Q$ -RTM is better than viscoacoustic RTM in approximating the inverse of viscoacoustic RTDM.

Forward propagate the source wavefield $S_{i} (x, t)$ with $Q$ -compensation by solving
$\frac{1}{c^{2}} \frac{\partial^{2} S_{i} (x, t)}{\partial t^{2}} - η L S_{i} (x, t) + τ \frac{\partial}{\partial t} H S_{i} (x, t) = δ (x - x_{i}) f (t) .$
(20)
Backward propagate the receiver wavefield $R_{i} (x, t)$ with $Q$ -compensation by solving
$\frac{1}{c^{2}} \frac{\partial^{2} R_{i} (x, t)}{\partial t^{2}} - η L R_{i} (x, t) + τ \frac{\partial}{\partial t} H R_{i} (x, t) = 0,$
(21)
with the boundary condition $R_{i} (x_{r}, t) = d_{i} (x_{r}, t)$ ⁠.
Apply the imaging condition (equation 15).

We propose to replace the original viscoacoustic RTM

A^{*}

with

Q

-RTM

A_{c}^{*}

as the backward operator. The true model defined in equation can be equivalently expressed as

m = {(A_{c}^{*} A)}^{- 1} A_{c}^{*} d .

(22)

Additionally, an RTM image may contain low-frequency noise, which can be efficiently removed by a Laplacian filter (Zhang and Sun, 2009). We propose to cascade the

Q

-RTM operator with a Laplacian filter to help with the least-squares inversion and speed up the convergence rate. Correspondingly, the inverted model is expressed as

m = {({LA}_{c}^{*} A)}^{- 1} {LA}_{c}^{*} d,

(23)

where

L

denotes the Laplacian operator. Because the operator

{LA}_{c}^{*} A

is closer to an identity operator than

A^{*} A

⁠, the iterative inversion process will converge faster. Equation can be viewed as the solution of the preconditioned (weighted) least-squares system that seeks to minimize

J_{p} (m) = \frac{1}{2} {‖ P (Am - d) ‖}_{2}^{2},

(24)

which leads to the solution

m = (A^{*} P^{*} PA)^{- 1} A^{*} P^{*} P d .

(25)

Instead of looking for the preconditioner $P$ ⁠, we simply replace $A^{*} P^{*} P$ with ${LA}_{c}^{*}$ ⁠. Note that, theoretically, the inverted matrix in equation 25 is Hermitian. However, the new formulation (equation 23) makes the inverted matrix numerically non-Hermitian. One complication with equation 23 is that because the square matrix being inverted is no longer Hermitian, iterative methods for Hermitian positive-definite matrices are not optimal (Saad, 2003). Therefore, we implement a complex-valued restarted GMRES(m) with m iterations between restarts, which solves a least-squares system by searching for the vector in the Krylov subspace with minimum residual (Saad and Schultz, 1986). We refer to the method of solving equation 23 by GMRES(m) as $Q$ -compensated LSRTM or $Q$ -LSRTM. As demonstrated in the numerical examples of the next section, $Q$ -LSRTM is capable of achieving a significantly faster convergence rate than conventional LSRTM and, in practice, can produce the desired image within only a few iterations.

NUMERICAL EXAMPLES

To test the convergence rate of $Q$ -LSRTM, we use a portion of the BP 2004 velocity model (Billette and Brandsberg-Dahl, 2004) and the corresponding $Q$ model suggested by Zhu et al. (2014) (Figure 1). The model features a low-velocity, low- $Q$ area, which is assumed to be caused by the presence of a gas chimney. The model has a spatial sampling rate of $12.5 m$ along vertical and horizontal directions. A total of 31 shots with a spacing of 162.5 m have been modeled with attenuation, and the source is a Ricker wavelet with 22.5 Hz peak frequency (Figure 2). Performing RTM without compensating for amplitude loss, i.e., using the dispersion-only operator, leads to an image corresponding to Figure 3b, which suffers from poor illumination below the gas chimney. In contrast, $Q$ -RTM appears capable of recovering the amplitude at deeper reflectors (Figure 3c), but the image still exhibits some differences from the true reflectivity. Note that the dispersion-only RTM image and $Q$ -RTM image have the same phase but differ in amplitude. Next, we perform LSRTM (equation 19) and $Q$ -LSRTM (equation 23) through a number of iterations. To test the separate effect of applying a Laplacian filter without compensating for attenuation, we also perform LSRTM with a Laplacian operator that removes low-frequency artifacts. For fairness of comparison, all three methods are driven by the GMRES method and because the tested model is small enough, they were not run in a restarted fashion. Using the original LSRTM (Figure 4), the inversion process attempts to remove low-frequency noise and improve the illumination of deeper reflectors. However, at 30th iteration, the reflector amplitude and sharpness beneath the attenuating area still have not been recovered. LSRTM with a Laplacian filter (Figure 5) achieves a somewhat sharper image, but because a Laplacian filter boosts high-frequency components in the image, the reflectors beneath the attenuating zone remain poorly illuminated. In contrast, the proposed $Q$ -LSRTM method (Figure 6) produces sharper reflectors with well-balanced illumination, especially in the area beneath the gas chimney using the same number of iterations. Note that the color scales used in all the three cases are kept the same as that of the true model (Figure 3a). Figure 7 compares the image traces extracted at $X = 2500 m$ from the 30th iteration results against the true model. Clearly, the result obtained by the proposed $Q$ -LSRTM best represents the true reflectivity, especially at deeper parts beneath the gas chimney (below 800 m depth).

To measure the convergence rate, we calculate the model residual as the

L_{2}

norm of the misfit between the model calculated at each iteration

m_{k}

and the true model

m^{*}

normalized by the

L_{2}

norm of the true model:

r = \frac{{‖ m_{k} - m^{*} ‖}_{2}^{2}}{{‖ m^{*} ‖}_{2}^{2}} .

(26)

Figure 8 shows the comparison of convergence rates. With the help of a Laplacian filter, LSRTM is able to achieve a slightly faster convergence rate at early iterations than the original LSRTM. The proposed

Q

-LSRTM, on the other hand, converges significantly faster than the other two methods. Convergence is achieved by

Q

-LSRTM within approximately 50 iterations, whereas the other two methods have not converged even after 100 iterations. The fast convergence is an important property because for large-scale 3D seismic imaging problems, only a few iterations can be afforded in practice.

DISCUSSION

In this work, the GMRES method is used to invert the non-Hermitian matrices. Similar to the CG method, the full GMRES method (with no restarts) converges in no more than $n$ steps, where $n$ is the total size of the model. However, the GMRES method requires additional memory to store the previous stepping directions. In the numerical examples presented above, because the model space is small enough, no restarts are needed. However, for large 3D models, restarts might be required for a large number of iterations, which could compromise global optimality. Fortunately, the proposed method is designed to achieve a satisfying result within only a small number of iterations. In practical applications where each iteration consumes large computing resources, only a small number of iterations is usually afforded.

The goal of preconditioning LSRTM in viscoacoustic media using $Q$ -RTM is to alleviate the computational burden on iterative inversion by compensating for attenuation explicitly in wave propagation. Therefore, the iterations can be spent on removing migration artifacts and compensating irregularities in subsurface illumination caused by other reasons, such as acquisition footprint. Due to attenuation, the events of the reflectors beneath the attenuating zone yield a smaller amplitude compared with unattenuated events and approximately correspond to smaller eigenvalues of the forward operator (Blanch and Symes, 1994). Inversion routines based on Krylov subspace methods, such as CG and GMRES, will favor large eigenvalues, which approximately correspond to shallower and unattenuated reflectors. This leads to the observed behavior of LSRTM without $Q$ -compensation, which first focuses on improving shallow reflectors above the attenuation zone. Blanch and Symes (1994) suggest a simple way of assigning more weights to deeper reflectors, by postconditioning the seismic record with an increasing function of time. The proposed method is similar in spirit but more accurate, in that $Q$ -compensation removes the true effect of attenuation in the gradient by accurately compensating for attenuation along the entire wavepath.

CONCLUSIONS

We have introduced a novel way of preconditioning least-squares RTM to achieve a faster convergence rate in viscoacoustic media. The data-space preconditioner is defined by the $Q$ -compensated RTM operator, with the goal of recovering amplitude loss due to attenuation and removing low-frequency artifacts in the gradient. Because the square matrix to be inverted becomes numerically non-Hermitian, we adopt the GMRES algorithm to perform iterative inversion. Our synthetic examples show that the proposed $Q$ -LSRTM is capable of producing an accurate $Q$ -compensated image within significantly fewer iterations than LSRTM, and thus it is preferable in application to accurate seismic imaging in attenuating media.

ACKNOWLEDGMENTS

We thank J. Blanch, A. Stanton, and one anonymous reviewer for their suggestions that helped to improve the paper. We thank M. Sen and L. Ying for constructive discussions. We thank the sponsors of the Texas Consortium for Computation Seismology (TCCS) for financial support. The first author is supported additionally by the Statoil Fellows Program at the University of Texas at Austin, and the third author is supported by the Jackson School Distinguished Postdoctoral Fellowship at the University of Texas at Austin. We thank Texas Advanced Computing Center (TACC) for providing computational resources used in this study.

DERIVATION OF THE LOW-RANK PSPI AND NSPS OPERATORS

Let

p (x, t)

be the seismic wavefield at location

x

and time

t

⁠, with the spatial Fourier transform denoted by

P (k, t)

⁠. The wavefield at the next time step

t + Δ t

can be approximated by the Fourier integral operator (Wards et al., 2008; Fomel et al., 2013):

p (x, t + Δ t) = \int P (k, t) e^{i ϕ (x, k, Δ t) + i x \cdot k} d k,

(A-1)

where

ϕ (x, k, Δ t)

is the phase function. The mixed-domain operator in equation is also referred to as the one-step wave extrapolation operator (Zhang and Zhang, 2009; Sun and Fomel, 2013). Because the wave extrapolation matrix is complex, it propagates a complex wavefield with the imaginary part being the Hilbert transform of the real part,

P (k, t) = P_{r} \pm F [H (p_{r} (x, t))],

(A-2)

where

P

and

P_{r}

⁠, respectively, denote the complex wavefield and real wavefield, and

F

denotes spatial Fourier transform (Zhang and Zhang, 2009).

Converting the dual-domain expression into the space domain, we obtain

p (x, t + Δ t) = \int p (y, t) \int e^{i ϕ (x, k, Δ t) + i k \cdot (x - y)} d k d y .

(A-3)

The adjoint form of operator can be written as

P (k, t) = \int p (x, t + Δ t) e^{- i \overline{ϕ} (x, k, Δ t) - i x \cdot k} d x .

(A-4)

where

\overline{ϕ}

denotes the complex conjugate of

ϕ

⁠. The

- i \overline{ϕ}

term in equation indicates stepping backward in time. Expressing the dual-domain operator in the space domain and stepping forward in time, we arrive at a different operator:

p (x, t + Δ t) = \int p (y, t) \int e^{i ϕ (y, k, Δ t) + i k \cdot (x - y)} d k d y .

(A-5)

The phase function in equation depends on the output space

x

and thus represents a kind of nonstationary combination filter (Margrave, 1998). In comparison, the phase function appearing in equation depends on the input space

y

and leads to a kind of nonstationary convolution filter. Operators and apply the same wave-propagation phase function

ϕ (x, k, Δ t)

⁠; the one-step low-rank wave extrapolation operator applies the phase function in the wavenumber domain after forward Fourier transform, whereas the new operator applies the phase function in the space domain before the inverse Fourier transform. The essential difference between the two is that the nonstationary convolution has the physical interpretation of scaled, linear superposition of the nonstationary filter impulse responses, as suggested by Huygens’ principle, whereas nonstationary combination filters do not have such implications (Margrave, 1998). The mixed-domain operator is an analog of the limiting case of PSPI method (Gazdag and Squazzero, 1984; Kesinger, 1992), which has been a popular choice for one-way wavefield extrapolators. The proposed operator is analogous to the nonstationary phase shift (NSPS) method (Margrave, 1998; Margrave and Ferguson, 1999) for one-way wave extrapolation.

The low-rank algorithm introduced by Fomel et al. (2013) is a separable approximation that selects a set of

N

representative spatial locations and

M

representative wavenumbers, which correspond to rows and columns from the original wave-propagation matrix. The low-rank one-step wave extrapolation uses low-rank decomposition to approximate the mixed-domain phase function in equation :

p (x, t + Δ t) \approx \sum_{m = 1}^{M} W (x, k_{m}) (\sum_{n = 1}^{N} a_{m n} (\int e^{i x \cdot k} W (x_{n}, k) P (k, t) d k)),

(A-6)

whose computational cost effectively equals that of applying

N

inverse fast Fourier transforms per time step, where

N

is the approximation rank and is typically a number less than ten.

With the help of low-rank decomposition, the computational effort for the new NSPS method can be made identical to that of the low-rank PSPI wave extrapolation, by approximating the wave propagation operator appearing in equation with

P (k, t + Δ t) \approx \sum_{n = 1}^{N} W (x_{n}, k) (\sum_{m = 1}^{M} a_{m n} (\int e^{- i x \cdot k} W (x, k_{m}) p (x, t) d x)) .

(A-7)

Note that, for simplicity of notation, equations and do not include an operation of the forward and inverse Fourier transforms in place between

p (x, t)

and

P (k, t)

⁠.

In this appendix, we have derived the adjoint operator to low-rank one-step wave extrapolation. Because the derivation of the adjoint operator is discrete, i.e., using the low-rank approximation matrices instead of state and adjoint state equations, the result presented in the appendix has wider applications not limited to the scope of this paper. For example, in full-waveform inversion applications, where the adjoint state equation is difficult to obtain, the discrete adjoint operator described in this appendix can be efficiently applied to calculate the adjoint state variable. This strategy is usually referred to as discretize then optimize (Betts, 2010).

© The Authors.Published by the Society of Exploration Geophysicists. All article content, except where otherwise noted (including republished material), is licensed under a Creative Commons Attribution 4.0 Unported License (CC BY-SA). See http://creativecommons.org/licenses/by/4.0/. Distribution or reproduction of this work in whole or in part commercially or noncommercially requires full attribution of the original publication, including its digital object identifier (DOI). Derivatives of this work must carry the same license.

The Authors

Freely available online through the SEG open-access option.