paper.tex

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

we need to take into account the connection between the traveltime
derivatives in~(\ref{eqn:ampint}) and the geometry of the reflector.
As follows directly from the trigonometry of the incident and
reflected rays triangle (Figure \ref{fig:ocoray}),
\begin{eqnarray}
h & = & {r-s \over 2}=
D\,{{\cos{\alpha}\,\sin{\gamma}\,\cos{\gamma}} \over
{\cos^2{\alpha}-\sin^2{\gamma}}}\;,
\label{eqn:geoh}\\
y & = & {r+s \over 2}=
x+D\,{{\cos^2{\alpha}\,\sin{\alpha}} \over
{\cos^2{\alpha}-\sin^2{\gamma}}}\;,
\label{eqn:geoy}\\
y_0 & = & x+D\,\sin{\alpha}\;,
\label{eqn:geoy0}
\end{eqnarray}
where $D$ is the length of the normal ray. Let $\tau_0=2\,D/v$ be the
zero-offset reflection traveltime. Combining
equations~(\ref{eqn:geoh}) and~(\ref{eqn:geoy0}) with
(\ref{eqn:length}), we can get the following relationship:
\begin{equation}
a={\tau_n\over\tau_0}={{\cos{\alpha}\,\cos{\gamma}}\over
\left(\cos^2{\alpha}-\sin^2{\gamma}\right)^{1/2}}=
\left(1+{{\sin^2{\alpha}\,\sin^2{\gamma}}\over
{\cos^2{\alpha}-\sin^2{\gamma}}}\right)^{1/2}=
{h\,\over\sqrt{h^2-\left(y-y_0\right)^2}}\;,
\label{eqn:A}
\end{equation}
which describes the ``DMO smile''~(\ref{eqn:smilezott}) found by
\cite{GPR29-03-03740406} in geometric terms.
Equation~(\ref{eqn:A}) allows for a convenient change of variables
in equation~(\ref{eqn:ampint}). Let the reflection angle $\gamma$ be a
parameter monotonically increasing along a time ray. In this case,
each time ray is uniquely determined by the position of the reflection
point, which in turn is defined by the values of $D$ and $\alpha$.
According to this change of variables, we can
differentiate~(\ref{eqn:A}) along a time ray to get
\begin{equation}
{{d\tau_n}\over\tau_n}=-{{\sin^2{\alpha}}\over
{2\,\cos^2{\gamma}\,\left(\cos^2{\gamma}-\sin^2{\alpha}\right)}}\,
d\left(\cos^2{\gamma}\right)\;.
\label{eqn:dt2tg}
\end{equation}
Note also that the quantity $h\,\left(\tau_n\,{\partial \tau_n \over
    \partial h} \right)^{-1}$ in equation~(\ref{eqn:ampint}) coincides
exactly with the time ray invariant $C_3$ found in
equation~(\ref{eqn:abc}). Therefore its value is constant along each
time ray and equals
\begin{equation}
h\,\left(\tau_n\,{\partial \tau_n \over \partial h}\right)^{-1}=
-{v^2 \over 4\, \sin^2{\alpha}}\;.
\label{eqn:c3}
\end{equation}
Finally, as shown in Appendix~\ref{chapter:deriv},
\begin{equation}
\tau_n\,\left({\partial^2 \tau_n \over \partial y^2}-
{\partial^2 \tau_n \over \partial h^2}\right)=4\,
{\cos^2{\gamma}\over v^2}\,\left({\sin^2{\alpha}+DK}\over
{\cos^2{\gamma}+DK}\right)\;,
\label{eqn:curve}
\end{equation}
where $K$ is the reflector curvature at the reflection point.
Substituting (\ref{eqn:dt2tg}),~(\ref{eqn:c3}), and~(\ref{eqn:curve})
into~(\ref{eqn:ampint}) transforms the integral to the form
\begin{eqnarray*}
 \int_{t_o}^{t_n}\left[
h\,\left({\partial^2 \tau_n \over \partial y^2}-
{\partial^2 \tau_n \over \partial h^2}\right)\,
\left(\tau_n\,{\partial \tau_n \over \partial h} \right)^{-1}\right]\,
d\tau_n = 
\end{eqnarray*}
\begin{equation}
 = -{1 \over 2}\,\int_{\cos^2{\gamma_0}}^{\cos^2{\gamma}}
\left({1 \over {\cos^2{\gamma'}-\sin^2{\alpha}}}-
{1 \over
{\cos^2{\gamma'}+DK}}\right)\,d\left(\cos^2{\gamma'}\right) 
\label{eqn:int2g}
\end{equation}
which we can evaluate analytically. The final equation for the
amplitude transformation is
\begin{eqnarray}
A_n & = & A_0\,{\sqrt{\cos^2{\gamma}-\sin^2{\alpha}}\over
\sqrt{\cos^2{\gamma_0}-\sin^2{\alpha}}}\,\left(
{\cos^2{\gamma_0}+DK}\over{\cos^2{\gamma}+DK}\right)^{1/2}=
\nonumber \\
& = & A_0\,{{\tau_0\,\cos{\gamma}}\over{\tau_n\,\cos{\gamma_0}}}\,
\left(
{\cos^2{\gamma_0}+DK}\over{\cos^2{\gamma}+DK}\right)^{1/2}\;. 
\label{eqn:ampcurve}
\end{eqnarray}
In case of a plane reflector, the curvature $K$ is zero, and
equation~(\ref{eqn:ampcurve}) coincides with~(\ref{eqn:ampplane}).
In the general case can be rewritten as
\begin{equation}
A_n={{c\,\cos{\gamma}}\over{\tau_n\,\sqrt{\cos^2{\gamma}+DK}}}\;,
\label{eqn:ampray}
\end{equation}
where $c$ is constant along each time ray (it may vary with the reflection point
location on the reflector but not with the offset). We should compare equation
(\ref{eqn:ampray}) 
with the known expression for the reflection wave amplitude of the leading
ray series term in 2.5-D media \cite[]{cwp}:
\begin{equation}
A={{C_R(\gamma) \Psi}\over G}\;,
\label{eqn:amptrue}
\end{equation}
where $C_R$ stands for the angle-dependent reflection coefficient, $G$ is the
geometric spreading
\begin{equation}
G=v \tau {\sqrt{\cos^2{\gamma}+DK}\over \cos{\gamma}}\;,
\label{eqn:GS}
\end{equation}
and $\Psi$ includes other possible factors (such as the source
directivity) that we can either correct or neglect in the preliminary
processing.  It is evident that the curvature dependence of the
amplitude transformation (\ref{eqn:ampray}) coincides completely with
the true geometric spreading factor (\ref{eqn:GS}) and that the angle
dependence of the reflection coefficient is not accounted for the offset
continuation process. If the wavelet shape of the reflected wave on
seismic sections [$R_n$ in equation~(\ref{eqn:raymethod})] is
described by the delta function, then, as follows from the known
properties of this function,
\begin{equation}
A\,\delta\left(t-\tau(y,h)\right)=\left|{{dt_n} \over {dt}}\right|\,
A\,\delta\left(t_n-\tau_n(y,h)\right) =
{t \over t_n}\,A\,\delta\left(t_n-\tau_n(y,h)\right)\;, 
\label{eqn:deltafun} 
\end{equation}
which leads to the equality
\begin{equation}
A_n=A\,{t \over t_n}\;.
\label{eqn:a2an} 
\end{equation}
Combining equation~(\ref{eqn:a2an}) with equations~(\ref{eqn:amptrue})
and~(\ref{eqn:ampray}) allows us to evaluate the amplitude after
continuation from some initial offset $h_0$ to another offset $h_1$,
as follows:
\begin{equation}
A_1={{C_R(\gamma_0) \Psi_0}\over G_1}\;.
\label{eqn:ampfin}
\end{equation}
According to equation~(\ref{eqn:ampfin}), the OC process described by
equation (\ref{eqn:OCequation}) is amplitude-preserving in the sense
that corresponds to the definition of Born DMO
\cite[]{born,GEO56-02-01820189}. This means that the geometric spreading
factor from the initial amplitudes is transformed to the true
geometric spreading on the continued section, while the reflection
coefficient stays the same. This remarkable dynamic property allows
AVO (amplitude versus offset) analysis to be performed by a dynamic
comparison between true constant-offset sections and the sections
transformed by OC from different offsets.  With a simple trick, the
offset coordinate is transferred to the reflection angles for the AVO
analysis.  As follows from~(\ref{eqn:A}) and~(\ref{eqn:length}),
\begin{equation}
{\tau_n^2 \over \tau\,\tau_0}=\cos{\gamma}\;.
\label{eqn:AVO}
\end{equation}
If we include the ${t_n^2 \over t\,t_0}$ factor in the DMO operator
(continuation to zero offset) and divide the result by the DMO section
obtained without this factor, the resultant amplitude of the reflected
events will be directly proportional to $\cos{\gamma}$, where the
reflection angle $\gamma$ corresponds to the initial offset. Of
course, this conclusion is rigorously valid for constant-velocity
2.5-D media only.

\cite{GEO58-01-00470066} suggest a definition of true-amplitude DMO different
from that of Born DMO. The difference consists of two important
components:
\begin{enumerate}
\item {\em True-amplitude DMO addresses preserving the peak amplitude
    of the image wavelet instead of preserving its spectral density.}
  In the terms of this paper, the peak amplitude corresponds to the
  pre-NMO amplitude $A$ from formula~(\ref{eqn:amptrue}) instead of
  corresponding to the spectral density amplitude $A_n$.  A simple
  correction factor $t \over t_n$ would help us take the difference
  between the two amplitudes into account. Multiplication by $t \over
  t_n$ can be easily done at the NMO stage.
\item {\em Seismic sections are multiplied by time to correct for the
    geometric spreading factor prior to DMO (or, in our case, offset
    continuation) processing.} \end{enumerate} As follows
from~(\ref{eqn:GS}), multiplication by $t$ is a valid geometric
spreading correction for plane reflectors only.  It is the
amplitude-preserving offset continuation based on the OC
equation~(\ref{eqn:OCequation}) that is able to correct for the
curvature-dependent factor in the amplitude. To take into account the
second aspect of Black's definition, we can consider the modified
field $\hat{P}$ such that
\begin{equation}
\hat{P}\left(y,h,t_n\right)=t\,P\left(y,h,t_n\right)\;.
\label{eqn:p2pt}
\end{equation}
Substituting~(\ref{eqn:p2pt}) into the OC equation~(\ref{eqn:OCequation}) transforms the
latter to the form
\begin{equation}
h \, \left( {\partial^2 \hat{P} \over \partial y^2} - {\partial^2 \hat{P} \over
\partial h^2} \right) \, = \, t_n \, {\partial^2 \hat{P} \over {\partial t_n \,
\partial h}}\, -{\partial \hat{P} \over \partial h}\; . 
\label{eqn:BSZequation} 
\end{equation}
Equations~(\ref{eqn:BSZequation}) and~(\ref{eqn:OCequation}) differ
only with respect to the first-order damping term $\partial \hat{P}
\over \partial h$.  This term affects the amplitude behavior but not
the traveltimes, since the eikonal-type equation~(\ref{eqn:eikonal})
depends on the second-order terms only.  Offset continuation operators
based on~(\ref{eqn:BSZequation}) conform to Black's definition of
true-amplitude processing.

\cite{menorm} describe an alternative approach to
confirming the kinematic and amplitude validity of the offset
continuation equation. Applying equation~(\ref{eqn:OCequation})
directly on the Kirchhoff model of prestack seismic data shows that
the equation is satisfied to the same asymptotic order of accuracy
as the Kirchhoff modeling approximation \cite[]{haddon,norm}.

\section{Integral offset continuation operator}
%%%%%%%%%%%%%%%%%%%%%%
Equation~(\ref{eqn:OCequation}) describes a continuous process of
reflected wavefield continuation in the time-offset-midpoint domain.
In order to find an integral-type operator that performs the one-step
offset continuation, I consider the following initial-value
problem for equation~(\ref{eqn:OCequation}):

{\em Given a post-NMO constant-offset section at half-offset $h_1$}
\begin{equation}
\left.P(t_n,h,y)\right|_{h=h_1}=P^{(0)}_1(t_n,y) 
\label{eqn:bound0} 
\end{equation}
{\em and its first-order derivative with respect to offset}
\begin{equation}
\left.\partial P(t_n,h,y)\over \partial h\right|_{h=h_1}=P^{(1)}_1(t_n,y)\;, 
\label{eqn:bound1} 
\end{equation}
{\em find the corresponding section $P^{(0)}(t_n,y)$ at offset $h$.}

Equation~(\ref{eqn:OCequation}) belongs to the hyperbolic type, with
the offset coordinate $h$ being a ``time-like'' variable and the
midpoint coordinate $y$ and the time $t_n$ being ``space-like''
variables.  The last condition~(\ref{eqn:bound1}) is required for the
initial value problem to be well-posed \cite[]{kurant}. From a physical
point of view, its role is to separate the two different wave-like
processes embedded in equation~(\ref{eqn:OCequation}), which are
analogous to inward and outward wave propagation. We will associate
the first process with continuation to a larger offset and the second
one with continuation to a smaller offset.  Though the offset
derivatives of data are not measured in practice, they can be
estimated from the data at neighboring offsets by a finite-difference
approximation. Selecting a propagation branch explicitly, for example
by considering the high-frequency asymptotics of the continuation
operators, can allow us to eliminate the need for
condition~(\ref{eqn:bound1}). In this section, I discuss the exact
integral solution of the OC equation and analyze its asymptotics.

The integral solution of problem~(\ref{eqn:bound0}-\ref{eqn:bound1})
for equation~(\ref{eqn:OCequation}) is obtained in with the help of
the classic methods of mathematical physics
\cite[]{me,Fomel.sepphd.107}. It takes the explicit form
\begin{eqnarray}
P(t_n,h,y) & = &
\int\!\!\int P^{(0)}_1(t_1,y_1)\,G_0(t_1,h_1,y_1;t_n,h,y)\,dt_1\,dy_1
\nonumber \\ 
+ & & 
\int\!\!\int P^{(1)}_1(t_1,y_1)\,G_1(t_1,h_1,y_1;t_n,h,y)\,dt_1\,dy_1\;,
\label{eqn:integral} 
\end{eqnarray}
where the Green's functions $G_0$ and $G_1$ are expressed as
\begin{eqnarray}
G_0(t_1,h_1,y_1;t_n,h,y) & = & \mbox{sign}(h-h_1)\,{H(t_n) \over \pi}\,
{\partial \over \partial t_n}\,\left\{
H(\Theta) \over 
\sqrt{\Theta}\right\}\;,
\label{eqn:green0} \\
G_1(t_1,h_1,y_1;t_n,h,y) & = & \mbox{sign}(h-h_1)\,
{H(t_n) \over \pi}\,h\,
{t_n \over t_1^2}\,\left\{
H(\Theta) \over 
\sqrt{\Theta}\right\}\;,
\label{eqn:green1} 
\end{eqnarray}
and the parameter $\Theta$ is
\begin{equation}
\Theta(t_1,h_1,y_1;t_n,h,y)  = 
\left(h_1^2/t_1^2-h^2/t_n^2\right)\,\left(t_1^2-t_n^2\right)-
\left(y_1-y\right)^2\;.
\label{eqn:gamma}
\end{equation}
$H$ stands for the Heaviside step-function. 
 
From equations~(\ref{eqn:green0}) and~(\ref{eqn:green1}) one can see
that the impulse response of the offset continuation operator is
discontinuous in the time-offset-midpoint space on a surface defined
by the equality
\begin{equation}
\Theta(t_1,h_1,y_1;t_n,h,y)  = 0\;,
\label{eqn:conoid}
\end{equation}
which describes the ``wavefronts'' of the offset continuation process.
In terms of the theory of characteristics \cite[]{kurant}, the surface
$\Theta=0$ corresponds to the characteristic conoid formed by the
bi-characteristics of equation~(\ref{eqn:OCequation}) -- time rays