\lefthead{Guerra and Valenciano}
\righthead{Multiple attenuation prior to inversion}
%\email{\footnotemark[1]guerrac@stanford.edu \footnotemark[2]valencia@stanford.edu}
\footer{SEP-134}
\title{Reconciling processing and inversion: Multiple attenuation prior to wave-equation inversion}
\author{Claudio Guerra and Alejandro Valenciano}

\maketitle

\begin{abstract}
Seismic inversion is very sensitive to the presence of noise. In an inversion scheme, noise is any event in the data not predicted by the forward modeling because of either inaccurate physics or inadequate model parameterization. Therefore, noise exists in the residual space and, if coherent, slows convergence to an acceptable result, or in the worst case, dominates the whole process, inhibiting the efficacy of inversion. The ideal solution is to incorporate the modeling of the noise into the forward modeling. 
A more practical  approach is to make the data agree with the physical assumptions of the inversion scheme. In this context, noise attenuation is a pre-processing step before inversion. Here, we illustrate the problem by applying one-way wave-equation inversion to a portion of the well known Sigsbee2b data. In the present example, noise takes the form of multiples, not modeled by the one-way wave-equation. After characterizing the noise in the migrated data, we use a dip filter to estimate it and a non-stationary adaptive filter technique to subtract it from the migrated data. 
\end{abstract}

\section{Introduction} 
Migration with the full wave equation is, in principle, capable of creating an image of the subsurface by reversing the propagation of all wave modes \cite[]{wapenaar10}. It is a computationally demanding process, which is why the wave equation is approximated by the one-way wave equation. Migration with the one-way approximation is formulated such that, ideally, the only seismic events present in the input data are primary reflections. When this is not the case, because the pre-processing was unable to completely eliminate noise (multiple reflections, for instance), the final image presents migrated events that do not correspond to the geology. In other words, these events are not modeled by the physical theory on which the migration method is based. They are frequently called migration noise or migration artifacts. 
\par
As is widely known, inversion problems are very sensitive to noise. In seismic, several authors have shown that coherent noise degrades the performance of the velocity inversion and yields erroneous results \cite[]{chauris, shin01, li01}. Therefore it is usually necessary to include some kind of a priori information about the model as an additional regularization equation in the inversion or to perform noise suppression prior to inversion.
\par
Noise also affects least-squares inversion of migrated images into reflectivity. In offshore data, multiples are the main coherent noise. \cite{Valenciano135} addresses the problem of migrated multiple reflections in inversion by applying a pre-processing step to attenuate them. Here, using the Sigsbee2b dataset \cite[]{paffenholz}, we characterize the multiples in the subsurface-offset domain, detail the pre-processing step and show the impact it has on inversion.

\section{Linear least-squares inversion}
Linear least-squares inversion provides a theoretical approach to compensate for experimental deficiencies (e.g., limited acquisition geometry) and complexities of the overburden, while maintaining consistency with the acquired data. For seismic imaging, it can be summarized as follows. Given a linear modeling operator ${\bf L}$, compute synthetic data {\bf d} using $ {\bf d}={\bf L}{\bf m}$, where {\bf m} is a reflectivity model.

The quadratic cost function,
\begin{equation}
 S({\bf m})=\Vert {\bf d} - {\bf d}_{obs}  \Vert^2
        =\Vert {\bf L}{\bf m} - {\bf d}_{obs}  \Vert^2, 
\end{equation}
 is formed, where ${\bf d}_{obs}$ denotes the recorded data.

The reflectivity model $\hat{{\bf m}}$ that minimizes $S({\bf m})$ is given by
\begin{equation}
\hat{{\bf m}}=({\bf L}'{\bf L})^{-1}{\bf L}' {\bf d}_{obs}
= {\bf H}^{-1} {\bf m}_{mig},
\label{eqn:tara}
\end{equation}
where ${\bf L}'$ (the migration operator) is the adjoint of the linear modeling operator ${\bf L}$, ${\bf m}_{mig}$ is the migrated image, ${\bf m}_{mig}={\bf L}' {\bf d}_{obs}$, and ${\bf H}={\bf L}'{\bf L}$ is the Hessian of $S({\bf m})$.

%Two possible ways to obtain the inverse image are suggested by equation \ref{eqn:tara}. One solves the linear system by iteratively evaluating the forward operator ${\bf L}$ and its adjoint  ${\bf L}'$ (implicit method), and the other solves it by explicitly computing the Hessian matrix.  
%The choice between the two methods generally depends on the problem. In seismic imaging, the forward operator and its adjoint can be represented by full matrices, making their evaluation a computationally expensive process. Each iteration of the implicit method using a solver such as conjugate gradients involves one evaluation of the forward and one of the adjoint operator. Because, according to theory, the conjugate-gradient method converges in a number of iterations roughly equal to the size of the model space, iteration until convergence with the implicit method is impossible for today's computer resources \cite[]{Mclapp122}. Still, \cite[]{Mclapp122} proved that after a small number of iterations, the image resulting from the inversion is better that the image resulting from the migration.

%The explicit method has a big advantage over the implicit method: after the Hessian matrix is computed, different inversion schemes can be used without having to incur the cost of wavefield propagation. Another advantage is in the parameterization of the inversion. In practice, inversion always requires parameter adjustments, which can be very expensive in the implicit case, since every time an inversion parameter is changed, wavefied propagation will be needed.
%But the explicit method also has problems.

Since the model space can be large, computing the inverse of the Hessian matrix is a big challenge for most geophysical imaging problems. For this reason, it is often more feasible to compute the inverse image as the solution of the linear system of equations,
\begin{equation}
{\bf H} \hat{{\bf m}}={\bf m}_{mig},
\label{eqn:lse}
\end{equation}
by using an iterative inversion algorithm. In this approach to the inverse problem, only one matrix-vector multiplication of the Hessian matrix with the model vector is necessary per iteration, in contrast with other methods \cite[]{Mclapp122} that require one migration and one modeling every iteration. Still, explicit computation of the Hessian for the entire model space is too expensive in practice. \cite{Valenciano135} discusses how exploiting the structure of the Hessian matrix and the localization in the model space makes this problem tractable.

By using a priori information about the covariance of the model (model regularization), one can add regularization to solve the otherwise ill-posed inversion problem. A more customary regularization for the inversion in the poststack image domain is to add a damping factor that penalizes an increase of the values of the model. This regularization makes no use of any physical knowledge we might have about the seismic reflectors. It is implemented by adding a small value to the diagonal of ${\bf H}({\bf x},{\bf x}')$ in equation \ref{eqn:lse}:
\begin{equation}
 \left({\bf H+\varepsilon I}\right)\hat{{\bf m}} - {\bf m}_{mig} = {\bf r} \approx 0 ,
\label{eqn:damping}
\end{equation}
where ${\bf I}$ is the identity operator, ${\bf r}$ is the residuals vector, and $\varepsilon$ is an scalar parameter that governs the strength of the regularization.

\section{Noise Characterization}
The occurrence of shadow zones below the salt body is strongly conditioned by the concavity of its base. In the Sigsbee2b model, the comparison of the diagonal of the Hessian matrix in Figure \ref{fig:diag-Sis} with the reflectivity model in Figure \ref{fig:refle-new} shows that where the base of the salt is concave down, the transmitted energy is focused and high illumination values are produced; where the base of the salt is concave up, the transmitted energy is spread along to divergent paths, producing low illumination zones. Unfortunately, for reflected energy and also for multiple energy, this geometry acts in the opposite way, concentrating energy where the base of the salt is concave up and dispersing energy where it is concave down. In the shadow zones it is very likely that the strongest remaining amplitudes correspond to multiple reflections, while in the well-illuminated zones these events produce little or no spurious interference with the primary energy (Figure \ref{fig:mig-stack}). Consequently, during one-way wave-equation inversion, multiple energy can dominate the residuals, surpassing primaries in the shadow zones.
\sideplot{diag-Sis}{width=3.5in}{Illumination pattern given by the diagonal of the Hessian matrix of the Sigsbee2B model. Dark gray represents low illumination and light gray represents high illumination. [{\bf ER}]}
\plot{refle-new}{width=4in}{Reflectivity model of Sigsbee2B model.  [{\bf ER}]}
\plot{mig-stack}{width=4in}{Migration of Sigsbee2b data -- zero-subsurface-offset section. [{\bf ER}]}

One solution to this problem is to attenuate the multiples (events that, in inversion, do not fit the model) from the migration. This requires that the multiples be separated from the primaries or modeled in such a way that during subtraction the primaries are minimally affected. 
\par
Traditional data-space demultiple schemes dealing either with periodicity or differential moveout become ineffective in sub-salt settings for two main reasons: a) complex ray paths make periodicity assumptions fail, and b) salt focusing effects concentrates sub-salt reflections in the near offset range, making Radon-type demultiple schemes ineffective. \cite{sava:V10} and \cite{alvarez:V97} demonstrated that primaries and multiples have different behavior in image space (subsurface-offset or reflection-angle), and that multiples can be adaptively subtracted from  the migration in the image domain without significantly affecting the primaries.
\par
The Sigsbee2B model \cite[]{paffenholz} has a high-reflectivity water-bottom. In contrast with Sigsbee2a, this characteristic generates strong peg-leg multiples, which bounce between the salt edge and the water-bottom. In Figure \ref{fig:ZOtime} we use part of the zero-offset section of Sigsbee2b data (non-free-surface data) to highlight these multiples and the salt limits in the time domain. The diffractions labeled 1 correspond to the undulating top of the salt; that labeled 2 is the reflection from its base; $1^{st}$-order peg-leg diffracted multiples, concurrent with the base of the salt, are  labeled 3; the $1^{st}$-order peg-leg multiple from the base of the salt is identified with label 4; and internal multiples generated inside the salt body are labeled 5.

\sideplot{ZOtime}{width=3.5in}{Sigsbee2b zero-offset time section. The labels correspond to 1) top of the salt; 2) base of the salt; 3) $1^{st}$-order peg-leg diffracted multiples; 4) internal multiples originated at the salt limits; and, 5) $1^{st}$-order peg-leg multiple from the base of the salt. [{\bf ER}]}
\par
After prestack depth-migration, multiple reflections are distorted, because they are downward propagated beyond their actual reflection depth, using the incorrect velocities of the time-concurrent primaries. If their reflection traveltimes are coincident with sub-salt reflections, their propagation directions are severely deviated; therefore, the shape they assume shows a strong imprint of the velocity model. 
\par
In Figure \ref{fig:sig2b01} the migration result of Sigsbee2b is shown in a magnified view. The frame on the left is the zero-subsurface-offset section, and the one on the right is a subsurface offset gather -- ODCIG -- taken at CMP position 38825 ft. The $1^{st}$-order peg-leg diffracted multiples, labeled 1, show the effect of over-migration, because they are migrated with salt velocity; labels 2 and 4 stand for migrated peg-leg multiples from the base of the salt, and immediately above them, labeled 3 and 5, are the migrated internal multiples. Some reflectors show focusing close to zero-offset. Notice the strong correlation between the shape of migrated multiples and the base of the salt caused by the deviation of the propagation direction mentioned above. For comparison, Figure \ref{fig:sig2a01} shows the migration of Sigsbee2a. In this case, multiples related to the water-bottom (labels 1, 2 and 4) are much weaker or absent, because the water-bottom in this model is characterized by a ``soft'' interface.

\plot{sig2b01}{width=6.0in,height=3.0in}{Sigsbee2b shot-profile migration using the cross-correlation imaging condition. The labels correspond to migrated events from 1) $1^{st}$-order peg-leg diffracted multiples; 2) $1^{st}$-order peg-leg multiple from the base of the salt; 3) internal multiples originated at the salt limits; 4) $1^{st}$-order peg-leg multiple from the base of the salt; and, 5) internal multiples originated at the salt limits. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at  $x=38600$ ft. [{\bf ER}]}

\plot{sig2a01}{width=6.0in,height=2.75in}{Sigsbee2a shot-profile migration using the cross-correlation imaging condition. The labels correspond to migrated events from 1) $1^{st}$-order peg-leg diffracted multiples; 2) $1^{st}$-order peg-leg multiple from the base of the salt; 3) internal multiples originated at the salt limits; 4) $1^{st}$-order peg-leg multiple from the base of the salt; and, 5) internal multiples originated at the salt limits. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at  $x=38600$ ft. Compare with Figure \ref{fig:sig2b01}. [{\bf ER}]}

Multiples are more evident in regions of low illumination. In Sigsbee data, these regions are associated with a concave-up base of the salt, in which downgoing energy is defocused and upgoing energy is focused. Much of the remaining energy in these regions corresponds to multiples. Figure \ref{fig:migoff1} exemplifies this problem. The ODCIG on the right frame was selected from a low illumination area. The horizontal to upward-curved events are migrated multiples, which contaminate the entire sub-salt section. Because of the low illumination, some reflectors appear as dipping events below 14000 ft. Notice that, in the zero-subsurface-offset section, because of the imprint of the salt velocity, all the multiple modes show different dips than the reflectors. 
\par
The different dip behavior of primaries and multiples can be used as a criterion to separate them using a dip filter in the $k_x-k_h$ wavenumber domain in a pre-processing step. After discrimination, and prior to subtraction from the original data, amplitude and phase of the estimated multiples must be adjusted. We perform this correction by using non-stationary filters, according to the strategy of \cite{alvarez125}. They advocate that the best adjustment between the estimated multiples and the multiples in the data is achieved by simultaneously matching the estimation of primaries and multiples to the data containing both. The strategy uses small overlapping patches of the input data to compute local filters in a least-squares inverse problem.
%
\par
%
The next section presents the results of pre-processing and its impact on the inversion output.

\plot{migoff1}{width=6.0in,height=3.0in}{Sigsbee2b shot-profile
  migration (subsurface-offset) using the cross-correlation imaging
  condition. The front face corresponds to the zero
  subsurface-offset image, and the side face corresponds to the
  subsurface-offset gather at $x=36800$ ft. Poorly illuminated areas
  are dominated by migrated multiples. [{\bf ER}]}
%by the presence of migrated multiples. [{\bf ER}]}


\section{Results}
The Hessian matrix for the Sigsbee2B model (diagonal shown in Figure \ref{fig:diag-Sis}) and the migration (Figure \ref{fig:mig-stack})  were computed using one-way wave-equation propagators. As result, the multiples were neither included in the modeling nor propagated with the correct velocity/direction during migration. For this reason, the multiples will be treated as noise for the inversion. In the following, we will show how they prevent the inversion fail from recovering the primary signal in the shadow zones. 
\par
Figure \ref{fig:migoff1} shows the migration result in the subsurface-offset domain. The front face corresponds to the image at zero subsurface offset, and the side face corresponds to the subsurface-offset gather at a position inside the shadow zone ($x=36800$ ft). Since the primaries are illuminated at few reflection angles and from a dominant slanted wave-propagation direction, their signature in the subsurface-offset domain is a slanted line. Because of the multiple bounces on dipping interfaces (15$^{\circ}$ watter-bottom and 30$^{\circ}$ salt flank), multiple energy is spread horizontally in the ODCIG, resembling a zero-reflection-angle event.
Figure \ref{fig:migang1} shows the migration result transformed to the reflection-angle domain. The front face corresponds to the image at the $0^\circ$ reflection angle, and the side face corresponds to the reflection-angle gather at a position inside the shadow zone ($x=36800$ ft). Notice how multiples  are concentrated around the $0^\circ$ reflection angle, and primaries are mapped at higher reflection angles. This behavior of the multiples and primaries can also be found in field datasets \cite[]{Valenciano135}.

\plot{migang1}{width=6.0in,height=2.75in}{Sigsbee2b shot-profile migration (reflection angle) using the cross-correlation imaging condition. The front face corresponds to the image at the $0^\circ$ reflection angle, and the side face corresponds to the reflection-angle gather at $x=36800$ ft. [{\bf ER}]}

We used a three-stage pre-processing strategy to attenuate the multiple energy before inversion. It includes discriminating the multiple energy in the $k_x-k_h$ plane (to generate a model of the multiples); performing amplitude and phase correction of the discriminated multiples in a least-square sense; and subtracting them from the original data. 

The discrimination of the multiple energy relies on the differences in dip patterns as described previously: multiples show up in the $k_x-k_h$ plane at low $k_h$ and high $k_x$, and primaries at high $k_h$ and low $k_x$. Therefore, we build a model of the primaries by submitting the migrated data to a $k_x-k_h$ filter at every depth step. Figure \ref{fig:not-filtered} shows a depth slice at $z=14150$ ft of the migration in the $k_x-k_h$ plane, and Figure \ref{fig:filtered} shows the result of filtering. The model for the multiples (Figure \ref{fig:diff-filtered}) is obtained after subtraction of the filtered result (Figure \ref{fig:filtered}) from the migration (Figure \ref{fig:not-filtered}). 

\sideplot{not-filtered}{width=3.5in}{Depth slice  at $z=14150$ ft of the migration in the $k_x-k_h$ plane. [{\bf ER}]}
%
\sideplot{filtered}{width=3.5in}{Primaries after bandpassing the migration (Figure \ref{fig:not-filtered})  in the $k_x-k_h$ plane. [{\bf ER}]}
%
\sideplot{diff-filtered}{width=3.5in}{The model for the multiples in the $k_x-k_h$ plane, obtained after subtraction of the filtered result (Figure \ref{fig:filtered}) from the migration (Figure \ref{fig:not-filtered}). [{\bf ER}]}

Since the separation of primaries and multiples is not perfect, some of the multiples are present in the model of primaries, and some primary energy leaks to the model of multiples. Because of this cross-talk between the estimates of primaries and multiples, it is desirable to control the amount of attenuation. Therefore, the impact on the primaries of the multiple-attenuation process can be ameliorated if the amplitude and phase of the multiples model are adjusted, minimizing cross-talk and the differences from the multiples present in the original data in a least-squares sense. To perform this adjustment, we use simultaneous adaptive matching of primaries and multiples, formulated by \cite{alvarez125}.

\plot{migoff1-filt}{width=6.0in,height=2.75in}{Sigsbee2b shot-profile migration (subsurface-offset) using the cross-correlation imaging condition, after multiple attenuation. The front face corresponds to the image at zero subsurface-offset, and the side face corresponds to the subsurface-offset gather at  $x=36800$ ft. Compare with Figure \ref{fig:migoff1}. [{\bf ER}]}

After pre-processing, the multiples have been largely attenuated, making the migration image more suitable for inversion. Figure \ref{fig:migoff1-filt} shows the filtered migration in the subsurface-offset domain. As in Figure \ref{fig:migoff1}, the front face corresponds to the image at zero subsurface offset, and the side face corresponds to the subsurface-offset gather at a position inside the shadow zone ($x=36800$ ft). The angle domain corroborates this hypothesis (Figure \ref{fig:migang1-filt}). The front face corresponds to the image at the $0^\circ$ reflection angle, and the side face corresponds to the reflection-angle gather at a position inside the shadow zone ($x=36800$ ft).  

\plot{migang1-filt}{width=6.0in,height=2.75in}{Sigsbee2b shot-profile migration (reflection angle) using the cross-correlation imaging condition, after multiple attenuation. The front face corresponds to the image at the $0^\circ$ reflection angle, and the side face corresponds to the reflection-angle gather at  $x=36800$ ft. Compare with Figure \ref{fig:migang1}. [{\bf ER}]}

Inversion was computed using equation \ref{eqn:lse} for two different right-hand-side vectors. One is the  migrated image at zero subsurface offset without filtering (front panel, Figure \ref{fig:migoff1}), and the other is the  migrated image at zero subsurface offset after filtering (front panel, Figure \ref{fig:migoff1-filt}). 

The results of the inversion are shown in Figures \ref{fig:invz} and \ref{fig:invz-filt} (unfiltered and filtered input), and should be compared with the migrations shown in Figures \ref{fig:migoff1} and \ref{fig:migoff1-filt}, respectively. The inversion with the unfiltered migration is more unstable. This should require the use of a high value of the regularization parameter, which reduces the effectiveness of the inversion outside the shadow zones. 

\sideplot{invz}{width=3.5in}{Inversion in the poststack image domain with the unfiltered migration as the input (equation \ref{eqn:lse})  after seven iterations  of a conjugate-gradient iterative solver and no regularization. [{\bf ER}]}

\sideplot{invz-filt}{width=3.5in}{Inversion in the poststack image domain with a filtered migration as the input to equation \ref{eqn:lse}. [{\bf ER}]}

 The inversion of the filtered migration shows a better convergence behavior. Again, no regularization was applied. Notice the well-collapsed diffractor at (35000,17000), in contrast to the unfiltered version (Figure \ref{fig:invz}). For the latter, conjugate-gradient iterations focused on reducing the multiple energy in the space of the residuals, decreasing the overall efficacy of inversion.

Figures \ref{fig:resinvz} and \ref{fig:resinvz-filt} show the residuals of the inversion with the unfiltered migration and with the filtered migration, respectively, after seven iterations of the conjugate-gradient solver.  Except for the residuals in the salt body, the highest residual amplitudes of the unfiltered inversion correspond to the multiples in the shadow zone. Notice that some amount of the diffractor energy (35000,17000) is still present in the residuals. On the other hand, the residuals of the inversion with the filtered data do not show the diffractor energy. Additionally, the residuals of the fault close to the diffractor are smaller than that of the unfiltered version. Unfortunately the multiples were not completely removed by filtering, so again the inversion procedure  increases their amplitudes.  

Overall, the inversion results after pre-processing have more balanced amplitudes, allowing the continuation of the reflector inside the shadow zones with much-improved kinematics. The reflectors  also gain more vertical and horizontal resolution, particularly seen at the two faults present in the reflectivity model (Figure \ref{fig:refle-new}) and at the diffractor mentioned above. 

\sideplot{resinvz}{width=3.5in}{Residuals of the inversion in the poststack image domain with the unfiltered migration as the input after seven iterations. [{\bf ER}]}
\sideplot{resinvz-filt}{width=3.5in}{Residuals of the inversion in the poststack image domain with a filtered migration as the input after seven iterations. [{\bf ER}]}

\section{Conclusions}

We showed how different modes of multiple reflections present in Sigsbee2b data behave after migration. The multiple reflections we face are not surface-related; therefore, they can still persist after multiple attenuation during processing. As the model used to compute this dataset is based on a realistic geological setting, the same behavior is likely to happen in real data. We explored the separability of primaries from multiples after migration by applying simple $k_x-k_h$ filters. Once discriminated, primaries and multiples were submitted to a simultaneous, adaptive, non-stationary filtering to adjust amplitude and phase and decrease the cross-talk between them. The adapted multiples were then subtracted from the migrated data.
\par
The inversion results show that, if the physics of the wave-propagation is not adequate to account for all the propagation modes in the migrated data, coherent noise -- in our case multiples -- will dominate the residual space, decreasing the efficiency of inversion. It is imperative, therefore, to make data satisfy the physical approximations by means of pre-processing that, in some sense reconciles processing and inversion.



\bibliographystyle{seg}
\bibliography{SEP,ref134}
