
\def\dtdhxx{\frac{\partial t}{\partial h_{x_1}}}
\def\dtdxz {\frac{\partial t}{\partial x_3}}
\def\dtdhx {\frac{\partial t}{\partial h_{x_1}}}
\def\dtdhz {\frac{\partial t}{\partial h_{x_3}}}
\def\dtdhex{\frac{\partial t}{\partial h_{\xi_1}}}
\def\dtdhez{\frac{\partial t}{\partial h_{\xi_3}}}

\def\dtdxx{\frac{\partial t}{\partial x_1}}
\def\dtdxz{\frac{\partial t}{\partial x_3}}
\def\dtdex{\frac{\partial t}{\partial \xi_1}}
\def\dtdez{\frac{\partial t}{\partial \xi_3}}
\def\dzdh{\frac{\partial z}{\partial h}}


% e over x
\def\dexdxz{\frac{\partial \xi_1}{\partial x_3}}
\def\dexdxx{\frac{\partial \xi_1}{\partial x_1}}
\def\dezdxz{\frac{\partial \xi_3}{\partial x_3}}
\def\dezdxx{\frac{\partial \xi_3}{\partial x_1}}

% x over e
\def\dxxdex{\frac{\partial x_1}{\partial \xi_1}}
\def\dxzdex{\frac{\partial x_3}{\partial \xi_1}}
\def\dxxdez{\frac{\partial x_1}{\partial \xi_3}}
\def\dxzdez{\frac{\partial x_3}{\partial \xi_3}}

% hx over he
\def\dhxxdhex{\frac{\partial h_{x_1}}{\partial h_{\xi_1}}}
\def\dhxxdhez{\frac{\partial h_{x_1}}{\partial h_{\xi_3}}}
\def\dhxzdhex{\frac{\partial h_{x_3}}{\partial h_{\xi_1}}}
\def\dhxzdhez{\frac{\partial h_{x_3}}{\partial h_{\xi_3}}}

% he over hx
\def\dhexdhxx{\frac{\partial h_{\xi_1}}{\partial h_{x_1}}}
\def\dhexdhxz{\frac{\partial h_{\xi_1}}{\partial h_{x_3}}}
\def\dhezdhxx{\frac{\partial h_{\xi_3}}{\partial h_{x_1}}}
\def\dhezdhxz{\frac{\partial h_{\xi_3}}{\partial h_{x_3}}}

% hx over e
\def\dhxxdex{\frac{\partial h_{x_1}}{\partial {\xi_1}}}
\def\dhxxdez{\frac{\partial h_{x_1}}{\partial {\xi_3}}}
\def\dhxzdex{\frac{\partial h_{x_3}}{\partial {\xi_1}}}
\def\dhxzdez{\frac{\partial h_{x_3}}{\partial {\xi_3}}}

% he over x
\def\dhexdxx{\frac{\partial h_{\xi_1}}{\partial {x_1}}}
\def\dhexdxz{\frac{\partial h_{\xi_1}}{\partial {x_3}}}
\def\dhezdxx{\frac{\partial h_{\xi_3}}{\partial {x_1}}}
\def\dhezdxz{\frac{\partial h_{\xi_3}}{\partial {x_3}}}

% hx over x
\def\dhxxdxx{\frac{\partial h_{x_1}}{\partial {x_1}}}
\def\dhxxdxz{\frac{\partial h_{x_1}}{\partial {x_3}}}
\def\dhxzdxx{\frac{\partial h_{x_3}}{\partial {x_1}}}
\def\dhxzdxz{\frac{\partial h_{x_3}}{\partial {x_3}}}

% he over e
\def\dhexdex{\frac{\partial h_{\xi_1}}{\partial {\xi_1}}}
\def\dhexdez{\frac{\partial h_{\xi_1}}{\partial {\xi_3}}}
\def\dhezdex{\frac{\partial h_{\xi_3}}{\partial {\xi_1}}}
\def\dhezdez{\frac{\partial h_{\xi_3}}{\partial {\xi_3}}}

% x over he
\def\dxxdhex{\frac{\partial x_1}{\partial h_{\xi_1}}}
\def\dxxdhez{\frac{\partial x_1}{\partial h_{\xi_3}}}
\def\dxzdhex{\frac{\partial x_3}{\partial h_{\xi_1}}}
\def\dxzdhez{\frac{\partial x_3}{\partial h_{\xi_3}}}

% x over hx
\def\dxxdhxx{\frac{\partial x_1}{\partial h_{x_1}}}
\def\dxxdhxz{\frac{\partial x_1}{\partial h_{x_3}}}
\def\dxzdhxx{\frac{\partial x_3}{\partial h_{x_1}}}
\def\dxzdhxz{\frac{\partial x_3}{\partial h_{x_3}}}

% e over he
\def\dexdhex{\frac{\partial \xi_1}{\partial h_{\xi_1}}}
\def\dexdhez{\frac{\partial \xi_1}{\partial h_{\xi_3}}}
\def\dezdhex{\frac{\partial \xi_3}{\partial h_{\xi_1}}}
\def\dezdhez{\frac{\partial \xi_3}{\partial h_{\xi_3}}}

% rm font
\def\rtan{{\rm tan}}
\def\rcos{ {\rm cos} }
\def\rsin{ {\rm sin} }
\def\rsinh{ {\rm sinh} }
\def\rcosh{ {\rm cosh} }
\def\w{\omega}

\chapter{Generalized-coordinate ADCIGs} 
%\author{Jeff Shragge}

\lefthead{Shragge}
%\righthead{ADCIGs in generalized coordinates}
%\footer{SEP-134}


%\maketitle

%\author{Jeff Shragge}

%\begin{abstract}
The theory of angle-domain common-image gathers (ADCIGs) is extended
to migrations performed in generalized 2D coordinate systems.  I develop an expression linking the definition of reflection opening angle to generalized coordinate various geometric factors.  I demonstrate that generalized coordinate ADCIGs can be calculated directly using Fourier-based offset-to-angle approaches for coordinate systems satisfying the Cauchy-Riemann conditions of differentiability.  The canonical examples of tilted Cartesian, polar, and elliptic coordinates are used to illustrate the ADCIG theory.  I  compare analytically and numerically generated image volumes for a set of elliptically shaped reflectors.  Experiments with the BP velocity synthetic data set illustrate that elliptic-coordinate ADCIGs better-resolve the reflection opening angles of steeply dipping structure relative to Cartesian image volumes,  due to better large-angle propagation and improved sensitivity to steep structural dips afforded by coordinate transformations.
%\end{abstract}

\section{Introduction}
%
Angle-domain common image gathers (ADCIGs) are used increasingly in seismic imaging to examine migration velocity model accuracy \cite[]{3DSI}.  The key idea is that migrating with the correct velocity model leads to flat angle gathers that shift neither vertically nor horizontally as a function of reflection opening angle. Migrating with an incorrect velocity, though, leads to inconsistent angle-domain reflectivity and generates residual curvature in the ADCIG volume.  ADCIGs are thus an effective velocity analysis tool and have been incorporated in wave-equation-based inversion schemes to update velocity profiles \cite[]{SavaWEMVA1, SavaWEMVA2}.    
%
\par
%
Wave-equation imaging techniques generate ADCIGs in straightforward manners for both shot-profile \cite[]{deBruin,Rickett,SavaADCIG} and shot-geophone \cite[]{PruchaADCIG, MosherADCIG} migration approaches.  In shot-profile migration, one first generates a subsurface offset axis at each depth step by correlating the source and receiver wavefields at a number of subsurface shifts.  The second step involves computing an offset-to-angle domain transformation, for example, using post-imaging Fourier-based operators \cite[]{SavaADCIG}. 
%
\par
%
Conventional ADCIG theory usually assumes horizontally and uniformly spaced wavefield shifts, largely because wavefield extrapolation and  imaging are most commonly performed in Cartesian coordinates.  However, a number of studies have noted that ADCIG results degrade for steeply dipping structure, such as salt flanks \cite[]{BiondiShan2002,BiondiSymes}.  Although this is partially due to problems associated with inaccurate large-angle extrapolation, ADCIGs calculated using horizontal wavefield shifts become decreasingly sensitive for increasingly steep structural dips.  This problem can be partially addressed by generating ADCIGs with vertical offset gathers; however, this approach is less desirable computationally because it requires storing larger wavefield volumes in memory.
%
\par
%
The introduction of shot-profile migration in more general coordinate systems, including tilted Cartesian \cite[]{SEG2004sgj} and elliptic meshes \cite[]{ShraggeGEORWE}, presents an opportunity to circumvent problems associated with generating ADCIGs for steeply dipping structure.  In particular, coordinate systems can be oriented such that geologic structures appearing steep in Cartesian meshes become relatively small angle, thus improving the robustness of the ADCIG calculation.  Developing an ADCIG theory able to handle more arbitrary geometries and wavefield shifts, though, requires properly handling the effects of non-Cartesian geometry.  For example, extrapolation in non-Cartesian coordinate systems induces local stretches, rotations and/or shearing of wavefields \cite[]{ShraggeGEONO}.  Similarly, non-uniform wavefield sampling can lead to anisotropic angle-domain stretching.  These effects can be accounted for using Jacobian change-of-variable transformations.   
%
\par
%
The goal of this paper is to extend ADCIG theory to non-Cartesian geometries.  I demonstrate that ADCIG theory, as developed in a differential sense \cite[]{SavaADCIG}, remains valid for arbitrary geometry provided that the corresponding derivative operators are properly specified.  Non-Cartesian coordinates do, though, introduce space-domain geometric factors that can render Fourier-based offset-to-angle methods unsuitable. However, I show that ADCIGs can be calculated directly in the Fourier domain for all coordinate systems that satisfying the Cauchy-Riemann differentiability criteria \cite[]{Arfken}.
%
\par
%
I begin by discussing how to generate subsurface offsets and ADCIGs in Cartesian coordinates.  I then provide an extension to generalized coordinate systems based on Jacobian change-of-variable arguments.  I examine two canonical coordinate systems, tilted Cartesian and elliptic meshes, where the reflection angle can be explicitly calculated using Fourier-based methods, and a third, polar coordinates, where it cannot.  I test the generalized ADCIG theory analytically using a set of elliptic reflectors, and demonstrate how computing angle gathers in elliptic coordinates can lead to improvements in ADCIG image volumes relative to Cartesian coordinates, especially for steeply dipping structure.
% 
\section{ADCIG theory}
%
The ADCIG theory presented in this section draws from that presented in \cite{SavaADCIG}.  In the ensuing development, $\boldsymbol{x}=[x_1,x_3]$ denotes the Cartesian variables and  $\boldsymbol{\xi}=[\xi_1,\xi_3]$ represents a generalized Riemannian coordinate system.   I also use a convention where the extrapolation axis is oriented in the $x_3$ ($\xi_3$) direction for Cartesian (Riemannian) coordinates.  Coordinates $x_2=\xi_2=0$ are assumed throughout.
%
\subsection{Generating  subsurface offsets axes}
%
Shot-profile migration in Cartesian coordinates consists of iteratively completing a two-step procedure.  The first step involves propagating the source and receiver wavefields, $S$ and $R$, from depth level $x_3-\Delta x_3$ to  $x_3$ using an extrapolation operator $E_{x_3}[\cdot]$
\begin{eqnarray}
E_{x_3}[S(x_3-\Delta x_3,x_1|\omega)] &=& S(x_3,x_1|\omega), \nonumber \\
E^*_{x_3}[R(x_3-\Delta x_3,x_1|\omega)] & = & R(x_3,x_1|\omega),
\end{eqnarray}
where $*$ denotes complex conjugate, $\omega$ is angular frequency, and $\Delta x_3$ is the depth step. A subsurface image, $I$, is subsequently computed at each $x_3$ extrapolation step by evaluating an imaging condition
\begin{equation} \label{eqn:subsurface1}
I(x_3,x_1,h_{x_1}) = \sum_\omega R^*(x_3,x_1+h_{x_1}|\omega) S(x_3,x_1-h_{x_1}|\omega).
\end{equation}
where the subsurface offset axis, $h_{x_1}$, in equation~\ref{eqn:subsurface1} is generated by correlating the source and receiver wavefields at various relative shifts in the $x_1$ direction.  Finally, the ADCIG volume is computed using an offset-to-angle transformation operator, ${\bf T}_{h_{x_1} \rightarrow \gamma}$
\begin{equation} \label{eqn:subsurface2}
I(x_3,x_1,\gamma) = {\bf T}_{h_{x_1}\rightarrow \gamma} I(x_3,x_1,h_{x_1}),
\end{equation} 
where $\gamma$ is the reflection opening angle discussed below.
%
\par
%
Imaging in generalized coordinate systems follows the same two-step procedure.   However, because of the different migration geometry in the $\boldsymbol{\xi}$-coordinate system, new extrapolation operators, $E_{\xi_3}[\cdot]$, must be used to propagate wavefields.  I specify these operators using Riemannian wavefield extrapolation (RWE).  RWE is not discussed herein, and readers interested in additional information are referred to \cite{Sava.geo.rwepub} and \cite{ShraggeGEONO}.  
%
\par
%
The first generalized coordinate imaging step is performing wavefield extrapolation
\begin{eqnarray}
E_{\xi_3}[S(\xi_3-\Delta \xi_3,\xi_1|\omega)] &=& S(\xi_3,\xi_1|\omega) \nonumber \\
E^*_{\xi_3}[R(\xi_3-\Delta \xi_3,\xi_1|\omega)] &=& R(\xi_3,\xi_1|\omega).
\end{eqnarray}
Generalized coordinate images then are constructed by evaluating an imaging condition
\begin{equation} \label{eqn:subsurface}
I(\xi_3,\xi_1,h_{\xi_1}) = \sum_\omega R^*(\xi_3,\xi_1 + h_{\xi_1}|\omega) S(\xi_3,\xi_1 - h_{\xi_1}|\omega). 
\end{equation}
where $h_{\xi_1}$ is the $\boldsymbol{\xi}$-coordinate equivalent of Cartesian subsurface offset axis $h_{\xi_1}$.   The generalized coordinate ADCIG volume is generated by applying an offset-to-angle transformation ${\bf T}_{h_{\xi}\rightarrow\gamma}$
\begin{equation} \label{eqn:subsurface5}
I(\xi_3,\xi_1,\gamma) = {\bf T}_{h_{\xi_1}\rightarrow \gamma} I(\xi_3,\xi_1,h_{\xi_1}).
\end{equation}
Conventional ADCIG volumes can be recovered by sinc interpolation of the $\boldsymbol{\xi}$-coordinate image in equation~\ref{eqn:subsurface5} to the final Cartesian coordinate  volume.  
%
\par
%
Figure~\ref{fig:RC2} illustrates this process  using the elliptic coordinate system.  Panel~\ref{fig:RC2}a shows the BP synthetic velocity model \cite[]{bpdata} with an elliptic mesh overlain.  Note that the salt flanks to the right-side of the model are nearly vertical in Cartesian coordinates.   Panel~\ref{fig:RC2}b shows the velocity model in panel \ref{fig:RC2}a interpolated to the elliptic coordinate system.   Importantly, the salt flanks in the elliptic coordinate system are nearly horizontal, which should lead to ADCIG calculations more robust than in Cartesian coordinates.   However, proving this requires understanding the differences, if any, between operators ${\bf T}_{h_{x_1}\rightarrow \gamma}$ and ${\bf T}_{h_{\xi_1}\rightarrow \gamma}$ in equations~\ref{eqn:subsurface2}~and~\ref{eqn:subsurface5}.
%
\plot{RC2}{width=\textwidth}{Prestack migration test in elliptic  coordinates.  a) Benchmark synthetic velocity model with an overlying elliptic coordinate system. b) Effective slowness model in the transformed elliptic coordinate system in a).  c) Benchmark synthetic velocity model with a different overlying elliptic coordinate system.  d) Effective elliptic coordinate slowness model for the coordinate system in c). } 
%
\subsection{Cartesian coordinate ADCIGs}
%
For constant velocity media in Cartesian coordinates, a straightforward link exists between differential changes in the travel time, $t$, of rays connecting the source-reflector and reflector-receiver paths to changes in the subsurface offset, $h_{x_1}$, and depth, $x_3$, coordinates.  Figure~\ref{fig:rays}a shows the geometry of these variables.    
%
\plot{rays}{width=\textwidth}{Cartoon illustrating the geometry of the ADCIG calculation.  a) Cartesian coordinates.  b) Generalized coordinates.  Adapted from \cite{SavaADCIG}. }   
%
Mathematically, these relationships are
\begin{eqnarray}\label{eqn:dtdh}
\left. \left[ \begin{array}{c}
\dtdhxx \\ \dtdxz
\end{array}  \right] \right|_{t,x_1}
 = 2 \,s \, \rcos \, \alpha
\left[ \begin{array}{c}
\rsin \, \gamma \\ \rcos \, \gamma
\end{array}  \right],
\end{eqnarray}
where $s$ is slowness, $\alpha$ is reflector dip, and $\gamma$ is the reflection opening angle.  The right-hand-side of equations~\ref{eqn:dtdh} are derived by \cite{SavaADCIG}.  Equations~\ref{eqn:dtdh} can be rewritten as 
\begin{equation} \label{eqn:Angle}
 \dxzdhxx  = - \left. \dtdhxx \right/ \dtdxz = \rtan \, \gamma, 
\end{equation} 
where the negative signs arise through use of the implicit functions theory \cite[]{ImplicitFunction}.    Note that the insensitivity of Cartesian ADCIGs is caused situations where $\dtdxz \rightarrow 0$, - i.e. for steeply dipping structures.
%
\par
%
Finally, because equation~\ref{eqn:Angle} has no explicit $\boldsymbol{x}$ dependence, Fourier-based methods can be used to calculate the reflection opening angle $\gamma$ directly in the frequency domain
\begin{equation} \label{eqn:Fangle} 
\rtan \gamma = - \frac{k_{h_{x_1}}}{k_{x_3}},
\end{equation}
where $k_{h_{x_1}}$ and $k_{x_3}$ are the wavenumbers in the $h_{x_1}$ and ${x_3}$ directions, respectively. 
%
\subsection{Generalized coordinate ADCIGs} 
%
Figure~\ref{fig:rays}b illustrates a scenario similar to that in Figure~\ref{fig:rays}a, but for a more general coordinate system.  The reflection opening angle,  $\gamma$, and the reflector dip, $\alpha$, obviously remain unchanged in the subsurface; however, the orientations of the $h_{\xi_1}$ and $\xi_3$ axes used to estimate $\gamma$ now differ. The key question is which quantities in the ADCIG calculation are affected by this change of variables?  
%
\par
%
To answer these questions, I first assume that generalized coordinate systems are related to the Cartesian variables through a bijection (i.e., one-to-one mapping) 
\begin{equation} \label{eqn:mapping}
x_1 = f(\xi_1,\xi_3) \quad {\rm and} \quad x_3 = g(\xi_1,\xi_3)
\end{equation}
with a non-vanishing Jacobian of coordinate transformation, $J_{\boldsymbol{\xi}}$.  The bijection between a generalized and Cartesian coordinate system allows us to rewrite equations \ref{eqn:dtdh} as  \cite[]{Widder}
\begin{equation} \label{eqn:dtdh2}
\frac{\partial t}{\partial {x_1}} = \frac{1}{J_{\boldsymbol{\xi}}}\frac{\partial (t,x_2) }{ \partial
 ({\xi_1},{\xi_3}) }
\quad {\rm and} \quad 
\frac{\partial t}{\partial x_3} = \frac{1}{J_{\boldsymbol{\xi}}}\frac{\partial (x_1,t) }{ \partial
 (\xi_1,\xi_3) }.
\end{equation} 
Expanding out the Jacobian notation leads to  
\begin{eqnarray}\label{eqn:PD}
%\left. 
\left[ \begin{array}{c}
\dtdex \dxzdez - \dtdez\dxzdex \\
\dtdez \dxxdex - \dtdex\dxxdez
\end{array} \right] 
%\right|_{\xi_1,t} 
= 
2\, J_{\boldsymbol{\xi}} \,s\, \rcos\, \gamma
\left[ \begin{array}{c}
\rsin \, \alpha \\
\rcos \, \alpha 
\end{array} \right].
\end{eqnarray}
The right-hand-sides of equations~\ref{eqn:PD} are analogous to those derived by \cite{SavaADCIG}.  Cross-multiplying the expressions by factors $\frac{\partial
  x_1}{\partial \xi_3}$ and $\frac{\partial x_3}{\partial \xi_3}$ 
\begin{eqnarray}\label{eqn:PD2}
%\left. 
\left[ \begin{array}{c}
\frac{\partial
  x_1}{\partial \xi_3}
\left(\dtdex \dxzdez - \dtdez\dxzdex\right) \\
\frac{\partial x_3}{\partial \xi_3}
\left(\dtdez \dxxdex - \dtdex\dxxdez\right)
\end{array} \right] 
%\right|_{\xi_1,t} 
= 
2\, J_{\boldsymbol{\xi}} \,s\, \rcos\, \gamma
\left[ \begin{array}{c}
\rsin \,\alpha \\
\rcos \,\alpha 
\end{array} \right]
\end{eqnarray}
and adding the two expressions results in 
\begin{equation} \label{eqn:PD3}
\dtdez \left( \dxzdez \dxxdex - \dxxdez \dxzdex \right) = 
2\, J_{\boldsymbol{\xi}} \,s\, \rcos\, \gamma \left( \dxxdez \rsin  \,\alpha + \dxzdez \rcos \,
\alpha \right).
\end{equation}
%
\par
%%%%%%%%%%%%%%%%%%%%%%%%
%
A similar argument can be used to construct the equations for the subsurface offset axis.  The bijection between the generalized coordinate and Cartesian subsurface offsets allows for equations~\ref{eqn:dtdh} to be rewritten as 
\begin{equation}
\frac{\partial t}{\partial h_{x_1} } = \frac{1}{J_{\boldsymbol{h}}}\frac{\partial
 (t,h_{x_2}) }{ \partial 
 (h_{\xi_1},h_{\xi_3}) }
\quad {\rm and} \quad 
\frac{\partial t}{\partial h_{x_3}} = \frac{1}{J_{\boldsymbol{h}}}\frac{\partial
 (h_{x_1},t) }{ \partial 
 (h_{\xi_1},h_{\xi_3}) },
\end{equation}
where $J_{\boldsymbol{h}}$ is the subsurface-offset Jacobian of transformation. Expanding out the Jacobian notation leads to  
\begin{eqnarray}\label{eqn:PDH}
%\left. 
\left[ \begin{array}{c}
\dtdhex \dhxzdhez - \dtdhez\dhxzdhex \\
\dtdhez \dhxxdhex - \dtdhex\dhxxdhez
\end{array} \right] 
%\right|_{\xi_1,t}
= 
2\, J_{\boldsymbol{h}} \,s\, \rsin\, \gamma
\left[ \begin{array}{c}
\rcos\, \alpha \\
\rsin \, \alpha 
\end{array} \right].
\end{eqnarray}
The right-hand-side of equations~\ref{eqn:PDH} are again analogous to those given by \cite{SavaADCIG}.  Cross-multiplying the expressions by factors $\frac{\partial h_{x_1}}{\partial h_{\xi_1}}$ and $\frac{\partial h_{x_3}}{\partial h_{\xi_1}}$  
\begin{eqnarray}\label{eqn:PDH2}
%\left. 
\left[ \begin{array}{c}
\frac{\partial h_{x_1}}{\partial h_{\xi_3}}
\left(\dtdhex \dhxzdhez - \dtdhez\dhxzdhex\right) \\
\frac{\partial h_{x_3}}{\partial h_{\xi_3}}
\left(\dtdhez \dhxxdhex - \dtdhex\dhxxdhez\right)
\end{array} \right] 
%\right|_{\xi_1,t} 
= 
2\, J_{\boldsymbol{h}} \,s\, \rsin\, \gamma
\left[ \begin{array}{c}
\rcos \, \alpha \\
\rsin \, \alpha 
\end{array} \right],
\end{eqnarray}
and adding the two expressions above yields
\begin{equation} \label{eqn:PDH3}
\dtdhex \left( \dhxxdhex \dhxzdhez - \dhxxdhez \dhxzdhex \right) = 
2\, J_{\boldsymbol{h}} \,s\, \rsin\, \gamma \left( \dhxxdhex \rcos \alpha - \dhxzdhex \rsin
\alpha \right).
\end{equation}
An expression for ADCIGs can be obtained by dividing equation~\ref{eqn:PDH3} by equation~\ref{eqn:PD3}
\begin{equation} \label{eqn:final1}
\frac{\dtdhex}{\dtdez} \frac
{\left( \dhxxdhex \dhxzdhez - \dhxxdhez  \dhxzdhex \right)}
{\left( \dxzdez   \dxxdex   - \dxxdez    \dxzdex   \right)}
 =  \rtan
\, \gamma \, \frac{J_{\boldsymbol{h}}}{J_{\boldsymbol{\xi}}} \, \frac
{\left( \dhxxdhex \rcos \,\alpha - \dhxzdhex \rsin \,\alpha \right)}
{\left( \dxzdez   \rcos \,\alpha + \dxxdez   \rsin \,\alpha \right)}.
\end{equation}
%
\par
%
One question arising from the geometric factors in equation~\ref{eqn:final1} is what do the terms $\dhxxdhex$, $\dhxzdhex$, $\dhxxdhez$ and $\dhxzdhez$ represent? In the examples below, I assume that the subsurface offset axes are geometrically equivalent to the coordinate axes
\begin{eqnarray}
\left[\begin{array}{c}
h_{x_1}  \\ 
h_{x_3}  \\
h_{\xi_1}\\
h_{\xi_1}
\end{array} 
\right] = \left[ \begin{array}{c}
x_1   \\ 
x_3   \\
\xi_1 \\
\xi_3
\end{array}\right] 
\quad {\rm such} \,\,\, {\rm that} \quad
\left[\begin{array}{c}
\dhxxdhex \\
\dhxzdhex \\
\dhxxdhez \\
\dhxzdhez
\end{array} 
\right] = \left[ \begin{array}{c}
\dxxdex \\
\dxzdex \\
\dxxdez \\
\dxzdez
\end{array}\right].
\end{eqnarray} 
This assumption is realized by uniform wavefield shifting.  Using these identities in equation~\ref{eqn:PD3} reduces equation~\ref{eqn:final1} to 
\begin{equation} \label{eqn:final3}
\dezdhex = - \left. \dtdhex \right/\dtdez  =  - \rtan\, \gamma \, \frac
{\left(\dxxdex \rcos\,\alpha - \dxzdex \rsin \,\alpha \right)} 
{\left( \dxzdez \rcos\,\alpha + \dxxdez \rsin \,\alpha \right)},
\end{equation}
where the two Jacobian transformations are equivalent (i.e. $J_{\boldsymbol{\xi}}=J_{\boldsymbol{h}}$). 
One important observation from equation~\ref{eqn:final3} is that if  the coordinate system obeys the Cauchy-Riemann conditions of differentiability \cite[]{ComplexAnalysis}
\begin{equation} \label{eqn:CRcriteria}
\dxxdex = \dxzdez \quad \hbox{and} \quad \dxzdex = - \dxxdez,
\end{equation}
then the following holds
\begin{equation} \label{eqn:final2}
\dezdhex =  - \rtan \, \gamma.
\end{equation}
This is the generalized coordinate equivalent of the Cartesian expression in equation~\ref{eqn:dtdh}.  Similar to Cartesian coordinates, elliptic coordinate ADCIGs become insensitive where structural dips cause $\dtdez \rightarrow 0$.  However, this insensitivity can be minimized when using many generalized coordinate systems, because structural dips appear at different angles in different translated elliptic meshes.   Figures~\ref{fig:RC2}c-d illustrate this showing a different coordinate shift than in panels~\ref{fig:RC2}a-b.  Note the changes in structural dip in the right-hand-side of the elliptic coordinate panels. Thus, while ADCIGs calculated on one elliptic mesh can be insensitive to certain structure, coordinate translations will ensure that ADCIGs become sensitive for other elliptic coordinate systems.  
%
\par
%
Finally, one may calculate reflection opening angles in the Fourier domain for coordinate systems satisfying equations~\ref{eqn:CRcriteria} 
\begin{equation} \label{eqn:STEP3}
\rtan \, \gamma = - \frac{k_{h_{\xi_1}}}{k_{\xi_3}},
\end{equation}
where $k_{h_{\xi_1}}$ and $k_{\xi_3}$ are the wavenumbers in the $h_{\xi_1}$
and ${\xi_3}$ directions, respectively.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Canonical Examples} 
%
This section presents three canonical examples that illustrate the generalized ADCIG theory: tilted Cartesian, polar, and elliptic  coordinate systems.  Figure~\ref{fig:COORDS} presents schematic examples of these three coordinate systems.
%
\plot{COORDS}{width=\textwidth}{Canonical coordinate system examples. a) Tilted Cartesian coordinates. b) Polar coordinates. c) Elliptic coordinates. }  
%
\subsection{Tilted Cartesian coordinates} 
%
Tilted Cartesian coordinates are a useful generalized coordinate system for performing migration (see Figure~\ref{fig:COORDS}a). \cite{SEG2004sgj} use this mesh in a plane-wave migration scheme where the coordinate system is oriented toward the plane-wave take-off angle.  A tilted Cartesian mesh is defined by  
\begin{eqnarray}
\left[ \begin{array}{c}
x_1\\
x_3 
\end{array}  \right] = 
\left[ \begin{array}{cc}
\rcos \, \theta &-\rsin \,\theta \\
\rsin \, \theta & \rcos \,\theta
\end{array}\right]
\left[ \begin{array}{c}
\xi_1 \\
\xi_3
\end{array} \right],
\end{eqnarray}
where $\theta$ is the tilt angle.  The partial derivative transform matrix is  
\begin{eqnarray}
\left[ \begin{array}{cc}
\dxxdex & \dxxdez \\
\dxzdex & \dxzdez
\end{array} \right] = 
\left[ \begin{array}{cc}
\rcos \, \theta &-\rsin \, \theta \\
\rsin \, \theta & \rcos \, \theta
\end{array}\right],
\end{eqnarray}
which leads to the following ADCIG equation:
\begin{equation}
\dezdhex  =  - \rtan \, \gamma \, \frac{\left( \rcos \,\theta\, \rcos \, \alpha  +\rsin \,\theta\, \rsin \, \alpha \right) }  {\left(   \rcos \,\theta\, \rcos \, \alpha  +\rsin \,\theta\, \rsin \, \alpha \right) } = - \rtan \, \gamma.
\end{equation}
Thus, calculating ADCIGs in tilted Cartesian coordinates will directly recover the correct reflection opening angle.   Note that setting tilt angle $\theta=0^\circ$ recovers the Cartesian expression in equation~\ref{eqn:Angle}.
%
\subsection{Polar coordinates}
%
The polar coordinate system (see Figure~\ref{fig:COORDS}b), where the extrapolation direction is oriented along the radial direction, is appropriate for generating 2D Green's function estimates.  The polar coordinate system is defined by 
\begin{eqnarray}
\left[ \begin{array}{c}
x_1\\
x_3 
\end{array}  \right] = 
\left[ \begin{array}{c}
a\, \xi_3 \,\rcos \,\xi_1 \\
a\,\xi_3\, \rsin \,\xi_1 \\
\end{array} \right].
\end{eqnarray} 
The partial derivative transformation matrix is 
\begin{eqnarray} \label{eqn:testy}
\left[ \begin{array}{cc}
\dxxdex & \dxxdez \\
\dxzdex & \dxzdez
\end{array} \right] = 
\left[ \begin{array}{cc}
-a \,\xi_3\, \rsin\, \xi_1  & a \, \rcos\, \xi_1 \\
 a \,\xi_3 \,\rcos \,\xi_1  & a \, \rsin\, \xi_1
\end{array}\right],
\end{eqnarray}
which leads to the following ADCIG equation:
\begin{equation}
\dezdhex  =  -\xi_3 \, \rtan \, \gamma \, \frac{ \left(  -\rsin\,\xi_1\,\rcos\,\alpha - \rcos\,\xi_1\,\rsin\,\alpha  \right) } {\left(  \rsin\,\xi_1\rcos\,\alpha+\rcos\,\xi_1 \,\rsin,\alpha  \right) }  = \xi_3 \,\rtan \, \gamma.
\end{equation}
Thus, one cannot calculate ADCIGs directly with Fourier-based methods in polar coordinates because of the spatial geometric dependence on $\xi_3$.  
%
\subsection{Elliptic coordinates} 
%
Elliptic coordinates (see Figure~\ref{fig:COORDS}c) are a useful coordinate system for performing 2D shot-profile migration \cite[]{ShraggeGEONO}.  An elliptic mesh is defined by
\begin{eqnarray}
\left[ \begin{array}{c}
x_1\\
x_3 
\end{array}  \right] = 
\left[ \begin{array}{c}
a\, {\rm cosh}\, \xi_3 \, \rcos \,\xi_1\\
a\, {\rm sinh}\, \xi_3 \, \rsin \,\xi_1
\end{array} \right].
\end{eqnarray}
The partial derivative transformation matrix is
\begin{eqnarray}
\left[ \begin{array}{cc}
\dxxdex & \dxxdez \\
\dxzdex & \dxzdez
\end{array} \right] = 
a \left[ \begin{array}{cc}
 {\rm cosh} \,\xi_3 \, \rsin  \,\xi_1 & {\rm sinh} \,\xi_3 \, \rcos \,
 \xi_1\\  
-{\rm sinh} \,\xi_3 \, \rcos \,\xi_1 & {\rm cosh} \,\xi_3 \, \rsin
 \,\xi_1 
\end{array}\right],
\end{eqnarray}
which leads to the following ADCIG equation:
\begin{equation} \label{eqn:ELLIP}
\dezdhex  = - \rtan
\, \gamma \, \frac{
\left( 
 {\rm cosh} \,\xi_3 \, \rsin \,\xi_1 \, \rcos \,\alpha 
- {\rm sinh} \,\xi_3 \, \rcos \,\xi_1 \, \rsin \,\alpha
\right) } {\left( 
  {\rm cosh} \,\xi_3 \, \rsin \,\xi_1 \, \rcos \,\alpha 
- {\rm sinh} \,\xi_3 \, \rcos \,\xi_1 \, \rsin \,\alpha 
 \right) } = -\rtan \, \gamma.
\end{equation}

Thus, calculating ADCIGs in elliptic coordinates with Fourier-based methods will directly recover the true reflection opening angle.  

\section{Numerical Examples} 
%
This section presents numerical tests of the generalized theory by comparing the ADCIG volumes for elliptic and Cartesian coordinate systems.  The results were generated using a shot-profile migration algorithm altered to account for elliptic geometry.  Replacing the Cartesian extrapolator, $E_{x_3}[\cdot]$  with an elliptic coordinate operator, $E_{\xi_3}[\cdot]$, represented the only significant code alteration.  In both coordinate systems, the implemented extrapolation operators were accurate to roughly $80^\circ$ \cite[]{GEO50-10-16341637} relative to the extrapolation axis.  The reader is referred to \cite{ShraggeGEORWE} for further implementation details on shot-profile migration in elliptic coordinates.
%
\par
%
I calculated ADCIG image volumes for each shot-profile following a three-step procedure: 1) generate image volume $I(\xi_1,\xi_3,h_{\xi_1})$ by computing 64 subsurface shifts at each point in every extrapolation step; 2) calculate ADCIG volume $I(\xi_1,\xi_3,\gamma)$ using the procedure described in \cite{SavaADCIG}; and 3) output the image $I(x_1,x_3,\gamma)$ by sinc-based interpolation of single-shot ADCIGs to the global volume.   Steps 1-3 were repeated for all shot-profile sections contributing to the final image.
%
\subsection{Test 1: Elliptic Isochrons}
%
The first test involves imaging reflectors with an elliptic geometry from two offset point sources in a constant velocity medium.  Figure~\ref{fig:Theory} shows the experimental geometry.  The source wavefield impulse is at time $t=0$~s and offset $h=-1.12$~km, while the receiver wavefield impulses are at $t$ = 2.0, 2.5, 3.0, and 3.5~s and $h=1.12$~km.  Correlating these two wavefields leads to four elliptic isochrons whose positions are given by
\begin{eqnarray}
x_1 (t,h) & = & \frac{vt}{4} \left( \rcos \, \theta_1 + \rcos \, \theta_2 \right) \\
x_3 (t,h) & = & \frac{vt}{4} \left(1 - \frac{h^2}{v^2t^2} \right) 
\left( \rsin \, \theta_1 + \rsin \, \theta_2 \right),
\end{eqnarray}
where, given a shooting angle $\theta_1$, the second angle $\theta_2$ is obtained by
\begin{equation}
\theta_2 = {\rm cos}^{-1} \left( \frac{vt}{2h}  \left[ \frac{1}{\frac{2}{1-\frac{4h^2}{v^2t^2}} - \frac{1}{1-\frac{2h \rcos\,\theta_1}{vt}} } -1 \right]  \right). 
\end{equation}
The reflection opening angle, given by $\gamma = 90 - \frac{(\theta_1 +\theta_2)}{2}$, is shown color-coded on the scatterplot in Figure~\ref{fig:Theory}. 
%
\plot{Theory}{width=\textwidth}{Theoretical results for an elliptic isochron for four travel times in a constant velocity medium.  The elliptic surfaces are color-coded according to reflection opening angle.}
%
\par
%
Figure~\ref{fig:EllipticTest} shows the ADCIG volumes calculated in elliptic and Cartesian coordinates.   Panel \ref{fig:EllipticTest}a presents the elliptic coordinate (EC) image extracted at -18.5$^\circ$ from the ADCIG volume.   The analytical ADCIG locations are represented by black bullets.  The breadth of the analytic ADCIG response is due to coarse binning in the angular domain.  Note that the analytically and numerically generated results are well matched.  Figure~\ref{fig:EllipticTest}a also shows three vertical lines indicating the locations of the three ADCIGs in panels \ref{fig:EllipticTest}b-d. Again, the analytic and numerical ADCIGs are well matched, though less so at shallower  depths.
%
\par
%
The Cartesian coordinate (CC) ADCIGs are presented in panels~\ref{fig:EllipticTest}e-h.  Panel \ref{fig:EllipticTest}e shows the Cartesian image again extracted at -18.5$^\circ$ in the angle domain.  Panels \ref{fig:EllipticTest}f-h present three ADCIGs at the same locations as in panels \ref{fig:EllipticTest}b-d.  Note that the Cartesian image volumes are well-matched to the elliptic coordinate examples.  Overall, good agreement between the theoretical results and the wavefield volume is observed in both images and energy is focused in the neighborhood of the correct locations. The angle gathers, though, are not always centered relative to the true location, which is again more noticeable at shallower depths where the numerical ADCIG volumes tend to overestimate the reflection opening angle.  
%
\plot{EllipticTest}{width=\textwidth,height=4in}{Elliptic reflector comparison tests between analytically (black bullets) and numerically generated ADCIG volumes.  Panels a-d are in elliptic coordinates (EC), while panels e-h are in Cartesian coordinates (CC).  a) EC image extracted at the -18.5$^\circ$ reflection angle.  b) EC angle gather at -1.5~km. c) EC gather at 0.75~km. d) EC gather at 0.0~km.  e) CC image extracted at the -18.5$^\circ$ reflection angle.  f) CC angle gather at -1.5~km. g) CC gather at 0.75~km. h) CC gather at 0.0~km.} 
%
\par
%
A second observation from Figure~\ref{fig:EllipticTest} is that while elliptic coordinate ADCIGs focus at the correct reflection opening angle, they tend to have greater residual curvature away from the correct image position than in the Cartesian coordinate images.  Figure~\ref{fig:AngleExplanation} illustrates why this occurs using near- and far-offset shift examples for Cartesian, tilted Cartesian, and elliptic coordinate systems.  Each panel has two sets of  perpendicular axes, the intersection of each indicating the two points being cross-correlated.  The red lines show the length of the $\xi_1$ coordinate axis, while the blue lines indicate the Euclidean distance between the two correlated points.   
%
\sideplot{AngleExplanation}{width=\textwidth}{Sketches illustrating wavefield shift differences between elliptic and (tilted) Cartesian coordinates for near- and far-offset shifts.  The intersection points of the black and gray perpendiculars indicate the points being correlated.  Red lines indicate the length along the coordinate system axes, while blue lines show the Euclidean distance between the two points.  Red and blue lines are equal for Cartesian and tilted Cartesian coordinate systems for both near- and far-offset shifts; however, they are unequal at far offsets for the elliptic coordinates generating residual curvature.}  
%
\par
%
In the Cartesian and tilted Cartesian panels, the red and blue line lengths are the same for both near- and far-offsets.  However, the blue and red lines in the elliptic panel are not the same length and become increasingly disparate at greater offsets.  The distance along the $\xi_1$ axis is thus greater than the Euclidean distance between the points, leading to a greater apparent  velocity and downward residual ADCIG curvature.  However, as the next tests demonstrate, this residual curvature is largely stacked out in the final volumes and does not contaminate the ADCIG results.
%
\subsection{Test 2: BP velocity model}
%
The second test compares elliptic and Cartesian coordinate ADCIG volumes computed for the BP synthetic velocity model.  Figure~\ref{fig:GOOD} shows slices from the corresponding elliptic and Cartesian ADCIG image volumes.  
%
\plot{GOOD}{width=\textwidth,height=3.5in}{Vertical elliptic and Cartesian ADCIGs slices using the correct migration velocity model.  a) Elliptic coordinate image with three vertical lines showing the locations of ADCIG gathers in panels b-d.  e) Cartesian coordinate image with three vertical lines showing the locations of ADCIG gathers in panels f-h.}
%
Panel~\ref{fig:GOOD}a shows an elliptic coordinate image with three vertical lines that indicate the angle gather locations in Figures~\ref{fig:GOOD}b-d. Panel~\ref{fig:GOOD}b shows mostly one-sided reflectivity, save for the second salt reflection in depth exhibiting energy to both sides.  I attribute this to the reversal of source and receiver wavefield orientations within the salt.  The image in panel~\ref{fig:GOOD}d has a wide reflection zone between 3.75-4.25~km in depth, which occurs because the shown angle-gather is a vertical slice through the vertical salt-flank reflector.  Panel~\ref{fig:GOOD}e shows the Cartesian image for the same location as panel~\ref{fig:GOOD}a, while panels~\ref{fig:GOOD}f-h are extracted from the same locations as panels~\ref{fig:GOOD}b-d.  The Cartesian gathers look similar to the corresponding elliptic coordinate gathers.  The angle gathers from steep reflectors are, except for the salt flanks to the right-hand-side of panel~\ref{fig:GOOD}h.
%
\par
%
Figure~\ref{fig:GOODZ} shows the image volume using horizontal slices that better resolve the vertical salt flank.  Panel~\ref{fig:GOODZ}a shows the elliptic coordinate image with three horizontal lines showing the ADCIG slice locations.  The right-hand sides of panels~\ref{fig:GOODZ}b-d show the vertical salt-flank reflector that is well focused.   This demonstrates the robustness of the ADCIG calculation in elliptic coordinates.
%
\plot{GOODZ}{width=\textwidth,height=7.5in}{Horizontal elliptic and Cartesian ADCIGs slices using the correct  migration velocity model.  a) Elliptic coordinate image with three horizontal lines showing the locations of horizontal ADCIG gathers in panels b-d.  e) Cartesian coordinate image with three horizontal lines showing the locations of horizontal ADCIG gathers in panels f-h.}
%
Panel~\ref{fig:GOODZ}e shows the Cartesian coordinate image with three horizontal lines showing the locations of the ADCIG slices.  The right-hand salt-flank reflector in panel~\ref{fig:GOODZ}f is similarly well-resolved, largely because the structural dip is relatively low.  However, the salt-flank reflector images in panels~\ref{fig:GOODZ}g-h are somewhat blurred out.  I attribute this to the combined effects of inaccurate large-angle extrapolation and insensitivity of the ADCIG calculation to steep structural dip.
%
\par
%
An additional test is examining how the ADCIG volumes change when introducing an incorrect migration velocity profile.  Figure~\ref{fig:BAD} presents ADCIG volumes similar to those shown in Figures~\ref{fig:GOOD} after using a migration velocity profile rescaled by 95\%.  
%
\plot{BAD}{width=\textwidth,height=3.5in}{Vertical elliptic and Cartesian ADCIGs slices using an incorrect  migration velocity model.  a) Elliptic coordinate image with three vertical lines showing the locations of vertical ADCIG gathers in panels b-d.  e) Cartesian coordinate image with three vertical lines showing the locations of vertical ADCIG gathers in panels f-h.}
%  
Note that both images are not well focused and have residual curvature indicating an incorrect migration velocity.  Because the reflectors are near vertical, though, the residual moveout is relatively small.  The situation becomes much clearer when examining the horizontal slices shown in Figure~\ref{fig:BADZ}.
%
\plot{BADZ}{width=\textwidth,height=7.5in}{Horizontal elliptic and Cartesian ADCIGs slices using an incorrect migration velocity model.  a) Elliptic coordinate image with three horizontal lines showing the locations of horizontal ADCIG gathers in panels b-d.  e) Cartesian coordinate image with three horizontal lines showing the locations of horizontal ADCIG gathers in panels f-h.}
%
Note that the angle gathers in panels~\ref{fig:BADZ}b-d, and especially to the right-hand side in panel~\ref{fig:BADZ}d, show much greater residual curvature.  Importantly, the imaging enhancements afforded by elliptic coordinates  will improve any migration velocity analysis approach based on residual curvature in ADCIGs.
%%%%%%%%%%%%%%%%%%%%%%%
%
%
\section{Conclusions}
%
I extend the Cartesian ADCIG theory to 2D generalized coordinate systems.  The expressions for ADCIGs contain additional coordinate factors describing mesh geometry and wavefield  shifts.  I show that these geometric expressions cancel out for coordinate systems satisfying the Cauchy-Riemann criteria of differentiability,  which include tilted Cartesian and elliptic meshes.  The method for calculating ADCIGs in elliptic coordinates is thus very close to the Cartesian coordinate procedure.  I confirm the validity of the approach by comparisons of analytically and numerically generated ADCIG domains, and with tests on the BP synthetic data set.  Calculated ADCIGs volumes are more robust when computed in elliptic coordinates than in Cartesian coordinate due to improved large-angle propagation and improved sensitivity to steep structural dips.  Finally, the imaging enhancements afforded by elliptic coordinates  will improve any migration velocity analysis approach based on residual curvature in ADCIGs.
%
\section{Acknowledgments}
%
I would like to thank Biondo Biondi, Alejandro Valenciano, Antoine Guitton and Brad Artman for helpful discussions, and BP for releasing the synthetic data used in the examples.
%
%
%\bibliographystyle{seg}
%\bibliography{SEP,SEG,MISC,GEOTLE,allrefs}
 
\cleardoublepage
