\lefthead{Cunha \& Muir}
\righthead{Separation of converted modes}
\footer{SEP-61}
\title{Separation of converted modes in marine data}
\author{Carlos A. Cunha Filho \& Francis Muir}
\ABS{ Two methods are tested to separate
the converted modes from the dominant P wavefield. The first one  
uses the critical slowness for P waves as the discriminant 
factor and is ineffective for deep reflectors.
A better result is obtained with another approach, 
in which the upcoming wavefield is divided into
zones sharing a common range of ray parameter; a filtering process
using the horizontal slowness as the discriminant factor is then applied
to each zone. Although more efficient, this 
method is more sensitive to contamination by multiples.
For this reason the suppression of multiples and peg-legs 
prior to the separation process becomes of vital importance.}
\INT
Wave propagation in a marine environment can be separated into two
parts: the propagation through the water, in which the acoustic wave equation 
applies, and the propagation in the subsurface strata, in which
the elastic wave equation applies.
Propagation in the water is described by the scalar pressure
wavefield recorded by the hydrophones; in solid layers 
the displacement vector, or 
elastic field separation should be used.
The coupling of the two potentials 
can be achieved by the application of the continuity of the normal
 components of stress and displacement at the liquid-solid interface.

The conversion of the scalar recorded field into the elastic
wavefield at the ocean floor and the subsequent separation of
the converted wavefield can be achieved by 
the application of successive steps.

First the recorded pressure field is decomposed into the upcoming and 
downgoing wavefields (Claerbout, 1971), and since the downgoing wave
at the cable is just the upcoming wave reflected at the water surface
(or ghost), it is possible to express the upcoming wavefield as a
function of the recorded wavefield.
Then, the principle of survey
sinking (Schultz and Sherwood, 1980; Claerbout, 1985) is used to 
extrapolate the upcoming pressure field down to the ocean floor.
The next step is the conversion of the scalar pressure field
into the displacement amplitude field of the P waves in the water;  
this conversion is accomplished through a simple relation in the $w-k_x$ domain.
Finally, a slowness filter is used to obtain the displacement 
amplitudes of P and SV waves below the water bottom interface. 
This filtering process shows to be more effective if a prior
separation of the data into \em P wave Snell rays domains \rm is
performed, so that a different filter can be used for each domain.

\mhead{UPCOMING PRESSURE FIELD AT SEA FLOOR}

As shown in Figure~\ref{fig1}, the recorded pressure field can be 
decomposed into the downgoing  and upcoming wavefields at the cable
depth
\begin{equation}
  \phi(x,z,t) = \acute{\phi}(x,z,t) + \grave{\phi}(x,z,t) ,
\end{equation}
\vspace{0.3cm} and, assuming a perfect reflection at the water's surface,
$$ \phi(x,z=z_0,t) =\acute{\phi}(x,z=z_0,t) - \acute{\phi}(x,z=-z_0,t) .$$

\plot{fig1}{3.5in}{/r0/carlos/data/vplottex/PUD}
{U/D decomposition of the recorded pressure wavefield}
\vspace{0.5cm}
Let $ {\cal L}(z) $ be the extrapolation operator for upcoming waves

$$ \acute{\phi}(x,z=-z_0,t) = {\cal L}(-2z_0) \acute{\phi}(x,z=z_0,t) , $$
then

$$\acute{\phi}(x,z=z_0,t)= [1 - {\cal L}(-2z_0) ]^{-1} \phi(x,z=z_0,t) . $$
\newline
and since \hspace{1.5cm} 
$ \acute{\phi}(x,z,t)= {\cal L}(z-z_0) \acute{\phi}(x,z=z_0,t)$ \hspace{1.5cm}
we get

\begin{equation}
\acute{\phi}(x,z,t) = {\cal H}(z,z0) \: \phi(x,z=z_0,t) ,
\label{eq:phi}
\end{equation}
\vspace{0.5cm}
where
\begin{equation}
{\cal H}(z,z_0) = {\cal L}(z-z_0) [ 1 - {\cal L}(-2z_0) ]^{-1} .
\label{eq:Hzz0}
\end{equation}

\vspace{0.5cm}
The operator defined in equation~(\ref{eq:Hzz0}) performs 
two distinct processes: the conversion
of the recorded pressure field into the upcoming pressure field, and its
downward continuation to the sea bottom.

In the $ \omega$-$ k_x $ domain, the downward extrapolation operator
for upcoming waves can be expressed by

 $$ {\cal L}(z) = e^{i k_z z} = e^{-i \sqrt{{\omega^2 \over v^2}-k_x^2}
\: {\textstyle z}} , $$
\vspace{0.5cm} which can be substituted into the Fourier transform 
of equation (\ref{eq:phi}) to obtain

\begin{equation}
\acute{\Phi}(k_x,z,\omega) = {\cal H} \; \Phi(k_x,z=z_0,\omega) ,
\end{equation}
\vspace{0.5cm}
where
\begin{equation}
{\cal H}(z,z_0,k_x,\omega)={e^{i k_z z} \over 2 \; i \; \sin{k_z z_0}} .
\label{eq:Hkxw}
\end{equation}

It is evident from equations~(\ref{eq:Hzz0})~and~(\ref{eq:Hkxw})
 that the operator $ {\cal H}$ has poles
at $k_z \: z_0 = n \pi $, corresponding to vertical wavelengths of
$ 2 \: z_0 \over n$ .

For waves with these wavelengths, the downgoing and upcoming fields will
cancel each other at the cable depth (assuming perfect reflection 
at the water's surface), and the result is a null recorded wavefield. Of course
no information can be obtained for the upcoming field in this case, but
since it represents only a couple of lines \footnote{For $ \triangle t=4$ms
and $z_0=10$m, only n=0 and n=1.} in the $\omega$-$k_x$ plane
we can neglect its contribution for the final, downward continued 
wavefield.

\mhead{SEPARATION OF P AND S WAVEFIELDS}

To apply scattering theory we need to transform the pressure
field into a displacement field. The pressure wavefield obeys
the scalar wave equation

$$\nabla^2 \phi(x,z,\omega)= -{\omega^2 \over v^2} \phi(x,z,\omega) . $$

The P wave displacement $ \vec{p} $ is related to the pressure by the
following expression:

$$ \phi = - K \: \nabla \! \cdot \! \vec{p} , 
\mbox{ \hspace{1.2cm} where \hspace{0.6cm} $K$ \hspace{0.6cm} is
the bulk modulus of water. } $$

Recalling that the medium (water) is homogeneous and that the displacement 
field is irrotational, we can use the relation 
$ \nabla^2 = \nabla \cdot \nabla $ 
and combine the two previous equations to get

$$ \vec{p} =  {v^2 \over \omega^2 K} \nabla \phi .$$

This  equation relates the displacement vector field to the scalar
pressure field and can be rewritten in the following form:

$$ \pmatrix{p_x \cr p_z }
=  {v^2 \over \omega^2 K} \pmatrix{{\partial \over \partial x}
\cr {\partial \over \partial z}} \phi(x,z,\omega) . $$

And after a double spatial Fourier transform, it will be expressed by 

$$ \pmatrix{P_x(k_x,k_z,\omega) \cr P_z(k_x,k_z,\omega) } =
i { v^2 \over \omega^2 K} \pmatrix{ k_x \cr k_z} \: \Phi(k_x,k_z,\omega) . $$

Finally, the displacement amplitude $\acute{P}$ of the upcoming P wave 
on the water can be expressed as a function of the upcoming pressure
field:

$$ \acute{P} = i {v^2 \over \omega^2 K} \sqrt{k_x^2 + k_z^2} \: \acute{\Phi} $$
or

$$ \acute{P} =  {v \over -i \omega K} \acute{\Phi} .$$
\newline
and, as one should expect, the pressure field differs from the velocity field 
(time derivative of the displacement) by just a scale factor 
$ v \over K $ .

\plot{scatter}{2.5in}{/r0/carlos/data/vplottex/SCAT}
{Scattering of the upcoming and downgoing wavefields at the ocean floor}

Figure~\ref{scatter} shows how the upcoming P and S waves are scattered
at the ocean bottom interface. The displacement amplitude of the upcoming 
pressure field just above the sea floor can be expressed as functions 
of the displacement amplitudes of all the incident waves:
\begin{equation}
 \acute{P_1} = [ \acute{P} \acute{P} ]_{2-1} \acute{P_2} +
[ \acute{S} \acute{P} ]_{2-1} \acute{S_2} +
[ \grave{P} \acute{P} ]_{1-1} \grave{P_1} ,      \label{eq:SPupPup}
\label{eq:Scatt}
\end{equation}
where $ [ \acute{P} \acute{P} ]_{2-1} $ is the transmission coefficient
for upcoming P into P waves, $ [ \acute{S} \acute{P} ]_{2-1} $
is the transmission coefficient for upcoming S into P waves
, and $ [ \grave{P} \acute{P} ]_{1-1} $
is the reflection coefficient for downgoing P into 
P waves (Aki and Richards, 1980).
The last term in the previous equation corresponds to the ocean floor
primary reflection and all its associated multiples and peg-legs. 
The conditions most suitable for the generation of converted waves in marine 
surveys are the same ones that are responsible for the occurrence of
strong multiple contamination. While the ocean floor 
primary will not interfere with any converted wave, the multiples
and peg-legs must be removed in order to isolate the PSSP wavefield.
The importance of this step and the method used for its accomplishment
are the subjects of Appendix A.

\plot{teta1}{2.0in}{/r0/carlos/data/vplottex/teta1}
{The relative recorded amplitudes of the several modes as a function
of the angle of incidence, for a specific model in which velocity
increases with depth ($ V_p = 3000 $ m/s for the bottom layer).}

The two remaining terms in equation~(\ref{eq:Scatt}) 
(which correspond to the P and S wavefields)
can be reasonably separated by the use of the critical angle for P waves at
the first layer as a discriminant factor (Cunha, 1989), 
as illustrated on Figure~\ref{teta1}.
Waves of the type PPPP, PPSP and PSPP will be found only with horizontal 
slownesses smaller than the slowness corresponding to the ocean floor 
critical angle for P waves (30 degrees on the figure),  
while only the PSSP mode will be found at higher horizontal slownesses. 
As a first approach we consider that only the dominant mode (PPPP)
is present for pre-critical horizontal slownesses.

 Under these assumptions, a ``P-sensible receiver" coupled 
at the sea bottom would record a displacement

\begin{equation}
\acute{P_2} ={1 \over [ \acute{P} \acute{P} ]_{2-1}} \acute{P_1}
\mbox{\hspace{1.5cm}for \hspace{1.5cm}$p = {k_x \over \omega} < {1 \over V_p}$}
 ,
\label{eq:Pup}
\end{equation}
\newline
\vspace{0.5cm} whereas an ``S-sensible receiver" would record
\begin{equation}
\acute{S_2} ={1 \over [ \acute{S} \acute{P} ]_{2-1}} \acute{P_1}
\mbox{\hspace{1.5cm}for \hspace{1.5cm}$p={k_x \over \omega} \geq {1 \over V_p}$}
 ,
\label{eq:Sup}
\end{equation}
\newline
\vspace{1.0cm} where $ V_p $ is the velocity of P waves at the first layer.


\mhead{DIVIDING THE DATA INTO SNELL ZONES}

Although the use of the critical horizontal slowness to perform the separation
between P and S waves is effective for shallow reflectors, the
finite size of the cable restricts its application to deep reflectors 
because a converted wave with horizontal slowness larger than the P wave
critical slowness will be received at very large offsets.

A possible way to overcome the limitations  associated with
the finiteness of the field aperture is to use a criterion  
that somehow takes into account the depth of the reflector.
Transforming the data into the \em $ \tau $-$ p $ \rm domain provides
the necessary flexibility for use of a variable slowness cutoff. 
Tatham and Goolsbee (1984) showed also that better
results can be achieved when a hyperbolic velocity filtering is used
during the $ \tau $-$ p $ transform to limit the range
of reasonable stacking velocities.

As we will see, the separation of the data into ranges of \em Snell rays \rm
would seem to be another appropriate way to perform the
slowness filtering with a variable-slowness cutoff.

\plot{snellray}{3.6in}{/r0/carlos/data/vplottex/SNELLRAY}
{\em (a) \rm A ray is propagating and has been reflected at each interface.
 \em (b) \rm The reflections associated with the same ray in the 
\em x-t \rm domain.}
	
A Snell ray (Ottolini, 1982) can be defined for a plane-layered
earth as a ray that keeps a constant horizontal slowness (obeys the Snell
law) while it propagates through the subsurface,  
as illustrated in Figure~\ref{snellray}. The reflection
points corresponding to a Snell ray with ray parameter (or horizontal
slowness) $p$ are defined by 

$$ {\partial x \over \partial t} = p v^2(z) , 
\mbox{ \hspace{1.2cm} where \hspace{0.6cm} $v(z)$ \hspace{0.6cm} is
the interval velocity at depth $z$. } $$

\sideplot{rayplot}{3.6in}{/r0/carlos/data/vplottex/rayplot}
{ Separation of the data into regions limited by a
set of Snell trajectories (reflection points associated with a Snell ray).}
	
For the simple case of a plane-layered earth, the horizontal slowness of a
P wave will always be lower than the horizontal slowness of a converted
wave recorded at the same position in the \em x-t \rm domain.
It is possible then to choose a set of ray parameters and 
divide the data into regions whose P waves Snell rays are
limited by two adjacent values (Figure~\ref{rayplot}).
Since the converted waves inside each region will have higher horizontal 
slownesses than the P waves, a different filter can be applied
to each region; the cutoff value is given by the ray parameter
of the Snell ray that defines the upper limit of that region.


\mhead{SYNTHETIC MODEL AND REALITY}

\plot{model}{7.5in}{/r0/carlos/data/vplottex/sldGM0}
{\em (a) \rm Elastic model. All the layers (except water) have 
density equal to 2 g/cc and Vp/Vs ratio =1.83. 
\em (b) \rm Elastic synthetic seismogram generated by the model in \em a \rm .}

The method has been applied on a 1D elastic model, generated
by a program that uses Haskel-Thompson propagation matrices.
Figure~\ref{model} shows the original synthetic seismogram
and the model, which contains a total of seven layers overlaid by water,
 each of them with constant density and Poisson's ratio. The
group interval is 50 m with 48 channels and a maximum offset
of 2650 m.

\plot{raytrace}{7.5in}{/r0/carlos/data/vplottex/raytrace}
{\em (a) \rm Ray-trace synthetic with primary arrival of PSPP-PPSP and the 
sea bottom \em (b) \rm  The same as \em a \rm but for PSSP mode.
The amplitude variations were not taken into account. }

\plot{sldGMP}{7.5in}{/r0/carlos/data/vplottex/sldGMP}
{\em (a) \rm Seismogram of an acoustic equivalent model \em (b) \rm
Seismogram corresponding to a full space (air replaced by water). The
ghost, multiples, and peg-legs associated with the
water layer are absent in \em b \rm .}

To facilitate the identification of each of the several superposed 
events, four other auxiliary synthetics were also generated. 
Two of them are simple ray-tracing synthetics for PSPP-PPSP primary
arrivals (Figure~\ref{raytrace}{\em a \rm}) and for PSSP primary
arrivals (Figure~\ref{raytrace}{\em b \rm}); the sea-bottom reflection is
also included in both seismograms. The other two correspond to an 
acoustical equivalent seismogram using the same P wave 
velocities as the elastic model (Figure~\ref{sldGMP}{\em a \rm}),  
and a seismogram free of \bf ``water-layer multiples" 
\rm (Figure~\ref{sldGMP}{\em b \rm}), where air was replaced by water.

All of the three most basic converted modes (PPSP, PPSP and PSPP) 
are present on Figure~\ref{sldGMP}{\em b \rm} (without the interference
of multiples), but not on Figure~\ref{sldGMP}{\em a\rm}; for example,  
the events that reach the farthest offset of 
Figure~\ref{sldGMP}{\em b \rm} at 2.4s (PSPP) and 3.0s (PSSP) or the nearest
offset at 0.9s (PSPP), 1.05(PSSP) and 1.8(PSPP) and are all absent
on Figure~\ref{sldGMP}{\em a \rm} and hard to distinguish 
on Figure~\ref{model}{\em b \rm}.

\plot{sldGMrev}{7.5in}{/r0/carlos/data/vplottex/sldGMrev}
{\em (a) \rm Seismogram of the synthetic model \em .(b) \rm 
The same seismogram after application of the multiple
suppression algorithm. }

A comparison between the model before and after 
the application of the multiples-suppression algorithm 
is showed in Figures~\ref{sldGMrev}{\em a \rm} and \ref{sldGMrev}{\em b. \rm} 
It is evident that the method fails to remove the multiples in the near traces. 
The reasons for that result are discussed along with the description of 
the method in Appendix A. Whereas the suppression's efficiency at 
the small offsets is not relevant for the present goals, the
algorithm's performance at large offsets is crucial for 
the isolation of the PSSP wavefield. Several converted-waves
arrivals are more evident after the multiple's attenuation, 
such as the reflections that cross the farthest offset
at 3.0s and 4.05s on section \em b \rm of 
Figures~\ref{sldGMrev}~and~\ref{sldGMP} and are absent
on section \em a \rm of both figures. Meanwhile, some
further refinement is still required for the multiples-removal 
process at large offsets, as it can be seen by
the presence of a PPSP water-bottom multiple at 2.7s 
in the last traces of the deconvolved data.

\plot{sldGMsplit}{7.5in}{/r0/carlos/data/vplottex/sldGMsplit}
{\em (a) \rm  Seismogram of the synthetic model after application of a
constant slowness filter and \em (b) \rm after the application
of a variable slowness filter. }

When a constant horizontal-slowness filter is used to perform the
S wavefield separation, the result exhibits an ``artificial appearance" 
due to the restricted slowness range accepted
by it. In contrast, the splitting method based on the separation
of the data into different Snell domains preserves a most
suitable range of stepouts in the S wavefield section. 
Figures~\ref{sldGMsplit}{\em a \rm} and \ref{sldGMsplit}{\em b \rm} 
provide a good comparison between the two methods. All the events are
doubtless more clearly discernible  in \em b \rm 
than in \em a \rm. A correlation between Figure~\ref{raytrace}
(which contains the expected positions of the converted
primary reflections) and Figure~\ref{sldGMsplit}{\em b \rm} (which
corresponds to the separated S field) shows that, although
most of the desired events are present, some multiples
associated with converted waves are also visible, such as
the ones that intercept the last trace at 3.5s and 4.5s.
These events correspond to multiples that were not
satisfactorily eliminated in the suppression process.

Since multiples and primaries have the same horizontal slownesses,
 the multiples of P waves will be also eliminated 
during the filtering process when  
a constant slowness cutoff is used. However, when the data 
is divided into different ranges of horizontal slowness so that
a variable filter can be applied, the multiples and primaries
inside a region will have different horizontal slownesses and
the contamination with multiples will became critical.

\plot{sldGMvel}{7.5in}{/r0/carlos/data/vplottex/sldGMvel}
{\em (a) \rm Velocity panel of the original model. \em (b) \rm 
Velocity panel after multiples-suppression. \em (c) \rm 
Velocity panel corresponding to the constant slowness filter. 
\em (d) \rm Velocity panel corresponding to the variable
slowness filter. }

An overall idea of the differences between the two methods
of mode separation and the contribution of the multiples-suppression
can be achieved if we contrast the velocity
analysis panels regarding each step 
(Figure~\ref{sldGMvel}{\em a \rm}-{\em d \rm}).
The panel in \em a \rm refers to the original model while
the velocity analysis in \em b \rm corresponds to the
output of the multiples-suppression process. The differences
of interest appear only on the shallow events, like the one
with velocity 1350m/s at 0.75s that is present only on 
\em b \rm or the one with velocity 1300m/s at 2.2s that
is stronger in \em b \rm. In addition, the velocity panel
corresponding to the constant-slowness splitting is
displayed in \em c \rm while the last panel (\em d \rm)
refers to the variable-slowness filtering applied over
different Snell zones. 
It is important to notice that whereas the shallow events
are easy to identify in \em c \rm and most of the
energy associated with multiples has been suppressed, 
an undesirable trend corresponding to aliased energy strongly
contaminates the whole panel. The use of a variable
ray parameter filter, however, provides a clean panel
with a definite trend of converted-waves' stacking velocities.

\plot{GMft25}{7.5in}{/r0/carlos/data/vplottex/GMft25}
{\em (a) \rm  Shot gather of a Barents Sea survey. \em (b) \rm  The same gather
 after the separation process. \em (c) \rm Velocity analysis
of \em a. (d) \rm Velocity analysis of \em b \rm. }

The same procedures were applied on data recorded by GECO in
Barents Sea (offshore Norway). The data is composed of 48
channels with a group interval of 50 m and
maximum offset of 2650 m. The water depth is close to
300 m and the ocean floor is ``semi-hard" (P velocity
around 1850 m/s). A shot gather of this data is shown in 
Figure~\ref{GMft25}{\em a \rm} and the same gather, after 
multiples-suppression and separation of the converted wavefield by
the variable slowness algorithm, is displayed in \em b \rm.
The velocity panels corresponding to the two seismograms
are shown in \em c \rm and \em d \rm. Although a 
nonambiguous interpretation of the velocity trend is not 
feasible, it is possible to delimit with reasonable
confidence the range of acceptable stacking velocities
for the converted waves. 

The application of the method in a more appropriate dataset
is still required, because the poor quality of the present
data and the inadequacy of the survey parameters (maximum
offset and group interval) strongly restrict the resolution
of any prestack procedure and compromise a critical
evaluation of its efficiency.


\mhead{CONCLUSION}

The critical point in the attainment of an elastic wavefield
in marine environments is the process of separation of
the several modes that cohabit in the recorded data.
Particularly important is the step regarding the 
deconvolution of water-bottom multiples
and peg-legs, that will share a common
range of horizontal slownesses with the converted
modes associated with deep reflectors.
The use of a constant slowness filter in the \em $ \omega $-$ k_x $
\rm domain proves to be ineffective for deep reflectors, for which
even the converted modes have a small ray parameter value.
The usual solution is to perform the splitting process in
the \em $ \tau $-$ p $ \rm domain, in which an appropriate filter
can be designed for each particular case. In this paper 
an alternative method was
introduced: it makes use of Snell rays to delimit zones of
different P wave horizontal slownesses in the data, and 
applies a different filter to each region. 
Although a better treatment of the multiples is
still required, the method shows to be effective even 
in a conventional survey.
The next step is the application of the whole
process in a more appropriate dataset (long cable,
small group interval, etc.). 


\mhead{ACKNOWLEDGMENTS}

I would like to thank Dave Nichols for the discussions
and suggestions that contributed to this work. 

\mhead{REFERENCES}

\reference{Aki, K.I., and Richards, P.G., 1980, Quantitative seismology: W. H. Freeman and Co.}

\reference{Cunha Filho, C.A., 1989, Energy partition on marine converted waves: SEP-{\bf 60}, 243-252.}

\reference{Claerbout, J.F., 1971, Toward a unified theory of reflector mapping: Geophysics, {\bf 36}, 467-481.}

\reference{Claerbout, J.F., 1985, Imaging the Earth's Interior: Blackwell Scientific Pub. Inc.}

\reference{Morley, L., 1982, Predictive techniques for marine multiple suppression: Ph.D. thesis, Stanford University, also SEP-{\bf29}.}

\reference{Ottolini, R.A., 1982, Angle-midpoint migration: Ph.D. thesis, Stanford University, also SEP-{\bf33}.}

\reference{Schultz, P.S. and Sherwood, J.W.C., 1980, Depth migration before stack: Geophysics, {\bf 45}, 376-393.}

\reference{Tatham, R. H., and Goolsbee, D. V., 1984, Separation of S-wave and P-wave reflections offshore western Florida: Geophysics, {\bf 49}, 493-508.}

\reference{Wiggins, J. W., 1988, Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction: Geophysics, {\bf 53}, 1527-1539.}
%

\mhead{APPENDIX A}
As stated before the conditions most suitable for the generation of
converted waves are also appropriate for the generation of
water-bottom multiples and peg-legs. This appendix 
describes the method used to attenuate this kind of
multiples before the separation process is applied.

\shead{Attenuation of Multiples}
The third term on equation~(\ref{eq:SPupPup}) corresponds to the multiples
and peg-legs associated with the water bottom and must be removed
so that equations~(\ref{eq:Pup})~and~(\ref{eq:Sup}) can be used to separate
the two upcoming waves. The process used is a simplified version of 
the method developed by Morley (1982) and described by Wiggins (1988) 
for the case of a nearly horizontal ocean floor.

\plot{multmethod}{4.5in}{/r0/carlos/data/vplottex/MULTMETH}
{For an horizontal ocean floor, 
the extrapolation of the upcoming wave by twice the water depth 
is equivalent to the actual propagation that generates
water-bottom multiples, except
for an angle-dependent factor.}

As illustrated in Figure~\ref{multmethod}, the ocean-bottom multiples 
associated with the upcoming wavefield can be predicted using the following
process:
\begin{enumerate}
\item
Replacement of the air layer by a water layer, above the water surface.

\item
Upward continuation of the wavefield by a distance of two water depths.

\item
Multiplication by a factor that corresponds to the product 
of the reflection coefficients of the water bottom and water surface.

\end{enumerate}

Similarly, the downgoing wavefield (or ghost), can have its multiples
predicted by the same process, except that the wavefield must be 
continued downward instead of upward.
Since the operators to downward-continue a downgoing wave or
upward-continue an upcoming wave in a homogeneous media are identical,
 the predicted wavefields corresponding to the multiples can be
expressed by

$$ m(x,t) = \alpha(p) {\cal L}(x,t,z) d(x,t) ,$$
\newline
where $ \alpha(p)$ is the product of the two mentioned reflection
coefficients at horizontal slowness \em $p$, \rm and \em $ z $ \rm 
is the water depth.
The determination of \em $ \alpha(p)$ \rm and \em $ z $ \rm is achieved
by an optimization algorithm that minimizes the objective function 

$$ \delta(\alpha_{(p)},z) = \sum_x \sum_t \parallel d(x,t) -
m(\alpha,z,x,t) \parallel^2 .$$

This optimization is easily performed in the \em $\omega$-$k_x$ \rm
domain because the extrapolation operator has the form of a simple
phase shift, and the decomposition into different values of slowness
is accomplished by simple radial slices within the transform plane. 
The objective function on this domain assumes the form

$$ \delta(\alpha_{(p)},z) = \sum_k \sum_{\omega} \parallel d(k,\omega) -
m(\alpha,z,k,\omega) \parallel^2 .$$

\plot{alpha}{3.0in}{/r0/carlos/data/vplottex/alpha}
{ The continuous line corresponds to $ \alpha(p) $, the water-bottom 
reflection coefficient obtained by the optimization
process, while the dashed line is the theoretical curve
for the same model. }

The comparison between the computed curve \em $ \alpha(p)$ \rm 
and the theoretical curve for the ocean-floor reflection 
coefficient (Figure~\ref{alpha}) shows that the computed reflection coefficient
diverges from the theoretical at small and large values
of slowness. The reason for this result is that no primaries are
recorded at small angles (because of the initial offset gap), and no
multiples are recorded with large incident angles. The primaries
that generated the multiples that are present in the near traces 
were not recorded, and so, can not be used to predict their
multiples. Another way to see the problem is to recall that the 
upward continuation process corresponds to a propagation of
energy from small to large offsets and in the direction of
increasing time. Therefore, the near traces in the model
will not have enough energy to properly represent the multiples.
It is important to notice also that in real data, the source's  
directiveness will have a major influence on $ \alpha(p) $.

Although the efficiency of the method at small offsets is
not critical for the present goals of recovering the converted
waves where they are sufficiently strong (large offsets), a
possible way to improve the performance at near
offsets is the use of an extrapolation process to fill in
the gap before the upward continuation.

\end{document}
