Improved wavefield extrapolation can be achieved in many ways. First, one may improve the high-angle accuracy of an operator while retaining a Cartesian computational grid. Examples of this include incorporating higher-order terms in the expansion of Fourier domain operators (e.g. Fourier finite-difference Ristow and Ruhl (1994), generalized screen propagator de Hoop et al. (2000)), or using tilted Cartesian coordinate systems Etgen (2002); Shan and Biondi (2004); Zhang and McMechan (1997) that extend the accuracy of high-angle propagation. Second, seismic wavefields may be spatially partitioned into more manageable sections and then independently extrapolated in preferred directions. For example, decomposing data sections to form local beams for extrapolating along tubes of finite thickness Albertin et al. (2001); Brandsberg-Dahl and Etgen (2003); Gray et al. (2002); Hill (2001).
A third option is to abandon the strictures of Cartesian coordinates altogether and represent the physics of one-way wavefield extrapolation in a generalized coordinate system that obeys the tenets of differential geometry Guggenheimer (1977). In particular, one could use a basis (or coordinate system) that conforms to where the wavefronts propagate Sava and Fomel (2004). In this reference frame low-angle operators remain applicable, and the extrapolation procedure is of high fidelity, even at arbitrarily large angles to depth axis. The strategy espoused in this paper is the latter: it is more prudent to adjust the coordinate system to conform better with the physics than to force the physics to work in Cartesian coordinates, or on an a priori spatial partition of the data or model space.
One judicious choice of non-Cartesian coordinate system is a basis derived from a suite of rays. In this approach, the natural wavefield extrapolation direction is travel-time along a ray, with orthogonal coordinates directed across the rayfront at a constant time step. However, unlike for Cartesian coordinates, the distance between adjacent rays may freely expand or contract according to the lateral variations in the velocity model. Thus, properly defining the coordinate metric requires additional parameters that account for the Jacobian-like coordinate spreading. Given a rayfield and the associated Jacobian parameters, the solution to the corresponding one-way acoustic wave-equation is generated in an ordinary fashion Sava and Fomel (2004). The ray-coordinate wavefield solution is then interpolated back to a Cartesian mesh through a simple mapping operation.
One of the practical difficulties of ray-coordinate-based wavefield extrapolation is developing a robust procedure for handling triplicating rayfields that naturally arise due to wavefield multipathing. In particular, we need to prevent numerical instabilities from arising when calculating coordinate Jacobian spreading and related parameters that require computing finite-difference derivatives at the ray-crossing locations. Sava and Fomel (2004) apply a regularization parameter that prevents division by zero. This procedure, though, can lead to anomalous extrapolation amplitudes, which motivates us to seek out new methods for calculating rayfields and circumventing the ubiquitous problem of ray-coordinate triplication.
Underlying ray-coordinate systems may be generated by assuming that the rayfield is frequency-independent, and computing the solution to the wide-band eikonal equation Cervený (2001). However, this approximation can be inappropriate for complex geology where a stationary ray-coordinate system inadequately describes monochromatic wave propagation over a range of frequencies. One example is the significant frequency-dependence of rayfield illumination across salt-sediment interfaces characterized by large impedance contrasts and rugose topography. Hence, an important question is how does one expect to maintain a sufficient and consistent wavefield illumination when the underlying rayfield is itself strongly dependent on frequency? Hence, a frequency-dependent ray-coordinate systems should be an invaluable tool for enhancing imaging practice in complex media.
This paper presents a procedure for constructing a frequency-dependent ray-coordinate system in an adaptive manner directly from the wavefield. The key idea is that the rayfront vectors at any given step are directly calculable from the phase-gradient of previous wavefield solution steps. This naturally leads to a bootstrapping procedure where one alternates between calculating the coordinate system for the next step, and the corresponding wavefield solution at that step. The methodology is similar to the Riemannian wavefield extrapolation (RWE) technique presented by Sava and Fomel (2004), where rayfields are traced through a smoothed velocity model using a Huygens' wavefront tracer Sava and Fomel (2001). The method differs, though, in that frequency-stationarity of the rayfield is not assumed, and the rayfield is instead calculated on-the-fly from the monochromatic wavefield. This method also differs from Shragge and Biondi (2003) in that an initial wavefield is not required as a precondition for solution. Also included in this report is a companion paper, Shragge and Biondi (2004), that discusses the strategy of using wavefield solutions precomputed on a background velocity model to train an updated ray-coordinate system using phase-rays.
We begin this paper with a review of phase-ray theory, frequency-dependent coordinate system generation, and ray-coordinate wavefield extrapolation. Then, we introduce the bootstrap procedure by which the ray-coordinate system and accompanying wavefield solutions are computed. Next, we show examples of wavefields extrapolated in adaptive phase-ray coordinates, and conclude with a discussion of the complications posed by triplicating coordinate systems. A more general formulation involving the oriented wave equation Fomel (2003) has the potential to address this problem in a robust theoretical framework, although such opportunity remains subject to future research.