paper.tex

国外免费地震资料处理软件包

\plot{sincint}{width=6in}{Sinc interpolant (left) and its spectrum
  (right).}

Function~(\ref{eq:w-cft}) is well-known as the Shannon sinc
interpolant. According to the sampling theorem
\cite[]{kotelnikov,shannon}, it provides an optimal interpolation for
band-limited signals.  A known problem prohibiting its practical
implementation is the slow decay with $(x - n)$, which results in a
far too expensive computation. This problem is solved in practice with
heuristic tapering \cite[]{Hale.sep.25.39}, such as triangle tapering
\cite[]{Harlan.sep.30.103}, or more sophisticated taper windows
\cite[]{wolberg}. One popular choice is the Kaiser window \cite[]{kaiser},
which has the form
\begin{equation}
  \label{eqn:kais}
  W (x, n) = \left\{\begin{array}{lcr} \displaystyle
      \frac{\sin \left[\pi (x - n) \right]}{\pi (x - n)}\,
      \frac{I_0\left(a\,\sqrt{1-\left(\frac{x-n}{N}\right)^2}\right)}
      {I_0(a)} & \mbox{for} & n - N < x < n + N \\
      0, & \mbox{otherwise} &
\end{array}\right.
\end{equation}
where $I_0$ is the zero-order modified Bessel function of the first
kind.  The Kaiser-windowed sinc interpolant~(\ref{eqn:kais}) has the
adjustable parameter $a$, which controls the behavior of its
spectrum. I have found empirically the value of $a=4$ to provide
a spectrum that deviates from 1 by no more than 1\% in a relatively
wide band.

While the function $W$ from equation (\ref{eq:w-cft}) automatically
satisfies properties (\ref{eq:sum}) and (\ref{eq:intwdx}), where both
$x$ and $n$ range from $-\infty$ to $\infty$, its tapered version may
require additional normalization.

\inputdir{chirp}

Figure~\ref{fig:cubkai} compares the interpolation error of the
8-point Kaiser-tapered sinc interpolant with that of cubic convolution
on the example from Figure~\ref{fig:chirp}. The accuracy improvement
is clearly visible.

\sideplot{cubkai}{height=2.5in}{Interpolation error of the
  cubic-convolution interpolant (dashed line) compared to that of an
  8-point windowed sinc interpolant (solid line).}

The differences among the described forward interpolation methods are
also clearly visible from the discrete spectra of the corresponding
interpolants.  The left plots in Figures~\ref{fig:speclincub}
and~\ref{fig:speccubkai} show discrete interpolation responses: the
function $W(x,n)$ for a fixed value of $x=0.7$. The right plots
compare the corresponding discrete spectra. Clearly, the spectrum gets
flatter and wider as the accuracy of the method increases.

\sideplot{speclincub}{height=2.5in}{Discrete interpolation responses
  of linear and cubic convolution interpolants (left) and their
  discrete spectra (right) for $x=0.7$.}

\sideplot{speccubkai}{height=2.5in}{Discrete interpolation responses
  of cubic convolution and 8-point windowed sinc interpolants (left)
  and their discrete spectra (right) for $x=0.7$.}

\subsection{Discrete Fourier basis}

\inputdir{Math}

Assuming that the range of the variable $x$ is limited in the interval
from $-N$ to $N$, the discrete Fourier basis (\emph{Fast Fourier
  Transform}) employs a set of orthonormal periodic functions
\begin{equation}\label{eq:fft}
  \psi_k (x) = \frac{1}{\sqrt{2N}} e^{i \pi \frac{k}{N} x} \;,
\end{equation}
where the discrete frequency index $k$ also ranges, according to the
Nyquist sampling criterion, from $-N$ to $N$. The interpolation
function is computed from equation~(\ref{eq:solution}) to be
\begin{eqnarray}\label{eq:w-fft}
  W (x, n) & = & \frac{1}{2 N} \sum_{k=-N}^{N-1} e^{i \pi \frac{k}{N} (x-n)} =
\frac{1}{2 N} e^{- i \pi (x-n)} \left[
  1 + e^{i \pi \frac{x-n}{N}} + \cdots +  e^{i \pi \frac{2N-1}{N} (x-n)}
\right] =  \nonumber \\
& & \frac{1}{2 N} e^{- i \pi (x-n)}
\frac{e^{2i \pi (x-n)} - 1}{e^{i \pi \frac{x-n}{N}} - 1} =
\frac{1}{2 N} e^{- i \pi \frac{x-n}{2 N}}
\frac{e^{i \pi (x-n)} - e^{-i \pi (x-n)}}{e^{i \pi \frac{x-n}{2 N}} -
e^{- i \pi \frac{x-n}{2 N}}} = \nonumber \\
& & e^{-i \pi \frac{x-n}{2 N}}
\frac{\sin \left[\pi (x - n)\right]}{2N \sin\left[\pi (x - n)/2N\right]} \;.
\end{eqnarray}
An interpolation function equivalent to (\ref{eq:w-fft}) has been
found by Muir \cite[]{Lin.sep.79.255,Popovici.sep.79.261,muir}. It can
be considered a tapered version of the sinc interpolant
(\ref{eq:w-cft}) with smooth tapering function
\[
\frac{\pi (x - n)/2N}{\tan\left[\pi (x - n)/2N\right]}\;.
\]
Unlike most other tapered-sinc interpolants, Muir's interpolant
(\ref{eq:w-fft}) satisfies not only the obvious property
(\ref{eq:wnn}), but also properties (\ref{eq:sum}) and
(\ref{eq:intwdx}), where the interpolation function $W (x,n)$ should
be set to zero for $x$ outside the range from $n - N$ to $n+N$. The
form of this function is shown in Figure \ref{fig:ma-sinc}.

\plot{ma-sinc}{width=6.0in}{The left plots show the sinc interpolation
  function.  Note the slow decay in $x$. The middle shows the
  effective tapering function of Muir's interpolation; the right is
  Muir's interpolant. The top is for $N=2$ (5-point interpolation);
  the bottom, $N=6$ (13-point interpolation).}

The development of the mathematical wavelet theory \cite[]{wavelet} has
opened the door to a whole universe of orthonormal function bases,
different from the Fourier basis. The wavelet theory should find many
useful applications in geophysical data interpolation, but exploring
this interesting opportunity would go beyond the scope of the present
work.

The next section carries the analysis to the continuum and compares
the mathematical interpolation theory with the theory of seismic
imaging.

\section{Continuous case and seismic imaging}
Of course, the linear theory is not limited to discrete grids.  It is
interesting to consider the continuous case because of its connection
to the linear integral operators commonly used in seismic imaging.
Indeed, in the continuous case, linear decomposition (\ref{eq:basis})
takes the form of the integral operator
\begin{equation}\label{eq:cont}
 f (y) = \int m (x) G (y; x) d x \;,
\end{equation}
where $x$ is a continuous analog of the discrete coefficient $k$ in
(\ref{eq:basis}), the continuous function $m (x)$ is analogous to the
coefficient $c_k$, and $G (y; x)$ is analogous to one of the basis
functions $\psi_k (x)$. The linear integral operator in
(\ref{eq:cont}) has a mathematical form similar to the form of
well-known integral imaging operators, such as Kirchhoff migration or
``Kirchhoff'' DMO. Function $G (y; x)$ in this case represents the
Green's function (impulse response) of the imaging operator.  Linear
decomposition of the data into basis functions means decomposing it
into the combination of impulse responses (``hyperbolas'').

In the continuous case, equation~(\ref{eq:solution}) transforms to
\begin{equation}\label{eq:wyn}
  W (y, n) = \int\!\!\int \Psi^{-1} (x_1, x_2) G (y;x_1) \bar{G} (n;x_2)
  dx_1\,dx_2\;,
\end{equation}
where $\Psi^{-1} (x_1, x_2)$ refers to the inverse of the ``matrix''
operator
\begin{equation}\label{eq:psi}
  \Psi (x_1, x_2) = \int G (y;x_1) \bar{G} (y;x_2) dy\;.
\end{equation}
When the linear operator, defined by equation~(\ref{eq:cont}), is
\emph{unitary},
\begin{equation}\label{eq:delta}
\Psi^{-1} (x_1, x_2) = \delta (x_1 - x_2)\;,
\end{equation}
and equation~(\ref{eq:wyn}) simplifies to the single integral
\begin{equation}\label{eq:wxn}
  W (y, n) = \int G (y;x) \bar{G} (n;x) dx \;.
\end{equation}
With respect to seismic imaging operators, one can recognize in the
interpolation operator (\ref{eq:wxn}) the generic form of azimuth moveout
\cite[]{Biondi.sep.93.15}, which is derived either as a cascade of adjoint
($\bar{G}(n;y)$) and forward ($G (x;y)$) DMO or as a cascade of migration
($\bar{G} (n;y)$) and modeling ($G (x;y)$)
\cite[]{Fomel.sep.84.25,GEO63-02-05740588}. In the first case, the
intermediate variable $y$ corresponds to the space of zero-offset data cube.
In the second case, it corresponds to a point in the subsurface.

\subsection{Asymptotically pseudo-unitary operators as orthonormal bases}

It is interesting to note that many integral operators routinely used
in seismic data processing have the form of operator (\ref{eq:cont})
with the Green's function
\begin{equation}
  \label{eq:green}
  G (t,\mathbf{y};z,\mathbf{x}) = \left|\frac{\partial}{\partial t}\right|^{m/2}
  A (\mathbf{x};t,\mathbf{y})
  \delta \left(z - \theta(\mathbf{x};t,\mathbf{y}) \right)\;.
\end{equation}
where we have split the variable $x$ into the one-dimensional
component $z$ (typically depth or time) and the $m$-dimensional
component $\mathbf{x}$ (typically a lateral coordinate with $m$ equal
$1$ or $2$). Similarly, the variable $y$ is split into $t$ and
$\mathbf{y}$.  The function $\theta$ represents the \emph{summation
  path}, which captures the kinematic properties of the operator, and
$A$ is the amplitude function. In the case of $m=1$, the fractional
derivative $\left|\frac{\partial}{\partial t}\right|^{m/2}$ is defined
as the operator with the frequency response $(i\,\omega)^{m/2}$, where
$\omega$ is the temporal frequency \cite[]{samko}.

The impulse response (\ref{eq:green}) is typical for different forms
of Kirchhoff migration and datuming as well as for velocity transform,
integral offset continuation, DMO, and AMO. Integral operators of that
class rarely satisfy the unitarity condition, with the Radon transform
(slant stack) being a notable exception. In an earlier paper
\cite[]{Fomel.sep.92.267}, I have shown that it is possible to define
the amplitude function $A$ for each kinematic path $\theta$ so that
the operator becomes \emph{asymptotically pseudo-unitary}. This means
that the adjoint operator coincides with the inverse in the
high-frequency (stationary-phase) approximation.  Consequently,
equation (\ref{eq:delta}) is satisfied to the same asymptotic order.

Using asymptotically pseudo-unitary operators, we can apply formula
(\ref{eq:wxn}) to find an explicit analytic form of the interpolation
function $W$, as follows:
\begin{eqnarray}
  \label{eq:apu}
  W (t, \mathbf{y}; t_n, \mathbf{y}_n) =  \int\!\!\int
  G (t, \mathbf{y}; z,\mathbf{x})\, G(t_n,\mathbf{y}_n;z,\mathbf{x})\,
  d z \, d \mathbf{x} =
  \nonumber \\
  \left|\frac{\partial}{\partial t}\right|^{m/2}
  \left|\frac{\partial}{\partial t_n}\right|^{m/2} \int
  A (\mathbf{x};t,\mathbf{y}) \, A (\mathbf{x};t_n,\mathbf{y}_n)\,
  \delta \left(\theta(\mathbf{x};t  ,\mathbf{y}  ) -
               \theta(\mathbf{x};t_n,\mathbf{y}_n) \right) \,
             d \mathbf{x}\;.
\end{eqnarray}
Here the amplitude function $A$ is defined according to the general
theory of asymptotically pseudo-inverse operators as
\begin{equation}
A = {\frac{1}{\left(2\,\pi\right)^{m/2}}}\,
{\left|F\,\widehat{F}\right|^{1/4}\,
\left|\frac{\partial \theta}{\partial t}\right|^{(m+2)/4}}\;,
\label{eq:weight}
\end{equation}
where
\begin{eqnarray}
F & = & \frac{\partial \theta}{\partial t}\,
\frac{\partial^2 \theta}{\partial \mathbf{x}\, \partial \mathbf{y}} -
\frac{\partial \theta}{\partial \mathbf{y}}\,
\frac{\partial^2 \theta}{\partial \mathbf{x}\, \partial t}\;,
 \\
\widehat{F} & = & \frac{\partial \widehat{\theta}}{\partial z}\,
\frac{\partial^2 \widehat{\theta}}{\partial \mathbf{x}\, \partial \mathbf{y}} -
\frac{\partial \widehat{\theta}}{\partial \mathbf{x}}\,
\frac{\partial^2 \widehat{\theta}}{\partial \mathbf{y}\, \partial z}\;,
\end{eqnarray}
and $\widehat{\theta} (\mathbf{x};t,\mathbf{y})$ is the \emph{dual}
summation path, obtained by solving equation $z=\theta(x;t,y)$ for $t$
(assuming that an explicit solution is possible).

For a simple example, let us consider the case of zero-offset time
migration with a constant velocity $v$. The summation path $\theta$
in this case is an ellipse
\begin{equation}
  \label{eq:ellips}
  \theta(\mathbf{x};t,\mathbf{y}) = \sqrt{t^2 -
    \frac{(\mathbf{x}-\mathbf{y})^2}{v^2}}\;,
\end{equation}
and the dual summation path $\widehat{\theta}$ is a hyperbola
\begin{equation}
  \label{eq:hyper}
  \widehat{\theta}(\mathbf{y};z,\mathbf{x}) = \sqrt{z^2 +
    \frac{(\mathbf{x}-\mathbf{y})^2}{v^2}}\;.
\end{equation}
The corresponding pseudo-unitary amplitude function is found from
formula (\ref{eq:weight}) to be \cite[]{Fomel.sep.92.267}
\begin{equation}
  A = {\frac{1}{\left(2\,\pi\right)^{m/2}}}\,
  {\frac{\sqrt{t/z}}{v^m z^{m/2}}}\;.
\label{eq:zomig}
\end{equation}
Substituting formula (\ref{eq:zomig}) into (\ref{eq:apu}), we derive
the corresponding interpolation function
\begin{equation}
  \label{eq:pumig}
  W (t, \mathbf{y}; t_n, \mathbf{y}_n)
  = \frac{1}{\left(2\,\pi\right)^{m}} \,
  \left|\frac{\partial}{\partial t}\right|^{m/2}
  \left|\frac{\partial}{\partial t_n}\right|^{m/2} \int
  \frac{\sqrt{t\,t_n}}{v^{2m} z^{m+1}}\,
  \delta (z - z_n) \,d \mathbf{x}\;,
\end{equation}
where $z = \theta(\mathbf{x};t,\mathbf{y})$, and $z_n =
\theta(\mathbf{x};t_n,\mathbf{y}_n)$. For $m=1$ (the two-dimensional
case), we can apply the known properties of the delta function to
simplify formula (\ref{eq:pumig}) further to the form
\begin{equation}
  W
  = {\frac{v}{\pi}}\,
  {\left|\frac{\partial}{\partial t}\right|^{1/2}
    \left|\frac{\partial}{\partial t_n}\right|^{1/2}
    \frac{\sqrt{t\,t_n}}{\sqrt{
	\left[(\mathbf{y}-\mathbf{y}_n)^2 - v^2 (t - t_n)^2\right]
	\left[v^2 (t + t_n)^2 - (\mathbf{y}-\mathbf{y}_n)^2\right]
  }}}\;.
  \label{eq:2dmig}
\end{equation}
The result is an interpolant for zero-offset seismic sections.  Like
the sinc interpolant in equation~(\ref{eq:w-cft}), which is based on
decomposing the signal into sinusoids, equation~(\ref{eq:2dmig}) is
based on decomposing the zero-offset section into hyperbolas. 

While opening a curious theoretical possibility, seismic imaging
interpolants have an undesirable computational complexity. Following
the general regularization framework of Chapter~\ref{chapter:funda}, I
shift the computational emphasis towards appropriately chosen
regularization operators discussed in Chapter~\ref{chapter:regul}.
For the forward interpolation method, all data examples in this
dissertation use either the simplest nearest neighbor and linear
interpolation or a more accurate B-spline method, described in the
next section.

\section{Interpolation with convolutional bases}

\cite{unser1} noticed that the basis function idea has an
especially simple implementation if the basis is convolutional and
satisfies the equation
\begin{equation}