paper.tex

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

six coefficients. It consists of two columns, each column having three
coefficients and the second column being a reversed copy of the first
one. When filter~(\ref{eqn:2dpef2}) is used in data regularization
problems, it can occasionally cause undesired high-frequency
oscillations in the solution, resulting from the near-Nyquist zeroes
of the polynomial $B(Z_t)$. The oscillations are easily removed in
practice with appropriate low-pass filtering.
\par
In the next section, I address the problem of estimating the local
slope $\sigma$ with filters of form~(\ref{eqn:2dpef2}). Estimating
the slope is a necessary step for applying the finite-difference
plane-wave filters on real data.

\section{Slope estimation}

Let us denote by $\mathbf{C}(\mathbf{\sigma})$ the operator of convolving the
data with the 2-D filter $C(Z_t,Z_x)$ of equation~(\ref{eqn:2dpef2}),
assuming the local slope $\mathbf{\sigma}$ is known. In order to determine 
the slope, we can define the least-squares goal
\begin{equation}
  \label{eqn:ls}
  \mathbf{C}(\mathbf{\sigma}) \, \mathbf{d} \approx 0\;,
\end{equation}
where $\mathbf{d}$ is the known data and the approximate equality
implies that the solution is found by minimizing the power of the
left-hand side. Equations~(\ref{eqn:b3}) and~(\ref{eqn:b5}) show that
the slope $\mathbf{\sigma}$ enters in the filter coefficients in an
essentially non-linear way. However, one can still apply the linear
iterative optimization methods by an analytical linearization of
equation~(\ref{eqn:ls}). The linearization (also known as the Gauss-Newton
iteration) implies solving the linear system
\begin{equation}
  \label{eqn:linit}
  \mathbf{C}'(\mathbf{\sigma}_0) \, \Delta \mathbf{\sigma} \,
  \mathbf{d}  + \mathbf{C}(\mathbf{\sigma}_0) \, \mathbf{d} \approx 0
\end{equation}
for the slope increment $\Delta \mathbf{\sigma}$. Here $\mathbf{\sigma}_0$
is the initial slope estimate, and $\mathbf{C}'(\mathbf{\sigma})$ is a
convolution with the filter, obtained by differentiating the filter
coefficients of $\mathbf{C}(\mathbf{\sigma})$ with respect to
$\mathbf{\sigma}$. After system~(\ref{eqn:linit}) is solved, the initial
slope $\mathbf{\sigma}_0$ is updated by adding $\Delta \mathbf{\sigma}$ to
it, and one can solve the linear problem again. Depending on the
starting solution, the method may require several non-linear
iterations to achieve an acceptable convergence. 
\par
The slope $\sigma$ in equation~(\ref{eqn:linit}) does not have to be
constant. We can consider it as varying in both time and space
coordinates.  This eliminates the need for local windows but may lead
to undesirably rough (oscillatory) local slope estimates.  Moreover,
the solution will be undefined in regions of unknown or constant data,
because for these regions the local slope is not constrained.  Both
these problems are solved by adding a regularization (styling) goal to
system~(\ref{eqn:linit}). The additional goal takes the form
\begin{equation}
  \label{eqn:regs}
  \epsilon \mathbf{D} \, \Delta \mathbf{\sigma} \approx 0\;,
\end{equation}
where $\mathbf{D}$ is an appropriate roughening operator and $\epsilon$
is a scaling coefficient. For simplicity, I chose $\mathbf{D}$ to be the
gradient operator. More efficient and sophisticated helical
preconditioning techniques are available 
\cite[]{GEO63-05-15321541,Fomel.sepphd.107,fandc}.
\par
In theory, estimating two different slopes $\mathbf{\sigma}_1$ and
$\mathbf{\sigma}_2$ from the available data is only marginally more
complicated than estimating a single slope. The convolution operator
becomes a cascade of $\mathbf{C}(\mathbf{\sigma}_1)$ and
$\mathbf{C}(\mathbf{\sigma}_2)$, and the linearization yields
\begin{equation}
  \label{eqn:lin2}
  \mathbf{C}'(\mathbf{\sigma}_1) \, \mathbf{C}(\mathbf{\sigma}_2) \, 
  \Delta \mathbf{\sigma}_1\, \mathbf{d} + \mathbf{C}(\mathbf{\sigma}_1) \, 
  \mathbf{C}'(\mathbf{\sigma}_2) \,
  \Delta \mathbf{\sigma}_2 \, \mathbf{d} + \mathbf{C}(\mathbf{\sigma}_1) \,
  \mathbf{C}(\mathbf{\sigma}_2) \, \mathbf{d} \approx 0\;.
\end{equation}
The regularization condition should now be applied to both $\Delta
\mathbf{\sigma}_1$ and $\Delta \mathbf{\sigma}_2$:
 \begin{eqnarray}
  \label{eqn:reg1}
  \epsilon \mathbf{D} \, \Delta \mathbf{\sigma}_1 & \approx & 0\;; \\
  \label{eqn:reg2}
  \epsilon \mathbf{D} \, \Delta \mathbf{\sigma}_2 & \approx & 0\;.
\end{eqnarray}
The solution will obviously depend on the initial values of
$\mathbf{\sigma}_1$ and $\mathbf{\sigma}_2$, which should not be equal to
each other. System~(\ref{eqn:lin2}) is generally underdetermined,
because it contains twice as many estimated parameters as equations:
The number of equations corresponds to the grid size of the data
$\mathbf{d}$, while characterizing variable slopes $\sigma_1$ and
$\sigma_2$ on the same grid involves two gridded functions.  However,
an appropriate choice of the starting solution and the additional
regularization~(\ref{eqn:reg1}-\ref{eqn:reg2}) allow us to arrive at a
practical solution.
\par
The application examples of the next section demonstrate that when the
system of equations~(\ref{eqn:linit}-\ref{eqn:regs})
or~(\ref{eqn:lin2}-\ref{eqn:reg2}) are optimized in the least-squares
sense in a cycle of several linearization iterations, it leads to
smooth and reliable slope estimates. The regularization
conditions~(\ref{eqn:regs}) and~(\ref{eqn:reg1}-\ref{eqn:reg2}) assure
a smooth extrapolation of the slope to the regions of unknown or
constant data.

\section{Application examples}

In this section, I examine the performance of the finite-difference
plane-destruction filters on several test applications. The general
framework for applying these filters consists of the two steps:
\begin{enumerate}
\item Estimate the dominant local slope (or a set of local slopes)
  from the data. This step follows the least-squares optimization
  embedded in equations~(\ref{eqn:linit}) or~(\ref{eqn:lin2}). Thanks
  to the general regularization technique of
  equations~(\ref{eqn:regs} ) and~(\ref{eqn:reg1}-\ref{eqn:reg2}),
  locally smooth slope estimates are obtained without any need for
  breaking the data into local windows. Of course, local windows can
  be employed for other purposes (parallelization, memory management,
  etc.) Selecting appropriate initial values for the local slopes can
  speed up the computation and steer it towards desirable results.
  It is easy to incorporate additional constraints on the local 
  slope values.
\item Using the estimated slope, apply non-stationary plane-wave
  destruction filters for the particular application purposes. In the
  fault detection application, we simply look at the output of
  plane-wave destruction.  In the interpolation application, the
  filters are used to constrain the missing data. In the noise
  attenuation application, they characterize the coherent 
  signal and noise components in the data.
\end{enumerate}
A description of these particular applications follows next.

\subsection{Fault detection}

\inputdir{lomo}

The use of prediction-error filters in the problem of detecting local
discontinuities was suggested by \cite{SEG-1994-1572,gee}, and
further refined by \cite{schwabSEG} and
\cite{Schwab.sepphd.99}. \cite{SEG-1998-0653} used simple
plane-destruction filters in a similar setting to compute coherency
attributes.
\par
To test the performance of the improved plane-wave destructors, I
chose several examples from \cite{gee}.
Figure~\ref{fig:sigmoid-txr} introduces the first example. The left
plot of the figure shows a synthetic model, which resembles
sedimentary layers with a plane unconformity and a curvilinear fault.
The model contains 200 traces of 200 samples each.
The right plot shows the corresponding \emph{texture}
\cite[]{EAE-1999-1009}, obtained by convolving a field of random
numbers with the inverse of plane-wave destruction filters. The
inverses are constructed using helical filtering techniques
\cite[]{GEO63-05-15321541,Fomel.sepphd.107}. Texture plots allow us to
quickly access the ability of the destruction filters to characterize
the main locally plane features in the data.  The dip field was
estimated by the linearization method of the previous section. The dip
field itself and the prediction residual [the left-hand side of
equation~(\ref{eqn:ls})] are shown in the left and right plots of
Figure~\ref{fig:sigmoid-dip} respectively. We observe that the
texture plot does reflect the dip structure of the input data, which
indicates that the dip field was estimated correctly. The fault and
unconformity are clearly visible both in the dip estimate and in the
residual plots. Anywhere outside the slope discontinuities and the
boundaries, the residual is close to zero.  Therefore, it can be used
directly as a fault detection measure.  Comparing the residual plot in
Figure~\ref{fig:sigmoid-dip} with the analogous plot of
\cite{SEG-1994-1572,gee}, reproduced in
Figure~\ref{fig:sigmoid-clae}, establishes a superior performance of
the improved finite-difference destructors in comparison with that of
the local $T$-$X$ prediction-error filters.

\plot{sigmoid-txr}{width=6in,height=3.5in}{Synthetic sedimentary
  model. Left plot: Input data. Right plot: Its texture. The texture is
  computed by convolving a field number with the inverse of plane-wave destruction 
  filters. It highlights the position of estimated local plane waves.}
\plot{sigmoid-dip}{width=6in,height=3.5in}{Synthetic sedimentary
  model. Left plot: Estimated dip field. Right plot: Prediction
  residual. Large absolute residual indicates the location of faults.}  
\sideplot{sigmoid-clae}{width=2.7in,height=3.5in}{Prediction
  residual of the 11-point prediction-error filter estimated in local
  20x6 windows (reproduced from \cite[]{gee}).To be compared with the
  right plot in Figure~\ref{fig:sigmoid-dip}.}
\par
The left plot in Figure~\ref{fig:conflict-txr} introduces a simpler
synthetic test. The model is composed of linear events with two
conflicting slopes. A regularized dip field estimation attempts to
smooth the estimated dip in the places where it is not constrained by
the data (the left plot of Figure~\ref{fig:conflict-dip}.) The effect
of smoothing is clearly seen in the texture image (the right plot in
Figure~\ref{fig:conflict-txr}). The corresponding residual (the right
plot of Figure~\ref{fig:conflict-dip}) shows suppressed linear events
and highlights the places of their intersection. Residuals are large at 
intersections because a single dominant dip model fails to adequately
represent both conflicting dips.

\plot{conflict-txr}{width=6in,height=3.5in}{Conflicting dips
  synthetic. Left plot: Input data. Right plot: Its texture.}
\plot{conflict-dip}{width=6in,height=3.5in}{Conflicting dips
  synthetic. Left plot: Estimated dip field. Right plot: Prediction
  residual.Large absolute residual indicates the location of
  conflicting dips.}

%\begin{comment}
\par
The left plot in Figure~\ref{fig:yc27-txr} shows a real shot gather:
a portion of \cite{yilmaz} data set 27. The initial dip in the dip
estimation program was set to zero. Therefore, the texture image (the
right plot in Figure~\ref{fig:yc27-txr}) contains zero-dipping plane
waves in the places of no data. Everywhere else the dip is accurately
estimated from the data. The data contain a missing trace at about
0.7~km offset and a slightly shifted (possibly mispositioned) trace at
about 1.1~km offset. The mispositioned trace is clearly visible in the
dip estimate (the left plot in Figure~\ref{fig:yc27-dip}), and the
missing trace is emphasized in the residual image (the right plot in
Figure~\ref{fig:yc27-dip}). Additionally, the residual image reveals
the forward and back-scattered surface waves, hidden under more
energetic reflections in the input data.

\plot{yc27-txr}{angle=90,totalheight=8.7in,width=5.075in}{Real 
  shot gather. Left plot: Input data. Right plot: Its texture.}
\plot{yc27-dip}{angle=90,totalheight=8.7in,width=5.075in}{Real 
  shot gather. Left plot: Estimated dip field. Right plot: 
  Prediction residual. The residual highlights surface waves 
  hidden under dominant reflection events in the original data.}
\par
%\end{comment}

Figure~\ref{fig:dgulf-txr} shows a stacked time section from the Gulf of
Mexico and its corresponding texture. The texture plot demonstrates
that the estimated dip (the left plot of Figure~\ref{fig:dgulf-dip})
reflects the dominant local dip in the data. After the plane waves
with the dominant dip are removed, many hidden diffractions appear in the
residual image (the right plot in Figure~\ref{fig:dgulf-dip}.) The
enhanced diffraction events can be used, for example, for 
estimating the medium velocity \cite[]{GEO49-11-18691880}.

\plot{dgulf-txr}{angle=90,totalheight=8.7in,width=5.075in}{Time
  section from the Gulf of Mexico. Left plot: Input data. Right plot:
  Its texture. The texture plot shows dominant local dips estimated
  from the data.}
\plot{dgulf-dip}{angle=90,totalheight=8.7in,width=5.075in}{Time
  section from the Gulf of Mexico. Left plot: Estimated dip field.
  Right plot: Prediction residual.The residual highlights diffraction events
  hidden under dominant reflections in the original data. }
\par
Overall, the examples of this subsection show that the
finite-difference plane-wave destructors provide a reliable tool for
enhancement of discontinuities and conflicting slopes in seismic
images. The estimation step of the fault detection procedure produces
an image of the local dominant dip field, which may have its own
interpretational value. An extension to 3-D is possible, as outlined
by \cite{Schwab.sepphd.99}, \cite{Clapp.sepphd.106}, and
\cite{Fomel.sepphd.107}.