Introduction

Imaging complex geology is one of the main challenges of today's seismic processing. Of the many seismic imaging methods available, downward continuation Claerbout (1985) has proven to be accurate, robust, and capable of handling models with large and sharp velocity variations. In addition, such methods naturally handle the multipathing which occurs in complex geology and provide a band-limited solution to the seismic imaging problem. Furthermore, as computational power increases, such methods are gradually moving into the mainstream of seismic processing.

However, migration by downward continuation imposes strong limitations on the dip of reflectors that can be imaged since, by design, it favors energy which is propagating mainly in the downward direction. Upward propagating energy, e.g., overturning waves, can be imaged in principle using downward continuation methods Hale et al. (1992), although the procedure is difficult, particularly for prestack data. In contrast, Kirchhoff-type methods based on ray-traced traveltimes can image steep dips and handle overturning waves, although those methods are far less reliable in complex environments given their high-frequency asymptotic nature.

The steep-dip limitation of downward continuation techniques has been addressed in several ways:

• A first option is to increase the angular accuracy of the extrapolation operator, for example by employing methods from the Fourier finite-difference (FFD) family Biondi (2002); Ristow and Ruhl (1994), or the Generalized Screen Propagator (GSP) family Huang and Wu (1996); de Hoop (1996). The enhancements brought about by these methods come at a price, since they increase the cost of extrapolation without guaranteeing unconditional stability.
• A second option is to perform the wavefield extrapolation in tilted coordinate systems Etgen (2002), or by designing sources which favor illumination of particular regions of the image Chen et al. (2002); Rietveld and Berkhout (1994). We can thus increase angular accuracy, although those methods favor a subset of the model (a salt flank, for example), while potentially decreasing the accuracy in other regions. In complex models it is also not obvious what is an optimal tilt angle for the extrapolation grid.
• A third possibility is hybridization of wavefield and ray-based techniques, either in the form of Gaussian beams Gray et al. (2002); Hill (1990, 2001), coherent states Albertin et al. (2001, 2002), or beam-waves Brandsberg-Dahl and Etgen (2003). Such techniques are quite powerful, since they couple wavefield methods with multipathing and band-limited properties, with ray methods, which deliver arbitrary directions of propagation, even overturning. The main strength of those techniques, is also their main weakness, i.e. they are formulated in terms of decoupled beams. Beams may leave shadow zones in various parts of the model, which hamper their imaging abilities. Furthermore, beams have limited size, which in turn limits the extent of the diffractions created by sharp features in the model to that of any particular beam, no matter how accurate the extrapolator within each beam is. In addition, the narrow extrapolation domain poses serious beam superposition problems, such as beam boundary effects.

The main idea of our paper is to couple the beams together and extrapolate within all of them at once. We, therefore, cannot talk about beams anymore, but instead we need to talk about continuously changing coordinate systems. We extend downward continuation in a regular Cartesian space to wavefield extrapolation in distorted coordinates, known as Riemannian spaces , thus the name of our method. We formulate the theory in arbitrary 3-D semi-orthogonal Riemannian spaces, e.g., ray coordinates, although those coordinate systems do not necessarily need to have a physical meaning as long as they fulfill the semi-orthogonality condition. Examples of such coordinate systems include, but are not limited to, fans of rays emerging from a source point, or bundles of rays initiated by plane waves of arbitrary initial dips at the source. A special case of our method is represented by the polar/spherical coordinate system Nichols and Palacharla (1994); Nichols (1994, 1996).

Our method can be seen as a finite-difference solution to the acoustic wave-equation in ray coordinates. In this respect, it is closely related to Huygens wavefront tracing Sava and Fomel (2001), which represents a finite-difference solution to the eikonal equation in ray coordinates.

Another idea related to our method is that of wave-equation in ray-centered coordinates Cervený (2001); Yedlin (1981). However, our method is different since the ray-centered coordinate system is orthogonal in 3-D and defined around an individual ray. In contrast, our method is formulated in ray coordinates which are parameterizations of the wavefields at the source, and which are non-orthogonal in 3-D and defined globally for an entire family of rays.

The upside of our method is that the coordinate system may follow the waves, and can even overturn, such that we can use one-way extrapolators to image diving waves (Figure 1).

 

overturned Figure 1 Ray coordinate systems are superior to tilted coordinate systems for imaging overturning waves using one-way wavefield extrapolators. Overturning reflected energy may become evanescent in tilted coordinate systems (a), but stays non-evanescent in ray coordinate systems (b).

We can also use extrapolators with small angle accuracy (e.g. $15^\circ$), since, in principle, we are never too far from the actual direction in which waves propagate. We are also not confined to the extent of any individual beam, therefore we can track diffractions for their entire spatial extent.