% Started 09/15/97%\shortnote\lefthead{Fomel \& Claerbout}\righthead{Implicit extrapolation}\title{Exploring three-dimensional implicit wavefield extrapolation  with the helix transform}\email{sergey@sep.stanford.edu, jon@sep.stanford.edu}%\keywords{depth migration, post-stack, finite-difference, velocity continuation, helix}\author{Sergey Fomel and Jon F. Claerbout}\maketitle\begin{abstract}  Implicit extrapolation is an efficient and unconditionally stable  method of wavefield continuation. Unfortunately, implicit wave  extrapolation in three dimensions requires an expensive solution of  a large system of linear equations. However, by mapping the  computational domain into one dimension via the helix transform, we  show that the matrix inversion problem can be recast in terms of an  efficient recursive filtering. Apart from the boundary conditions,  the solution is exact in the case of constant coefficients (that is,  a laterally homogeneous velocity.) We illustrate this fact with an  example of three-dimensional velocity continuation and discuss  possible ways of attacking the problem of lateral variations.\end{abstract}\section{Introduction}Implicit finite-difference wavefield extrapolation played anexceptionally important role in the early development of seismicmigration methods. Using limited-degree approximations to the one-waywave equation, implicit schemes have provided efficient andunconditionally stable numerical wave extrapolation operators\cite[]{Godfrey.sep.16.83,Claerbout.blackwell.85}. Unfortunately, theadvantages of \emph{implicit} methods were lost with the developmentof three-dimensional seismic exploration. While the cost of 2-Dimplicit extrapolation is linearly proportional to the mesh size, thesame approach, applied in the 3-D case, leads to a nonlinearcomputational complexity. Primarily for this reason, implicitextrapolators were replaced in practice by \emph{explicit} ones,capable of maintaining linear complexity in all dimensions. A numberof computational tricks \cite[]{GEO56-11-17701777} allow the commonlyused explicit schemes to behave stably in practical cases.  However,their stability is not unconditional and may break in unusualsituations \cite[]{SEG-1994-1266}.\parIn this paper, we present an approach to three-dimensionalextrapolation, based on the helix transform of multidimensionalfilters to one dimension \cite[]{Claerbout.gem.97}. The traditionalapproach involves an inversion of a banded matrix (tridiagonal in the2-D case and blocked-tridiagonal in the 3-D case). With the help ofthe helix transform, we can recast this problem in terms of inverserecursive filtering.  The coefficients of two-dimensional filters on ahelix are obtained by one-dimensional spectral factorization methods.As a result, the complexity of three-dimensional implicitextrapolation is reduced to a linear function of the computationalmesh size. This approach doesn't provide an exact solution in thepresence of lateral velocity variations. Nevertheless, it can be usedfor preconditioning iterative methods, such as those described by\cite{Nichols.sep.70.31}.  In this paper, we demonstrate thefeasibility of 3-D implicit extrapolation on the example of laterallyinvariant velocity continuation and, in the final part, discusspossible strategies for solving the problem of lateral variations.\parThe main application of finite-difference wave extrapolation is\emph{post-stack} depth migration. An application of similar methodsfor \emph{prestack} common-shot migration is constrained by thelimited aperture of commonly used seismic acquisition patterns.Recently developed acquisition methods, such as the vertical cabletechnique \cite[]{SEG-1993-1376}, open up new possibilities for 3-D waveextrapolation applications. An alternative approach is common-azimuthmigration \cite[]{Biondi.sep.80.109,Biondi.sep.93.1}. Other interestingapplications include finite-difference data extrapolation in offset\cite[]{Fomel.sep.84.179}, migration velocity \cite[]{Fomel.sep.92.159},and anisotropy \cite[]{Alkhalifah.sep.94.tariq3}.\section{Implicit versus Explicit extrapolation}The difference between implicit and explicit extrapolation is bestunderstood through an example. Following \cite{Claerbout.blackwell.85},let us consider, for instance, the diffusion (heat conduction) equationof the form\begin{equation}  {\frac{\partial T}{\partial t}} = {a (x)\,{\frac{\partial^2 T}{\partial x^2}}}\;.\label{eqn:heat}\end{equation}Here $t$ denotes time, $x$ is the space coordinate, $T (x,t)$ is thetemperature, and $a$ is the heat conductivity coefficient.Equation (\ref{eqn:heat}) forms a well-posed boundary-value problem ifsupplied with the initial condition\begin{equation}  \label{eqn:heatinit}  \left.T\right|_{t=0} = T_0 (x)\end{equation}and the appropriate boundary conditions. Our task is to build adigital filter, which transforms a gridded temperature $T$ from onetime level to another.\parIt helps to note that when the conductivity coefficient $a$ isconstant and the space domain of the problem is infinite (or periodic)in $x$, the problem can be solved in the wavenumber domain. Indeed,after the Fourier transform over the variable $x$, equation(\ref{eqn:heat}) transforms to the ordinary differential equation\begin{equation}  {\frac{d \hat{T}}{d t}} = {- a k^2\, \hat{T}}\;,  \label{eqn:heatk}	\end{equation}which has the explicit analytical solution\begin{equation}  \label{eqn:heatsol}  \hat{T} (k,t) = \hat{T}_0 (k) e^{- a k^2 t}\;,\end{equation}where $\hat{T}$ denotes the Fourier transform of $T$, and $k$ standsfor the wavenumber. Therefore, the desired filter in thewavenumber domain has the form\begin{equation}  \label{eqn:heatf}  H (k) = e^{- a k^2}\;,\end{equation}where for simplicity the coefficient $a$ is normalized for the timestep $\triangle t$ equal to $1$.\parReturning now to the time-and-space domain, we can approach the filterconstruction problem by approximating the space-domain response offilter (\ref{eqn:heatf}) in terms of the differential operators$\frac{\partial^2}{\partial x^2} = - k^2$, which can be approximatedby finite differences. An \emph{explicit} approach would amount toconstructing a series expansion of the form\begin{equation}  \label{eqn:heatexpl}  H_{\mbox{ex}} (k) \approx a_0 + a_1 k^2 + a_2 k^4 + \ldots\;,\end{equation}and selecting the coefficients $a_j$ to approximate equation(\ref{eqn:heatf}). For example, the three-term Taylor series expansionaround the zero wavenumber yields\begin{equation}  \label{eqn:heattayl}  H_{\mbox{ex}} (k) = 1 - a\,{k^2}  +  {\frac{{{a }^2}\,{k^4}}{2}} \;.\end{equation}The error of approximation (\ref{eqn:heattayl}) as a function of $k$for two different values of $a$ is shown in the left plot of Figure\ref{fig:error}.\inputdir{Math}\plot{error}{width=6in}{Errors of second-order explicit and implicit  approximations for the heat extrapolation.}\parAn \emph{implicit} approach also approximates the ideal filter(\ref{eqn:heatf}), but with a rational approximation of the form\begin{equation}  \label{eqn:heatpade}  H_{\mbox{im}} (k) \approx \frac{b_0 + b_1 k^2 + b_2 k^4 + \ldots}  {1 + c_1 k^2 + c_2 k^4 + \ldots}\;.\end{equation}One way of selecting the coefficients $b_i$ and $c_i$ is to apply anappropriate Pad\'{e} approximation \cite[]{pade}\footnote{If the  denominator and the numerator have the same order, Pad\'{e}  approximants are equivalent to the corresponding continuous  fraction expansions.}.  For example the $[2/2]$ Pad\'{e}approximation is\begin{equation}  \label{eqn:heatcrank}  H_{\mbox{im}} (k) =  \frac{1 - \frac{a}{2}\,k^2}{1 + \frac{a}{2}\,k^2}  \;.\end{equation}This approximation corresponds to the famous Crank-Nicolson implicitmethod \cite[]{cn}. The error of approximation (\ref{eqn:heatcrank}) asa function of $k$ for different values of $a$ is shown in the rightplot of Figure \ref{fig:error}. Not only is it significantly smallerthan the error of the same-order explicit approximation, but it alsohas a negative sign. It means that the high-frequency numerical noisegets suppressed rather than amplified. In practice, this propertytranslates into a stable numerical extrapolation.\parThe second derivative operator $-k^2$ can be approximated in practiceby a digital filter. The most commonly used filter has the$Z$-transform $D_2 (Z) = -Z^{-1} + 2 - Z$, and the Fourier transform\begin{equation}  \label{eqn:d2k}  D_2 (k) = e^{-ik} - 2 + e^{-ik} = 2 (\cos{k} - 1) = -4  \sin^2{\frac{k}{2}}\;.\end{equation}Formula (\ref{eqn:d2k}) approximates $-k^2$ well only for smallwavenumbers $k$. As shown in Appendix A, the implicit scheme allowsthe accuracy of the second-derivative filter to be significantlyimproved by a variation of the ``1/6-th trick''\cite[]{Claerbout.blackwell.85}. The final form of the implicitextrapolation filter is\begin{equation}  \label{eqn:heatfk}   H_{\mbox{im}} (k) =   \frac{1 + \frac{a+\beta}{2}\,D_2 (k)}{1 - \frac{a-\beta}{2}\,D_2 (k)}   \;,\end{equation}where $\beta$ is a numerical constant, found in Appendix A.\inputdir{heat}\plot{heat}{width=6in,height=2.5in}{Heat extrapolation with explicit  and implicit finite-different schemes. Explicit extrapolation  appears stable for $a=2/3$ (left plot) and unstable for $a=4/3$  (middle plot). Implicit interpolation is stable even for larger  values of $a$ (right plot).}\parA numerical 1-D example is shown in Figure \ref{fig:heat}. The initialtemperature distribution is given by a step function. Thediscontinuity at the step gets smoothed with time by the heatdiffusion. The left plot shows the result of an explicit extrapolationwith $a=2/3$, which appears stable. The middle plot is an explicitextrapolation with $a=4/3$, which shows a terribly unstable behavior:the high-frequency numerical noise is amplified and dominates thesolution. The right plot shows a stable (though not perfectlyaccurate) extrapolation with the implicit scheme for the larger value of$a=2$.\parThe difference in stability between explicit and implicit schemes iseven more pronounced in the case of \emph{wave extrapolation}. Forexample, let us consider the ideal depth extrapolation filter in theform of the phase-shift operator\cite[]{GEO43-07-13421351,Claerbout.blackwell.85}\begin{equation}  \label{eqn:gazdag}  W (k) = e^{i \sqrt{a^2 - k^2}}\;,\end{equation}where $a = \omega / v$, $\omega$ is the time frequency, and $v$ is theseismic velocity (which may vary spatially); we assume for simplicitythat both the depth step $\triangle z$ and the space sampling$\triangle x$ are normalized to $1$.  A simple implicit approximationto filter (\ref{eqn:gazdag}) is\begin{equation}  \label{eqn:wave45}   W_{\mbox{im}} (k) = e^{i a}\,   \frac{1 -4 a^2 + i a\,k^2}{1 - 4 a^2 - i a\,k^2} = e^{i \phi}\;,\end{equation}where $\phi = a - 2 \arctan{\frac{a\,k^2}{4 a^2-1}}$. We can seethat approximation (\ref{eqn:wave45}) is again a pure phase shiftoperator, only with a slightly different phase. For that reason, theoperator is unconditionally stable for all values of $a$: the totalwave energy from one depth level to another is preserved. Operator(\ref{eqn:gazdag}) corresponds to the Crank-Nicolson scheme for the45-degree one-way wave equation \cite[]{Claerbout.blackwell.85}. Itsphase error as a function of the dip angle $\theta =\arcsin{\frac{k}{a}}$ for different values of $a$ is shown in Figure\ref{fig:phase}.\inputdir{Math}