paper.tex

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

\long\def\HideThis#1{}
%\def\SEPCLASSLIB{../../../../sepclasslib}
\def\CAKEDIR{.}
\def\RMS{{{\sc rms}}} 
\def\RMS{{\rm RMS}}

\title{Moveout, velocity, and stacking}
\author{Jon Claerbout}
\maketitle
\label{paper:vela}

%{\today. \em  This chapter is owned by JFC.}

In this chapter we handle data as though the earth had no dipping reflectors.
The earth model is one of stratified layers
with velocity a (generally increasing) function of depth.
We consider reflections from layers,
which we process by normal moveout correction (NMO).
The NMO operation is an interesting example
of many general principles of linear operators and numerical analysis.
Finally, using NMO, we estimate the earth's velocity with depth
and we stack some data,
getting a picture of an earth with dipping layers.
This irony, that techniques developed for a stratified earth
can give reasonable images of non-stratified reflectors,
is one of the ``lucky breaks'' of seismic processing.
We will explore the limitations of this phenomenon
in the chapter on dip-moveout.

\par
First, a few words about informal language.
The inverse to velocity arises more frequently
in seismology than the velocity itself.
This inverse is called the ``slowness.''
In common speech, however, the word ``velocity'' is a catch-all,
so what is called a ``velocity analysis''
might actually be a plane of slowness versus time.

\section{INTERPOLATION AS A MATRIX}
Here we see how general principles of linear operators
are exemplified by linear interpolation.
Because the subject matter is so simple and intuitive,
it is ideal to exemplify abstract mathematical concepts
that apply to all linear operators.
\par
Let an integer $k$ range along a survey line,
and let data values $x_k$ be packed into a vector $\bold x$.
(Each data point $x_k$ could also be a seismogram.)
Next we resample the data more densely,
say from 4 to 6 points.
For illustration, I follow a crude
\bx{nearest-neighbor interpolation}
\sx{interpolation, nearest-neighbor}
scheme by sprinkling ones along the diagonal of a rectangular matrix
that is
\begin{equation}
\bold y \eq  \bold B \, \bold x
\label{eqn:interp1}
\end{equation}
where
\begin{equation}
 \left[ 
  \begin{array}{c}
   y_1 \\ 
   y_2 \\
   y_3 \\
   y_4 \\
   y_5 \\
   y_6
  \end{array}
 \right] 
\eq
 \left[ 
  \begin{array}{cccc}
  1 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 1 \\
  0 & 0 & 0 & 1
  \end{array}
 \right] \;
 \left[ 
  \begin{array}{c}
   x_1 \\ 
   x_2 \\
   x_3 \\
   x_4
  \end{array}
 \right] 
\label{eqn:interp2}
\end{equation}
The interpolated data is simply
$\bold y = (x_1, x_2,x_2,x_3, x_4,x_4)$.
The matrix multiplication~(\ref{eqn:interp2})
would not be done in practice.
Instead there would be a loop running over the space of the
outputs $\bold y$ that picked up values from the input.

\subsection{Looping over input space}
The obvious way to program a deformation
is to take each point from the {\em  input} space
and find where it goes on the output space.
Naturally, many points could land in the same place,
and then only the last would be seen.
Alternately, we could first erase the output space,
then add in points, and finally divide by the number of points
that ended up in each place.
The biggest aggravation is that some places could end up with no points.
This happens where the transformation \bx{stretch}es.
There we need to decide whether to interpolate the missing points,
or simply low-pass filter the output.

\subsection{Looping over output space}
The alternate method that is usually preferable
to looping over input space
is that
our program have a loop over the space of the {\em  outputs,}
and that each output find its input.
The matrix multiply of (\ref{eqn:interp2})
can be interpreted this way.
Where the transformation \bx{shrink}s is a small problem.
In that area
many points in the input space are ignored,
where perhaps they should somehow be averaged with their neighbors.
This is not a serious problem unless we are contemplating
iterative transformations back and forth between the spaces.

\par
We will now address interesting questions
about
the reversibility of these deformation transforms.

\subsection{Formal inversion}
We have thought of equation~(\ref{eqn:interp1})
as a formula for finding $\bold y$ from $\bold x$.
Now consider the opposite problem, finding $\bold x$ from $\bold y$.
Begin by multiplying equation (\ref{eqn:interp2})
by the \bx{transpose matrix}
to define a new quantity $\tilde\bold x$:
\begin{equation}
 \left[ 
  \begin{array}{c}
   \tilde x_1 \\ 
   \tilde x_2 \\
   \tilde x_3 \\
   \tilde x_4
  \end{array}
 \right] 
\eq
 \left[ 
  \begin{array}{ccccccc}
  1 & 0 & 0 & 0 & 0 & 0\\
  0 & 1 & 1 & 0 & 0 & 0\\
  0 & 0 & 0 & 1 & 0 & 0\\
  0 & 0 & 0 & 0 & 1 & 1
  \end{array}
 \right] \;
 \left[ 
  \begin{array}{c}
   y_1 \\ 
   y_2 \\
   y_3 \\
   y_4 \\
   y_5 \\
   y_6
  \end{array}
 \right] 
\end{equation}
$\tilde \bold x$ is not the same as $\bold x$,
but these two vectors have the same dimensionality
and in many applications
it may happen that
$\tilde \bold x$ is a good approximation to $\bold x$.
In general, 
$\tilde \bold x$ may be called an ``image'' of $\bold x$.
Finding the image is the first step of finding $\bold x$ itself.
Formally, the problem is
\begin{equation}
 \bold y  \eq \bold B  \, \bold x
						\label{eqn:formalproblem}
\end{equation}
And the formal solution to the problem is
\begin{equation}
\bold x \eq ({\bf B'\, B})^{-1} \, {\bf B'} \, \bold y
							\label{eqn:interpinv}
\end{equation}
Formally, we verify this solution by substituting
(\ref{eqn:formalproblem}) into
(\ref{eqn:interpinv}).
\begin{equation}
\bold x \eq ( {\bf B' \, B} )^{-1} \, ({\bf B'} \, \bold B) \,\bold x
		   \eq \bold I \, \bold x  \eq \bold x
\end{equation}
In applications,
the possible nonexistence of an inverse for the matrix $( {\bf B' \, B} )$
is always a topic for discussion.
For now we simply examine this matrix for the interpolation problem.
We see that it is diagonal:
\begin{equation}
\bold B' \, \bold B
\eq
 \left[ 
  \begin{array}{ccccccc}
  1 & 0 & 0 & 0 & 0 & 0\\
  0 & 1 & 1 & 0 & 0 & 0\\
  0 & 0 & 0 & 1 & 0 & 0\\
  0 & 0 & 0 & 0 & 1 & 1
  \end{array}
 \right] \;
 \left[ 
  \begin{array}{cccc}
  1 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 1 \\
  0 & 0 & 0 & 1
  \end{array}
 \right] \;
 \eq
 \left[ 
  \begin{array}{cccc}
  1 & 0 & 0 & 0 \\
  0 & 2 & 0 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 2
  \end{array}
 \right] \;
\end{equation}
So, ${\bf \tilde x}_1 =   \bold x_1$; but
    ${\bf \tilde x}_2 = 2 \bold x_2$.
To recover the original data,
we need to divide ${\bf \tilde x}$ by the diagonal matrix $\bold B'\,\bold B$.
Thus, matrix inversion is easy here.

\par
Equation
(\ref{eqn:interpinv})
has an illustrious reputation, which
arises in the context of ``least squares.''
\bxbx{Least squares}{least squares}
is a general method for solving sets of equations
that have more equations than unknowns.

\par
Recovering $\bold x$ from $\bold y$ using equation
(\ref{eqn:interpinv})
presumes the existence of the inverse of
$\bold B'\,\bold B$.
As you might expect, this matrix is nonsingular when ${\bf B}$
{\em  stretches} the data,
because then a few data values are distributed
among a greater number of locations.
Where the transformation {\em  squeezes} the data,
$\bold B'\,\bold B$
must become singular,
since returning
uniquely to the uncompressed condition is impossible.
%chapter \LS\ explains how to handle singular matrices,
%enabling us to decompress losing only the high frequencies.

\par
We can now understand why an adjoint operator is often an
approximate inverse.
This equivalency happens in proportion to the nearness of the matrix
$\bold B'\,\bold B$
to an identity matrix.
The interpolation example we have just examined is one in which
$\bold B'\,\bold B$
differs from an identity matrix merely by a scaling.

\section{THE NORMAL MOVEOUT MAPPING}
\inputdir{stretch}
Recall the traveltime equation~(\ref{wvs/eqn:hyper}).
\begin{eqnarray}
v^2 \, t^2
&=&
z^2 \ +\  x^2
\\
t^2
&=&
\tau^2 \ +\ { x^2   \over  v^2 } 
\end{eqnarray}
where $\tau$ is traveltime depth.
This equation gives either time from a surface source
to a receiver at depth $\tau$,
or it gives time to a surface receiver
from an image source at depth $\tau$.

\par
A seismic \bx{trace} is a signal $d(t)$
recorded at some constant $x$.
We can convert the trace
to a ``vertical propagation'' signal
$m(\tau)=d(t)$
by stretching $t$ to $\tau$.
This process is called
``\bx{normal moveout} correction'' (NMO).
Typically we have many traces at different $x$ distances
each of which theoretically produces the same
hypothetical zero-offset trace.
Figure~\ref{fig:stretch} shows a marine shot profile
before and after NMO correction at the water velocity.
You can notice that the wave packet reflected from the ocean bottom
is approximately a constant width on the raw data.
After NMO, however,
this waveform broadens considerably---a phenomenon known
as ``NMO stretch."

\sideplot{stretch}{width=3.0in}{
  Marine data moved out with water velocity.
  Input on the left, output on the right.
%  Press button for \bx{movie} sweeping through
%  velocity (actually through slowness squared).
}

\par
The  NMO  transformation ${\bf N}$ is representable as a square matrix.
The matrix  ${\bf N}$  is a $(\tau,t)$-plane containing all zeros
except an interpolation operator centered along the hyperbola.
The dots in the matrix below are zeros.
The input signal $d_t$ is put into the vector $\bold d$.
The output vector $\bold m$---i.e.,~the NMO'ed signal---is simply
$(d_6,d_6,d_6, d_7,d_7, d_8,d_8, d_9, d_{10}, 0)$.
In real life examples such as Figure~\ref{fig:stretch}
the subscript goes up to about one thousand instead of
merely to ten.

\begin{equation}
{\bf m\eq Nd} \eq
 \left[ 
  \begin{array}{c}
   m_1 \\ 
   m_2 \\
   m_3 \\
   m_4 \\
   m_5 \\
   m_6 \\
   m_7 \\
   m_8 \\
   m_9 \\
   m_{10}
  \end{array}
 \right] 
\eq
 \left[ 
  \begin{array}{cccccccccc}
   .&.&.&.&.&1&.&.&.&. \\
   .&.&.&.&.&1&.&.&.&. \\
   .&.&.&.&.&1&.&.&.&. \\
   .&.&.&.&.&.&1&.&.&. \\
   .&.&.&.&.&.&1&.&.&. \\
   .&.&.&.&.&.&.&1&.&. \\
   .&.&.&.&.&.&.&1&.&. \\
   .&.&.&.&.&.&.&.&1&. \\
   .&.&.&.&.&.&.&.&.&1 \\
   .&.&.&.&.&.&.&.&.&.
  \end{array}
 \right] \;
 \left[ 
  \begin{array}{c}
   d_1 \\ 
   d_2 \\
   d_3 \\
   d_4 \\
   d_5 \\
   d_6 \\
   d_7 \\
   d_8 \\
   d_9 \\
   d_{10}
  \end{array}
 \right] 
 \label{eqn:NMOarray}
\end{equation}

You can think of the matrix as having a horizontal $t$-axis and
a vertical $\tau$-axis.
The 1's in the matrix are arranged on the hyperbola
$t^2=\tau^2+x_0^2/v^2$.
The transpose matrix defining some ${\bf \tilde d}$
from $\bold m$ gives synthetic data ${\bf \tilde d}$ from the zero-offset
(or stack) model $\bold m$, namely,