Appendix A: Derivation of the lowrank PSPI and NSPS operators

Let $p(\mathbf{x},t)$ be the seismic wavefield at location $\mathbf{x}$ and time $t$, with the spatial Fourier transform denoted by $P(\mathbf{k},t)$. The wavefield at the next time step $t+ \Delta t$ can be approximated by the Fourier integral operator (Fomel et al., 2013; Wards et al., 2008):

\begin{displaymath}
p(\mathbf{x},t+\Delta t) = \int P(\mathbf{k},t) e^{i \p...
...f{k},\Delta t) + i\mathbf{x} \cdot \mathbf{k}} d\mathbf{k}\;,
\end{displaymath} (27)

where $\phi (\mathbf{x},\mathbf{k},\Delta t)$ is the phase function. The mixed-domain operator in equation 27 is also referred to as the one-step wave extrapolation operator (Sun and Fomel, 2013; Zhang and Zhang, 2009). Because the wave extrapolation matrix is complex, it propagates a complex wavefield with the imaginary part being the Hilbert transform of the real part:
\begin{displaymath}
P(\mathbf{k},t)=P_r \pm F[H(p_r(\mathbf{x},t))]\;,
\end{displaymath} (28)

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 27 into the space domain, we obtain

\begin{displaymath}
p(\mathbf{x},t+\Delta t) = \int p(\mathbf{y},t) \int e^{i...
...i\mathbf{k} \cdot \mathbf{(x-y)}} d\mathbf{k} d\mathbf{y}\;.
\end{displaymath} (29)

The adjoint form of operator 27 can be written as:
\begin{displaymath}
P(\mathbf{k},t) = \int p(\mathbf{x},t+\Delta t) e^{-i \b...
...{k},\Delta t) - i\mathbf{x} \cdot \mathbf{k}} d\mathbf{x}\;.
\end{displaymath} (30)

where $\bar{\phi}$ denotes the complex conjugate of $\phi$. The $-i\bar{\phi}$ term in equation 30 indicates stepping backward in time. Expressing the dual domain operator 30 in the space domain and stepping forward in time, we arrive at a different operator:
\begin{displaymath}
p(\mathbf{x},t+\Delta t) = \int p(\mathbf{y},t) \int e^{i...
...i\mathbf{k} \cdot \mathbf{(x-y)}} d\mathbf{k} d\mathbf{y}\;.
\end{displaymath} (31)

The phase function in equation 29 depends on the output space $\mathbf{x}$, and thus represents a kind of non-stationary combination filter (Margrave, 1998). In comparison, the phase function appearing in equation 31 depends on the input space $\mathbf{y}$, and leads to a kind of non-stationary convolution filter. Both operators 29 and 31 apply the same wave-propagation phase function, $\phi (\mathbf{x},\mathbf{k},\Delta t)$; the one-step low-rank wave extrapolation operator 29 applies the phase function in the wavenumber domain after forward Fourier transform, whereas the new operator 31 applies the phase function in the space domain before the inverse Fourier transform. The essential difference between the two is that the non-stationary convolution has the physical interpretation of scaled, linear superposition of the non-stationary filter impulse responses, as suggested by Huygens' principle, whereas non-stationary combination filters do not have such implications (Margrave, 1998). The mixed-domain operator 29 is an equivalent to the most accurate limiting case of phase shift plus interpolation (PSPI) method (Kesinger, 1992; Gazdag and Squazzero, 1984), which has been a popular choice for one-way wavefield extrapolators. The proposed operator 31 is analogous to the non-stationary phase shift (NSPS) method (Margrave and Ferguson, 1999; Margrave, 1998) 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 29:

    $\displaystyle p(\mathbf{x},t+\Delta t) \approx$ (32)
    $\displaystyle \sum\limits_{m=1}^M W(\mathbf{x},\mathbf{k}_m) \left( \sum\limits...
...{k}} W(\mathbf{x}_n,\mathbf{k}) P(\mathbf{k},t) d\mathbf{k} \right) \right)\; ,$  

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 31 with

    $\displaystyle P(\mathbf{k},t+\Delta t) \approx$ (33)
    $\displaystyle \sum\limits_{n=1}^N W(\mathbf{x}_n,\mathbf{k}) \left( \sum\limits...
...{k}} W(\mathbf{x},\mathbf{k}_m) p(\mathbf{x},t) d\mathbf{x} \right) \right)\; .$  

Note that, for simplicity of notation, equations 32 and 33 do not include an operation of the forward and inverse Fourier transforms in place between $p(\mathbf{x},t)$ and $P(\mathbf{k},t)$.

In this appendix we have presented the adjoint operator to lowrank one-step wave extrapolation. Because the derivation of the adjoint operator is discrete, i.e., using the lowrank 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 calculated the adjoint state variable. This strategy is usually refered to as discretize then optimize (Betts, 2010).


2019-05-03