% This texfile was generated by executing Latify on file paper_act.tex
\documentclass[12pt]{sep2e}
\begin{document}

\cleardoublepage
\thispagestyle{septwo}
\label{.start}
\def\papername{}
\def\figdir{./Fig}
\chapter{Image space multiple prediction}\label{chapter:DepthMult}
\section{summary}
A very important aspect of removing multiples from seismic data is 
accurate prediction of their kinematics.  I cast the multiple 
prediction problem as an operation in the image space parallel
to the conventional surface-related multiple prediction 
methodology.  Though developed in the image domain, the technique 
shares the data-driven strengths of data-domain surface-related
multiple elimination (SRME) by being independent of the Earth 
(velocity) model.  Also, the data are used to predict the kinematics
of the multiples 
exactly so that a Radon transform need not be designed to separate 
the two types of events.  The cost of the prediction is approximately 
the same as that of data-space methods, though it can be computed during 
the course of migration.  The additional cost is not significant 
compared to that incurred by shot-profile migration.
\par
Image-space multiple predictions are generated by convolving 
the traces in each shot-gather at every depth level during the course 
of a shot-profile migration.  The prediction in the image 
domain is equivalent to that produced by migrating the data-space 
convolutional prediction.  Adaptive subtraction of the prediction from
the image is required.  Subtraction in the image domain, however,
provides the advantages of focused energy in a smaller domain since
extrapolation removes some of the imperfections of the input data.
%%\end{abstract}

\section{Introduction}
Removing multiples from seismic data is often an
imperative to producing interpretable images of the subsurface. Multiple
attenuation has a rich history in the geophysical literature ranging
from methods that predict the multiples 
via convolution of the data \cite{GPR14.04.03890426} to methods that use 
a differential characteristic between
primaries and multiples as a discriminator for
separating the two types of events \cite{Hampson86,TLE18-01-00400044}.
\par
Kinematic prediction of multiple reflections by convolution 
leads to amplitude and bandwidth inconsistencies.  Therefore,
direct subtraction of such predictions is not possible.  \longcite{GEO50-09-14101425} 
suggests modeling the waveform
of the multiples to subtract events at times calculated by
convolution from the data.  The surface-related multiple elimination (SRME)
method \cite{GEO57-09-11661177,GEO62-05-15861595} convolves traces within
shot-gathers to predict multiples (surface-related multiple
prediction, SRMP) followed by an iterative subtraction scheme to eliminate
them from the data.  Alternatively, after predicting multiples, via
convolution or filter-based methods, \longcite{antmult} uses a
pattern-based subtraction technique that resembles the match filter
application described in \longcite{SEG-2001-13011304}. 
\par
While the above references all operate in the data domain, authors
have also suggested removing multiples in the image space, after migration
\cite{Sava.segab1.2003,Sava.geo.mat,Alvarez.segab1.2004}.  These authors
capitalize on kinematic differentiation of the primaries and multiples and 
separate the two via a Radon transform tailored to residual moveout or
the analytic expression of over-migrated multiples, respectively.  There are
several motivators for attacking multiples in the image space.  First,
the image space is usually much smaller than the data space.  Second,
given a reasonably accurate velocity model, the kinematics 
of the image domain are simplified.  Appropriately migrated primary events
have little to no residual moveout, and multiples, migrated with
velocities that are too fast, have predictable moveout in both angle-
and offset-domain common-image gathers.  Alternatively, one could
migrate the data-domain SRMP volume to produce an image-space prediction 
of the multiples.
\par
I extend the SRMP approach to the image domain through the
commutability of wavefield extrapolation and convolution to produce a
multiple prediction in the image domain without needing to migrate two
data volumes.  My approach is directly analogous to SRMP, though the
prediction is calculated during the course of a shot-profile
migration.  The added cost of image-space surface-related
multiple prediction (IS-SRMP) is only a second imaging 
condition.  Because extrapolation dominates the cost of
shot-profile migration, IS-SRMP does not significantly increase 
its cost.   
\par
\longcite{muijs} develops an analogous strategy to IS-SRMP.  That work
requires substantial preprocessing, as its input is a previous
migration result.  Furthermore, the imaging procedure uses a
multi-dimensional Fourier domain (transformed over source position)
deconvolutional imaging condition that is iteratively performed over
frequency and position that would seem very computationally
burdensome.  Generating a multiple model in the image domain during
survey-sinking migrations was also discussed in \longcite{alison}.
Because the up-coming and down-going energy are not explicitly separate,
as in shot-profile algorithms, the authors require extra work that is
not necessary for this method.  That extra work effectively doubles
the cost and memory requirements of the (DSR) migration.  In contrast
to our method, it implicitly acts as a layer-stripping methodology
that is also capable of predicting intrabed multiples. 
\par
The multiple prediction produced in the image space with
this technique is mathematically equivalent to migrating
the multiple prediction produced by SRMP convolutions.  The predicted
multiples are then removed from the data via adaptive subtraction or
pattern-matching techniques.  The need for adaptive subtraction is 
two-fold.  First, convolution squares the wavelet in the data 
which reduces the effective bandwidth, squares the amplitude, and destroys the 
polarity of the output.  Therefore, even if a true-amplitude migration 
algorithm is available, the multiple prediction will not share the 
amplitude characteristics of the multiples in the migrated image.  Second, 
adaptive subtraction can account for imperfections including kinematic
errors and the presence of residual primaries or higher-order multiples in the
prediction.
\par 
As wave-equation migration becomes more of a commodity, iterative migration
as part of estimating the Earth's velocity has become more common.  The 
presence of multiples in the data makes velocity analysis more difficult since
events with conflicting velocity properties are present at the same time 
or depth.  For this reason, there will always be a need to remove multiples 
in the data domain.  However, if a shot-profile migration strategy is planned,
it will be shown that simultaneously producing a multiple prediction adds no
significant cost, and does not require the 
presence of multiples in the data being migrated.  Therefore, the prediction 
can be generated even after perfect elimination as a comparative volume for 
interpreters to use as an example of events that are not geologic.
\par
In this chapter, I develop the image-space surface-related multiple prediction
(IS-SRMP) technique by combining the SRMP approach with a 
wave-equation shot-profile depth migration algorithm.   \longcite{antmult} 
shows that pattern-based and adaptive subtraction of multiple models
work better when higher dimensionalities can be exploited by the
subtraction algorithm.  The various imaging conditions I present can
all produce the image as a function of subsurface offset (and from
there reflection angle) to facilitate better subtraction.  A simple 
flat-layer synthetic and the Sigsbee2B synthetic data set are used to
show the efficacy of the prediction and its subtraction from 
the migrated image.  Multiple prediction and subtraction is also
presented with a Gulf of Mexico data set acquired by Western-Geco in
the Mississippi Canyon lease area.  

\section{Theory}
The image space is the output of migration, which
I produce with a shot-profile depth migration 
algorithm.  Shot-profile wave-equation 
depth migration \cite{GEO36-03-04670481} is the cascade of two
component operations: extrapolation and imaging.  Extrapolation is
carried out with an anticausal wave-equation operator applied
to the up-going wavefield, $U$, and a causal operator 
applied to the down-going wavefield,
$D$.  $U$ is a shot record with traces placed along a wavefield axis
${\bf x}$.  $D$ is a zero valued wavefield, also defined along the
${\bf x}$-axis, where a source function is placed at the location of
the shot being migrated, ${\bf x}_s$.  The wavefields are recursively
extrapolated to all depths $z$ using one-way mixed-domain solutions to
the wave-equation 
\begin{eqnarray}\label{extr}
U_{z+1}({\bf x};{\bf x}_s,\omega)&=&E^+({\bf x},\omega)
U_{z}({\bf x};{\bf x}_s,\omega)   \\
&\mbox{and} & \nonumber  \\ 
D_{z+1}({\bf x};{\bf x}_s,\omega)&=&E^-({\bf x},\omega)
D_{z}({\bf x};{\bf x}_s,\omega) \;.
\end{eqnarray}
In this work, the form of the extrapolator $E$ used to propagate 
the wavefields is the split-step Fourier plus interpolation (SSF-PI)
\cite{SEG-1992-0917} with multiple reference velocities, though the
degree of complexity of the operator does not change the discussion
herein.  The importance of these equations is that the operator that
extrapolates a wavefield from one depth level to the next is
commutable with convolution.  Using the simple, though less accurate,
bulk phase-shift operator of the form $\mbox{e}^{i\omega s \Delta z}$ for
$E$, makes this commutativity more obvious.  The extrapolation can
also be viewed as repeatedly redatuming the shot gathers into the
image space.
\par
The correlation based multi-offset imaging condition for shot-profile
migration is \cite{GEO67-03-08830889} 
\begin{equation} \label{spic}
i_z({\bf x},{\bf h})=
\sum_{{\bf x}_s}\sum_\omega 
U_z  ({\bf x}+{\bf h};{\bf x}_s,\omega)
D_z^*({\bf x}-{\bf h};{\bf x}_s,\omega)\;
\end{equation}
where the $*$ represents conjugation, and $\bf{h}$ is subsurface 
offset.  Extraction of the zero lag of the
correlations, by summation over $\omega$, combines the energy in the
two wavefields that is collocated at each depth level.  Overlapping
acquisition patches from the individual shots are stacked by the
sum over ${\bf x}_s$.
\par
Acknowledging the approximation of convolving raw
data traces in lieu of only 
primaries, the prediction of multiples in the data space (SRMP) 
can be written in the Fourier domain 
(Berkhout and Verschuur, equation 13f,  1997)
\begin{equation}\label{issrme1}
M({\bf x}_g;{\bf x}_s,\omega)=\sum_{{\bf x}_a} 
R({\bf x}_g;{\bf x}_a,\omega)
R({\bf x}_a;{\bf x}_s,\omega)\;,
\end{equation}
where $R$ is the data-space volume of shot-gathers defined at geophone
and source locations on the acquisition surface.  Equation~\ref{issrme1} is a
trace-by-trace operation to produce the multiple prediction with any
geophone-source, $({\bf x}_g,{\bf x}_s)$, combination by convolving each
trace of every shot gather with all the others followed by summing over the
convolution index ${\bf x}_a$.  Note however, the similarity of the
SRMP equation to the imaging condition of shot-profile migration,
equation~\ref{spic}.  
\par
Wave-equation extrapolation is performed on wavefields where data and
source-functions are used as initial conditions to propagate energy
into the subsurface.  To begin, traces at locations ${\bf x}_g$ are
inserted into a zero-valued wavefield defined along the axis ${\bf
  x}$.  Although data-space SRMP is a 
trace-by-trace operation, equation~\ref{issrme1} can be redefined in
terms of the wavefield $U({\bf x},{\bf x}_s)$.  With null-traces at
locations ${\bf x}$ where no data were collected, a multiple
prediction can also be written 
\begin{equation}\label{issrme2}
M_{z=0}({\bf x};{\bf x}_s,\omega)=\sum_{{\bf x}_a} 
U_{z=0}({\bf x};{\bf x}_a,\omega)
U_{z=0}({\bf x}_a;{\bf x}_s,\omega)\;,
\end{equation}
where the resultant volume has been regularized along ${\bf x}$ by
adding zero-traces, and I have added the specification that the
operation is being performed at the recording surface $z=0$.
\par
Using the principle of reciprocity between the receiver and source
locations (first and second arguments of the wavefields respectively),
the multiple prediction becomes
\begin{equation}\label{issrme3}
M_{z=0}({\bf x};{\bf x}_s,\omega)=\sum_{{\bf x}_a} 
U_{z=0}({\bf x};{\bf x}_a,\omega)
U_{z=0}({\bf x}_s;{\bf x}_a,\omega)\;.
\end{equation}
Here, the subscript $s$ on the RHS represents any different
receiver location (since it is the first argument of the wavefield),
and the dummy index ${\bf x}_a$ is recognized as a summation over
source location.  Therefore, using arbitrary index subscripts $c,d$
and restoring the significance of source location to subscript $s$
\begin{equation}\label{issrme4}
M_{z=0}({\bf x}_c;{\bf x}_d,\omega)=\sum_{{\bf x}_s} 
U_{z=0}({\bf x}_c;{\bf x}_s,\omega)
U_{z=0}({\bf x}_d;{\bf x}_s,\omega)\;.
\end{equation}
Finally, I define the dummy indices $c,d$ in terms of physically
significant variables location and half-offset,
${\bf x}_c={\bf x}+{\bf h}$ and ${\bf x}_d={\bf x}-{\bf h}$, such that
\begin{equation}\label{issrme}
M_{z=0}({\bf x},{\bf h},\omega)=
\sum_{{\bf x}_s} 
U_{z=0}({\bf x}+{\bf h};{\bf x}_s,\omega)
U_{z=0}({\bf x}-{\bf h};{\bf x}_s,\omega)\;.
\end{equation}
Thus reconfigured, equation~\ref{issrme} is now of parallel
construction to the shot-profile imaging condition,
equation~\ref{spic}, lacking only the summation over frequency.
\par
Extrapolation by $E({\bf x},\omega)$, in 
equation~\ref{extr}, simply redatums the
shot-gather $U$, after which the conventional SRMP volume is
calculated by convolution.  Image-space SRMP (IS-SRMP) is the application of a
second imaging condition evaluated at each subsurface depth level
during the migration that images only multiples.  It is the chain of
multiple prediction (convolution) and zero-time extraction (summation
over frequency).  The image-space multiple prediction, as a function of
sub-surface offset, is therefore
\begin{equation} \label{spmic}
m_z({\bf x},{\bf h})=
\sum_{{\bf x}_s}\sum_\omega 
U_z({\bf x}+{\bf h};{\bf x}_s,\omega)
U_z({\bf x}-{\bf h};{\bf x}_s,\omega)\;.
\end{equation}
\par
There are two important ramifications associated with
the equation for predicting multiples with the imaging condition 
above.  The first is that this operation is intrinsically a shot-domain
manipulation of the data.  After sorting to midpoint-offset coordinates, 
the source and receiver coordinates are mixed
in such a way as to make IS-SRMP more difficult for survey-sinking style
migration algorithms \cite{alison}.  Second, because reciprocity was
invoked to derive equation~\ref{spmic}, off-end (marine) acquisition
geometries will need to have split-spread gathers generated via
reciprocity.  The split-spread gathers will include the ray-paths from
multiples that emerge in front of the receiver spread (boat) which
need to be included in the shot-gathers to predict all possible
multiple events.
\par
Further understanding of the IS-SRMP imaging condition for those
familiar with shot-profile migration algorithms can be elicited by
defining the down-going wavefield in equation~2 as
$D\equiv U$.  Therefore, equations~1~-~3
become 
\begin{eqnarray}\label{extrU}
%%U_{z+1}({\bf k}_x;{\bf x}_s,\omega)&=&
%%U_{z}({\bf k}_x;{\bf x}_s,\omega) \mbox{exp}(+i\omega k_z)  \\
%%&\mbox{and} & \nonumber  \\ 
%%\hat{U}_{z+1}({\bf k}_x;{\bf x}_s,\omega)&=&
%%\hat{U}_{z}({\bf k}_x;{\bf x}_s,\omega) \mbox{exp}(-i\omega k_z) 
%%
U_{z+1}({\bf x};{\bf x}_s,t)&=&E^-({\bf x},t)
U_{z}({\bf x};{\bf x}_s,t) \\
&\mbox{and} & \nonumber  \\ 
\hat{U}_{z+1}({\bf x};{\bf x}_s,t)&=&E^+({\bf x},t)
\hat{U}_{z}({\bf x};{\bf x}_s,t) 
\end{eqnarray}
with the the imaging condition 
\begin{equation} \label{spmicswitch1}
i_z({\bf x},{\bf h})=
\sum_{{\bf x}_s}\sum_\omega 
 U_z({\bf x}+{\bf h};{\bf x}_s,\omega)
\hat{U}_z^*({\bf x}-{\bf h};{\bf x}_s,\omega)\;,
%=m_z({\bf x},{\bf h})\;,
\end{equation}
where $\hat{}$ denotes that after extrapolation in different
directions, the wavefields are no longer identical.  Because the 
conjugation of $D$ in the imaging condition
of equation~\ref{spic} can be commuted with the causality of the
extrapolation operator, it is not necessary to
extrapolate $U$ in two different directions.  Instead, the second
extrapolation above can be ignored and the imaging condition becomes
convolutional rather than correlational.  In this case,
equation~\ref{spmicswitch1} is exactly the multiple prediction,
equation~\ref{spmic}. 
\par
Cast in this manner, the
migration shows similarity with reverse-time migration
\cite{baysal:1514} and using multiples to migrate primaries
\cite{Shan.segab2.2004}.  The difference is that the data used as source
function is not first time-reversed.   IS-SRMP uses the shot gathers
as both an areal source function and data to image multiples.  Conversely,
time-reversing the data will use the primaries as areal source
functions (delayed from time zero) to image the subsurface with the
multiple reflections.  Using both an impulsive source function at
${\bf x}_s$ and time-reversed primaries as a source function 
would image Earth structure with the primaries and the multiples, but still
include the multiples in the image.

\subsection{Analytic example}
In 1D, let a trace be represented in the Fourier domain by the
expression
\begin{equation}\label{trace}
R(\omega)=
\e^{-i\phi_{w}}+
\e^{-i\phi_{e}}-
\e^{-i2\phi_{w}}-
2\e^{-i(\phi_{w}+\phi_{e})}-
\e^{-i2\phi_{e}}
= P + M_1. 
\end{equation}
The trace has primary reflections, $P$, at phase delays, $\phi$,
representing a water-bottom and a subsurface event (subscripts $w,e$).  Also included are
first order multiples, $M_1$, which are the water-bottom multiple at
$2\phi_{w}$, two peg-leg multiples at $\phi_{w}+\phi_{e}$, and the
event multiple at $2\phi_{e}$.  SRMP dictates the autoconvolution of
$P$ to derive $M_1$, which is clearly true.  Including the events
$M_1$ in the autoconvolution will predict the higher-order multiples
as well.  
\par
Extrapolating trace $R$ to a deeper depth applies a common phase
shift, say $\e^{-i\phi_z}$, to all terms in equation~\ref{trace}.
The trace then becomes
\begin{equation}\label{extrap}
R_z(\omega)=e^{-i\phi_z}(
\e^{-i\phi_{w}}+
\e^{-i\phi_{e}}+
\e^{-i2\phi_{w}}+
2\e^{-i(\phi_{w}+\phi_{e})}+
\e^{-i2\phi_{e}})
= e^{-i\phi_z}P + e^{-i\phi_z}M_1. 
\end{equation}
This equation shows that the extrapolation of data without multiple
subtraction produces the superposition of the redatumed primaries and
the redatumed multiples.  The extraction of the zero-time lag in the imaging
condition of migration states that energy in the wavefield should be
mapped into the image domain only when the extrapolation phasor
$\phi_z$, is equal to the time delay of the event in the data.  Thus
the water bottom primary is imaged when $\phi_z=\phi_w$, the water-bottom
multiple is imaged when $\phi_z=2\phi_w$, etc.   Whether the data $R$ is
first separated into its constituent parts, $P$ and $M_i$, or not, it
can be seen where in the image domain the various events in the above
example will be placed.  However, by first squaring the trace, 
implementing data-space SRMP, the water-bottom primary is
mapped into the image domain when $\phi_z=2\phi_w$.  This is the same
phasor that maps the recorded water-bottom multiple into the image
domain when migrating data contaminated with multiples.
\par
Figure~\ref{fig:sketch} is a cartoon depicting the 
generation of a 1D image-space surface-related multiple prediction
without first convolving the gathers.  The scenario drawn is for 
the simple case of a water-bottom reflector and its multiple.  The trace
$U$ denotes a zero-offset recording from a shot gather, and $D$ is the 
modeled source function used for shot-profile migration.  The trace 
$I$ represents the result of the conventional imaging condition, while
$M$ is the multiple prediction in the image space generated by 
autoconvolution.  The superscript $*$ denotes conjugation.  Both 
imaging conditions extract energy only at $t=0$.
%%\sideplot{GEO-2006-0199-R1-fig1}{width=1.8in,height=7in}{One-dimensional
\activesideplot{sketch}{width=1.8in,height=7in}{NR}{One-dimensional
  example of IS-SRMP during shot-profile migration at three
  extrapolation levels, $z=0,1,2$.  Trace $U$ represents data with a
  primary and multiple event.  Trace $D$ is the modeled impulsive
  source wavefield.  Trace $I$ is the conventional image ($UD^*$).
  Trace $M$ is the multiple prediction  ($UU$).  Final image volumes
  for each depth level are produced by extracting the values of $I$
  and $M$ at $t=0$.} 
\par
The first panel, $z=0$, shows the initial conditions of $U$ and $D$
and the fact that both $I|_{t=0}$ and $M|_{t=0}$ are devoid of 
events.  Since correlation subtracts the time to energy in the 
source trace $D$, currently zero, from the receiver trace, $U$, the 
image trace $I$ is the same as $U$.  The time for events on trace $UU$ 
are doubled from the initial condition, which results in the multiple 
on trace $U$ mapping below the time interval in the cartoon.
\par
The second panel depicts
the situation when the wavefields have been extrapolated to the 
depth of the water column.  Now energy in $D$ and $U$ is 
collocated and the primary is imaged during correlation, $UD^*$.  The
multiple prediction is still zero valued.  The last panel depicts
the wavefields and imaging results after the wavefields have been 
propagated to twice the depth of the water column.  Now the source
energy is collocated with the multiple, and the correlation results
in a negative polarity event.  Simultaneously, the primary is
now at $t=0$ which maps energy into the IS-SRMP volume, $M$.  However,
even though the kinematics of the multiple in $UU$ are correct, the 
event has the wrong polarity.  Also notice that the multiple is 
now at the original time of the primary, which sets up prediction 
of second-order multiples from the convolution of the first-order
multiples.  Importantly, even if the velocity used 
to calculate the extrapolation phasors is incorrect, both $I$ and $M$
will share the same error.  Also, $U$ consisting of
only primaries will correctly predict the location of multiples 
even if no multiples exist in the data.
\section{Examples}
\subsection{Simple synthetic}
The left panel of Figure~\ref{fig:syndat} is a synthetic shot
gather with two reflections and three multiples.  The velocity to the
first flat layer at depth $z=400$~m is 1500 m/s.  The velocity is 2500
m/s to the second flat layer at depth $z=1200$~m.  The traveltimes were
computed analytically and convolved with a wavelet.  The white
events are primaries, and the multiples have opposite polarity.  The
three multiple events are the simple multiples to both events, and the
asymmetric peg-leg.  The intrabed multiples between the two layers were
not included since they are not predicted by SRMP.  Receiver and source
spacing was 20~m.  The right panel is a subsurface-offset common-image 
gather produced by shot-profile migration of 350 shot-gathers identical 
to the data shown.  Because the layers are flat,
every location ${\bf x}$, away from the edges, is the same.  Only 
$1/6^{th}$ of the surface offsets are required to capture the 
moveout of the multiples in the subsurface-offset common image gathers
after migration. 
\par
%%\sideplot{GEO-2006-0199-R1-fig2}{width=2.25in}{Left:  Synthetic data
\activesideplot{syndat}{width=2.25in}{NR}{Left:  Synthetic data
  with two primaries and three multiples.  White events are 
  primaries.  Two simple and a peg-leg multiple were 
  modeled.  Right:  Subsurface-offset common-image gather produced 
  by shot-profile migration.}
The two panels in Figure~\ref{fig:synoffprim} are offset-domain 
common-image gathers produced by migrating only the primaries in the 
data above.  The first panel shows energy tightly focused at zero 
offset for the two reflectors.  The second panel is the multiple 
prediction.  Three clear events are imaged with concave-up 
moveout.  The middle event is the peg-leg multiple and the others are 
the simple multiples to the two reflectors.  At zero-offset, the events 
are analytically calculated to arrive at depths $z=1066,1866$ and 
$2666$~m, respectively. 
%%\sideplot{GEO-2006-0199-R1-fig3}{width=2in}{Offset-domain common-image gathers from
\activesideplot{synoffprim}{width=2in}{NR}{Offset-domain common-image gathers from
  migrating only the primaries in Figure 2.  Left: Conventional
  imaging condition.  Right: IS-SRMP.}
\par
The two panels in Figure~\ref{fig:synoffni} are
offset-domain common-image gathers produced by migrating all the
events in Figure~\ref{fig:syndat}.  Autoconvolution
predicts multiples from primaries.  Convolution of primaries with the
first-order multiples produces second-order multiples.  Autoconvolution 
of first-order multiples predicts third-order multiples.  The left panel 
now has the two focused primary events and the three multiples predicted 
by IS-SRMP in Figure~\ref{fig:synoffprim}.  These multiples 
are again predicted in the IS-SRMP gather on the right of this 
figure.  However, higher order multiples are now included as well.  The 
broad event at $z=1733$~m is the multiple that makes three trips through 
the water column.  The faint event at $z=2400$~m is the double
water-column multiple.  The last new event is at $z=2533$~m.  This
multiple has two trips in the water column and one trip from the
surface to the second reflector. 
%%\plot{GEO-2006-0199-R1-fig4}{width=2in}{Offset-domain common-image gathers from
\activesideplot{synoffni}{width=2in}{NR}{Offset-domain common-image gathers from
  migrating all events in Figure 2.  Left: Conventional 
  imaging condition.  Right: IS-SRMP with higher order multiples.}

\subsection{Sigsbee2{B} synthetic}
The Sigsbee2B data-set was designed to model strong surface-related
multiples from an offshore acquisition.  Two datasets were generated
with a 2D finite difference algorithm: One with the perfectly
reflecting free surface, and one
without\footnote{http://www.delphi.tudelft.nl/SMAART/S2Breadme.htm}.
Therefore, the direct subtraction of the two data volumes yields a
nearly perfect multiple model without the need for SRMP.  There are
slight differences in the source and receiver ghost effects between
the two data sets, so their subtraction retains faint residuals of
primary energy.  Though the data were modeled with an off-end
acquisition strategy, split-spread gathers were computed via
reciprocity for all of the examples below.  All of the images herein
were produced with four reference velocities in a SSF-PI shot-profile
migration code.  
\par
Figure~\ref{fig:sig1} shows zero subsurface-offset images of the Sigsbee2B 
data set.  The top panel, generated with the conventional imaging condition,
contains primaries and multiples.  The bottom panel has multiples 
and some migration artifacts.  The worst artifacts in the prediction are 
above the first multiple, especially near steep salt flanks.  This noise 
can be easily muted before subtraction.  The IS-SRMP image was produced using 
the data with primaries only.  The data was modeled such that the details of 
the bottom left corner have simple kinematic differences between primaries 
and multiples:  Events dipping up-right are multiples, and those dipping
down-right are primaries.
%%\plot{GEO-2006-0199-R1-fig5}{width=6.0in}{Migration of the Sigsbee2B data 
\activeplot{sig1}{height=7in}{NR}{Migration of the Sigsbee2B data 
containing multiples and the IS-SRMP result using data containing only 
primaries.}
\par
Figure~\ref{fig:mult} shows the bottom third of the
image produced with the Sigsbee2B data sets.  The top image used the
data modeled without the reflecting free surface and contains only
primaries.  The second image migrated the data with the 
free-surface and contains multiples.  The third panel is the
image produced by migrating the difference between the two modeled data
volumes (only multiples).  The complex multiples in this deep section
quickly overwhelm the primary events and could easily be mistaken for
primaries in some instances.  Notice the faint primaries in the top 
left corner and the basement reflector due to the imperfect subtraction.  The
bottom panel is a zoomed-in version of the multiple prediction in 
Figure~\ref{fig:sig1}.  No residual primaries are present, but some edge 
effects are visible at $z=6000$~m on the left side.  The bottom two panels 
are effectively identical which demonstrates the commutability of SRMP
convolutions and shot-profile migration.
%%\plot{GEO-2006-0199-R1-fig6}{width=6.0in}{Images of the bottom third
\activeplot{mult}{width=6.0in}{NR}{Images of the bottom third
  of Sigsbee2B modeled data.  Top to bottom:  Migration of primaries
  only.  Migration of multiples and primaries.  Migration 
  of multiples only.  IS-SRMP image using data containing 
  only primaries.}
\subsection{Gulf of {M}exico data}
The left panel of Figure~\ref{fig:gulfp} is a shot gather from a data 
set acquired
in 1997 by Western-Geco in the Mississippi Canyon lease area of the Gulf of 
Mexico.  Split-spread coverage was generated via reciprocity.  The data 
set contains 1096 shots and recorded to 10~s.  The source and receiver 
sampling was 27~m.  The right panel is a migrated angle gather from 
the same location.  The raw data contains 367 traces, while the angle gather 
contains only 60 traces.  Primaries are flat, and multiples are 
characteristically concave-down.  Before migration, the data were regularized, 
and bandpassed from 3-65~Hz.  The data were migrated with the same SSF-PI 
shot-profile migration algorithm with four reference velocities.  To produce 
only the zero-subsurface offset of the pre-stack migration, the cost of the 
imaging condition, including the 
multiple prediction, was 1.4\% of the cost of extrapolation.  Calculating
60 subsurface offsets approximately equalizes the costs of extrapolation 
(with four reference velocities) 
and imaging (including multiple prediction and conventional image 
generation).  This number of offsets is about $1/6^{th}$ of the traces in
the split-spread gather, and more than sufficient to capture the unfocused
multiple energy at all depths in the subsurface-offset gathers.
%%\plot{GEO-2006-0199-R1-fig7}{height=4in}{Left:  Shot gather from a Gulf 
\activeplot{gulfp}{height=4in}{NR}{Left:  Shot gather from a Gulf 
of Mexico data set.  Near-offset traces are null.  Split-spread
gathers were created via reciprocity.  Right:  Subsurface-offset 
gather from the migrated image at the same location.  Subsurface 
gather contains $1/6^{th}$ of the number of traces in the shot gather.}
\par
Figure~\ref{fig:gulf} shows the resulting migration and image-space
multiple prediction.  The multiples have the opposite polarity as 
expected.  The quantity of multiple events in the prediction is 
overwhelming.  Several 
events in the deep section of the image that have correspondence in the
multiple prediction, look very much like primaries.  The event that runs 
off the image at $z,x=6,2$~km could very easily be misinterpreted as a 
primary. 
%% Below the salt bottom and above the first water-bottom
%%multiple, several intrabed multiples can be identified in the image
%%that are not predicted in this context.
%%\plot{GEO-2006-0199-R1-fig8}{width=6in}{Mississippi Canyon, Gulf of 
\activeplot{gulf}{height=7in}{NR}{Mississippi Canyon, Gulf of 
Mexico zero-offset image and the IS-SRMP computed during the course
of shot-profile migration of 1096 shots.}

\section{Adaptive Subtraction}
After predicting the kinematics of the multiples, they must be
adaptively subtracted from the image due to the amplitude and
bandwidth problems associated with convolving (squaring) the
traces.  Another potential problem with generating a multiple
prediction by convolving data instead of only the primaries is the
prediction of higher-order multiples, even if they were not recorded
in the time interval of the original data.  This was noted
above in the IS-SRMP predictions with the flat-layer
synthetic, Figure~\ref{fig:synoffni}.  The adaptive
subtraction must be able to ignore events in the multiple prediction
that are completely absent from the data.  Also, the multiple
prediction has more artifacts than the image.  This is especially true
at shallow depths before the first multiple.  These artifacts must not
be introduced to the estimate of the primaries when removing the
multiples. 
\par
In this section I present subtraction results from the previous
predictions using adaptive subtraction.  Adaptive subtraction proceeds
by inverting for match filters connecting two similar objects, convolving
the filters with its object, and subtracting this result from the
second object.  This can be a complicated and subtle art that I make
no claims of performing at the optimum level.  All subtraction results
presented operated in two dimensions, and could likely be improved 
by incorporating more dimensions in the process \cite{antmult}.
\par
I pose the subtraction of the predicted multiples from the data as
the following linear inversion problem:
\begin{eqnarray}\label{adapt}
\mathbf{M}\mathbf{f}&\approx &\mathbf{d} \\
\epsilon\mathbf{A}\mathbf{f}&\approx&\mathbf{0},
\end{eqnarray}
where $\mathbf{M}$ is the convolutional matrix of the multiples (a matrix
whose columns contain shifted versions of a vector of the multiple 
prediction, $\mathbf{m}$), and
$\mathbf{f}$ is a bank of non-stationary filters acting on patches
\cite{GEE1} of the data vector, $\mathbf{d}$, which contains primaries
and multiples.  The matrix $\mathbf{A}$ is a regularization operator
(in this case a Laplacian) and $\epsilon$ controls the strength of 
regularization.  The
result of this linear inversion is a multiple model, $\mathbf{Mf}$, that 
matches the amplitude and wavelet of the data.

\subsection{Simple synthetic}
Figure~\ref{fig:m.ni} details the steps in the
subtraction process for an angle-domain common-image gather from the
flat-layer synthetic.  The first panel is the conventional image.  The
second panel is the IS-SRMP gather.  The last panel is the matched
version of the IS-SRMP gather, $\mathbf{Mf}$, ready to subtract from 
the image gather.  Some energy from the water-bottom primary has leaked into
the matched multiple prediction.  Where the first multiple crosses the
second primary, the filters have difficulty separating the two
events.  The higher-order multiples at the bottom of the initial prediction
have been removed.
%%\plot{GEO-2006-0199-R1-fig9}{width=3.0in}{Angle-domain image gathers from
\activeplot{m\DOT ni}{width=3.0in}{NR}{Angle-domain image gathers from
a flat-layer model.  Left to right: Conventional image, multiple 
prediction, prediction matched to data.}
\par
Figure~\ref{fig:p.ni} contains the final subtraction 
results.  The first panel is the subtraction of the first and third
panels in the previous figure, $\mathbf{d-Mf}$.  This is a reasonable result, 
but suffers from some artifacts.  The last panel is the original image
gather, $\mathbf{d}$.  The improved subtraction result in the center is produced
with a second application of the the match-filter technique described
in equation~\ref{adapt}.  For this application, the convolutional
matrix $\mathbf{M}$ is made from the subtraction result (first panel),
and the data vector $\mathbf{d}$ is once again the conventional image
gather.  The result shown is therefore the subtraction result matched 
to the input gather.  This two-stage matching methodology is a partial 
application of the much more thorough approach presented in 
\longcite{Alvarez.sep.125.gabriel1}.
%%\plot{GEO-2006-0199-R1-fig10}{width=3.0in}{Angle-domain image gathers from
\activeplot{p\DOT ni}{width=3.0in}{NR}{Angle-domain image gathers from
a flat-layer model.  Left to right:  Subtraction of matched multiple
prediction from data, subtraction matched to data, original data.}
\par
Missing data, either source or receiver positions, results in
incomplete multiple prediction for SRMP.  However, extrapolation
spreads energy across such gaps after a few propagation steps.  This
is commonly referred to as wave-front healing.  Collecting near-offset
traces is always difficult in the field.  Figures
\ref{fig:m.0n}~and~\ref{fig:p.0n} were produced to test the IS-SRMP
algorithm when data does not contain near-offset information.  Ten
null-traces surrounding the source were substituted in the gather
shown in Figure~\ref{fig:syndat}, which corresponds to a gap of
200~m.  There is some dimming near
zero angle in the primary events, and a small deviation in the
continuity of the curvature of the predicted multiples.
%%\plot{GEO-2006-0199-R1-fig11}{width=3.0in}{Angle-domain image gathers from
\activeplot{m\DOT 0n}{width=3.0in}{NR}{Angle-domain image gathers from
a flat-layer model with missing near offsets.  Left to right: Conventional 
image, multiple prediction, prediction matched to data.}
\par
In the interest of direct comparison, the adaptive subtraction 
parameters were kept constant for both Figures
\ref{fig:p.ni}~and~\ref{fig:p.0n}.  The faint energy in the center
panel of Figure~\ref{fig:p.0n} above the second reflector indicates
that the same parameters for the adaptive subtraction are not as
appropriate for this image volume.  However, the remainder of multiple
energy is easily removed from the final result when the subtraction is
better tuned.
%%\plot{GEO-2006-0199-R1-fig12}{width=3.0in}{Angle-domain image gathers from
\activeplot{p\DOT 0n}{width=3.0in}{NR}{Angle-domain image gathers from
a flat-layer model with missing near offsets.  Left to right:  Subtraction 
of matched multiple prediction from data, subtraction matched to data, 
original data.}
\subsection{Sigsbee2{B}}
The top panel in Figure~\ref{fig:sigsub} shows the bottom third of the 
multiple-contaminated image produced with the Sigsbee2B 
synthetic data.  The bottom image was produced with the data volume 
without multiples.  The center image is the result of predicting the 
multiples during shot-profile migration of data containing 
multiples and removing them from the top image.  The point diffractors 
masked in the top image are much better realized in the subtraction 
result.  The multiples in the sedimentary section before $x=10$~km are
almost perfectly removed except for some up-right dipping energy at 
$x=7$~km emerging from the basement reflector.  Several events between 
$x=10-19$~km are much more continuous and interpretable.  
%%\plot{GEO-2006-0199-R1-fig13}{width=6in}{Bottom third of the zero-offset image
\activeplot{sigsub}{height=7in}{NR}{Bottom third of the zero-offset image
produced with the Sigsbee2B synthetic.  Top to bottom:  Image with primaries 
and multiples, subtraction result, image produced from data without 
multiples.}

\subsection{Gulf of {M}exico data}
The top panel in Figure~\ref{fig:gsub} shows the bottom half of the 
multiple contaminated image produced with the Gulf of Mexico data.  The
center panel is the multiple prediction generated during shot-profile
migration.  The bottom panel is the subtraction result.  The subtraction
still contains some residual multiples associated with the 3D nature
of the salt body.  It is possible that the prediction could improve
with more diligent adaptive subtraction since the prediction captures
much of the detail of these multiples.  Many diffracted multiples
associated with the sharp tips of the salt have been removed around
$x=7$~and $20$~km.  Primaries at $z=4.5$~km on the left edge of the
image have been brought out and may suggest an anticlynal 
structure.  Most of the migration diffraction energy associated with
the rugose salt surfaces has also been attenuated.

%%\plot{GEO-2006-0199-R1-fig14}{width=6in}{Bottom half of the zero-offset image
\activeplot{gsub}{height=7in}{NR}{Bottom half of the zero-offset image
produced with the Mississippi Canyon, Gulf of Mexico, data provided by
Western-Geco.  Top to bottom:  Image with primaries 
and multiples, image-space multiple prediction, subtraction result.}

\section{Discussion}
In many geophysical applications, improved results are generated 
by performing operations in the pre-stack domain.  While more 
expensive, operating in higher dimensions can often produce better 
images.  Migration 
performs a sum over frequency, shot position, and (possibly) 
reflection angle to produce an image volume.  I suggest that 
the amplitude and bandwidth normalization performed in adaptive 
subtraction can be implemented at the imaging condition using 
deconvolutional variants \cite{lee:1786,EAGE-2006-P265}.  To 
normalize the frequency content and collapse the wavelet, the
imaging condition in equation~\ref{spic} becomes 
\begin{equation} \label{decspic}
i_z({\bf x},{\bf h})=
\sum_{{\bf x}_s}\sum_\omega    \frac{
  U_z({\bf x}+{\bf h};{\bf x}_s,\omega)
D_z^*({\bf x}-{\bf h};{\bf x}_s,\omega) } 
{\ll 
D_z ({\bf x}-{\bf h};{\bf x}_s,\omega)
D_z^*({\bf x}-{\bf h};{\bf x}_s,\omega)
\gg}\;,
\end{equation}
where the denominator is smoothed across horizontal coordinates
(angle braces) for stability.  Similarly, the IS-SRMP,
equation~\ref{spmic}, becomes 
\begin{equation} \label{decspmic}
m_z({\bf x},{\bf h})=
\sum_{{\bf x}_s}\sum_\omega    \frac{
U_z  ({\bf x}+{\bf h};{\bf x}_s,\omega)
U_z  ({\bf x}-{\bf h};{\bf x}_s,\omega) }  {\ll 
U_z  ({\bf x}-{\bf h};{\bf x}_s,\omega)
U_z^*({\bf x}-{\bf h};{\bf x}_s,\omega)\gg }\;.
\end{equation}
\par
Our experiments to date have not shown this approach particularly 
effective.  There are minor improvements in both zero-offset images, 
but adaptive subtraction of the multiple prediction is still 
required.  The stability of the 
multi-offset images is very sensitive to the smoothing parameters 
selected.  The deconvolutional variants are substantially more expensive,
and have not, so far, proven to be worth the computational 
burden.
\par
While IS-SRMP, by definition, produces only surface-related multiples,
the technique could be manipulated to address strong multiple
generators in the subsurface.  A layer-stripping type approach 
derived in parallel to the DSR multiple prediction \cite{alison} 
could probably be implemented at the cost of extrapolating a third 
(in addition to the up-coming and down-going) wavefield.
\section{Conclusion}
Image space SRMP produces a multiple prediction by convolving the data
with itself at every subsurface depth level during the course of a
shot-profile migration.  The result is the same as migrating
the conventional data-space multiple prediction (SRMP).  This method 
is most convenient with shot-profile migration strategies 
since the convolution operation must operate on distinct gathers rather 
than the combinations thereof produced by resorting to CMP 
coordinates.  The simplicity of this
approach can immediately be leveraged to generate the image-space
multiple prediction directly from any shot-profile migration 
program.  The method is immediately applicable to 3D, and non-zero
subsurface-offset and angle.  Any migration algorithm that maintains
separate up-coming and down-going wavefields and uses a combinatory
imaging condition (e.g. plane-wave and reverse-time migrations) can
be easily modified to produce an IS-SRMP volume. 
\par
Given a reasonably accurate velocity 
model for migration, it is only necessary to compute $O(10)$ subsurface
offsets.  This results in many fewer traces involved in calculating
the multiple prediction than the $O(1000)$ offsets collected at  
the surface.  This savings will be reduced however by the need to 
convolve the traces at every depth level, $O(100)$, of the image space rather
than just at the surface.  Whatever the balance of floating-point 
operations for a particular survey, the convenience of being
able to generate the multiple prediction during another required 
processing step can save file manipulation, sorting, and overhead 
costs.  Furthermore, this technique can also be implemented in a target 
oriented fashion simply by not calculating the multiple prediction 
at shallow depths where it is not required.
\par
Importantly, split-spread gathers must be pre-computed via reciprocity
for data collected with off-end acquisition geometries.  Off-end
gathers will not contain (nor therefore predict) emerging 
rays which pierce the acquisition surface in front of the boat.  This may 
increase the size of the computational domain used for propagating 
each individual shot-record.  The cost increase by performing two
imaging conditions is not severe, as the cost of calculating an
imaging condition with in-line offsets is usually a fraction of
the cost of a shot-profile migration.  Therefore, whenever a shot-profile
migration is being performed, it may be advantageous to generate the
IS-SRMP even if a data-space elimination effort has already been
performed, especially when the adaptive subtraction leaves some
multiple energy in the final result.
\par
The quality of the multiple prediction produced in the image space with this
technique is independent of the accuracy of the velocity model used
during the migration.  The multiple prediction, propagated with the same
velocity model, must be kinematically accurate with the location of
any multiples in the migrated image.  Though both the image and the
multiple prediction may be mis-migrated, events in the volumes
always correspond when constructed with the same extrapolation
operators and velocity model.   
\par
Due to the squaring of the wavelet when convolving the data, the multiple 
prediction can not be directly subtracted from the data or the 
image.  The adaptive subtraction results above show high quality estimations
of the primaries for both the Sigsbee2B synthetic and the Gulf of Mexico
data examples.  Substantial improvement to the estimation of primaries
was achieved by re-using the match filter inversion process to match 
the subtraction result to the original image.  This second iteration 
mitigates some of the residual energy from the subtraction, 
restores some of the primary energy inadvertently removed during 
subtraction, and attenuates some of the artifacts associated with migration
and the transform from subsurface offset to angle domain common-image
gathers.

\section{Acknowledgments}
This chapter was written for publication in Geophysics (accepted 2006) 
with two coauthors.  Gabriel Alvarez provided expertise with the
adaptive subtraction efforts.  Ken Matson provided valuable insight into 
the world of multiple prediction (with which I have little to no 
experience).  During the review process for publication, reviewers 
C.P.A Wapenaar, D.J. Verschuur, and Paul Sava all made 
insightful comments to drafts of the manuscript that substantially 
improved its quality.  
\clearpage
\label{.end}
\clearpage

\end{document}
