The fast marching method: an effective tool for tomographic imaging and tracking multiple phases in complex layered media

The accurate prediction of seismic traveltimes is required in many areas of seismology, including the processing of seismic reflection profiles, earthquake location, and seismic tomography at a variety of scales. In this paper, we present two seismic applications of a recently developed grid-based numerical scheme for tracking the evolution of monotonically advancing interfaces, via finite-difference solution of the eikonal equation, known as the fast marching method (FMM). Like most other practical grid-based techniques, FMM is only capable of locating the first-arrival phase in continuous media; however, its combination of unconditional stability and rapid computation make it a truly practical scheme for velocity fields of arbitrary complexity. The first application of FMM that we present focuses on the prediction of multiple reflection and refraction phases in complex 2D layered media. By treating each layer that the wavefront enters as a separate computational domain, we show that sequential application of FMM can be used to track phases comprising any number of reflection and transmission branches in media of arbitrary complexity. We also show that the use of local grid refinement in the source neighbourhood, where wavefront curvature is high, significantly improves the accuracy of the scheme with little extra computational expense. The second application of FMM that we consider is in the context of 3D teleseismic tomography, which uses relative traveltime residuals from distant earthquakes to image wavespeed variations in the Earth’s crust and upper mantle beneath a seismic array. Using teleseismic data collected in Tasmania, we show that FMM can rapidly and robustly calculate two-point traveltimes from an impinging teleseismic wavefront to a receiver array located on the surface, despite the presence of significant lateral variations in wavespeed in the intervening crust and upper mantle. Combined with a rapid subspace inversion method, the new FMM based tomographic scheme is shown to be extremely efficient and robust. INTRODUCTION One of the classic problems in seismology is to accurately and robustly predict the traveltime and path of seismic energy between two points within a laterally heterogeneous 2D or 3D medium. Traditionally, this has been solved using geometric ray tracing based on a shooting or bending approach. Shooting methods of ray tracing (e.g., Julian and Gubbins, 1977; Sambridge and Kennett, 1990; Rawlinson et al., 2001) formulate the ray equation as an initial value problem that allows a complete ray to be traced if the source trajectory is specified. The boundary value problem of locating the required two-point path is then solved using an iterative update procedure. The bending method of ray tracing (e.g., Julian and Gubbins, 1977; Um and Thurber, 1987; Grechka and McMechan, 1996) iteratively adjusts the geometry of an initial arbitrary path that joins source and receiver until it becomes a true ray path (i.e., it satisfies Fermat’s principle). The principal drawbacks of ray tracing are related to robustness, speed, and ray selection. In the presence of even small velocity variations, both shooting and bending methods may fail to converge; this lack of robustness increases with the complexity of the medium. Ray tracing can also be a time-consuming process, particularly in the presence of a large number of sources and/or receivers. The final difficulty, that of ray selection, results from the potential existence of multiple two-point paths. Both shooting and bending methods do not necessarily converge to a globalminimum solution (i.e., the first arrival), and often it can be difficult to ascertain which arrival has been located. A more recently developed and increasingly popular class of method, particularly in the exploration industry, for predicting traveltimes in complex media is to seek finite-difference solutions to the eikonal equation throughout a gridded velocity field (e.g., Vidale, 1988; van Trier and Symes, 1991; Hole and Zelt, 1995; Buske and Kästner, 2004). Although this class of scheme is restricted to locating first arrivals only, the complete traveltime field can usually be computed extremely rapidly, which allows twopoint traveltimes, ray paths, and wavefront geometry to be easily extracted. The main drawback of finite-difference eikonal solvers is that they often suffer from stability problems; in particular, the progressive integration of traveltimes along an expanding square, which is commonly used to compute the traveltime field, has the potential to breach causality in the presence of large velocity gradients (Qin et al., 1992). However, recent developments with essentially non-oscillatory (ENO) finite difference schemes (e.g., Kim and Cook, 1999) have helped to address this shortcoming. The problem of calculating traveltimes to every point on a grid can be posed in terms of tracking the evolution of a monotonically advancing interface (e.g., a first-arrival seismic wavefront) that propagates throughout the medium from the source. The need to track an advancing interface is not limited to seismic wavefronts; many other areas of science require this problem to be solved. A recently introduced technique called the fast marching method or FMM (Sethian, 1996; Sethian and Popovici, 1999) was developed Research School of Earth Sciences Australian National University Canberra ACT 0200 Australia Phone: (02) 6125 0339 Facsimile: (02) 6257 2737 Email: [email protected] Presented at the 17th ASEG Geophysical Conference & Exhibition, August 2004. Revised manuscript received 9 June, 2005. 341 Exploration Geophysics (2005) Vol 36, No. 4 Exploration Geophysics (2005) 36, 341–350