
\sepinput{macros}
% Started 09/22/99

%\shortnote

\lefthead{Clapp \& Brown}
\righthead{Multiple suppression}
\footer{SEP--102}

\title{Applying SEP's latest tricks to the multiple suppression problem}

\email{bob@sep.stanford.edu,morgan@sep.stanford.EDU}
\keywords{multiples, Huber, signal-noise separation}

\author{Robert G. Clapp and Morgan Brown}

\ABS{
Two methods for suppressing multiples are proposed.
In  the first, 
multiple suppression is expressed as a signal to noise separation
problem.  The problem is  solved in the time domain using space-varying
Prediction Error Filters (PEFs).  
The second method shows how greater separation between primaries
and multiples can be obtained by velocity space inversion using
a Huber functional rather than the standard $L_2$ functional.
Early results of both methods are encouraging.
}

\section{INTRODUCTION}
Multiple suppression is one of the biggest problems facing the seismic industry.
Methods that have proven effective in 2-D  are either
cost-prohibitive or not easily extendible into 3-D \cite{berkhout.geo.97,verschuur.geo.97,Sun.sep.100.yalei1}.
\longcite{TLE18-1-00550058}  
proposed forming the multiple suppression as a signal-to-noise separation
in the frequency domain, but this method suffered from stability problems.
\par
Until recently, an equivalent time domain method was not possible. 
Claerbout\shortcite{helix} discovered that multi-dimensional time filters 
can be mapped into 1-D, therefore making it possible to do inverse filtering.
\longcite{Crawley.sep.97.sean2} 
showed how non-stationary filters could more accurately
predict seismic data.  \longcite{Fomel.sn} demonstrated how 
Spitz's method could
be changed to work with time domain PEFs.
\par
In the first section of the paper, we
perform time domain multiple suppression by a two step method. We first
estimate a space-varying PEF  from data (a CMP gather) and a noise model
(an estimate of the multiples obtain by   
downward continuing through the water column twice).  
We then separate  out the signal (primaries)
from the noise (multiples) by a simple inversion scheme. 
\par
In the second portion of the paper we present a better way to separate
multiples in velocity space.
\longcite{Lumley.sep.82.25}  described a CMP gather as a  sum of hyperbolic
events.  They then inverted  this velocity-space transform
into  ($\tau$,$v$) space, muted multiples,
and transformed back into ($t,h$) space.  \longcite{Guitton.sep.100.antoine1}
showed that a Huber functional \cite{Huber:73} produces a velocity scan where
reflection energy is better behaved.
We invert into ($\tau,v$) using the Huber functional,
rather than an $L_2$ functional.  We show that the Huber method provides
more separation between primary and multiple trends, and therefore improved
multiple suppression.

\par

\section{MULTIPLE SUPPRESSION USING SIGNAL-NOISE SEPARATION}
Signal to noise separation has a long history at SEP  \cite{ Harlan.sepphd.47,Kostov.sepphd.63,Claerbout.sep.71.270}.
The method we use is similar to Abma's \shortcite{Abma.sepphd.88} formation.
\longcite{Abma.sepphd.88} 
proposed solving the set of equations:
\beqa
  \npef \noise &\approx\ \zero \nonumber \\
  \epsilon \spef \signal &\approx\ \zero   \label{eq:abma} \\
    \mbox{subject to} &\leftrightarrow& \bf d = s+n \nonumber 
\eeqa
  where the operators $\bf N$ and $\bf S$ represent $t-x$ domain convolution with
(PEF's) which decorrelate the unknown noise $\noise$ and signal $\signal$,
 respectively, and the factor $\epsilon$ balances the energies of the residuals.

For his problem he assumed that the noise was uncorrelated, therefore $\npef$
becomes the identity and $\spef$ is  the PEF that best predicted the
data in a given window [patching approach \cite{Claerbout.sep.73.391,Schwab.sep.84.271}].
\par
In the multiple problem the noise is not uncorrelated so we must find
another way to find $\npef$.  \longcite{TLE18-1-00550058} proposed defining
$\spef$ as $\spef =  \dpef \npef^{-1}$ where $\dpef$ is a filter that
      characterizes the data rather than the signal.  Using this new definition
we get a new set of fitting goals:
\beqa
\npef \signal &\approx& \npef \data \nonumber \\
\epsilon \dpef \npef^{-1} &\approx& \zero  .
\eeqa
Following \longcite{Fomel.sep.95.sergey1} we can set 
up the conversion  by reformulating it as a preconditioned problem
by a simple change of variables ($\pvar=\dpef \npef^{-1}$)
\beqa
\npef \npef \dpef^{-1} \pvar &\approx& \bf \npef \data \nonumber \\
\epsilon \pvar &\approx& \zero \label{eq:signoi} ,
\eeqa
where $\pvar$ is just a dummy preconditioning variable.
\par
Instead of using patching we followed the methodology of
\longcite{Crawley.sep.97.sean2} and constructed and estimated
a space varying filter.
\beqa
\zero &\approx& \dpef \prec \pvar \nonumber \\
\zero &\approx& \pvar \label{eq:destroy}
\eeqa
where $\reg$ is a radial smoother \cite{Clapp.sep.100.bob3}.
For $\npef$ we follow a similar procedure assuming an a priori  model for the
noise.

\subsection{Synthetic example}
At this early stage we decided to show how the method works on
a simple synthetic. Figure \ref{fig:bee-cmp}  shows a sea floor
with a series of water-bottom multiples.  The amplitude of the
reflector increases as a function of angle, something that
frequency methods have had difficult time handling.
The multiple model was
constructed by downward continuing \cite{Bevc.sepphd.87} the sea floor reflection,
 Figure \ref{fig:bee-mult}.  
As a result the
amplitude information is incorrect, but the kinematics is correct.

\sideplot{bee-cmp}{width=3.0in}{A simple CMP gather with multiples.}
\sideplot{bee-mult}{width=3.0in}{CMP gather in Figure \ref{fig:bee-cmp} 
datumed to the surface and back to the seafloor.
Note how the amplitude information is
wrong but the general kinematics is correct.}
\par
Once we have our data and noise model we estimated a space varying  filter
for each by applying fitting goal (\ref{eq:destroy}). 
To conserve memory we  put a new filter every 15th point in time
and third point in offset (Figure~\ref{fig:pef}).
These two filters were then used, and
fitting goal (\ref{eq:signoi}) were applied.  Figure \ref{fig:signoi} shows
the result of the separation. We can see some residuals of
the filter patches but generally we have done a good job in removing
the multiple while preserving the signal.  

\plot{signoi}{width=5.0in}{The estimated primaries and mulitpes from
Figure \ref{fig:bee-cmp}.}

\plot{pef}{width=3.5in}{Space varying filter composition. A different
filter is placed inside each patch.  The filter estimation problem
is done globally, with the filter coefficient smoothed in a radial
direction.}

\subsection{Problems}
At this preliminary stage there are still some problems to overcome.
The separation technique is much more reliant on a 
stable filter than the interpolation\cite{Crawley.sep.102.sean1} or
feature identification problem\cite{Brown.sep.102.morgan2}.  
Some preliminary work on filter stability was done in \cite{Sava.sep.100.paul3,Sava.sep.97.paul1}, but the problems associated with non-stationary filters
are only beginning to be addressed \cite{Rickett.sep.102.james1}.
To be effective in 3-D, the problem of spatial aliasing
has to be dealt with. 


\section{HUBER VELOCITY SPACE  MULTIPLE ELIMINATION}
Transforming seismic data into another domain and then muting is
a common multiple suppression technique.
The common method is the  parabolic
radon transform \cite{GEO57-03-03860395}.
This method has the advantage of having a analytic inverse (and is therefore 
faster), but involves approximating moveouts by parabolas.
\longcite{Lumley.sep.82.25} used  the more expensive hyperbolic transform and
went a step further by forming it as an inversion problem,
\beq
\data \approx {\bf H} {\bf V} \model \label{eq:scan} ,
\eeq
where:
\begin{description}
\item [$\data$] is the CMP gather, 
\item [$\bf H$] is a half derivative operator \cite{Prucha.sep.102.marie1},
\item [$\bf V$] is a velocity transform operator.
\end{description}
They then used an $L_2$ conjugate gradient algorithm to estimate
a ($\tau,v$) model and then muted out the multiples.
Unfortunately both methods suffer from multiple and primary energy overlapping.

\subsection{Huber separation}
To get a better ($\tau,v$) model we decided to replace the linear iterative
solver used by \longcite{Lumley.sep.82.25} with the Fletcher-Reeves non-linear
conjugate gradient conditions \cite{opbook} and  the Dennis-Schnabel
line search method \cite{den}.
We replaced the $L_2$ function, with a Huber functional \cite{Huber:73} 
that is
less sensitive to large outliers.
The Huber
functional is $L_2$ until some cutoff value and then smoothly switches
to $L_1$ (Figure~\ref{fig:huber}).  
The idea is  compromise between the convergence speed of
$L_2$ and the less sensitive nature to outliers with the $L_1$.
Guitton and Symes\shortcite{Guitton.sep.100.antoine1} showed that the Huber 
functional does a better job of localizing energy in ($t,v$) space.  For
multiples this mean that the primary and multiple trains are
better separated. 
\par
\sideplot{huber}{width=3in}{The $L_2$, Huber, and $L_1$ functionals.}

To show the advantage of the Huber method over a straight $L_2$ we 
took a multiple contaminated CMP gather, Figure \ref{fig:mobil}, and iterated
on fitting goals (\ref{eq:scan}).  
Figure \ref{fig:compare} shows the envelope of the
($\tau,v$) representation of both the Huber and $L_2$ approach.  Note
how the Huber result shows a more compact representation of the
primary and multiple trains. 

\sideplot{mobil}{width=2.5in,height=4.0in} { A multiple infested gather from the Mobil AVO dataset.}
\plot{compare}{width=6.0in}{The envelope of the tau-velocity space representation
of the Mobil AVO gather (Figure \ref{fig:mobil}) using
both an $L_2$, left, and a Huber, right, functional.  
Note how the Huber gather shows
more energy and better isolation of the primary train.  Further, the $L_2$
approach  show significantly more energy at high, unreasonable, velocities.}


\bibliographystyle{sep}
\bibliography{SEP,MISC,bob}

%\APPENDIX{A}

%\plot{name}{width=6in,height=}{caption}
%\sideplot{name}{height=1.5in,width=}{caption}
%
%\begin{equation}
%\label{eqn:}
%\end{equation}
%
%\ref{fig:}
%(\ref{eqn:})
