paper.tex

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

\begin{equation}
  \label{eqn:system}
  \left(\mathbf{B} + \mathbf{D}\right) \mathbf{m} = \mathbf{d}\;,
\end{equation}
where $\mathbf{m}$ is the vector of extrapolated wavefield, and
$\mathbf{d}$ is an appropriate righthand side. The helix transform
provides us with the operator $\mathbf{B}^{-1}$, which we can use to
precondition system (\ref{eqn:system}). Introducing the change of
variables
\begin{equation}
  \label{eqn:prec}
  \mathbf{m} = \mathbf{B}^{-1} \mathbf{x}\;,
\end{equation}
we can transform the original system (\ref{eqn:system}) to the form
\begin{equation}
  \label{eqn:system2}
  \mathbf{d} = \left(\mathbf{B} + \mathbf{D}\right) \mathbf{B}^{-1} \mathbf{x} = 
  \mathbf{x} + \mathbf{D}\,\mathbf{B}^{-1} \mathbf{x}\;.
\end{equation}
When the velocity perturbation is small\footnote{In the
  linear-algebraic sense, this means that the spectral radius of
  operator $\mathbf{D}\,\mathbf{B}^{-1}$ is strictly less than one.}, even
the simple iteration
\begin{eqnarray}
  \label{ralg}
  \mathbf{x}_0 & = & \mathbf{d}\;; \\
  \mathbf{x}_k & = & \mathbf{d} - \mathbf{D}\,\mathbf{B}^{-1} \mathbf{x}_{k+1}\;
\end{eqnarray}
will converge rapidly to the desired solution. This interesting
possibility needs thorough testing.
\par
The third untested possibility (Papanicolaou, personal communication)
is to implement a clever patching in the velocity domain, applying a
constant-velocity filter locally inside each patch. Recently developed
fast wavelet transform techniques \cite[]{wavelet}, in particular the
\emph{local cosine transform}, provide a formal framework for that
approach.

\section{Conclusions}

The feasibility of multidimensional deconvolution, proven by the
helix transform, allows us to revisit the problem of implicit
wavefield extrapolation in three dimensions. The attraction of
implicit finite-difference methods lies in their unconditional
stability, a property invaluable for practical applications.
\par
We have shown that at least in the constant coefficient case (that is,
laterally invariant velocity), it is possible to implement an
extremely efficient implicit extrapolation by a recursive inverse
filtering in the helix-transformed computational model. Unfortunately,
the case of lateral velocity variations still presents a difficult
problem that may not have an exact solution. We are currently
exploring different roads to that goal.

\section{Acknowledgments}

We thank Biondo Biondi for useful discussions and for sharing his
three-dimensional expertise.

\bibliographystyle{seg}
\bibliography{SEP2,SEG,paper}

\append{The 1/6-th trick}
\inputdir{Math}
Given the filter $D_2 (k)$, defined in formula (\ref{eqn:d2k}), we can
construct an accurate approximation to the second derivative operator
$-k^2$ by considering a filter ratio (another Pad\'{e}-type
approximation) of the form
\begin{equation}
  \label{eqn:ratio}
  -k^2 \approx \frac{D_2(k)}{1 + \beta D_2 (k)}\;,
\end{equation}
where $\beta$ is an adjustable constant \cite[]{Claerbout.blackwell.85}.
The actual Pad\'{e} coefficient is $\beta=1/12$.  As pointed out by
Francis Muir, the value of $\beta = 1/4 - 1/\pi^2 \approx 1/6.726$
gives an exact fit at the Nyquist frequency $k = \pi$. Fitting the
derivative operator in the $L_1$ norm yields the value of $\beta
\approx 1/8.13$. All these approximations are shown in Figure
\ref{fig:sixth}.

\sideplot{sixth}{width=2.5in}{The second-derivative operator in the
  wavenumber domain and its approximations.}

\append{Constructing an ``isotropic'' Laplacian operator}
\inputdir{Math}
The problem of approximating the Laplacian operator in two dimensions
not only inherits the inaccuracies of the one-dimensional
finite-difference approximations, but also raises the issue of
azimuthal asymmetry. For example, the usual five-point filter
\begin{equation}
\label{eqn:usual}
F_5 =
\begin{array}{|r|r|r|}
\hline
0 & 1 & 0 \\ \hline
1 & -4 & 1 \\ \hline
0 & 1 & 0  \\ \hline 
\end{array}
\end{equation}
exhibits a clear difference between the grid directions and the
directions at a 45-degree angle to the grid. To overcome this
unpleasant anisotropy, we can consider a slightly larger filter of the
form
\begin{equation}
\label{eqn:new}
F_9 = \begin{array}{|r|r|r|}
\hline
\alpha & \gamma & \alpha \\ \hline
\gamma & -4 (\alpha+\gamma) & \gamma \\ \hline
\alpha & \gamma & \alpha  \\ \hline 
\end{array}
\end{equation}
where the constants $\alpha$ and $\gamma$ are to be defined. The
Fourier-domain representation of filter (\ref{eqn:new}) is
\begin{equation}
  \label{eqn:l9k}
  F_9 (k_x,k_y) = 4\,\alpha\,[\cos{k_x}\,\cos{k_y} - 1] +
  2\,\gamma\,[\cos{k_x}+\cos{k_y}-2]\;,
\end{equation}
and the isotropic filter that we can try to approximate is defined
analogously to its one-dimensional equivalent, as follows:
\begin{equation}
  \label{eqn:lk}
  F (k_x,k_y) = 2 (\cos{k} -1) = 2 (\cos{\sqrt{k_x^2+k_y^2}} - 1)\;.
\end{equation}
Comparing equations (\ref{eqn:l9k}) and (\ref{eqn:lk}), we notice that
they match exactly, when either of the wavenumbers $k_x$ or $k_y$ is
equal to zero, provided that 
\begin{equation}
  \label{eqn:alpha}
  \alpha = \frac{1-\gamma}{2}\;.
\end{equation}
Therefore, we can reduce the problem to estimating a single
coefficient $\gamma$. Another way of expressing this conclusion is to
represent filter $F_9$ in equation (\ref{eqn:l9k}) as a linear
combination of filter $F_5$ from equation (\ref{eqn:l9k}) and its
rotated version \cite[]{cole}, as follows:
\begin{equation}
\label{eqn:new2}
F_9 = \gamma\, 
\begin{array}{|r|r|r|}
\hline
0 & 1 & 0 \\ \hline
1 & -4  & 1 \\ \hline
0 & 1 & 0  \\ \hline 
\end{array}
+ (1-\gamma)\,
\begin{array}{|r|r|r|}
\hline
1/2 & 0 & 1/2 \\ \hline
0 & -2  & 0 \\ \hline
1/2 & 0 & 1/2  \\ \hline 
\end{array}
\end{equation}
With the value of $\gamma=0.5$, filter $F_9$ takes the value
\begin{equation}
\label{eqn:mc}
F_9 = \begin{array}{|r|r|r|}
\hline
1/4 & 1/2 & 1/4 \\ \hline
1/2 & -3 & 1/2 \\ \hline
1/4 & 1/2 & 1/4  \\ \hline 
\end{array}
\end{equation}
and corresponds precisely to the nine-point McClellan filter
\cite[]{mc,GEO56-11-17781785}. On the other hand, the value of
$\gamma=2/3$ gives the least error in the vicinity of the zero
wavenumber $k$. In this case, the filter is
\begin{equation}
\label{eqn:my}
F_9 = \begin{array}{|r|r|r|}
\hline
1/6 & 2/3 & 1/6 \\ \hline
2/3 & -10/3 & 2/3 \\ \hline
1/6 & 2/3 & 1/6  \\ \hline 
\end{array}
\end{equation}
Errors of different approximations are plotted in Figure
\ref{fig:laplace}\footnote{Another way of constructing
  circular-symmetric filters is suggested by the rotated McClellan
  transform \cite[]{Biondi.sep.77.27}.}

\plot{laplace}{width=6in}{The numerical anisotropy error of different
  Laplacian approximations. Both the five-point Laplacian (plot a) and
  its rotated version (plot b) are accurate along the axes, but
  exhibit significant anisotropy in between at large wavenumbers. The
  nine-point McClellan filter (plot c) has a reduced error, while the
  filter with $\gamma=2/3$ (plot d) has the flattest error around the
  origin.}
\par

\inputdir{laplace}
Under the helix transform, a filter of the general form
(\ref{eqn:new}) becomes equivalent to a one-dimensional filter with
the $Z$ transform
\begin{eqnarray}
\label{eqn:f9z}
F_9 (Z) & = & \alpha\,Z^{-N_x-1} + \gamma\,Z^{-N_x} + \alpha\,Z^{-N_x+1} + \gamma\,Z^{-1} 
-4\,(\alpha + \gamma) \nonumber \\ & & + 
\gamma\,Z + \alpha\,Z^{N_x-1} + \gamma\,Z^{N_x} + \alpha\,Z^{N_x+1}\;,
\end{eqnarray}
where $N_x$ is the helix period (the number of grid points in the $x$
dimension). To find the inverse of a convolution with filter
(\ref{eqn:f9z}), we factorize the filter into the causal minimum-phase
component and its adjoint:
\begin{equation}
\label{eqn:factor}
F_9 (Z) = P (Z) P (1/Z)\;.
\end{equation}
To find the coefficients of the filter $P$, any one-dimensional
spectral factorization method can be applied. It is important to point
out that the result of factorization (neglecting the numerical errors)
does not depend on $N_x$. Another approach is to define a residual
error vector for the coefficients of Z in equation (\ref{eqn:factor})
and minimize it for some particular norm. For example, minimizing the
$\mathcal{L}_1$ norm when $F_9$ is the McClellan filter
(\ref{eqn:mc}), we discover that the filter $P$, after transforming
back to two dimensions, takes the form
\begin{equation}
\label{eqn:l1long}
\begin{array}{|r|r|r|r|r|r|r|r|}
\hline
     &         &           &        &     -1.6094 & 0.4293 & 0.05157 & 0.017406 \\ \hline
0.01428 & 0.033513 & 0.0808 & 0.2543 & 0.3521  & 0.1553 &   &          \\ \hline
\end{array}
\end{equation}
The results of applying a recursive deconvolution with filter
(\ref{eqn:l1long}) are shown in Figure \ref{fig:inv-laplace}.  An
essentially similar procedure, only with a different set of filters,
works for implicit wavefield extrapolation.
  
\plot{inv-laplace}{width=4in,height=4in}{Inverting the Laplacian
  operator by a helix deconvolution.  The top left plot shows the
  input, which contains a single spike and the causal minimum-phase
  filter $P$. The top right plot is the result of inverse filtering.
  As expected, the filter is deconvolved into a spike, and the spike
  turns into a smooth one-sided impulse. After the second run, in the
  backward (adjoint) direction, we obtain a numerical solution of
  Laplace's equation! In the two bottom plots, the solution is shown
  with grayscale and contours.}

%\activeplot{name}{width=6in,height=}{}{caption}
%\activesideplot{name}{height=1.5in,width=}{}{caption}
%
%\begin{equation}
%\label{eqn:}
%\end{equation}
%
%\ref{fig:}
%(\ref{eqn:})