paper.tex

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

  \label{eqn:conv}
  \psi_k (x) = \beta (x-k)\;.
\end{equation}
In other words, the basis is constructed by integer shifts of a single
function $\beta(x)$. Substituting expression~(\ref{eqn:conv}) into
equation~(\ref{eq:basis}) yields
\begin{equation}\label{eqn:beta}
  f (x) = \sum_{k \in K} c_k \beta (x - k)\;.
\end{equation}
Evaluating the function $f(x)$ in equation~(\ref{eqn:beta}) at an
integer value $n$, we obtain the equation
\begin{equation}\label{eqn:dig}
  f (n) = \sum_{k \in K} c_k \beta (n-k)\;,
\end{equation}
which has the exact form of a discrete convolution. The basis function
$\beta(x)$, evaluated at integer values, is digitally convolved with
the vector of basis coefficients to produce the sampled values of the
function $f(x)$. We can invert equation~(\ref{eqn:dig}) to obtain the
coefficients $c_k$ from $f(n)$ by inverse recursive filtering
(deconvolution). In the case of a non-causal filter $\beta(n)$, an 
appropriate spectral factorization will be
needed prior to applying the recursive filtering.
\par

\inputdir{XFig}

According to the convolutional basis idea, forward interpolation
becomes a two-step procedure. The first step is the direct inversion
of equation~(\ref{eqn:dig}): the basis coefficients $c_k$ are found by
deconvolving the sampled function $f(n)$ with the factorized filter
$\beta(n)$. The second step reconstructs the continuous (or arbitrarily
sampled) function $f(x)$ according to formula~(\ref{eqn:beta}). The
two steps could be combined into one, but usually it is more
convenient to apply them separately. I show a schematic relationship
among different variables in Figure~\ref{fig:scheme}.

\sideplot{scheme}{height=1.5in}{Schematic relationship among
  different variables for interpolation with a convolutional basis.}

\subsection{B-splines}

\inputdir{Math}

B-splines represent a particular example of a convolutional basis.
Because of their compact support and other attractive numerical
properties, B-splines are a good choice of the basis set for the
forward interpolation problem and related signal processing problems
\cite[]{unser}.  According to \cite{handbook}, they exhibit superior
performance for any given order of accuracy in comparison with other
methods of similar efficiency.
\par
B-splines of the order 0 and 1 coincide with the nearest neighbor and
linear interpolants~(\ref{eq:bin}) and~(\ref{eq:lin}) respectively.
B-splines $\beta^n(x)$ of a higher order $n$ can be defined by a
repetitive convolution of the zeroth-order spline $\beta^0(x)$ (the
box function) with itself:
\begin{equation}
\label{eqn:bdef}
\beta^n(x) = 
\underbrace{\beta^0(x) \ast \cdots \ast \beta^0(x)}_{(n+1)\quad 
\mbox{times}}\;.
\end{equation}
There is also the explicit expression
\begin{equation}
\label{eqn:expl}
\beta^n(x) = 
\frac{1}{n!}\,\sum_{k=0}^{n+1} C_k^{n+1} (-1)^k 
(x + \frac{n+1}{2} - k)_{+}^n\;,
\end{equation}
which can be proved by induction. Here $C_k^{n+1}$ are the binomial
coefficients, and the function $x_{+}$ is defined as follows:
\begin{equation}
  \label{eqn:xp}
  x_{+} = \left\{\begin{array}{lcr}
x, & \mbox{for} & x > 0 \\
0, & \mbox{otherwise} &
\end{array}\right.
\end{equation}
As follows from formula~(\ref{eqn:expl}), the most commonly used cubic
B-spline $\beta^3(x)$ has the expression
\begin{equation}
  \label{eqn:beta3}
  \beta^3(x) = \left\{\begin{array}{lcr}
\displaystyle \left(4-6 |x|^2+3 |x|^3\right)/6, & 
\mbox{for} & 1 > |x| \geq 0 \\
\displaystyle (2-|x|)^3/6, & \mbox{for} & 2 > |x| \geq 1 \\
0, & \mbox{elsewhere} &
\end{array}\right.
\end{equation}
The corresponding discrete filter $\beta^3(n)$ is a centered 3-point
filter with coefficients 1/6, 2/3, and 1/6. According to the
traditional method, deconvolution with this filter is performed as a
tridiagonal matrix inversion \cite[]{deBoor}. One can, however,
accomplish the same task more efficiently by spectral factorization
and recursive filtering \cite[]{unser1}. The recursive filtering
approach generalizes straightforwardly to B-splines of higher orders.
\par
Both the support length and the smoothness of B-splines increase with
the order. In the limit, B-splines converge to the Gaussian function.
Figures~\ref{fig:splint3} and~\ref{fig:splint7} show the third- and
seventh-order splines $\beta^3(x)$ and $\beta^7(x)$, respectively, and
their continuous spectra.

\plot{splint3}{width=6in}{Third-order B-spline $\beta^3(x)$ (left)
  and its spectrum (right).}

\plot{splint7}{width=6in}{Seventh-order B-spline $\beta^7(x)$ (left)
  and its spectrum (right).}
\par
It is important to realize the difference between B-splines and the
corresponding interpolants $W(x,n)$, which are sometimes called
\emph{cardinal splines}.  An explicit computation of the cardinal
splines is impractical, because they have infinitely long support.
Typically, they are constructed implicitly by the two-step
interpolation method outlined above. The cardinal splines of orders 3
and 7 and their spectra are shown in Figures~\ref{fig:crdint3}
and~\ref{fig:crdint7}. As B-splines converge to the Gaussian function,
the corresponding interpolants rapidly converge to the sinc
function~(\ref{eq:w-cft}). Good convergence is achieved with the help
of the implicitly-generated long support, which results from
recursive filtering at the first step of the interpolation procedure.

\plot{crdint3}{width=6in}{Effective third-order B-spline interpolant
  (left) and its spectrum (right).}

\plot{crdint7}{width=6in}{Effective seventh-order B-spline interpolant
  (left) and its spectrum (right).}

\inputdir{chirp}

In practice, the recursive filtering step adds only marginally to the
total interpolation cost. Therefore, an $n$-th order B-spline
interpolation is comparable in cost with any other method that uses an
$(n+1)$-point interpolant. The comparison in accuracy usually turns
out in favor of B-splines. Figures~\ref{fig:cubspl4}
and~\ref{fig:kaispl8} compare interpolation errors of B-splines and
other similar-cost methods on the example from Figure~\ref{fig:chirp}.

\sideplot{cubspl4}{height=2.5in}{Interpolation error of the
  cubic-convolution interpolant (dashed line) compared to that of the
  third-order B-spline (solid line).}

\sideplot{kaispl8}{height=2.5in}{Interpolation error of the 8-point
  windowed sinc interpolant (dashed line) compared to that of the
  seventh-order B-spline (solid line).}
\par
Similarly to the comparison in Figures~\ref{fig:speclincub}
and~\ref{fig:speccubkai}, we can also compare the discrete responses
of B-spline interpolation with those of other methods. The right plots
in Figures~\ref{fig:speccubspl4} and~\ref{fig:speckaispl8} show that the
discrete spectra of the effective B-spline interpolants are genuinely
flat at low frequencies and wider than those of the competitive
methods. Although the B-spline responses are infinitely long because
of the recursive filtering step, they exhibit a fast amplitude decay.

\sideplot{speccubspl4}{height=2.5in}{Discrete interpolation responses
  of cubic convolution and third-order B-spline interpolants (left)
  and their discrete spectra (right) for $x=0.7$.}

\sideplot{speckaispl8}{height=2.5in}{Discrete interpolation responses of
  8-point windowed sinc and seventh-order B-spline interpolants (left)
  and their discrete spectra (right) for $x=0.7$.}

\subsection{2-D example}

\inputdir{chirp2}

For completeness, I include a 2-D forward interpolation example.
Figure~\ref{fig:chirp2} shows a 2-D analog of the function in
Figure~\ref{fig:chirp} and its coarsely-sampled version.

\plot{chirp2}{width=6in,height=3in}{Two-dimensional test function
  (left) and its coarsely sampled version (right).}
\par
Figure~\ref{fig:plcbinlin} compares the errors of the 2-D nearest
neighbor and 2-D linear (bi-linear) interpolation. Switching to
bi-linear interpolation shows a significant improvement, but the error
level is still relatively high. As shown in
Figures~\ref{fig:plccubspl} and~\ref{fig:plckaispl}, B-spline
interpolation again outperforms other methods with comparable cost. In
all cases, I constructed 2-D interpolants by orthogonal splitting.
Although the splitting method reduces computational overhead, the main
cost factor is the total interpolant size, which is squared when the
interpolation goes from one to two dimensions.

\plot{plcbinlin}{width=5.5in,height=3.2in}{2-D Interpolation errors of
  nearest neighbor interpolation (left) and linear interpolation
  (right).  The top graphs show 1-D slices through the center of the
  image. Bi-linear interpolation exhibits smaller error and therefore
  is more accurate.}

\plot{plccubspl}{width=5.5in,height=3.2in}{2-D Interpolation errors of
  cubic convolution interpolation (left) and third-order B-spline
  interpolation (right).  The top graphs show 1-D slices through the
  center of the image. B-spline interpolation exhibits smaller error
  and therefore is more accurate.}

\plot{plckaispl}{width=5.5in,height=3.2in}{2-D Interpolation errors of
  8-point windowed sinc interpolation (left) and seventh-order
  B-spline interpolation (right).  The top graphs show 1-D slices
  through the center of the images. B-spline interpolation exhibits
  smaller error and therefore is more accurate.}

\subsection{Beyond B-splines}

It is not too difficult to construct a convolutional basis with more
accurate interpolation properties than those of B-splines, for example
by sacrificing the function smoothness. The following piece-wise cubic
function has a lower smoothness than $\beta^3(x)$ in
equation~(\ref{eqn:beta3}) but slightly better interpolation behavior:
\begin{equation}
  \label{eqn:mu3}
  \mu^3(x) = \left\{\begin{array}{lcr}
\displaystyle \left(10-13 |x|^2+6 |x|^3\right)/16, & 
\mbox{for} & 1 > |x| \geq 0 \\
\displaystyle (2-|x|)^2 (5-2 |x|)/16, & \mbox{for} & 2 > |x| \geq 1 \\
0, & \mbox{elsewhere} &
\end{array}\right.
\end{equation}
\par

\begin{comment}
Figures \ref{fig:splmom4} and \ref{fig:specsplmom4} compare the test
interpolation errors and discrete responses of methods based on the
B-spline function $\beta^3(x)$ and the lower smoothness function
$\mu^3(x)$. The latter method has a slight but visible performance
advantage and a slightly wider discrete spectrum.

%\sideplot{splmom4}{height=2.5in}{Interpolation error of the
  third-order B-spline interpolant (dashed line) compared to that of
  the lower smoothness spline interpolant (solid line).}
%\sideplot{specsplmom4}{height=2.5in}{Discrete interpolation responses
  of third-order B-spline and lower smoothness spline interpolants
  (left) and their discrete spectra (right) for $x=0.7$. A slight but
  visible difference in the interpolation responses accounts for a
  small improvement in accuracy.}
\end{comment}

\par
\cite{mom} have developed a general approach for constructing
non-smooth piece-wise functions with optimal interpolation properties.
However, the gain in accuracy is often negligible in practice. In the
rest of the dissertation, I use the classic and better tested B-spline
method.

\section{Seismic applications of forward interpolation}

\inputdir{stolt}

For completeness, I conclude this section with two simple examples of
forward interpolation in seismic data processing.
Figure~\ref{fig:stolt} shows a 3-D impulse response of Stolt migration
\cite[]{GEO43-01-00230048}, computed by using 2-point linear
interpolation and 8-point B-spline interpolation. As noted by
\cite{Ronen.sep.30.95} and \cite{Harlan.sep.30.103},
inaccurate interpolation may lead to spurious artifact events in
Stolt-migrated images. Indeed, we see several artifacts in the image
with linear interpolation (the left plots in Figure~\ref{fig:stolt}).
The artifacts are removed if we use a more accurate interpolation
method (the right plots in Figure~\ref{fig:stolt}).

\plot{stolt}{width=5.5in}{Stolt-migration impulse response. Left: using
  linear interpolation. Right: using seventh-order B-spline
  interpolation. Migration artifacts are removed by a more accurate
  forward interpolation method.}

\inputdir{radial}

\par
Another simple example is the radial trace transform
\cite[]{Ottolini.sepphd.33}. Figure~\ref{fig:radialdat} shows a land
shot gather contaminated by nearly radial ground-roll. As discussed by
\cite{Claerbout.sep.35.43}, \cite{SEG-1999-12041207,SEG-2000-21112114}, and
\cite{Brown.sep.103.morgan1,SEG-2000-21152118}, one can effectively
eliminate ground-roll noise by applying a radial trace transform
followed by high-pass filtering and the inverse radial transform.
Figure~\ref{fig:radial} shows the result of the forward radial
transform of the shot gather in Figure~\ref{fig:radialdat} in the
radial band of the ground-roll noise and the transform error after we
go back to the original domain. Comparing the results of using linear
and third-order B-spline interpolation, we see once again that the
transform artifacts are removed with a more accurate interpolation
scheme.

\sideplot{radialdat}{height=2.5in}{Ground-roll-contaminated shot
  gather used in a radial transform test}

\plot{radial}{width=5.5in}{Radial trace transform results. Top: radial
  trace domain. Bottom: residual error after the inverse transform.
  The error should be zero in a radial band from 0 to 0.65 km/s radial
  velocity. Left: using linear interpolation. Right: using third-order
  B-spline interpolation.}

\section{Acknowledgments}

A short conversation with Dave Hale led me to a better understanding
of different forward interpolation methods. Tamas Nemeth helped me 
better understand the general interpolation theory.

\bibliographystyle{segnat}
\bibliography{SEG,SEP2,forwd}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: