\sideplot{phase}{width=2.5in}{The phase error of the  implicit depth extrapolation with the Crank-Nicolson method.}\parThe unconditional stability property is not achievable with theexplicit approach, though it is possible to increase the stability ofexplicit operators by using relatively long filters\cite[]{GPR36-02-00990114,GEO56-11-17701777}.\section{Spectral factorization and three-dimensional extrapolation}In this section, we continue our review of extrapolation methods toreveal the principal difficulties of three-dimensional extrapolation.We then describe a new, helix-transform approach to this old andfascinating problem.\subsection{Inverse filter factorization}The conventional way of applying implicit finite-difference schemesreduces to solving a system of linear equations with a sparse matrix.For example, to apply the scheme of equation (\ref{eqn:heatfk}), wecan put the filter denominator on the other side of the extrapolationequation, writing it as\begin{equation}  \label{eqn:heattris}  \left(\mathbf{I} - \frac{a-\beta}{2}\,\mathbf{D}_2\right) \mathbf{T}_{t+1} =  \left(\mathbf{I} + \frac{a+\beta}{2}\,\mathbf{D}_2\right) \mathbf{T}_t\;,\end{equation}where $\mathbf{I}$ is the identity matrix, $\mathbf{D}$ is the convolutionmatrix for filter (\ref{eqn:d2k}), and $\mathbf{T}_t$ is the vector oftemperature distribution at time level $t$. In the case oftwo-dimensional extrapolation, the matrix on the left side of equation(\ref{eqn:heattris}) takes the tridiagonal form\begin{equation}\label{eqn:matrix}  \mathbf{A} = \left(\mathbf{I} -c\,\mathbf{D}_2\right) =  \left[\begin{array}{cccccc}  1+2 c_{1}   & -c_{1}     &  0     & \cdots &        & 0      \\  -c_{2}     & 1+2c_{2}   & -c_{2}     & 0      & \cdots &        \\  0      & -c_{3}     & \cdots & \cdots &        & \cdots \\  \cdots & 0      & \cdots &        &        &        \\         & \cdots &        &        & \cdots & -c_{n-1}     \\  0      &        &        &        & -c_{n}     & 1 + 2c_{n}  \end{array}\right]\;,\end{equation}where $c = \frac{a-\beta}{2}$, and where, for simplicity, we assumezero-slope boundary conditions. Like any positive-definite tridiagonalmatrix, matrix $\mathbf{A}$ can be inverted recursively by an $LU$decomposition into two bidiagonal matrices. The cost of inversion isdirectly proportional to the number of vector components.  The sameconclusion holds for the case of depth extrapolation [equation(\ref{eqn:wave45})] with the substitution $c = \frac{\beta + i  a}{1-4\,a^2}$.\parIn the case of a laterally constant coefficient $a$, we can take adifferent point of view on the tridiagonal matrix inversion. In thiscase, the matrix $\mathbf{A}_2$ represents a convolution with asymmetric three-point filter $1-c\,D_2(k)$. The $LU$decomposition of such a matrix is precisely equivalent to filter\emph{factorization} into the product of a causal minimum-phase filterwith its adjoint. This conclusion can be confirmed by the easilyverified equality\begin{equation}\label{eqn:d2d1}  1 + c (Z^{-1} - 2 + Z) = \frac{(1+b)^2}{4}\, \left(1 + \frac{1-b}{1+b} Z\right)  \,\left(1 + \frac{1-b}{1+b} Z^{-1}\right)\;,\end{equation}where $b = \sqrt{1+ 4\,c}$. The inverse of the causal minimum-phasefilter $1 + \frac{1-b}{1+b} Z$ is a recursive inverse filter.Correspondingly, the inverse of its adjoint pair, $1 + \frac{1-b}{1+b}Z^{-1}$, is the same inverse filtering, performed in the adjoint mode(backwards in space).  In the next subsection, we show how thisapproach can be carried into three dimensions by applying the helixtransform.\subsection{Helix and multidimensional deconvolution}The major obstacle of applying an implicit extrapolation in threedimensions is that the inverted matrix is no longer tridiagonal. If weapproximate the second derivative (Laplacian) on the 2-D plane withthe commonly used five-point filter $Z_x^{-1} + Z_y^{-1} -4 + Z_x +Z_y$, then the matrix on the left side of equation(\ref{eqn:heattris}), under the usual mapping of vectors from atwo-dimensional mesh to one dimension, takes the infamousblocked-tridiagonal form \cite[]{birk}\begin{equation}\label{eqn:matrix2}  \tilde{\mathbf{A}} = \left(\mathbf{I} -c\,\mathbf{D}_2\right) =  \left[\begin{array}{cccccc}  \mathbf{A}_1   & -c_1 \,\mathbf{I}    &  0     & \cdots &        & 0      \\  -c_2 \,\mathbf{I}    & \mathbf{A}_2   & -c_2 \,\mathbf{I}    & 0      & \cdots &        \\  0      & -c_3 \,\mathbf{I}     & \cdots & \cdots &        & \cdots \\  \cdots & 0      & \cdots &        &        &        \\         & \cdots &        &        & \cdots & -c_{n-1} \,\mathbf{I}    \\  0      &        &        &        & -c_{n} \,\mathbf{I}    &  \mathbf{A}_n  \end{array}\right]\;.\end{equation}Inspecting this form more closely, we see that the main diagonal of$\tilde{\mathbf{A}}$, as well as the two offset diagonals formed by thescaled identity matrices, remains continuous, while the second top andbottom diagonals are broken. Therefore, even for constant $c$, theinverted matrix does not have a simple convolutional structure, andthe cost of its inversion grows nonlinearly with the number of gridpoints.\parA \emph{helix transform}, recently proposed by one of the authors\cite[]{Claerbout.sep.95.jon1}, sheds new light on this old problem.Let us assume that the extrapolation filter is applied by sliding italong the $x$ direction in the $\{x,y\}$ plane. The diagonaldiscontinuities in matrix $\tilde{\mathbf{A}}$ occur exactly in theplaces where the forward leg of the filter slides outside thecomputational domain. Let's imagine a situation, where the leg of thefilter that went to the end of the $x$ column, would immediatelyappear at the beginning of the next column. This situation defines adifferent mapping from two computational dimensions to the onedimension of linear algebra. The mapping can be understood as thehelix transform, illustrated in Figure \ref{fig:helix1} and explainedin detail by \cite{Claerbout.sep.95.jon1}. According to thistransform, we replace the original two-dimensional filter with a longone-dimensional filter. The new filter is partially filled with zerovalues (corresponding to the back side of the helix), which can besafely ignored in the convolutional computation.\plot{helix}{width=5in,bb=210 155 630 390}{The helix transform of  two-dimensional filters to one dimension. The two-dimensional filter  in the left plot is equivalent to the one-dimensional filter in the  right plot, assuming that a shifted periodic condition is imposed on  one of the axes.}\parThis is exactly the helix transform that is required to make all thediagonals of matrix $\tilde{\mathbf{A}}$ continuous. In the case oflaterally invariant coefficients, the matrix becomes strictly Toeplitz(having constant values along the diagonals) and represents aone-dimensional convolution on the helix surface. Moreover, thissimplified matrix structure applies equally well to largersecond-derivative filters ( such as those described in Appendix B),with the obvious increase of the number of Toeplitz diagonals.Inverting matrix $\tilde{\mathbf{A}}$ becomes once again a simpleinverse filtering problem.  To decompose the 2-D filter into a pairconsisting of a causal minimum-phase filter and its adjoint, we canapply spectral factorization methods from the 1-D filtering theory\cite[]{Claerbout.blackwell.76,Claerbout.blackwell.92}, for example,Kolmogorov's highly efficient method \cite[]{kolmog}. Thus, in the caseof a laterally invariant implicit extrapolation, matrix inversionreduces to a simple and efficient recursive filtering, which we needto run twice: first in the forward mode, and second in the adjointmode.\inputdir{heat}\plot{heat3d}{width=6in,height=2in}{Heat extrapolation in two  dimensions, computed by an implicit scheme with helix recursive  filtering. The left plot shows the input temperature distributions;  the two other plots, the extrapolation result at different time  steps.  The coefficient $a$ is 2.}\parFigure \ref{fig:heat3d} shows the result of applying the helixtransform to an implicit heat extrapolation of a two-dimensionaltemperature distribution. The unconditional stability properties arenicely preserved, which can be verified by examining the plot ofchanges in the average temperature (Figure \ref{fig:heat-mean}).\sideplot{heat-mean}{width=2.5in}{Demonstration of the stability of  implicit extrapolation. The solid curve shows the normalized mean  temperature, which remains nearly constant throughout the  extrapolation time. The dashed curve shows the normalized maximum  value, which exhibits the expected Gaussian shape.}\parIn principle, we could also treat the case of a laterally invariantcoefficient with the help of the Fourier transform. Under whatcircumstances does the helix approach have an advantage over Fouriermethods? One possible situation corresponds to a very large input datasize with a relatively small extrapolation filter. In this case, the$O(N log N)$ cost of the fast Fourier transform is comparable with the$O(N_f N)$ cost of the space-domain deconvolution (where $N$corresponds to the data size, and $N_f$ is the filter size). Anothersituation is when the boundary conditions of the problem have anessential lateral variation. The latter case may occur in applicationsof velocity continuation, which we discuss in the next section.  Laterin this paper, we return to the discussion of problems associated withlateral variations.\section{Three-dimensional implicit velocity continuation}\inputdir{vc3}Velocity continuation is a process of navigating in the migrationvelocity space, applicable for time migration, residual migration, andmigration velocity analysis \cite[]{Fomel.sep.92.159}. In thezero-offset (post-stack) case, the velocity continuation process isdescribed by the simple partial differential equation\cite[]{Claerbout.sep.48.79,me}\begin{equation}  {\frac{\partial^2 P}{\partial v\,\partial t}} +  {v\,t\,\left({\frac{\partial^2 P}{\partial x^2}} + {\frac{\partial^2      P}{\partial y^2}}\right)} = 0\;,\label{eqn:velcon}\end{equation}where $t$ is the vertical time coordinate of the migrated image, $x$and $y$ are spatial (midpoint) coordinates, and $v$ is the migrationvelocity. Slightly different versions of two-dimensional implicitextrapolation with equation (\ref{eqn:velcon}) have been described by\cite{Li.sep.48.85} and \cite[]{Fomel.sep.92.159}. \plot{velcon}{width=6in}{Impulse responses of the velocity  continuation operator, computed by an implicit, unconditionally  stable extrapolation via the helix transform. The left plot  corresponds to continuation towards higher velocities (migration  mode); the right plot, smaller velocities (modeling mode).}\parThe helix approach has allowed us to modify the old code for threedimensions. Figure \ref{fig:velcon} shows impulse responses of animplicit helix-based three-dimensional velocity continuation.\sideplot{qdome}{width=3in}{Qdome synthetic model, used for testing  the 3-D velocity continuation program.}\parFigure \ref{fig:modmig} illustrates the velocity continuation processon the Qdome synthetic model \cite[]{Claerbout.gem.97}, shown in Figure\ref{fig:qdome}. Continuation backward in velocity corresponds to the``modeling'' mode, while forward continuation corresponds to the``migration'' mode. It is possible to balance the amplitudes of thetwo processes so that the finite-difference velocity continuationbehaves as a unitary operator\cite[]{Fomel.sep.92.159,Fomel.sep.92.267}.\plot{modmig}{width=6in}{Modeling (left) and migration (right) with  the Qdome synthetic model, obtained by running the 3-D velocity  continuation backward and forward in velocity.}\section{Depth extrapolation and the v(x) challenge}Can the constant-velocity result help us achieve the challenging goalof a stable implicit depth extrapolation through media with lateralvelocity variations?\parThe first idea that comes to mind is to replace the space-invarianthelix filters with a precomputed set of spatially varying filters,which reflect local changes in the velocity fields. This approachwould merely reproduce the conventional practice of explicit depthextrapolators, popularized by \cite{GPR36-02-00990114} and\cite{GEO56-11-17701777}. However, it hides the danger of losingthe property of unconditional stability, which is obviously the majorasset of implicit extrapolators.\parAnother route, partially explored by \cite{Nichols.sep.70.31}, isto implement the matrix inversion in the three-dimensional implicitscheme by an iterative method. In this case, the helix inversion mayserve as a powerful preconditioner, providing an immediate answer inconstant velocity layers and speeding up the convergence in the caseof velocity variations. To see why this might be true, one can writethe variable-coefficient matrix $\tilde{\mathbf{A}}$ in the form\begin{equation}  \label{eqn:golub}  \tilde{\mathbf{A}} = \mathbf{B} + \mathbf{D}\;,\end{equation}where matrix $\mathbf{B}$ corresponds to some constant average velocity,and $\mathbf{D}$ is the matrix of velocity perturbations. The system oflinear equations that we need to solve is then