Riemannian Wavefield Extrapolation: Non-orthogonal coordinate systems

Jeff Shragge Stanford Exploration Project, Mitchell Bldg., Department of Geophysics, Stanford University, Stanford, CA 94305-2215 (October 11, 2006)

ABSTRACT

Riemannian wavefield extrapolation (RWE) is a method for performing one-way wave propagation on generalized coordinate meshes. Previous RWE implementations assume that coordinate systems are either orthogonal or semi-orthogonal, which can lead to situations where meshes suffer from problematic bunching, singularities, or multivaluedness. This paper develops a procedure for avoiding these problems that is based on wavefield extrapolation on smoother, but generally non-orthogonal, coordinate system meshes. The corresponding extrapolation wavenumber is comprised of a number of mixed-domain fields, in addition to velocity, that encode coordinate system geometry. An extended split-step Fourier approximation of the extrapolation operator is developed that provides accurate results when multiple reference parameter sets are used. Two analytic coordinate system examples are presented to validate the theory. A procedure for generating triplication-free coordinate systems in complex media is developed, which is followed by calculations of 2D and 3D Green’s functions estimates in cylindrical and near-spherical geometries. An important consequence is that greater emphasis can be placed on generating smooth computational meshes for RWE rather than satisfying more restrictive semi-orthogonal criteria.

INTRODUCTION

A persistent goal of wave-equation migration research is to improve seismic imaging capabilities in complex geologic settings. Although ubiquitous velocity model uncertainty and uneven illumination can contribute greatly to image interpretation ambiguity in these contexts, extrapolation operator inaccuracy remains a significant problem. The central issues with conventional one-way wave-equation extrapolation operators are well documented: while naturally handling wavefield multipathing in the presence of lateral velocity variation, they are of limited accuracy at steep propagation angles and cannot propagate overturning waves by design (though extensions including phase-shift table lookup (Hale et al., 1992) or two-pass migration (Zhang et al., 2005) obviate some of these concerns). Propagation errors are subsequently manifest in migration images as defocused or misplaced reflectors or even by a complete absence of interpretable reflectivity. Accordingly, minimizing these deleterious effects should improve image quality and any subsequent interpretation based thereon.

One strategy for reducing extrapolation operator inaccuracy is to decompose the complete computational grid into subdomains oriented in the wave propagation direction (see, e.g., Gaussian beams (Hill, 2001), beam-waves (Brandsberg-Dahl and Etgen, 2003), coherent states (Albertin et al., 2001) or tilted Cartesian meshes (Etgen, 2002)). The key concept in each of these approaches is that a judicious choice of reference frame lowers the effective propagation angle, reducing the need for expensive extrapolation operators and enabling imaging of overturned waves. Sava and Fomel (2005) followed this approach in developing Riemannian wavefield extrapolation (RWE), a theory of one-way wavefield propagation through a generalized Riemannian space. This formulation specifies waveequation operators appropriate for wavefield extrapolation on generalized computational meshes. One important ramification is that the user is free to specify the degree to which wave-propagation effects are incorporated into the computational mesh. However, finding the optimal trade-off between computational mesh simplicity, how well the mesh conforms to the wavefield propagation direction, and the computational cost is not a straightforward task.

RWE was initially implemented to model high-quality Green’s functions. This process involved extrapolating wavefields on a point-source coordinate mesh comprised of a suite of rays traced beforehand through a smoothed version of the migration velocity model. Hence, RWE computational meshes explicitly were asserted to exhibit ray-field characteristics: semi-orthogonal geometry with an extrapolation direction (i.e. travel-time along a ray) orthogonal to the two other axes (i.e. shooting angles) that are not necessarily mutually orthogonal. This supposition results in a wave-equation dispersion relationship for semi-orthogonal meshes that contains a number of mixed-domain fields, in addition to velocity, that encoded coordinate system geometry. In most examples, modeled Green’s function estimates interpolated into the Cartesian domain are highly accurate at steep propagation angles; however, accuracy is compromised in certain situations exhibiting certain unfavorable mesh characteristics.

In general, the semi-orthogonality assertion is overly restrictive because it precludes using nonorthogonal computational meshes designed to be smoothly varying, with singularities and singlevalued. One problematic example occurs in coordinate system triplications (see upper panel of figure 1), that arise wherever a mesh is developed from a triplicating field of rays. This generates spatial singularities that lead to zero-division of the wavefield during extrapolation and cause infinite amplitudes. Although ray-coordinate triplications can be avoided by iterative velocity model smoothing, this less-than-ideal solution counters the goal of having a coordinate system conformal to the wavefield propagation direction. A second example of restrictive semi-orthogonal geometry is illustrated in Shragge and Sava (2005), who formulate a wave-equation migration from topography strategy that poses wavefield extrapolation directly in locally orthogonal meshes conformal to the acquisition surface. This approach successfully generates subsurface images beneath areas exhibiting longer wavelength and lower amplitude relief; however, imaging results in situations involving more rugged acquisition topography degrade due to the grid compression/extension demanded by semi-orthogonality geometry (see lower panel of figure 1b).

In this paper, I argue that the RWE approach can be made more broadly applicable by resolving the issues related to semi-orthogonal geometry using a combination of localized computational mesh smoothing and the judicious removal of all coordinate triplication branches. In doing so, though, the asserted semi-orthogonality mesh geometry is lost, and one must propagate wavefields on nonorthogonal grids. Fortunately, a non-orthogonal extrapolation wavenumber can be isolated from the general RWE dispersion relationship. This wavenumber is used in forming a basis for a one-way extrapolation operator that is used in the usual sense of wave-equation imaging.

The goal of this paper is develop and implement a one-way wave-equation extrapolation operator appropriate for RWE in 3D non-orthogonal coordinates. A second goal is to specify a procedure for generating unconditionally triplication-free computational meshes. The development presented herein naturally follows that of Sava and Fomel (2005); however, I recast the theory in a more compact notation that leads to a closer analytic connection of the generalized computation geometry with the underlying Cartesian grid. The paper begins with the formulation of the 3D Riemannian acoustic wave-equation and the corresponding non-orthogonal one-way wavefield extrapolation wavenumber. Appendix A presents an overview of the required differential geometry theory, while the split-step Fourier approximation corresponding to the extrapolation operator developed herein is presented in Appendix B. Two analytic 2D non-orthogonal coordinate system examples are then provided to validate the theory. The final sections detail the procedure for generating triplication-free coordinate systems, and present 2D and 3D Green’s functions estimates modeled in cylindrical and near-spherical coordinates, respectively.

ACOUSTIC WAVE-EQUATION IN 3D GENERALIZED RIEMANNIAN SPACES

Specifying the acoustic wave-equation in a 3D Riemannian space requires formulating the physics of wave-propagation in a generalized coordinate system framework. By definition, generalized Riemannian coordinates are related to the underlying Cartesian mesh by unique transformations (i.e. non-triplicating and one-to-one). In this paper, I use a notation where a generalized coordinate system ξ D {ξ_{1},ξ_{2},ξ_{3}} to mapped to a Cartesian grid x D {x_{1},x_{2},x_{3}} through transformation ξ_{i}(x_{j}) D f_{i}. Provided these conditions are met, the acoustic wave-equation for a wavefield, �, in a generalized Riemannian space is,

∇_{ξ}^{2}� D −ω^{2}s^{2}(ξ)�,  (1)

where ∇_{ξ}^{2} is the Laplacian operator applied in coordinates ξ, ω is frequency, and s is the propagation slowness.

Correctly formulating the wave-equation in coordinate system ξ requires that Laplacian operator ∇_{ξ}^{2} be specified by differential geometry relationships. (An overview of necessary differential geometry theory is provided in Appendix A.) The Laplacian operator in generalized coordinates is,

∇_{ξ}^{2}� D \frac{1}{\sqrt{|g|}} \frac{∂}{∂ξ_{i}} ( \sqrt{|g|}g^{ij} \frac{∂�}{∂ξ_{j}} ) D \frac{1}{\sqrt{|g|}} \frac{∂}{∂ξ_{i}} ( m^{ij} \frac{∂�}{∂ξ_{j}} ) , i,j D 1,2,3,  (2)

where g^{ij} is an element of the metric tensor g, |g| is the metric tensor discriminant, and m^{ij} D \sqrt{|g|}g^{ij} is weighted metric tensor element that enables a more compact notation. Unless otherwise stated, summation over all repeated indices (i.e. i,j D 1,2,3) is assumed throughout this paper. Substituting equation 2 into 1 generates a Helmholtz equation appropriate for propagating waves through a 3D Riemannian space,

\frac{1}{\sqrt{|g|}} \frac{∂}{∂ξ_{i}} ( m^{ij} \frac{∂�}{∂ξ_{j}} ) D −ω^{2}s^{2}�.  (3)

The first step in developing a generalized RWE wave-equation dispersion relationship is to expand the derivative terms in equation 3 and multiply through by \sqrt{|g|} to obtain,

\frac{∂m^{ij}}{∂ξ_{i}} \frac{∂�}{∂ξ_{j}} C m^{ij} \frac{∂^{2}�}{∂ξ_{i}∂ξ_{j}} D −\sqrt{|g|}ω^{2}s^{2}�.  (4)

The spatial derivative of the weighted metric tensor in the first term of equation 4 is written concisely using the following substitution,

n^{j} D \frac{∂m^{ij}}{∂ξ_{i}} D \frac{∂m^{1j}}{∂ξ_{1}} C \frac{∂m^{2j}}{∂ξ_{2}} C \frac{∂m^{3j}}{∂ξ_{3}} .  (5)

Fields n^{j} are interpreted as measures of the rates by which space expands, compresses and/or shears in the j^{th}direction and can be non-zero even for orthogonal coordinate systems. Using this substitution, equation 4 is rewritten,

n^{j} \frac{∂�}{∂ξ_{j}} C m^{ij} \frac{∂^{2}�}{∂ξ_{i}∂ξ_{j}} D −\sqrt{|g|}ω^{2}s^{2}�.  (6)

A wave-equation dispersion relation is developed by replacing the partial differential operators acting on wavefield � with their Fourier domain wavenumber duals (Claerbout, 1985),

(m^{ij} k_{ξ_{i}} − in^{j} ) k_{ξ_{j}} D \sqrt{|g|}ω^{2}s^{2},  (7)

where ik_{ξ_{i}} is the Fourier domain dual of differential operator \frac{∂}{∂ξ_{i}}.

Equation 7 represents the dispersion relationship required to propagate a wavefield through a generalized 3D Riemannian space. Quantity m^{ij} in the first term, m^{ij} k_{ξ_{i}}k_{ξ_{j}}, is a measure of the dot product between wavenumber vectors in the k_{ξ_{i}} and k_{ξ_{j}} directions (i.e. orthogonal wavenumbers will have a zero m^{ij} coefficient). Fields n^{j} in the second term, in^{j} k_{ξ_{j}}, represent a scaling of wavenumber k_{ξ_{j}} caused by local expansion, contraction and/or shearing of the coordinate system in the j^{th} direction.

Note that the expression in equation 7 reduces to the more familiar Cartesian expression when introducing n^{j} D 0 and m^{ij} D δ^{ij}, which leaves

k_{ξ_{i}}k_{ξ_{i}} D k_{ξ_{1}}^{2}C k_{ξ_{2}}^{2}C k_{ξ_{3}}^{2} D ω^{2}s^{2}.  (8)

Extrapolation wavenumber isolation

Specifying a one-way extrapolation operator requires isolating one of the wavenumbers in equation 7. In this paper, I associate the extrapolation direction with coordinate ξ_{3}. Expanding equation 7 by evaluating indices and rearranging the results yields,

m^{33}k_{ξ_{3}}^{2}C (2m^{13}k_{ξ_{1}}C 2m^{23}k_{ξ_{2}}− in^{3}) k_{ξ_{3}} D \sqrt{|g|}ω^{2}s^{2}C i(n^{1}k_{ξ_{1}}C n^{2}k_{ξ_{2}}) − m^{11}k_{ξ_{1}}^{2}− m^{22}k_{ξ_{2}}^{2}− 2m^{12}k_{ξ_{1}}k_{ξ_{2}}.(9)

An expression for wavenumber k_{ξ_{3}} is obtained through a complete-the-square transform,

m^{33}( k_{ξ_{3}}C^{\frac{(2m^{13}k_{ξ_{1}}C2m^{23}k_{ξ_{2}}−in^{3})}{2m^{33}}} ) ^{2} D \sqrt{|g|}ω^{2}s^{2}− k_{ξ_{1}}^{2}( m^{11}−^{\frac{(m^{13})^{2}}{m^{33}}}) −  (10)
k_{ξ_{2}}^{2}( m^{22}−^{\frac{(m^{23})^{2}}{m^{33}}}) − k_{ξ_{1}}k_{ξ_{2}}(2m^{12}− \frac{2m^{13}m^{23}}{m^{33}}) C ik_{ξ_{1}}(n^{1}− \frac{m^{13}n^{3}}{m^{33}}) C ik_{ξ_{2}}(n^{2}− \frac{m^{23}n^{3}}{m^{33}}) − \frac{(n^{3})^{2}}{m^{33}}.

Isolating the contributions of wavenumber k_{ξ_{3}} yields,

k_{ξ_{3}} D −a_{1}k_{ξ_{1}}− a_{2}k_{ξ_{2}}C ia_{3}± [a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}− a_{6}^{2}k_{ξ_{2}}^{2}− a_{7}k_{ξ_{1}}k_{ξ_{2}}C a_{8}ik_{ξ_{1}}C a_{9}ik_{ξ_{2}}− a_{10}^{2}] ^{\frac{1}{2}} ,  (11)

where a_{i} are non-stationary coefficients given by,

a D [ \frac{g^{13}}{g^{33}} \frac{g^{23}}{g^{33}} \frac{n^{3}}{2m^{33}} \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}} − ( \frac{g^{13}}{g^{33}} ) ^{2}} \sqrt{\frac{g^{22}}{g^{33}} − ( \frac{g^{23}}{g^{33}} ) ^{2}} [ \frac{2g^{12}}{g^{33}} − \frac{2g^{13}g^{23}}{(g^{33})^{2}} ] [ \frac{n^{1}}{m^{33}} − \frac{m^{13}n^{3}}{(m^{33})^{2}} ] [ \frac{n^{2}}{m^{33}} − \frac{m^{23}n^{3}}{(m^{33})^{2}} ] \frac{n^{3}}{m^{33}} ] ^{T} . (12)

Note that the coefficients contain mixed m^{ij} and g^{ij} terms and positive definite terms a_{4},a_{5},a_{6} and a_{10} are squared.

The special Cartesian case is again recovered from the two equations above by substitution n^{j} D 0 and g^{ij} D δ^{ij} (for all values i and j indicies) in the coefficients in equation 12

k_{ξ_{3}} D [s^{2}ω^{2}− k_{ξ_{1}}^{2}− k_{ξ_{2}}^{2}] ^{\frac{1}{2}} .  (13)

The dispersion relationship specified by equations 11 and 12 contains ten coefficients that represent mixed-domain fields. Similar to Cartesian-based wavefield extrapolation through homogeneous media, a constant-coefficient Fourier-domain (ω − k_{ξ}) phase-shift extrapolation scheme can be developed to recursively advance a wavefield from level ξ_{3} to level ξ_{3}C 1ξ_{3},

�(ξ_{3}C 1ξ_{3},k_{ξ_{1}},k_{ξ_{2}};ω) D �(ξ_{3},k_{ξ_{1}},k_{ξ_{2}};ω)e^{ik_{ξ_{3}}1ξ_{3}} .  (14)

If � represents a post-stack wavefield, an image � can be produced from the propagated wavefield through an imaging condition (Claerbout, 1985),

�(ξ_{3},ξ_{1},ξ_{2}) D ∑_{ω} �(ξ_{3},ξ_{1},ξ_{2};ω).  (15)

Situations where coefficients vary across an extrapolation step, though, required additional approximations. One straightforward approach is to use an extended split-step Fourier (SSF) method (Stoffa et al., 1990), which is detailed in Appendix B. This method uses Taylor expansions of the dispersion relation about a set of reference parameters to form a bulk phase-shift operator in the Fourier domain (ω − k_{ξ}). Differences between the reference and true parameters then form a correction term applied in the mixed ω − ξ domain. In addition, the phase-shift extrapolation scheme can be extended to incorporate multiple reference media followed by interpolation similar to the phase-shift plus interpolation (PSPI) technique of Gazdag and Sguazzero (1984).

Discussion

The accuracy of the extended SSF approach is directly related to the degree to which coefficients in equation 12 vary across a propagation step. At a first glance, one might expect that far too many expansions are required to make a PSPI approach practical. (For example, three reference expansions for each of the ten terms would seemingly require 3^{10} D 22599 separate wavefield extrapolations.)

However, three factors combine to greatly reduce the total number of required reference coefficient sets.

First, the a_{i} coefficients in equation 11 tend to be highly correlated because they are composed of similar metric tensor elements g^{ij}. Thus, the central issue is how accurately can we characterize vector coefficient fields. Fortunately, this problem is similar to the quantization problem in computer graphics: What is the fewest number of colors by which an image can be represented given a maximum allowable error? In this paper, I calculate reference coefficients using a multi-dimensional Lloyd’s algorithm developed by Tang and Clapp (2006). This iterative procedure allows for variable rate quantization that represents the multi-dimensional histogram of the coefficients with the sparsest number of points within a specified error tolerance.

Second, numerous situations exists where some coefficients are zero or negligible. For example, the mixed-domain fields for a 3D weakly non-orthogonal coordinate system (i.e. max|{g^{12},g^{13},g^{23}}| << max{g^{11},g^{22},g^{33}}) within a kinematic approximation (i.e. no imaginary terms) can be accurately approximated with only four coefficients. (Appendix C details situations where additional approximations are appropriate.) Third, one may apply algorithms that locally smooth the coordinate system mesh, which reduces the spatial variability of the coefficients and allows a more reliable representation of wavenumber k_{ξ_{3}}.

NUMERICAL MODELING EXAMPLES

2D sheared Cartesian coordinate system

An instructive analytic coordinate system to examine is a 2D sheared Cartesian grid formed by a uniform shearing action on a 2D Cartesian mesh (see left panel of figure 2). This coordinate system is uniquely specified by one additional degree of freedom and is related to an underlying Cartesian mesh through the following transformation,

[ \begin{matrix}x_{1} \\ x_{3} \end{matrix} ] D [ \begin{matrix}1 & cosθ \\ 0 & sinθ \end{matrix} ] [ \begin{matrix}ξ_{1} \\ ξ_{3} \end{matrix} ] ,  (16)

where θ is the shear angle of the coordinate system (θ D 90^{∘} is Cartesian). The metric tensor of this transformation is,

[g_{ij} ] D [ \begin{matrix}\frac{∂x_{k}}{∂ξ_{1}} \frac{∂x_{k}}{∂ξ_{1}} & \frac{∂x_{k}}{∂ξ_{1}} \frac{∂x_{k}}{∂ξ_{3}} \\ \frac{∂x_{k}}{∂ξ_{1}} \frac{∂x_{k}}{∂ξ_{3}} & \frac{∂x_{k}}{∂ξ_{3}} \frac{∂x_{k}}{∂ξ_{3}} \end{matrix} ] D [ \begin{matrix}g_{11} & g_{13} \\ g_{13} & g_{33} \end{matrix} ] D [ \begin{matrix}1 & cosθ \\ cosθ & 1 \end{matrix} ] ,  (17)

and has a discriminant |g| D sin^{2}θ and an associated metric tensor g^{ij} given by,

[g^{ij} ] D \frac{1}{sin^{2}θ} [ \begin{matrix}1 & −cosθ \\ −cosθ & 1 \end{matrix} ] .  (18)

Note that because the tensor in equation 18 is coordinate invariant, equation 6 simplifies to,

∇_{ξ}^{2}� D g^{ij} \frac{∂^{2}�}{∂ξ_{i}∂ξ_{j}} D −ω^{2}s^{2}�,  (19)

which generates the following dispersion relation,

g^{ij} k_{ξ_{i}}k_{ξ_{j}} D ω^{2}s^{2}.  (20)

Expanding out these terms leads to an expression for wavenumber k_{ξ_{3}},

k_{ξ_{3}} D −\frac{g^{13}}{g^{33}}k_{ξ_{1}}± \sqrt{\frac{s^{2}ω^{2}}{g^{33}} − ( \frac{g^{11}}{g^{33}} − ( \frac{g^{13}}{g^{33}}) ^{2}) k_{ξ_{1}}^{2}} .  (21)

Substituting the values of the associated metric tensor in equation 18 into equation 21 yields,

k_{ξ_{3}} D cosθ k_{ξ_{1}}± sinθ\sqrt{s^{2}ω^{2}− k_{ξ_{1}}^{2}} ,  (22)

which is appropriate for performing RWE on the 2D sheared Cartesian coordinate system shown in left panel in figure 2.

The right panel of figure 2 shows the results of extrapolating plane waves in a Cartesian coordinate system sheared at 25^{∘} from vertical. The background velocity model is 1500 m/s and the zero-offset data consist of 4 flat plane-waves at times t D 0.2,0.4,0.6 and 0.8 s. Zero-offset migration results generated by equation 15 show migrated reflectors at the expected depths of z=300, 600, 900, and 1200 m. Note that the propagation creates explainable boundary artifacts: reflections on the left are due to a truncated coordinate system and hyperbolic diffractions on the right are caused by truncated plane-waves.

Polar Ellipsoidal Coordinates

A second example is an ellipsoidal polar coordinate system (see upper left panel in figure 3) appropriate for migrating overturning plane-waves. A polar ellipsoidal coordinate system is specified by,

[ \begin{matrix}x_{1} \\ x_{3} \end{matrix} ] D [ \begin{matrix}a(ξ_{3})ξ_{1}cosξ_{3} \\ a(ξ_{3})ξ_{1}sinξ_{3} \end{matrix} ] ,  (23)

where parameter coordinate ξ_{1} is the radius from the center focus, ξ_{3} is polar angle, and is a D a(ξ_{3}) is a smooth function that controls coordinate system ellipticity and has curvature parameters b D \frac{∂a}{∂ξ_{3}} and c D \frac{∂^{2}a}{∂ξ_{3}^{2}}. The metric tensor g_{ij} for the polar ellipsoidal coordinate system defined in equation 23 is,

[g_{ij} ] D [ \begin{matrix}a^{2} & ξ_{1}a b \\ ξ_{1}a b & ξ_{1}^{2}(b^{2}C a^{2}) \end{matrix} ] ,  (24)

and has a metric discriminant given by |g| D a^{4}ξ_{1}^{2}. The associated and weighted associated metric tensors are given by,

[g^{ij} ] D [ \begin{matrix}\frac{b^{2}Ca^{2}}{a^{4}} & −\frac{b}{a^{3}ξ_{1}} \\ −\frac{b}{a^{3}ξ_{1}} & \frac{1}{a^{2}ξ_{1}^{2}} \end{matrix} ] and [m^{ij} ]D [ \begin{matrix}\frac{ξ_{1}(b^{2}Ca^{2})}{a^{2}} & −\frac{b}{a} \\ −\frac{b}{a} & \frac{1}{ξ_{1}} \end{matrix} ] .  (25)

Tensors g^{ij} and m^{ij} are used to form the extrapolation wavenumber appropriate for one-way wavefield propagation on a 2D polar ellipsoidal mesh. However, because the computational mesh is nonstationary, we must also compute the n^{i} fields: n^{1} D \frac{a^{2}C2b^{2}−ac}{a^{2}} and n^{3} D 0. Inserting these values leads to the following extrapolation wavenumber expression (see equation C-7),

k_{ξ_{3}} D \frac{ξ_{1}b}{a}k_{ξ_{1}}± \sqrt{a^{2}ξ_{1}^{2}s^{2}ω^{2}− ξ_{1}^{2}k_{ξ_{1}}^{2}− ik_{ξ_{1}}ξ_{1}( \frac{a^{2}C 2b^{2}− ac}{a^{2}} )} .  (26)

The kinematic approximation of equation 26 is given by,

k_{ξ_{3}} D ξ_{1}[ \frac{b}{a}k_{ξ_{1}}± \sqrt{a^{2}s^{2}ω^{2}− k_{ξ_{1}}^{2}} ] ,  (27)

and the orthogonal polar case (i.e. a D 1 and b D 0) yields,

k_{ξ_{3}} D ±ξ_{1}\sqrt{s^{2}ω^{2}− k_{ξ_{1}}^{2}} ,  (28)

which is examined in Nichols (1994).

Figure 3 shows a wavefield extrapolation example for an ellipsoidal polar coordinate system in equation 23 defined by ellipticity parameter a(ξ_{3}) D 1 C 0.2ξ_{3}− 0.05ξ_{3}^{2}. The upper and lower panels of figure 3 correspond to velocity/coordinate and wavefield domains, respectively. Similarly, the left and right panels represent the Cartesian and Riemannian domains. Note that wavefield interpolation between the latter two domains is possible because of the established mapping relationships. The upper left panel shows the polar ellipsoidal coordinate system mesh overlying a linear v(z) D 1500C 0.35z m/s velocity function. The upper right panel presents the velocity model mapped into the RWE domain under the transformations defined in equation 23. The data used in this test consist of ten plane-waves defined by a ray parameter p_{x} D −0.5 s/km. At greater depths the waves propagate faster and go through a turning point before propagating upward to the left. The tops of the waves, though, still travel through slower material and have not yet overturned. One interesting observation is that if a propagating plane-wave wavepath is well represented by a single ellipticity parameter a D a(ξ_{3}), then a polar elliptical mesh can form the coordinate basis for a plane-wave migration strategy that is potentially more dynamic than tilted Cartesian coordinates.

GENERATING TRIPLICATION-FREE COORDINATE MESHES

A RWE computational mesh design challenge is finding a coordinate system that is fairly conformal to the wavefield propagation direction, but is unconditionally triplication-free. This problem is illustrated by the panels in figure 4. The top panel shows a v(z) velocity model with three Gaussian anomaly inclusions overlain by a ray-coordinate system calculated by the Huygens rayfront tracing method (Sava and Fomel, 2001). Note that the anomalies cause both mesh triplications to the left and right of the model, as well as a grid rarefaction directly beneath the shot-point.

The center panel shows a first-arrival Eikonal equation solution for the same shot-point presented in the top panel. Plotted over-top are single-valued Eikonal solution isochrons. Note that isochron generally conform to the propagation direction and can be used to construct the extrapolation steps of a RWE computational mesh. The first mesh generation step is to extract initial and final isochron surfaces from the Eikonal equation solution to form the inner and outer mesh boundaries. The mesh domain is then enclosed by interpolating between the edges of the inner and outer bounding surfaces. The interior mesh can then be formed through bi-linear interpolation methods, such as blending functions (Liseikin, 2004; Shragge, 2006).

A corresponding triplication-free mesh is presented in the bottom panel of figure 4. The grid is regularly-spaced on the outer isochron and has a couple of dimples at the locations of removed triplications. These discontinuities have been reduced by applying a smoothing operator to the Eikonal equation solution before calculating the mesh. Importantly, smoothing does not greatly affect propagation accuracy because the coordinate system mesh forms only the skeleton on which wavefield extrapolation occurs.

2D Green’s Function Generation

The third test uses RWE to model 2D Green’s functions on coordinate systems constructed by the aforementioned approach. Figure 5 presents the slice of the SEG-EAGE salt velocity model chosen for the test. Importantly, the contrast between the salt body and sediment velocities causes complicated wavefield propagation including triplication and multipathing. The upper left panel shows the velocity model with an overlain coordinate system generated by meshing procedure discussed in the section above. The velocity model in the RWE domain is illustrated in the upper right panel.

The lower right panel shows the impulse response tests in the RWE domain. The impulses conform fairly well to the travel-time steps, except where they enter the salt body in the lower left of the image. The migration results mapped back to Cartesian space are shown in the lower left panel. The complicated wavefield left of the shot point advances through the salt body and is subsequently refracted upward. Note that wide-angle reflections off the top-salt/sediment boundary are also strongly present. Figure 6 presents a comparison test between RWE and Cartesian extrapolation. Beneath and to the right of the shot-point the wavefields are fairly similar except for a 90^{∘} phase-change associated with differences between modeling a Cartesian point-source and a RWE plane-wave. However, the wavefields are significantly different to the left of the shot-point because Cartesian-based extrapolation cannot propagate wavefields laterally with the same accuracy and upward at all. Hence, this energy is absent from the propagating wavefield in the lower panel.

3D Point Source Coordinates

The procedure discussed above can be extended to 3D point-source coordinate systems. However, any wavefield extrapolation procedure in near-spherical coordinates must account for the apparent singularity at the poles. One approach is to rotate the poles such that the small arc-length regions are at the Earth’s surface (Schneider, 1995). A second approach is to avoid sampling at the existing polar singularity (Fowler, 1995). Both solutions, though, have highly variable spatial sampling because regular spherical angle discritization leads to highly variable sampling in Cartesian space.

Another solution is to choose a discritization of the spherical surface that has a fairly uniform distribution in Cartesian domain. Herein, I use a 3D point source gridding technique based on the Winkel-Tripel cartographic projection (Bugayevskiy and Snyder, 1995). Note that he mapping of a spherical surface illustrated in figure 7 does not eliminate distortions in area, direction or distance; rather, it tries to minimize the sum of all three. The gridding equations for a Winkel-Tripel surface are the following,

x_{1} D ξ_{3}(2ξ_{1}/π C 2w αcosξ_{2}sin(ξ_{1}/2))/2
x_{2} D ξ_{3}(ξ_{2}C 2α sinξ_{2})/2  (29)
x_{3} D \sqrt{ξ_{3}^{2}− x_{1}^{2}− x_{2}^{2}}

where ξ_{3} is in the radial direction, ξ_{1} and ξ_{2} are a latitude parallel and longitude meridian, respectively, and,

α D cos^{−1} (cosξ_{2}cos(ξ_{1}/2))
w D \{ \begin{matrix}0,  sinα D 0 \\ (sinα)^{−1} , otherwise \end{matrix}  (30)

Note that the 3D grid is created by stepping outward in the radial direction ξ_{3} (or equally by the Eikonal equation solution method described above).

Figure 8 presents an example of a wavefield extrapolated on 3D point source coordinate mesh. The Elf-IFP-CGG synthetic velocity model used in this example is inspired from field data observations from North Sea Block L7d. The sediments located beneath the shot point form a mostly v(z) profile, while the salt body to the right generates a more complex geologic setting. Overlain are a number of impulse responses associated with different travel times. Note that the wavefield does not exhibit azimuthal anisotropy found in 3D impulse response tests based on finite-difference splitting methods. Also, the coordinate system boundaries are restricted at steep angles, which generates boundary reflection artifacts. Overall, though, this extrapolation results demonstrates the stability of wavefield propagation on a 3D grid defined by the Winkel-Tripel spatial discritization.

CONCLUDING REMARKS

This paper extends the implementation of Riemannian wavefield extrapolation theory to 3D nonorthogonal coordinate systems. The extrapolation wavenumber decouples from the other wavenumbers enabling the use an extended split-step Fourier approximate solution of the one-way wavefield extrapolation operator. Also documented is a methodology for generating computational meshes that are unconditionally triplication-free. The provided examples illustrate that wavefields can be extrapolated on non-orthogonal coordinate meshes, which opens up a range of possibilities including Riemannian plane-wave migration in elliptical polar coordinates.

ACKNOWLEDGMENTS

I acknowledge the contributions of Paul Sava and Sergey Fomel in laying the groundwork for the current theory and for ongoing RWE conversations, and I thank Biondo Biondi, Bob Clapp, Brad Artman, Paul Fowler, Tom Dickens, and Peter Traynor for enlightening discussions. The sponsors of the SEP consortium provided funding for this work.

REFERENCES

Albertin, U., D. Yingst, and H. Jaramillo, 2001, Comparing common-offset Maslov, Gaussian beam, and coherent state migrations: 71st Ann. Internat. Mtg., Soc. of Expl. Geophys., Expanded Abstracts, 913–916.

Brandsberg-Dahl, S. and J. Etgen, 2003, Beam-wave imaging: Soc. of Expl. Geophys., Expanded

Abstracts, 977–980. Bugayevskiy, L. M. and J. P. Snyder, 1995, Map projections - A Reference Manual: Taylor & Francis

Inc., Philadelphia, PA. Claerbout, J., 1985, Imaging the Earth’s Interior: Stanford University. Etgen, J., 2002, Waves, beams and Dimensions: An illuminating if incoherent view of the future of

migration: Soc. of Expl. Geophys., invited presentation. Fowler, P. J., 1995, Finite-difference solutions of the 3-D Eikonal equation in spherical coordinates:

65^{th} Annual SEG Convention, SEG Expanded Abstracts. Gazdag, J. and P. Sguazzero, 1984, Migration of seismic data by phase-shift plus interpolation: Geo-

physics, 69, 124–131. Guggenheimer, H., 1977, Differential Geometry: Dover Publications, Inc., New York. Hale, D., N. Hill, and J. Stefani, 1992, Imaging salt with turning seismic waves: Geophysics, 57,

1453–1462. Hill, N. R., 2001, Prestack Gaussian-beam depth migration: Geophysics, 66, 1240–1250.

Liseikin, V., 2004, A Computational Differential Geometry Approach to Grid Generation: Springer-Verlag, Berlin.

Nichols, D., 1994, Imaging complex structures using band-limited Green’s functions: Ph.D. thesis, Stanford University.

Sava, P. and S. Fomel, 2001, 3-D traveltime computation using Huygens wavefront tracing: Geophysics, 66, 883–889.

Sava, P. and S. Fomel, 2005, Riemannian wavefield extrapolation: Geophysics, 70, T45–T56.

Schneider, W. A., 1995, Robust and efficient upwind finite-difference traveltime calculations in three dimensions: Geophysics, 60, 1108–1117.

Shragge, J. and P. Sava, 2005, Wave-equation migration from topography: Soc. of Expl. Geophys., Expanded Abstracts, 1842–1845.

Shragge, J., 2006, Differential mesh generation: Soc. of Expl. Geophys., Expanded Abstracts, 2206– 2209.

Stoffa, P. L., J. T. Fokkema, R. M. de Luna Freire, and W. P. Kessinger, 1990, Split-step Fourier migration: Geophysics, 55, 410–421.

Tang, Y. and R. Clapp, 2006, Selection of reference-anisotropy parameters for wavefield extrapolation by lloyd’s algorithm: Soc. of Expl. Geophys., Expanded Abstracts, 189–192.

Zhang, Y., G. Zhang, and N. Bleistein, 2005, Theory of true-amplitude one-way wave equations and true-amplitude common-shot migration: Geophysics, 70, E3–E10.

APPENDIX A

Geometry in a generalized 3D Riemannian space is described by a symmetric metric tensor, g_{ij} D g_{ji}, that relates the geometry in a general non-orthogonal coordinate system, {ξ_{1},ξ_{2},ξ_{3}}, to an underlying Cartesian mesh, {x_{1},x_{2},x_{3}} (Guggenheimer, 1977). In matrix form, the metric tensor is written,

[g_{ij} ] D [ \begin{matrix}g_{11} & g_{12} & g_{13} \\ g_{21} & g_{22} & g_{23} \\ g_{31} & g_{32} & g_{33} \end{matrix} ] D [ \begin{matrix}g_{11} & g_{12} & g_{13} \\ g_{12} & g_{22} & g_{23} \\ g_{13} & g_{23} & g_{33} \end{matrix} ] ,  (A-1)

where g_{11}, g_{12}, g_{22}, g_{13}, g_{23} and g_{33} are functions linking the two coordinate systems through,

g_{11} D \frac{∂x_{k}}{∂ξ_{1}} \frac{∂x_{k}}{∂ξ_{1}}, g_{12} D \frac{∂x_{k}}{∂ξ_{1}} \frac{∂x_{k}}{∂ξ_{2}}, g_{22} D \frac{∂x_{k}}{∂ξ_{2}} \frac{∂x_{k}}{∂ξ_{2}},

g_{13} D \frac{∂x_{k}}{∂ξ_{1}} \frac{∂x_{k}}{∂ξ_{3}}, g_{23} D \frac{∂x_{k}}{∂ξ_{2}} \frac{∂x_{k}}{∂ξ_{3}}, g_{33} D \frac{∂x_{k}}{∂ξ_{3}} \frac{∂x_{k}}{∂ξ_{3}}.  (A-2)

The associated (or inverse) metric tensor, g^{ij}, is defined by g_{ij} D |g|g^{ij}, where |g| is metric tensor matrix discriminant. The associated metric tensor is given by,

[g^{ij} ] D \frac{1}{|g|} [ \begin{matrix}g_{22}g_{33}− g_{23}^{2} & g_{13}g_{23}− g_{12}g_{33} & g_{12}g_{23}− g_{13}g_{22} \\ g_{13}g_{23}− g_{12}g_{33} & g_{11}g_{33}− g_{13}^{2} & g_{12}g_{13}− g_{11}g_{23} \\ g_{12}g_{23}− g_{13}g_{22} & g_{12}g_{13}− g_{11}g_{23} & g_{11}g_{22}− g_{12}^{2} \end{matrix} ] ,  (A-3)

and has the following metric discriminant,

|g| D g_{33}(g_{11}g_{22}− g_{12}^{2}) [1 − \frac{g_{11}g_{23}^{2}C g_{22}g_{13}^{2}− 2g_{12}g_{23}g_{13}}{g_{33}(g_{11}g_{22}− g_{12}^{2})} ] .  (A-4)

A weighted metric tensor, m^{ij} D \sqrt{|g|}g^{ij}, is also used throughout the paper.

APPENDIX B

The extrapolation wavenumber defined in equations 11 and 12 cannot generally be implemented exactly in the Fourier domain due to a simultaneous spatial dependence (i.e. a function of both x_{1} and k_{ξ_{1}}). This can be addressed using an extended version of the split-step Fourier approximation (Stoffa et al., 1990) that uses Taylor expansions to separate k_{ξ_{3}} into two parts: k_{ξ_{3}} ≈ k_{ξ_{3}}^{P S} C k_{ξ_{3}}^{SSF}. Wavenumbers k_{ξ_{3}}^{P S} and k_{ξ_{3}}^{SSF} represent a pure Fourier (ω−k_{x}) domain phase-shift and a mixed (ω−x) domain split-step correction, respectively.

The phase-shift term is given by,

k_{ξ_{3}}^{P S} D −b_{1}k_{ξ_{1}}− b_{2}k_{ξ_{2}}C ib_{3}± [b_{4}^{2}ω^{2}− b_{5}^{2}k_{ξ_{1}}^{2}− b_{6}^{2}k_{ξ_{2}}^{2}− b_{7}k_{ξ_{1}}k_{ξ_{2}}C b_{8}i k_{ξ_{1}}C b_{9}i k_{ξ_{2}}− b_{10}^{2}] ^{\frac{1}{2}} , (B-1)

where b_{i} D b_{i}(ξ_{3}) are reference values of a_{i} D a_{i}(ξ_{1},ξ_{2},ξ_{3}). The split-step approximation is developed by performing a Taylor expansion about each coefficient a_{i} and evaluating the results at stationary reference values b_{i}. Assuming that the stationary values of k_{ξ_{1}} and k_{ξ_{2}} are zero, the split-step correction is as follows,

k_{ξ_{3}}^{SSF} D \frac{∂k_{ξ_{3}}}{∂a_{3}}|_{0}(a_{3}− b_{3}) C \frac{∂k_{ξ_{3}}}{∂a_{4}}|_{0}(a_{4}− b_{4}) C \frac{∂k_{ξ_{3}}}{∂a_{10}} |_{0}(a_{10}− b_{10}) ,  (B-2)

where “0” denotes "with respect to a reference medium". The partial differential expressions in equation B-2 are,

\frac{∂k_{ξ_{3}}}{∂a_{3}}|_{0} D b_{3}, \frac{∂k_{ξ_{3}}}{∂a_{4}}|_{0} D \frac{b_{4}ω^{2}}{\sqrt{b_{4}^{2}ω^{2}− b_{10}^{2}}} , \frac{∂k_{ξ_{3}}}{∂a_{10}} |_{0} D −\frac{b_{10}}{\sqrt{b_{10}^{2}ω^{2}− b_{10}^{2}}} ,  (B-3)

resulting in the following split-step Fourier correction wavenumber,

k_{ξ_{3}}^{SSF} D i b_{3} (a_{3}− b_{3}) C \frac{b_{4}ω^{2} (a_{4}− b_{4})}{\sqrt{b_{4}^{2}ω^{2}− b_{10}^{2}}} − \frac{b_{10} (a_{10}− b_{10})}{\sqrt{b_{4}ω^{2}− b_{10}^{2}}} .  (B-4)

APPENDIX C

The extrapolation wavenumber developed in equation 11 is appropriate for any non-orthogonal Riemannian geometry. However, there are a number of situations where symmetry or partial orthogonality are present. Moreover, one may wish to make a kinematic approximation where all of the imaginary components of the wavenumber are ignored. These situations are discussed in this Appendix.

3D Semi-orthogonal Coordinate Systems - Semi-orthogonal coordinate systems occur where one coordinate (ξ_{3}) is orthogonal to the other two coordinates (ξ_{1} and ξ_{2}) (Sava and Fomel, 2005). 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,

k_{ξ_{3}} D ia_{3}± [a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}− a_{6}^{2}k_{ξ_{2}}^{2}− a_{7}k_{ξ_{1}}k_{ξ_{2}}C ia_{8}k_{ξ_{1}}C ia_{9}k_{ξ_{2}}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-1)

where,

a D [ 0 0 \frac{n^{3}}{2m^{33}} \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}}} \sqrt{\frac{g^{22}}{g^{33}}} \frac{2g^{12}}{g^{33}} \frac{n^{1}}{m^{33}} \frac{n^{2}}{m^{33}} \frac{n^{3}}{m^{33}} ] ^{T}.  (C-2)

which are identical to the coefficients recovered by Sava and Fomel (2005).

3D Kinematic Coordinate Systems - Wave-equation migration amplitudes are generally inexact in laterally variant media - even in Cartesian-based systems. Hence, one beneficial approximation that reduces computational cost is to consider only the kinematic terms in equation 11. This approximate generates the following extrapolation wavenumber,

k_{ξ_{3}} D a_{1}k_{ξ_{1}}C a_{2}k_{ξ_{2}}± [a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}− a_{6}^{2}k_{ξ_{2}}^{2}− a_{7}k_{ξ_{1}}k_{ξ_{2}}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-3)

where,

a D [ −\frac{g^{13}}{g^{33}} − \frac{g^{23}}{g^{33}} 0 \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}} − ( \frac{g^{13}}{g^{33}}) ^{2}} \sqrt{\frac{g^{22}}{g^{33}} − ( \frac{g^{23}}{g^{33}}) ^{2}} \frac{2g^{12}}{g^{33}} − \frac{2g^{13}g^{23}}{(g^{33})^{2}} 0 0 \frac{n^{3}}{m^{33}} ] ^{T}. (C-4)

3D Kinematic Semi-orthogonal coordinate systems - Combining the two above restrictions yields the following extrapolation wavenumber,

k_{ξ_{3}} D ±[a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}− a_{6}^{2}k_{ξ_{2}}^{2}− a_{7}k_{ξ_{1}}k_{ξ_{2}}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-5)

where,

a D [ 0 0 0 \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}}} \sqrt{\frac{g^{22}}{g^{33}}} \frac{2g^{12}}{g^{33}} ; 0 0 \frac{n^{3}}{m^{33}} ] ^{T}.  (C-6)

2D Non-orthogonal coordinate systems - Two-dimensional situations are handled by identifying ξ_{2} D 0. All derivatives in the associated metric tensor g^{ij} with respect coordinate ξ_{2} are identically zero, and the resulting 2D non-orthogonal coordinate system wavenumber is,

k_{ξ_{3}} D a_{1}k_{ξ_{1}}C ia_{3}± [a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}C a_{8}ik_{ξ_{1}}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-7)

where,

a D [ −\frac{g^{13}}{g^{33}} 0 \frac{n^{3}}{2m^{33}} \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}} − ( \frac{g^{13}}{g^{33}}) ^{2}} 0 0 \frac{n^{1}}{m^{33}} 0 \frac{n^{3}}{m^{33}} ] ^{T}.  (C-8)

2D Non-orthogonal Kinematic Coordinate Systems - Two-dimensional kinematic situations are handled through identity ξ_{2} D 0. Again, all derivatives in the associated metric tensor g^{ij} with respect coordinate ξ_{2} are identically zero, and the 2D non-orthogonal kinematic extrapolation wavenumber is

k_{ξ_{3}} D a_{1}k_{ξ_{1}}± [a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-9)

where,

a D [ −\frac{g^{13}}{g^{33}} 0 0 \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}} − ( \frac{g^{13}}{g^{33}}) ^{2}} 0 0 0 0 \frac{n^{3}}{m^{33}} ] ^{T}.  (C-10)

2D Orthogonal Coordinate Systems - Two-dimensional situations are handled with ξ_{2} D g_{13} D 0. Accordingly, all derivatives in the associated metric tensor g^{ij} with respect coordinate ξ_{2} are identically zero, and the 2D non-orthogonal coordinate system is represented by

k_{ξ_{3}} D ia_{3}± [a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}C ia_{8}k_{ξ_{1}}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-11)

where,

a D [ 0 0 \frac{n_{3}}{2m^{33}} \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}}} 0 0 \frac{n^{1}}{m^{33}} 0 \frac{n^{3}}{m^{33}} ] ^{T}.  (C-12)

2D Orthogonal Kinematic Coordinate Systems - The above two approximations can be combined to yield the following extrapolation wavenumber for 2D orthogonal kinematic coordinate systems,

k_{ξ_{3}} D ±[a_{4}^{2}ω^{2}− a_{5}^{2}k_{ξ_{1}}^{2}− a_{10}^{2}] ^{\frac{1}{2}} ,  (C-13)

where,

a D [ 0 0 0 \frac{s}{\sqrt{g^{33}}} \sqrt{\frac{g^{11}}{g^{33}}} 0 0 0 0 \frac{n^{3}}{m^{33}} ] ^{T}.  (C-14)

LIST OF FIGURES

1 Illustration of problems with RWE computational grids. Top panel: Ray-coordinate system triplications. Bottom panel: Grid bunching for topographically conformal coordinate system.

2 Sheared Cartesian coordinate system test. Coordinate system (left panel) shear angle and velocity are θ D 25^{∘} and 1500 ms^{−1} , respectively. Zero-offset data consist of 4 flat plane-wave impulses at t D 0.2,0.4,0.6 and 0.8 s that are correctly imaged (right panel) at depths z D 300,600,900, and 1200 m.

3 Ellipsoidal polar coordinate system test example. Upper left: v(z) D 1500 C 0.35z velocity function overlain by a polar ellipsoidal coordinate system defined by parameter a D 1 C 0.2ξ_{3}− 0.05ξ_{3}^{2}. Upper right: Velocity model mapped in the RWE domain. Bottom right: Imaged reflectors in RWE domain. Bottom left: the RWE domain image mapped to a Cartesian mesh.

4 Example of triplication-free mesh generation. Top panel: Velocity model with three Gaussian velocity perturbations. Overlain is a coordinate mesh generated from ray-tracing. Note the triplication to either side of the shot-point, as well as the spreading beneath the shot point. Middle panel: Velocity model overlain by isochrons of an Eikonal equation solution for same shot-point. Bottom panel: Triplication-free computational mesh generated by through the approach specified in this section.

5 Example of wave-equation generated Green’s Functions on structured non-orthogonal mesh for a slice through the SEG-EAGE salt velocity model. Top left: Salt model in physical space with an overlain ray-based mesh. Top right: Velocity model in the RWE domain. Bottom right: Wavefield propagated in ray-coordinates through velocity model shown in the upper right. Bottom left: Wavefield in bottom right interpolated back to Cartesian space.

6 Results of a comparison test between Riemannian wavefield extrapolation (top panel) and Cartesian-based extrapolation (bottom panel).

7 Starting surface for Winkel-Tripel projection of unit hemisphere. Note that the cells on the unit hemisphere surface have nearly consistent area indicating fairly uniform sampling, and that not all parts of the unit hemisphere are sampled with this approach.

8 Wavefield extrapolation example through Elf-IFP-CGG synthetic velocity model.

Figure 1. 26 Figure 2.

27 Figure 3.

28 Figure 4.

29 Figure 5.

30 Figure 6.

31 Figure 7.

32 Figure 8.

33