paper.tex

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

%An aspect of the computation that was hidden in the pseudo code
%That you can see in the real code
%is that we must also handle negative frequencies.
%When doing so, we use the signum function
%$\sgn(\omega ) =\omega /|\omega |$ to choose the sign of $k_z$
%so that $k_z$
%reverses polarity when $\omega$ reverses polarity.
%\progdex{signum}{sign function}

\par
Conjugate migration (modeling) proceeds in much the same way.
Beginning from an upcoming wave that is zero at great depth,
the wave is marched upward in steps
by multiplication with  $\exp ( i k_z \Delta z )$.
As each level in the earth is passed,
exploding reflectors from that level are added into the upcoming wave.
Pseudo code for modeling the upcoming wave  $u$  is
\par\noindent
$$\vbox{
\begin{tabbing}
Image$(k_x, z) = FT[$image$(x, z)]$    \\
For \= all $\omega$ and all $k_x$    \\
    \> $U(\omega, k_x) = 0$.         \\
For all $\omega$ \{                  \\
For all $k_x$ \{                     \\
For $z = z_{\rm max}, z_{\rm max} -\Delta z, z_{\rm max} -2\Delta z
         $, \ldots, $0$ \{           \\
    \> $C = \exp(+i\Delta z\omega\sqrt{v^{-2} -
            {k_x}^2/\omega^2})$          \\
    \> $U(\omega, k_x) = U(\omega, k_x) \ast C $     \\
    \> $U(\omega, k_x) = U(\omega, k_x)  + $Image$(k_x, z)$     \\
    \> \} \} \}                      \\
$u(t, x) = FT[U(\omega, k_x)]$       \\
\end{tabbing}
}$$
\par\noindent
%Some real code for this job is in subroutine {\tt phasemod()}.
%\progdex{phasemod}{diffraction}

\par
The positive sign in the complex exponential is a combination of two negatives,
the {\em  up}
coming wave and the {\em  up}ward extrapolation.
In principle,
the three loops on  $\omega$, $k_x$, and  $z$  are interchangeable,
however, since this tutorial program uses
a velocity $v$ that is a constant function of depth,
I speeded it by a large factor by putting the $z$-loop on the inside
and pulling the complex exponential out of the inner loop.
%Figure~\ref{fig:stormhole} was made with subroutine~\texttt{phasemod()} \vpageref{/prog:phasemod}.

\subsection{Kirchhoff versus phase-shift migration}
\inputdir{sigmoid}
In chapter~\ref{krch/paper:krch}, we were introduced to the Kirchhoff
migration and modeling method by means of subroutines
\texttt{kirchslow()} \vpageref{/prog:kirchslow} and \texttt{kirchfast()} \vpageref{/prog:kirchfast}.
From chapter~\ref{ft1/paper:ft1} we know that these routines should be
supplemented by a $\sqrt{-i\omega}$
filter such as subroutine \texttt{halfdifa()} \vpageref{/prog:halfdifa}.
Here, however,
we will compare results of the unadorned subroutine \texttt{kirchfast()} \vpageref{/prog:kirchfast}
with our new programs, \texttt{phasemig()} \vpageref{/prog:phasemig} and \texttt{phasemod()} \vpageref{/prog:phasemod}.
Figure~\ref{fig:comrecon} shows the result of modeling data and then migrating it.
Kirchhoff and phase-shift migration methods both work well.
As expected, the Kirchhoff method lacks some of the higher frequencies
that could be restored by $\sqrt{-i\omega}$.
Another problem is the irregularity of the shallow bedding.
This is an operator aliasing problem
addressed in chapter~\ref{trimo/paper:trimo}.

\plot{comrecon}{width=6.00in,height=2.25in}{
        Reconstruction after modeling.
        Left is by the nearest-neighbor Kirchhoff method.
        Right is by the phase shift method.
        }
%\newslide

Figure~\ref{fig:phaspec} shows the temporal spectrum of the original sigmoid model,
along with the spectrum of the reconstruction via phase-shift methods.
We see the spectra are essentially identical
with little growth of high frequencies
as we noticed with the Kirchhoff method
in Figure~\ref{krch/fig:kirspec}.
\sideplot{phaspec}{width=3.00in,height=1.5in}{
        Top is the temporal spectrum of the model.
        Bottom is the spectrum of the reconstructed model.
        }
%\newslide

Figure~\ref{fig:commod} shows the effect of coarsening the space axis.
Synthetic data is generated from an increasingly subsampled model.
Again we notice that the phase-shift method of this chapter
produces more plausible results than
the simple Kirchhoff programs of chapter~\ref{krch/paper:krch}.

\plot{commod}{width=6.00in,height=6.75in}{
        Modeling with increasing amounts of lateral subsampling.
        Left is the nearest-neighbor Kirchhoff method.
        Right is the phase-shift method.
        Top has 200 channels,
        middle has 100 channels,
        and bottom has 50 channels.
        }
%\newslide


\subsection{Damped square root}
\sx{damped square root}
The definition of $k_z$ as $k_z=\sqrt{\omega^2/v^2 - k_x^2}$
obscures two aspects of $k_z$.
First, which of the two square roots is intended,
and second, what happens when $k_x^2 > \omega^2/v^2$.
For both coding and theoretical work we
need a definition of $ik_z$ that is valid
for both positive and negative values of $\omega$
and for all $k_x$.
Define a function $R=ik_z(\omega,k_z)$ by
\begin{equation}
        R \eq ik_z \eq
        \sqrt{(-i \omega + \epsilon)^2+k_x^2}
\end{equation}
It is important to know that for any $\epsilon >0$,
and any real $\omega$ and  real $k_x$ that
the real part $\Re R>0$ is positive.
This means we can extrapolate waves safely
with $e^{-Rz}$ for increasing $z$ or
with $e^{+Rz}$ for decreasing $z$.
To switch from downgoing to upcoming we use
the complex conjugate $\overline{R}$.
Thus we have disentangled the damping from the direction of propagation.
%In life, waves damp as they propagate, but in imaging
%we want damping even when we extrapolate backwards.
%Imaging is not inversion.

\begin{comment}
\par
Let us see why $\Re R>0$ is positive
for all real values of $\omega$ and $k_x$.
Recall that for $\omega$ ranging between $\pm\infty$,
$e^{i\omega\Delta t}$ rotates around the unit circle
in the complex plane.
Examine Figure \ref{fig:francis}
which shows the complex functions:
\newslide
\begin{enumerate}
\item       $f(\omega)= \epsilon - i\omega$,
\item       $ -i\hat\omega = (1+\epsilon) -e^{i\omega\Delta t}$,
\item       $(-i \hat\omega)^2$,
\item       $(ik_z)^2 =(-i \hat\omega)^2+k_x^2$, and
\item       $ik_z = [(-i \hat\omega)^2+k_x^2]^{1/2}$
\end{enumerate}
\activeplot{francis}{width=6.0in,height=1.15in}{ER} {
        Some functions in the complex plane.
        }
\newslide
The first two panels are explained by the first two functions.
The first two functions and the first two panels look different
but they become the same in the practical limit
of $\epsilon\rightarrow 0$ and $\Delta t \rightarrow 0$.
The left panel represents a time derivative in continuous time,
and the second panel likewise
in sampled time is for
a ``causal finite-difference operator''
representing a time derivative.
Notice that the graphs look the same near $\omega = 0$.
As we sample seismic data with increasing density,
$\Delta t \rightarrow 0$,
the frequency content shifts further away from the Nyquist frequency.
Measuring $\omega$ in radians/sample,
in the limit
$\Delta t \rightarrow 0$, the physical energy is all near
$\omega = 0$.

\par
The third panel in Figure \ref{fig:francis}
shows $(-i \hat\omega)^2$ which is a cardioid that
wraps itself close up to the negative imaginary axis without touching it.
(To understand the shape near the origin, think about the square
of the leftmost plane.  You may have seen examples
of the negative imaginary axis being a branch cut.)
In the fourth panel a small positive quantity $k_x^2$ is added which
shifts the cardioid to the right a bit.
Taking the square root gives the last panel
which shows the curve in the right half plane
thus proving
the important result we need,
that $\Re\; ik_z(\omega,k_z)>0$ for all real $\omega$.
It is also positive for all real $k_x$ because
any $k_x^2>0$ shifts the cardioid to the right.
The additional issue of time causality in forward modeling
is covered in IEI.
\par
Luckily the Fortran {\tt csqrt()} function assumes the phase of
the argument is between $\pm 180^\circ$
exactly as we need here.
Thus the square root itself will have a phase between $\pm 90^\circ$
as we require.
In applications, $\epsilon$
would typically be chosen proportional to the maximum time on the data.
Thus the mathematical expression $-i\omega + \epsilon$
might be rendered in Fortran as 
{\tt cmplx(qi,-omega)}
where
{\tt qi=1./tmax}
and the whole concept implemented as in function~\texttt{eiktau()} \vpageref{/prog:eiktau}.
Do not set {\tt qi=0} because then the {\tt csqrt()} function
cannot decipher positive from negative frequencies.

\progdex{eiktau}{exp ikz}
\end{comment}

\par
Finally, you might ask, why bother with all this careful theory
connected with the damped square root.
Why not simply abandon the evanescent waves?
%as done by the ``{\tt if}'' statement in subroutines
%{\tt phasemig()} and {\tt phasemod()}?
There are several reasons:
\begin{enumerate}
\item
The exploding reflector concept fails for evanescent waves
(when $\omega^2 < v^2k_x^2$).
Realistic modeling would have them damping with depth.
Rather than trying to handle them correctly we will make a choice,
either (1) to abandon evanescent waves effectively setting them to zero,
or (2) we will take them to be damping.
(You might notice that when we switch from downgoing to upgoing,
a damping exponential switches to a growing exponential,
but when we consider the adjoint of applying a damped exponential,
that adjoint is also a damped exponential.)
\par
I'm not sure if there is a practical difference between
choosing to damp evanescent waves or simply to set them to zero,
but there should be a noticable difference on synthetic data:
When a Fourier-domain amplitude drops abruptly
from unity to zero, we can expect a time-domain signal
that spreads widely on the time axis,
perhaps dropping off slowly as $1/t$.
We can expect a more concentrated pulse
if we include the evanescent energy, even though it is small.
I predict the following behavior:
Take an impulse; diffract it and then migrate it.
When evanescent waves have been truncated, I predict
the impulse is turned into a ``butterfly'' whose wings
are at the hyperbola asymptote.
Damping the evanescent waves, I predict,
gives us more of a ``rounded'' impulse.
\item
In a later chapter we will handle the $x$-axis by finite differencing
(so that we can handle $v(x)$.
There a stability problem will develop unless we begin
from careful definitions as we are doing here.
\item
Seismic theory includes an abstract mathematical concept
known as branch-line integrals.
Such theory is most easily understood beginning from here.
\end{enumerate}




\subsection{Adjointness and ordinary differential equations}

It is straightforward to adapt the simple programs
\texttt{phasemig()} \vpageref{/prog:phasemig} and
\texttt{phasemod()} \vpageref{/prog:phasemod}
to depth variable velocity.
As you might suspect, the two processes are adjoint to each other,
and for reasons given at the end of chapter~\ref{conj/paper:conj}
it is desirable to code them to be so.
With differential equations and their boundary conditions,
the concept of adjoint is more subtle than previous examples.
Thus, I postponed till here the development of adjoint code
for phase-shift migration.
This coding is a little strenuous,
but it affords a review of many basic concepts,
so we do so here. (Feel free to skip this section.)
It is nice to have a high quality code for this fundamental operation.

\par
Many situations in physics are expressed by the differential equation
\begin{equation}
{du \over dz} \ - \  i \alpha \, u \eq s(z)
\end{equation}
In the migration application,
$u(z)$ is the up-coming wave,
$\alpha=-\sqrt{\omega^2 /v^2 - k_x^2}$,
$s(z)$ is the exploding-reflector source.
We take the medium to be layered
($v$ constant in layers)
so that $\alpha$ is constant in a layer,
and we put the sources at the layer boundaries.
Thus within a layer we have $du / dz - i \alpha \, u = 0$
which has the solution 
\begin{equation}
u(z_k+\Delta z ) \eq u(z_k)\  e^{i \alpha \Delta z}
\end{equation}
For convenience, we use the ``\bx{delay operator}''
in the $k$-th layer
$Z_k= e^{-i \alpha \Delta z}$
so the delay of upward propagation is expressed by
$ u(z_k) = Z_k \, u(z_k+\Delta z )$.
(Since $\alpha$ is negative for upcoming waves,
$Z_k= e^{-i \alpha \Delta z}$
has a positive exponent which represents delay.)
Besides crossing layers, we must cross layer boundaries
where the (reflection) sources add to the upcoming wave.
Thus we have
\begin{equation}
u_{k-1} \eq Z_{k-1} u_k + s_{k-1}
\label{eqn:modeling}
\end{equation}
Recursive use of equation~(\ref{eqn:modeling}) across a medium
of three layers is expressed in matrix form as
\begin{equation}
\bold M \, \bold u \eq
\left[
        \begin{array}{cccc}
                1  & -Z_0  & .    &  .    \\
                .  &  1    & -Z_1 &  .    \\
                .  &  .    & 1    & -Z_2  \\
                .  &  .    &  .   &  1
        \end{array}
        \right] \ 
\left[
        \begin{array}{c}
                u_0 \\
                u_1 \\
                u_2 \\
                u_3
        \end{array}
        \right]
\eq
\left[
        \begin{array}{c}
                s_0 \\
                s_1 \\
                s_2 \\
                s_3
        \end{array}
        \right]
\eq \bold s
\label{eqn:matrecur}
\end{equation}

A recursive solution begins at the bottom
with $u_3=s_3$ and propagates upward.

\par
The adjoint (complex conjugate) of the delay operator $Z$
is the time advance operator $\bar Z$.
The adjoint of equation~(\ref{eqn:matrecur})
is given by

\begin{equation}
\bold M' \, \tilde {\bold s}  \eq
\left[
        \begin{array}{cccc}
                1         & .        &     .    &  .  \\