By Alan Schiemenz
Nuclear test ban treaty monitoring, earthquake early warning systems, and knowledge of the earth’s deep interior are all made possible by seismic wave analysis. Numerical modeling of seismic waves is also a well-practiced art within the oil industry, where construction of high-resolution seismic velocity models improves seismic imaging workflows, thus reducing drilling risk and enhancing recovery of hydrocarbons. Full-waveform inversion (FWI) [8] is now the state-of-the-art seismic velocity model-building algorithm, with a growing research community supporting commercial development at many oil and oil service companies. FWI is a nonlinear optimization problem which iteratively updates the velocity model to reduce misfit between recorded and synthesized seismic data via the adjoint method. The problem was originally conceived in the early 1980s [7], although it was not until the late 2000s [6] that computers caught up to the intense demands of industrial application. Although FWI resides broadly in the geophysical domain, applied mathematicians stand to make the most meaningful contributions to the community’s unsolved problems. Among these contributions are improved wave equation solvers, accuracy of the inversion with respect to starting model, and robust strategies for multi-parameter inversion.
Seismic Waves
Elastic wave propagation is modeled by the equation of motion
\[\rho \: \partial^2_t\textbf{u} = \bigtriangledown \: \cdot \: \textbf{T} + \textbf{f}, \qquad (1) \]
where \(\textbf{u}\) is displacement of the seismic waveform, \(\rho\) is mass density, \(\textbf{f}\) represents body forces (e.g. earthquakes, explosions, air guns), and \(\textbf{T}\) is the elastic stress tensor. Hooke’s law prescribes a linear stress-strain constitutive relationship: \(\textbf{T} = \textbf{c} : \: \bigtriangledown \textbf{u}. \) The velocity model parameters are embedded within the tensor \(\textbf{c}.\) In the simplified case of an isotropic medium (i.e. \(\textbf{c}\) is invariant to rotation), there are two independent parameters: P-wave \((V_p)\) and S-wave \((S_p)\) velocity. Although the “speed of sound” in a medium is commonly conceptualized as a single scalar quantity, seismologists have long observed multiple types of wavefronts with greatly differing arrival time and amplitude behaviors. Primary (P) waves and secondary (S) waves are so named due to the order of their arrival \((V_p > V_s).\) S-waves do not exist in fluids, which led global seismologists to discover the fluid nature of the earth’s outer core. Marine exploration seismologists tend to use acoustic modeling engines (i.e. P-waves only), which are far cheaper and typically of sufficient accuracy, as the ocean layer removes much sensitivity to S-wave structure anyway. After acoustic velocity, anisotropic parameters are generally of greatest concern, and are crucial to validating field models with vertical profiles obtained from well logs.
Numerical Modeling
The acoustic approximation to (1) is a scalar equation and can be written in a semi-discrete form as
\[M \partial^2_t \phi = S \phi + f,\qquad (2) \]
where the mass and stiffness matrices \((M,S)\) are spatial-differential operators applied to the potential \(\phi,\) with \(\textbf{u} = \rho^{-1}\bigtriangledown\phi.\) These operators are defined by the choice of spatial discretization deployed. Finite difference methods have for decades been the workhorse of the oil industry, owing to ease of implementation and efficient parallelization. (2) is simulated for each of the seismic sources f, multiple times per FWI iteration, typically totaling 105-107 wave equation solves for a standard industry acquisition. In recent years, the seismology community has investigated more exotic schemes, such as the continuous finite element [5] and discontinuous Galerkin [2] methods. Element-based methods allow for local mesh refinement and adaptivity to irregular geological structure. Improving mesh construction remains a prime challenge for such methods in the context of FWI, where updated seismic velocity models necessitate a new grid for wave propagation.
Multi-step, explicit time integration (e.g. classic Runge-Kutta) schemes are the industry standard approach to time evolution. Alternatively, one may apply a Fourier transform to (2), yielding the sparse linear system \(B_\omega \tilde{\phi}(\textbf{x}, \omega) = \tilde{f}(\textbf{x}, \omega),\) where \(B_\omega \) is the discrete impedance matrix for frequency \(\omega.\) When a direct solve approach is tractable, the frequency-domain approach is superior, as all seismic sources can be cheaply simulated after inverting \(B_\omega.\) However, large-scale, three-dimensional problems typically exceed memory capacities, necessitating iterative methods instead.
Inversion of Model Parameters
Conventional FWI uses a least-squares objective function, where the velocity model \(\tilde{m}\) is derived by minimizing data misfit:
\[ \tilde{m} = \min\limits_{m} \frac {1}{2} \Sigma_{source \: i}\!\parallel \!\: d_i^{syn}(m) -d_i^{rec} \parallel ^2 : = \min\limits_{m} \chi (m).\qquad (3) \]
The synthetic data \(d_i^{syn}\) is computed by modeling (2) with a starting velocity model \(m_0,\) while \(d_i^{rec}\) is recorded data from a field survey. Gradient descent methods are used to calculate a model perturbation, where the steepest descent direction \(-\frac {\partial_\chi}{\partial(m)} \bigg|_{m=m_0} \) is computed from the cross-correlation of forward and adjoint-wavefield solutions to (2). Additional wave equation solves are used to compute a step length search, deriving a new model for the next iteration. As is common with gradient-based methods, convergence of FWI is sensitive to local minima in the objective function \(\chi.\) Denoted as “cycle-skipping” in the FWI community, the inversion will fail to converge to the global minimum when features of the synthetic waveform are more than half a wavelength out of phase with their recorded-data counterpart.
Figure 1. Results after applying standard acoustic FWI workflow to the Marmousi2 velocity model. Initial models are constructed by smoothing the true model. Model dimensions are 17 km width x 3.5 km depth. Final model accuracy depends strongly on the accuracy of the starting model.
Although FWI delivers highly accurate results when applied to a suitable starting model, FWI models derived from poor-quality starting models contain spurious cycle-skipping artifacts (see Figure 1). Better signal-to-noise ratio within the data (particularly at low frequencies) and longer acquisition geometries reduce the requisite starting model accuracy, but the fundamental problem remains. Quasi-Newton optimization schemes [4] may help by accounting for second-order scattering effects, at the expense of additional wave equation solves per iteration. Model-space strategies [e.g. 3] post-process the update by steering towards models of greater geological plausibility, while data-space strategies [e.g. 9] precondition the input to the objective function. Alternative objective functions [1] have likewise been explored in search of an error measure less sensitive to cycle-skipping.
Applied mathematics has played an integral role in the advancement of FWI technology, from the development of high-powered wave equation solvers to the innate nature of the optimization problem itself. The field is a promising one for applied mathematicians in search of multidisciplinary application, with a well-established foundation and many remaining open problems. Advancement beyond single-parameter inversion is a commercially lucrative proposition, but is currently restricted by the difficulty of mitigating parameter crosstalk and the cost involved to model additional physics. Avoiding convergence to local minima in the presence of poor data quality and/or starting model accuracy remains the top challenge in the field.
References
[1] Bozda, E., Trampert, J., & Tromp, J. (2011). Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements. Geophys. J. Int., 185(2), 845-870.
[2] Dumbser, M., & Käser, M. (2006). An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes–II. The three-dimensional isotropic case. Geophys. J. Int., 167, 319–336.
[3] Lewis, W., Amazonas, D., Vigh, D., & Coates, R. (2014). Geologically constrained full-waveform inversion using an anisotropic diffusion based regularization scheme: Application to a 3D offshore Brazil dataset. SEG Technical Program Expanded Abstracts, 1083-1088.
[4] Métivier, L., Brossier, R., Virieux, J., & Operto, S. (2013). Full waveform inversion and the truncated Newton method. SIAM J. Sci. Comput., 35(2), B401–B437.
[5] Peter, D., Komatitsch, D., Luo, Y., Martin, R., Le Goff, N., Casarotti, E.,…Tromp, J. (2011). Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes. Geophys. J. Int., 186, 721-739.
[6] Sirgue, L., Barkved, O.I., Dellinger, J., Etgen, J., Albertin, U., & Kommedal, J.H. (2010.) Full waveform inversion: the next leap forward in imaging at Valhall. First Break, 28, 65-70.
[7] Tarantola, A. (1984.) Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49, 1259–1266.
[8] Virieux, J., & Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics, 74, WCC1-WCC26.
[9] Warner M., & Guasch L. (2014.) Adaptive Waveform Inversion–FWI Without Cycle Skipping-Theory. 76th EAGE Conference and Exhibition 2014.
Alan Schiemenz is a data scientist at Crown Castle International Corp. He has enjoyed a diversity of career experiences enabled by a foundation of applied mathematics, including time as a member of the seismic modeling and inversion team at Schlumberger, and prior postdoctoral appointments in the computational earth sciences.