\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} 
Angle-domain common image gathers (ADCIGs), created from downward-continuation or reverse time migration, can provide useful lithological and velocity information \cite{Prucha.sep.100.marie1}. In geologically complex areas, poor illumination causes undesirable kinematic effects and amplitude variations along the angle axis \cite{Prucha.sep.103.marie1,Valenciano.sep.124.alejandro1}. 
\par
The subsurface-offset-to-angle transformation consists of a radial trace transform in the Fourier space with some regularization in the angle direction \cite{Sava.sep.103.paul2} or slant-stack in the physical space plus an additional transformation from offset ray-parameter to reflection angle \cite{Prucha.sep.100.marie1}. The regularization, to some extent, can diminish the amplitude variation caused by poor illumination. The more accurate solution to the illumination problem, however, is achieved by computing a regularized least-squares inverse image \cite{Mclapp122} rather than the simply the adjoint (migration). The inverse image problem can be solved either by computing the Hessian implicitly \cite{Mclapp122} or explicitly \longcite{Valenciano.sep.124.alejandro1}, preferentially, in the reflection-angle domain or, without any physically meaningful regularization direction, in the subsurface-offset domain.
\par
In the reflection-angle domain, the inverse image problem which explicitly computes the Hessian can be performed according to two different strategies. First, by computing the angle-domain Hessian. \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 lacked the resolution to determine which angles were more illuminated. Recently, \longcite{fomel01} introduced the theoretical framework of the oriented wave equation, under which computing the angle-domain Hessian could be promising. In the other approach, the angle-domain Hessian can be evaluated by chaining the offset-to-angle operator and the subsurface-offset Hessian \cite{Valenciano.sep.123.alejandro1}. \longcite{Valenciano.sep.131.alejandro1}, in this report, shows good results obtained by applying this strategy in Sigsbee dataset.
\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 a weighted transformation from subsurface-offset to reflection-angle domain, and finally illustrate the technique 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 subsurface-offset domain common image gathers (SODCIGs), $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, $z$ is depth, $\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 that 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 whose value depends on the reflection coefficient, illumination and focusing, $f$ is the depth domain representation of the seismic pulse, and $z_r$ is the reflector depth. The fact that $A$ and $z_r$ are functions of $h$ accommodates the focusing of reflector amplitudes at nonzero-subsurface offsets because of inaccuracies in migration velocity and problems in illumination. A SODCIG containing several reflectors can be described by the superposition of individual reflectors, each described by 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 in 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 stationary-phase method. 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$. The term $sgn(\phi''(c))$ corresponds to the signal of the 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 \cite{Tang.sep.129.yaxun1}, there will be as many stationary points as 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 point contributions. The other special case is when all the energy is focused at zero subsurface-offset, indicating 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 two images migrated with slightly different weights. \longcite{tygel} applied the same ideas to what they called a multiple-weight diffraction stack to obtain the stationary point location that in turn, along with source and receiver position, specifies the reflection ray. For instance, to estimate reflector dips one can compute two different migrated images, $M_a$ and $M_1$, using two distinct migration-weighting functions, say the migration angle ($M_a$) and simply a constant value of one ($M_1$). For $M_a$, the resulting amplitudes are weighted by the migration angles around the stationary point. In this region the migration operator and reflectors are tangent. Consequently the local 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
The phase behavior of the offset-to-angle transformation shows that the main contribution for the image in the angle domain comes from the vicinity of the stationary point. Therefore, 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 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, $H(z,h)$ ---  defined in the subsurface-offset domain. The amplitudes of the resulting ADCIGs, according to the stationary phase results in the previous section, can be represented by
\begin{equation}
\tilde{Q}(z,\gamma) \sim H(z,h^{\star}) \frac{A(h^{\star})} {\sqrt{\left | \Phi''(h^{\star})  \right |}} F(z-\Phi(\gamma,h^{\star})) .
\label{eq:stat4}
\end{equation}
where $H(z,h^{\star})$ are the averaged values of the subsurface-offset Hessian diagonals in the vicinity of the stationary point.
The ADCIGs, $\tilde{Q}(z,\gamma)$, are to be divided by the non-weighted transformed results from equation \ref{eq:stat3}, $Q(z,\gamma)$, using a regularization term, $\epsilon$, in the denominator to avoid division by small numbers. Finally, a median filter in the $(z,x)$-plane (represented by $\langle \rangle _{(z,x)}$) for every angle section is applied to remove spurious amplitudes, thus providing an estimate of subsurface-offset parameter in the angle domain. 
The angle-transformed Hessian diagonals, $H(z,\gamma)$ are computed according to equation \ref{eq:hesang}:
\begin{equation}
H(z,\gamma) \sim \langle \frac{\tilde{Q}Q}{Q^2+\epsilon}\rangle_{(z,x)} .
\label{eq:hesang}
\end{equation}
\par
The general formula of the subsurface-offset 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_s$ denotes Green's function from the source point, $\bf{x_s}$, to the image point, $\bf{x}$, and $G_r$ the Green's function from the image point to the receiver point, $\bf{x_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.
\par
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 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. However, this procedure does not correct for kinematic errors of the migrated image and, depending on the complexity of the illumination pattern, only the least-squares inverse image may be able to 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}
In this case, the diagonals just represent the cross-correlation between Green's functions computed for a specific image point but shifted by different subsurface-offsets. The off-diagonal terms of the subsurface-offset Hessian are an expression of how much illumination for a specific subsurface-offset is conditioned by the illumination in another subsurface-offset. For the ideal case of infinite cable length, infinite frequency bandwidth and constant velocity, the subsurface-offset Hessian is the identity operator meaning that subsurface-offsets are linearly independent. In this extreme situation, all the energy in a SODCIG will be concentrated at zero-subsurface offset. Therefore, in general, the main diagonal represents an estimate of how much illumination for a specific subsurface-offset is not conditioned by the illumination in another subsurface-offset. Consequently, it is straightforward to consider the main diagonal as the natural candidate to be transformed to angle domain.
\par
Although, in principle, any diagonal of the subsurface-offset Hessian can be transformed to the angle domain by the proposed approach, at present I have conceived a direct application only for the transformed main diagonal, which is to use it as a weight to balance the amplitudes of the ADCIGs. In the next section I show 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}
To test this methodology I applied it on the well known Sigsbee synthetic dataset. This dataset presents illumination problems due to an irregular salt body shape, which results in unbalanced amplitude patterns in the seismic section (Figure \ref{fig:Sigsbee}). The small rectangle in Figure \ref{fig:Sigsbee} highlights the target area. The off-end acquisition geometry consists of 348 receivers, 75 ft apart, resulting in 26025 ft maximum offset. As source coordinates are smaller than receiver coordinates, the source-receiver offsets are positive. Therefore, the energy is mainly distributed at positive reflection angles.
\plot{Sigsbee}{width=3.0in}{Shot profile migration of part of Sigsbee dataset --- zero-subsurface offset. The small box highlights the target area.}
Figures \ref{fig:Ojoin16} and \ref{fig:Ojoin50} show two different SODCIGs and their respective subsurface-offset Hessian main diagonal, located at CMP coordinates 33200 ft and 35500 ft, respectively. The illumination problem gets more severe as we approach the dipping salt flank. In this work, all figures related to illumination show high-illumination values in dark gray and low illumination values in light gray. Both SODCIGs clearly exhibit the effects of poor illumination represented by horizontal and dipping ($\sim$ 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 35500 ft, events curving upward are the expression of multiples. 
\par
\plot{Ojoin16}{width=3.0in}{SODCIG and diagonal of the subsurface-offset Hessian at CMP coordinate 33200 ft. Note the effects of poor illumination represented by horizontal and dipping straight events in the SODCIG.}
\plot{Ojoin50}{width=3.0in}{SODCIG and diagonal of the subsurface-offset Hessian at CMP coordinate 33200 ft. Upward curved events correspond to  multiples. Note the effects of poor illumination represented by horizontal and dipping straight events in the SODCIG.}

\subsection{Transformation of the Main diagonal}
Figures~\ref{fig:join16} and ~\ref{fig:join50} show that, for ADCIGs at CMP coordinates 33200 ft and 35500 ft, the energy is mainly focused at reflection angles around 0$^{\circ}$ and 40$^{\circ}$ -- 50$^{\circ}$. 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 only the diagonal of the Hessian is being used, the kinematic artifacts remain. 
\plot{join16}{width=6.0in}{ADCIGs and diagonal of the transformed angle-domain Hessian at CMP coordinate 33200 ft. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c). Note the improved amplitude balance in the angle direction.}
\plot{join50}{width=6.0in}{ADCIGs and diagonal of the transformed angle-domain Hessian at CMP coordinate 35500 ft. 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 strength of the events in the ADCIGs. However, as shown in the next example, it yields useful information about illumination. Figures~\ref{fig:join00}, ~\ref{fig:join15} and ~\ref{fig:join30} show 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. However, 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, to some extent, that the proposed approach can yield reliable information about illumination despite the presence of high-amplitude events not predicted in the computation of the Green's functions.
 
\plot{join00}{width=6.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=6.0in}{15$^{\circ}$-angle section. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c). Notice the improved amplitude balance along the CMP direction.}
\plot{join30}{width=6.0in}{30$^{\circ}$-angle section. Before amplitude compensation (a); diagonal of the Hessian in the angle domain (b); and after amplitude compensation (c).}

Figure~\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. Unfortunately, however, the amplitudes of internal multiples are also increased.
\plot{Stk}{width=6.0in,height=6.0in}{Stack along the angle axis. Before amplitude compensation (a) and after (b).}

\subsection{Transformation of the off-diagonal terms}
In the final example, I show, for the zero-angle section, off-diagonal terms after the transformation to the angle domain. As previously mentioned, it is still not clear the way this results can be used to correct for illumination problems. According to the initial interpretation of the results shown 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. All the transformed diagonals were scaled to the same value. Actually, the rms value of the main diagonal is approximately 20 and 180 times higher than the rms value of the 5th and 15th off-diagonal, respectively.

\plot{jDoff01}{width=6.0in,height=8.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 the subsurface-offset domain using what I call weighted offset-to-angle, particularly subsurface-offset Hessian diagonals. The proposed approach provides useful information, which 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, it is still not clear how to use these transformed off-diagonal terms in inversion schemes.
\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} 
