\def\figdir{./Fig/}
\lefthead{Shragge}
\righthead{Non-orthogonal RWE}
\footer{SEP--124}
\title{Non-orthogonal Riemannian wavefield extrapolation}
\email{jeff@sep.stanford.edu}
\keywords{}
\author{Jeff Shragge}

\input{mainmacros_act}
\input{macros_act}

\ABS{ This paper generalizes wavefield extrapolation to
non-orthogonal Riemannian spaces.  Central to the development is the 
specification of an extrapolation wave-number for propagating wavefields
on non-orthogonal grids. The expression contains additional
non-stationary coefficients, and generates further complexity to
existing fields.  However, I show that an extended split-step Fourier
approach provides an acceptable approximation.  Importantly, rather
than attempting to satisfy semi-orthogonal grid criteria, this
development allows the user to focus on generating smoother and more
stationary coordinate meshes. }
%
\section{Introduction}
%
Riemannian wavefield extrapolation (RWE) \cite{Sava.geo.rwe}
generalizes the wavefield extrapolation component of wave-equation
migration to non-Cartesian coordinate systems. The original
formulation assumes that coordinate systems are (semi-)orthogonal with
an extrapolation direction orthogonal to the other two axes.   This
supposition resulted in a wave-equation 
dispersion relationship for the extrapolation wavenumber that 
contained non-stationary coefficients (additional to propagation slowness)
describing the coordinate system geometry.  A wavefield extrapolation
scheme could then be developed using a multi-coefficient version of
the split-step Fourier approach \cite{GEO55-04-04100421}.  However, in
many cases semi-orthogonal geometry turns out to be a fairly strong
restriction for two reasons: i) generating semi-orthogonal coordinate
meshes with appropriate characteristics is often onerous; and ii)
developed semi-orthogonal meshes generated greatly varying
non-stationary coefficients. 
%
\par
%
RWE was designed initially for 'dynamic' applications where a 
coordinate system orientation is conformal to the wavefield propagation
direction.  In many situations, this approach generated high-quality Green's
functions using ray-traced coordinate systems originating at the location of an
impulsive source point.  However, numerical instability (i.e. zero
divisions) occurred wherever ray coordinate systems went through
triplications.  \longcite{Sava.geo.rwe} addressed this issue by
iteratively smoothing the velocity models until all coordinate system
triplications disappeared.  This {\it ad hoc} solution, though, was
somewhat self-defeating because it countered the original purpose of
coordinate systems and propagation direction conformality.    
%
\activeplot{Gauss\DOT ps}{width=6in}{NR}{Gaussian slowness anomaly illustrating the
  problem with triplicating coordinate systems}
%
\par
%
Another application of RWE in performing imaging to and from surfaces
with irregular geometry.  \longcite{Shragge.sep.117.jeff1} formulated
a wave-equation migration from topography strategy that posed
wavefield extrapolation directly in locally orthogonal coordinates
orthogonal to the acquisition surface.  Though successful in areas
with longer wavelength and lower amplitude topography, the imaging
results degraded in situations involving more rugged acquisition
topography.  A more significant observation, though, is that
coordinate system compression/extension demanded by
(semi)-orthogonality caused the degradation. 
%
\activeplot{Foothills\DOT coords}{width=6in,height=3in}{CR}{Topographic coordinate system
  from the Canadian Rocky Mountain Foothills.  Note the
  non-stationarity in the coordinate system compression and extension
  caused by the orthogonality constraint.}
%
\par
%
A solution to both of the above problems is to extend RWE to
non-orthogonal coordinate systems.  This can be achieved by
eliminating the semi-orthogonal constraint and applying the same
technique to a generalized Riemannian space.  Non-orthogonality
introduces two additional terms in the RWE dispersion relationship,
and makes existing coefficient terms slightly more complicated.
However, the resulting dispersion relationship and extrapolation
wavenumber can be handled with the same extended split-step Fourier
approach.  Importantly, this solution affords greater flexibility in
coordinate system design while facilitating more rapid coordinate
system computation.  In addition, instead of focusing on local and/or
global orthogonality constraints, greater emphasis on optimizing grid
quality can help control grid clustering and generate more stationary
an smoother coefficients. 
%
\par
%
This paper develops the wave-equation dispersion relationship
appropriate for 3D non-orthogonal RWE, following the approach of
\longcite{Sava.geo.rwe}.  I first discuss generalized Riemannian
geometry and then show how the acoustic wave-equation can be
formulated in a non-orthogonal Riemannian space.  Subsequently,
expressions for one-way wavefield extrapolation wavenumber and its
split-step Fourier approximation are developed.  I then  present two
analytic 2D non-orthogonal coordinate systems and verify developed
extrapolation wave-number expressions with numerical tests.    
%
\section{Generalized Riemannian Geometry}
%
Geometry in a generalized Riemannian space is described by a symmetric
metric tensor, $g_{ij}=g_{ji}$, that relates the geometry in a
non-orthogonal coordinate system, $\{\xi_1,\xi_2,\xi_3\}$, to an
underlying Cartesian system,  
$\{x_1,x_2,x_3\}$  \cite{Guggenheimer}.  The metric tensor is written
in the following matrix form,   
\begin{equation}
\label{eqn:metric}
\left[g_{ij}\right] 
= 
\left[\begin{array}{ccc}
g_{11} & g_{12} & g_{13} \\
g_{21} & g_{22} & g_{23} \\
g_{31} & g_{32} & g_{33} 
\end{array}\right] 
= 
\left[\begin{array}{ccc}
g_{11} & g_{12} & g_{13} \\
g_{12} & g_{22} & g_{23} \\
g_{13} & g_{23} & g_{33} 
\end{array}\right] \end{equation}
where $\EE$, $\FF$, $\GG$, $\HH$, $\LL$ and $\AA$ are functions
relating the two geometries through the following partial differential
operations, 
\beqa \label{eqn:pderels}
\EE=\hxxdef, \;\;\;\; \FF=\hxydef, \;\;\;\; \GG = \hyydef,\nonumber \\
\HH=\hxzdef, \;\;\;\; \LL=\hyzdef, \;\;\;\; \AA = \hzzdef. 
\eeqa
The convention used in equation~\ref{eqn:pderels}, and in the
following equations, is that repeated indicies indicate summation
(i.e. $g_{ii} = g_{11}+g_{22}+g_{33}$).  Metric tensor, $g_{ij}$, is
related to the associated (inverse) metric tensor, $g^{ij}$ through
$g_{ij}=\,|\mathbf{g}|\,g^{ij}$, where $|\mathbf{g}|$ is metric tensor
matrix determinant.  The associated metric tensor is  given by,  
\begin{equation}
\label{eqn:ametric}
  \left[g^{ij}\right] =\frac{1}{\left|\mathbf{g}\right|}
  \left[\begin{array}{ccc}
      \GG\AA -\LL^2 & \HH\LL-\FF\AA & \FF\LL-\HH\GG \\
      \LL\HH-\FF\AA & \EE\AA-\HH^2  & \HH\FF-\EE\LL \\
      \FF\LL-\GG\HH   & \FF\HH-\EE\LL   & \EE\GG-\FF^2 
    \end{array}\right]\;,\end{equation}
where $\PP$ is given by
\begin{equation}
\PP = \left[1 - \frac{\EE\LL^2+\GG\HH^2-2\FF\LL\HH}
{\AA\JJ}\right].\end{equation}
The metric determinant takes the form
\begin{equation}
\label{eqn:det}
  |\mathbf{g}| = \AA\,\JJ\,\PP^2\;.\end{equation}
Combining the expressions in equations~\ref{eqn:ametric}-\ref{eqn:det}
leads to the following expression for the associated metric tensor, 
\begin{equation}
\label{eqn:ametric2}
  \left[g^{ij}\right] =
  \left[\begin{array}{ccc}
      \frac{\GG\AA -\LL^2}{\gdet} & 
      \frac{\HH\LL-\FF\AA}{\gdet} & 
      \frac{\FF\LL-\HH\GG}{\gdet} \\
      & & \\
      \frac{\LL\HH-\FF\AA}{\gdet} &
      \frac{\EE\AA-\HH^2 }{\gdet} &
      \frac{\HH\FF-\EE\LL}{\gdet} \\
      & & \\
      \frac{\FF\LL-\GG\HH}{\gdet} &
      \frac{\FF\HH-\EE\LL}{\gdet} &
      \frac{\JJ          }{\gdet}
    \end{array}\right]\;.\end{equation}
A weighted metric matrix, $\mathbf{m}=\gdetf \, \mathbf{g}$,
\beq\label{eqn:ametric3}
  \left[m^{ij}\right] =
  \left[\begin{array}{ccc}
      \frac{\GG\AA -\LL^2}{\gdetf} & 
      \frac{\HH\LL-\FF\AA}{\gdetf} & 
      \frac{\FF\LL-\HH\GG}{\gdetf} \\
      & & \\
      \frac{\LL\HH-\FF\AA}{\gdetf} &
      \frac{\EE\AA-\HH^2 }{\gdetf} &
      \frac{\HH\FF-\EE\LL}{\gdetf} \\
      & & \\
      \frac{\FF\LL-\GG\HH}{\gdetf} &
      \frac{\FF\HH-\EE\LL}{\gdetf} &
      \frac{\JJ          }{\gdetf}
    \end{array}\right]\;.\end{equation}
is also useful in for developing an expression for the extrapolation wavenumber. 
%
\section{Acoustic wave-equation in 3D generalized Riemannian spaces}
%
The acoustic wave-equation for an acoustic wavefield, $\W$, in a
generalized Riemannian space is given by,
\begin{equation}
\label{eqn:helm}
  \dd \W = - \ww^2\ss^2 \lp \mathbf{x} \rp \W\end{equation}
where the $\ww$ is frequency, $\ss$ is the propagation slowness, and
$\dd$ is the Laplacian operator given by, 
\begin{equation}
\label{eqn:laplac2}
  \Delta \W = \frac{1}{\sqrt{|\mathbf{g}|}}\,
    \frac{\partial}{\partial \xi_i}\,
    \left(\,m^{ij}\,\frac{\partial \W}{\partial \xi_j}\right).\end{equation}
Substituting equation~\ref{eqn:laplac2} into equation~\ref{eqn:helm} generates a
Helmholtz equation appropriate for propagating waves through a 3D Riemannian space:
\begin{equation}
\label{eqn:we1}
\frac{1}{\gdetf} \de\qi \lp \mmm\ii\jj \done\W\qj \rp = - \ww^2\ss^2 \W\end{equation}
Expanding the derivative terms on the left side of
equation~\ref{eqn:we1} and multiplying through by $\gdetf$ yields,  
\begin{equation}
\label{eqn:helm1}
\done{\mmm\ii\jj}\qi \done\W\qj + \mmm\ii\jj \mtwo\W\qi\qj
= - \gdetf \ww^2 \ss^2 \W\end{equation}
Defining $n_j$ as,
\begin{equation}
n_j=\done{\mmm\ii\jj}\qi = 
\done{m^{1j}}\qx + 
\done{m^{2j}}\qy +
\done{m^{3j}}\qz,\end{equation}
leads to a more compact notation of equation~\ref{eqn:helm1}
\begin{equation}
\label{eqn:helm2}
n_j \done\W\qj + \mmm\ii\jj \mtwo\W\qi\qj
= -\gdetf \ww^2 \ss^2 \W\end{equation}
Developing a wave-equation dispersion relation is achieved by
replacing the partial differential operators acting on the wavefield
$\W$ with the Fourier domain operator dual.  This yields the
following relationship, 
\begin{equation}
\label{eqn:disp}
i n_j \kj -m^{ij} \ki\kj = - \gdetf \ww^2\ss^2\end{equation}
where $\ki$ is the Fourier domain dual of $\frac{\partial}{\partial \xi_i}$.
%
\subsection{One-way wavefield extrapolation}
%
In the following, coordinate $\xi_3$ is assumed to be the
extrapolation is direction.  Hence, obtaining a extrapolation
direction wavenumber requires isolating the contribution of $\kz$.
This is done by expanding out 
equation~\ref{eqn:disp} to yield,
\beqa \label{eqn:disp2}
m^{11}\kx^2+m^{22}\ky^2+m^{33}\kz^2+2m^{12}\kx\ky+2m^{13}\kx\kz+2m^{23}\ky\kz
\nonumber \\
=  \gdetf \ww^2 \ss^2 + i \lp n_1\kx+n_2\ky+n_3\kz \rp.
\eeqa
Reordering the terms in equation~\ref{eqn:disp2} yields,
\beqa \label{eqn:disp2}
m^{33}\kz^2 + \lp 2m^{13}\kx + 2m^{23}\ky - i\,n_3 \rp \kz =
\nonumber \\
\gdetf \ww^2 \ss^2 + i \lp n_1\kx + n_2\ky \rp - m^{11}\kx^2 -
m^{22}\ky^2 - 2m^{12}\kx\ky 
\eeqa
The expression for wavenumber $\kz$ is isolated by completing
the square. 
%(i.e. use the fact that, 
%\beq
%ax^2+bx+c = a \lp x+\frac{b}{2a}\rp^2 + \lp c-\frac{b^2}{4a} \rp
%\eeq
%where, in this case, 
%%% = m^{33}, \;\;\;\;\;\; 
%b = (2m^{13}\kx+2m^{23}\ky - i\, n_3), \;\;\;\;\;\;\; 
%%c = 0.
%\eeq%
%Hence, 
%\beqa \label{eqn:cts}
%m^{33}\kz^2 + \lp 2m^{13}\kx + 2m^{23}\ky - i \,n_3 \rp \kz = \\
%\, 
%m^{33} 
%\lp 
%\kz + 
%\frac{ \lp 2m^{13}\kx + 2m^{23}\ky- i\,n_3 \rp }{ 2m^{33} } 
%\rp^2
%- \frac{ \lp 2m^{13}\kx+2m^{23}\ky- i\,n_3 \rp^2 }{4m^{33}} \nonumber
%\eeqa
Completion of this procedure yields, 
\beqa \label{eqn:disp3}
m^{33} 
\lp 
\kz + \frac{ \lp 2m^{13}\kx + 2m^{23}\ky-i \,n_3 \rp }{ 2m^{33} }  \rp^2
= & \, \\
\gdetf\ww^2\ss^2+
i \, \lp n_1\kx+n_2\ky\rp-
m^{11}\kx^2-
m^{22}\ky^2-
\, &
2m^{12}\kx\ky+
\frac{ \lp 2m^{13}\kx+2m^{23}\ky-i \, n_3 \rp^2 }{4m^{33}} \nonumber
\eeqa
Expanding out the last term in equation~\ref{eqn:disp3} leads to the
following arrangement of terms, 
\beqa
m^{33} \lp \kz + 
\frac{ \lp 2 m^{13}\kx + 2m^{23}\ky- i \,n_3 \rp }{ 2m^{33} } 
\rp^2 =  \, \\
\gdetf\ww^2\ss^2 -
\kx^2 \lp m^{11}-\frac{ \lp m^{13} \rp^2 }{ m^{33}} \rp -
\ky^2 \lp m^{22}-\frac{ \lp m^{23} \rp^2 }{ m^{33}} \rp -  \, \nonumber \\
\kx\ky \lp 2 m^{12}-\frac{ 2 m^{13}m^{23}}{ m^{33}} \rp +
i \, \kx \lp n_1 - \frac{m^{13}\,n_3}{ m^{33} } \rp + i \, \ky \lp n_2
- \frac{m^{23}\,n_3}{ m^{33} }\rp - \frac{n_3^2 }{ m^{33}} 
\,\nonumber
\eeqa
Finally, solving for wavenumber $\kz$ yields,
\begin{equation}
\label{eqn:kzfinal}
\kz =
- a_1 \kx
- a_2 \ky 
+ i a_3  
\pm
\lb
  a_4^2 \ww^2 
- a_5^2 \kx^2 
- a_6^2 \ky^2
- a_7 \kx\ky 
+ a_8 i\,\kx 
+ a_9 i\,\ky  
- a_{10}^2  
\rb^{\frac{1}{2}}\end{equation}
where non-stationary coefficients $a_i$ are given by,
\begin{equation}
\label{eqn:ametric2}
\left[
\begin{array}{c}
a_1 \\
a_2 \\
a_3 \\
a_4 \\
a_5 \\
a_6 \\
a_7 \\
a_8 \\
a_9 \\
a_{10}
\end{array}
\right] 
=
\left[
\begin{array}{c}
g^{13}/ g^{33} \\
g^{23}/ g^{33} \\
n_3 / 2 m^{33} \\
\ss / \sqrt{g^{33}} \\
\sqrt{g^{11}/ g^{33} - \lp g^{13}/g^{33}  \rp^2 } \\
\sqrt{g^{22}/ g^{33} - \lp g^{23}/g^{33}  \rp^2 } \\
2\,g^{12}/g^{33} -2\, g^{13}\,g^{23}/ \lp g^{33} \rp^2 \\
n_1/m^{33} - m^{13}\, n_3/ (m^{33})^2 \\
n_2/m^{33} - m^{23}\, n_3/ (m^{33})^2 \\
n_3 / m^{33}
\end{array}
\right]\;.\end{equation}
Note that there is a mixture of $m^{ij}$ and $g^{ij}$ terms, and that
coefficients $a_4, a_5, a_6$ and $a_{10}$ are squared.  This latter
point is required if the familiar expression of the split-step Fourier
correction is to be recovered.  Importantly, even though there are
more coefficients in the case of non-orthogonality, concentration can
be placed on making coefficients $a_i$ as smooth as possible, which 
should require fewer sets of reference parameters.  For example, the
coefficients for a weakly non-orthogonal coordinate system
(i.e. $\{g^{12},g^{13},g^{23}\} << \{g^{11},g^{22},g^{33} \}$) within a
kinematic approximation reduce to,
\begin{equation}
\label{eqn:ametric3}
\left[
\begin{array}{c}
a_1 \\
a_2 \\
a_3 \\
a_4 \\
a_5 \\
a_6 \\
a_7 \\
a_8 \\
a_9 \\
a_{10}
\end{array}
\right] 
\approx
\left[
\begin{array}{c}
0 \\
0 \\
0 \\
\ss / \sqrt{g^{33}} \\
\sqrt{g^{11}/ g^{33} } \\
\sqrt{g^{22}/ g^{33}} \\
0 \\
0 \\
0 \\
n_3 / m^{33}
\end{array}
\right]\;.\end{equation}
%
\subsection{Split-Step Fourier Approximation}
%
The extrapolation wavenumber in equation~\ref{eqn:kzfinal} cannot be
implemented as given because there are mixed domain terms (i.e. a
function of $\x$ and $\kx$ simultaneously).  The split-step Fourier
approach \cite{GEO55-04-04100421} uses Taylor expansions to separate this
wavenumber into two parts: a phase-shift, $\kz^{PS}$, applied in the
pure Fourier domain, and a split-step correction, $\kz^{SSF}$, applied
in the mixed $\omega-\mathbf{x}$ domain (i.e. $\kz 
\approx \kz^{PS} + \kz^{SSF}$).  These phase-shift term is given by
the following relationships, 
\begin{equation}
\kz^{PS} = 
- b_1 \kx
- b_2 \ky 
+ i b_3  
\pm
\lb
  b_4^2 \ww^2 
- b_5^2 \kx^2 
- b_6^2 \ky^2
- b_7 \kx\ky 
+ b_8 i\,\kx 
+ b_9 i\,\ky  
- b_{10}^2  
\rb^{\frac{1}{2}}\end{equation}
where $b_i=b_i(\xi_3)$ are reference values of $a_i = a_i(\xi_1,\xi_2,\xi_3)$.
The split-step approximation is developed by performing a Taylor
expansion about each coefficient $a_i$ and then evaluating the results
at stationary reference values (i.e. $b_i$).  The stationary values of
$\kx$ and $\ky$ are assumed to be zero.  This leads to the following
split-step correction, 
\begin{equation}
\label{eqn:kzssf}
\kz^{SSF}=
\left. \frac{\partial \kz}{\partial a_3} \right|_0 \lp a_3-b_3 \rp +
\left. \frac{\partial \kz}{\partial a_4} \right|_0 \lp a_4-b_4 \rp + 
\left. \frac{\partial \kz}{\partial a_{10}} \right|_0 \lp a_{10}-b_{10} \rp\end{equation}
where ``$0$'' denotes with respect to a reference medium.  The partial
differential expressions in equation~\ref{eqn:kzssf} are, 
\begin{equation}
\left. \frac{\partial \kz}{\partial a_3} \right|_0 
= b_3,
\;\;\;\;\;
\left. \frac{\partial \kz}{\partial a_4} \right|_0
= \frac{b_4 \,\ww^2}{\sqrt{b_4^2 \, \ww^2 - b_{10}^2 }},
\;\;\;\;\;
\left. \frac{\partial \kz}{\partial a_{10}} \right|_0
= - \frac{b_{10}}{\sqrt{b_{10}^2 \, \ww^2 - b_{10}^2 }},
\;\;\;\;\;\end{equation}
and the total split-step Fourier correction is, 
\begin{equation}
\kz^{SSF} = i\,b_3 \, \lp a_3-b_3 \rp +
\frac{b_4 \, \ww^2 \, \lp a_4-b_4 \rp}{\sqrt{b_4^2 \, \ww^2 - b_{10}^2 }}  -
\frac{b_{10} \, \lp a_{10}-b_{10} \rp}{\sqrt{b_4\,\ww^2 - b_{10}^2 }}.\end{equation}
Note that we could use many reference media followed by interpolation
similar to the phase-shift plus interpolation (PSPI) technique of
\longcite{Gazdag.1984}. 
%
\section{Example 1 - 2D Sheared Cartesian Coordinates}
%
An interesting non-orthogonal coordinate system example is one formed
by a 'shearing' action on a Cartesian mesh.  This coordinate system is
illustrated in figure~\ref{fig:2Dexamp}.  
\activeplot{2Dexamp}{width=6in,height=2.5in}{NR}{2D Sheared Cartesian
  Coordinates.  Left panel: Physical domain represented by sheared
  Cartesian coordinates defined by ${x_1,x_3}$; Right panel:
  Canonical domain represented by a unit square where coordinate 
  $\xi_3$ is the extrapolation direction and $\xi_1$ is the orthogonal
  direction.} 
A sheared Cartesian coordinate system is defined by the following
equations,   
\begin{equation}
\label{eqn:2Dexamp}
  \left[\begin{array}{c}
  	\xone \\
  	\xtwo \\
  	\xthr
  	\end{array}\right] 
  	= 
  \left[\begin{array}{ccc}
      1 & 0 & \rm{cos} \theta \\
      0 & 0 & 0      \\
      0 & 0 & \rm{sin} \theta
  \end{array}\right]
  \left[\begin{array}{c}
    \qx \\
    \qy \\
    \qz \\
  \end{array}\right]\end{equation}
where $\theta$ is defined as the dip angle of the coordinate system
(i.e. $\theta=90^{\circ}$ is Cartesian).  This system can be reduced
to a more workable 2D set of equations.
\begin{equation}
\label{eqn:2Dred}
  \left[\begin{array}{c}
  	\xone \\
  	\xthr
  	\end{array}\right] 
  	= 
  \left[\begin{array}{cc}
      1 & \rm{cos} \theta      \\
      0 & \rm{sin} \theta
  \end{array}\right]
  \left[\begin{array}{c}
    \qx \\
    \qz \\
  \end{array}\right]\end{equation}
This leads to a metric tensor $g_{ij}$ that has the following form,
\begin{equation}
\label{eqn:2Dmet}
\left[g_{ij}\right]
=
\left[\begin{array}{cc}
    \EE & \HH \\
    \HH & \AA
\end{array}\right] 
= 
\left[\begin{array}{cc}
    1                & \rm{cos}\theta \\
    \rm{cos}\,\theta & 1  
\end{array}\right]\end{equation}
The determinant of $\mathbf{g}$ is given by,
\begin{equation}
\gdet = \rm{sin}^2\,\theta\end{equation}
The associated metric tensor (i.e. the inverse of $g_{ij}$) is given
by,
\begin{equation}
\label{eqn:2Dmet2}
\left[g^{ij}\right]
=
\left[\begin{array}{cc}
    \frac{1}{\rm{sin}^2 \theta}  &
    - \frac{ \rm{cos}\,\theta}{\rm{sin}^2\,\theta }\\
    - \frac{ \rm{cos}\,\theta}{\rm{sin}^2\,\theta } & 
    \frac{1}{\rm{sin}^2 \theta}  
\end{array}\right]\end{equation}
Note that equation~\ref{eqn:2Dmet2} is independent of all
coordinates. Hence, in this special case we are able to use a
simplified version of equation~\ref{eqn:laplac}.
\begin{equation}
\label{eqn:laplac}
  \Delta \W = g^{ij} \,\frac{\partial^2 \W}{\partial \qi \partial \qj}
  = - \ww^2\ss^2\W\end{equation}
which allows for the following dispersion relation,
\begin{equation}
g^{ij}\ki\kj = \ww^2\ss^2\end{equation}
Expanding out all of the terms leads to an equation similar to the
above 2-D nonorthogonal coordinate systems discussed above 
\begin{equation}
\label{eqn:kz2Dex}
\kz =
-\frac{ g^{13} }{ g^{33} } \kx
\pm
\sqrt{\frac{\ss^2\ww^2}{g^{33}}- \lp\frac{g^{11}}{g^{33}}-\lp
  \frac{g^{13}}{g^{33}} \rp^2\rp \kx^2}. \end{equation}
Substituting the values of the associated metric tensor in
equation~\ref{eqn:2Dmet2} into equation~\ref{eqn:kz2Dex} yields, 
\begin{equation}
\kz =
\rm{cos} \, \theta \;\kx
\pm \rm{sin} \,\theta
\sqrt{\ss^2\ww^2 - \kx^2}.  \end{equation}
%
\par
%
A numerical test was used on a Cartesian coordinate system sheared at
25$^{\circ}$ from the vertical axis.  The coordinate system is shown
in figure~\ref{fig:Rays0}.  The background velocity was set at
1500~ms$^{-1}$.  Zero-offset data consist of 4 planar impulses with
zero dip at t=0.2, 0.4, 0.6 and 0.8 s.  Zero-offset migration should
image reflectors at z=300, 600, 900, and 1200 m in depth.
%
\activesideplot{Rays0}{width=4in}{ER}{Sheared Cartesian coordinate system with
  a shear angle of 25$^{\circ}$.  Background velocity is 1500
  ms$^{-1}$.  Zero-offset data consist of 4 planar impulses at t=0.2,
  0.4, 0.6 and 0.8 s.  Zero-offset migration should place
  reflectors at depths z=300, 600, 900, and 1200 m.}   
%
The zero-offset migration section is presented in
figure~\ref{fig:CC0test}.  As expected in a constant velocity medium,
the zero-offset migration has imaged reflectors at depths z=300, 600,
900, and 1200 m. 
\activeplot{CC0test}{width=5in}{ER}{Zero-offset migration in a constant velocity
  medium (1500~ms$^{-1}$).  As expected in a constant velocity medium,
the zero-offset migration has imaged reflectors at depths z=300, 600,
900, and 1200 m.}
Note that the limits of the migrated area are visible, indicating that
the data have been successfully translated to depth.
%
\section{Example 2 - Polar Ellipsoidal}
%
A second interesting non-orthogonal coordinate system is one formed by
stretching a polar coordinate system.  For example, a wavefield may be
extrapolated in increments of $\Delta \theta$ in the $\xi_3$
direction.  The other axis is formed in the radial direction $\xi_1$.
The polar ellipsoidal coordinate system is illustrated in
figure~\ref{fig:2Dex2}. 
%
\activeplot{2Dex2}{width=5in}{NR}{Illustration of the polar ellipsoidal
coordinate system} 
%
A polar ellipsoidal coordinate system can formed by the following equations:
\begin{equation}
\label{eqn:2DPE1}
  \left[\begin{array}{c}
  	\xone \\
  	\xthr
  	\end{array}\right] 
  	= 
  \left[\begin{array}{c}
    a(\xi_3)\,\xi_1\, \rm{cos}\,\xi_3 \\
    a(\xi_3)\,\xi_1\, \rm{sin}\,\xi_3 \\
  \end{array}\right]\end{equation}
Here parameter $a=a(\xi_3)$ is a smooth function that controls the
ellipticity of the coordinate system.  The following usesSchematic of a pure polar coordinate
  system with $a=1$.
$b=\frac{\partial a}{\partial \xi_3}$ and $c=\frac{\partial^2
  a}{\partial \xi_3^2}$ for notational simplicity.  Metric tensor
$g_{ij}$ has the following form, 
\begin{equation}
\label{eqn:2DPE2}
\left[g_{ij}\right]
=
\left[\begin{array}{cc}
    a^2 & \xi_1\, a\, b \\
    \xi_1 \, a \, b & \xi_1^2(b^2+a^2) 
\end{array}\right]\end{equation}
The determinant of $\mathbf{g}$ is given by,
\begin{equation}
\gdet = a^4 \xi_1^2\end{equation}
The associated metric tensor (i.e. the inverse of $g_{ij}$) is given
by,
\begin{equation}
\label{eqn:2DPE3}
\left[g^{ij}\right]
=
\left[\begin{array}{cc}
    \frac{b^2+a^2}{a^4}   & -\frac{b}{a^3\,\xi_1} \\
    -\frac{b}{a^3\,\xi_1} & \frac{1}{a^2\, \xi_1^2}
\end{array}\right]\end{equation}
Matrix $m^{ij}=\gdetf \, g^{ij}$ is given by,
\begin{equation}
\label{eqn:2DPE4}
\left[g^{ij}\right]
=
\left[\begin{array}{cc}

    \frac{\xi_1(b^2+a^2)}{a^2} & -\frac{b}{a} \\
    -\frac{b}{a} & \frac{1}{\xi_1}
\end{array}\right]\end{equation}
We may use tensor $m^{ij}$ in the 2D non-orthogonal coordinate system
in equation~\ref{eqn:2DNO}.  In this case, we must compute $n_1$ and
$n_3$ since the coordinate system varies as a function of space.
These are given by expressions $n_1=\frac{a^2+2b^2-ac}{a^2}$ and
$n_3=0$, respectively.  Inserting the values for $m^{ij}$ and $n_j$
yields the following extrapolation wavenumber $\kz$, 
\begin{equation}
\label{eqn:2Dex5}
\kz = \frac{\xi_1 b}{a} \kx \pm \sqrt{ a^2 \xi_1^2 s^2 \omega^2 -
  \xi^2_1 \kx^2 - i \kx \xi_1 \frac{a^2+2b^2-ac}{a^2} } \end{equation}
A kinematic version of equation~\ref{eqn:2Dex5} is,
\begin{equation}
\label{eqn:2Dex6}
\kz = \xi_1 \lb \frac{b}{a} \kx \pm \sqrt{ a^2 s^2 \omega^2 -
  \kx^2 }\rb \end{equation}
Note that for the orthogonal polar case where $a=1$, we recover the
following, 
\begin{equation}
\label{eqn:2Dex7}
\kz = \pm \xi_1 \sqrt{ s^2 \omega^2 -\kx^2 } \end{equation}
%
\par
%
The ellipsoidal polar coordinate system is illustrated in the next few
figures.  Figure~\ref{fig:Rays1} presents a coordinate system with
ellipticity factor $a=1$ (i.e. pure polar coordinates).
%
\activesideplot{Rays1}{width=3in}{ER}{Schematic of a pure polar coordinate
  system with $a=1$.}
%

%
\activesideplot{Vel1}{width=4in,height=3in}{ER}{Velocity model projected into
  Riemannian coordinates.  The velocity model in Cartesian
  coordinates is a $v=1500+0.2z ms^{-1}$.}
%

%
\activesideplot{Data1}{width=4in,height=3in}{ER}{Synthetic zero-offset data.}
%

%
\activeplot{CC1test}{width=7in}{ER}{Zero-offset migration results.  As expected
  the plane-waves with zero-dip propagate downwards and do not begin
  to overturn.}
%
\newpage
\bibliographystyle{sep}
\bibliography{SEP,MISC,GEOTLE,refs,refs2}

\newpage

\APPENDIX{A}
%
The extrapolation wavenumber developed in equation~\ref{eqn:kzfinal}
is appropriate for any non-orthogonal Riemannian geometry.  However,
there are a number of situations where symmetry or partial
orthogonality are present in coordinate systems.  Moreover, one may
also wish to make a kinematic approximation where all of the the
imaginary components of the wavenumber are ignored.  These situations
are discussed in the following section.
%
\subsection{3D Semi-orthogonal Coordinate Systems}
%
Semi-orthogonal coordinate systems occur  where one coordinate is orthogonal
to the other two coordinates \longcite{Sava.geo.rwe}.  In these cases
the $m^{13}$ and $m^{23}$ components of the weighted metric tensor are
identically zero, which leads to the following extrapolation
wavenumber,  
\begin{equation}
\kz = \frac{\rm{i} n_3}{2 m^{33}}
\pm   
\lb
  a_0 \ww^2 
- a_1 \kx^2 
- a_2 \ky^2
- a_3 \kx\ky 
+ a_4 \rm{i}\kx 
+ a_5 \rm{i}\ky  
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where, $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$, $a_1 = \frac{m^{11}}{
  m^{33} }$, $a_2 = \frac{m^{22} }{ m^{33}}$,  
$a_3 = \frac{ 2 m^{12} }{ m^{33} }$, $a_4 = \frac{n_1}{m^{33}}$, $a_5
= \frac{n_2}{m^{33}}$, and $a_6 = \lp \frac{n_3}{m^{33}} \rp^2$, which
are the coefficients recovered by \longcite{Sava.geo.rwe}.
%
\subsection{3-D Kinematic Coordinate Systems}
%
Wave-equation migration amplitudes are generally incorrect in
laterally variant media - even in a Cartesian based system.  Hence,
one beneficial approximation that reduces computational cost
is to consider only the kinematic terms in equation~\ref{eqn:kzfinal}.
This approximate generates the following extrapolation wavenumber,  
\begin{equation}
\kz =
-\frac{ m^{13} }{ m^{33} }\kx
-\frac{ m^{23} }{ m^{33} }\ky 
\pm
\lb
  a_0 \ww^2 
- a_1 \kx^2 
- a_2 \ky^2
- a_3 \kx\ky 
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where, $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$,$a_1 = \frac{m^{11} }{
  m^{33} } - \lp \frac{ m^{13} }{ m^{33} } \rp^2$, $a_2 = \frac{m^{22}
}{ m^{33} } - \lp \frac{ m^{23} }{ m^{33} }\rp^2$, $a_3 = \frac{ 2
  m^{12} }{ m^{33} } -\frac{ 2 m^{13} m^{23}} {\lp m^{33}\rp^2 }$, and
$a_6 = \lp \frac{n_3}{m^{33}}\rp^2$.
%
\subsection{3-D Kinematic Semi-orthogonal coordinate systems}
%
Combining the two above restrictions yields the following
extrapolation wavenumber, 
\begin{equation}
\kz =
\pm
\lb
  a_0 \ww^2 
- a_1 \kx^2 
- a_2 \ky^2
- a_3 \kx\ky 
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where, $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$, $a_1=\frac{m^{11} }
{m^{33}}$, $a_2 = \frac{   m^{22} }{ m^{33} }$, $a_3 = \frac{ 2 m^{12}
}{ m^{33} }$, and $a_6 = \lp \frac{n_3}{m^{33}} \rp^2$.  Note that the 
expression $\PP=1$, and that components of the metric tensor are
significantly simplified.  
% 
\subsection{2-D Non-orthogonal coordinate systems}
%
Two-dimensional situations are handled by identifying $\qy=0$.  Hence,
all derivatives in the associated metric tensor $g^{ij}$ with respect
coordinate $\qy$ are identically zero. Hence, the 2-D non-orthogonal
coordinate system can be represented by  
\begin{equation}
\label{eqn:2DNO}
\kz =
-\frac{ m^{13} }{ m^{33} }\kx
+\frac{ \rm{i} n_3 }{ 2 m^{33} } 
\pm
\lb
  a_0 \ww^2 
- a_1 \kx^2 
+ a_4 \rm{i}\kx 
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$, $a_1 = \frac{m^{11}
}{ m^{33}}-\lp \frac{ m^{13} }{ m^{33} } \rp^2$, $a_4 =
\frac{n_1}{m^{33}}-\frac{m^{13} n_3}{(m^{33})^2}$, and $a_6 =
\lp \frac{n_3}{m^{33}} \rp^2$.
%
\subsection{2-D Non-orthogonal Kinematic Coordinate Systems}
%
Two-dimensional kinematic situations are handled through identity
$\qy=0$.  Again, all derivatives in the associated metric tensor
$g^{ij}$ with respect coordinate $\qy$ are identically zero, and
the 2-D non-orthogonal kinematic extrapolation wavenumber is 
\begin{equation}
\kz =
-\frac{ m^{13} }{ m^{33} }\kx
\pm
\lb
  a_0 \ww^2 
- a_1 \kx^2 
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$, $a_1 = \frac{   m^{11}
}{ m^{33} } - \lp \frac{ m^{13} }{ m^{33} } \rp^2$, and $a_6 =
\lp \frac{n_3}{m^{33}} \rp^2$. 
%
\subsection{2-D Orthogonal Coordinate Systems}
%
Two-dimensional situations are handled through the following identity
$\qy=$Constant and $\HH=0$.  Hence, all derivatives in the associated
metric tensor $g^{ij}$ with respect coordinate $\qy$ are identically
zero. Hence, the 2-D non-orthogonal coordinate system can be
represented by 
\begin{equation}
\kz =
\frac{ \rm{i} n_3 }{ 2 m^{33} } 
\pm
\lb
  a_0 \ww^2 
- a_1 \kx^2 
+ a_4 \rm{i}\kx 
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where, $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$, $a_1 =
\frac{m^{11}}{m^{33}}$, $a_4=\frac{n_1}{m^{33}}$, and $a_6 = \lp
\frac{n_3}{m^{33}} \rp^2$. 
%
\subsection{2-D Orthogonal Kinematic Coordinate Systems}
%
The above two approximations can be combined to yield the following
extrapolation wavenumber for 2D orthogonal kinematic coordinate
systems, 
\begin{equation}
\kz =
\pm
\lb
  a_0 \ww^2 
- a_1 \kx^2 
- a_6  
\rb^{\frac{1}{2}}\end{equation}
where $a_0 = \frac{ \gdetf \ss^2 }{ m^{33} }$,
$a_1=\frac{m^{11}}{m^{33} }$ and $a_6=\lp \frac{n_3}{m^{33}} \rp^2$.
%
