paper.tex

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

An important transformation in exploration geophysics
takes data as a function of shot-receiver offset
and transforms it to data as a function of apparent velocity.
Data is summed along hyperbolas of many velocities.
This important industrial process is adjoint to another that may
be easier to grasp:  data is synthesized by a superposition
of many hyperbolas.
The hyperbolas have various asymptotes (velocities) and various
tops (apexes).
Pseudocode for these transformations is
\par\noindent

\vbox{\begin{tabbing}
indent \= char \= char \= char \= char \= lotsofcharshereIhope \= \kill
\>do $v$ \{		\\
\>do $\tau$ \{ \\
\>do $x$  \{\\
\>	\>	\>$t = \sqrt{ \tau^2 + x^2/v^2 }$ \\
\>	\>	\>if hyperbola superposition \\
\>	\>	\>\>\> data$(t,x)$ = data$(t,x)$ + vspace$(\tau,v)$ \\
\>	\>	\>else if velocity analysis \\
\>	\>	\>\>\> vspace$(\tau,v)$ = vspace$(\tau,v)$ + data$(t,x)$\\
\>	\>	\}\}\}
\end{tabbing}}
\par\noindent
This pseudocode transforms one plane to another using the equation
$t^2 = \tau^2 +x^2/v^2$.  This equation relates four variables,
the two coordinates of the data space $(t,x)$
and the two of the model space $(\tau,v)$.
Suppose a model space is all zeros except for an impulse at $(\tau_0, v_0)$.
The code copies this inpulse to data space everywhere where
$t^2 = \tau_0^2 + x^2/v_0^2$.   In other words, the impulse
in velocity space is copied to a hyperbola in data space.
In the opposite case an impulse at a point in data space $(t_0,x_0)$
is copied to model space everywhere that satisfies the equation
$t_0^2 = \tau^2 + x_0^2/v^2$.  
Changing from velocity space to
slowness space this equation
$t_0^2 = \tau^2 + x_0^2 s^2$
has a name.   In $(\tau, s)$-space it is an ellipse
(which reduces to a circle when $x_0^2=1$.

%XX
\par
Look carefully in the model spaces of
Figure~\ref{fig:velvel} and
Figure~\ref{fig:mutevel}.
Can you detect any ellipses?
For each ellipse,
does it come from a large $x_0$ or a small one?
Can you identify the point $(t_0,x_0)$ causing the ellipse?



\par
We can ask the question, if we transform data to velocity space,
and then return to data space,
will we get the original data?
Likewise we could begin from the velocity space,
synthesize some data, and return to velocity space.
Would we come back to where we started?
The answer is yes, in some degree.
Mathematically, the question amounts to this:
Given the operator $\bold A$, is $\bold A'\bold A$ approximately
an identity operator, i.e.~is $\bold A$ nearly a unitary operator?
It happens that $\bold A'\bold A$ defined by the pseudocode above
is rather far from an identity transformation,
but we can bring it much closer
by including some simple scaling factors.
It would be a lengthy digression here to derive all these weighting factors
but let us briefly see the motivation for them.
One weight arises because waves lose \bx{amplitude} as they spread out.
Another weight arises because some angle-dependent effects should be taken
into account.  A third weight arises because in creating a velocity space,
the near offsets are less important than the wide offsets
and we do not even need the zero-offset data.
A fourth weight is a frequency dependent one
which is explained in chapter~\ref{ft1/paper:ft1}.
Basically, the summations in the velocity transformation are like integrations,
thus they tend to boost low frequencies.
This could be compensated by scaling
in the frequency domain
with frequency as $\sqrt{-i\omega}$
with subroutine \texttt{halfdifa()} \vpageref{/prog:halfdifa}.
% To remove reference to halfdifa if Section 6.5.2 is removed - Biondo 4/96

\par
The weighting issue will be examined in more detail later.
Meanwhile, we can see nice quality examples
from very simple programs
if we include the weights
in the physical domain, $w= \sqrt{1/t}\; \sqrt{x/v}\; \tau /t $.
(Typographical note:  Do not confuse
the weight $w$ (double you) with omega $\omega$.)
To avoid the coding clutter of the frequency domain weighting
$\sqrt{-i\omega}$ I omit that,
thus getting smoother results than theoretically preferable.
Figure~\ref{fig:velvel} illustrates this smoothing by starting
from points in velocity space, transforming to offset,
and then back and forth again.

\plot{velvel}{width=5.50in}{
  Iteration between spaces.
  Left are model spaces.
  Right are data spaces.
  Right derived from left.
  Lower model space derived from upper data space.
}

\par
There is one final complication relating to weighting.
The most symmetrical approach is to put
$w$ into both $\bold A$ and $\bold A'$.
%This is what subroutine \texttt{velsimp()} \vpageref{/prog:velsimp} does.
Thus, because of the weighting by $\sqrt{x}$,
the synthetic data in Figure~\ref{fig:velvel} is
nonphysical.
An alternate view is to {\em  define} $\bold A$
(by the pseudo code above, or by some modeling theory)
and then for reverse transformation
use $w^2\bold A'$.

%\progdex{velsimp}{velocity spectra}

\par
An example %of applying subroutine \texttt{velsimp()} \vpageref{/prog:velsimp}
%to field data 
is shown in Figure~\ref{fig:mutvel}.

\plot{mutvel}{width=6.00in,height=3.6in}{
  Transformation of data as a function of offset (left)
  to data as a function of slowness (velocity scans)
  on the right
  using subroutine {\tt velsimp()}.
}

\subsection{Velocity picking}
\sx{velocity!picking}
For many kinds of data analysis,
we need to know the velocity of the earth as a function of depth.
To derive such information
we begin from Figure~\ref{fig:mutvel}
and draw a line through the maxima.
In practice this is often a tedious manual process,
and it needs to be done everywhere we go.
There is no universally accepted way to automate
this procedure, but we will consider one
that is simple enough that it can be fully described here,
and which works well enough for these demonstrations.
(I plan to do a better job later.)

\par
Theoretically we can define the velocity or slowness
as a function of traveltime depth by the moment function.
Take the absolute value of the data scans and smooth
them a little on the time axis to make something like an unnormalized
probability function, say $p(\tau ,s)>0$.
Then the slowness $s(\tau )$ could be defined by the moment function, i.e.,
\begin{equation}
s(\tau ) \eq { \sum_s \ s\  p(\tau ,s) \over \sum_s \ p(\tau ,s) }
\label{eqn:mommy}
\end{equation}
The problem with defining slowness $s(\tau )$ by the moment is that it is 
strongly influenced by noises away from the peaks,
particularly water velocity noises.
Thus, better results can be obtained if the sums in equation~(\ref{eqn:mommy})
are limited to a range about the likely solution.
To begin with, we can take the likely solution to be defined
by universal or regional experience.
It is sensible to begin from a one-parameter equation
for velocity increasing with depth where the form of the equation
allows a ray tracing solution
such as equation~(\ref{wvs/eqn:Vrms}).
Experience with Gulf of Mexico data shows that
$\alpha\approx 1/2\  {\rm sec}^{-1}$ is reasonable there
for equation~(\ref{wvs/eqn:Vrms}),
and that is the smooth curve
in Figure~\ref{fig:slowfit}.

\par
Experience with moments,
equation~(\ref{eqn:mommy}),
shows they are reasonable when
the desired result is near the guessed center of the range.
Otherwise, the moment is biased towards the initial guess.
This bias can be reduced in stages.
At each stage we shrink the width of the zone used to compute the moment.
%This procedure is used in subroutine \texttt{slowfit()} \vpageref{/prog:slowfit}
%which
%after smoothing to be described,
%gives the oscillatory curve you see in Figure~\ref{fig:slowfit}.

%\progdex{slowfit}{velocity est.}

A more customary way to view velocity space
is to square the velocity scans
and normalize them by the sum of the squares of the signals.
This has the advantage that the remaining information
represents velocity spectra
and removes variation due to seismic \bx{amplitude}s.
Since in practice, reliability seems somehow proportional to \bx{amplitude}
the disadvantage of normalization
is that reliability becomes more veiled.

\begin{comment}
\par
An appealing visualization of velocity is shown in the right side
of Figure~\ref{fig:slowfit}.
This was prepared from the absolute value of left side,
followed by filtering spatially with an antisymmetric
leaky integral function.
(See PVI page 57).
An example is shown on the right side of Figure~\ref{fig:slowfit}.

%\plot{slowfit}{width=6.00in,height=3.6in}{
  Left is the slowness scans.
  Right is the slowness scans after absolute value,
  smoothing a little in time,
  and antisymmetric leaky integration over slowness.
  Overlaying both is the line of slowness picks.
}
\end{comment}

\sideplot{fit}{width=3.00in,height=3.6in}{
  Slowness scans.
%  Right is the slowness scans after absolute value,
%  smoothing a little in time,
%  and antisymmetric leaky integration over slowness.
  Overlaying is the line of slowness picks.
}

%
% The plot below is not stable enough to include.
%
%\aKtivesideplot{vstack}{width=3.00in}{vscan}{
%	Stacking velocity.
%	Repeating the analysis of Figure~\protect\FIG{slowfit}
%	for hundreds of midpoint $y$-values gives this slowness
%	model, $s(t,y)$.
%	The velocity obviously varies with depth.
%	The validity of lateral variation is doubtful here,
%	though we will examine it more carefully later.
%	}

\subsection{Stabilizing RMS velocity}
\par
\sx{velocity!RMS}
\sx{velocity!interval}
With velocity analysis, we estimate the \RMS\ velocity.
Later we will need both the \RMS\ velocity and the \bx{interval velocity}.
(The word ``interval'' designates an interval between two reflectors.)
Recall from chapter~\ref{wvs/paper:wvs} equation~(\ref{wvs/eqn:vrmshyp})
$$
t^2 \eq \tau^2 + \frac{4h^2}{V^2(\tau)}  \nonumber
$$
\par
%Routine~\texttt{vint2rms()} \vpageref{/prog:vint2rms}
%converts from interval velocity to \RMS\ velocity
%and vice versa.
%\progdex{vint2rms}{interval to/from RMS vel}
The forward conversion follows 
in straightforward steps: square, integrate, square root.
The inverse conversion, like an adjoint,
retraces the steps of the forward transform
but it does the inverse at every stage.
There is however,
a messy problem with nearly all field data
that must be handled along the inverse route.
The problem is that the observed \RMS\ velocity function
is generally a rough function,
and it is generally unreliable over a significant portion of its range.
To make matters worse,
deriving an \bx{interval velocity} begins as does a derivative,
roughening the function further.
We soon find ourselves taking square roots of negative numbers,
which requires judgement to proceed.
\begin{comment}
The technique used in \texttt{vint2rms()} \vpageref{/prog:vint2rms}
is to average the squared interval velocity
in ever expanding neighborhoods until there are no longer
any negative squared interval velocities.
As long as we are restricting $v^2$ from being negative,
it is easy to restrict it to be above some allowable velocity,
say {\tt vminallow}.
Figures~\ref{fig:rufsmo} and~\ref{fig:vrmsint}
were derived from the velocity scans in Figure~\ref{fig:slowfit}. %
\activeplot{rufsmo}{width=5.00in,height=2.5in}{ER}{
	Left is the raw \RMS\ velocity.
	Right is a superposition of \RMS\ velocities,
	the raw one,
	and one constrained to have realistic interval velocities.
	}%
Figure~\ref{fig:rufsmo} shows the \RMS\ velocity before and after
a trip backward and forward through routine~\texttt{vint2rms()} \vpageref{/prog:vint2rms}.
The interval velocity associated with the smoothed velocity
is in figure~\ref{fig:vrmsint}.
	

\activesideplot{vrmsint}{width=2.50in,height=2.5in}{ER}{
	Interval velocity associated with the smoothed
	\RMS\ velocity of Figure~\protect\ref{fig:rufsmo}.
	Pushbutton allows experimentation with {\tt vminallow}.
	}
\end{comment}
\plot{wgvel1}{width=\textwidth}{Left is a superposition of \RMS\
  velocities, the raw one, and one constrained to have realistic
  interval velocities. Right is the nterval velocity.}