\chapter{Delayed-shot migration in TEC coordinates} %\author{Jeff Shragge}
%\date{} % delete this line to display the current date
%\lefthead{Shragge and Shan}
%\righthead{TEC Coordinates}
%\renewcommand{\thefootnote}{\fnsymbol{footnote}}
%\newcommand{\boldvec}[1]{\boldsymbol{\mathrm{#1}}}
%\let\vec\boldvec
\def\w{\omega}
\def\bX{\boldsymbol{x}}
\def\bE{\boldsymbol{\xi}}

%\maketitle

%\author{Jeff Shragge and Guojian Shan}

%\begin{abstract}
This chapter extends the analytical RWE approach discussed in Chapter~3 to 3D coordinate systems.  I formulate an inline delayed-shot migration procedure in tilted elliptical-cylindrical (TEC) coordinate systems.  When inline coordinate tilt angles are well-matched to the inline source ray parameters, the TEC coordinate extension affords accurate propagation of both steep-dip and turning-wave components.  I show that wavefield extrapolation in TEC coordinates is no more complicated than propagation in elliptically anisotropic media.  Impulse response tests illustrate the accuracy and lack of numerical anisotropy of the implemented scheme.  I apply this approach to a realistic 3D wide-azimuth synthetic derived from a field Gulf of Mexico data set.  The resulting images demonstrate the imaging advantages made possible through 3D RWE implementations, including the improved imaging of steeply dipping salt flanks, potentially at a reduced computational cost.  Narrow-azimuth migration results demonstrate the applicability of the approach to typical Gulf of Mexico field data.
%\end{abstract}
%
\section{Introduction} 
%
Wave-equation migration (WEM) methods routinely generate accurate seismic images in areas of complex geology.  One common class of WEM approaches is shot-profile migration using one-way wavefield extrapolation.  The first shot-profile migration step is to specify source and receiver wavefields that consist of modeled point sources and an individual shot profile, respectively.  The migration algorithm propagates these two wavefields through the velocity model and correlates them at each extrapolation step to form an image.  Although this procedure generates high-quality migration results, two drawbacks make shot-profile migration a less-than-ideal strategy.  The first issue is that each individual shot migration requires a large aperture to propagate energy to wide offsets.  The second drawback is that one migrates each shot record individually, which can be computationally expensive for large surveys with a high shot density.
%one drawback is that individual profiles routinely cover only a small portion of the total survey %area and many independent shot-profile migrations are required to infill the full image volume.  
%
\par
%
One way to make the shot-profile style of WEM more efficient is to migrate a reduced number of composite source and receiver profiles each covering a broader aperture.  For example, one can image a number of shot profiles simultaneously on the same migration domain.  The key idea is that one makes a computationally advantageous trade-off of a broader migration aperture for a reduced number of shots.   Shot-profile migration with composite wavefields, though, leads to the mixing of information from different shots and generates image crosstalk.  A number of authors address this problem using a variety of phase-encoding migration approaches \cite[]{Morton1998, Jing2000, Romero2000, Sun2002}, that minimize the deleterious crosstalk effects. 
%
\par
%
Plane-wave migration (PWM) is one technique for reducing total migration cost using composite wavefields \cite[]{whitmorephd, Mosher1998, Duquet2001, Zhang2003, Liu2004, LiuPW}.  As originally demonstrated by \cite{whitmorephd}, the key idea is to synthesize from the full wavefield volume the set of composite receiver wavefields that would have been recorded were a planar source function used.  Generally, the number of synthesized wavefields is fewer than the corresponding number of shot profiles.  One generates PWM images by propagating the modeled planar source and composite receiver wavefields through the velocity and computing a (weighted) correlation.  \cite{LiuPW} and \cite{Duquet2006} demonstrate that PWM is equivalent to shot-profile migration in the limit where one uses many plane waves with well-sampled plane-wave dip spectra. \cite{LiuPW} also prove that 3D PWM is equivalent to conical-wave migration of individual sail lines synthesized as inline composite wavefields.  The approach is termed conical wave because the source wavefronts form conic sections (in constant media) for non-zero inline plane-wave ray parameters.    
%
\par
%
The migration of plane- and conical-wave data, though more efficient than shot-profile migration, is similarly restricted in accuracy by one-way wavefield extrapolation assumptions.  The most common limitation is a difficulty in propagating waves at large angles and turning waves by design, both of which are important for accurately imaging salt flanks in complex geologic areas.  \cite{SEG2004sgj} circumvent this problem by implementing 3D PWM in tilted Cartesian meshes.  This coordinate system effectively orients the wavefield extrapolation axis toward the plane-wave take-off vector, enabling more accurate bulk propagation of plane-wave energy.  One logistical complication of performing fully 3D PWM is that it requires propagating image-space-sized data volumes on a number of meshes tilting in both the inline and cross-line directions.  This leads to a number of computational issues associated with the significant memory footprint.  
%
\par
%
This chapter presents an alternate phase-encoding approach that similarly uses alternative coordinate systems. The key differences between the work herein and that of \cite{SEG2004sgj} are two-fold.  The first difference is that I phase encode only according to the inline source coordinate, leading to the inline delayed-shot migration algorithm.  This leads to a straightforward coarse-grain parallelization of the migration tasks across individual sail lines, where each migration has a significantly smaller aperture than the corresponding image-space-sized PWM volumes.  A second efficiency gain over PWM is a reduction in the total number of migrations, because the number of sail lines is quite often fewer than the required number of cross-line plane waves.  Thus, the inline-delayed shot approach has attractive computational advantages over the 3D PWM technique.
%
\par
%
The second difference is that I perform the migration in tilted elliptical-cylindrical (TEC) coordinates, rather than tilted Cartesian meshes.  TEC coordinate systems, formed by concatenating a set of the 2D elliptical coordinates discussed in Chapter~3 along the invariant third axis, are well-suited for migrating individual sail lines.  I orient the invariant and elliptical surfaces in the inline and cross-line directions, respectively.  I extrapolate the inline delay-shot synthesized wavefield volumes outward on a series of elliptical-cylindrical shells.  This allows source and receiver wavefields with zero inline dip to overturn in the cross-line direction, if necessary.  I introduce an extra degree of freedom that permits the coordinate system to tilt along the invariant inline axis, thus enabling the propagation of turning waves inline. Consequently, inline delayed-shot migration in TEC coordinates allows wavefields with most non-zero dips to propagate and overturn to all azimuths as appropriate.
%
\par
%
The chapter begins by examining 3D full-plane-wave and inline delayed-shot migration theory.  I then introduce the TEC coordinate geometry and develop the corresponding wavenumber that forms the basis of the TEC wavefield extrapolation operator.  I discuss the finite-difference extrapolation implementation and present the 3D impulse response.  I apply the technique to a 3D wide-azimuth synthetic data set derived from real Gulf of Mexico velocity model to demonstrate the imaging advantages of 3D RWE migration.  I then discuss the numerical costs associated with performing inline delayed-shot migration in TEC coordinates relative to Cartesian meshes.  The chapter concludes with narrow-azimuth migration results demonstrate the applicability of the approach to typical Gulf of Mexico field data.
%
\section{3D plane-wave migration}
%
The full plane-wave and inline delayed-shot migration theory discussed herein draws largely from \cite{LiuPW}.  I restate a number of key points for completeness, though with a slightly different notation.   As in previous chapters, I define Cartesian coordinates by $\boldsymbol{x}=[x_1,x_2,x_3]$ and a generalized coordinate system by $\boldsymbol{\xi}=[\xi_1,\xi_2,\xi_3]$.
%
\subsection{Full plane-wave phase-encoding migration}
%
Performing 3D plane-wave migration is similar in many respects to 3D shot-profile migration.  The main differences derive from how the composite source and receiver wavefield volumes, $\overline{S}$ and  $\overline{R}$, are re-synthesized from individual source and receiver profiles, $S_{jk}$ and $R_{lm}$, prior to imaging. The complete wavefields are generated by filtering the source and receiver profiles by a  function dependent on the inline and cross-line plane-wave ray parameters, $\boldsymbol{p_\xi}=[p_{\xi_1}, p_{\xi_2}]$.  These wavefields are then propagated through the migration domain  to generate the full source and receiver wavefield volumes 
\begin{eqnarray} \label{eqn:ImVolss}
\overline{S(\bE|\w)}= \sum_{j=1}^A \sum_{k=1}^B S_{jk}(\bE|\w) f(\w)
	 {\rm e}^{ i\w [p_{\xi_1} \Delta \xi_1 (j-p) + p_{\xi_2}
	    \Delta \xi_2 (k-q) ]}, \\ 
\overline{R(\bE|\w)}= \sum_{l=1}^A \sum_{m=1}^B R_{lm}(\bE|\w) f(\w)
	 {\rm e}^{ i\w [p_{\xi_1} \Delta \xi_1 (l-p) + p_{\xi_2}
	 	\Delta \xi_2 (m-q) ]},  
\end{eqnarray}
where $f(\w)$ is a frequency filter to be discussed below, $\Delta \xi_1$ and $\Delta \xi_2$ are the inline and cross-line sampling intervals, $p$ and $q$ are reference spatial indices in the inline and cross-line directions, $j$ and $k$ are indices fixing the inline and crossline source position, $l$ and $m$ are indices fixing the inline and cross-line receiver position, and $A$ and $B$ are the number of inline and cross-line source records, respectively.  The phase encoding, implemented at the surface independent of wavefield extrapolation, is valid for any generalized coordinate system.  Note that the wavefield propagation throughout the migration volume in equations~5.1~and~5.2 is understood, and assumed to be governed by the wavefield propagation techniques described in Chapter 2.
% 
\par
%
An image volume $I(\bE)$ is formed from a series of individual full plane-wave migration images, $I^{PW}(\bE| \boldsymbol{p_{\xi}})$, by correlating the composite plane-wave source and receiver wavefields and stacking the results over frequency.  The plane-wave migration kernel mixes source and receiver wavefield energy, $S_{jk}(\bE|\w)$ and $R_{lm}(\bE|\w)$, according to   
\begin{eqnarray} \label{eqn:PWimage1}
I(\bE) & = &\sum_{p_{\xi_1}} \sum_{p_{\xi_2}} \sum_{j,l=1}^A
\sum_{k,m=1}^B I^{PW}_{jklm}(\bE| \boldsymbol{p_\xi}) \, \\ \, & = &
\sum_{p_{\xi_1}} \sum_{p_{\xi_2}} \sum_{j,l=1}^A \sum_{k,m=1}^B
\sum_{\w} |f(\w)|^2 \, S_{jk}^*(\bE|\w) \,R_{lm}(\bE|\w)
{\rm e}^{ i \w[ p_{\xi_1} \Delta \xi_1 (j-l) + p_{\xi_2} \Delta
    \xi_2 (k-m) ]}, \nonumber 
\end{eqnarray}
where $^*$ indicates complex conjugate.
%
\par
%
Generally, mixing wavefields of differing $S_{jk}$ and $R_{lm}$ indices introduces image crosstalk.  A plane-wave migration image will be crosstalk-free, though, in the following limits: 
\begin{eqnarray}\label{eqn:Dapprox}
\lim_{N_{p_{\xi_1}} \rightarrow \infty} \sum_{\alpha =
  -N_{p_{\xi_1}}}^{N_{p_{\xi_1}}}{\rm e}^{ 
  i \w\alpha \Delta p_{\xi_1} \Delta {\xi_1} (j-l) } & = & |\w
|^{-1} \delta_{jl}, \nonumber \\
\lim_{N_{p_{\xi_2}} \rightarrow \infty} \sum_{\alpha =
  -N_{p_{\xi_2}}}^{N_{p_{\xi_2}}}{\rm e}^{ 
  i \w\alpha \Delta p_{\xi_2} \Delta {\xi_2} (k-m) } &=& |\w
|^{-1} \delta_{km}.
\end{eqnarray}
where $N_{p_{\xi_1}}$ and $N_{p_{\xi_2}}$ are the number of plane waves in the $\xi_1$ and $\xi_2$ directions.  Assuming that equation~\ref{eqn:Dapprox} approximately is valid (i.e., for large values of $N_{p_{\xi_1}}$ and $N_{p_{\xi_2}}$), I rewrite equation~\ref{eqn:PWimage1} as 
\begin{equation}
I(\bE) \approx \sum_{j=1}^A \sum_{k=1}^B \sum_{\w} |f(\w)|^2 |\w|^{-2}
S_{jk}^*(\bE|\w) R_{jk}(\bE|\w),
\end{equation}
which, by defining $|f(\w)|^2=|\w|^2$, generates the following expression:
\begin{equation}
I(\bE) \approx \sum_{j=1}^M \sum_{k=1}^N \sum_{\w} S_{jk}^*(\bE|\w)
R_{jk}(\bE|\w).
\end{equation}
This demonstrates the equivalence between plane-wave and shot-profile migration \cite[]{LiuPW}. 
% 
\subsection{Inline delayed-shot migration}
%
An alternate 3D migration formulation is to phase-encode individual sail lines for a given ray parameter, $p_{\xi_1}$, solely according to inline source position.  This phase-encoding approach is related to conical-wave migration, which requires $j-l=0$ in equation~\ref{eqn:PWimage1}.  However, I choose to not make this restriction because it is realized only by straight sail lines and non-flip-flop sources \cite[]{LiuPW}.  Rather, I present an alternative theory of inline delayed-shot migration that allows more general crossline source and receiver distribution.
%
\par
%
Inline delayed-shot wavefields, propagated through the migration domain to generate the full source and receiver wavefield volumes, are defined by
\begin{eqnarray} \label{eqn:ImVols}
\overline{S(\bE|\w)}= \sum_{l=1}^A \sum_{j=1}^B S_{jl}(\bE|\w) f(\w)
	 {\rm e}^{ i\w  [p_{\xi_1} \Delta \xi_1 (j-p)]}, \\ 
\overline{R(\bE|\w)}= \sum_{l=1}^A \sum_{k=1}^B R_{kl}(\bE|\w) f(\w)
	 {\rm e}^{ i \w [p_{\xi_1} \Delta \xi_1 (k-p)]},  
\end{eqnarray}
where $j$ and $k$ are the source and receiver inline position, respectively, $B$ is the number of inline records, $l$ is the sail line index out of a total of $A$ sail lines, and $p$ is a reference inline index.
%
\par
%
An image volume $I(\bE)$ is generated from a series of inline delayed-shot migration images, $I_l^{DS}(\bE| p_{\xi_1})$, formed by correlating the composite inline source and receiver wavefields and stacking the results over frequency.  The inline delayed-shot migration kernel mixes source and receiver wavefield energy, $S_{jl}(\bE|\w)$ and $R_{kl}(\bE|\w)$, according to   
\begin{eqnarray} \label{eqn:PWimage}
I(\bE) & = & \sum_{l=1}^A \sum_{p_{\xi_1}} \sum_{j=1}^B
\sum_{k=1}^B I^{DS}_{jkl}(\bE| p_{\xi_1}) \, \\ \, & = &
\sum_{l=1}^A \sum_{p_{\xi_1}} \sum_{j=1}^B \sum_{k=1}^B
\sum_{\w} |f(\w)|^2 \, S_{jl}^*(\bE|\w) \,R_{kl}(\bE|\w)
{\rm e}^{ i \w[ p_{\xi_1} \Delta \xi_1 (j-k)]}, \nonumber 
\end{eqnarray}
Similar to plane-wave migration, mixing wavefields of differing $S_{jl}$ and $R_{kl}$ indices will introduce crosstalk into the image volume.  However, inline delayed-shot migration will be crosstalk-free in the following limit: 
\begin{equation}\label{eqn:Dapprox2}
\lim_{N_{p_{\xi_1}} \rightarrow \infty} \sum_{\alpha = -N_{p_{\xi_1}}}^{N_{p_{\xi_1}}}
{\rm e}^{ i \w\alpha \Delta p_{\xi_1} \Delta {\xi_1} (j-k) } = | \w |^{-1} \delta_{jk},
\end{equation} 
Defining $|f(\w)|^2=|\w|$ and using the approximation in equation~\ref{eqn:Dapprox2}, I rewrite 
\begin{equation} \label{eqn:IPWimage}
I_l^{DS}(\bE) \approx \sum_{j=1}^B  \sum_{\w}S_{jl}^*(\bE|\w) R_{jl}(\bE|\w).
\end{equation}
Stacking over all inline delayed-shot sail-line migration results
yields the full image volume, 
\begin{equation} \label{eqn:CWimage2}
I(\bE) \approx \sum_{l=1}^A I_l^{DS}(\bE)  \approx\sum_{l=1}^A \sum_{j=1}^B
\sum_{\w}S_{jl}^*(\bE|\w) R_{jl}(\bE|\w). 
\end{equation} 
This proves the equivalence of inline delayed-shot and shot-profile migration.  
%
\section{Tilted elliptical-cylindrical coordinates}
%
One question to be addressed is what coordinate system geometry optimally conforms to the impulse response of a conical wavefield?  I assert that the best geometry is that of the TEC coordinate system shown in Figures~\ref{fig:TECgeom}~and~\ref{fig:TEC}.  One advantage is that the breadth of the first extrapolation step at the surface allows multiple streamers of a single sail line to be positioned directly on a single mesh.  Hence, this geometry is applicable to both narrow- and wide-azimuth acquisition.  A second advantage is that one direction of large-angle propagation can be handled by coordinate system tilting, while the other is naturally handled by the ellipticality of the mesh.  (Note that the geometry of another natural mesh - cylindrical polar coordinates - would not be a judicious choice for because the geometry permits migration of only single-streamer data and has singular points located on the surface at the first extrapolation step.)
%
\par
%
I set up the migration geometry of the elliptical-cylindrical mesh as follows:
\begin{itemize} 
\item $\xi_3 \in [0,\infty]$ is the extrapolation direction, where surfaces of constant $\xi_3$ form concentric elliptical cylinders, shown in Figure~\ref{fig:TECgeom}a.
\item $\xi_2 \in [0,2\pi)$ is the crossline direction, where surfaces of constant $\xi_2$ are folded hyperbolic planes, shown in Figure~\ref{fig:TECgeom}b; and 
\item $\xi_1 \in [-\infty,\infty]$ is the inline direction, where surfaces of constant $\xi_1$ are 2D elliptical coordinate meshes, shown in Figure~\ref{fig:TECgeom}c;\end{itemize}
%
\plot{TECgeom}{height=0.8\textheight}{Constant surfaces of the elliptical-cylindrical coordinate system (with zero inline tilt).  Cartesian coordinate axes are given by the vector diagram. a) constant $\xi_3$ surfaces forming confocal elliptical-cylindrical shells that represent the direction of extrapolation direction.  b) constant $\xi_2$ surfaces representing folded hyperbolic planes.  c) constant $\xi_1$ surfaces representing 2D elliptical meshes. {\bf NR}}
%
\par%
The mapping relationship between the two coordinate systems, adapted from \cite{Arfken}, is 
\begin{equation}
\label{eqn:2DPE1}
  \left[\begin{array}{c}
  	x_1 \\
  	x_2 \\
  	x_3 
  	\end{array}\right] 
  	= 
  \left[\begin{array}{c}
\xi_1 \,{\rm cos}\,\theta - a \, \rm{sinh}\,\xi_3\,\rm{sin}\,\xi_2 \,{\rm sin}\,\theta \\
a \,\rm{cosh} \, \xi_3 \, \rm{cos}\,\xi_2  \\ 
\xi_1\, {\rm sin}\, \theta + a \,\rm{sinh} \, \xi_3 \,\rm{sin}\,\xi_2\, {\rm cos}\, \theta\\  
  \end{array}\right],  
\end{equation}
where $\theta$ is the inline tilt angle of the coordinate system and parameter $a$ controls the coordinate system breadth.  Panels~\ref{fig:TEC}a and \ref{fig:TEC}b show the TEC coordinate system at $0^\circ$ and  $25^\circ$ tilt angles, respectively. 
%
\plot{TEC}{width=\textwidth}{Four extrapolation steps in $\xi_3$ of an TEC coordinate system, where the $\xi_1$  and $\xi_2$ coordinate axes are oriented in the inline and crossline directions, respectively. a) 0$^\circ$ tilt angle.  b) 25$^\circ$ tilt angle.  {\bf NR}}   
%
\subsection{TEC extrapolation wavenumber}
%
A metric tensor $g_{jk}$ can be specified from the mapping relationship given in equations~\ref{eqn:2DPE1}:   
\begin{equation}
\left[g_{jk}\right] =
\left[\begin{array}{ccc}
    1 & 0 & 0\\
    0 & A^2 & 0 \\
    0 & 0 & A^2 \\
\end{array}\right],
\end{equation}
where $A=a \sqrt{ {\rm sinh}^2 \,\xi_3 + {\rm sin}^2 \,\xi_2 }$.  The determinant of the metric tensor is: $|\boldsymbol{g}| = A^4$.  The corresponding inverse weighted metric tensor, $m^{jk}$ as developed in Chapter 2, is given by:  
\begin{equation}
\left[m^{jk}\right]
=
\left[\begin{array}{ccc}
    A^2 & 0 & 0 \\
    0 & 1 & 0 \\
    0 & 0 & 1 \\
\end{array}\right].
\end{equation}
Note that even though the metric of the TEC coordinate system varies spatially, the local curvature parameters ($n^j=\frac{\partial   m^{jk}}{\partial\xi_k}$) remain constant: $n^1=n^2=n^3=0$.   The corresponding extrapolation wavenumber, $k_{\xi_3}$, can be generated by inputting tensor $m^{jk}$ and fields $n^j$ into the general wavenumber expression for 3D non-orthogonal coordinate systems in equation \ref{eqn:kzfinal} 
%
\begin{equation} \label{eqn:EWN}
k_{\xi_3} = \pm \sqrt{ A^2 s^2 \w^2 - A^2 k_{\xi_1}^2 -  k_{\xi_2}^2 },
\end{equation}
where $s$ is the slowness (reciprocal of velocity), $k_{\xi_3}$ is the extrapolation wavenumber, and $k_{\xi_1}$ and $k_{\xi_2}$ are the inline and crossline wavenumbers, respectively.
%
\par
%
The wavenumber specified in equation~\ref{eqn:EWN} is central to the inline delayed-shot migration algorithm.  The first step is to extrapolate the  source and receiver wavefields
\begin{eqnarray}
E_{\xi_3}[S_{jl} (\xi_3,\xi_1,\xi_2|\w) ] 
& = & 
S_{jl} (\xi_3+\Delta \xi_3, \xi_1,\xi_2|\w), \\ 
E^*_{\xi_3} [R_{kl} (\xi_3,\xi_1,\xi_2|\w) ]
&=&
R_{kl} (\xi_3+\Delta \xi_3,\xi_1,\xi_2|\w),
\end{eqnarray}
where $E_{\xi_3}[\cdot]$ and $E^*_{\xi_3}[\cdot]$ are the extrapolation operator and its conjugate, respectively.  The results herein were computed using the $\omega-\xi$ finite-difference extrapolators discussed below.  The second step involves summing the individual inline delayed-shot images contributions, $I^{DS}_{jk} (\bE)$, into the total image volume, $I (\bE)$ according to equation~\ref{eqn:CWimage2}.
%
\section{3D Implicit finite-difference extrapolation}
%
One obvious concern is whether the dispersion relationship in equation~\ref{eqn:EWN} can be implemented accurately and efficiently in a wavefield extrapolation scheme.  I address this question by comparing the elliptical-cylindrical dispersion relationship to that for elliptically anisotropic media in Cartesian coordinates.  By defining an effective slowness $s_A=A s$ and rewriting equation~\ref{eqn:EWN} as  
\begin{equation} \label{eqn:FDstyle}
\frac{k_{\xi_3}}{\omega s_A} = 
\sqrt{ 1 - A^2\frac{k_{\xi_1}^2}{\omega^2 s_A^2} - \frac{
    k_{\xi_2}^2 }{\omega^2 s_A^2}}, 
\end{equation}
the TEC coordinate dispersion relationship resembles that of elliptically anisotropic media \cite[]{Tsvankin1996}.   More specifically, extrapolation in TEC coordinates is related to a special case where the Thomsen parameters \cite[]{Thomsen1986} obey  $\epsilon=\delta$:
\begin{equation} \label{eqn:Aniso}
\left. \frac{k_{x_3}}{\omega s} \right|_{\epsilon=\delta} = 
\left. \sqrt{
\frac{1 - (1+2\epsilon) \frac{k_{x_1}^2+k_{x_2}^2}{\omega^2 s^2}}
{1-2(\epsilon-\delta) \frac{k_{x_1}^2+k_{x_2}^2}{\omega^2 s^2}}
}\right|_{\epsilon=\delta}=
\sqrt{ 1-(1+2\epsilon) \frac{k_{x_1}^2}{\omega^2 s^2}-(1+2\epsilon)
\frac{k_{x_2}^2}{\omega^2 s^2}}.
\end{equation}
%
From equation~\ref{eqn:Aniso} we see that equation~\ref{eqn:EWN} is no more complex than the dispersion relationship for propagating waves in elliptically anisotropic media, which is now routinely handled with finite-difference approaches \cite[]{Zhang2001,Baumstein2003, Shan2005}.  
%
\par
%
A general approach to 3D implicit finite-difference propagation is to approximate the square-root  by a series of rational functions \cite[]{Ma}
\begin{equation} \label{eqn:expansion}
S_{\xi_3} = \sqrt{1 - A^2 S_{\xi_1}^2 - S_{\xi_2}^2} \approx
\sum_{j=1}^{n} \frac{ a_j S_r^2}{1-b_j S_r^2}, 
\end{equation}
where $S_{\xi_j}= \frac{k_{\xi_j}}{\omega s_A}$ and $S_r^2=  A^2S_{\xi_1}^2 + S_{\xi_2}^2$,  for $j=1,2,3$, and $n$ is the order of the coefficient expansion.  
%
\par
%
An optimal set of coefficients can be found by solving an optimization problem \cite[]{Shan2005}, 
\begin{equation}
E(a_j,b_j) = {\rm min} \int_{0}^{{\rm sin} \phi} 
\left[ \sqrt{1-S_r^2} - \sum_{j=1}^{n} \frac{a_j S_r^2}{1 - b_j S_r^2}
  \right]^2 {\rm d}S_r, 
\end{equation}
where $\phi$ is the maximum optimization angle.  I generated the following results using a 4th-order approximation and coefficients found in Table~\ref{tab:coeffs} \cite[]{GEO50-10-16341637}. 
%
\begin{table}
   \centering
   \begin{tabular}{|c|c|c|}
      \hline
      Coeff. order $j$ & Coeff. $a_j$ &  Coeff. $b_j$ \\
      \hline
      1     & 0.040315157 &  0.873981642\\
      2      & 0.457289566 & 0.222691983 \\
     \hline
   \end{tabular}
   \caption{Coefficients used in 3D implicit finite-difference wavefield extrapolation.}  
   \label{tab:coeffs}
\end{table} 
%
\par
%
Specifying a finite-difference extrapolator operator using the 4th-order approximation is equivalent to solving a cascade of partial differential equations \cite[]{Shan2005} 
\begin{eqnarray}
\frac{\partial}{\partial \xi_3} U_{\xi_3+\Delta\xi_3/3}  = & i \omega s U_{\xi_3}, \nonumber \\ 
\frac{\partial}{\partial \xi_3} U_{\xi_3+2\Delta\xi_3/3}  = & i \omega s 
\left[
\frac{ \frac{a_1}{\omega^2s^2} \frac{\partial^2}{\partial
    \xi_1^2} }{1+\frac{b_1}{\omega^2s^2}
  \frac{\partial^2}{\partial \xi_1^2} }+\frac{
  \frac{a_1}{\omega^2s_A^2} \frac{\partial^2}{\partial \xi_2^2}
}{1+\frac{b_1}{\omega^2s_A^2} \frac{\partial^2}{\partial \xi_2^2}
} 
\right] U_{\xi_3+\Delta\xi_3/3}, \\
\frac{\partial}{\partial \xi_3} U_{\xi_3+\Delta\xi_3}  = & i \omega s 
\left[
\frac{ \frac{a_2 }{\omega^2s^2} \frac{\partial^2}{\partial
    \xi_1^2} }{1+\frac{b_2 }{\omega^2s^2}
  \frac{\partial^2}{\partial \xi_1^2} }+\frac{
  \frac{a_2}{\omega^2s_A^2} \frac{\partial^2}{\partial \xi_2^2}
}{1+\frac{b_2}{\omega^2s_A^2}\frac{\partial^2}{\partial \xi_2^2} }  
\right] U_{\xi_3+2\Delta\xi_3/3}. \nonumber
\end{eqnarray}
I solve these equations implicitly at each extrapolation step by a finite-difference splitting approach that alternatively advances the wavefield in the $\xi_1$ and $\xi_2$ directions.  Splitting methods allow the direct application of the $A$ scaling factor in equation~\ref{eqn:expansion} by introducing the original slowness model, $\frac{s_A}{A}=s$, for the $\xi_1$ direction split. 
%
\par
%
One drawback to finite-difference splitting methods is that they commonly generate numerical anisotropy.  To minimize these effects, I apply a Fourier-domain phase-correction filter $L[\cdot]$ \cite[]{LiFilter}
\begin{equation}
L[U] = U{\rm e}^{ i \Delta \xi_3 k_L},
\end{equation}
 where
 \begin{equation}
 k_L = 
 \sqrt{ 1-\frac{k_{\xi_1}^2}{(\omega s_1^r)^2}
 -\frac{k_{\xi_2}^2}{(\omega s_2^r)^2}} - 
 \left[ 1 - \sum_{j=1}^{2}
 \left( \frac{ a_j ( \frac{k_{\xi_1}}{\omega s_1^r})^2}{1-b_j (
 \frac{k_{\xi_1}}{\omega s_1^r})^2} - 
 \frac{ a_j ( \frac{k_{\xi_2}}{\omega s_2^r})^2}{1-b_j (
 \frac{k_{\xi_2}}{\omega s_2^r})^2}  
 \right) \right],
 \end{equation}  
 and $s_1^r$ and $s_2^r$ are reference slownesses chosen to be the
 mean value of $s^A_{eff}$ and $s$ defined above, respectively.  Note that while this phase-shift correction is explictly correct for v($\xi_3$) media, the Li filter in v($\xi_3,\xi_1,\xi_2$) media is only approximate and will introduce error.
 %
\subsection{Impulse response tests}
%
I conducted impulse response tests on a 500x400x400 cube in a homogeneous medium of slowness $s=0.0005$~sm$^{-1}$.  The initial wavefield consisted of three horizontally smoothed point sources at $t$=0.5, 0.75, and 1~s. The impulse responses are expected to consist of three hemispherical surfaces of radii $r$=1000, 1500, and 2000~m. 
%
\par
%
Figures~\ref{fig:CrossIn}a and \ref{fig:CrossIn}b show the inline and crossline responses.   The three lines overlying the analytic curves show the correct impulse response locations.  Note that the impulse responses are restricted at large angles both by the coordinate system boundaries and by the 50 sample cosine-taper function along the edges of the TEC mesh. Figure~\ref{fig:Depth1300} shows a  impulse response slice extracted at 1300~m depth.  The symmetric response indicates that the numerical anisotropy from the numerical splitting is accounted for by the Li phase-correction filter. 
%
\plot{CrossIn}{width=4.5in,height=7.2in}{Elliptical-cylindrical-coordinate
  impulse-response tests.  a) Inline section. b) Crossline section.  {\bf CR}} 
%
\sideplot{Depth1300}{width=3in,height=3in}{Elliptical-cylindrical impulse response at 1300~m depth.  Note the circular symmetry of the impulse response indicating little-to-no numerical anisotropy. {
\bf CR}}  
%
\section{3D Wide-azimuth synthetic tests}
%
This section presents the inline delayed-shot migration algorithm test results on a wide-azimuth synthetic data set generated from a realistic 3D Gulf of Mexico velocity model.  Figure~\ref{fig:VELCUT} presents some depth slices and sections through the model.  The velocity model is comprised of typical Gulf of Mexico sedimentary profile with a velocity gradient of  approximately $0.2$~s$^{-1}$, with a number of salt bodies of complex 3D geometry characterized by smoothly varying salt tops and steep flanks below overhangs.  Key imaging targets include the steep salt flanks and the on-lapping sedimentary units that comprise the likely exploration areas.
%
\plot{VELCUT}{width=0.7\textwidth,height=0.875\textheight}{Depth section and inline crossline sections of the Gulf of Mexico velocity model through complex 3D salt bodies.  Top: 3900 m depth slice.  Middle: 33000 m inline section. Top: 16000 m crossline section.  {\bf ER}}
%
\par
%
Table~\ref{tab:DSvals} summarizes the acquisition geometry of the data set.  The data used for migration consisted of 72 sail lines separated 250 m apart.  Each sail line consists of 100 shots sampled at a 250 m shot interval.  The receiver pattern for each shot record contains 321 inline samples with a maximum offset of $\pm$8000 m computed at a 50 m interval, and 161 crossline samples with a maximum offset of $\pm$4000 m at a 50 m interval.  
%
\begin{table}
   \centering
   \begin{tabular}{|c|c|c|c|}
      \hline
      Source Parameter & Value  &  Receiver Parameter & Value \\
      \hline
      Number of sail lines & 72 & Max. inline offset (m) & $\pm$8000 \\
      Sail line interval (m) & 250 &  Max. crossline offset (m) & $\pm$4000 \\
      Shots per sail line & 100 & Inline receiver interval (m)& 50 \\
      Shot interval (m)& 250 & Crossline receiver interval (m)& 50\\
     \hline
   \end{tabular}
   \caption{Parameters associated with the 3D synthetic data set.}  
   \label{tab:DSvals}
\end{table} 
%
\par
%
A total of 192 frequencies were selected for migration starting at 1.42~Hz at a sampling rate of 0.075~Hz.  Filtered data from each sail line data were transformed into a plane-wave data set by phase-encoding over a range of inline ray parameters, $p_{\xi_1}$.   I selected a total of 101 inline ray parameters between $\pm$8.33x10$^{-4}$~sm$^{-1}$ at a sampling rate of 8.33x10$^{-6}$~sm$^{-1}$.  Given the $1500$~ms$^{-1}$ water velocities at the surface, the maximum values correspond to a surface take-off angle of $\pm 38.7^\circ$.
%
\par
%
I applied the inline delayed-shot migration technique to the plane-wave data on a sail-line by sail-line basis, which allowed for a coarse-grain computational parallelism at a scripting level.   (The  migration code was also OMP-enabled, which led to a second level of coarse-grain parallelism over the frequency axis.)  Migration runs were conducted for Cartesian coordinate (CC) and TEC geometries with both tilting and non-tilting meshes.  For CC migrations, the data volumes were zero-padded with 40 samples on each inline side and 95 samples on each crossline side.  The data volume for TEC migrations were padded with 40 samples on the inline sides, but only one sample on each crossline side because the coordinate system aperture expands naturally in the crossline direction. 
%
\par
%
Figure~\ref{fig:TESTY15400} presents the 15400~m cross section from the 24500~m sail-line migration image (for 101 plane-waves) for the TEC (top panel) and the CC (bottom panel) geometries.  The gently dipping sedimentary reflections in both sections are imaged across a ~6000~m swath.  The TEC migration, relative to that in CC geometry, shows a significant improvement in the vertical salt flank on the right-hand-side of the image.  Although the salt-flank is weakly present in the CC image under strong clipping, it is mis-positioned due to the ~$80^\circ$  limit of extrapolation operator accuracy.
%
\par
%
%Figure~\ref{fig:TESTX22250} presents the 22250~m inline section for the 22250~m sail line migration image.  Note that because the coordinate system has not been tilted in the inline direction in this example, the EC coordinate (top panel) and the CC  (bottom panel) images are very similar.  Both images show well-imaged salt tops and lower flank to the left.  However, the near-vertical salt flanks are absent.
%
\plot{TESTY15400}{width=\textwidth,height=0.35\textheight}{Sections for the 24500~m sail line at the 15400m crossline coordinate.  Top: Elliptical-cylindrical coordinate image.  Bottom: Cartesian coordinate image. {\bf NR}}
%
%\plot{TESTX22250}{width=\textwidth,height=0.35\textheight}{Sections for the 22250~m sail line at the 22250m inline coordinate.  Top: Elliptical-cylindrical coordinate image.  Bottom: Cartesian coordinate image.  }
%
\par
%
Figures~\ref{fig:FIG3}~and~\ref{fig:FIG4} present crossline sections from the full TEC and CC image volumes.   Figure~\ref{fig:FIG3} presents the EC and Cartesian crossline sections at the 33700~m inline coordinate in the upper and lower panels, respectively.  The TEC image has an improved left-hand salt flank (marked A) that is more correctly positioned relative to the CC image.  Similarly, the right-hand salt flank (marked B) is more accurately positioned and forms a more continuous reflector. 
%
\plot{FIG3}{width=\textwidth}{Crossline sections through the velocity model and full image volumes at inline coordinate 33700~m.  Top panel: Elliptical-cylindrical coordinate image.  Bottom panel: Cartesian coordinate image.  The imaging improvements for the left-hand salt flank are denoted by the oval marked A.  The oval marked B illustrates a more continuous and correctly placed reflector in the TEC coordinate system.  {\bf NR}}
%
Figure~\ref{fig:FIG4} presents the 15100~m crossline sections extracted from the two image volumes.  Note the differences in the vertical right-hand salt flank (marked A) between the ECC (top panel) and CC (bottom panel) images.  The TEC image exhibits a stronger reflector that is better positioned than that in the CC image (again because of the high-angle limits of the extrapolation operator).  
%
\plot{FIG4}{width=\textwidth}{Crossline sections through the velocity model and full image volumes at inline coordinate 15100~m.  Top panel: Elliptical-cylindrical coordinate image.  Bottom panel: Cartesian coordinate image.  The oval marked A indicates the location of the vertical salt flank that is better imaged in TEC coordinates.  {\bf NR}}
%
\par
%
Figure~\ref{fig:FIG1} shows the 21750~m inline section through the complete TEC (upper panel) and CC (lower panel) image volumes.  The left-hand salt flank (marked A) is more accurately located and continuous in the TEC image.  The right-hand salt flank (marked B), again, is more continuous in the 
TEC image.  Another observation is that the TEC image (and in Figures~\ref{fig:FIG3}-\ref{fig:FIG4}) does not contain the same spatial frequency content as the CC images (see below).
%
\plot{FIG1}{width=\textwidth}{Inline sections through the velocity model and full image volumes at crossline coordinate 21750~m.  Top panel: Velocity section.  Middle panel: Elliptical-cylindrical coordinate image.  Bottom panel: Cartesian coordinate image.  The left-hand salt flank, shown in oval A, is more accurately positioned in the TEC coordinate image, while the right-hand flank, marked by oval B, is similarly more accurately positioned and continuous.  {\bf NR}}
%
\par
%
Figure~\ref{fig:FIG5} presents slices extracted at 6150~m depth from the TEC (top panel) and CC (bottom panel) images.  The images are again fairly similar, though are there slight differences that correspond to amplitude differences between the weakly imaged steep flank reflectors.  Examples include the regions marked A and C that corresponds to the salt flanks in Figure~\ref{fig:FIG1} and Figure~\ref{fig:FIG3}, respectively.  Finally, the migration algorithm has well-imaged the set of channels denoted in region B in both coordinate system images.
%
\plot{FIG5}{height=0.8\textheight}{Depth slices through the velocity model and image volumes at 6150~m depth.  Top panel: Elliptical-cylindrical coordinate image.  Bottom panel: Cartesian coordinate image.  Oval A illustrates the improved TEC image for the vertical salt flank shown in Figure~\ref{fig:FIG4}.  Oval B demarcates a region where some of the smaller-scale fractures are well imaged in both images.  Oval C shows the region where the near-vertical flank shown in TEC coordinate image in Figure~\ref{fig:FIG3} is better imaged.  {\bf NR}}
%
\subsection{Discussion}
%
Relative computational cost is one important metric to consider when comparing the migration algorithms in different coordinate systems.  In the above tests, padding in the crossline direction tended to be the most important factor in determining the migration run time.  One benefit of the TEC geometry is its naturally outward-expanding mesh in the elliptical direction that effectively increasing the migration aperture.  Thus, TEC migrations usually require less zero-padding in the crossline direction relative to CC geometries.  I performed the TEC migrations on meshes with inline-by-crossline-by-depth grids of 720x324x400.  Migrations in CC geometries required a 720x512x400 mesh in order to achieve similar crossline aperture, which resulted in a fairly significant additional computational overhead.  
%
\par
%
Table~\ref{tab:runtimes} shows the comparative costs for various TEC and CC migration runs for both the shot-profile and delayed-shot migration styles.  I used 72 data points in specifying each median runtime times for the four different migration runs.  The test migrations indicate that the TEC geometry migrations were faster than the those in Cartesian tests (for equivalent effective aperture), with 29\% and 21\% computational cost reduction for the shot-profile and inline delayed-shot migration strategies, respectively.
%
%\sideplot{Runtimes}{width=\textwidth}{Run-time comparisons for the shot-profile and delayed-shot migration styles in the elliptical-cylindrical and Cartesian coordinate systems.}
\par
\begin{table}
   \centering
   \begin{tabular}{|c|c|c|}
      \hline
      Migration & Coordinate  &  Mean run \\
       style &  system &  time (hrs) \\
      \hline
      Plane-wave & Tilted elliptical cylindrical & 37.2  \\
      Plane-wave & Cartesian & 45.0 \\
      Shot-profile & Tilted elliptical cylindrical & 15.5 \\
      Shot-profile & Cartesian & 20.0 \\ 
      \hline
   \end{tabular}
   \caption{Run-time comparisons for the delayed-shot migration and shot-profile styles in the tilted elliptical-cylindrical and Cartesian coordinate systems.}  
   \label{tab:runtimes}
\end{table} 
%

%
One question worth addressing is how far can the TEC sampling be reduced before imaging artifacts become apparent?  As one moves outward between successive extrapolation surfaces, the TEC geometry expands  at increasingly larger step sizes.  Fortunately, most realistic velocity models have velocity increasing with depth, causing the wavelengths of the propagated waves to lengthen.  This phenomenon acts as a natural wavefield filter that, in most cases, prevents wavenumbers from aliasing (except near-surface in the grid extremities).  A good rule-of-thumb is that one must ensure that the grid point of TEC coordinate system mesh does not go below one grid point for every two CC grid points in each direction; however, maintaining this relationship throughout the image volume is not a straightforward task.  Additional work on the craft of 3D coordinate-system interpolation is necessary and would likely help restore some of the absent high frequency information.
%
\par
%
An additional consideration of parameter choice is the interpolation window over which the surface wavefields are injected onto the TEC coordinate mesh.  Not using a sinc-based interpolation over the near-surface depth axis can lead to significant artifacts; however, choosing too large of a window will blend information from different extrapolation steps leading to smoother and lower frequency images.  Figures~\ref{fig:FIG3}-\ref{fig:FIG5} show the result of a somewhat overcautious parameter choice (interpolating wavefields three additional depth steps) that led to the lower spatial wavenumber content of the TEC images relative to the CC images.  I assert that his effective low-pass filtering can be reduced by interpolating only one or two additional steps in depth.
%
\section{Narrow-azimuth field data test}
%
This section presents the results of applying the inline delayed-shot imaging procedure to a 3D Gulf of Mexico narrow-azimuth data set provided by ExxonMobil.  The velocity model, shown in Figure~\ref{fig:XOM-VEL}, consists of typical sediment-controlled $v(\xi_3)$ velocity structure, save for the salt body intruding in the center of the block, and offset associated with throw along the dipping fault plane observable in Figure~\ref{fig:XOM-VEL}.  By agreement with ExxonMobil, the depths shown in all figures differ from the true values.
%
\plot{XOM-VEL}{width=\textwidth}{Velocity model example for Gulf of Mexico field data set.  {\bf ER}}
%
Key imaging targets in this model include the steep salt flanks around the salt structure and the onlapping sediments.  
%
Previous imaging work in this area indicates that the sediments surrounding the salt body exhibit moderate-to-strong degrees of anisotropy.  \cite{Bear} estimated the vertical velocity and anisotropy parameters (assuming VTI media) using a joint inversion technique that combined surface seismic and borehole constraints. \cite{ShanThesis}, using a 3D tilted Cartesian coordinate plane-wave migration algorithm for transversely isotropic (TTI) media, demonstrated that accounting for anisotropy greatly improves migrated image quality for this data set.   
%
\par
%
The migration strategy presented herein differs from that in \cite{ShanThesis} in a number of respects.  First, I perform migration using only isotropic sediment flood velocity model.  Second, I use a multi-streamer data set for imaging, rather than the more optimally regularized single-streamer data formed through azimuthal move-out preprocessing \cite[]{Biondi.3DSI}.  
%
\par
%
Table~\ref{tab:FieldParams} summarizes the acquisition geometry of the data set.  The data used for migration consisted of 54 sail lines separated roughly 160~m apart, each sail line consists of approximately 300 shots acquired every 50~m. 
%
\begin{table}
   \centering
   \begin{tabular}{|c|c|c|c|}
      \hline
      Source Parameter & Value  &  Receiver Parameter & Value \\
      \hline
      Number of sail lines & 54 & Max. inline offset (m) & $\pm$3725 \\
      Sail line interval (m) & 160 &  Max. crossline offset (m) & $\pm$500 \\
      Shots per sail line & 300 & Inline receiver interval (m)& 25 \\
      Shot interval (m)& 50 & Crossline receiver interval (m)& 80\\
     \hline
   \end{tabular}
   \caption{Approximate acquisition parameters associated with the 3D Gulf of Mexico field data set.}  
   \label{tab:FieldParams}
\end{table} 
%
I binned the sources in 25~m and 80~m intervals in the inline and crossline, respectively.  Figure~\ref{fig:XOM-SRC} shows the source distribution, and illustrates the sail line direction, herein chosen to be the inline direction.  Figure~\ref{fig:XOM-RCV} shows the chosen offset distribution.  The receiver points fall to both positive and negative offsets, as the sail lines were acquired in two directions.  The gap in offset coverage between offsets of $\pm$2500-2750~m arises due to a corrupted data tape.  Receivers were binned at 25~m in both the inline and crossline directions.
%
\sideplot{XOM-SRC}{width=\textwidth}{Chosen source distribution for the field data set. {\bf CR}}
%
\sideplot{XOM-RCV}{width=\textwidth}{Chosen receiver distribution for the field data set.  The missing data between offsets $\pm$2500-2750~m is due to a corrupted data tape. {\bf CR}}
%f
\par
%
I prepared the data for migration by applying an inline delay-shot phase-encoding algorithm according to the inline source position.  A total of 54 plane-wave sub-volumes were generated from the total 5D shot record volume, each consisting of 41 plane-waves equally sampled between $\pm 20^\circ$.  I chose a total of 244 frequencies between 3~Hz and 25~Hz for migration.  The data were imaged on migration grids with dimensions of 800x350x300 samples. Migrations in TEC coordinates were performed using tilt angles between $\pm 20^\circ$ at $1^\circ$ increments, where the tilt angle was matched with the same plane-wave take-off slowness.
%
\par
%
Figures~\ref{fig:RFIG1}-\ref{fig:RFIG5} present comparative slices from the 3D Gulf of Mexico migration images computed in the TEC and Cartesian coordinate systems.  Figure~\ref{fig:RFIG1} presents an inline section taken at the constant 8750~m crossline coordinate for the TEC (top panel) and CC (bottom panel) images corresponding to the front face of Figure~\ref{fig:XOM-VEL}.  The top of salt body is well-imaged in both images; however, the near-vertical salt-flanks to the right are nearly entirely absent.  Oval A shows the imaging improvements in TEC coordinates for the left-hand flank.
%
\plot{RFIG1}{width=\textwidth}{Inline sections through the migration images taken at the 8750~m crossline coordinate location.  Top: TEC coordinate migration results.  Bottom: Cartesian coordinate migration results.  Oval A shows the imaging improvements in TEC coordinates for the left-hand flank.  {\bf CR}}
%
Figure~\ref{fig:RFIG3} presents crossline sections for the TEC (top panel) and CC (bottom panel) images.  Both images are subject to fair amounts of near-surface aliasing a sali line contribution every 160~m in the crossline direction, as well as the artifaces due to using a non-regularized data set with locations where no data are present.  The outlines of the salt body reflector, though, are imaged.  Oval A shots an example of an area where the TEC coordinate image is better than the Cartesian image.
%
\plot{RFIG3}{width=\textwidth}{Crossline sections through the migration images taken at the 7100~m crossline coordinate location.  Top: TEC coordinate migration results.  Bottom: Cartesian coordinate migration results.  Oval A shows an example of an area where the TEC coordinate image is better than the Cartesian image.  {\bf CR}}
%
\par
%
Figure~\ref{fig:RFIG5} presents a depth slice extracted from the TEC (top panel) and CC (bottom panel) image volumes.  The annular ring, showing the location of the salt body, is apparent in both images; however, the image is sharper in the TEC image indicating improved focussing of energy.  Oval A shows an example of an area where the TEC image is better than that generated in Cartesian, including two parts of the right-hand salt flank.  Oval B shows the TEC coordinate image improvements in the crossline direction.
%
\plot{RFIG5}{width=\textwidth,height=6.5in}{Migration results for the 3D Gulf of Mexico field data set through the sedimentary section.  Top: TEC coordinate migration results.  Bottom: Cartesian coordinate migration results.  Oval A shows an example of an area where the TEC image is better than that generated in Cartesian, including two parts of the right-hand salt flank.  Oval B shows the TEC coordinate image improvements in the crossline direction.  {\bf CR}}
%
\par
%
The results of the 3D field data application likely could have been improved in a number of aspects.  First, a migration velocity model incorporating anisotropy information to emphasize horizontal wave propagation could have been used instead of the isotropic velocity profile.  Although this would affect the vertical location of the flat-lying sedimentary reflectors, it likely would have led to more accurate horizontal propagation and imaging of waves reflecting off the target salt flanks.  Second, if additional computational resources were made available, migrating the full data set (i.e. every 80~m in crossline source position rather than every 160~m) with a higher frequency content would have led to a more infilled and higher resolution image.  Third, extending the generalized RWE theory to incorporate TTI anisotropy likely would have enabled a more consistent imaging of the steep salt flanks.  This extension is likely to be a subject for future research.
%
\section{Conclusions}
%
This chapter discusses an inline delayed-shot migration technique in tilted elliptical-cylindrical coordinates.  I argue that migration approach, relative to the full 3D plane-wave technique, offers both lower memory requirements (due to small migration aperture), as well as a potential reduction in the number of total migrations needed (by migrating fewer sail lines than crossline plane waves).  I demonstrate that the impulse response of inline-source delayed-shot wavefields are well-matched to TEC geometry, and that corresponding extrapolation wavenumber is no more complicated than that of elliptically anisotropic media.  This leads to an accurate 3D finite-difference splitting algorithm that both accurately propagates wavefields and handles the associated numerical anisotropy.  The 3D synthetic Gulf of Mexico data tests demonstrate the migration technique's ability to generate improved images of steeply dipping structure, relative to Cartesian coordinate migration, at reduced computational cost.  Field data tests illustrate the utility of the 3D migration approach in exploration practice.
%
\section{Acknowledgments}
%
I thank Biondi Biondi, Bob Clapp, Alejandro Valenciano, Ben Witten, and Bill Curry for helpful conversations.  John Etgen is greatly thanked for donating his wisdom and time on the BP research computer cluster, where the 3D synthetic imaging results were generated.  I acknowledge BP for the release of the synthetic data set, and ExxonMobil for releasing the field data set.
%
\cleardoublepage
%\bibliographystyle{seg}
%\bibliography{SEP,SEG,MISC,GEOTLE,allrefs}

