paper.tex

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

both input and output are taken
to be three-dimensional arrays.
Externally, however, either could be a scalar, vector, plane, or cube.
For example,
the internal array  {\tt xx(n1,1,n3)}
could be externally the matrix  {\tt map(n1,n3)}.
When
module \texttt{spraysum} is given
the input dimensions and output dimensions stated below,
the operations stated alongside are implied.

\begin{tabular}{p{1em}lll}
&{\tt (n1,n2,n3)}  &{\tt (1,1,1)}     &Sum a cube into a value.      \\
&{\tt (1,1,1)}     &{\tt (n1,n2,n3)}  &Spray a value into a cube.\\
&{\tt (n1,1,1)}    &{\tt (n1,n2,1)}   &Spray a column into a matrix.\\
&{\tt (1,n2,1)}    &{\tt (n1,n2,1)}   &Spray a row into a matrix.\\
&{\tt (n1,n2,1)}   &{\tt (n1,n2,n3)}  &Spray a plane into a cube.\\
&{\tt (n1,n2,1)}   &{\tt (n1,1,1)}    &Sum rows of a matrix into a column.\\
&{\tt (n1,n2,1)}   &{\tt (1,n2,1)}    &Sum columns of a matrix into a row.\\
&{\tt (n1,n2,n3)}  &{\tt (n1,n2,n3)}  &Copy and add the whole cube.
\end{tabular}

If an axis is not of unit length on either input or output,
then both lengths must be the same; otherwise, there is an error.
Normally, after (possibly) erasing the output,
we simply loop over all points on each axis, adding the input to the output.
This implements either a copy or an add, depending on the {\tt add} parameter.
It is either a spray, a sum, or a copy,
according to the specified axis lengths.
\opdex{spraysum}{sum and spray}{42}{56}{user/gee}

%\end{notforlecture}

\subsection{Causal and leaky integration}
\sx{causal integration}
\sx{leaky integration}
\sx{integration ! leaky}
\sx{integration ! causal}
\inputdir{causint}
Causal integration is defined as
\begin{equation}
y(t) \eq \int_{-\infty}^t \ x(\tau )\ d\tau 
\end{equation}
Leaky integration is defined as
\begin{equation}
y(t) \eq \int_0^\infty \ x(t-\tau)\ e^{-\alpha\tau} \ d\tau
\end{equation}
As $\alpha \rightarrow 0$, leaky integration becomes causal integration.
The word ``leaky'' comes from electrical circuit theory
where the voltage on a capacitor would be the integral of the current
if the capacitor did not leak electrons.

\par
Sampling the time axis gives a matrix equation that
we should call causal summation,
but we often call it causal integration.
Equation (\ref{eqn:leakinteq}) represents causal integration for $\rho=1$
and leaky integration for $0<\rho <1$.

\begin{equation}
\bold y \eq
  \left[
        \begin{array}{c}
                y_0 \\
                y_1 \\
                y_2 \\
                y_3 \\
                y_4 \\
                y_5 \\
                y_6 
        \end{array}
  \right]
 \quad = \quad
  \left[
        \begin{array}{cccccccc}
             1 &     0 &     0 &     0 &     0 &     0 &     0 \\
        \rho   &     1 &     0 &     0 &     0 &     0 &     0 \\
        \rho^2 &\rho   &     1 &     0 &     0 &     0 &     0 \\
        \rho^3 &\rho^2 &\rho   &     1 &     0 &     0 &     0 \\
        \rho^4 &\rho^3 &\rho^2 &\rho   &     1 &     0 &     0 \\
        \rho^5 &\rho^4 &\rho^3 &\rho^2 &\rho   &     1 &     0 \\
        \rho^6 &\rho^5 &\rho^4 &\rho^3 &\rho^2 &\rho   &     1 
        \end{array}
  \right]
  \ \ 
  \left[
        \begin{array}{c}
                x_0 \\
                x_1 \\
                x_2 \\
                x_3 \\
                x_4 \\
                x_5 \\
                x_6 
        \end{array}
  \right]
\eq \bold C \bold x
\label{eqn:leakinteq}
\end{equation}
(The discrete world is related to the continuous by
$ \rho = e^{-\alpha\Delta\tau}$ and in
some applications, the diagonal is 1/2 instead of 1.)
Causal integration is
the simplest prototype of a recursive operator.
\sx{recursion ! integration}
The coding is trickier than
that for the 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 = \rho y_4 + x_5$.
Thus equation~(\ref{eqn:leakinteq}) is more efficiently thought of as
the recursion
\begin{equation}
y_t \eq \rho\; y_{t-1} + x_t
\quad
\quad
\quad t\ {\rm increasing}
\label{eqn:myrecur}
\end{equation}
(which may also be regarded as a numerical representation
of the \bx{differential equation} $dy/dt+y (1-\rho)/\Delta t=x(t)$.)
\par
When it comes time to think about the adjoint, however,
it is easier to think of equation~(\ref{eqn:leakinteq}) than of~(\ref{eqn:myrecur}).
Let the matrix of equation~(\ref{eqn:leakinteq}) 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, as in
\begin{equation}
\tilde x_{t-1} \eq \rho \tilde x_{t} + y_{t-1}
\quad
\quad
\quad t \ {\rm decreasing}
\label{eqn:backrecur}
\end{equation}
Thus the adjoint of causal integration
is \bx{anticausal integration}.
\par
A module to do these jobs is \texttt{leakint}.
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.
Figure \ref{fig:causint} illustrates the program for $\rho = 1$.
\opdex{causint}{causal integral}{35}{46}{filt/lib}
\par
%\activesideplot{causint}{width=3.00in,height=3.2in}{}
\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 a negative box.
  {\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{notforlecture}
\begin{exer}

\item
Consider the matrix
\begin{equation}
  \left[
        \begin{array}{cccccccc}
             1 &     0 &     0 &     0 &     0 &     0 &     0 \\
             0 &     1 &     0 &     0 &     0 &     0 &     0 \\
             0 &  \rho &     1 &     0 &     0 &     0 &     0 \\
             0 &     0 &     0 &     1 &     0 &     0 &     0 \\
             0 &     0 &     0 &     0 &     1 &     0 &     0 \\
             0 &     0 &     0 &     0 &     0 &     1 &     0 \\
             0 &     0 &     0 &     0 &     0 &     0 &     1 
        \end{array}
  \right]
\end{equation}
and others like it with $\rho$ in other locations.
Show what combination of these matrices will represent
the leaky integration matrix
in equation (\ref{eqn:leakinteq}).  What is the adjoint?

\item
Modify the calculation in Figure~\ref{fig:causint} so that there is
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 \texttt{leakint} \vpageref{/prog:leakint} by changing the diagonal to contain
1/2 instead of 1.
Notice how time alignment changes in Figure~\ref{fig:causint}.
\end{exer}
%\end{notforlecture}




\subsection{Backsolving, polynomial division and deconvolution}

\par
Ordinary differential equations often lead us to the backsolving operator.
For example, the damped harmonic oscillator leads to
a special case of equation
(\ref{eqn:polydiv})
where $(a_3,a_4,\cdots)=0$.
There is a huge literature on finite-difference solutions
of ordinary differential equations
that lead to equations of this type.
Rather than derive such an equation on the basis
of many possible physical arrangements,
we can begin from the filter transformation $\bold B$ in (\ref{eqn:contran1})
but put the matrix on the other side of the equation
so our transformation can be called one of inversion or backsubstitution.
Let us also force the matrix $\bold B$ to be a square matrix
by truncating it with
$\bold T = \left[ \bold I \quad \bold 0 \right]$, say
$\bold A =\left[ \bold I \quad \bold 0 \right] \bold B = \bold T \bold B$.
To link up with applications in later chapters,
I specialize to 1's on the main diagonal and insert some bands of zeros.
\begin{equation}
\bold A \bold y
\eq
\left[ 
\begin{array}{ccccccc}
   1  & 0   & 0    & 0   & 0   & 0   & 0  \\
  a_1 &  1  & 0    & 0   & 0   & 0   & 0  \\
  a_2 & a_1 &  1   & 0   & 0   & 0   & 0  \\
  0   & a_2 & a_1  &  1  & 0   & 0   & 0  \\
  0   & 0   & a_2  & a_1 &  1  & 0   & 0  \\
  a_5 & 0   & 0    & a_2 & a_1 &  1  & 0  \\
  0   & a_5 & 0    &   0 & a_2 & a_1 & 1
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  y_0 \\ 
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6
  \end{array} \right]
\eq
\left[ 
\begin{array}{c}
  x_0 \\ 
  x_1 \\ 
  x_2 \\ 
  x_3 \\ 
  x_4 \\ 
  x_5 \\ 
  x_6
  \end{array} \right] 
\eq
\bold x
\label{eqn:polydiv}
\end{equation}
Algebraically, this operator goes under the various names,
``backsolving'',
``polynomial division'', and
``deconvolution''.
The leaky integration transformation
(\ref{eqn:leakinteq})
is a simple example of backsolving
when $a_1=-\rho$ and $a_2=a_5=0$.
To confirm this, you need to verify that the matrices in
(\ref{eqn:polydiv})
and
(\ref{eqn:leakinteq})
are mutually inverse.

\par
A typical row in equation (\ref{eqn:polydiv}) says
\begin{equation}
x_t \eq y_t \ +\ \sum_{\tau > 0} \ a_\tau\ y_{t-\tau}
\label{eqn:conrecur}
\end{equation}
Change the signs of all terms in equation
(\ref{eqn:conrecur})
and
move some terms to the opposite side
\begin{equation}
y_t \eq x_t \ -\ \sum_{\tau > 0} \ a_\tau\ y_{t-\tau}
\label{eqn:deconrecur}
\end{equation}
Equation (\ref{eqn:deconrecur}) is a recursion
to find $y_t$ from the values of $y$ at earlier times.

\par
In the same way that
equation (\ref{eqn:contran1})
can be interpreted as $Y(Z)=B(Z)X(Z)$,
equation (\ref{eqn:polydiv})
can be interpreted as $A(Z) Y(Z) = X(Z)$
which amounts to
$Y(Z) = X(Z)/A(Z)$.
Thus, convolution is amounts to polynomial multiplication
while the backsubstitution we are doing here
is called deconvolution, and it amounts to polynomial division.




\par
A causal operator is one that uses its present and past inputs
to make its current output.
Anticausal operators use the future but not the past.
Causal operators are generally associated with lower triangular matrices
and positive powers of $Z$, whereas
anticausal operators are associated with upper triangular matrices
and negative powers of $Z$.
A transformation like equation
(\ref{eqn:polydiv})
but with the transposed matrix would require us
to run the recursive solution the opposite direction in time,
as we did with leaky integration.


\par
A module to backsolve~(\ref{eqn:polydiv}) is \texttt{recfilt}.
\opdex{recfilt}{deconvolve}{48}{76}{filt/lib}

\par
The more complicated an operator, the more complicated is its adjoint.
Given a transformation from $\bold x$ to $\bold y$ that is
$\bold T\bold B \bold y = \bold x$,
we may wonder if the adjoint transform really is
$(\bold T\bold B)' \hat{\bold x}= \bold y$.
It amounts to asking if the adjoint of
$\bold y = (\bold T\bold B)^{-1}\bold x$ is
$\hat{\bold x} = ((\bold T\bold B)')^{-1}\bold y$.
Mathematically we are asking if
the inverse of a transpose is the transpose of the inverse.
This is so because in
$\bold A\bold A^{-1} = \bold I = \bold I' = (\bold A^{-1})'\bold A'$
the parenthesized object must be the inverse of its neighbor $\bold A'$.



\par
The adjoint has a meaning which is nonphysical.
It is like the forward operator except that we must
begin at the final time and revert towards the first.
The adjoint pendulum damps as we compute it backward in time,
but that, of course, means that the adjoint pendulum diverges
as it is viewed moving forward in time.


\subsection{The basic low-cut filter}
Many geophysical measurements contain
very low-frequency noise called ``drift.''
For example, it might take some months to survey the depth of a lake.
Meanwhile, rainfall or evaporation could change the lake level so that
new survey lines become inconsistent with old ones.
Likewise, gravimeters are sensitive to atmospheric pressure,
which changes with the weather.
A magnetic survey of an archeological site would need to contend
with the fact that the earth's main magnetic field is changing randomly
through time while the survey is being done.
Such noises are sometimes called ``secular noise.''

\par
The simplest way to eliminate low frequency noise is
to take a time derivative.
A disadvantage is that the derivative
changes the waveform
from a pulse to a doublet (finite difference).
Here we examine the most basic low-cut filter.
It preserves the waveform at high frequencies;
it has an adjustable parameter
for choosing the bandwidth of the low cut;
and it is causal (uses the past but not the future).

\par
We make our causal lowcut filter (highpass filter) by