\lefthead{Guerra}
\righthead{Residual multiple attenuation}
\email{guerrac@stanford.edu}
\footer{SEP--129}
\title{Residual multiple attenuation using AVA modeling}
\author{Claudio Guerra}

\ABS{Residual multiple attenuation is a common process applied towards the end of the seismic processing flow.
Its main objective is to decrease the energy of migrated residual multiples. This process can be time-consuming. 
I present a method which uses amplitude-versus-angle (AVA) 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 data set 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} 
Conventional time and depth imaging considers that the input data is made up of only primary reflections. Because of this basic assumption, much effort is spent attenuating events considered as noise. In the marine case, multiple reflections are the most important noise. If not successfully attenuated, multiple reflections can bias velocity estimation, and after migration, give rise to images contaminated with residual multiples which can make interpretation difficult. For example, such contamination may severely interfere with AVA analysis.
\par
In marine data processing, it is a common practice to run residual multiple attenuation as one of the later steps --- 
after migration, for instance. 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, sometimes, a surprising amount of residual-multiple energy is left and is spread by migration over the entire migrated data set. 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, it can allow residual multiples to appear, manifesting as crosshatched patterns and high-frequency migration noise resembling overmigrated events.
\par
Another very common process after pre-stack migration is automatic residual-moveout correction to mitigate the effects of velocity inaccuracies, which aims to align reflectors horizontally to improve the stacking power and better constrain the data for AVA 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. However, 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 appropriate domains where multiples  can be filtered out \cite{SEG-1987-BP3.CSEG}. Since, regardless of the domain, primaries and multiples normally overlap, particularly at small reflection angles and offsets, the filtering is not perfect; therefore, some of the primary energy can be attenuated, or some of the multiple energy is left out (crosstalk). Prediction methods simulate multiples by auto-convolution of shot gathers and then subtract them from the original data using nonlinear adaptive subtraction \cite{GEO.57.9.11661177}. As prediction methods became popular, filtering methods inherited the nonlinear adaptive subtraction to attack the crosstalk problem. So, filtering methods can be used to estimate multiples by exploring the same differences between them and primaries as before, but instead of simply filtering them out, a possible approach is to subtract the estimated multiples from the original data using nonlinear adaptive-subtraction schemes. This provides more flexibility to solve specific problems.
\par 
Generally, residual-multiple attenuation is performed by filtering methods, by zeroing the residual-multiple energy (along with some primary energy) or by using nonlinear adaptive subtraction. 
\longcite{choo} devised a method to perform residual-multiple attenuation in a different manner. They apply AVA inversion to simulate the primaries and then estimate the residual multiples using nonlinear adaptive subtraction. Although many details of their work are unpublished, the method I present here is, in essence, very similar to the one proposed by them. The method simulates the primaries according to an AVA curve, adjusts amplitude and phase of the simulated primaries to the amplitude and phase of the original data by applying filters computed in a nonlinear adaptive strategy, subtracts the adjusted primaries from the original data to get an estimate of the residual multiples, and, finally, obtains the attenuated data by nonlinear adaptive subtraction. This approach solves the AVA modeling problem globally. 
\par
The method is applied on angle-domain common-image gathers (ADCIGs) of a Gulf of Mexico (GOM) 2D data subjected to multiple attenuation after migration. The data set 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 AVA modeling. The aim of this work is not
 to invert the AVA 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 AVA 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 AVA 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 AVA parameters in equation (\ref{eq:aki}). One
 problem that arises from this local solution is that the AVA parameters are not constrained to be smooth along depth. This can lead to anomalous amplitudes in the simulated primaries for a certain depth.
\par
In the present approach, the AVA 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}{ccc}
{\bold 1}_1&{\bold S}_1&{\bold T}_1\\
{\bold 1}_2&{\bold S}_2&{\bold T}_2\\
{\bold 1}_3&{\bold S}_3&{\bold T}_3\\
\vdots&\vdots&\vdots\\
{\bold 1}_{n_1}&{\bold S}_{n_1}&{\bold T}_{n_1}\\
        \end{array}
  \right]
  \ \ 
\label{eq:Loper}
\end{equation}
which has dimension $n_{_1}n_{_2} \times 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 the number of traces in the ADCIG, respectively. ${\bold 1}_j$, ${\bold S}_j$ and ${\bold T}_j$ are $n_{_2} \times n_{_1}$ matrices consisting of $1$, $(\sin^2 \theta_i)$ and $(\tan^2 \theta_i)$ along their $j$-th column entries, respectively. $i$ stands for the reflection-angle indexes.
\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 which applies the appropriate internal and external mutes.
\par
If unrealistic variations of the AVA parameters with depth are an issue, the following regularization goal can be introduced:
\begin{equation}
\bold 0 \approx \epsilon \bold D \bold m ,\\
\label{eq:fit03}
\end{equation}
where $\bold D$ is the derivative operator along depth and $\epsilon$ is the regularization parameter. As usual, care must be taken when choosing $\epsilon$ not to destroy 
any recoverable residual moveout information and not to spread simulated primaries to angles at
 which they originally do not occur.
\par
To accelerate the solution, the problem can be solved with preconditioning \cite{gee}, using the transformation $\bold m = \bold C \bold p$, where $\bold C = D^{-1}$ and $\bold p$ is the preconditioned variable. Finally, the fitting goals reduce to
\begin{equation}
\begin{array}{c}
\bold 0 \approx \bold M (\bold {LC} \bold p - \bold d)\\
\bold 0 \approx \epsilon \bold p .\\
\end{array}
\label{eq:fit04}
\end{equation}

The final model is obtained with $\bold m = \bold C \bold p$. I use conjugate-gradients to solve the inverse problem.
\par
After the determination of the three AVA parameters for every depth step, primaries are simulated by computing the reflection coefficient for all reflection angles using equation (\ref{eq:aki}).
 Of course, the method relies on the flatness of the reflectors in the CIG to correctly extract the 3 
parameters of the AVA 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 preconditioning equations, $\bold s$ represents the primaries, $\bold n$ the multiples; $\bold N^{-1}$ and $\bold S^{-1}$ represent deconvolution with PEFs for multiples and primaries, respectively; $\bold {m_n}$ is the multiple model content and $\bold {m_s}$ is the primary model content. Finally, the estimated primaries are obtained by $\bold {\hat s} = \bold d - \bold N^{-1} \bold {m_n}$. All PEFs are computed with conjugate-gradients.


\section{Examples}
The method is applied on ADCIGs of a Gulf of Mexico (GOM) 2D data. The data have been migrated with shot-profile algorithm and the transformation from subsurface-offset to angle was performed after imaging. The ADCIGs 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 without multiple-attenuation, shows complex
 patterns of multiples, especially below the edges of the salt, corresponding to 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 is the one given in equation (\ref{eq:fit04}).

\plot{migstk}{width=6.00in,height=3.00in}{Migrated data without multiple attenuation -- stacked data. Notice strong migrated multiples, especially below the edges of the salt, corresponding to migrated diffracted multiples.}
 
The results are shown in CIGs near CMP locations 2000, 12000, and 20000 m, in zero-angle sections, and in stacked sections. Figures \ref{fig:cig02000}, \ref{fig:cig12000} and \ref{fig:cig20000} illustrate the ADCIGs. In all of the figures are shown in panel (a) the original ADCIG, panel (b) the residual-multiple-attenuated, panel (c) the estimated residual multiples, and panel (d) 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 the primary information. In this particular case, the velocity model seems to be less accurate for deeper reflectors. The estimated primary panels show significant residual-multiple attenuation, specially at smaller reflection angles.

\plot{cig02000}{width=5.00in,height=8.00in}{ADCIG --- CMP location 2000 m. (a): multiple-attenuated data; (b): residual-multiple-attenuated data; (c): estimated multiples; and (d): simulated primaries. Notice significant residual-multiple attenuation, specially at smaller reflection angles.}

\plot{cig12000}{width=5.00in,height=8.00in}{ADCIG --- CMP location 12000 m. (a): multiple-attenuated data; (b): residual-multiple-attenuated data; (c): estimated multiples; and (d): simulated primaries. Notice significant residual-multiple attenuation, specially at smaller reflection angles.}

\plot{cig20000}{width=5.00in,height=8.00in}{ADCIG --- CMP location 20000 m. (a): multiple-attenuated data; (b): residual-multiple-attenuated data; (c): estimated multiples; and (d): simulated primaries. Notice, in the estimated residual-multiples panel, the apex of diffracted multiples (around depth = 3700 m) located at reflection angles = $\pm 20^{\circ}$.}
 
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 and the primaries. In sub-salt regions, this problem is especially prominent in ADCIGs because of the poor illumination of larger 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 may 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 a low-frequency content of the residual multiples remains after the process. Figure \ref{fig:zzo} shows details of the zero-angle sections. As can be seen by inspecting the Figures \ref{fig:zzo}a and \ref{fig:zzo}b, after the residual-multiple attenuation, some reflectors have been revealed. This final result obtained for the near-angle section has great impact in AVA 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 Figures \ref{fig:zzo}c and \ref{fig:zzo}d, the residual multiples are not as strong as before.
\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 (around depth = 3700 m and reflection angles = $\pm 20^{\circ}$) is out of the internal mute (roughly $10^{\circ}$) applied. Additionally, their tops are relatively flat. Therefore, they give rise to ``false'' simulated primaries and are not attenuated at all.
\par
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
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. A zoomed version is presented in Figure \ref{fig:zstk}. 

\section{Conclusions}
I present a simple and fast method, based on AVA 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 data set 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, Guojian Shan and Madhav Vyas for discussions and suggestions, and Antoine Guitton for the program to perform adaptive subtraction. I also thank WesternGeco for providing the data set.
\vspace{.5cm}
\plot{or00deg}{width=6.00in,height=3.00in}{Multiple-attenuated data --- zero-angle section. Strong migrated multiples make interpretation of primaries impossible.}

\plot{pest00deg}{width=6.00in,height=3.00in}{Residual-multiple-attenuated data --- zero-angle section. Residual multiples have been attenuated.}

\plot{zzo}{width=6.00in,height=6.00in}{Details of Figures \ref{fig:or00deg} and  \ref{fig:pest00deg} --- zero-angle sections. (a) and (c): multiple-attenuated data. (b) and (d): residual-multiple-attenuated data.}
  
\plot{orstk}{width=6.00in,height=3.00in}{Multiple-attenuated data --- stacked data. Residual multiples dominate the deeper portion of the data.}

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

\plot{zstk}{width=6.00in,height=6.00in}{Details of Figures \ref{fig:orstk} and  \ref{fig:peststk} at the same portion as shown in Figure \ref{fig:zzo}  --- stacked data. (a) and (c): multiple-attenuated data. (b) and (d): residual-multiple-attenuated data. Notice persistent multiples in the original data on panels (a) and (c) and how attenuated they are on panels (b) and (d).}

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