paper.tex

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

\end{itemize}

\subsection{Factorization examples}

\inputdir{helix}

The first simple example of helical spectral factorization is shown in
Figure~\ref{fig:autowaves}. A minimum-phase factor is found by
spectral factorization of its autocorrelation. The result is
additionally confirmed by applying inverse recursive filtering, which
turns the filter into a spike (the rightmost plot in
Figure~\ref{fig:autowaves}.)

\plot{autowaves}{width=4.5in, height=4.5in}{Example of 2-D Wilson-Burg
  factorization. Top left: the input filter. Top right: its
  auto-correlation. Bottom left: the factor obtained by the Wilson-Burg
  method. Bottom right: the result of deconvolution.}

A practical example is depicted in Figure~\ref{fig:laplac}.  The
symmetric Laplacian operator is often used in practice for
regularizing smooth data. In order to construct a corresponding
recursive preconditioner, we factor the Laplacian autocorrelation
(the biharmonic operator) using the Wilson-Burg algorithm.
Figure~\ref{fig:laplac} shows the resultant filter. The minimum-phase
Laplacian filter has several times more coefficients than the original
Laplacian. Therefore, its application would be more expensive in a
convolution application. The real advantage follows from the
applicability of the minimum-phase filter for inverse filtering
(deconvolution). The gain in convergence from recursive filter
preconditioning outweighs the loss of efficiency from the longer
filter.  Figure~\ref{fig:thin42} shows a construction of the smooth
inverse impulse response by application of the $\mathbf{C} = \mathbf{P
  P}^T$ operator, where $\mathbf{P}$ is deconvolution with the
minimum-phase Laplacian. The application of $\mathbf{C}$ is equivalent
to a numerical solution of the biharmonic equation, discussed in the
next section.

\inputdir{laplac}

\plot{laplac}{width=4.5in, height=4.5in}{Creating a minimum-phase
  Laplacian filter.  Top left: Laplacian filter. Top right:
  its auto-correlation (bi-harmonic filter). Bottom left: factor obtained by the Wilson-Burg method
  (minimum-phase Laplacian). Bottom right: the result of deconvolution.}

\plot{thin42}{width=6in, height=2in}{2-D deconvolution with the
  minimum-phase Laplacian. Left: input. Center: output of
  deconvolution.  Right: output of deconvolution and adjoint
  deconvolution (equivalent to solving the biharmonic differential
  equation).}

\section{Application of spectral factorization: \newline
  Regularizing smooth data with splines in tension}

\hyphenation{Woi-now-sky}

The method of minimum curvature is an old and ever-popular approach
for constructing smooth surfaces from irregularly spaced data
\cite[]{GEO39-01-00390048}. The surface of minimum curvature corresponds
to the minimum of the Laplacian power or, in an alternative
formulation, satisfies the biharmonic differential equation.
Physically, it models the behavior of an elastic plate. In the
one-dimensional case, the minimum curvature method leads to the
natural cubic spline interpolation \cite[]{deBoor}. In the
two-dimensional case, a surface can be interpolated with biharmonic
splines \cite[]{sandwell} or gridded with an iterative finite-difference
scheme \cite[]{swain}.  We approach the gridding (data regularization)
problem with an iterative least-squares optimization scheme.
\par
In most of the practical cases, the minimum-curvature method produces
a visually pleasing smooth surface. However, in cases of large changes
in the surface gradient, the method can create strong artificial
oscillations in the unconstrained regions. Switching to lower-order
methods, such as minimizing the power of the gradient, solves the
problem of extraneous inflections, but also removes the smoothness
constraint and leads to gradient discontinuities
\cite[]{galilee}. A remedy, suggested by \cite{schweikert},
is known as \emph{splines in tension}. Splines in tension are
constructed by minimizing a modified quadratic form that includes a
tension term. Physically, the additional term corresponds to tension
in elastic plates \cite[]{timoshenko}. \cite{GEO55-03-02930305}
developed a practical algorithm of 2-D gridding with splines in
tension and implemented it in the popular GMT software
package.
\par
In this section, we develop an application of helical preconditioning
to gridding with splines in tension. We accelerate an iterative data
regularization algorithm by recursive preconditioning with
multidimensional filters defined on a helix \cite[]{GEO68-02-05770588}. The
efficient Wilson-Burg spectral factorization constructs a
minimum-phase filter suitable for recursive filtering.

We introduce a family of 2-D
minimum-phase filters for different degrees of tension.  The filters
are constructed by spectral factorization of the corresponding
finite-difference forms. In the case of zero tension (the original
minimum-curvature formulation), we obtain a minimum-phase version of
the Laplacian filter. The case of infinite tension leads to spectral
factorization of the Laplacian and produces the \emph{helical
  derivative} filter \cite[]{iee}.
\par
The tension filters can be applied not only for data regularization
but also for preconditioning in any estimation problems with smooth
models. Tomographic velocity estimation is an obvious example of such
an application \cite[]{SEG-1998-1218}.

\subsection{Mathematical theory of splines in tension}

The traditional minimum-curvature criterion implies seeking a
two-dimensional surface $f(x,y)$ in region $D$, which corresponds to
the minimum of the Laplacian power:
\begin{equation}
  \label{eqn:l2}
  \iint\limits_{D} \left|\nabla^2 f(x,y)\right|^2\,dx\,dy\;,
\end{equation}
where $\nabla^2$ denotes the Laplacian operator: $ \nabla^2 =
\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$.
\par
Alternatively, we can seek $f(x,y)$ as the solution of the biharmonic
differential equation
\begin{equation}
  \label{eqn:bi}
  (\nabla^2)^2 f(x,y) = 0\;.
\end{equation}
\cite{fung} and \cite{GEO39-01-00390048} derive
equation~(\ref{eqn:bi}) directly from~(\ref{eqn:l2}) with the help of
the variational calculus and Gauss's theorem.
\par
Formula~(\ref{eqn:l2}) approximates the strain energy of a thin
elastic plate \cite[]{timoshenko}. Taking tension into account modifies
both the energy formula~(\ref{eqn:l2}) and the corresponding
equation~(\ref{eqn:bi}). \cite{GEO55-03-02930305} suggest the
following form of the modified equation:
\begin{equation}
  \label{eqn:bit}
  \left[(1-\lambda) (\nabla^2)^2 - \lambda (\nabla^2)\right] f(x,y) = 0\;,
\end{equation}
where the tension parameter $\lambda$ ranges from 0 to 1. The
corresponding energy functional is
\begin{equation}
  \label{eqn:l2lam}
  \iint\limits_{D} \left[(1-\lambda)\,\left|\nabla^2 f(x,y)\right|^2
\;+\;
\lambda\,\left|\nabla f(x,y)\right|^2\right]\,dx\,dy\;.
\end{equation}
Zero tension leads to the biharmonic equation~(\ref{eqn:bi}) and
corresponds to the minimum curvature construction. The case of
$\lambda=1$ corresponds to infinite tension. Although infinite tension
is physically impossible, the resulting Laplace equation does have the
physical interpretation of a steady-state temperature distribution. An
important property of harmonic functions (solutions of the Laplace
equation) is that they cannot have local minima and maxima in the free
regions. With respect to interpolation, this means that, in the case
of $\lambda=1$, the interpolation surface will be constrained to have
its local extrema only at the input data locations.

Norman Sleep (2000, personal communication) points out that if the
tension term $\lambda \nabla^2$ is written in the form $\nabla \cdot
(\lambda \nabla)$, we can follow an analogy with heat flow and
electrostatics and generalize the tension parameter~$\lambda$ to a
local function depending on $x$ and $y$. In a more general form,
$\lambda$ could be a tensor allowing for an anisotropic smoothing in
some predefined directions similarly to the steering-filter method
\cite[]{SEG-1998-1851}.

To interpolate an irregular set of data values, $f_k$ at points
$(x_k,y_k)$, we need to solve equation~(\ref{eqn:bit}) under the
constraint
\begin{equation}
  \label{fk}
  f(x_k,y_k) = f_k\;.
\end{equation}
We can accelerate the solution by recursive filter preconditioning. If
$\mathbf{A}$ is the discrete filter representation of the differential
operator in equation~(\ref{eqn:bit}) and we can find a minimum-phase
filter $\mathbf{D}$ whose autocorrelation is equal to $\mathbf{A}$, then
an appropriate preconditioning operator is a recursive inverse
filtering with the filter $\mathbf{D}$. The preconditioned formulation
of the interpolation problem takes the form of the least-squares system 
\cite[]{iee}
\begin{equation}
\mathbf{K}\, \mathbf{D}^{-1} \mathbf{p} \approx  \mathbf{f}_k\;,
\label{eqn:prec2}
\end{equation}
where $\mathbf{f}_k$ represents the vector of known data, $\mathbf{K}$ is
the operator of selecting the known data locations, and $\mathbf{p}$ is
the preconditioned variable: $\mathbf{p} = \mathbf{D\, f}$. After
obtaining an iterative solution of system~(\ref{eqn:prec2}), we
reconstruct the model $\mathbf{f}$ by inverse recursive filtering:
$\mathbf{f} = \mathbf{D}^{-1}\,\mathbf{p}$. Formulating the problem in
helical coordinates \cite[]{helix0,GEO63-05-15321541} enables both the spectral
factorization of $\mathbf{A}$ and the inverse filtering with $\mathbf{D}$.

\subsection{Finite differences and spectral factorization}

\inputdir{tension}

In the one-dimensional case, one finite-difference representation of
the squared Laplacian is as a centered 5-point filter with
coefficients $(1,-4,6,-4,1)$. On the same grid, the Laplacian operator
can be approximated to the same order of accuracy with the filter
$(1/12,-4/3,5/2,-4/3,1/12)$.  Combining the two filters in accordance
with equation~(\ref{eqn:bit}) and performing the spectral
factorization, we can obtain a 3-point minimum-phase filter suitable
for inverse filtering.  Figure~\ref{fig:otens} shows a family of
one-dimensional minimum-phase filters for different values of the
parameter $\lambda$.  Figure~\ref{fig:int} demonstrates the
interpolation results obtained with these filters on a simple
one-dimensional synthetic. As expected, a small tension value
($\lambda=0.01$) produces a smooth interpolation, but creates
artificial oscillations in the unconstrained regions around sharp
changes in the gradient. The value of $\lambda=1$ leads to linear
interpolation with no extraneous inflections but with discontinuous
derivatives. Intermediate values of $\lambda$ allow us to achieve a
compromise: a smooth surface with constrained oscillations.

\sideplot{otens}{width=3in,height=2.5in}{One-dimensional minimum-phase
  filters for different values of the tension parameter $\lambda$. The
  filters range from the second derivative for $\lambda=0$ to the first
  derivative for $\lambda=1$.}

\plot{int}{width=5.5in,height=7.33in}{Interpolating a simple
  one-dimensional synthetic with recursive filter preconditioning for
  different values of the tension parameter $\lambda$. The input data are
  shown on the top. The interpolation results range from a natural
  cubic spline interpolation for $\lambda=0$ to linear interpolation for
  $\lambda=1$.}

\inputdir{Math}

To design the corresponding filters in two dimensions, we define the
finite-difference representation of operator~(\ref{eqn:bit}) on a
5-by-5 stencil. The filter coefficients are chosen with the help of
the Taylor expansion to match the desired spectrum of the operator
around the zero spatial frequency.  The matching conditions lead to
the following set of coefficients for the squared Laplacian:
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline 
-1/60 & 2/5 & 7/30 & 2/5 & -1/60 \\
\hline
2/5 & -14/15 & -44/15 & -14/15 & 2/5 \\
\hline