paper.tex

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

\label{eqn:dcp}
\end{equation}
\par
It is natural to wonder where in real life
we would encounter a \bx{Snell wave} that we could downward continue
with equation~(\ref{eqn:dcp}).
The answer is that any wave from real life
can be regarded as a sum of waves propagating in all angles.
Thus a field data set should first be decomposed
into Snell waves of all values of $p$,
and then equation~(\ref{eqn:dcp})
can be used to downward continue each $p$,
and finally the components for each $p$ could be added.
This process akin to Fourier analysis.
We now turn to Fourier analysis as a method of downward continuation
which is the same idea but
the task of decomposing data into Snell waves
becomes the task of decomposing data into sinusoids along the $x$-axis.

\subsection{Downward continuation with Fourier transform}

\sx{Fourier downward continuation}
\par
One of the main ideas in Fourier analysis
is that an impulse function
(a delta function)
can be constructed by the superposition of sinusoids
(or complex exponentials).
In the study of time series this construction is used for the
{\em 
impulse response
}
of a filter.
In the study of functions of space,
it is used to make a physical point source
that can manufacture the downgoing waves
that initialize the reflection seismic experiment.
Likewise observed upcoming waves can be Fourier transformed over $t$ and $x$.

\par
Recall in chapter~\ref{wvs/paper:wvs}, a plane wave carrying
an arbitrary waveform, specified by
equation~(\ref{wvs/eqn:mvwv}).
Specializing the arbitrary function to be
the real part of the function  $ \exp [ - i \omega (t-t_0 ) ]$  gives
\begin{equation}
\hbox{moving cosine wave} \ \ =\ \  \cos
\left[ \  \omega \left(  {x \over v }\  \sin \, \theta \ +\ 
{z \over v }\  \cos \, \theta \ -\  t \   \right) \   \right]
\label{eqn:2iei5}
\end{equation}

Using Fourier integrals on time functions we encounter the
{\em  Fourier kernel}
$\exp(-i\omega t)$.
To use Fourier integrals on the
space-axis  $x$  the spatial angular frequency must be defined.
Since we will ultimately encounter many space axes
(three for shot, three for geophone, also the midpoint and offset),
the convention will be to use a
subscript on the letter  $k$  to denote the
axis being Fourier transformed.
So  $ k_x $  is the angular spatial frequency on
the $x$-axis and  $ \exp ( i k_x x ) $  is
its Fourier kernel.
For each axis and Fourier kernel there is the question of the 
sign before the $i$.
The sign convention used here is the one used in most physics books,
namely, the one that agrees with equation~(\ref{eqn:2iei5}).
Reasons for the choice are given in chapter~\ref{ft1/paper:ft1}.
With this convention, a wave moves in the
{\em  positive}
direction along the space axes.
Thus the Fourier kernel for  $(x , z , t)$-space
will be taken to be
\begin{equation}
\hbox{Fourier kernel} \ =\  
%\ \ =\ \ 
e^{ i \, k_x x} \  e^{ i \, k_z z} \  e^{ - \, i \omega t} 
\ = \ 
\exp [ i ( k_x x  \ +\  k_z z \ -\  \omega t ) ]
\label{eqn:2iei6}
\end{equation}
\par
Now for the whistles, bells, and trumpets.
Equating (\ref{eqn:2iei5}) to the real part of (\ref{eqn:2iei6}),
physical angles and velocity are related to Fourier components.
The Fourier kernel has the form of a plane wave.
These relations should be memorized!
\begin{equation}
\vbox{\offinterlineskip
  \def\tablerule{\noalign{\hrule}}
  \hrule
  \halign {&\vrule#&
    \strut\quad\hfil#\quad\cr
    height10pt &\multispan3 &\cr
    &\multispan3 {\hfil\rm Angles and Fourier Components}\hfil&\cr
    height10pt &\multispan3 &\cr
    \tablerule
    height10pt &\omit    &&\omit  &\cr
    &$\sin \theta =\displaystyle {\strut v\ k_x \over \omega}$
          &&$\cos \theta = \displaystyle {\strut v\ k_z \over \omega}$ &\cr
   height10pt &\omit    &&\omit  &\cr
  } 
  \hrule}
\label{eqn:sincos}
\end{equation}
\par\noindent
A point in $(\omega,k_x,k_z)$-space is a plane wave.
The one-dimensional Fourier kernel extracts frequencies.
The multi-dimensional Fourier kernel extracts (monochromatic) plane waves.

\par
Equally important is what comes next.
Insert the angle definitions into the familiar
relation  $ \sin^2 \theta + \cos^2  \theta = 1 $.
This gives a most important relationship:
\begin{equation}
k_x^2 \ +\  k_z^2 \ \ =\ \  { \omega^2   \over  v^2 }
\label{eqn:wavedisp}
\end{equation}
The importance of (\ref{eqn:wavedisp}) is that it enables
us to make the distinction between
an arbitrary function and a chaotic function
that actually is a wavefield.
Imagine any function  $u(t,x,z)$.
Fourier transform it to  $U(\omega , k_x , k_z ) $.
Look in the  $ ( \omega , k_x , k_z ) $-volume for any nonvanishing
values of  $U$.
You will have a wavefield if and only if
all nonvanishing $U$ have coordinates that satisfy (\ref{eqn:wavedisp}).
Even better,
in practice the $(t,x)$-dependence at $z=0$ is usually known,
but the $z$-dependence is not.
Then the $z$-dependence is found by assuming $U$ is a wavefield,
so the $z$-dependence is inferred from (\ref{eqn:wavedisp}).
\par
Equation~(\ref{eqn:wavedisp}) also achieves fame as the
``dispersion relation of the scalar \bx{wave equation},''
a topic developed more fully in IEI.


\par
Given any $f(t)$ and its Fourier transform $F(\omega)$
we can shift $f(t)$ by $t_0$
if we multiply $F(\omega)$ by $e^{i\omega t_0}$.
This also works on the $z$-axis.
If we were given $F(k_z)$ we could shift it from the earth surface $z=0$
down to any $z_0$ by multiplying by
$e^{ik_z z_0}$.
Nobody ever gives us $F(k_z)$,
but from measurements on the earth surface $z=0$
and double Fourier transform, we can compute $F(\omega,k_x)$.
If we assert/assume that we have measured a wavefield, then we have
$k_z^2 = \omega^2/v^2 - k_x^2$,
so knowing $F(\omega,k_x)$ means we know $F(k_z)$.
Actually, we know $F(k_z,k_x)$.
Technically, we also know $F(k_z,\omega)$, but
we are not going to use it in this book.

\par
We are almost ready to extrapolate waves from the surface into the earth
but we need to know one more thing --- which square root do
we take for $k_z$?
That choice amounts to the assumption/assertion of upcoming or
downgoing waves.
With the exploding reflector model we have no downgoing waves.
A more correct analysis has two downgoing waves to think about:
First is the spherical wave expanding about the shot.
Second arises when upcoming waves hit the surface and reflect back down.
The study of multiple reflections requires these waves.


\subsection{Linking Snell waves to Fourier transforms}
To link \bx{Snell wave}s to Fourier transforms we merge
equations~(\ref{wvs/eqn:5iei4a}) and~(\ref{wvs/eqn:5iei4b})
with equations~(\ref{eqn:sincos})
\begin{eqnarray}
{k_x \over \omega} \eq
{\partial t_0 \over \partial x} \ \ \ &=&\ \ \  { \sin \, \theta  \over v }
                   \eq p
\label{eqn:kxofvp}
\\
{k_z \over \omega} \eq
{\partial t_0 \over \partial z} \ \ \ &=&\ \ \  { \cos \, \theta  \over v }
                   \eq { \sqrt{1-p^2 v^2} \over v}
\label{eqn:kzofvp}
\end{eqnarray}

\inputdir{.}
The basic downward continuation equation for upcoming waves in Fourier space
follows from equation~(\ref{eqn:dcp}) by eliminating $p$ by using
equation~(\ref{eqn:kxofvp}).
For analysis of real seismic data
we introduce a minus sign because 
equation~(\ref{eqn:kzofvp}) refers to downgoing waves
and observed data is made from up-coming waves.
\begin{equation}
U( \omega , k_x ,z+\Delta z)
 \eq
U( \omega , k_x ,z) \ 
 \exp \left( \, -\  {i \omega \Delta z \over v} \ 
\sqrt{ 1  -\ {v^2k_x^2  \over \omega^2 }\  } \   \right)
\label{eqn:dckxw}
\end{equation}
In Fourier space we delay signals by multiplying by
$e^{i\omega \Delta t}$,
analogously, equation~(\ref{eqn:dckxw}) says
we downward continue signals into the earth
by multiplying by $e^{i k_z \Delta z}$.
Multiplication in the Fourier domain
means convolution in time which can be depicted by the engineering diagram
in Figure~\ref{fig:inout}.
\plot{inout}{width=5.0in}{
	Downward continuation of a downgoing wavefield.
	}
%\newslide
\par
Downward continuation is a product relationship
in both the $\omega$-domain and the $k_x$-domain.
Thus it is a convolution in both time and $x$.
What does the filter look like in the time and space domain?
It turns out like a cone, that is,
it is roughly an impulse function
of  $x^2+z^2 - v^2 t^2$.
More precisely, it is the Huygens secondary wave source
that was exemplified by ocean waves entering a gap through a storm barrier.
Adding up the response of multiple gaps in
the barrier would be convolution over  $x$.
\inputdir{phasemod}
\par
A nuisance of using Fourier transforms in migration and modeling
is that spaces become periodic.
This is demonstrated in Figure~\ref{fig:diag}.
Anywhere an event exits the frame at a side, top, or bottom boundary,
the event immediately emerges on the opposite side.
In practice, the unwelcome effect of periodicity
is generally ameliorated by padding zeros around the data and the model.

\plot{diag}{width=6.0in,height=2.5in}{
        A reflectivity model on the left
        and synthetic data using a Fourier method on the right.
        }
%\newslide

\section{PHASE-SHIFT MIGRATION}
\sx{migration!phase-shift}
The phase-shift method of migration
begins with a two-dimensional Fourier transform (2D-FT) of the dataset.
(See chapter \ref{ft1/paper:ft1}.)
This transformed data is downward continued
with  $\exp(ik_z z)$  and subsequently evaluated
at $t=0$ (where the reflectors explode).
Of all migration methods,
the phase-shift method
most easily incorporates depth variation in velocity.
The phase angle and obliquity function are correctly included,
automatically.
Unlike Kirchhoff methods,
with the phase-shift method there is no danger of aliasing the operator.
(Aliasing the data, however, remains a danger.)

\par
Equation~(\ref{eqn:dckxw}) referred to upcoming waves.
However in the reflection experiment,
we also need to think about downgoing waves.
With the exploding-reflector concept of a zero-offset section,
the downgoing ray goes along the same path as the upgoing ray,
so both suffer the same delay.
The most straightforward way of converting one-way propagation
to two-way propagation is to multiply time everywhere by two.
Instead, it is customary to divide velocity everywhere by two.
Thus the Fourier transformed data values,
are downward continued to a depth $\Delta z$ by multiplying by
\begin{equation}
e^{ i\,k_z \Delta z }
\eq
\exp \left(  \  - \ i \, {2\,\omega \over v }\ \sqrt{
1 \ -\  { v^2\,k_x^2   \over 4\,\omega^2 } }
\  \Delta z \right)
\label{eqn:3.4}
\end{equation}
Ordinarily the time-sample interval $\Delta \tau$  for the output-migrated
section is chosen equal to the time-sample rate
of the input data (often 4 milliseconds).
Thus, choosing the depth  $ \Delta z = (v/2) \Delta \tau $,
the downward-extrapolation operator for a single time step $\Delta \tau$ is
\begin{equation}
C\eq
\exp  \left( \  - \,i \, \omega \  \Delta \tau\,\sqrt{
1 \ -\  { v^2\,k_x^2   \over 4\,\omega^2 } }
\ \ \right)
\label{eqn:3.5}
\end{equation}
Data will be multiplied many times by  $C$,  thereby downward
continuing it by many steps of  $\Delta \tau$.

\subsection{Pseudocode to working code}
Next is the task of imaging.
Conceptually, at each depth an inverse Fourier transform is followed
by selection of its value at  $t = 0$.
(Reflectors explode at  $t=0$).
Since only the Fourier transform at one point,
$t = 0$,
is needed,
other times need not be be computed.
We know the $\omega =0$ Fourier component
is found by the sum over all time,
analogously, the $t=0$ component
is found as the sum over all $\omega$.
(This is easily shown by substituting $t=0$
into the inverse Fourier integral.)
Finally,
inverse Fourier transform  $ k_x $  to  $x$.
The migration process,
computing the image from the upcoming wave  $u$,
may be summarized in the following pseudo code:
\par\noindent
\def\eq{\quad =\quad}
$$\vbox{
\begin{tabbing}
  $U(\omega, k_x, \tau = 0) = FT[u(t, x)]$                   \\
  For \= $\tau = \Delta\tau, 2\Delta\tau$,
         \ldots, end of time axis on seismogram   \\
      \> For all \= $k_x$                        \\
      \>         \> For all \= $\omega$            \\
      \>         \>        \> $C =\ $exp$(-i\omega \Delta
            \tau \sqrt{1 - v^2 k_x^2/ 4\omega^2})$   \\
      \>         \>        \> $U(\omega, k_x, \tau) = 
                               U(\omega, k_x, \tau - \Delta\tau ) \ast C$  \\
      \> For all \= $k_x$                         \\
      \>         \> Image$ (k_x, \tau) = 0$.        \\
      \>         \> For all \= $\omega$            \\
      \>         \>        \> Image$(k_x, \tau) = \ 
                              $Image$(k_x, \tau) + U(\omega, k_x, \tau)$   \\
      \> image$(x, \tau) = FT[$Image$(k_x, \tau)]$    \\
\end{tabbing}
}$$
\par\noindent
This pseudo code Fourier transforms a wavefield
observed at the earth's surface $\tau =0$,
and then it marches that wavefield down into the earth ($\tau>0$)
filling up a three-dimensional function, $U(\omega, k_x, \tau)$.
Then it selects $t=0$, the time of the exploding reflectors
by summing over all frequencies $\omega$.
(Mathematically,
this is like
finding the signal at $\omega =0$ by summing over all $t$).

\par
Turning from pseudocode to real code,
an important practical reality
is that computer memories are not big enough
for the three-dimensional function $U(\omega, k_x, \tau)$.
But it is easy to intertwine the downward continuation
with the summation over $\omega$
so a three-dimensional function need not be kept in memory.
%This is done in the real code in subroutine {\tt phasemig()}.
%\progdex{phasemig}{migration}