%\section{THEORY}
%\subsection{Regularization}
In general, geophysical inverse problems (inverting for some
model ($\model$), given some ($\data$),
while applying some operator ($\lop$))
are ill-posed.
To obtain a pleasing result we impose some type of regularization
criteria  such as
adding an additional model regularization term to our
objective function,

\beq
\min_{m} | \lop \model - \data|^2 +\epsilon |\reg \model|^2  
\eeq

or, alternately, described by the fitting goal
\beqa
\data &\approx& \lop \model \\ \nonumber
\zero &\approx& \epsilon  \reg \model
\label{eq:basic}
\eeqa

Ideally, $\reg$ should be the inverse of the model covariance matrix
\cite{tarantola-ipt}.
Unfortunately, we don't have an obvious mean of getting the model covariance
matrix.
      The lack of knowledge about the model often leads to the Laplacian
or some other isotropic
operator being used for $\reg$.
As a result, we fill the null-space of the model with isotropic features, that,
while explaining the data,
may be unreasonable when judging the results
from a geologic perspective.
Fortunately, we often do have other informational sources that 
can give us at least a first-order approximation to the model covariance matrix.

%NOW WE NEED TO GET INTO PRECONDITIONING
%\subsection{Preconditioning}
Another important consideration is the speed of convergence of the problem.  
The size of most geophysical
problems make direct matrix inversion methods impractical. 
An appealing alternative for linear problems is the family of
conjugate gradient methods.  Unfortunately, 
the operators used in seismic reflection  problems
are often computationally expensive.
As a result it is important to minimize the number of steps it takes
to get to a reasonable solution.
One of the reasons for the slow convergence is that our regularization operator
($\reg$)  is small, which means information will be spread only
a short distance with each iteration, slowing down convergence.
One way that we can reduce the number of iterations is by reformulating the 
problem in terms of some new variable ($\pvar$),
with an inexpensive preconditioning
operator ($\bold B$) that smoothes long distances in each iteration.

Now the question becomes what to use for $\bold B$.  The ideal
choice for $\bold B$ would be $\prec$, because if $\bold B = \prec$,
the  fitting goals (\ref{eq:basic}) become 
\beqa
\data &\approx& \lop \prec \pvar \\ \nonumber
\zero &\approx& \reg \model = \reg \prec \pvar =  \iop \pvar .
\eeqa
This new set of fitting goals is advantagous 
because they are 
equivalent to our initial goals {\bf and} should converge 
more quickly.
The speed up is due not only to this magical
preconditioner spreading information long-distances
with every ittteration but also  to our regularization
equation now being the identity matrix, so the
conjugate algorthym does not have to spend any
energy trying to invert it\cite{gee}.


But, before we can implement these fitting goals we must 
\begin{itemize}
\item find a
regularization operator that will steer our model to follow geologic dip and
\item find, and apply, efficiently the inverse of this operator.
\end{itemize}


%\subsection{Steering filters}
Schwab\shortcite{Schwab.sep.94.matt2} showed that plane waves with a given slope
can be predicted on a discrete grid by compact filters.  Using these same
compact filters (Figure~\ref{fig:steering}) as a regularization operator we 
can steer our solution towards a given dip (Figure \ref{fig:sweep}). 
By changing the number of coefficients we use to describe the dip, we
can change the dip procession of our plane waves. This representation is
also easily adaptible.  Figure~\ref{sweep} shows that by
changing the location of our filter coeficients
we can destroy different angled plane-waves. And by changing the number
of coeficients we use to describe a given angle, we can create more issotropic
features in our data (Figure~\ref{fig:size}).

Unfortunately, 
for anything but the simplest of regimes, the velocity gradient can not
be described by a single dip, but instead varies spatially.  As a result
we must introduce one more level of sophistication, a {\it space varying}
regularization operator(Figure~\ref{fig:filters}.)  

\sideplot{steering}{height=1.0in}{An example plane-wave anhilation filter that will anhilate waves at 45 degree slope.}

\plot{sweep}{width=4.0in}{The result of using plane-wave
anhilators  as a regularization operator.  The resulting model has
preferentially orriented along the direction at which the plane-wave
anhilator was constructed.}

Given early migrated reflector position, a geologists guess on structure,
or even well log dip measurements we can make some guess on the general
velocity gradient throughout our model space.  Then, for each model point,
we can construct a small plane wave annihilation ( or {\it steering}) filter, 
oriented along our estimate at 
model gradient.   The combination of all of these filters becomes our
regularization operator (Figure \ref{fig:filters}.

\plot{filters}{width=4.0in}{Composite steering filter operator.  Each model
point has its one steering filter orriented along the desired velocity gradient.}

To give you an example of the ability of this operator to steer a solutions,
lets look at a simple missing data problem.
We will start with a slice form a rather complex synthetic velocity model 
(Figure~\ref{fig:qdome-geophysics}.  that has a dome -feature overlaying
an angular unconformity.
We will chose, some representative velocity boundaries,
calculate their dip at each point along the boundary, and then smoothly spread
this dip information to every point in the model space, Figure~\ref{fig:qdome-geophysics}.
We simulate nine well logs along the survey 
by subsampling the velocity field.
Next we set up our fitting goals for our problem.
We want to honor our data, so we won't allow the model to change at known
locations ($\bold J$ will  be a diagonal masking operator with 
zeros at known  model locations, ones otherwise), and we want to regularize
with our steering filter operator ($\reg$),

\plot{qdome-geophysics}{width=6.0in}{Left panel represent the input velocity
field with some representative reflectors highlighted.  The center panel
is the result of calculating the dip of these reflectors and smoothing
into the entire model space.  The right panel is the input data to our
interpolation problem.}

\beqa
\data &\approx& (\iop - \bold J) \model \\ \nonumber
\zero &\approx& \reg \model .
\eeqa

After 120 iterations, you can see that we did an excellent job recovering
our initial velocity model (Figure~\ref{fig:sds}), 
for comparison I have included the result
of using the Laplacian as our operator.

%\subsection{Precondioner from the helix}
To invert this  multi-dimensional operator, 
involves a little trick that image processors
have known for a long time (cite jon's image procesing book)
it has only been recently used
in geophysics 
\cite{helix}.  Given a multi-dimensional filter,
such as Figure \ref{fig:filter}, you first 
map it onto the coordinate space of the data (panel a of Figure~\ref{fig:helix}). Then imagine wrapping the data arround a cylinder with the end of column one
touch the begining of column 3.  Finally, unwind the data into a single
string of numbers and you have converted your multi-dimensional filter into
a 1-D filter!
If this new, one-dimensional
filter is causal and minimum phase, which our steering filters
our, you can apply the inverse of the filter, Claerbout\shortcite{fgdp} shows
that  cheaply through polynomial
division.
\plot{helix}{width=6.0in}{ Filtering on a helix.
        The same filter coefficients overlay the same data values
        if the 2-D coils are unwound into 1-D strips.}
Now that we have the inverse filter, we can change our regularized 
interpolation problem, into a preconditioned problem,
\beqa
\data &\approx& (\iop - \bold J) \prec \pvar \\ \nonumber
\zero &\approx& \pvar .
\eeqa
Applying this new set of equations to the same problem, we get
convergence, to approximately the same model in 12, rather than 120 
iterations, a two orders of magnitude speed up!
