paper.tex

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

\bold 0 \eq {\partial Q \over \partial \bold x'}  \eq
\bold F' (\bold F\bold x - \bold d)
\label{eqn:normeq}
\end{equation}
The result is merely the complex form of
our earlier result (\ref{eqn:fitorth1}).
Therefore,
differentiating by a complex vector
is an abstract concept,
but it gives the same set of equations
as differentiating by each scalar component,
and it saves much clutter.
%\end{notforlecture}

%\begin{notforlecture}
\subsection{From the frequency domain to the time domain}
Equation~(\ref{eqn:z4}) is a frequency-domain quadratic form
that we minimized by varying a single parameter,
a Fourier coefficient.
Now we will look at the same problem in the time domain.
We will see that the time domain offers flexibility with
boundary conditions, constraints, and weighting functions.
The notation will be that a filter $f_t$ has input $x_t$ and output $y_t$.
In Fourier space this is $Y=XF$.
There are two problems to look at,
unknown filter $F$ and unknown input $X$.

\subsubsection{Unknown filter}
When inputs and outputs are given,
the problem of finding an unknown filter appears to be overdetermined,
so we write $\bold y \approx \bold X \bold f$
where the matrix $\bold X$ is a matrix of downshifted columns like
(\ref{ajt/eqn:contran2}).
Thus the quadratic form to be minimized
is a restatement of equation~(\ref{eqn:6-1-23})
with filter definitions:
\begin{equation}
Q(\bold f', \bold f) \eq
(\bold X\bold f - \bold y)'
(\bold X\bold f - \bold y)
\end{equation}
The solution $\bold f$ is found just as we found
(\ref{eqn:normeq}),
and it is the set of simultaneous equations
$ \bold 0 = \bold X'(\bold X\bold f - \bold y)$.

\subsubsection{Unknown input: deconvolution with a known filter}
\sx{deconvolution}
For solving the unknown-input problem,
we put the known filter $f_t$ in a matrix of downshifted columns $\bold F$.
Our statement of wishes is now to find $x_t$ so that
$\bold y \approx \bold F \bold x$.
We can expect to have trouble finding unknown inputs $x_t$
when we are dealing with certain kinds of filters,
such as \bx{bandpass filter}s.
If the output is zero in a frequency band,
we will never be able to find the input in that band
and will need to prevent $x_t$ from diverging there.
We do this by the statement that we wish
$\bold 0\approx\epsilon\,\bold x$,
where $\epsilon$ is a parameter that is small
and whose exact size will be chosen by experimentation.
Putting both wishes into a single, partitioned matrix equation gives
\begin{equation}
 \left[ 
  \begin{array}{c}
   \bold 0 \\ 
   \bold 0
  \end{array}
 \right] 
\quad\approx\quad
 \left[ 
  \begin{array}{c}
   \bold r_1 \\ 
   \bold r_2
  \end{array}
 \right] 
\eq
 \left[ 
  \begin{array}{c}
   \bold F \\ 
   \epsilon \  \bold I
  \end{array}
 \right] 
 \ 
 \bold x
\quad - \quad
 \left[ 
  \begin{array}{c}
   \bold y \\ 
   \bold 0
  \end{array}
 \right] 
\end{equation}
To minimize the residuals $\bold r_1$ and $\bold r_2$,
we can minimize the scalar
$\bold r' \bold r = \bold {r'}_1\bold r_1 + \bold {r'}_2\bold r_2$.
This is
\begin{eqnarray}
Q(\bold x', \bold x) &=& (\bold F  \bold x  - \bold y)' (\bold F\bold x-\bold y)
                                        + \epsilon^2 \bold x' \bold x
                                                                \nonumber \\
                     &=& (\bold x' \bold F' - \bold y') (\bold F\bold x-\bold y)
                                        + \epsilon^2 \bold x' \bold x
\label{eqn:knownfilt}
\end{eqnarray}
We solved this minimization
in the frequency domain
(beginning from equation~(\ref{eqn:z4})).
\par
Formally the solution is found just as with equation (\ref{eqn:normeq}),
but this solution looks unappealing in practice
because there are so many unknowns and because
the problem can be solved much more quickly
in the Fourier domain.
To motivate ourselves to solve this problem in the time domain,
we need either to find an approximate solution method that is
much faster, or to discover that
constraints or time-variable weighting functions
are required in some applications.
This is an issue we must be continuously alert to,
whether the cost of a method is justified by its need.







%\item
%Try other lags in~\EQN{2conv} such as
%$(0,1,0)'$ and 
%$(0,0,1)'$.
%Which works best?
%Why?

%\end{notforlecture}

\begin{exer}
\item
\sx{zero absolute temperature}
In 1695, 150 years before Lord Kelvin's absolute temperature scale,
120 years before Sadi Carnot's PhD thesis,
40 years before Anders Celsius,
and 20 years before Gabriel Farenheit,
the French physicist Guillaume
\bx{Amontons},
deaf since birth,
took a mercury manometer (pressure gauge) and
sealed it inside a glass pipe (a constant volume of air).
He heated it to the boiling point of water at $100^\circ$C.
As he lowered the temperature to freezing at $0^\circ$ C,
he observed the pressure dropped by 25\% .
He could not drop the temperature any further
but he supposed that if he could drop it further by a factor of three,
the pressure would drop to zero (the lowest possible pressure)
and the temperature would have been the lowest possible temperature.
Had he lived after Anders Celsius he might have calculated
this temperature to be $-300^\circ$ C (Celsius).
Absolute zero is now known to be $-273^\circ$ C.
\par
It is your job to be Amontons' lab assistant.
Your $i$th measurement of temperature
$T_i$ you make with Issac Newton's thermometer and
you measure pressure $P_i$ and volume $V_i$ in the metric system.
Amontons needs you to fit his data with the regression
$
0\approx \alpha(T_i-T_0)-P_i V_i
$
and calculate the temperature shift $T_0$ that Newton should have made
when he defined his temperature scale.
%In the theory of least squares fitting, what is the data?
%What are the fitting functions?
%and what are the two model parameters?
Do not solve this problem!
Instead, cast it in the form of equation (\ref{eqn:wish1}),
identifying the data $\bold d$ and the two column vectors
$\bold f_1$ and
$\bold f_2$
that are the fitting functions.
Relate the model parameters $x_1$ and $x_2$
to the physical parameters $\alpha$ and $T_0$.
Suppose you make ALL your measurements at room temperature,
can you find $T_0$?  Why or why not?

\end{exer}

\section{KRYLOV SUBSPACE ITERATIVE METHODS}
The \bx{solution time} for simultaneous \bx{linear equations}
grows cubically with the number of unknowns.
There are three regimes for solution;
which one is applicable
depends on the number of unknowns $m$.
For $m$ three or less, we use analytical methods.
We also sometimes use analytical methods on matrices of size $4\times 4$
when the matrix contains enough zeros.
Today in year 2001,
a deskside workstation, working an hour solves about a
$4000\times 4000$ set of simultaneous equations.
A square image packed into a 4096 point vector is a $64\times 64$ array.
The computer power for linear algebra to give us solutions that
fit in a $k\times k$ image is thus proportional
to $k^6$, which means that even though computer power grows rapidly,
imaging resolution using ``exact numerical methods'' hardly
grows at all from our $64\times 64$ current practical limit.

\par
The retina in our eyes captures an image of size about $1000\times 1000$
which is a lot bigger than $64\times 64$.
Life offers us many occasions where final images exceed the 4000
points of a $64\times 64$ array.
To make linear algebra (and inverse theory) relevant to such problems,
we investigate special techniques.
A numerical technique known as the
``\bx{conjugate-direction method}''
works well for all values of $m$ and is our subject here.
As with most simultaneous equation solvers,
an exact answer (assuming exact arithmetic)
is attained in a finite number of steps.
And if $n$ and $m$ are too large to allow enough iterations,
the iterative methods can be interrupted at any stage,
the partial result often proving useful.
Whether or not a partial result actually is useful
is the subject of much research;
naturally, the results vary from one application to the next.

\subsection{Sign convention}
On the last day of the survey, a storm blew up,
the sea got rough, and the receivers drifted further downwind.
The data recorded that day
had a larger than usual difference
from that predicted by the final model.
We could call
$(\bold d-\bold F\bold m)$
the {\it \bx{experimental error}}.
(Here
$\bold d$ is data,
$\bold m$ is model parameters, and
$\bold F$ is their linear relation).

\par
The alternate view is that our theory was too simple.
It lacked model parameters for the waves and the drifting cables.
Because of this model oversimplification
we had a {\it \bx{modeling error}} of the opposite polarity
$(\bold F\bold m-\bold d)$.
\par
A strong experimentalist prefers to think of the error
as experimental error, something for him or her to work out.
Likewise a strong analyst likes to think
of the error as a theoretical problem.
(Weaker investigators might be inclined to take the opposite view.)

\par
Regardless of the above, and opposite to common practice,
I define the \bx{sign convention} for the error (or residual) as
$(\bold F\bold m-\bold d)$.
When we choose this sign convention,
our hazard for analysis errors will be reduced
because $\bold F$ is often complicated and formed by combining many parts.
\begin{comment}
So in this book
we see positive signs on operators
and we see residuals initialized by the negative of the data,
often with subroutine \texttt{negcopy()}.
\progdex{negcopy}{copy and negate}
\end{comment}
\par
\boxit{
        Beginners often feel disappointment
        when the data does not fit the model very well.
        They see it as a defect in the data
        instead of an opportunity
        to design a stronger theory.
        }
\par

\subsection{Method of random directions and steepest descent}
\sx{random directions} \sx{steepest descent}
\par
Let us minimize the sum of the squares of the components
of the \bx{residual} vector given by
%\begin{eqnarray}
%{\rm residual}
%&=&
%{\rm 
%\hbox{transform} \ \ \  \ \hbox{model space}
%\ \ -\ \  \
%\ \hbox{data space}
%}
%        \\
%\left[
%\matrix { \matrix {   \cr \cr \bold r  \cr \cr \cr } }
%\right]
%       &=&
%       \ \ 
%\left[
%  \matrix {
%    \matrix { \cr \cr         \cr \cr \cr }
%    \matrix { \cr \cr         \cr \cr \cr }
%    \matrix { \cr \cr         \cr \cr \cr }
%    \matrix { \cr \cr \bold F \cr \cr \cr }
%    \matrix { \cr \cr         \cr \cr \cr }
%    \matrix { \cr \cr         \cr \cr \cr }
%    \matrix { \cr \cr         \cr \cr \cr }
%    }
%  \right]
% \ \ \ \
%\left[
%  \matrix {
%    \matrix {  \cr \bold x \cr \cr }
%    }
%  \right]
%  \ \ \ 
%\ \ \ \ \ -\ \ \ \ \ 
%\left[
%  \matrix { 
%    \matrix {  \cr  \cr \bold d \cr  \cr  \cr}
%    }
%  \right]
%\label{eqn:cg1b}
%\end{eqnarray}
%
%
\begin{eqnarray}
{\rm residual}
&=&
{\rm 
\quad \hbox{transform} \quad \quad \hbox{model space}
\quad - \quad
\hbox{data space}
}
        \\
\left[
\begin{array}{c}
	  \\
	  \\
	  \\
	 \bold r \\
	  \\
	  \\
	 \  
\end{array}
\right]
       &=&
\left[
  \begin{array}{ccc}
     \quad & \quad   &  \quad \\
     \quad & \quad   &  \quad \\
     \quad & \quad   &  \quad \\
     \quad & \bold F &  \quad \\
     \quad & \quad   &  \quad \\
     \quad & \quad   &  \quad \\
     \quad & \quad   &  \quad 
  \end{array}
\right]
          \quad
\left[
  \begin{array}{c}
	  \\
	  \\
	 \bold x \\
	  \\
	 \ 
  \end{array}
\right]
	\quad - \quad
\left[
  \begin{array}{c}
	  \\
	  \\
	  \\
	 \bold d \\
	  \\
	  \\
	 \