By Bill Symes
Bill Symes, Rice University.
The following is a short introduction to an invited lecture to be presented at the upcoming 2018 SIAM Annual Meeting (AN18) in Portland, Ore., from July 9-13. Look for feature articles by other AN18 invited speakers introducing the topics of their talks in future issues.
Most of what we know about Earth’s interior comes from measuring physical fields on or near the surface and inferring subsurface structures from these measurements. Apart from inspecting outcrops, the optical band of the electromagnetic spectrum is useless—you cannot actually see through rock, despite the title of this talk—but there are many other choices for probing fields. Seismic fields (roughly speaking, elastic waves in Earth) have a few key advantages as subsurface probes. They cover relatively long distances—measurable waves of the lowest frequency excite all of Earth—and combine this range with resolution, propagating a few hundred wavelengths before disappearing into Earth’s background seismic noise. To some extent, they exhibit space-time locality and directionality (that is, locality in phase space) and definite velocity, which makes them waves.
Whether elastodynamics or a more complex model of wave physics is adequate to predict measured seismic data, and what information about Earth’s structure might result from doing so, is an inverse problem. One can only answer these questions affirmatively by exhibiting a distribution of elastic (or other) parameters for which the solution of the associated partial differential equations (PDEs) actually predicts the seismic measurements to satisfactory accuracy. Joe Keller proposed this “find the question given the answer” definition of “inverse problem” in his 1976 American Mathematical Monthly article. However, Keller’s examples do not clarify how instances of such problems involving nonlinear relations between billions of data points and billions of parameters, mediated by complex systems of PDEs, might be wrestled into computational practicality.
As I will explain in my 2018 SIAM Annual Meeting talk, answers to this question began taking shape in the 1980s, as Albert Tarantola and many others elaborated the “outer loop” point of view mentioned in Bruce Hendrickson’s November 2017 SIAM News piece. The outer loop perspective posits that posing data fit as an optimization problem over choices of elastic parameters invites selection by iteration over solutions of the elastodynamic (or similar) wave equations. Tarantola’s classic text, Inverse Problem Theory and Methods for Model Parameter Estimation, republished by SIAM in 2005, shows how tools from control theory make the computation of artifacts needed for local (descent-based) optimization more feasible.
The first practicable algorithms required another 20 years of computational and mathematical advances and inspired a very active research area known as full waveform inversion, with branches in industrial and academic seismology. The computational needs were obvious consequences of the world’s three-dimensionality. The mathematical issues were more surprising — already evident in the 1980s, they stem from the simple observation that perturbation of wave speed effectively differentiates the waveform. Combined with the nonlinearity of the elastic wavefield as a function of the coefficients in the (linear) wave equation, this hypersensitivity with respect to wave velocities—which are functions of the coefficients—can cause local optimization to stagnate far from a useful elastic model of Earth, no matter how many Pflops are expended. I will review some of the numerous ideas advanced to overcome this obstacle, and show why they might work via simple examples.