paper.tex

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

\sx{convolution ! two-dimensional}
Convolution in two dimensions is just like convolution
in one dimension except that convolution is done on two axes.
The input and output data are planes of numbers
and the filter is also a plane.
A two-dimensional filter \sx{filter ! two-dimensional} is
a small plane of numbers that
is convolved over a big data plane of numbers.
\par
Suppose the data set is a collection of seismograms
uniformly sampled in space.
In other words, the data is numbers in a $(t,x)$-plane.
For example, the following filter
destroys any wavefront
aligned along the direction of a line containing both the ``+1''
and the ``$-1$''.
\begin{equation}
   \begin{array}{cc}
        -1     &\cdot  \\
        \cdot  &\cdot  \\
        \cdot  &1      \end{array}
\label{eqn:leftwave}
\end{equation}
The next filter destroys a wave with a slope
in the opposite direction:
\begin{equation}
   \begin{array}{cc}
        \cdot  &1     \\
        -1     &\cdot  \end{array}
\label{eqn:ritewave}
\end{equation}
To convolve the above two filters,
we can reverse either on (on both axes) and correlate them,
so that you can get
\begin{equation}
        \begin{array}{ccc}
              \cdot &  -1  &\cdot      \\
                 1  &\cdot &\cdot      \\
              \cdot &\cdot &  1        \\
              \cdot &  -1  &\cdot      \end{array}
\label{eqn:twowave}
\end{equation}
which destroys waves of both slopes.
\par
A \bx{two-dimensional filter}
\sx{filter ! two-dimensional}
that can be a \bx{dip-rejection filter} like (\ref{eqn:leftwave}) or (\ref{eqn:ritewave}) is
\sx{filter ! dip-rejection}
\begin{equation}
   \begin{array}{cc}
        a & \cdot  \\
        b & \cdot  \\
        c &   1    \\
        d & \cdot  \\
        e & \cdot
        \end{array}
%       a     &b     &c    &d     &e        \\
%       \cdot &\cdot &1    &\cdot &\cdot
\label{eqn:onedip}
\end{equation}
where the coefficients
$(a,b,c,d,e)$
are to be estimated by least squares in order
to minimize the power out of the filter.
(In the filter table,
the time axis runs vertically.)
\par
Fitting the filter to two neighboring traces
that are identical but for a time shift, we see that
the filter coefficients $(a,b,c,d,e)$ should turn out to be
something like $(-1,0,0,0,0)$ or
$(0,0,-.5,-.5, 0)$,
depending on the dip (stepout) of the data.
But if the two channels are not fully coherent, we expect to see
something like
$(-.9,0,0,0,0)$ or
$(0,0,-.4,-.4,0)$.
To find filters such as (\ref{eqn:twowave}),
we adjust coefficients to minimize the power out
of filter shapes, as in
\begin{equation}
        \begin{array}{ccccc}
                v  & a & \cdot \\
                w  & b & \cdot \\
                x  & c &   1    \\
                y  & d & \cdot \\
                z  & e & \cdot 
%                v  &  w   &   x  &   y  &   z      \\
%                a  &  b   &   c  &   d  &   e     \\
%             \cdot &\cdot &   1  &\cdot &\cdot
        \end{array}
\end{equation}
\par
With 1-dimensional filters,
we think mainly of power spectra,
and with 2-dimensional filters
we can think of temporal spectra and spatial spectra.
What is new, however,
is that in two dimensions we can think of dip spectra
(which is when a 2-dimensional spectrum has a particularly common form,
namely when energy organizes on radial lines in the $(\omega,k_x)$-plane).
% A line of energy in the $(\omega,k_x)$-plane
% when seen as a function of $\omega$ for a fixed $k_x$ is an impulse.
% Fourier transforming $k_x$ to $x$ and then at fixed $x$
% at a function of $\omega$ we do not see sinusoids.
As a short (three-term) 1-dimensional filter can devour a sinusoid,
we have seen that simple 2-dimensional filters can devour
a small number of dips.

\subsection{PEF whiteness proof in 2-D}
\inputdir{whitepruf}
\par
\sx{whiteness ! multidimensional}
A well-known property (see FGDP or PVI)
of a 1-D PEF is that its energy clusters immediately after the
impulse at zero delay time.
Applying this idea to
the helix in Figure %\CHAPFIG{hlx}{sergey-helix}
shows us that we can consider a 2-D PEF
to be a small halfplane like
%\CHAPEQN{hlx}{2dpef}
with an impulse along a side.
These shapes are what we see here in
Figure~\ref{fig:whitepruf}.

\sideplot{whitepruf}{height=2.0in}{
  A 2-D whitening filter template, and itself lagged.
  At output locations ``A'' and ``B,'' the filter coefficient
  is constrained to be ``1''.
  When the semicircles are viewed as having infinite radius,
  the B filter is contained in the A filter.
  Because the output at A is orthogonal to all its inputs,
  which include all inputs of B,
  the output at A is orthogonal to the output of B.
}

\par
Figure~\ref{fig:whitepruf} shows the input plane with a 2-D filter on top
of it at two possible locations.
The filter shape is a semidisk,
which you should imagine being of
infinitely large radius.
Notice that semidisk A includes all the points in B.
The output of disk A will be shown to be orthogonal to the output
of disk B.
Conventional least squares theory says that the coefficients of the filter
are designed so that the output of the filter
is orthogonal to each of the inputs to that filter
(except for the input under the ``1,''
because any nonzero signal cannot be orthogonal to itself).
Recall that if a given signal is orthogonal to each in a given group of signals,
then the given signal is orthogonal
to all linear combinations within that group.
The output at B is a linear combination of members
of its input group,
which is included in the input group of A,
which are already orthogonal to A.
Therefore the output at B is orthogonal to the output at A.
In summary,
\par            % don't remove this line
\begin{tabular}{lll}
residual     & $\perp$ &  fitting function \\
output at A  & $\perp$ &  each input to A \\
output at A  & $\perp$ &  each input to B \\
output at A  & $\perp$ &  linear combination of each input to B \\
output at A  & $\perp$ &  output at B
\end{tabular}
\par            % don't remove this line
\noindent
The essential meaning is that
a particular lag of the output \bx{autocorrelation} function vanishes.

\par
Study Figure~\ref{fig:whitepruf} to see for what lags
all the elements of the B filter are wholly contained in the A filter.
These are the lags
where we have shown the output autocorrelation to be vanishing.
Notice another set of lags where we have proven nothing
(where B is moved to the right of A).
Autocorrelations are centrosymmetric,
which means that the value at any lag
is the same as the value at the negative of that lag,
even in 2-D and 3-D where the lag is a vector quantity.
Above we have shown that a halfplane of autocorrelation values vanishes.
By the centrosymmetry, the other half must vanish too.
Thus the autocorrelation of the PEF output is an impulse function,
so its 2-D spectrum is white.
\sx{white 2-D spectrum}
\sx{spectrum ! white 2-D}

\par
The helix tells us why the proper filter form
is not a square with the ``1'' on the corner.
Before I discovered the helix, I understood it another way
(that I learned from John P. Burg):
For a spectrum to be white,
{\it all}
nonzero autocorrelation lags must be zero-valued.
If the filter were a quarter-plane,
then the symmetry of autocorrelations
would only give us vanishing in another quarter,
so there would be two remaining quarter-planes
where the autocorrelation was not zero.

\par
Fundamentally,
the white-output theorem requires a
one-dimensional ordering to the values in a plane or volume.
The filter must contain a halfplane of values
so that symmetry gives the other half.

\par
You will notice some nonuniqueness.
We could embed the helix
with a $90^\circ$ rotation
in the original physical application.
Besides the difference in side boundaries,
the 2-D PEF would have a different orientation.
Both PEFs should have an output that tends to whiteness as
the filter is enlarged.
It seems that we could design whitening autoregression filters
for $45^\circ$ rotations also,
and we could also design them for hexagonal coordinate systems.
In some physical problems,
you might find the nonuniqueness unsettling.
Does it mean the ``final solution'' is nonunique?
Usually not, or not seriously so.
Recall even in one dimension, the time reverse of a PEF
has the same spectrum as the original PEF.
When a PEF is used for regularizing a fitting problem,
it is worth noticing that the quadratic form minimized
is the PEF times its adjoint so the phase drops out.
Likewise, a missing data restoration also amounts to minimizing
a quadratic form so the phase again drops out.

\subsection{Examples of modeling and deconvolving with a 2-D PEF }
\inputdir{morgan}
Here we examine elementary signal-processing applications of
2-D prediction-error filters (PEFs)
on both everyday 2-D textures and on seismic data.
Some of these textures are easily modeled with
prediction-error filters (PEFs) while others are not.
All figures used the same $10\times 10$ filter shape.
No attempt was made to optimize filter size
or shape or any other parameters.

\par
Results in Figures
\ref{fig:granite}-\ref{fig:WGstack}
are shown with various familiar textures\footnote{
	I thank Morgan Brown for finding these textures.
	}
on the left
as training data sets.
From these training data sets,
a prediction-error filter (PEF) is estimated
using module \texttt{pef} \vpageref{lst:pef}.
The center frame is simulated data made by deconvolving
(polynomial division) random numbers by the estimated PEF.
The right frame is the more familiar process,
convolving the estimated PEF on the training data set.

\plot{granite}{width=6.0in,height=1.75in}{
  Synthetic granite matches the training image quite well.
  The prediction error (PE) is large at grain boundaries
  so it almost seems to outline the grains.
  % The PE might be more interesting if I plotted its absolute value.
}

\par
Theoretically, the right frame tends towards a white spectrum.
Earlier you could notice
the filter size by knowing that the output
was taken to be zero where the filter is only partially on the data.
This was annoying on real data where we didn't want to throw
away any data around the sides.
Now the filtering is done without a call to the boundary module
so we have typical helix wraparound.

\plot{wood}{width=6.0in,height=1.75in}{
  Synthetic wood grain has too little white.
  This is because of the nonsymmetric brightness histogram of natural wood.
  Again, the PEF output looks random as expected.
}

\plot{herr}{width=6.0in,height=1.75in}{
  A banker's suit (left).  A student's suit (center).  My suit (right).
  The prediction error is large where the weave changes direction.
}

\plot{fabric}{width=6.0in,height=1.70in}{
  Basket weave.
  The simulated data fails to segregate the two dips into a checkerboard pattern.
  The PEF output looks structured perhaps because the filter is too small.
}

\plot{brick}{width=6.0in,height=1.70in}{
  Brick.
  Synthetic brick edges are everywhere
  and do not enclose blocks containing a fixed color.
  PEF output highlights the mortar.
}

\plot{ridges}{width=6.0in,height=1.70in}{
  Ridges.
  A spectacular failure of the stationarity assumption.
  All dips are present but in different locations.
  Never-the-less,
  the ridges have been sharpened by the deconvolution.
}

\plot{WGstack}{width=6.0in,height=3.5in}{
  Gulf of Mexico seismic section, modeled, and deconvolved.
  Do you see any drilling prospects in the simulated data?
  In the deconvolution, the strong horizontal layering
  is suppressed giving a better view of the hyperbolas.
  The decon filter is the same $10\times 10$ used on the everyday textures.
}


\par
Since a PEF tends to the inverse of the spectrum of its input,
results similar to these could probably be found
using Fourier transforms, smoothing spectra, etc.
We used PEFs because of their flexibility.
The filters can be any shape.
They can dodge around missing data, or we can use them
to estimate missing data.
We avoid periodic boundary assumptions inherent to FT.
The PEF's are designed only internal to known data, not off edges
so they are readily adaptible to nonstationarity.
Thinking of these textures as seismic time slices,
the textures could easily be required to pass thru specific
values at well locations.  



\subsection{Seismic field data examples}
\inputdir{pefex}
\par
Figures~\ref{fig:specdecon}-\ref{fig:zof}
are based on exploration seismic data from the Gulf of Mexico deep water.
A ship carries an air gun and tows a streamer with some hundreds of geophones.
First we look at a single pop of the gun.
We use all the geophone signals to create a single 1-D PEF for the time axis.
This changes the average temporal frequency spectrum
as shown in Figure~\ref{fig:specdecon}.
\plot{specdecon}{width=6in,height=2in}{
  $\omega$ spectrum of a shot gather
  of Figure~\ref{fig:decon0}
  before and after 1-D decon
  with a 30 point filter.
}
Signals from 60 Hz to 120 Hz are boosted substantially.
The raw data has evidently been prepared
with strong filtering against signals below about 8 Hz.
The PEF attempts to recover these signals, mostly unsuccessfully,
but it does boost some energy near the 8 Hz cutoff.
Choosing a longer filter would flatten the spectrum further.
The big question is, ``Has the PEF improved the appearance of the data?''

\par
The data itself from the single pop,
both before and after PE-filtering is shown in
Figure~\ref{fig:decon0}.
For reasons of esthetics of human perception
I have chosen to display a mirror image of the PEF'ed data.
To see a blink movie of superposition of
before-and-after images you need the electronic book.
We notice that signals of high temporal frequencies
indeed have the expected hyperbolic behavior in space.
Thus, these high-frequency signals are wavefields, not mere random noise.

\plot{decon0}{width=6in,height=8.4in}{
  Raw data with its mirror.