paper.tex

来自「国外免费地震资料处理软件包」· TEX 代码 · 共 1,940 行 · 第 1/5 页

  y_5 \\ 
  y_6 
  \end{array} \right] 
\eq
\left[ 
\begin{array}{ccc}
  x_3 & x_2 & x_1  \\
  x_4 & x_3 & x_2  \\
  x_5 & x_4 & x_3  \\
  x_6 & x_5 & x_4 
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  b_1 \\ 
  b_2 \\ 
  b_3 \end{array} \right]
\label{eqn:trunccontran2}
\end{equation}
The adjoint is
\begin{equation}
\left[ 
\begin{array}{c}
\hat  b_1 \\ 
\hat  b_2 \\ 
\hat  b_3
  \end{array} \right] 
\eq
\left[ 
\begin{array}{ccccc}
  x_3 & x_4 & x_5 & x_6  \\
  x_2 & x_3 & x_4 & x_5  \\
  x_1 & x_2 & x_3 & x_4 
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6 
  \end{array} \right]
\label{eqn:ajtintcontran2}
\end{equation}
The difference between (\ref{eqn:ajtintcontran2})
and (\ref{eqn:ajtcontran2}) is that here 
the adjoint crosscorrelates a fixed portion of {\em  output}
across a variable portion of {\em  input},
whereas with (\ref{eqn:ajtcontran2})
the adjoint crosscorrelates a fixed portion of {\em  input}
across a variable portion of {\em  output}.
\par
In practice we typically allocate equal space for input and output.
Because the output is shorter than the input,
it could slide around in its allocated space,
so its location is specified by an additional parameter called its {\tt lag}.
%The statement {\tt x=y-b+lag} in module
%\texttt{icaf1}
%says that the
%output time {\tt y} aligns with the
% input time {\tt x} for the filter point {\tt b=lag}.
\opdex{icaf1}{convolve internal}{45}{50}{user/gee}
The value of {\tt lag} always used in this book is {\tt lag=1}.
For {\tt lag=1} the module \texttt{icaf1}
implements not equation (\ref{eqn:trunccontran2})
but (\ref{eqn:padconvin}):
\begin{equation}
\left[ 
\begin{array}{c}
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6 
  \end{array} \right] 
\eq
\left[ 
\begin{array}{ccc}
  0   & 0   & 0    \\
  0   & 0   & 0    \\
  x_3 & x_2 & x_1  \\
  x_4 & x_3 & x_2  \\
  x_5 & x_4 & x_3  \\
  x_6 & x_5 & x_4 
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  b_1 \\ 
  b_2 \\ 
  b_3 \end{array} \right]
\label{eqn:padconvin}
\end{equation}
It may seem a little odd to put the required zeros at the beginning
of the output, but filters are generally designed so that their
strongest coefficient is the first, namely {\tt bb(1)} so
the alignment of input and output in equation (\ref{eqn:padconvin})
is the most common one.

\inputdir{conv}
\par
The
\bx{end effect}s
of the convolution modules are summarized
in Figure~\ref{fig:conv}.%
\sideplot{conv}{width=3in}{
        Example of convolution end-effects.
        From top to bottom:
        input;
        filter;
        output of {\tt tcai1()};
        output of {\tt icaf1()} also with ({\tt lag=1}).
        }
%\end{notforlecture}

\subsection{Zero padding is the transpose of truncation}
Surrounding a dataset by zeros
(\bx{zero pad}ding)
is adjoint to throwing away the extended data
(\bx{truncation}).
Let us see why this is so.
Set a signal in a vector $\bold x$, and
then to make a longer vector $\bold y$,
add some zeros at the end of $\bold x$.
This zero padding can be regarded as the matrix multiplication
\begin{equation}
\bold y\eq
 \left[ 
  \begin{array}{c}
   \bold I \\ 
   \bold 0
  \end{array}
 \right] 
 \ 
 \bold x
\end{equation}
The matrix is simply an identity matrix $\bold I$
above a zero matrix $\bold 0$.
To find the transpose to zero-padding, we now transpose the matrix
and do another matrix multiply:
\begin{equation}
\tilde {\bold x} \eq
 \left[ 
  \begin{array}{cc}
   \bold I & \bold 0
  \end{array}
 \right] 
\ 
\bold y
\end{equation}
So the transpose operation to zero padding data
is simply {\em  truncating} the data back to its original length.
Module \texttt{zpad1} below
pads zeros on both ends of its input.
Modules for two- and three-dimensional padding are in the
library named {\tt zpad2()} and {\tt zpad3()}.
\opdex{zpad1}{zero pad 1-D}{27}{31}{user/gee}


\subsection{Adjoints of products are reverse-ordered products of adjoints}
Here we examine an example of the general idea that
adjoints of products are reverse-ordered products of adjoints.
For this example we use the Fourier transformation.
No details of \bx{Fourier transformation} are given here
and we merely use it as an example of a square matrix $\bold F$.
We denote the complex-conjugate transpose (or \bx{adjoint}) matrix
with a prime,
i.e.,~$\bold F'$.
The adjoint arises naturally whenever we consider energy.
The statement that Fourier transforms conserve energy is 
$\bold y'\bold y=\bold x'\bold x$ where $\bold y= \bold F \bold x$.
Substituting gives $\bold F'\, \bold F = \bold I$, which shows that
the inverse matrix to Fourier transform
happens to be the complex conjugate of the transpose of $\bold F$.
\par
With Fourier transforms,
\bx{zero pad}ding and \bx{truncation} are especially prevalent.
Most modules transform a dataset of length of $2^n$,
whereas dataset lengths are often of length $m \times 100$.
The practical approach is therefore to pad given data with zeros.
Padding followed by Fourier transformation $\bold F$
can be expressed in matrix algebra as
\begin{equation}
{\rm Program} \eq
\bold F \ 
 \left[ 
  \begin{array}{c}
   \bold I \\ 
   \bold 0
  \end{array}
 \right] 
\end{equation}
According to matrix algebra, the transpose of a product,
say $\bold A \bold B = \bold C$,
is the product $\bold C' = \bold B' \bold A'$ in reverse order.
So the adjoint routine is given by
\begin{equation}
{\rm Program'} \eq
 \left[ 
  \begin{array}{cc}
   \bold I & \bold 0
  \end{array}
 \right] 
\
\bold F'
\end{equation}
Thus the adjoint routine
{\em  truncates} the data {\em  after} the inverse Fourier transform.
This concrete example illustrates that common sense often represents
the mathematical abstraction
that adjoints of products are reverse-ordered products of adjoints.
It is also nice to see a formal mathematical notation
for a practical necessity.
Making an approximation need not lead to collapse of all precise analysis.


\subsection{Nearest-neighbor coordinates}
\sx{nearest neighbor coordinates}
In describing physical processes,
we often either specify models as values given on a uniform mesh
or we record data on a uniform mesh.
Typically we have
a function $f$ of time $t$ or depth $z$
and we represent it by {\tt f(iz)}
corresponding to $f(z_i)$ for $i=1,2,3,\ldots, n_z$
where $z_i = z_0+ (i-1)\Delta z$.
We sometimes need to handle depth as
an integer counting variable $i$
and we sometimes need to handle it as
a floating-point variable $z$.
Conversion from the counting variable to the floating-point variable
is exact and is often seen in a computer idiom
such as either of \begin{verbatim}
            do iz= 1, nz {   z = z0 + (iz-1) * dz 
            do i3= 1, n3 {  x3 = o3 + (i3-1) * d3
\end{verbatim}
%{\tt
%   \begin{tabbing}  indent \= \kill
%       \> do iz= 1, nz \{   z = z0 + (iz-1) * dz  \\
%       \> do i3= 1, n3 \{  x3 = o3 + (i3-1) * d3
%   \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 = 0.5 + ( z - z0) / dz;
        i3 = 0.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 C function {\tt floor()}.

\subsection{Data-push binning}
\sx{nearest neighbor binning}
\sx{data-push binning}
\inputdir{galilee}
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 module \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.
To capture this idea in an equation,
let $\bold B$ denote the linear operator
in which the bin value is sprayed to the data values.
The inverse operation,
in which the data values in the bin are summed
and divided by the number in the bin, is represented by
\begin{equation}
\bold m \eq {\bf diag}( \bold B' \bold 1)^{-1} \bold B' \bold d
\label{eqn:dpbin}
\end{equation}
Empty bins, of course, leave us a problem.
That we'll address in chapter \ref{iin/paper:iin}.
In Figure \ref{fig:galbin}, the empty bins contain zero values.
\plot{galbin}{width=6in,height=4.0in}{
  Binned depths of the Sea of Galilee.
}

\subsection{Linear interpolation}
\par
\inputdir{XFig}
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 bin centers 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,
the adjoint distributes the data ordinate value
between the two nearest bins 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}
   .7 & .3 &  .  & .  & .  & .  \\
   .  & .  &  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
Module \texttt{lint1} does linear interpolation and its adjoint.
In chapters \ref{iin/paper:iin} and \ref{mda/paper:mda}
we build inverse operators.
\opdex{lint1}{linear interp}{45}{59}{user/gee}

%\begin{notforlecture}
\subsection{Spray and sum : scatter and gather}
\sx{spray}
\sx{scatter}
\sx{gather}
Perhaps the most common operation is the summing of many values
to get one value.
Its adjoint operation takes a single input value
and throws it out to a space of many values.
The \bx{summation operator} is a row vector of ones.
Its adjoint is a column vector of ones.
In one dimension
this operator is almost too easy for us to bother showing a routine.
But it is more interesting in three dimensions,
where we could be summing or spraying on any of three subscripts,
or even summing on some and spraying on others.
In module \texttt{spraysum},