h\,t\,\left({4 \over {v^2}} + \left({{\partial t} \over {\partial h}}\right)^2\,-\left({{\partial t} \over {\partial x}}\right)^2\right)\;.\label{eq:OCeikonal}\end{equation}Differentiating (\ref{eq:length}) with respect to the velocity $v$ yields\begin{equation}- v^2\,{{\partial t} \over {\partial v}} = 2\,h\,{\cos{\alpha} \over \sin{\theta}}\;.\label{eq:dtdv} \end{equation}Finally, dividing (\ref{eq:dtdv}) by (\ref{eq:snells3}) produces\begin{equation}v\,{{\partial z} \over {\partial v}} = {h \over {\cos{\theta}\,\sin{\theta}}}\;.\label{eq:dzdv} \end{equation}Equation (\ref{eq:dzdv}) can be written in a variety of ways with the helpof an explicit geometric relationship between the half-offset $h$ andthe depth $z$, \begin{equation}h = z\,{{\sin{\theta}\,\cos{\theta}} \over{\cos^2{\alpha}-\sin^2{\theta}}}\;, \label{eq:z2h} \end{equation}which follows directly from the trigonometry of the triangle in Figure\ref{fig:vlcray} \cite[]{ofcon}. For example, equation (\ref{eq:dzdv}) canbe transformed to the form obtained by \cite{GEO60-01-01420153}:\begin{equation}v\,{{\partial z} \over {\partial v}} = {z \over{\cos^2{\alpha}-\sin^2{\theta}}} ={z \over{\cos{\alpha_1}\,\cos{\alpha_2}}}\;.\label{eq:liu} \end{equation}In order to separate different factors contributing to the velocitycontinuation process, one can transform this equation to the form\begin{eqnarray}\nonumberv\,{{\partial z} \over {\partial v}} & = & {z \over {\cos^2{\alpha}}} +{{h^2} \over z}\,\left(1-\tan^2{\alpha}\,\tan^2{\theta}\right) \\& = & z\,\left(1 + \left({{\partial z} \over {\partial x}}\right)^2\right) +{{h^2} \over z}\,\left(1-\left({{\partial z} \over {\partial x}}\right)^2\,\left({{\partial z} \over {\partial h}}\right)^2\right)\;.\label{eq:zeikonal} \end{eqnarray}Rewritten in terms of the vertical traveltime $\tau = z/v$, it furthertransforms to equation \begin{equation}{{\partial \tau} \over {\partial v}} = v\,\tau\,\left({{\partial \tau} \over {\partial x}}\right)^2 +{{h^2} \over {v^3\,\tau}}\,\left(1 - v^4\,\left({{\partial \tau} \over {\partial x}}\right)^2\,\left({{\partial \tau} \over {\partial h}}\right)^2\right)\;,\label{eq:eikonal0} \end{equation}equivalent to equation~(\ref{eq:eikonal}) in the main text. Yetanother form of the kinematic velocity continuation equation followsfrom eliminating the reflection angle $\theta$ from equations(\ref{eq:dzdv}) and (\ref{eq:z2h}). The resultant expression takes thefollowing form:\begin{equation}v\,{{\partial z} \over {\partial v}} = {{2\,(z^2 + h^2)} \over{\sqrt{z^2 + h^2 \sin^2{2\,\alpha}} + z\,\cos{2\,\alpha}}} ={z \over {\cos^2{\alpha}}} + {{2\,h^2} \over{\sqrt{z^2 + h^2 \sin^2{2\,\alpha}} + z}}\;.\label{eq:notheta} \end{equation}\append{Derivation of the residual DMO kinematics}This appendix derives the kinematical laws for the residual NMO+DMOtransformation in the prestack offset continuation process.The direct solution of equation (\ref{eq:DMONMOeikonal}) isnontrivial. A simpler way to obtain this solution is to decomposeresidual NMO+DMO into three steps and to evaluate their contributionsseparately. Let the initial data be the zero-offset reflection event$\tau_0(x_0)$. The first step of the residual NMO+DMO is the inverseDMO operator. One can evaluate the effect of this operator by means ofthe offset continuation concept \cite[]{ofcon}. According to thisconcept, each point of the input traveltime curve $\tau_0(x_0)$travels with the change of the offset from zero to $h$ along a specialtrajectory, which I call a {\em time ray}. Time rays are paraboliccurves of the form\begin{equation}x\left(\tau\right) =  x_0+{{\tau^2-\tau_0^2\left(x_0\right)} \over{\tau_0\left(x_0\right)\,\tau_0'\left(x_0\right)}}\;,\label{eq:xoftau}\end{equation}with the final points constrained by the equation\begin{equation}h^2 = \tau^2\,{{\tau^2-\tau_0^2\left(x_0\right)} \over{\left(\tau_0\left(x_0\right)\,\tau_0'\left(x_0\right)\right)^2}}\;,\label{eq:hoftau}\end{equation}where $\tau_0'\left(x_0\right)$ is the derivative of $\tau_0\left(x_0\right)$.The second step of the cumulative residual NMO+DMO process is theresidual normal moveout. According to equation (\ref{eq:ResNMO}), residualNMO is a one-trace operation transforming the traveltime $\tau$ to$\tau_1$ as follows:\begin{equation}\tau_1^2 = \tau^2 + h^2\,d\;,\label{eq:ResNMO2} \end{equation}where\begin{equation}d = \left({1 \over v_0^2} - {1 \over v_1^2}\right)\;.\label{eq:formers}\end{equation}The third step is dip moveout corresponding to the new velocity$v_1$. DMO is the offset continuation from $h$ to zerooffset along the redefined time rays \cite[]{ofcon}\begin{equation}x_2\left(\tau_2\right) =  x + {{h\,X} \over {\tau_1^2\,H}}\,\left(\tau_1^2-\tau_2^2\right)\;,\label{eq:xoftau2}\end{equation}where $H = {{\partial \tau_1} \over {\partial h}}$, and $X ={{\partial \tau_1} \over {\partial x}}$. The end points of the time rays (\ref{eq:xoftau2}) are defined by theequation\begin{equation}\tau_2^2 = - \tau_1^2\,{{\tau_1\,H} \over {h\,X^2}}\;.\label{eq:tauofh2}\end{equation}The partial derivatives of the common-offset traveltimes areconstrained by the offset continuation kinematic equation\begin{equation}h\,(H^2 - X^2) = \tau_1\,H\;,\label{eq:OCequation}\end{equation}which is equivalent to equation (\ref{eq:OCeikonal}) in AppendixA. Additionally, as follows from equations (\ref{eq:ResNMO2}) and the rayinvariant equations from \cite[]{ofcon},\begin{equation}\tau_1\,X = \tau\,{{\partial \tau} \over {\partial x}} = {{\tau^2\,\tau_0'\left(x_0\right)} \over {\tau_0\left(x_0\right)}}\;.\label{eq:invariant}\end{equation}Substituting~(\ref{eq:xoftau}-\ref{eq:formers}) and(\ref{eq:OCequation}-\ref{eq:invariant}) into equations(\ref{eq:xoftau2}) and (\ref{eq:tauofh2}) and performing the algebraicsimplifications yields the parametric expressions for velocityrays of the residual NMO+DMO process:\begin{equation}\left\{\begin{array}{rcl}x_2(d) & = & \displaystyle{x_0 + {{h^2\,\tau_0'(x_0)} \over T}\,\left(1 - {T^2 \over T_2^2(d)}\right)}\;,\\ \\\tau(d) & = & \displaystyle{{{\tau_1^2(d)} \over {T_2(d)}}}\;,\end{array}\right.\label{eq:ResDMO2}\end{equation}where the function$T\left(h,\tau_0(x_0),\tau_0'\left(x_0\right)\right)$ is defined by\begin{equation}T\left(h,\tau,\tau_x\right) = {{\tau + \sqrt{\tau^2 + 4\,h^2\,\tau_x^2}} \over 2}\;,\label{eq:CapT}\end{equation}\begin{equation}T_2(d) = \sqrt{T\left(h,\tau_1^2(d),\tau_0'\left(x_0\right)\,T\left(h,\tau_0(x_0),\tau_0'\left(x_0\right)\right)\right)}\;,\label{eq:CapT2}\end{equation}and\begin{equation}\tau_1^2(d) = \tau_0\,T + d\,h^2\;.\label{eq:tauofs}\end{equation}The last step of the cascade of inverse DMO, residual NMO, and DMO isillustrated in Figure \ref{fig:vlcvoc}. The three plots in the figure showthe offset continuation to zero offset of the inverse DMO impulseresponse shifted by the residual NMO operator. The middle plotcorresponds to zero NMO shift, for which the DMO step collapses thewavefront back to a point.  Both positive (top plot) and negative(bottom plot) NMO shifts result in the formation of the specifictriangular impulse response of the residual NMO+DMO operator. Asnoticed by \cite{Etgen.sepphd.68}, the size of the triangularoperators dramatically decreases with the time increase.  For largetimes (pseudo-depths) of the initial impulses, the operator collapsesto a point corresponding to the pure NMO shift.\inputdir{Math}\sideplot{vlcvoc}{width=0.9\textwidth}{Kinematic residual NMO+DMO  operators constructed by the cascade of inverse DMO, residual NMO,  and DMO. The impulse response of inverse DMO is shifted by the  residual NMO procedure. Offset continuation back to zero offset  forms the impulse response of the residual NMO+DMO operator. Solid  lines denote traveltime curves; dashed lines denote the offset  continuation trajectories (time rays). Top plot: $v_1/v_0 = 1.2$.  Middle plot: $v_1/v_0 = 1$; the inverse DMO impulse response  collapses back to the initial impulse. Bottom plot: $v_1/v_0 = 0.8$.  The half-offset $h$ in all three plots is 1 km.}\append{INTEGRAL VELOCITY CONTINUATION AND KIRCHHOFF MIGRATION}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%The main goal of this appendix is to prove the equivalence between theresult of zero-offset velocity continuation from zero velocity andconventional post-stack migration. After solving the velocitycontinuation problem in the frequency domain, I transform the solutionback to the time-and-space domain and compare it with the conventionalKirchhoff migration operator \cite[]{GEO43-01-00490076}. The frequency-domainsolution has its own value, because it forms the basis for an efficientspectral algorithm for velocity continuation \cite[]{second}.  Zero-offset migration based on velocity continuation is the solutionof the boundary problem for equation (\ref{eq:POMequation2}) with theboundary condition\begin{equation}\left.P\right|_{v=0} = P_0\;,\label{eq:POMbound} \end{equation}where $P_0(t_0,x_0)$ is the zero-offset seismic section, and$P(t,x,v)$ is the continued wavefield. In order to find the solutionof the boundary problem composed of (\ref{eq:POMequation2}) and(\ref{eq:POMbound}), it is convenient to apply the functiontransformation $R(t,x,v) = t\,P(t,x,v)$, the time coordinatetransformation $\sigma = t^2/2$, and, finally, the double Fouriertransform over the squared time coordinate $\sigma$ and the spatialcoordinate $x$:\begin{equation}\widehat{R}(v) = \int \int\,P(t,x,v)\,\exp(i \Omega \sigma - i k x )\,t^2\,dt\,dx\;.\label{eq:FTK} \end{equation}With the change of domain, equation (\ref{eq:POMequation2}) transformsto the ordinary differential equation\begin{equation}{{d\,\widehat{R}} \over {d\,v}} = i\,{k^2 \over \Omega}\,v\,\widehat{R}\;,\label{eq:ODE} \end{equation}and the boundary condition (\ref{eq:POMbound}) transforms to the initialvalue condition\begin{equation}\widehat{R}(0) = \widehat{R}_0\;, \label{eq:ODEbound} \end{equation}where \begin{equation}\widehat{R}_0 = \int \int\,P_0(t_0,x_0)\,\exp(i \Omega \sigma_0 - i k x_0 )\,t_0^2\,dt_0\,dx_0\;,\label{eq:FTK0}\end{equation}and $\sigma_0 = t_0^2/2$.  The unique solution of the initial value(Cauchy) problem (\ref{eq:ODE}) - (\ref{eq:ODEbound}) is easily found to be\begin{equation}\widehat{R}(v) = \widehat{R}_0\,\exp\left( i\,{{k^2} \over {2\,\Omega}}\,v^2\right)\;.\label{eq:ODEsolution} \end{equation}In the transformed domain, velocity continuation appears to be a unitaryphase-shift operator. An immediate consequence of this remarkable fact is thecascaded migration decomposition of post-stack migration\cite[]{GEO52-05-06180643}:\begin{equation}\exp\left( i\,{{k^2} \over {2\,\Omega}}\,(v_1^2 +  \cdots + v_n^2)\right) =\exp\left( i\,{{k^2} \over {2\,\Omega}}\,v_1^2\right)\,\cdots\,\exp\left( i\,{{k^2} \over {2\,\Omega}}\,v_n^2\right)\;.\label{eq:cascaded} \end{equation}Analogously, three-dimensional post-stack migration is decomposedinto the two-pass procedure \cite[]{GPR31-01-00340056}:\begin{equation}\exp\left( i\,{{k_1^2+k_2^2} \over {2\,\Omega}}\,v^2\right) =\exp\left( i\,{{k_1^2} \over {2\,\Omega}}\,v^2\right)\,\exp\left( i\,{{k_2^2} \over {2\,\Omega}}\,v^2\right)\;.\label{eq:two-pass}\end{equation}The inverse double Fourier transform of both sides of equality(\ref{eq:ODEsolution}) yields the integral (convolution) operator\begin{equation}P(t,x,v) = \int\int\,P_0(t_0,x_0)\,K(t_0,x_0;t,x,v)\,dt_0\,dx_0\;,\label{eq:convolution}\end{equation}with the kernel $K$ defined by\begin{equation}K = {{t_0^2/t} \over {(2\,\pi)^{m+1}}}\,\int\int\,\exp\left(i\,{{k^2} \over {2\,\Omega}}\,v^2 + ik\,(x - x_0) - {{i\Omega} \over 2}\,(t^2 - t_0^2)\right)\,dk\,d\Omega\;,\label{eq:kernel}\end{equation}where $m$ is the number of dimensions in $x$ and $k$ ($m$ equals $1$or $2$). The inner integral on the wavenumber axis $k$ in formula(\ref{eq:kernel}) is a known table integral \cite[]{grad}. Evaluating thisintegral simplifies equation (\ref{eq:kernel}) to the form\begin{equation}K = {{t_0^2/t} \over {(2\,\pi)^{m/2+1}\,v^m}}\,\int\,(i\Omega)^{m/2}\,\exp\left[{{i\Omega} \over 2}\,\left(t_0^2 - t^2 - {{(x - x_0)^2} \over v^2}\right)\right]\,d\Omega\;.\label{eq:skernel}\end{equation}The term $(i\Omega)^{m/2}$ is the spectrum of the anti-causalderivative operator ${d \over {d\sigma}}$ of the order $m/2$. Notingthe equivalence\begin{equation}\left({\partial \over {\partial \sigma}}\right)^{m/2} =\left({1 \over t}\,{\partial \over {\partial t}}\right)^{m/2} =\left({1 \over t}\right)^{m/2}\,\left({\partial \over {\partial t}}\right)^{m/2}\;,\label{eq:halfdif}\end{equation}which is exact in the 3-D case ($m=2$) and asymptotically correct inthe 2-D case ($m=1$), and applying the convolution theoremtransforms operator (\ref{eq:convolution}) to the form\begin{equation}P(t,x,v) = {1 \over {(2\,\pi)^{m/2}}}\,\int\,{{\cos{\alpha}} \over {(v\,\rho)^{m/2}}}\,\left(- {\partial \over {\partial t_0}}\right)^{m/2}P_0\left({\rho \over v},x_0\right)\,dx_0\;,\label{eq:Kirchhoff}\end{equation}where $\rho = \sqrt{v^2\,t^2 + (x - x_0)^2}$, and $\cos{\alpha} =t_0/t$. Operator (\ref{eq:Kirchhoff}) coincides with the Kirchhoff operatorof conventional post-stack time migration \cite[]{GEO43-01-00490076}.\newpage%%% Local Variables: %%% mode: latex%%% TeX-master: t%%% TeX-master: t%%% TeX-master: t%%% TeX-master: t%%% TeX-master: t%%% TeX-master: t%%% TeX-master: t%%% End: