paper.tex

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

undersampling. The bottom plots in Figure \ref{fig:fft-inv} correspond
to increasing the grid size by a factor of three. Some of the
artifacts have been suppressed, at the expense of dealing with a
larger grid.

\plot{fft-inv}{width=6in}{The left plots show the reconstruction of the original data 
  after transforming back from the $t^2$ grid to the original $t$
  grid.  The right plots show the difference with the original model.
  Top: using the original grid size ($N_t = 200$). Bottom: increasing
  the grid size by a factor of three.}

To perform an accurate transform of the grid, I adopted the following
method, inspired by \cite[]{Claerbout.sep.48.347}. Let $d_{\mbox{\tiny
    new}}$ denote the data on the new grid and $d_{\mbox{\tiny old}}$
be the data on the old grid. If $L$ is the interpolation operator,
defined on the new grid, then the optimal least-square transformation
is
\begin{equation}
  \label{interp}
  d_{\mbox{\tiny new}} = (L^T L)^{-1}\,L\,d_{\mbox{\tiny old}}\;,
\end{equation}
where $L^T$ denotes the adjoint interpolation operator. The operator
$(L^T L)^{-1}$ provides a proper scaling of the result. If we use
simple linear interpolation for the $L$ operator, then $L^T L$ is a
tridiagonal matrix, which can be easily inverted (in $8 N$
operations). If some parts in $d_{\mbox{\tiny new}}$ are not fully
constrained, then the tridiagonal matrix is not invertible. To obtain
a solution in this case, we can include a regularization operator $D$
in (\ref{interp}), as follows:
\begin{equation}
  \label{regul}
  d_{\mbox{\tiny new}} = (L^T L + \epsilon^2 D)^{-1}\,L\,d_{\mbox{\tiny
      old}}\;,
\end{equation}
A convenient choice for $D$ is a second derivative operator,
represented with the second-order finite-difference approximation.
This operator allows the selection of the smoothest possible function
$d_{\mbox{new}}$ while preserving the efficient tridiagonal structure
of $L^T L + \epsilon^2 D$. In this problem, the parameter $\epsilon$
can be chosen as small as possible, as long as it prevents the
inversion from getting unstable.
\par
\subsection{Suppressing wraparound artifacts of the Fourier method}

\inputdir{impls}
\par
The periodic boundary conditions both in the squared time $\sigma$ and
the spatial coordinate $x$, implied by the Fourier approach, are
artificial in the problem of velocity continuation. The artificial
periodicity is convenient from the computational point of view.
However, false periodic events (wraparound artifacts) should be
suppressed in the final output. A natural method for attacking this
problem is to apply zero padding in the physical space prior to
Fourier transform. Of course, this method involves an additional
expense of the grid size increase.
\par
\plot{fft-imp}{width=6in}{Impulse responses (Green's functions) of velocity continuation, 
  computed by the Fourier method. Top: without zero padding, bottom:
  with zero padding. The left plots correspond to continuation to a
  larger velocity ($+1$ km/sec); the right plots, smaller velocity,
  ($-1$ km/sec).}
\par
The top plots in Figure \ref{fig:fft-imp} show the numerical impulse
responses of velocity continuation, computed by the Fourier method.
The initial data contained three spikes, passed through a narrow-band
filter. Theoretically, continuation to larger velocity (the left plot)
should create three elliptical wavefronts, and continuation to smaller
velocity (right plot) should create three hyperbolic wavefronts
\cite[]{GEO50-01-01100126}. We can see that the results are largely
contaminated with wraparound artifacts. The result of applying zero
padding (the bottom plots in Figure \ref{fig:fft-imp}) shows most of
the artifacts suppressed.
\par
Chebyshev spectral method, discussed in the next section, provides a
spectral accuracy while dealing correctly with non-periodic data.
\par 
\section{Chebyshev Approach}
\inputdir{sigvc}
For an alternative spectral approach, I adopted the Chebyshev-$\tau$
method \cite[]{lanc,orszag}. The Chebyshev-$\tau$ velocity continuation
algorithm consists of the following steps:
\begin{enumerate}
\item Transform the regular grid in $t$ to Gauss-Lobato collocation
  points, required for the fast Chebyshev transform.  First, a new
  variable $\xi$ is introduced by the shift transform:
  \begin{equation}
    \label{t2x}
    \xi = 1 - \frac{2\,t^2}{T^2}
  \end{equation}
  so that the domain $0 \leq t \leq T$ is mapped into the domain $1
  \geq \xi \geq -1$. Second, the $\xi$ grid points are distributed
  regularly in the cosine projection: $\xi_j = \cos(\frac{\pi j}{N}), j
  = 0,1,2,\ldots,N$.
\item Transform the initial image $P_0 (x,t)$ into the Chebyshev space
  in $\xi$ and Fourier transform in $x$, using the FFT algorithm. The
  Chebyshev-Fourier representation of $P_0 (x,t)$ is
  \begin{equation}
    \label{cheb}
    P_0 (x,t) = \sum_{k=-N_x/2}^{N_x/2-1}\sum_{j=0}^{N_t} 
    \hat{P}_{kj} T_j (\xi) e^{i k x}\;, 
  \end{equation}
 where $T_j$ denotes the Chebyshev polynomial of degree $j$.
\item Apply equation (\ref{PDE}) to advance the image in velocity $v$.
 It is convenient to rewrite this equation in the form
 \begin{equation}
   \frac{\partial P}{\partial v} =  \frac{v\,T^2}{4}\,
   \int d\xi\,\frac{\partial^2 P}{\partial x^2}\;.
   \label{chebPDE} 
 \end{equation}
  In the Chebyshev-$\tau$ domain, the double differentiation in $x$ is
  performed by multiplying the Fourier transform of $P$ by $-k^2$, and
  integration in $\xi$ is performed
  as a direct operations on the Chebyshev coefficients. In
  particular, if $\sum_{j=0}^{N} a_j T_j (\xi)$ is the Chebyshev
  representation of the function $f (\xi)$, then the coefficients
  $b_j$ of $\int f (\xi) d\xi$ are defined by the relation
  \begin{equation}
    \label{int}
    2 j\,b_j = c_{j-1} a_{j-1} - a_{j+1}
  \end{equation}
  where $c_0 = 2$, $c_j = 0$ for $j < 0$, and $c_j = 1$ for $j > 0$.
  The constant of integration (and, correspondingly, the coefficient
  $b_0$) can be found at each velocity step from the boundary
  conditions (\ref{BC}), which are transformed to the form
  \begin{equation}
    \label{chebBC}
    \left.P\right|_{\xi=-1} = \sum_{j=0}^{N_t} \hat{P}_{kj} (-1)^j = 0\;.
  \end{equation}
  
  For the velocity advancement I used an implicit Crank-Nicolson scheme,
  which is unconditionally stable independent of the velocity step size.
  By writing equation (\ref{chebPDE}) in the matrix form
  \begin{equation}
    \frac{\partial \mathbf{P}}{\partial v} = \mathbf{A}\,\mathbf{P}\;,
    \label{matrix}
  \end{equation}
  the Crank-Nicolson advancement is represented by the equation
    \begin{equation}
    \label{CN}
    \mathbf{P}_{v+dv} = \left(\mathbf{I} - \mathbf{A}\,\frac{dv}{2}\right)^{-1}
    \left(\mathbf{I} + \mathbf{A}\,\frac{dv}{2}\right) \mathbf{P}_v\;,
  \end{equation}
  where $\mathbf{I}$ is the identity matrix. The inverted matrix
  $\left(\mathbf{I} - \mathbf{A}\,\frac{dv}{2}\right)$ has a tridiagonal
  structure, except for the first row, implied by the boundary
  condition (\ref{chebBC}). A careful treatment of the boundary
  condition by the matrix-bordering method \cite[]{faddeev,boyd} allows
  for an efficient inversion at a tridiagonal solver speed.
\item Transform the result of the velocity advancement back to the
  physical domain.
\item Transform the grid back to being regularly space in $t$.
\end{enumerate}
\par
\plot{cheb1}{width=6in}{Synthetic seismic data before (left) and 
  after (right) transformation to the Chebyshev grid in squared time.}
\par
\plot{cheb1-inv}{width=6in}{The left plots show the reconstruction of the original data 
  after transforming back from the Chebyshev grid to the original $t$
  grid.  The right plots show the difference with the original model.
  Top: using the original grid size ($N_t = 200$). Bottom: increasing
  the grid size by a factor of three.}
\par
The first advantage of the Chebyshev approach comes from the better
conditioning of the grid transform.  Figure \ref{fig:cheb1} shows the
synthetic data before and after the grid transform.  Figure
\ref{fig:cheb1-inv} shows a reconstruction of the original data after
transforming back from the Chebyshev grid (Gauss-Lobato collocation
points). The difference with the original image is negligibly small.
\par
\plot{cheb1-fft}{width=6in}{Left: Synthetic data after Chebyshev
  transform.  Right: the real part of the Fourier transform in the
  space coordinate.}
\par
The second advantage is the compactness of the Chebyshev
representation.  Figure \ref{fig:cheb1-fft} shows the data after the
decomposition into Chebyshev polynomials in $\xi$ and Fourier
transform in $x$. We observe a very rapid convergence of the Chebyshev
representation: a relatively small number of polynomials suffices to
represent the data.
\inputdir{impls}
\par
\plot{cheb-impl}{width=6in}{Impulse responses (Green's functions) of
  velocity continuation, computed by the Chebyshev-$\tau$ method. Top:
  without zero padding, bottom: with zero padding on the $x$ axis. The
  left plots correspond to continuation to a larger velocity ($+1$
  km/sec); the right plots, smaller velocity, ($-1$ km/sec).}
\par
The third advantage is the proper handling of the non-periodic boundary
conditions. Figure \ref{fig:cheb-impl} shows the velocity continuation
impulse responses, computed by the Chebyshev method. As expected, no
wraparound artifacts occur on the time axis, and the accuracy of the
result is noticeably higher than in the case of finite differences
(Figure \ref{fig:fd-imp}).
\par 
\section{Conclusions}
\inputdir{sigvc}
\par
I have applied two spectral methods for a numerical solution of the
velocity continuation problem.
\par
The Fourier method is attractive because of its numerical efficiency.
However, it requires additional computational effort to suppress
numerical artifacts: the inaccuracy of the grid transform and
the artificial periodicity in the physical space.
\par
The Chebyshev-$\tau$ method is free of most of these difficulties,
although its overall efficiency can be slightly inferior to that of
the Fourier method.
\par
Both methods possess a ``spectral'' accuracy, which is highly desired
if accuracy is a concern.
\par
\plot{mig-impl}{width=6in}{Top left: synthetic model (the ideal image). 
  Top right: synthetic data. 
  Bottom left: the result of velocity continuation 
  with the Fourier method. Bottom right: the result of 
  velocity continuation with the Chebyshev method.}
\par
Figure \ref{fig:mig-impl} compares the results of velocity
continuation with different methods. The top left plot shows an
implied subsurface model (an ``ideal image''). The top right plot is
the corresponding synthetic data. The bottom left plot is the output
of the Fourier method, and the bottom right plot is the output of the
Chebyshev method. The Fourier result shows a poor quality in the
shallow part (caused by subsampling in the $t^2$ grid). The wraparound
artifacts were suppressed by a zero-padding correction. The quality of
the Chebyshev result is noticeably higher. It is close to the best
possible accuracy, under the natural limitations of seismic resolution.

\bibliographystyle{seglike} \bibliography{SEG,SEP2,spec}

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