paper.tex

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

the waves emitted by both sides have the same polarity.
The physics of reflection coefficients says
reflections from opposite sides have opposite polarities.


\section{HYPERBOLA PROGRAMMING}
Consider an exploding reflector at the point $(z_0,x_0)$.
The location of a circular wave front at time $t$ is
$v^2t^2= (x-x_0)^2 + (z-z_0)^2$.
At the surface, $z=0$, we have the equation of the hyperbola
where and when the impulse arrives on the surface data plane $(t,x)$.
We can make a ``synthetic data plane'' by copying the explosive
source amplitude to the hyperbolic locations in the $(t,x)$ data plane.
(We postpone including the amplitude reduction caused
by the spherical expansion of the wavefront.)
Forward modeling amounts to taking every point from the $(z,x)$-plane
and adding it into the appropriate hyperbolic locations in 
the $(t,x)$ data plane.  Hyperbolas get added on top of hyperbolas.
\plot{yaxun}{width=5in}{
  Point response model to data and converse.
}
\par
Now let us think backwards.
Suppose we survey all day long and record no echos except for
one echo at time $t_0$ that we can record only at location $x_0$.
Our data plane is thus filled with zero values except the one
nonzero value at $(t_0,x_0)$.
What earth model could possibly produce such data?
\par
An earth model that is a spherical mirror with bottom at
$(z_0,x_0)$ will produce a reflection at only one point
in data space.  Only when the source is at the center
of the circle will all the reflected waves return to the source.
For any other source location, the reflected waves will
not return to the source.
The situation is summarized in Figure~\ref{fig:yaxun}.

\inputdir{beach}
\par
Above explains how an impulse at a point in image space
can transform to a hyperbola in data space,
likewise, on return, an impulse in data space can transform
to a semicircle in image space.
We can simulate a straight line in either space
by superposing points along a line.
Figure~\ref{fig:dip} shows
how points making up a line reflector diffract to a line reflection, and
how points making up a line reflection migrate to a line reflector.
\plot{dip}{width=6in}{
  Left is a superposition of many hyperbolas.
  The top of each hyperbola lies along a straight line.
  That line is like a reflector, but instead of using a continuous line,
  it is a sequence of points.
  Constructive interference gives an apparent reflection off to the side.
  Right shows a superposition of semicircles.
  The bottom of each semicircle lies along a line that
  could be the line of an observed plane wave.
  Instead the plane wave is broken into point arrivals,
  each being interpreted as coming from a semicircular mirror.
  Adding the mirrors yields a more steeply dipping reflector.
}

First we will look at the simplest,
most tutorial migration subroutine I could devise.
Then we will write an improved version
and look at some results.

\subsection{Tutorial Kirchhoff code}

Subroutine {\tt kirchslow()} below is the best tutorial
\bx{Kirchhoff migration}-modeling program I could devise.
A nice feature of this program is that it works OK
while the edge complications do not clutter it.
The program copies information from data space {\tt data(it,iy)}
to model space {\tt modl(iz,ix)} or vice versa.
Notice that of these four axes,
three are independent (stated by loops)
and the fourth is derived by the circle-hyperbola
relation $t^2=\tau^2+x^2/v^2$.
Subroutine {\tt kirchslow()}
for {\tt adj=0} copies
information from model space to data space,
i.e.~from the hyperbola top to its flanks.
For {\tt adj=1}, data summed over the hyperbola flanks
is put at the hyperbola top.
\opdex{kirchslow}{hyperbola sum}{44}{59}{user/gee}
Notice how this program has the ability to create a hyperbola
given an input impulse in $(x,z)$-space, and a circle
given an input impulse in $(x,t)$-space.

\par
The three loops in subroutine {\tt kirchslow()}
may be interchanged at will without changing the result.
To emphasize this flexibility, the loops are set at the same
indentation level.
We tend to think of fixed values of the outer two loops
and then describe what happens on the inner loop.
For example, if the outer two loops are those of the model space
{\tt modl(iz,ix)}, then
for {\tt adj=1} the program sums data along the hyperbola into
the ``fixed'' point of model space.
When loops are reordered,
we think differently
and opportunities arise for speed improvements.

\subsection{Fast Kirchhoff code}
\inputdir{sep73}
\par
Subroutine {\tt kirchslow()} can easily be speeded by a factor
that is commonly more than 30.
The philosopy of this book is to avoid minor optimizations,
but a factor of 30 really is significant,
and the analysis required for the speed up is also interesting.
Much of the inefficiency of {\tt kirchslow()} arises when
$x_{\max} \gg v t_{\max}$ because then many values of $t$
are computed beyond $t_{\max}$.
To avoid this,
we notice that for fixed offset ({\tt ix-iy}) and variable depth {\tt iz},
as depth increases,
time {\tt it} eventually goes beyond the bottom of the mesh
and, as soon as this happens,
it will continue to happen for all larger values of {\tt iz}.
Thus we can {\tt break} out of the {\tt iz} loop the first time
we go off the mesh to avoid computing anything beyond
as shown in subroutine {\tt kirchfast()}.
(Some quality compromises, limiting the aperture or the dip,
also yield speedup, but we avoid those.)
Another big speedup arises from reusing square roots.
Since the square root depends only on offset and depth,
once computed it can be used for all {\tt ix}.
Finally, these changes of variables have left us
with more complicated side boundaries,
but once we work these out,
the inner loops can be devoid of tests and
in {\tt kirchfast()}
they are in a form that is highly optimizable by many compilers.%
\opdex{kirchfast}{hyperbola sum}{45}{61}{user/gee}

\par
Originally the two Kirchhoff programs produced identical output,
but finally I could not resist adding an important feature
to the fast program,
scale factors $z/t=\cos\theta$
and $1/\sqrt{t}$ that are described elsewhere.
The fast program
allows for velocity variation with depth.
%(The {\tt break} should be conditioned
%on the maximum velocity followed by the {\tt if()} conditioned
%on the velocity at the current depth.)
When velocity varies laterally the story becomes much more complicated.

\par
Figure~\ref{fig:kfgood} shows an example.
The model includes dipping beds, syncline, anticline, fault,
unconformity, and buried focus.
The result is as expected with a ``bow tie'' at the buried focus.
On a video screen, I can see hyperbolic events originating
from the unconformity and the fault.
At the right edge are a few faint edge artifacts.
We could have reduced or eliminated these edge artifacts
if we had extended the model to the sides with some empty space.
\plot{kfgood}{width=6.00in,height=2.25in}{
  Left is the model.
  Right is diffraction to synthetic data.
}

\subsection{Kirchhoff artifacts}
\sx{artifacts}
Reconstructing the earth model with the adjoint option in \texttt{kirchfast()} \vpageref{/prog:kirchfast}
yields the result in Figure~\ref{fig:skmig}. %
\plot{skmig}{width=6.00in,height=2.25in}{
  Left is the original model.
  Right is the reconstruction.
}%
The reconstruction generally succeeds
but is imperfect in a number of interesting ways.
Near the bottom and right side, the reconstruction fades away,
especially where the dips are steeper.
Bottom fading results because in modeling the data
we abandoned arrivals after a certain maximum time.
Thus energy needed to reconstruct dipping beds near the bottom
was abandoned.
Likewise along the side we abandoned rays shooting off the frame.

\par
Difficult migrations are well known for producing semicircular reflectors.
Here we have controlled everything fairly well so none are obvious,
but on a video screen I see some semicircles.

\par
Next is the problem of the spectrum.
Notice in Figure~\ref{fig:skmig} that the reconstruction
lacks the sharp crispness of the original.
It is shown in chapter~\ref{ft1/paper:ft1}
that the spectrum of our reconstruction
loses high frequencies by a scale of $1/ | \omega |$.
Philosophically, we can think of the hyperbola summation
as integration, and integration boosts low frequencies.
Figure~\ref{fig:kirspec} shows the average over $x$
of the relevant spectra. %
\sideplot{kirspec}{width=3.00in}{
  Top is the spectrum of the the model, i.e.~the left side of
  Figure~\protect\ref{fig:skmig}.
  Bottom is the spectrum of the the reconstruction,
  i.e.~the right side of
  Figure~\protect\ref{fig:skmig}.
  Middle is the reconstruction times frequency $f$.
} %
First, notice the high frequencies are weak because
there is little high frequency energy in the original model.
Then notice that our cavalier approach to interpolation
created more high frequency energy.
Finally, notice that multiplying the spectrum of our
migrated model by frequency, $f$, brought the important
part of the spectral bands into agreement.
This suggests applying an $|\omega |$ filter to our reconstruction,
or $\sqrt{-i\omega}$ operator to both the modeling and the reconstruction,
an idea implemented in subroutine \texttt{halfdifa()} \vpageref{/prog:halfdifa}.

\par
Neither of these Kirchhoff codes addresses the issue of spatial \bx{alias}ing.
Spatial aliasing is a vexing issue of numerical analysis.
The Kirchhoff codes shown here do not work as expected
unless the space mesh size is suitably more refined than the time mesh.
Figure~\ref{fig:skmod} shows an example of forward modeling
with an $x$ mesh of 50 and 100 points. %
\plot{skmod}{width=6.00in}{
  Left is model.
  Right is synthetic data from the model.
  Top has 50 points on the $x$-axis,
  bottom has 100.
} %
(Previous figures used 200 points on space.
All use 200 mesh points on the time.)
Subroutine \texttt{kirchfast()} \vpageref{/prog:kirchfast} does interpolation by moving values
to the nearest neighbor of the theoretical location.
Had we taken the trouble to interpolate the two nearest points,
our results would have been a little better,
but the basic problem (resolved in chapter~\ref{trimo/paper:trimo}) would remain.

\subsection{Sampling and aliasing}
\inputdir{alias}
{\em 
Spatial \bx{alias}ing
}
means insufficient sampling of the data along the space axis.
This difficulty is so universal, that all migration methods
must consider it.
\par
Data should be sampled at more than two points per wavelength.
Otherwise the wave arrival direction becomes ambiguous.
Figure~\ref{fig:alias} shows synthetic data that is
sampled with insufficient density along the $x$-axis.
\sideplot{alias}{width=3.0in}{
  Insufficient spatial sampling of synthetic data.
  To better perceive the ambiguity of arrival angle,
  view the figures at a grazing angle from the side.
}
You can see that the problem becomes more acute at high frequencies 
and steep dips.
\par
There is no generally-accepted, automatic method for migrating
spatially aliased data.
In such cases, human beings may do better than machines,
because of their skill in recognizing true slopes.
When the data is adequately sampled however, computer migrations
give better results than manual methods.

\subsection{Kirchhoff migration of field data}
\inputdir{wgkirch}
Figure~\ref{fig:wgkirch} shows migrated field data.

\par
The on-line movie behind the figure shows the migration
before and after amplitude gain with time.
You can get a bad result if you gain up the data,
say with automatic gain or with $t^2$,
for display before doing the migration.
What happens is that the hyperbola flanks are then included
incorrectly with too much strength.
\par
The proper approach is to gain it first with $\sqrt{t}$
which converts it from 3-D wavefields to 2-D.
Then migrate it with a 2-D migration like {\tt kirchfast()},
and finally gain it further for display
(because deep reflectors are usually weaker).

\plot{wgkirch}{width=6.20in,height=8.5in}{
  Kirchhoff migration of Figure~\protect\ref{vela/fig:agcstack}.
  % Press button for movie comparing stack to migrated stack.
}

%\activeplot{wgkirch}{width=6.00in,height=3.6in}{ER}{
%       Kirchhoff migration of Figure~\protect\CHAPFIG{vela}{agcstack}.
%       Press button for movie comparing stack to migrated stack.
%       }

\todo{  \section{STRATIFIED V(z) FANTASY} 
        {\em Topics to cover}
        \begin{itemize}
        \item Put in the migrated WG data set!
        \item Normal raypaths from a dipping reflector in a stratified medium
        \item Which RMS velocity should you use (vertical or along normal ray)?
        \item Cheating: using Vrms(x) even though
                we initially assumed a stratified earth.
        \end{itemize}
        \section{MIGRATION VELOCITY ANALYSIS FANTASY}
        {\em Topics to cover:}
        \begin{itemize}
        \item When is it valid to use Vnmo
                (Vrms derived from stacking-velocity
                analysis) as a migration velocity?
        \item Industry practice: start with Vnmo,
                smooth it, and migrate with it.
                Then migrate again with $V_{nmo} \pm 10\%$
        \item Empirical definitions of ``overmigration'' and ``undermigration''.
        \end{itemize}
 }