paper.tex

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

\subsection{Internal boundaries to multidimensional convolution }
\sx{convolution ! two-dimensional}

\par
Sometimes we deal with small patches of data.
In order that boundary phenomena not dominate
the calculation intended in the central region,
we need to take care that
input data is not assumed to be zero
beyond the interval that the data is given.

\par
The two little triangular patches of zeros in
the convolution matrix in equation (\ref{eqn:exmiss})
describe end conditions where it is assumed that the data $y_t$
vanishes before $t=1$ and after $t=6$.
Alternately we might not wish to make that assumption.
Thus the triangles filled with zeros could be regarded as missing data.
In this one-dimensional example,
it is easy to see that the filter, say
\verb#yy->mis#
should be set to \texttt{true} at the ends so no output
would ever be computed there.
We would like to find a general multidimensional algorithm
to correctly specify
\verb#yy->mis#
around the multidimensional boundaries.
This proceeds like the missing data algorithm,
i.e. we apply a filter of all ones to a data space template that is taken
all zeros except ones at the locations of missing data,
in this case $y_0,y_{-1}$ and $y_7,y_8$.
This amounts to surrounding the original data set with some missing data.
We need padding the size of the filter on all sides.
The padded region would be filled
with ones (designating missing inputs).
Where the convolution output is nonzero, there
\verb#yy->mis# is set to \texttt{true}
denoting an output with missing inputs.

\inputdir{XFig}
\par
The two-dimensional case is a little more cluttered than the 1-D case
but the principle is about the same.
Figure \ref{fig:rabdomain} shows a larger input domain,
a $5\times 3$ filter, and a smaller output domain.
\sideplot{rabdomain}{width=1.5in}{
  Domain of inputs and outputs of a two-dimensional filter like a PEF.
}
There are two things to notice.
First, sliding the filter everywhere inside the outer box,
we get outputs (under the 1 location) only in the inner box.
Second, (the adjoint idea) crosscorrelating the inner and outer boxes gives
us the $3\times 5$ patch of information
we use to build the filter coefficients.
We need to be careful not to assume
that signals vanish outside the region where they are defined.
In a later chapter we will break data spaces into overlapping patches,
separately analyze the patches, and put everything back together.
We do this because crosscorrelations change with time
and they are handled as constant in short time windows.
There we must be particularly careful that zero signal values not be presumed
outside of the small volumes;
otherwise the many edges and faces of the many small volumes
can overwhelm the interior that we want to study.

\par
In practice, the input and output are allocated equal memory,
but the output residual is initialized to zero everywhere
and then not computed
except where shown in figure \ref{fig:rabdomain}.
Below is module \texttt{bound}
to build a selector for filter outputs that should 
never be examined or even computed
(because they need input data from outside the given data space).
Inputs are a filter \texttt{aa}
and the size of its cube \texttt{na = (na(1),na(2),...)}.
Also input are two cube dimensions,
that of the data last used by the filter \texttt{nold} and
that of the filter's next intended use \texttt{nd}.
(\texttt{nold} and \texttt{nd} are often the same).
Module \texttt{bound}
begins by defining a bigger data space with room for a filter
surrounding the original data space \texttt{nd} on all sides.
It does this by the line \texttt{nb=nd+2*na}.
Then we allocate two data spaces
\texttt{xx} and \texttt{yy} of the bigger size \texttt{nb}
and pack many ones
in a frame of width \texttt{na} around the outside of \texttt{xx}.
The filter \texttt{aa} is also filled with ones.
The filter \texttt{aa} must be regridded for the bigger \texttt{nb}
data space (regridding merely changes the lag values of the ones).
Now we filter the input \texttt{xx} with \texttt{aa} getting \texttt{yy}.
Wherever the output is nonzero,
we have an output that has been affected by the boundary.
Such an output should not be computed.
Thus we allocate the logical mask \verb#aa->mis#
(a part of the helix filter definition
in module \texttt{helix} \vpageref{lst:helix})
and wherever we see a nonzero value of \texttt{yy}
in the output,
we designate the output as depending on missing inputs by setting
\verb#aa->mis# to \texttt{true}.
\moddex{bound}{out of bounds dependency}{30}{74}{user/gee}

\par
In reality one would set up the boundary conditions
with module
\texttt{bound}
before identifying locations of missing data
with module
\texttt{misinput}.
Both modules are based on the same concept,
but the boundaries are more cluttered and confusing
which is why we examined them later.


\subsection{Finding the prediction-error filter}
\sx{filter ! prediction-error}

\par
The first stage of the least-squares estimation
is computing the
\bx{prediction-error filter}.
The second stage will be using it to find the missing data.
The input data space contains a mixture of known data values
and missing unknown ones.
For the first stage of finding the filter,
we generally have many more fitting equations than we need
so we can proceed by ignoring the fitting equations
that involve missing data values.
We ignore them everywhere that the missing inputs hit the filter.

\par
The codes here do not address the difficulty that maybe
too much data is missing so that all weights are zero.
To add stabilization we could supplement the data volume
with a ``training dataset''
or by a ``prior filter''.
With things as they are,
if there is not enough data to specify a prediction-error filter,
you should encounter the error exit from \texttt{cgstep()} 
\vpageref{lst:cgstep}.
\moddex{pef}{estimate PEF on a helix}{28}{38}{user/gee}

%Results are shown in Figures
%\ref{fig:wood-hole}- \ref{fig:brick-hole}.
%Here again, the PEF is carefully designed using
%module \texttt{bound} \vpageref{lst:bound}
%so that no filter outputs are used where the filter
%is only partway on the data.
%After the PEF was designed, it is applied with the
%more cavalier helix treatment of boundaries
%so although you cannot see the frame of excluded outputs,
%you can see what happens as the filter is climbing onto the data.
%\activeplot{wood-hole}{width=6.0in,height=1.75in}{} {
%	Hole in wood.
%	}
%\activeplot{herr-hole}{width=6.0in,height=1.75in}{} {
%	Hole in herringbone.
%	}
%\activeplot{brick-hole}{width=6.0in,height=1.75in}{} {
%	Hole in brick.
%	}



\section{TWO-STAGE LINEAR LEAST SQUARES}
\sx{Two-stage linear least squares}
\sx{two-stage ! least squares}
In Chapter \ref{iin/paper:iin} and
Chapter    \ref{prc/paper:prc}
we filled empty bins
by minimizing the energy output from the filtered mesh.
In each case there was arbitrariness in the choice of the filter.
Here we find and use the optimum filter, the PEF.

\par
%We will develop code in which \bx{missing data}
%in a training dataset
%is restored by a two-stage linear least squares process.
%Since we are using a helix,
%the training volume could be a line, a plane, a volume, or a hypervolume.
%In the first stage,
%we fit a prediction-error filter (PEF) to the given plane.
%Fitting equations that involve empty bins
%are weighted to zero
%so they have no effect on the answer, the PEF.
The first stage is that of the previous section,
finding the optimal PEF while
carefully avoiding using any regression equations
that involve boundaries or missing data.
For the second stage, we take the PEF as known and
find values for the empty bins so that
the power out of the prediction-error filter is minimized.
To do this we find missing data with
module \texttt{mis2()} \vpageref{lst:mis2}.

\par
This two-stage method avoids the nonlinear problem
we would otherwise face if we included the fitting
equations containing both free data values and free filter values.
Presumably, after
two stages of linear least squares
we are close enough to the final solution
that we could switch over to the full nonlinear setup
described near the end of this chapter.

\inputdir{hole}
\par
The synthetic data in Figure \ref{fig:hole}
is a superposition of two plane waves of different directions,
each with a random (but low-passed) waveform.
After punching a hole in the data,
we find that the lost data is pleasingly restored,
though a bit weak near the side boundary.
This imperfection could result
from the side-boundary behavior of the operator
or from an insufficient number of missing-data iterations.
\plot{hole}{width=6in,height=2.5in}{
  Original data (left), with a zeroed hole, restored,
  residual selector (right).
}
\par
The residual selector in Figure \ref{fig:hole}
shows where the filter output has valid inputs.
From it you can deduce the size and shape of the filter,
namely that it matches up with Figure \ref{fig:rabdomain}.
The ellipsoidal hole in the residual selector is larger
than that in the data because we lose regression equations
not only at the hole,
but where any part of the filter overlaps the hole.

\par
The results in Figure \ref{fig:hole} are essentially perfect
representing the fact that that synthetic example
fits the conceptual model perfectly.
Before we look at the many examples
in Figures
\ref{fig:herr-hole-fillr}-\ref{fig:WGstack-hole-fillr}
we will examine another gap-filling strategy.




\subsection{Adding noise (Geostat)}
\inputdir{morgan}
\sx{geostatistics}
In chapter \ref{iin/paper:iin} we restored missing data
by adopting the philosopy of minimizing the energy in filtered output.
In this chapter we learned about an optimum filter
for this task, the prediction-error filter (PEF).
Let us name this method the ``minimum noise'' method
of finding missing data.

\par
A practical problem with the minimum-noise method is evident
in a large empty hole such as 
in Figures \ref{fig:herr-hole-fillr}-
\ref{fig:brick-hole-fillr}.
In such a void the interpolated data diminishes greatly.
Thus we have not totally succeeded in the goal of
``hiding our data acquisition footprint''
which we would like to do if we are trying to make
pictures of the earth and not pictures of our
data acquisition footprint.

\par
What we will do next is useful in some applications but not in others.
Misunderstood or misused it is rightly controversial.
We are going to fill the empty holes
with something that looks like the original data but really isn't.
I will distinguish the words ``\bx{synthetic data}''
(that derived from a physical model)
from ``\bx{simulated data}'' (that manufactured from a statistical model).
We will fill the empty holes with simulated data
like what you see in the center panels of Figures
\ref{fig:granite}-\ref{fig:WGstack}.
We will add just enough of that ``wall paper noise'' to keep
the variance constant as we move into the void.

\par
Given some data $\bold d$, we use it in a filter operator $\bold D$,
and as described with equation (\ref{eqn:exmiss}) we build
a weighting function $\bold W$ that throws out the
broken regression equations (ones that involve missing inputs).
Then we find a PEF $\bold a$ by using this regression.
\begin{equation}
	\bold 0
\quad\approx\quad
	\bold r
\quad = \quad
	\bold W \bold D \bold a
\end{equation}
Because of the way we defined $\bold W$,
the ``broken'' components of $\bold r$ vanish.
We need to know the variance $\sigma$ of the nonzero terms.
It can be expressed mathematically in a couple different ways.
Let $\bold 1$ be a vector filled with ones and let
$\bold r^2$ be a vector containing the squares of the components of $\bold r$.
\begin{equation}
\sigma
\quad = \quad
	\sqrt{ {1\over N} \sum_i^N \,r_i^2 }
\quad = \quad
	\sqrt{
		\bold 1' \bold W \bold r^2
	\over
		\bold 1' \bold W \bold 1 
	}
\end{equation}
Let us go to a random number generator
and get a noise vector $\bold n$
filled with random numbers of variance $\sigma$.
We'll call this the ``added random noise''.
Now we solve this new regression for
the data space $\bold d$ (both known and missing)
\begin{equation}
\bold 0
\quad\approx\quad
\bold r
\quad = \quad
\bold A \bold d   \ - \  \bold n
\label{eqn:noiseeverywhere}
\end{equation}
\par\noindent
keeping in mind that known data is constrained
(as detailed in chapter \ref{iin/paper:iin}).

\par
To understand why this works,
consider first the training image, a region of known data.
Although we might think that the data defines the white noise
residual by $\bold r=\bold A\bold d$, we can also imagine that the white noise
determines the data by $\bold d=\bold A^{-1}\bold r$.
Then consider a region of wholly missing data.  This data
is determined by $\bold d=\bold A^{-1}\bold n$.
Since we want the data variance to be the same in known and unknown
locations, naturally we require the variance of $\bold n$
to match that of $\bold r$.

\par
A very minor issue remains.
Regression equations may have all of their required input data,
some of it, or none of it.
Should the $\bold n$ vector add noise to every regression equation?
First, if a regression equation has all its input data
that means there are no free variables so it doesn't matter
if we add noise to that regression equation because the constraints
will overcome that noise.
I don't know if I should worry about how
{\it many}
inputs are missing for each regression equation.

%We could choose to add noise only to those regression equations
%that have all their inputs.
%Then we'd be fitting
%$\bold 0 \approx\ \bold r = \bold A \bold d    - (\bold I -\bold W)  \bold n$.
%I feel, however, that equation
%(\ref{eqn:noiseeverywhere}) would be slightly better.
%The difference would be most apparent when the PEF was allocated
%more space than it required.
%Then we might see a halo of dimness just inside the rim of the void.

\par
It is fun making all this interesting ``wall paper''
noticing where it is successful and where it isn't.
We cannot help but notice that it seems to work better with
the genuine geophysical data than
it does with many of the highly structured patterns.
Geophysical data is expensive to acquire.
Regrettably, we have uncovered a technolo