paper.tex

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

%   \end{tabbing}
%}
The reverse conversion from the floating-point variable
to the counting variable is inexact.
The easiest thing is to place it at the nearest neighbor.
This is done by solving for {\tt iz}, then adding one half,
and then rounding down to the nearest integer.
The familiar computer idioms are:\begin{verbatim}
        iz =  .5 + 1 + ( z - z0) / dz
        iz = 1.5 +     ( z - z0) / dz
        i3 = 1.5 +     (x3 - o3) / d3
\end{verbatim}
A small warning is in order:
People generally use positive counting variables.
If you also include negative ones,
then to get the nearest integer,
you should do your rounding with the Fortran function {\tt NINT()}.

\subsection{Data-push binning}
\sx{nearest neighbor binning}
\sx{data-push binning}
Binning is putting data values in bins.
Nearest-neighbor binning is an operator.
There is both a forward operator and its adjoint.
Normally the model consists of values given on a uniform mesh,
and the data consists of pairs of numbers (ordinates at coordinates)
sprinkled around in the continuum
(although sometimes the data is uniformly spaced and the model is not).
\par
In both the forward and the adjoint operation,
each data coordinate is examined
and the nearest mesh point (the bin) is found.
For the forward operator,
the value of the bin is added to that of the data.
The adjoint is the reverse:
we add the value of the data to that of the bin.
Both are shown in two dimensions in subroutine \texttt{bin2()}.
\opdex{bin2}{push data into bin}{46}{55}{user/gee}
The most typical application requires an additional step, inversion.
In the inversion applications
each bin contains a different number of data values.
After the adjoint operation is performed,
the inverse operator divides the bin value
by the number of points in the bin.
It is this inversion operator that is generally called binning.
To find the number of data points in a bin,
we can simply apply the adjoint of \texttt{bin2()} to pseudo data of all ones.


\subsection{Linear interpolation}
\inputdir{XFig}
\par
The \bx{linear interpolation}
operator is much like the binning operator but a little fancier.
When we perform the forward operation, we take each data coordinate
and see which two model mesh points bracket it.
Then we pick up the two bracketing model values
and weight each of them
in proportion to their nearness to the data coordinate,
and add them to get the data value (ordinate).
The adjoint operation is adding a data value
back into the model vector;
using the same two weights,
this operation distributes the ordinate value
between the two nearest points in the model vector.
For example, 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,
as shown in Figure \ref{fig:helgerud},
we have the operator in (\ref{eqn:lintseq}).
\sideplot{helgerud}{width=3in}{
  Uniformly sampled model space
  and irregularly sampled data space corresponding
  to \protect(\ref{eqn:lintseq}).
}
\begin{equation}
\left[ 
\begin{array}{c}
  d_0 \\ 
  d_1 \\ 
  d_2 
  \end{array} \right] 
\quad \approx \quad
\left[ 
\begin{array}{rrrrrr}
   .8 & .2 &  .  & .  & .  & .  \\
   .  & .  &  1  & .  & .  & .  \\
   .  & .  &  .  & .  & .5 & .5 
  \end{array} \right] 
\left[ 
	\begin{array}{c}
	  m_0 \\ 
	  m_1 \\ 
	  m_2 \\ 
	  m_3 \\ 
	  m_4 \\ 
	  m_5
	\end{array}
\right] 
\label{eqn:lintseq}
\end{equation}
The two weights in each row sum to unity.
If a binning operator were used for the same data and model,
the binning operator would contain a ``1.'' in each row.
In one dimension (as here),
data coordinates are often sorted into sequence,
so that the matrix is crudely a diagonal matrix like equation (\ref{eqn:lintseq}).
If the data coordinates covered the model space uniformly,
the adjoint would roughly be the inverse.
Otherwise,
when data values pile up in some places and gaps remain elsewhere,
the adjoint would be far from the inverse.
\par
Subroutine \texttt{lint1()} does linear interpolation and its adjoint.
\opdex{lint1}{linear interp}{45}{59}{user/gee}


\subsection{Causal integration}
\sx{causal integration}
Causal integration is defined as
\begin{equation}
y(t) \eq \int_{-\infty}^t \ x(t)\ dt
\end{equation}
Sampling the time axis gives a matrix equation which
we should call causal summation, but we often call it causal integration.
\begin{equation}
  \left[
	\begin{array}{c}
		y_0 \\
		y_1 \\
		y_2 \\
		y_3 \\
		y_4 \\
		y_5 \\
		y_6 \\
		y_7 \\
		y_8 \\
		y_9 \\
	\end{array}
  \right]
 \quad = \quad
  \left[
	\begin{array}{ccccccccccc}
	1 &0 &0 &0 &0 &0 &0 &0 &0 &0 \\
	1 &1 &0 &0 &0 &0 &0 &0 &0 &0 \\
	1 &1 &1 &0 &0 &0 &0 &0 &0 &0 \\
	1 &1 &1 &1 &0 &0 &0 &0 &0 &0 \\
	1 &1 &1 &1 &1 &0 &0 &0 &0 &0 \\
	1 &1 &1 &1 &1 &1 &0 &0 &0 &0 \\
	1 &1 &1 &1 &1 &1 &1 &0 &0 &0 \\
	1 &1 &1 &1 &1 &1 &1 &1 &0 &0 \\
	1 &1 &1 &1 &1 &1 &1 &1 &1 &0 \\
	1 &1 &1 &1 &1 &1 &1 &1 &1 &1 
	\end{array}
  \right]
  \ \ 
  \left[
	\begin{array}{c}
		x_0 \\
		x_1 \\
		x_2 \\
		x_3 \\
		x_4 \\
		x_5 \\
		x_6 \\
		x_7 \\
		x_8 \\
		x_9 \\
	\end{array}
  \right]
\label{eqn:mytri}
\end{equation}
%Having 1/2 instead of 1 on the diagonal clutters our code a little
%but it is a worthwhile improvement. 
%With 1/2, the summation equivalent of
%$\int_{-\infty}^t + \int_t^{\infty} = \int_{-\infty}^\infty$
%remains exactly true.
%Simplifying the 1/2 to a 1 is also accurate if
%we agree to having the integral evaluated half way between two mesh points.
%For most applications, however, it is more convenient to have
%the integral represented on the same mesh as the function itself,
%so we need the 1/2.
(In some applications the 1 on the diagonal is replaced by 1/2.)
Causal integration is
the simplest prototype of a recursive operator.
\sx{recursion, integration}
The coding is trickier than operators we considered earlier.
Notice when you compute $y_5$ that it is the sum of 6 terms,
but that this sum is more quickly computed as $y_5 = y_4 + x_5$.
Thus equation~(\ref{eqn:mytri}) is more efficiently thought of as
the recursion
\begin{equation}
y_t \quad = \quad y_{t-1} + x_t
\quad
\quad
\quad {\rm for\ increasing\ } t
\label{eqn:myrecur}
\end{equation}
(which may also be regarded as a numerical representation
of the \bx{differential equation} $dy/dt=x$.)
\par
When it comes time to think about the adjoint, however,
it is easier to think of equation~(\ref{eqn:mytri}) than of~(\ref{eqn:myrecur}).
Let the matrix of equation~(\ref{eqn:mytri}) be called $\bold C$.
Transposing to get $\bold C '$ and applying it to $\bold y$
gives us something back in the space of $\bold x$,
namely $\tilde{\bold x} = \bold C' \bold y$.
From it we see that the adjoint calculation,
if done recursively,
needs to be done backwards like
\begin{equation}
\tilde x_{t-1} \quad = \quad \tilde x_{t} + y_{t-1}
\quad
\quad
\quad {\rm for\ decreasing\ } t
\label{eqn:backrecur}
\end{equation}
We can sum up by saying that the adjoint of causal integration
is anticausal integration.
\par
A subroutine to do these jobs is \texttt{causint()} \vpageref{/prog:causint}.
The code for anticausal integration is not obvious
from the code for integration and the adjoint coding tricks we
learned earlier.
To understand the adjoint, you need to inspect
the detailed form of the expression $\tilde{\bold x} = \bold C' \bold y$
and take care to get the ends correct.
\opdex{causint}{causal integral}{35}{46}{filt/lib}
\par
%\activesideplot{causint}{width=3.00in,height=3.2in}{ER}{
\inputdir{causint}
\sideplot{causint}{width=3.00in}{
	{\tt in1} is an input pulse.
	{\tt C in1} is its causal integral.
	{\tt C' in1} is the anticausal integral of the pulse.
	{\tt in2} is a separated doublet.
	Its causal integration is a box and its anticausal integration
	is the negative.
	{\tt CC in2} is the double causal integral of {\tt in2}.
	How can an equilateral triangle be built?
	}

\par
Later we will consider equations
to march wavefields up towards the earth's surface,
a layer at a time, an operator for each layer.
Then the adjoint will start from the earth's surface
and march down, a layer at a time, into the earth.
\begin{exer}
\item
Modify the calculation in Figure~\ref{fig:causint} to make
a triangle waveform on the bottom row.
%\item
%Notice that the triangle waveform is not time aligned
%with the input {\tt in2}.
%Force time alignment with the operator 
%${\bold C' \bold C}$ or
%${\bold C  \bold C'}$.
%\item
%Modify \GPROG{causint} by changing the diagonal to contain
%1/2 instead of 1.
%Notice how time alignment changes in Figure~\FIG{causint}.
\end{exer}


\section{ADJOINTS AND INVERSES}

Consider a model $\bold m$ and an operator $\bold F$ which creates some
theoretical data $\bold d_{\rm theor}$.
\begin{equation}
	\bold d_{\rm theor} \eq \bold F \bold m
\end{equation}
The general task of geophysicists is to begin from
observed data $\bold d_{\rm obs}$ and
find an estimated model $\bold m_{\rm est}$
that satisfies the simultaneous equations
\begin{equation}
	\bold d_{\rm obs} \eq \bold F \bold m_{\rm est}
\end{equation}
This is the topic of a large discipline variously called
``inversion'' or ``estimation''.
Basically, it defines a residual
$\bold r = \bold d_{\rm obs}-\bold d_{\rm theor}$
and then minimizes its length $\bold r \cdot \bold r$.
Finding $\bold m_{\rm est}$ this way is called
the \bx{least squares} method.
The basic result (not proven here) is that
\begin{equation}
\bold m_{\rm est} = (\bold F'\bold F)^{-1}\bold F'\bold d_{\rm obs}
\end{equation}
In many cases including all seismic imaging cases,
the matrix
$\bold F'\bold F$
is far too large to be invertible.
People generally proceed by a rough guess at an approximation
for $(\bold F'\bold F)^{-1}$.
The usual first approximation is
the optimistic one that $(\bold F'\bold F)^{-1}=\bold I$.
To this happy approximation, the inverse $\bold F^{-1}$
is the adjoint $\bold F'$.

\par
In this book we'll see examples where
$\bold F'\bold F\approx \bold I$
is a good approximation and other examples where it isn't.
We can tell how good the approximation is.
We take some hypothetical data and convert it to a model,
and use that model to make some reconstructed data
$\bold d_{\rm recon} = \bold F \bold F' \bold d_{\rm hypo}$.
Likewise we could go from a hypothetical model to some data and then
to a reconstructed model
$\bold m_{\rm recon} = \bold F' \bold F \bold m_{\rm hypo}$.
Luckily, it often happens that the reconstructed differs from
the hypothetical in some trivial way,
like by a scaling factor, or by a scaling factor
that is a function of physical location or time,
or a scaling factor that is a function of frequency.
It isn't always simply a matter of a scaling factor,
but it often is, and when it is, we often simply
redefine the operator to include the scaling factor.
Observe that there are two places for scaling functions (or filters),
one in model space, the other in data space.

\par
We could do better than the adjoint
%(better than assuming $\bold F'\bold F=\bold I$)
by iterative modeling methods (conjugate gradients)
that are also described elsewhere.
These methods generally demand that the adjoint be computed correctly.
As a result, we'll be a little careful about adjoints in
this book to compute them correctly
even though this book does not require them to be exactly correct.

\subsection{Dot product test}

\par
We define an adjoint when we write a program that computes one.
In an abstract logical mathematical sense, however,
every adjoint is defined by a \bx{dot product test}.
This abstract definition gives us no clues how to code our program.
After we have finished coding, however, this abstract definition
(which is actually a test) has considerable value to us.


\par
Conceptually, the idea of matrix transposition is simply ${a}_{ij}'=a_{ji}$.
In practice, however, we often encounter matrices far too large
to fit in the memory of any computer.
Sometimes it is also not obvious how to formulate the process at hand
as a matrix multiplication.
(Examples are differential equations and fast Fourier transforms.)
What we find in practice is that an application and its adjoint 
amounts to two subroutines. The first subroutine
amounts to the matrix multiplication $ \bold F \bold x$.
The adjoint subroutine computes $\bold F' \bold y$,
where $\bold F'$ is the conjugate-transpose matrix.
Most methods of solving inverse problems will fail
if the programmer provides an inconsistent pair of subroutines
for $\bold F$ and $\bold F'$.
The dot product test described next
is a simple test for verifying that the two 
subroutines really are adjoint to each other.

\par
The matrix expression
$\bold y' \bold F \bold x $
may be written with parentheses as either
$(\bold y' \bold F) \bold x $ or
$\bold y' (\bold F \bold x)$.
Mathematicians call this the ``associative'' property.
If you write matrix multiplication using summation symbols,
you will notice that putting parentheses around matrices simply
amounts to reordering the sequence of computations.
But we soon get a very useful result.
Programs for some linear operators are far from obvious,
for example \texttt{causint()} \vpageref{/prog:causint}.
Now we build a useful test for it.
%The associative property of linear algebra says that
%we do not need parentheses in a vector-matrix-vector product
%like $\bold y' \bold F \bold x $ because we get the same
%result no matter where we put the parentheses.

\begin{eqnarray}
\bold y' ( \bold F \bold x ) &=& ( \bold y' \bold F )  \bold x   \\
\bold y' ( \bold F \bold x ) &=& ( \bold F' \bold y )' \bold x 
\label{eqn:bilin2}
\end{eqnarray}
For the dot-product test,
load the vectors $\bold x$ and $\bold y$ with random numbers.
Compute the vector $\tilde \bold y = \bold F\bold x$
using your program for $\bold F$,
and compute
$\tilde \bold x = \bold F'\bold y$
using your program for $\bold F'$.
Inserting these into equation~(\ref{eqn:bilin2})
gives you two scalars that should be equal.
\begin{equation}
\bold y' ( \bold F \bold x ) \eq
\bold y' \tilde \bold y \eq \tilde \bold x ' \bold x
\eq ( \bold F' \bold y )' \bold x
\label{eqn:bilin}
\end{equation}
The left and right sides of this equation will be computationally equal
only if the program doing $\bold F'$ is indeed adjoint
to the program doing $\bold F$
(unless the random numbers do something miraculous).
Note that the vectors $\bold x$ and $\bold y$
are generally of different lengths.

\par
Of course passing the dot product test does not prove that a computer code
is correct, but if the test fails we know the code is incorrect.
More information about adjoint operators,
and much more information about inverse operators
is found in my other books,
Earth Soundings Analysis: Processing versus inversion (PVI) and
Geophysical Estimation by Example (GEE).