paper.tex

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

\begin{equation}
{\bf \tilde d \eq N' m } \eq
 \left[ 
  \begin{array}{c}
   \tilde d_1 \\ 
   \tilde d_2 \\
   \tilde d_3 \\
   \tilde d_4 \\
   \tilde d_5 \\
   \tilde d_6 \\
   \tilde d_7 \\
   \tilde d_8 \\
   \tilde d_9 \\
   \tilde d_{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}
   m_1 \\ 
   m_2 \\
   m_3 \\
   m_4 \\
   m_5 \\
   m_6 \\
   m_7 \\
   m_8 \\
   m_9 \\
   m_{10}
  \end{array}
 \right] 
\label{eqn:NMOTarray}
\end{equation}

A program for \bx{nearest-neighbor normal moveout} as defined by
\sx{normal moveout, nearest neighbor}
equations~(\ref{eqn:NMOarray}) and (\ref{eqn:NMOTarray})
is {\tt nmo0()}.
Because of the limited alphabet of programming languages,
I used the keystroke {\tt z} to denote $\tau$.%
\opdex{nmo0}{normal moveout}{52}{61}{user/gee}
A program is a ``pull'' program if the loop creating the output
covers each location in the output and gathers the input from wherever
it may be.
A program is a ``push'' program if it takes each input and
pushes it to wherever it belongs.
Thus this NMO program is a ``pull'' program
for doing the model building (data processing),
and it is a ``push'' program for the data building.
You could write a program that worked the other way around,
namely, a loop over $t$ with $z$ found
by calculation $z=\sqrt{t^2/v^2-x^2}$.
What is annoying is that if you want a push program going
both ways, those two ways cannot be adjoint to one another.

\par
Normal moveout is a linear operation.
This means that data can be decomposed into any two parts,
early and late,
high frequency and low,
smooth and rough,
steep and shallow dip, etc.;
and whether the two parts
are NMO'ed either separately or together, the result is the same.
The reason normal moveout is a linear operation
is that we have shown it is effectively a matrix multiply operation
and that operation fulfills ${\bf N(d_1+d_2) = Nd_1+Nd_2}$.

\section{COMMON-MIDPOINT STACKING}
\inputdir{stack}
Typically, many receivers record every shot, and there are many shots
over the reflectors of interest.
It is common practice to define 
the midpoint $y =(x_s+x_g)/2$ and then to sort the
seismic traces into ``\bx{common-midpoint} gathers''.
%\todo {
%	The reason for this is shown in Figure~\FIG{midpoint},
%	which shows that the traces in a common-midpoint
%	gather all correspond to reflection from same place on the reflector.
%	This would not be true on a common-shot gather,
%	nor would it be true if the earth model were not stratified
%	(as we shall see in a later chapter).
%	\activesideplot{midpoint}{width=2.2in}{}{
%		Raypaths for a common-midpoint gather.  Various shot/receiver
%		combinations that have the same shot-receiver midpoint
%		will correspond to reflection from the same point on
%		the reflector, regardless of the offset $2h$.
%		NOTE:  Need to create this figure (JLB)!!
%		}
%	}
After sorting,
each trace on a common-midpoint gather can be transformed by NMO
into an equivalent zero-offset trace
and the traces in the gather can all be added together.
This is often called ``common-depth-point (\bx{CDP}) stacking'' 
or, more correctly, ``\bx{common-midpoint stack}ing''.

\par
The adjoint to this operation is to begin from
a model that is identical to the zero-offset trace
and spray this trace to all offsets.
There is no ``official'' definition
of which operator of an operator pair is the operator itself
and which is the adjoint.
On the one hand, I like to think of the modeling operation itself
as {\em  the} operator.
On the other hand,
the industry machinery keeps churning away at many processes
that have well-known names,
so people often think of one of them as {\em  the} operator.
Industrial data-processing operators are typically
\bx{adjoint}s
to modeling operators.

\par
Figure~\ref{fig:stack} illustrates the operator pair, consisting
of spraying out a zero-offset trace (the model) to all offsets
and the adjoint of the spraying,
which is \bx{stack}ing.
The moveout and stack operations are in subroutine {\tt stack0()}.
\opdex{stack0}{NMO stack}{53}{56}{user/gee}%
Let $\bold S'$ denote NMO, and let the stack be
defined by invoking {\tt stack0()} with the {\tt adj=1} argument.
Then $\bold S$ is the modeling operation
defined by invoking {\tt stack0()} with the {\tt adj=0} argument.
Figure~\ref{fig:stack} illustrates both.%
\sideplot{stack}{width=3in}{
  Top is a model trace $\bold m$.
  Center shows the spraying to synthetic traces, $\bold S' \bold m$.
  Bottom is the stack of the synthetic data, ${\bf S S' m}$.
}
Notice the roughness on the waveforms caused by
different numbers of points landing in one place.
Notice also the increase of \bx{AVO}
(\bx{amplitude} versus \bx{offset})
as the waveform gets compressed into a smaller space.
Finally, notice that the stack is a little rough,
but the energy is all in the desired time window.

\par
We notice a contradiction of aspirations.
On the one hand,
an operator has smooth outputs if it
``loops over output space''
and finds its input where ever it may.
On the other hand,
it is nice to have modeling and processing
be exact adjoints of each other.
Unfortunately,
we cannot have both.
If you loop over the output space of an operator,
then the adjoint operator has a loop over input space
and a consequent roughness of its output.

\subsection{Crossing traveltime curves}
\inputdir{strat}
\sx{crossing traveltime curves}
\sx{traveltime curves, crossing}
\par
Since velocity increases with depth,
at wide enough offset a deep enough path
will arrive sooner than a shallow path.
In other words,
traveltime curves for shallow events
must cut across the curves of deeper events.
Where traveltime curves cross,
NMO is no longer a one-to-one transformation.
To see what happens to the \bx{stack}ing process
I prepared Figures~\ref{fig:nmoalfa0}-\ref{fig:nmoalfa0.5}
using a typical marine recording geometry
(although for clarity I used
larger $(\Delta t,\Delta x)$) and we will use
a typical Texas gulf coast average velocity,
$v(z)=1.5+\alpha z$ where $\alpha=.5$.

First we repeat the calculation of Figure~\ref{fig:stack}
with constant velocity $\alpha=0$ and more reflectors.
We see in Figure~\ref{fig:nmoalfa0} that the \bx{stack} reconstructs
the model except for two details:
(1) the \bx{amplitude} diminishes with time, and 
(2) the early waveforms have become rounded.

\sideplot{nmoalfa0}{width=3in}{
  Synthetic CMP gather for constant velocity earth and reconstruction.
}

Then we repeat the calculation
with the Gulf coast typical velocity gradient $\alpha=1/2$.
The polarity reversal on the first arrival of the wide offset trace
in Figure~\ref{fig:nmoalfa0.5}
is evidence that in practice traveltime curves do cross.
(As was plainly evident in Figures
\ref{wvs/fig:wzl.34},
\ref{wvs/fig:wzl.20} and
\ref{wvs/fig:wzl.32}
crossing traveltime curves are even more significant elsewhere in the world.)
Comparing Figure~\ref{fig:nmoalfa0} to Figure~\ref{fig:nmoalfa1}
we see that an effect of the velocity gradient
is to degrade the \bx{stack}'s reconstruction of the model.
Velocity gradient has ruined the waveform on the shallowest event,
at about 400ms.
If the plot were made on a finer mesh
with higher frequencies,
we could expect ruined waveforms a little deeper too.

\sideplot{nmoalfa1}{width=3in}{
  Synthetic CMP gather 
  for velocity linearly increasing with depth
  (typical of Gulf of Mexico)
  and reconstruction.
}

\begin{comment}
\par
Our NMO and stack subroutines can be used for modeling or for data processing.
In designing these programs we gave no thought to signal \bx{amplitude}s
(although results showed an interesting \bx{AVO} effect in Figure~\ref{fig:stack}.)
We could redesign the programs so that the modeling operator has the
most realistic \bx{amplitude} that we can devise.
Alternately, we could design the \bx{amplitude}s
to get the best approximation to
$\bold S'\bold S\approx \bold I$
which should result in ``{\tt Stack}''
being a good approximation to ``{\tt Model}.''
I experimented with various weighting functions until
I came up with
subroutines \texttt{nmo1()} \vpageref{/prog:nmo1}
and {\tt stack1()} (like \texttt{stack0()} \vpageref{/prog:stack0})
%\sx{subroutine!{\tt stack1}, weighted NMO stack}
which embodies the \bx{weighting function} $(\tau/t)(1/\sqrt{t})$
and which produces the result in Figure~\ref{fig:nmo1alfa.5}.

\activesideplot{nmo1alfa.5}{width=3in}{ER}{
	Synthetic CMP gather 
	for velocity linearly increasing with depth
	and reconstruction with weighting functions
	in subroutine {\tt nmo1()}.
	Lots of adjustable parameters here.
	}
The result in Figure~\ref{fig:nmo1alfa.5} is very pleasing.
Not only is the \bx{amplitude} as a function of time better preserved,
more importantly, the shallow wavelets are less smeared
and have recovered their rectangular shape.
The reason the reconstruction is much better is
the cosine weighting implicit in $\tau/t$.
It has muted away much of the energy in the shallow asymptote.
I think this energy near the asymptote is harmful
because the waveform stretch is so large there.
Perhaps a similar good result could be found by experimenting
with muting programs such as \texttt{mutter()} \vpageref{/prog:mutter}.
However, subroutine \texttt{nmo1()} \vpageref{/prog:nmo1} differs from {\tt mutter()}
in two significant respects:
(1) {\tt nmo1()} is based on a theoretical concept
whereas {\tt mutter()} requires observational parameters
and (2) {\tt mutter()} applies a weighting in
the coordinates of the $(t,x)$ input space,
while {\tt nmo1()} does that
but also includes the coordinate $\tau$ of the the output space.
With {\tt nmo1()} events from different $\tau$ depths see different mutes
which is good where a shallow event asymptote
crosses a deeper event far from its own asymptote.
In practice the problem of crossing traveltime curves is severe,
as evidenced by
Figures \ref{wvs/fig:wzl.34}-\ref{wvs/fig:wzl.32}
and both weighting during NMO and muting should be used.

\progdex{nmo1}{weighted NMO}

It is important to realize that the most accurate possible
physical \bx{amplitude}s are not necessarily those for which
$\bold S'\bold S\approx \bold I$.
Physically accurate amplitudes involve many theoretical
issues not covered here.
It is easy to include some effects
(spherical divergence based on velocity depth variation)
and harder to include others
(surface ghosts and arrays).
We omit detailed modeling here
because it is the topic of so many other studies.
\end{comment}

\subsection{Ideal weighting functions for stacking}
The difference between \bx{stack}ing as defined by
\texttt{nmo0()} \vpageref{/prog:nmo0} and by
\texttt{nmo1()} \vpageref{/prog:nmo1}
is in the weighting function $(\tau/t)(1/\sqrt{t})$.
This weight made a big difference in the resolution of the stacks
but I cannot explain
whether this weighting function is the best possible one,
or what systematic procedure
leads to the best weighting function in general.
To understand this better, 
notice that $(\tau/t)(1/\sqrt{t})$
can be factored into two weights, $\tau$ and $t^{-3/2}$.
One weight could be applied before NMO and the other after.
That would also be more efficient than weighting inside NMO,
as does {\tt nmo1()}.
Additionally, it is likely that
these weighting functions should take into account
data truncation at the cable's end.
Stacking is the most important operator in seismology.
Perhaps some objective measure of quality can be defined
and arbitrary powers of $t$, $x$, and $\tau$
can be adjusted until the optimum stack is defined.
Likewise, we should consider weighting functions in the spectral domain.
As the weights $\tau$ and $t^{-3/2}$
tend to cancel one another,
perhaps we should filter with opposing filters
before and after \bx{NMO} and stack.

\subsection{Gulf of Mexico stack and AGC}
\inputdir{vscan}
\par
Next we create a ``CDP stack'' of our the Gulf of Mexico data set.
Recall the moved out common-midpoint (CMP) gather
Figure~\ref{wvs/fig:nmogath}.
At each midpoint there is one of these CMP gathers.
Each gather is summed over its offset axis.
Figure~\ref{fig:wgstack} shows the result of stacking over offset,
at each midpoint.
The result is an image of a cross section of the earth.
\plot{wgstack}{width=6.00in,height=3.6in}{
  Stack done with a given velocity profile for all midpoints.
}

In Figure~\ref{fig:wgstack} the early signals are too weak to see.
This results from the small number of traces
at early times because of the mute function.
(Notice missing information at wide offset and early time
on Figure~\ref{wvs/fig:nmogath}.)
To make the stack properly, we should divide by the number of nonzero traces.
The fact that the mute function is tapered rather than cut off abruptly
complicates the decision of what is a nonzero trace.
In general we might like to apply a weighting function of offset.
How then should the stack be weighted with time to preserve
something like the proper signal strength?
A solution is to make constant synthetic data (zero frequency).
Stacking this synthetic data gives a weight that can be used
as a divisor when stacking field data.
I prepared code for such weighted stacking,
but it cluttered the NMO and stack program and required
two additional new subroutines,
so I chose to leave the clutter in the electronic book
and not to display it here.
Instead, I chose to solve the signal strength problem
by an old standby method, Automatic Gain Control (AGC).
A divisor for the data is found by smoothing the absolute
values of the data over a moving window.
To make Figure~\ref{fig:agcstack} I made the divisor by smoothing
in triangle shaped windows about a half second long.
To do this, I used subroutine \texttt{triangle()} \vpageref{/prog:triangle}.

\plot{agcstack}{width=6.00in,height=3.6in}{
  Stack of Figure~\protect\ref{fig:wgstack} after AGC.
}

%\subsection{Jon's stack weighting function fantasy}
%The code is sprinkled around in vela/Trueamp.
%If the figure looks OK,
%my cleanup will consist of:
%
%Merge and include nmow() and stack2(), renaming them suitably.
%
%Include waiter().
%
%\activeplot{stackw}{width=6.00in,height=3.6in}{}{
%	Stack with weighting function.
%	}

\section{VELOCITY SPECTRA}
\sx{velocity spectrum}