\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$. %