paper.tex

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

it simply jumped from a well value at one point
to a seismic value at a neighboring point.
The map would contain discontinuities around each well.
Our philosophy of finding an earth model $\bold m$
is that our earth map should contain no obvious 
``footprint'' of the data acquisition (well locations).
We adopt the philosopy that the difference
between the final map (extended wells)
and the seismic information $\bold x=\bold m-\bold s$ should be smooth.
Thus,
we seek the minimum residual $\bold r$
which is the roughened difference between the seismic data $\bold s$
and the map $\bold m$ of hypothetical omnipresent wells.
With roughening operator $\bold A$ we fit
\begin{equation}
\bold 0\quad\approx\quad \bold r \eq \bold A ( \bold m - \bold s )
	\eq \bold A \bold x
\label{eqn:unconmap}
\end{equation}
along with the constraint
that the map should match the wells at the wells.
We could write this as
$\bold 0 = (\bold I-\bold J) ( \bold m - \bold w )$.
We honor this constraint by initializing the map $\bold m = \bold w$
to the wells (where we have wells, and zero elsewhere).
After we find the gradient direction to suggest some changes
to $\bold m$, we simply will not allow those changes at well locations.
We do this with a mask.
We apply a "missing data selector" to the gradient.
It zeros out possible changes at well locations.
Like with the goal (\ref{eqn:misfitgoal2}),
we have
\begin{equation}
\bold 0\quad\approx\quad \bold r \eq
\bold A \bold J \bold x + \bold A \bold x_{\rm known}
\end{equation}
After minimizing $\bold r$ by adjusting $\bold x$,
we have our solution $ \bold m =  \bold x + \bold s $.

\par
Now we prepare some roughening operators $\bold A$.
We have already coded a 2-D gradient operator
\texttt{igrad2} \vpageref{lst:igrad2}.
Let us combine it with its adjoint to get the 2-D laplacian operator.
(You might notice that the laplacian operator is ``self-adjoint'' meaning
that the operator does the same calculation that its adjoint does.
Any operator of the form $\bold A'\bold A$ is self-adjoint because
$(\bold A'\bold A)'=\bold A'\bold A''=\bold A'\bold A$. )
\par
\opdex{laplac2}{Laplacian in 2-D}{34}{78}{user/fomels} 
Subroutine \texttt{lapfill()}
\vpageref{lst:lapfill} is the same idea as \texttt{mis1()}
\vpageref{lst:mis1} except that 
the filter $\bold A$ has been specialized to the 
laplacian
implemented by module \texttt{laplac2} \vpageref{lst:laplac2}. 


\moddex{lapfill}{Find 2-D missing data}{63}{76}{user/fomels}

\par
Subroutine \texttt{lapfill()}
can be used for each of our two problems,
(1) extending the seismic data to fill space, and
(2) fitting the map exactly to the wells and approximately to the seismic data.
When extending the seismic data,
the initially non-zero components $\bold s \ne \bold 0$ are fixed
and cannot be changed.
\begin{comment}
That is done by calling
\texttt{lapfill()} with \texttt{mfixed=(s/=0.)}.
When extending wells,
the initially non-zero components $\bold w \ne \bold 0$ are fixed
and cannot be changed.
That is done by calling
\texttt{lapfill()} with \texttt{mfixed=(w/=0.)}.
\end{comment}

\par
The final map is shown in Figure \ref{fig:finalmap}.
\plot{finalmap}{width=4.5in,height=3in}{
  Final map based on Laplacian roughening.
}

\par
Results can be computed with various filters.
I tried both $\nabla^2$ and $\nabla$.
There are disadvantages of each,
$\nabla$ being too cautious and
$\nabla^2$ perhaps being too aggressive.
Figure \ref{fig:diffdiff} shows the difference $\bold x$ between
the extended seismic data and the extended wells.
Notice that for $\nabla$ the difference shows
a localized ``tent pole'' disturbance about each well.
For $\nabla^2$ there could be large overshoot between wells,
especially if two nearby wells have significantly different values.
I don't see that problem here.
\par
My overall opinion is that the Laplacian does the better job in this case.
I have that opinion because in viewing the extended gradient
I can clearly see where the wells are.
The wells are where we have acquired data.
We'd like our map of the world to not show where we acquired data.
Perhaps our estimated map of the world cannot help but show where
we have and have not acquired data, but we'd like to minimize that aspect.
\par
\boxit{ A good image of the earth hides our data \bx{acquisition footprint}.  }


\plot{diffdiff}{width=5in,height=2.25in}{
  Difference between wells (the final map)
  and the extended seismic data.
  Left is plotted at the wells (with gray background for zero).
  Center is based on gradient roughening and shows
  tent-pole-like residuals at wells.
  Right is based on Laplacian roughening.
}

\par
To understand the behavior theoretically,
recall that in one dimension 
the filter $\nabla$ interpolates with straight lines
and $\nabla^2$ interpolates with cubics.
This is because the fitting goal
$\bold 0 \approx \nabla \bold m$,
leads to
${\partial \over\partial \bold m'} \bold m'\nabla'\nabla \bold m = \bold 0$
or $\nabla'\nabla \bold m = \bold 0$, whereas the fitting goal 
        $\bold 0 \approx \nabla^2 \bold m$
leads to
        $\nabla^4 \bold m = \bold 0$
which is satisfied by cubics.
In two dimensions, minimizing the output of $\nabla$
gives us solutions of Laplace's equation with sources at the known data.
It is as if $\nabla$ stretches a rubber sheet over poles at each well,
whereas $\nabla^2$ bends a stiff plate.
\par
Just because $\nabla^2$ gives smoother maps than  $\nabla$
does not mean those maps are closer to reality.
This is a deeper topic, addressed in Chapter~\ref{mda/paper:mda}.
It is the same issue we noticed when comparing
figures \ref{fig:mlines}-\ref{fig:mbest}.

\section{SEARCHING THE SEA OF GALILEE}
\sx{Sea of Galilee}
\sx{Galilee, Sea of}
\inputdir{galilee}
Figure \ref{fig:locfil} shows a bottom-sounding survey
of the Sea of Galilee\footnote{
        Data collected by Zvi \bx{ben Avraham}, TelAviv University.
        Please communicate with him {\tt [e-mail]}
        for more details or if you make something
        publishable with his data.
        }
at various stages of processing.
The ultimate goal is not only a good map of
the depth to bottom,
but images useful for the purpose
of identifying \bx{archaeological}, geological, or
geophysical details of the sea bottom.
The Sea of Galilee is unique
because it is a {\it fresh}-water lake {\it below} sea-level.
It seems to be connected to the great rift (pull-apart)
valley crossing east Africa.
We might delineate the Jordan River delta.
We might find springs on the water bottom.
We might find archaeological objects.
\par
\plot{locfil}{width=6in,height=8.5in}{
  Views of the bottom of the Sea of Galilee.
}
The raw data is 132,044 triples, $(x_i,y_i,z_i)$,
where $x_i$ ranges over about 12 km and
where $y_i$ ranges over about 20 km.
The lines you see in Figure \ref{fig:locfil}
are sequences of data points, i.e., the track of the survey vessel.
The depths $z_i$ are recorded to an accuracy of about 10 cm.
\par
The first frame in Figure~\ref{fig:locfil} shows simple binning.
A coarser mesh would avoid the empty bins but lose resolution.
As we refine the mesh for more detail,
the number of empty bins grows
as does the care needed in devising a technique
for filling them.
This first frame uses the simple idea from Chapter \ref{ajt/paper:ajt} of
spraying all the data values to the nearest bin
with \texttt{bin2()} \vpageref{lst:bin2}
and dividing by the number in the bin.
Bins with no data obviously need to be filled in some other way.
I used a missing data program like that in the recent section
on ``wells not matching the seismic map.''
Instead of roughening with a Laplacian, however,
I used the gradient operator \texttt{igrad2} \vpageref{lst:igrad2}
% I used the \bx{roughening operator} \texttt{rufftri2()} \vpageref{lst:rufftri2},
% which is an impulse minus a broad smoothing operator.
The solver is \texttt{grad2fill()}.
\moddex{grad2fill}{low cut missing data}{52}{60}{user/fomels}

%\par
%Having a roughener with an adjustable width in the smoothing operator
%allows us to choose a cutoff frequency.
%This is interesting for two reasons:
%(1) The output of the roughening operator
%is interesting to look at; and
%(2) Larger filters allow missing data iterations to converge more quickly.

\par
The output of the roughening operator is an image,
a filtered version of the depth,
a filtered version of something real.
Such filtering can enhance the appearance of interesting features.
For example,
scanning the shoreline of the roughened image
(after missing data was filled),
we see several ancient shorelines, now submerged.
\par
\boxit{The adjoint is the easiest image to build.
The roughened map is often more informative than the map itself.}
\par
The views expose several defects
of the data acquisition and of our data processing.
The impulsive glitches (St.~Peter's fish?)
need to be removed but we must be careful not to throw
out the sunken ships along with the bad data points.
Even our best image shows clear evidence of the recording vessel's tracks.
Strangely, some tracks are deeper than others.
Perhaps the survey is assembled from work done in different seasons
and the water level varied by season.
Perhaps some days the vessel was more heavily loaded
and the depth sounder was on a deeper keel.
As for the navigation equipment,
we can see that some data values are reported outside the lake!
\par
%\boxit{ A good image of the earth hides our data \bx{acquisition footprint}.  }
\par
We want the sharpest possible view
of this classical site.
A treasure hunt is never easy
and no one guarantees we will
find anything of great value
but at least the exercise is a good warm-up
for submarine petroleum exploration.


\section{INVERSE LINEAR INTERPOLATION}
\inputdir{invint}
In Chapter \ref{ajt/paper:ajt} we defined \bx{linear interpolation}
\sx{interpolation}
as the extraction of values from between mesh points.
In a typical setup (occasionally the role of data and model are swapped),
a model is given on a uniform mesh
and we solve the easy problem of extracting values
between the mesh points with subroutine \texttt{lint1()} \vpageref{lst:lint1}.
The genuine problem is the inverse problem, which we attack here.
Data values are sprinkled all around,
and we wish to find a function on a uniform mesh
from which we can extract that data by \bx{linear interpolation}.
The adjoint operator for subroutine {\tt lint1()}
simply piles data back into its proper location in model space
without regard to how many data values land in each region.
Thus some model values may have many data points added
to them while other model values get none.
We could interpolate by minimizing the energy in the model gradient,
or that in the second derivative of the model,
or that in the output of any other roughening filter
applied to the model.
\par
Formalizing now our wish
that data $\bold d$ be extractable by \bx{linear interpolation} $\bold F$,
from a model $\bold m$,
and our wish that application of a roughening filter
with an operator $\bold A$ have minimum energy, we write the fitting goals:
\begin{equation}
        \begin{array}{lll}
        \bold 0 &\approx & \bold F \bold m - \bold d \\
        \bold 0 &\approx & \bold A \bold m
        \label{eqn:invintwish}
        \end{array}
\end{equation}
Suppose we take the roughening filter to be the second difference operator
$(1,-2,1)$
scaled by a constant $\epsilon$,
and suppose we have a data point near each end of the model
and a third data point exactly in the middle.
Then,
for a model space 6 points long,
the fitting goal could look like
\def\E{\epsilon}
\begin{equation} { 
\left[ 
\begin{array}{rrrrrr}
   .8 & .2 & .  & .  & .  & .  \\
   .  & .  & 1  & .  & .  & .  \\
   .  & .  & .  & .  & .5 & .5 \\
   \hline
   \E & .  & .  & .  & .  & .  \\
  -2\E& \E & .  & .  & .  & .  \\
   \E &-2\E&  \E& .  & .  & .  \\
   .  & \E &-2\E&  \E& .  & .  \\
   .  & .  &  \E&-2\E&  \E& .  \\
   .  & .  & .  &  \E&-2\E&  \E\\
   .  & .  & .  & .  &  \E&-2\E\\
   .  & .  & .  & .  & .  &  \E
  \end{array} \right] 
\left[ 
        \begin{array}{c}
          m_0 \\ 
          m_1 \\ 
          m_2 \\ 
          m_3 \\ 
          m_4 \\ 
          m_5 
        \end{array}
\right] 
\ -\ 
\left[ 
\begin{array}{c}
  d_0 \\ 
  d_1 \\ 
  d_2 \\ 
  \hline
  0   \\
  0   \\
  0   \\
  0   \\
  0   \\
  0   \\
  0   \\
  0
  \end{array} \right] 
\eq
\left[ 
\begin{array}{c}
  \bold r_d \\ 
  \bold r_m 
  \end{array} \right] 
\quad \approx \ \bold 0
\label{eqn:tworegexam}
} \end{equation}
\par
The residual vector has two parts,
a data part $\bold r_d$ on top
and a model part $\bold r_m$ on the bottom.
The data residual
should vanish except where contradictory data values
happen to lie in the same place.
The model residual is the roughened model.

\begin{comment}
\par
Two fitting goals (\ref{eqn:invintwish}) are so common in practice
that it is convenient to adopt our least-square fitting
subroutine \texttt{solver_smp} \vpageref{lst:solver_smp} accordingly.
The modification
is shown in module \texttt{solver_reg} \vpageref{lst:solver_reg}.
In addition to specifying the ``data fitting'' operator $\bold F$
(parameter \texttt{Fop}),
we need to pass the ``model regularization'' operator $\bold A$
(parameter \texttt{Aop}) and the
size of its output (parameter \texttt{nAop}) for proper memory allocation.

%\begin{notforlecture}
\par
(When I first looked at module \texttt{solver_reg} I was appalled
by the many lines of code, especially all the declarations.
Then I realized how much much worse was Fortran 77 where
I needed to write a new solver for every pair of operators.
This one solver module works for all operator pairs
and for many optimization descent strategies
because these ``objects'' are arguments.
These more powerful objects require declarations that are more complicated
than the simple objects of Fortran 77.
As an author I have a dilemma:
To make algorithms compact (and seem simple) requires many careful definitions.
When these definitions put in the code, they are careful,
but the code becomes annoyingly verbose.
Otherwise, the definitions must go in the surrounding natural language
where they are not easily made precise.)