\lefthead{Guerra}
\righthead{Residual multiple attenuation}
\email{guerrac@stanford.edu}
\title{Residual multiple attenuation using AVO modeling}
\author{Claudio Guerra}

\ABS{Residual multiple attenuation is a common process applied toward the end of the seismic processing flow.
Its main objective is to decrease the energy of migrated residual multiples, which can be time-consuming. 
I present a method which uses AVO modeling to simulate primary events and adaptive subtraction to estimate residual multiples and residual-multiple-attenuated data. The method is applied
 on angle-domain common-image gathers (ADCIGs) from a 2D dataset from the Gulf of Mexico that is contaminated with residual multiples. This method is simple, fast and, depending on the 
 parameterization of the adaptive subtraction, can preserve primary reflection amplitudes.}

\vspace{.5cm}

\section{Introduction} 

In marine data processing, it is common practice to run residual multiple attenuation as one of the later steps --- 
after migration, for instance. Assuming that the appropriate strategy for multiple attenuation is used, the multiple-attenuation tests do not access every line of 3D data, and 
the evaluation of the efficacy of the attenuation is performed along control lines. Since the multiple-attenuation
 parameters are optimized for the control lines, unfortunately, a surprising amount of residual-multiple energy from other lines can be spread by migration over the entire migrated dataset. Sometimes, the evaluation of the 
multiple attenuation is performed by inspecting post-stack migrated data, or even stacked volumes. Since this ignores the
 effect of pre-stack migration on the multiple-attenuated data, is not accessed during the multiple-attenuation tests.
 Eventually, after pre-stack migration, residual multiples can appear like crosshatching
 patterns and high-frequency migration noise resembling overmigrated events, because of the faster primary velocities
 used to migrate them.
\par
Another very common process after pre-stack migration is automatic residual-moveout correction to correct for 
velocity inaccuracies, which aims to align reflectors horizontally to improve the stacking power and better constrain
 the data for AVO analysis. Consequently, a method which relies on the flatness of the gathers is a natural process to
 attenuate residual multiples after automatic residual-moveout correction, as long as it is not time consuming. The problem of residual moveout will not be addressed in this work.
\par
\longcite{wegtle} classified multiple-attenuation methods into two main categories: (1) filtering methods and (2) prediction methods. The most widely used filtering methods rely on differences in traveltime between primaries and multiples to separate them in an appropriate domain where multiples  can be zeroed out \cite{SEG-1987-BP3.CSEG}. Since, in spite of the domain, primaries and multiples normally overlap, the filtering is not perfect and some of the primary energy can be attenuated or some of the multiple energy is kept (crosstalk). Prediction methods simulate multiples by auto-convolution of shot gathers and then subtract them from the original data in a nonlinear adaptive subtraction fashion \cite{GEO.57.9.11661177}. As prediction methods became popular, filtering methods inherited the nonlinear adaptive subtraction to attack the crosstalk problem. Nowadays, filtering methods are being used to simulate the multiples by exploring the same differences between primaries and multiples as before but, instead of simply filter multiples out, they are subtracted from the original data using the nonlinear adaptive subtraction scheme. This provides more flexibility to solve specific problems.
\par 
Generally, residual-multiple attenuation is performed by filtering methods, zeroing the residual-multiple energy (along with some primary) or using nonlinear adaptive subtraction. 
\longcite{choo} devised a method to perform residual-multiple attenuation in a different manner. They apply AVO inversion to simulate the primaries and then estimate the residual multiples using nonlinear adaptive subtraction. In spite of the unpublished details of their work, the method I present here is, in essence, very similar to theirs. The
 method simulates the primaries according to an AVO curve, estimates residual multiples by nonlinear adaptive 
 subtraction, and, finally, estimates the primaries. 
This approach solves the AVO modeling problem globally. 
\par
The method is applied on a Gulf of Mexico (GOM) 2D dataset comprised of angle-domain common-image gathers (ADCIGs) which have been subjected to multiple attenuation after migration. The dataset contains very complex multiple patterns \cite{Alvarez.sep.115.gabriel1}.

\section{Description of the method}       
\longcite{aki} stated that, up to $40^{\circ}$, the {\it p-wave} reflection coefficients, $R(\theta)$, as a
 function of the incidence angle, $\theta$, can be approximated by 
\begin{equation}
R(\theta) = A + B sin^{2}{\theta} + C tan^2{\theta}, 
\label{eq:aki}
\end{equation}
\noindent
{where A, B and C are parameters related to the average and contrasts of the elastic properties of the limiting layers. Therefore, if it is
possible to determine these parameters, one can simulate the amplitudes of primary events at every depth or time step.
 I refer to the simulation of the primary amplitudes as AVO modeling. The aim of this work is not
 to invert the AVO parameters for rock and fluid properties, but just to obtain a reasonable estimate of the primaries to be further used in a residual-multiple-attenuation sequence.}
\par
The AVO modeling problem can be addressed by solving locally, for every different depth (or time) step, the normal equations $\bold {(G^{^T}G)} \bold {\hat m} = \bold {G^{^T}} \bold d$,
where $\bold {\hat m}$ are the estimated AVO parameters, $\bold G$ has the column vector form
 $[\bold 1 $ $\bold{sin^2 \theta_i} $ $\bold{tan^2 \theta_i}]$ and $\bold d$ is the ADCIG (or a CMP gather after NMO). The simulated primaries are obtained by using the estimated AVO parameters in equation (\ref{eq:aki}). One
 problem that arises from this local solution is that the AVO parameters are is not constrained to be smooth in the
 depth direction. This can lead to anomalous amplitudes in the simulated primaries for a certain depth.
\par
In the present approach, the AVO modeling problem consists of determining the three parameters A, B and C of equation
 (\ref{eq:aki}) by solving the following data fitting problem:
\begin{equation}
\bold L \bold m = \bold d ,\\
\label{eq:fit01}
\end{equation}
where $\bold L$ in matrix form is
\begin{equation}
  \left[
        \begin{array}{cccccccccccc}
1&0&\cdot \cdot \cdot&0&s_{_{11}}&0&\cdot \cdot \cdot&0&t_{_{11}}&0&\cdot \cdot \cdot&0\\
0&1&0&\cdot \cdot \cdot&0&s_{_{12}}&0&\cdot \cdot \cdot&0&t_{_{12}}&0&\cdot \cdot \cdot\\
0&0&\cdot&0&\cdot \cdot \cdot&0&\cdot&0&\cdot \cdot \cdot&0&\cdot&0\\
0&\cdot \cdot \cdot&0&1&0&\cdot \cdot \cdot&0&s_{_{1n_2}}&0&\cdot \cdot \cdot&0&t_{_{1n_2}}\\
1&0&\cdot \cdot \cdot&0&s_{_{21}}&0&\cdot \cdot \cdot&0&t_{_{21}}&0&\cdot \cdot \cdot&0\\
0&1&0&\cdot \cdot \cdot&0&s_{_{22}}&0&\cdot \cdot \cdot&0&t_{_{22}}&0&\cdot \cdot \cdot\\
0&0&\cdot&0&\cdot \cdot \cdot&0&\cdot&0&\cdot \cdot \cdot&0&\cdot&0\\
0&\cdot \cdot \cdot&0&1&0&\cdot \cdot \cdot&0&s_{_{2n_2}}&0&\cdot \cdot \cdot&0&t_{_{2n_2}}\\

1&0&\cdot \cdot \cdot&0&s_{_{n_11}}&0&\cdot \cdot \cdot&0&t_{_{n_11}}&0&\cdot \cdot \cdot&0\\
0&1&0&\cdot \cdot \cdot&0&s_{_{n_12}}&0&\cdot \cdot \cdot&0&t_{_{n_12}}&0&\cdot \cdot \cdot\\
0&0&\cdot&0&\cdot \cdot \cdot&0&\cdot&0&\cdot \cdot \cdot&0&\cdot&0\\
0&\cdot \cdot \cdot&0&1&0&\cdot \cdot \cdot&0&s_{_{n_1n_2}}&0&\cdot \cdot \cdot&0&t_{_{n_1n_2}}\\
        \end{array}
  \right]
  \ \ 
\label{eq:Loper}
\end{equation}
which has dimension $n_{_1}n_{_2} X 3n_{_1}$, $\bold m$ is the model-parameter vector ($3n_{_1}$ elements), and $\bold d$ is the data vector ($n_{_1}n_{_2}$ elements), 
where $n_{_1}$ and $n_{_2}$ are the number of samples in depth (or time) and number of traces in the ADCIG, respectively. 
The terms $s_{_{ij}}$ and $t_{_{ij}}$ stand for sine and tangent squared terms at depth $i$ and for reflection angle $j$.
\par
Generally, in ADCIGs, residual multiples are more persistent at near angles, because of insufficient moveout
 difference between 
them and the flattened primaries. Additionally, at the farthest angles, stretch occurs. To avoid these imperfections in the input data, the fitting goal is to minimize the residual,
\begin{equation}
\bold 0 \approx \bold M (\bold L \bold m - \bold d),\\
\label{eq:fit02}
\end{equation}
where $\bold M$ is a selector operator (mute).
\par
If unrealistic variations of the AVO parameters with depth are an issue, another fitting goal is required. The model 
can be regularized with
\begin{equation}
\bold 0 \approx \epsilon \bold D \bold m ,\\
\label{eq:fit03}
\end{equation}
where $\bold D$ is the derivative operator. As usual, care must be taken when choosing $\epsilon$ not to destroy 
any recoverable residual moveout information and to not spread simulated primaries to angles in
 which they originally do not occur.
\par
Finally, to accelerate the solution, the problem can be solved with preconditioning \cite{gee}, using 
$\bold m = \bold C \bold p$, where $\bold C$, the anti-causal integration, is the inverse of the causal
 differentiation, $\bold D$. After the determination of the three AVO parameters, the simulated primaries are obtained by substituting them into
 equation (\ref{eq:aki}). Of course, the method relies on the flatness of the reflectors to correctly extract the 3 
parameters of the AVO curve.
\par
To get the residual-multiple-attenuated data, the adaptive subtraction must be applied in two steps. The first one 
aims to obtain an estimate of the residual multiples. This is done by subtracting the adjusted version of the 
simulated primaries from the original data. Amplitude and phase adjustments are achieved by the convolution of the
 simulated primaries with a prediction-error Wiener filter.
\par
As the estimated residual multiples may contain some primary information (mainly at small reflection angles), I use the
 strategy proposed by \longcite{Guitton.sep.108.antoine2}, the so-called `subtraction method,' in which two different prediction-error filters (PEFs), 
which model primaries and multiples, are computed. Their method allows regularization,
 to decrease the crosstalk between multiples and primaries. The corresponding fitting goals are
\begin{center}
$\bold 0 \approx \bold m_n$\\
$\bold 0 \approx \epsilon \bold m_s$\\
subjected to  $\leftrightarrow \bold d = \bold S^{-1} \bold m_s + \bold N^{-1}\bold m_n$\\
\end{center}
after preconditioning, with $\bold s = \bold {S^{-1}m_s}$ and $\bold n = \bold {N^{-1}m_n}$. In the previous equations, 
$\bold s$ represents the primaries, $\bold n$ the multiples; $\bold S$ and $\bold N$ characterize convolution with 
PEFs that whiten primaries and multiples, respectively; $\bold {m_s}$ is the model of primaries and 
$\bold {m_n}$ the model of multiples. Finally, the estimated primaries are obtained by 
$\bold {\hat s} = \bold d - \bold N^{-1} \bold {m_n}$.\\

\section{Examples}
The dataset used to test the method consist of ADCIGs of a 2D line from the GOM. They have been subjected to multiple attenuation
 using the apex-shifted tangent-squared radon transform \cite{Alvarez.sep.115.gabriel1}. Figure (\ref{fig:migstk}), 
which corresponds to the stack of ADCIGs before the previously mentioned multiple-attenuation method, shows complex
 patterns of multiples, especially below the edges of the salt, corresponding to migrated diffracted multiples.
 Around a depth of 3500 m, the first-order multiple of the sea bottom is probably distorted due
 to the alternation between fast and slow velocities. The fitting goal I used to simulate the primaries corresponds to equation (\ref{eq:fit02}).
\par
\plot{migstk}{width=6.00in}{Migrated data with no multiple attenuation applied. Notice strong migrated multiples .}

The results are shown in CIGs near CMP locations 2000, 6000, 12000, and 20000 m, in zero reflection-
angle gathers, and in stacked sections. Figures (\ref{fig:cig02000}), 
(\ref{fig:cig06000}), (\ref{fig:cig12000}) and (\ref{fig:cig20000}) illustrate the ADCIGs. In all of the figures, the two top
 panels show the original ADCIG on the left and the residual-multiple-attenuated one on the right; the bottom left
 panels show the estimated residual multiples, and the bottom right show the simulated primaries. 
The residual-multiple panels tell us that not only multiples, but every non-flat event in the
 ADCIG is considered noise to be attenuated (e.g. kinematic effects due to
 illumination problems and artifacts due to the transformation from subsurface offset to reflection angle). Primaries
 with significant residual moveout may also be considered events to be attenuated, so it is important to
 consider the application of residual-moveout correction before the residual-multiple attenuation to preserve primary information.
 In this particular case, the velocity model seems to be less accurate for deeper reflectors.
\par
\plot{cig02000}{width=4.00in,height=8.00in}{ADCIG --- CMP location 2000 m. Top left: multiple-attenuated data; top right: residual-multiple-attenuated data; bottom left: estimated multiples; and bottom right: simulated primaries.}

\plot{cig06000}{width=4.00in,height=8.00in}{ADCIG --- CMP location 6000 m. Top left: multiple-attenuated data; top right: residual-multiple-attenuated data; bottom left: estimated multiples; and bottom right: simulated primaries.}

\plot{cig12000}{width=4.00in,height=8.00in}{ADCIG --- CMP location 12000 m. Top left: multiple-attenuated data; top right: residual-multiple-attenuated data; bottom left: estimated multiples; and bottom right: simulated primaries.}

\plot{cig20000}{width=4.00in,height=8.00in}{ADCIG --- CMP location 20000 m. Top left: multiple-attenuated data; top right: residual-multiple-attenuated data; bottom left: estimated multiples; and bottom right: simulated primaries.}

Another issue in methods that rely on the moveout difference between primaries and multiples is the persistence of 
multiples in the near angles (or offsets) due to insufficient moveout difference between them. In sub-salt regions, this problem is especially
 prominent in ADCIGs because of the poor illumination of farther angles. To guarantee that residual multiples
 in the near-angle range do not interfere with the computation of the simulated primaries, I applied an internal mute, 
roughly below depth=3000 m between angles $-10^{\circ}$ and $10^{\circ}$. At the CMP location 12000 m; what seem to be primaries
 in the near angle are actually residual multiples, as can be clearly seen in the zero-angle section (Figure 
(\ref{fig:or00deg})). The zero-angle section of the residual-multiple-attenuated data in Figure (\ref{fig:pest00deg}) shows that the low-frequency contents of the residual multiples remain after the process.
\par
\plot{or00deg}{width=6.00in}{Multiple-attenuated data --- zero-reflection-angle gather.}

\plot{pest00deg}{width=6.00in}{Residual-multiple-attenuated data --- zero-reflection-angle gather.}

Disregarding the imperfections which break the assumption of flatness, by comparing the original and the attenuated
 ADCIGs we can see that the results are promising. Most of the residual-multiple energy has been attenuated.
\par
Interestingly, below the edges of the salt body where diffracted multiples are more prominent, parts of the residuals 
of diffracted multiples, which look like migration smiles dipping to the right,  remain untouched 
(e.g. CMP location 20000 m). Conversely, diffracted migrated multiples dipping to the left are well attenuated. 
By inspecting the corresponding ADCIG, (Figure \ref{fig:cig20000}) one can
 see that the apex of the untouched diffracted multiples is out of the internal mute (roughly $10^{\circ}$) applied to simulate the primaries.
 Additionally, their tops are relatively flat. Therefore, they give rise to ``false'' simulated primaries and are not
 attenuated at all.
\par
In spite of the simplicity of the assumptions, the method provides a cleaner image for interpretation, as the 
 comparison of Figures (\ref{fig:orstk}) and (\ref{fig:peststk}) shows. Detailed versions of these figures can be 
found in Figure (\ref{fig:zstk}). Additionally, Figure (\ref{fig:zzo}) shows details of the 
 the zero-reflection-angle gathers, at the same locations as the stacked data. As can be seen by inspecting the upper frames in Figure 
(\ref{fig:zzo}), after the residual-multiple attenuation, some reflectors have been revealed. 
This final result obtained for the near-angle section has great impact in AVO analysis. Unfortunately, it is very difficult to observe a typical reflector
 below the edge of the salt body, in the zero-reflection-angle gathers, making it hard to evaluate the efficacy of the method 
in preserving primaries in such a complex geology. However, as can be observed in the bottom frames of Figure (\ref{fig:zzo}), 
the residual multiples are not as strong as before.

\plot{orstk}{width=6.00in}{Multiple-attenuated data. Residual multiples dominate the deeper portion of the data.}

\plot{peststk}{width=6.00in}{Data after residual-multiple attenuation. The final data is much cleaner than the 
multiple-attenuated data.}

\plot{zstk}{width=6.00in}{Details of Figures (\ref{fig:orstk}) and  (\ref{fig:peststk}). Notice persistent multiples 
in the original data on the left and how attenuated they are on the right panels.}

\plot{zzo}{width=6.00in}{Details of Figures (\ref{fig:or00deg}) and  (\ref{fig:pest00deg}) --- zero-reflection-angle data 
at the same portion as shown in the previous figure. Again, on the left panels is the input data and on the right, 
the residual-multiple-attenuated data.}

\section{Conclusions}
I present a simple and fast method, based on AVO modeling, to obtain simulated primaries for estimating
 residual multiples, in such a way that they can be further attenuated. The method relies on the flatness of the input 
CIGs. I apply the method to a GOM 2D dataset which has strong residual multiples with complex geometries. 
The method proves its efficiency in attenuating residual multiples, yielding cleaner images than the original ones. In principle, any residual noise which does not obey the flatness criteria in a CIG can be attenuated by the method. 

\section{Acknowledgments}
I thank Gabriel Alvarez for introducing me to the multiple-attenuated migrated data and for help and discussions, and
 Antoine Guitton for the program to perform adaptive subtraction. Thanks to Bob Clapp for helping me to run Guitton's 
program in parallel. I also thank the sponsors of the Stanford Exploration Project for their support and WesternGeco 
for providing the dataset.

\bibliographystyle{sep}
\bibliography{SEP,ini1}
