\title{Seismic reflection data interpolation with differential offsetand shot continuation}\author{Sergey Fomel}\maketitle\begin{abstract}I propose a finite-difference offset continuation filter for  interpolating seismic reflection data. The filter is constructed  from the offset continuation differential equation and is applied on  frequency slices in the log-stretch frequency domain.  Synthetic and  real data tests demonstrate that the proposed method succeeds in  structurally complex situations where more simplistic approaches  fail.\end{abstract}\section{Introduction}Data interpolation is one of the most important problems of seismicdata processing. In 2-D exploration, the interpolation problem arisesbecause of missing near and far offsets, spatial aliasing andoccasional bad traces. In 3-D exploration, the importance of thisproblem increases dramatically because 3-D acquisition almost neverprovides a complete regular coverage in both midpoint and offsetcoordinates \cite[]{Biondi.3dsi.99}. Data regularization in 3-D cansolve the problem of Kirchoff migration artifacts \cite[]{SEG-1994-1553},prepare the data for wave-equation common-azimuth imaging\cite[]{GEO61-06-18221832}, or provide the spatial coverage required for 3-Dmultiple elimination \cite[]{SEG-1998-1321}.\par\cite{Claerbout.blackwell.92,gee} formulates the following generalprinciple of missing data interpolation:\begin{quote}  A method for restoring missing data is to ensure that the restored  data, after specified filtering, has minimum energy.\end{quote}How can one specify an appropriate filtering for a given interpolationproblem? Smooth surfaces are conveniently interpolated with Laplacianfilters \cite[]{GEO39-01-00390048}. Steering filters help usinterpolate data with predefined dip fields \cite[]{SEG-1998-1851}.Prediction-error filters in time-space or frequency-space domainsuccessfully interpolate data composed of distinctive plane waves\cite[]{GEO56-06-07850794,gee}. Local plane waves are handled withplane-wave destruction filters \cite[]{GEO67-06-19461960}.  Becauseprestack seismic data is not stationary in the offset direction,non-stationary prediction-error filters need to be estimated, whichleads to an accurate but relatively expensive method with manyadjustable parameters \cite[]{SEG-1999-11541157}.  \par A simple modelfor reflection seismic data is a set of hyperbolic events on a commonmidpoint gather. The simplest filter for this model is the firstderivative in the offset direction applied after the normal moveoutcorrection. Going one step beyond this simple approximation requirestaking the dip moveout (DMO) effect into account\cite[]{FBR04-07-00070024}. The DMO effect is fully incorporated inthe offset continuation differential equation\cite[]{me,GEO68-02-07180732}.  \par Offset continuation is a processof seismic data transformation between different offsets\cite[]{GPR29-03-03740406,GPR30-06-08130828,GPR30-06-08290849}.Different types of DMO operators \cite[]{DMP00-00-01300130} can beregarded as continuation to zero offset and derived as solutions of aninitial-value problem with the revised offset continuation equation\cite[]{GEO68-02-07180732}. Within a constant-velocity assumption, this equationnot only provides correct traveltimes on the continued sections, butalso correctly transforms the corresponding wave amplitudes\cite[]{SEG-1996-1731}. Integral offset continuation operators havebeen derived independently by \cite{Chemingui.sep.82.117},\cite{GEO61-06-18461858}, and \cite{stovas}.  The 3-D analog is knownas azimuth moveout (AMO) \cite[]{GEO63-02-05740588}. In theshot-record domain, integral offset continuation transforms to shotcontinuation \cite[]{Schwab.sep.77.117,SEG-1993-0673,SEG-1996-0439}.Integral continuation operators can be applied directly for missingdata interpolation and regularization\cite[]{SEG-1994-1549,SEG-1999-19951998}. However, they don't behavewell for continuation at small distances in the offset space becauseof limited integration apertures and, therefore, are not well suitedfor interpolating neighboring records. Additionally, as all integral(Kirchoff-type) operators they suffer from irregularities in the inputgeometry. The latter problem is addressed by accurate but expensiveinversion to common offset \cite[]{Chemingui.sepphd.101}.  \par Inthis paper, I propose an application of offset continuation in theform of a finite-difference filter for Claerbout's method of missingdata interpolation.  The filter is designed in the log-stretchfrequency domain, where each frequency slice can be interpolatedindependently.  Small filter size and easy parallelization amongdifferent frequencies assure a high efficiency of the proposedapproach. Although the offset continuation filter lacks the predictivepower of non-stationary prediction-error filters, it is much simplerto handle and serves as a good \emph{a priori} guess of aninterpolative filter for seismic reflection data. I first test theproposed method by interpolating randomly missing traces in aconstant-velocity synthetic dataset. Next, I apply offset continuationand related shot continuation field to a real data example from theNorth Sea. Using a pair of offset continuation filters, operating intwo orthogonal directions, I successfully regularize a 3-D marinedataset. These tests demonstrate that the offset continuation canperform well in complex structural situations where more simplisticapproaches fail.\section{Offset continuation}A particularly efficient implementation of offset continuation resultsfrom a log-stretch transform of the time coordinate\cite[]{GPR30-06-08130828}, followed by a Fourier transform of thestretched time axis. After these transforms, the offset continuationequation from \cite[]{GEO68-02-07180732} takes the form\begin{equation}  h \, \left( {\partial^2 \tilde{P} \over \partial y^2} -     {\partial^2 \tilde{P} \over \partial h^2} \right) -   i\,\Omega \, {\partial \tilde{P} \over   {\partial h}} = 0 \;,  \label{eqn:OC} \end{equation}where $\Omega$ is the corresponding frequency, $h$ is the half-offset,$y$ is the midpoint, and $\tilde{P} (y,h,\Omega)$ is the transformeddata. As in other $F$-$X$ methods, equation~(\ref{eqn:OC}) can beapplied independently and in parallel on different frequency slices.\parWe can construct an effective offset-continuation finite-differencefilter by studying first the problem of wave extrapolation betweenneighboring offsets. In the frequency-wavenumber domain, theextrapolation operator is defined by solving the initial-value problemon equation~(\ref{eqn:OC}). The solution takes the following form\cite[]{GEO68-02-07180732}:\begin{equation}\widehat{\widehat{P}}(h_2) = \widehat{\widehat{P}}(h_1)\,Z_{\lambda}(kh_2)/Z_{\lambda}(kh_1)\;,\label{eqn:OKOC}\end{equation}where $\lambda = (1 + i \Omega)/2$, and $Z_\lambda$ is the specialfunction defined as\begin{eqnarray}\nonumberZ_{\lambda}(x) & = & \Gamma(1-\lambda)\,\left(x \over 2\right)^{\lambda}\,J_{-\lambda}(x)={}_0F_1\left(;1-\lambda;-\frac{x^2}{4}\right) \\& = &\sum_{n=0}^{\infty} {(-1)^n \over n!}\,{\Gamma(1-\lambda) \over \Gamma(n+1-\lambda)}\,\left(x \over 2\right)^{2n}\;,\label{eqn:z}\end{eqnarray}where $\Gamma$ is the gamma function, $J_{-\lambda}$ is the Besselfunction, and ${}_0F_1$ is the confluent hypergeometric limit function\cite[]{ab}. The wavenumber $k$ in equation~(\ref{eqn:OKOC}) correspondsto the midpoint $y$ in the original data domain.  In thehigh-frequency asymptotics, operator~(\ref{eqn:OKOC}) takes the form\begin{equation}\widehat{\widehat{P}}(h_2) = \widehat{\widehat{P}}(h_1)\,F(2 k h_2/\Omega)/F(2 k h_1/\Omega)\,\exp{\left[i\Omega\,\psi\left(2 k h_2/\Omega - 2 k h_1/\Omega\right)\right]}\;,\label{eqn:AOKOC}\end{equation}where\begin{equation}F(\epsilon)=\sqrt{{1+\sqrt{1+\epsilon^2}} \over{2\,\sqrt{1+\epsilon^2}}}\,\exp\left({1-\sqrt{1+\epsilon^2}} \over 2\right)\;,\label{eqn:F}\end{equation}and\begin{equation}\psi(\epsilon)={1 \over 2}\,\left(1 - \sqrt{1+\epsilon^2} +\ln\left({1 + \sqrt{1+\epsilon^2}} \over 2\right)\right)\;.\label{eqn:psi}\end{equation}Returning to the original domain, we can approximate the continuationoperator with a finite-difference filter with the $Z$-transform\begin{equation}\label{eqn:OCpass}\hat{P}_{h+1}(Z_y) = \hat{P}_{h} (Z_y) \frac{G_1(Z_y)}{G_2(Z_y)}\;.\end{equation}The coefficients of the filters $G_1(Z_y)$ and $G_2(Z_y)$ are found byfitting the Taylor series coefficients of the filter response aroundthe zero wavenumber.  In the simplest case of 3-pointfilters\footnote{An analogous technique applied to the case of  wavefield depth extrapolation with the wave equation would lead to  the famous 45-degree implicit finite-difference operator  \cite[]{Claerbout.blackwell.85}.}, this procedure uses four Taylorseries coefficients and leads to the following expressions:\begin{eqnarray}  \label{eqn:OCnum}  G_1(Z_y) & = & 1 - \frac{1 - c_1(\Omega) h_2^2 + c_2(\Omega) h_1^2}{6} +  \frac{1 - c_1(\Omega) h_2^2 + c_2(\Omega) h_1^2}{12}\,  \left(Z_y + Z_y^{-1}\right)\;, \\  \label{eqn:OCden}  G_2(Z_y) & = & 1 - \frac{1 - c_1(\Omega) h_1^2 + c_2(\Omega) h_2^2}{6} +  \frac{1 - c_1(\Omega) h_1^2 + c_2(\Omega) h_2^2}{12}\,  \left(Z_y + Z_y^{-1}\right)\;,\end{eqnarray}where \[c_1(\Omega) = \frac{3\,(\Omega^2 + 9 - 4  i\,\Omega)}{\Omega^2\,(3+i\,\Omega)}\]and \[c_2(\Omega) =\frac{3\,(\Omega^2 - 27 - 8 i\,\Omega)}{\Omega^2\,(3+i\,\Omega)}\;.\]Figure~\ref{fig:arg} compares the phase characteristic of thefinite-difference extrapolators~(\ref{eqn:OCpass}) with the phasecharacteristics of the exact operator~(\ref{eqn:OKOC}) and theasymptotic operator~(\ref{eqn:AOKOC}). The match between differentphases is poor for very low frequencies (left plot inFigure~\ref{fig:arg}) but sufficiently accurate for frequencies in thetypical bandwidth of seismic data (right plot inFigure~\ref{fig:arg}).Figure~\ref{fig:off-imp} compares impulse responses of the inverse DMOoperator constructed by the asymptotic $\Omega-k$ operator with thoseconstructed by finite-difference offset continuation. Neglectingsubtle phase inaccuracies at large dips, the two images look similar,which provides an experimental evidence of the accuracy of theproposed finite-difference scheme.When applied on the offset-midpoint plane of an individual frequencyslice, the one-dimensional implicit filter~(\ref{eqn:OCpass})transforms to a two-dimensional explicit filter with the2-D $Z$-transform \begin{equation}\label{eqn:gfilt}G(Z_y,Z_h) = G_1(Z_y) - Z_h G_2(Z_y)\;.\end{equation}Convolution with filter~(\ref{eqn:gfilt}) is the regularizationoperator that I propose to use for interpolating prestack seismic data.%I propose to adopt a finite-difference form of the differential%operator~(\ref{eqn:OC}) for the regularization operator $\mathbf{D}$. A%simple analysis of equation~(\ref{eqn:OC}) shows that at small%frequencies, the operator is dominated by the first term. The form%${\partial^2 P \over \partial y^2} - {\partial^2 P \over \partial%  h^2}$ is equivalent to the second mixed derivative in the source and%receiver coordinates. Therefore, at low frequencies, the offset waves%propagate in the source and receiver directions. At high frequencies,%the second term in~(\ref{eqn:OC}) becomes dominating, and the entire%method becomes equivalent to the trivial linear interpolation in%offset. The interpolation pattern is more complicated at intermediate%frequencies.  \inputdir{Math}