Introduction

Accurately imaging the Earth's subsurface using seismic wavefield data requires estimating realistic velocity profiles. Most information used to ascertain velocity models, however, derives from two sources: inverting a subset of the seismic wavefield (e.g., first-arrival tomography or waveform inversion), or from a measure of migration image quality (e.g., residual migration or migration velocity analysis). All approaches are nevertheless subject to the non-linearity of the seismic imaging problem and often require multiple iterations of velocity estimation, migration and image quality assessment before converging toward a well-resolved velocity model and seismic image.

One established data-domain approach for velocity estimation is ray-based, first-arrival, travel-time tomography, where earliest arrival times predicted by a linearized theory are matched to those picked from data. Travel-time discrepancies (or residuals) are back-projected along calculated ray-paths for all source and receiver combinations to obtain a velocity model update. Although iterative ray-based tomography is computationally efficient and can demonstrably recover longer-wavelength velocity profile structure, it has a limited ability to resolve finer-scale model components. Although this shortcoming is partially addressed by more "wave-like" corrections to the asymptotic ray theory, a more accurate (but more expensive) way to improve inverted velocity model resolution is using methods derived from finite-frequency approaches.

A popular alternative to ray-based tomography is waveform inversion based on wave-equation tomography. Assuming a linearized theory, this approach aims to match waveforms modeled by band-limited propagation operators to acquired seismic data. Importantly, higher-resolution velocity models can be recovered where a greater percentage of the seismic wavefield is used in the inversion. (In fact, velocity models developed from ray-based tomography often are the starting models for waveform inversion.) Waveform inversion can be implemented in either the time domain Bunks et al. (1985); Mora (1987); Shipp and Singh (2002); Tarantola (1984) or in the frequency domain Liao and McMechan (1996); Pratt and Worthington (1989); Sirgue and Pratt (2004) and for either acoustic and elastic wave equations. Frequency-domain approaches have the benefit that only a limited number of frequencies need be inverted Pratt and Worthington (1989), usually starting at lower frequencies then moving up a sparsely sampled spectrum. Moreover, because data are more linear with respect to model parameters at lower frequencies, convergence toward the global minimum is more likely Sirgue and Pratt (2004).

Two drawbacks to the general waveform inversion procedure are the memory requirements and computation complexity required to solve the 3D forward modeling problem. Most frequency-domain procedures employ finite-differences (FD) to solve the acoustic or elastic wave-equation, and often implement a LU decomposition of the correspondingly large, but relatively sparse, impedance matrix. The memory requirement for a LU decomposition of a 2D matrix is roughly proportional to nx⁴, where nx is the average model space dimension. This can be lowered to roughly nx³ using graph theoretic techniques able to exploit matrix sparsity Stekl and Pratt (1998). The memory requirement for 3D problem, however, theoretically rises to approximately nx⁶, which remains too costly for common exploration model sizes. However, novel iterative approaches may be able to reduce the cost of 3D FD modeling (see, e.g., Plessix (2006)).

An alternate way to lower the forward modeling costs is to replace FD modeling with an extrapolation scheme based on propagators derived from one-way wave-equations. Although one-way wavefield extrapolation is not as accurate, the numerical cost of a 3D implementation is roughly ns nx³ (log₂ (nx))² Biondi (2006), where ns is the number of shots. Moreover, the memory requirement is now of a similar or lower order of magnitude, and is no longer the impeding constraint. Accordingly, wavefield extrapolation potentially offers significant computational and memory savings and is worth examining in the context of waveform inversion. One-way extrapolation, nevertheless, comes with caveats. For example, long-offset seismic data usually contain turning waves, wide-angle reflections, and forward-scattering from shorter wavelength structure that carry important and complementary information about velocity structure. However, accurately modeling these waves with conventional wavefield extrapolation techniques remains difficult (though not impossible, see Zhang et al. (2006)). One approach is to use Riemannian wavefield extrapolation Sava and Fomel (2005) on a generalized coordinate mesh oriented in the general direction of turning-wave propagation. Because these grids are designed to incorporate the bulk of the turning-wave propagation directly into the coordinate system, it allows the user to implement lower-order extrapolation operators while achieving more accurate global propagation. In addition, one can avoid modeling at depths greater than those of the deepest turning waves, beyond which the transmission wavefield is insensitive to velocity variations Mulder and Plessix (2006).

This paper examines the use of one-way Riemannian wavefield extrapolation (RWE) operators in the forward modeling component of frequency-domain waveform inversion. The paper begins with a general review of the Pratt and Worthington (1989) approach. Subsequently, I describe the implementation of RWE forward modeling and present results of tests on the SMAART JV Pluto 1.5 dataset that indicate that RWE waveforms are fairly well matched with those from FD modeling at wider offsets. I then discuss results of a RWE waveform inversion scheme and demonstrate its ability to invert for a moderate (10$\%$) 1D velocity perturbation assuming an a priori constant velocity background.