\shortnote
\lefthead{Guerra}
\righthead{Weighted offset-to-angle}
\email{guerrac@stanford.edu}
\footer{SEP--131}
\title{Angle-domain parameters computed via weighted slant-stack}
\author{Claudio Guerra}

\vspace{.5cm}

\section{Introduction} 
It is well known that angle-domain common-image gathers (ADGIC) in geologically complex regions can show undesirable kinematic effects and amplitude variations along the angle axis due to illumination problems \cite{Prucha.sep.103.marie1,Valenciano.sep.124.alejandro1}.
 The subsurface-offset-to-angle transformation includes some regularization of the linear interpolation problem \cite{Prucha.sep.103.marie1,Sava.sep.103.paul2} along the angle axis which, in some extent, diminish the amplitude variations not related to the interpolation problem. However, the correct solution to the amplitude variation problem is achieved by computing a least-squares inverse image \cite{Mclapp122}. The inverse problem can be solved iteratively by applying the migration operator and its transpose every iteration, or according to the  approach proposed by \longcite{Valenciano.sep.124.alejandro1} in which the Hessian is computed in a target-oriented way and used an inversion problem solved iteratively. This problem can be solved both in the subsurface-offset domain or in the reflection-angle domain.
\par
To solve the inversion problem in the reflection-angle domain it is necessary to compute the Hessian in this domain or to chain the offset-to-angle operator and the subsurface-offset Hessian \cite{Valenciano.sep.123.alejandro1}. \longcite{Valencia125} proposed to obtain the angle-domain Hessian by applying the slant-stack technique to compute ADCIGs on the subsurface-offset Hessian. They noticed that the resulting angle-domain Hessian for a model with a Gaussian velocity anomaly "\textit{ ... lacks of resolution to be able to interpret which angles get more illumination.}" As Bob Clapp mentioned in one of our seminars, this lack of resolution can be explained by the fact that integration of a smooth function yields an even more smooth result. Recently, \longcite{fomel01} introduced the theoretical framework of the oriented wave equation under which to compute the angle-domain Hessian can be promissing.
\par
Here, I propose a general framework to map any information computed in the subsurface-offset domain to the angle-domain. The proposed approach relies on the asymptotic nature of the slant-stack transformation from subsurface-offset to angle domain. I first show the validity of the stationary phase assumption for the offset-to-angle transformation; then describe the weighted transformation from subsurface-offset to reflection-angle domain; and finally exemplify with the transformation of the diagonal of the Hessian in the subsurface-offset domain to the angle domain, yielding amplitude factors to compensate for illumination problems in ADCIGs. Additionally, I show the transformation of some off-diagonal terms which, at present, does not have a direct application in the amplitude correction problem.

\section{Phase behavior of the offset-to-angle transformation}

The offset-to-angle transformation can be expressed by the integration of the SODCIG, $P(z,h)$, along a certain slanted path, according to the equation  
\begin{equation}
Q(z,\gamma) = \int_{A} \varrho[P(z,h)] dh \arrowvert_{z = \zeta(\gamma,h) } ,
\label{eq:int1}
\end{equation}
where $Q(z,\gamma)$ is the output ADCIG, $\gamma$ is the aperture angle, $h$ is the subsurface offset, $\varrho$ is the {\textit {rho}}-filter which aims to yield the correct phase of the output ADCIG \cite{iei}, $A$ is the domain of integration which defines the range of subsurface offsets to be summed and $\zeta(\gamma,h)$ is the slanted path given by 
\begin{equation}
\zeta(\gamma,h) = z_0 + h \tan \gamma ,
\label{eq:slant}
\end{equation}
where $z_0$ is the depth coordinate at zero subsurface offset. A single reflector in a SODCIG can be represented by
\begin{equation}
P(z,h) = A(h) f(z - z_r(h)), 
\label{eq:reflector}
\end{equation}
where $A(h)$ is an amplitude term which value depends on the reflection coefficient, illumination and focusing, $f$ is the seismic pulse and $z_r$ is the reflector depth. The fact that $A$ and $z_r$ are function of $h$ accommodates reflector amplitudes being focused at subsurface offsets different from zero due to inaccuracies in migration velocity and problems in illumination. A SODCIG containing several reflectors can be described by the superposition of individual reflectors like the equation \ref{eq:reflector}.
\par
Equation \ref{eq:int1} Fourier transformed to the $k_z$ domain after inserting equation \ref{eq:reflector} reads 
\begin{equation}
\hat{Q}(k_z,\gamma) = \sqrt{\frac{-ik_z}{2\pi}} F(k_z) \int_{-h}^h A(h) e^{-ik_z \Phi(\gamma,h)}dh ,
\label{eq:int2}
\end{equation}
where $\Phi(\gamma,h) = \zeta(\gamma,h) - z_r(h)$ is the phase function. Assuming that $A(h)$ is not itself an oscillating function and considering the high-$k_z$ regime, the argument of the integral in equation \ref{eq:int2} rapidly oscillates yielding negligible amplitudes for integration over a full period, except for the case where the phase function, $\Phi(\gamma,h)$, remains stationary. This condition is achieved at the vicinity of a point -- the stationary point -- in the SODCIG with a certain subsurface offset, $h^{\star}$, where $\zeta(\gamma,h)$ is tangent to $z_r(h)$, or 
\begin{eqnarray*}
\zeta(\gamma,h) &=& z_r(h) \\
&&\\
\frac{\partial \zeta(\gamma,h)}{\partial h} &=& \frac{\partial z_r(h)}{\partial h} ,\\
\end{eqnarray*}
estimated in $h = h^{\star}$.
\par
Equation \ref{eq:int2} can be evaluated by the method of the stationary phase. According to \longcite{blei02}, under the assumption of a single stationary point in which the second derivative does not vanish, integrals like 
\begin{equation}
I(\lambda) = \int_A f(t) e^{i\lambda \phi(t)}dt ,
\label{eq:int3}
\end{equation}
where $f(t)$ is a smooth and compact function, can be asymptotically approximated by
\begin{equation}
I(\lambda) \sim e^{i(\lambda \phi(c) + sgn(\phi''(c))\frac{\pi}{4})} f(c) \sqrt{\frac{2\pi}{\lambda \left | \phi''(c)  \right | }} ,
\label{eq:stat1}
\end{equation}
if $\lambda \rightarrow \infty$. $sgn(\phi''(c))$ corresponds to the signal of second derivative of the phase function, $\phi(t)$, evaluated at the stationary point, $c$.
\par
It turns out that the stationary phase formula of equation \ref{eq:int2} is given by
\begin{equation}
\hat{Q}(k_z,\gamma) \sim \frac{A(h^{\star})} {\sqrt{\left | \Phi''(h^{\star})  \right |}} F(k_z) e^{-ik_z \Phi(\gamma,h^{\star})} .
\label{eq:stat2}
\end{equation}
Finally, the inverse Fourier transform of equation \ref{eq:stat2} gives
\begin{equation}
Q(z,\gamma) \sim \frac{A(h^{\star})} {\sqrt{\left | \Phi''(h^{\star})  \right |}} F(z-\Phi(\gamma,h^{\star})) .
\label{eq:stat3}
\end{equation}
Equation \ref{eq:stat3} shows that the main contribution to the amplitudes in the ADCIG comes from the vicinity of the stationary point. The second derivative of the phase function with respect to $h$ is basically the second derivative of $z_r$, as $\zeta(\gamma,h)$ is a straight line. If $z_r$ is a straight event in the SODCIG, meaning that just a very small range of angles has been illuminated, as in the case of a single shot and well-sampled receiver wavefield as discussed by \longcite{tang1}, there will be as much stationary points as the number of subsurface offsets. In this situation, the integration interval is divided in such a way that each new interval contains only one stationary point, and the final result is the sum of all individual stationary-points contribution. The other special case is when all the energy is focused at zero subsurface-offset corresponding to a situation of good illumination for all reflection angles and migration with the correct velocity. It is a generalization of the previous case and is solved in the same way for various illumination angles. 

\section{Weighted offset-to-angle transformation}
\longcite{blei01} describes a strategy to estimate parameters from the subsurface using different images migrated with two slightly different weights. \longcite{tygel} applied the same ideas to what they called multiple-weights diffraction stack to obtain the stationary point location that in turn, along with source and receiver position, specifies the reflection ray. For instance, if one desires to estimate reflector dips, a possible way to do that is by computing two different migrated images, $M_a$ and $M_1$ obtained with two distinct migration-weighting functions, say migration angle ($M_a$ result) and simply one ($M_1$ result). For $M_a$, the resulting amplitudes are weighted by the migration angles around the stationary point. As in this region the migration operator and reflectors are tangent, the average of the migration angle is an estimate of the reflector dip. So, the division $M_a/M_1$ results in an estimate of the dip of the reflectors.
\par
After analyzing the phase behavior of the offset-to-angle transformation in the previous section and taking into account that the main contribution for the image in the angle domain comes from the vicinity where the phase function remains stationary in the subsurface offset domain, the use of the weighted stacking strategy \cite{blei01} to map quantities computed in the subsurface-offset domain to the angle domain is straightforward. The mapping of certain attributes can be useful, for instance, to balance amplitudes in the angle domain.
\par
In the following, the aim of the weighted offset-to-angle transformation is to compute weights to be applied on ADCIGs in such a way that amplitude variations due to illumination problems are attenuated. The weighted offset-to-angle transformation is represented by the computation of ADCIGs from SODCIGs previously multiplied by some parameter -- in the present case, the subsurface-offset Hessian diagonals --  defined in the subsurface-offset domain. These ADCIGs are to be divided by the non-weighted transformed results, using a regularization term to avoid division by small numbers. Finally, a median filter is applied to filter spurious amplitudes out, thus providing an estimate of that parameter in the transformed domain. 
\par
The general formula of the Hessian in the prestack-inversion problem is   
\begin{equation}
\bf{H(x,h;x',h')} = \sum_{\omega} \sum_{x_{s}} \bf{G^{\ast}_s(x+h,x_s};\omega)\bf{G_s(x'+h',x_s};\omega) \sum_{x_r} \bf{G^{\ast}_r(x-h,x_r};\omega) \bf{G_r(x'-h',x_r};\omega)
\label{eq:Hes1}
\end{equation}
where $G$ denotes Green's functions from source -- $\bf{x_s}$ -- to the image point -- $\bf{x}$ -- (subscript $s$) and from the image point to receiver -- $\bf{x_r}$ -- (subscript $r$), $\bf{h}$ is the subsurface-offset, the prime indicates points in the image space in the vicinity of the image point and different subsurface-offsets and the $^{\ast}$ stands for the conjugate transpose of the Green's functions. In equation \ref{eq:Hes1}, the Green's functions, computed for every source and receiver, at different image points and for different subsurface-offsets are cross-correlated. The main diagonal of the Hessian, which is the Laplacian of the cost function related to the model parameters, contains the autocorrelation of the Green's functions and, generally, carries most of the information about the illumination. Sometimes, a good and cheap solution is just to approximate the Hessian by its main diagonal and apply its inverse to the migrated image. But, this procedure does not correct for kinematic errors of the migrated image and, depending on how complex is the illumination pattern, only the least-squares inverse image can provide reasonable results \cite{Mclapp122}. 
\par
Equation \ref{eq:Hes2} shows the structure of the subsurface-offset domain Hessian used in the examples. 
\begin{equation}
\bf{H(x,h;h')} = \sum_{\omega} \sum_{x_{s}} \bf{G^{\ast}_s(x+h,x_s};\omega)\bf{G_s(x+h',x_s};\omega) \sum_{x_r} \bf{G^{\ast}_r(x-h,x_r};\omega) \bf{G_r(x-h',x_r};\omega)
\label{eq:Hes2}
\end{equation}
As one can notice, the diagonals represent just the cross-correlation between Green's functions computed at the same image point for different subsurface-offsets. Although, in principle, any diagonal of the subsurface-offset Hessian of such a strucuture can be transformed to the angle domain by the proposed approach, at present I have only conceived a direct application for the transformed main diagonal. 
\par
In the next section it will be shown examples of the angle-domain transformed subsurface-offset Hessian diagonals, as well as the comparison of migrated images before and after the amplitude compensation with the transformed main diagonal, for a small portion of the Sigsbee dataset.

\section{Examples}
The Sigsbee synthetic dataset was used in this work. This well known dataset presents illumination problems due to irregular salt body shape which manifest as unbalanced amplitude patterns in the seismic section. The small rectangle in figure \ref{fig:Sig} highlights the target area. The off-end acquisition geometry consists of 348 receivers, 75ft apart, resulting in 26025ft maximum offset. According to the CMP coordinates in Figure \ref{fig:Sig}, the source coordinates are lower than receiver coordinates, defining positive source-receiver offsets. Therefore, the energy will be mainly distributed at positive reflection angles.
\par
Figures \ref{fig:Ojoin16} and \ref{fig:Ojoin50} shows two different SODCIGs and their respective subsurface-offset Hessian main diagonal, located at CMP coordinates 33200ft and 35500ft, respectively. The lower CMP position is closer to the tip of the salt and the higher is closer to the salt body. In this work, all figures related to illumination presents in dark grey high illumination values and in light grey low illumination values.
Both SODCIGs clearly exhibt the effects of poor illumination represented by horizontal and dipping ($\approx$ 40$^{\circ}$ -- 50$^{\circ}$) straight events. The energy smeared along these directions will be mapped to the reflection-angle domain according the dips observed in the subsurface-offset domain. In the SODCIG at CMP position 35500ft events curving upward are the expression of internal multiples originated at the top and base of the salt body. 
\par
\plot{Ojoin16}{width=2.0in}{SODCIG and diagonal of the subsurface-offset Hessian at CMP coordinate 33200ft. Note the effects of poor illumination represented by horizontal and dipping straight events in the SODCIG.}
\plot{Ojoin50}{width=2.0in}{SODCIG and diagonal of the subsurface-offset Hessian at CMP coordinate 33200ft. Upward curved events correspond to internal multiples originated at the top and base of the salt body. Note the effects of poor illumination represented by horizontal and dipping straight events in the SODCIG.}
\par
As already predicted, the ADCIGs at CMP coordinates 33200ft and 35500ft shows focusing of energy at reflection angle around 0$^{\circ}$ and 40$^{\circ}$ -- 50$^{\circ}$, as can be seen in Figures~\ref{fig:join16} and ~\ref{fig:join50}. Additionally, these figures show the original reflection-angle gather (a), the main diagonal of the Hessian transformed to the reflection-angle domain (b) and the amplitude compensated ADCIG (c). Notice how the amplitudes are better distributed along the reflection-angle axis after compensation by the inverse of the diagonal of the Hessian. However, as just the diagonal of the Hessian is being used, the kinematic artifacts remain. 
\plot{join16}{width=3.0in}{ADCIGs and diagonal of the transformed angle-domain Hessian at CMP coordinate 33200ft. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c). Note how better balanced the amplitude is in the angle direction.}
\plot{join50}{width=3.0in}{ADCIGs and diagonal of the transformed angle-domain Hessian at CMP coordinate 35500ft. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c).}
\par
The proposed approach seems to be dependent on the amplitude strenght of the events in the ADCIGs. However, as can be seen in the next example, it yields a useful information  about illumination. In Figures~\ref{fig:join00}, ~\ref{fig:join15} and ~\ref{fig:join30} can be seen angle sections of the original angle data (a), the main diagonal of the Hessian in the angle domain (b) and the amplitude balanced angle data (c). Again, the amplitude compensation proved to be effective. But, notice how for the zero-angle section the illumination computed in the angle domain is low at the right part of the section, in spite of the high amplitudes of the internal multiples. This confirms, in some extent, that the proposed approach can yield reliable information about illumination even though high amplitudes of events not predicted in the computation of the Green's functions are present.
 
\plot{join00}{width=3.0in}{Zero-angle section. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c). Notice the low illumination in the right part of the section, near the flank of the salt body.}
\plot{join15}{width=3.0in}{15$^{\circ}$-angle section. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c). Notice how better balanced is the amplitude along the CMP direction.}
\plot{join30}{width=3.0in}{30$^{\circ}$-angle section. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c).}

Figures~\ref{fig:Stk} shows the stacked section along the angle axis, before (a) and after(b) the amplitude compensation by the inverse of the diagonal of the Hessian in the angle domain. The dimming of the amplitudes at the right portion of the section is almost eliminated. However, unfortunately, the amplitudes of internal multiples are increased too.
\plot{Stk}{width=3.0in}{Stack along the angle axis. Before amplitude compensation (a) and after (b).}

In the last example, I show, for the zero-angle section, off-diagonal terms after the transformation to angle domain. In Figure~\ref{fig:jDoff01} it is clear that off-diagonal terms start gaining importance as we approach the flank of the salt body at the right of the section. 

\plot{jDoff01}{width=3.0in}{Zero-angle section(a); main diagonal (b);5th subsurface-offset domain off-diagonal transformed to angle (c); and 15th subsurface-offset domain off-diagonal transformed to angle (d).}

\section{Conclusions}
I showed how to estimate angle-domain parameters from subsurface-offset domain using what I call weighted offset-to angle, particularly subsurface-offset Hessian diagonals. The proposed approach provides useful information, that can be confirmed by the amplitude compensation results. The transformation of off-diagonal terms of the subsurface-offset Hessian indicates that the results are not strongly dependent on the amplitude distribution in the ADCIGs. However, how to use these transformed off-diagonal terms in some inversion scheme is still not clear. 
\section{Aknowledgements}
I would like to thank Alejandro Valenciano for providing me the subsurface-offset domain Hessian and common-image gathers and for the explanations and fruitful discussions.
 
\bibliographystyle{sep}
\bibliography{SEP,ref1}