Research paperSeisElastic2D: An open-source package for multiparameter full-waveform inversion in isotropic-, anisotropic- and visco-elastic media
Introduction
In recent decades, benefiting from the significant growth in power of high-performance computing (HPC) technology, and in seismic data acquisition technology, full-waveform inversion (FWI) methods are increasingly applied to exploration geophysics for oil and gas detection (Tarantola, 1984, Pratt et al., 1998, Virieux and Operto, 2009, Biondi and Almomin, 2014, Yuan and Simons, 2014). Multiparameter FWI holds the promise of providing high-resolution subsurface properties in terms of elasticity, anisotropy and attenuation, which are essential for lithology discrimination and hydrocarbon reservoir characterization (Pan et al., 2018b, Pan et al., 2019). However, field data application of multiparameter elastic FWI remains a challenging task.
When simultaneously determining multiple physical parameters, the interparameter trade-offs increase nonlinearity and uncertainty of the inverse problem and complicate multiparameter FWI significantly (Operto et al., 2013, Innanen, 2014, Biondi et al., 2017, Pan et al., 2018a). Another challenge is known as cycle-skipping arising from the highly nonlinear relationship between seismic data and model properties (Bunks et al., 1995, Leeuwen and Herrmann, 2013, Yuan and Simons, 2014, Chi et al., 2015, Li and Demanet, 2016, Liu et al., 2018, Yao et al., 2018, Yong et al., 2018, Wu et al., 2019). With poor initial models and lack of low frequencies, solving the inverse problem using local optimization methods always fails to converge to the global minimum. In large-scale seismic experiments, iterative solution of the inverse problem involves a huge number of forward and adjoint simulations (Castellanos et al., 2015, Zhang et al., 2020), which is extremely expensive especially for 3D elastic media (Vigh and Starr, 2008, Borisov and Singh, 2015, Biondi et al., 2018).
A range of open-source FWI packages have been released, including FWT2D (Sourbier et al., 2009a, Sourbier et al., 2009b), the SEISCOPE optimization toolbox (Métivier and Brossier, 2016), SeisFlows (Modrak et al., 2018), DENISE (Köhn, 2011), SimpleFWI (https://github.com/TristanvanLeeuwen/SimpleFWI), pysit (https://github.com/pysit/pysit), etc. In this paper, we present an open-source package SeisElastic2D, which consists of a set of FORTRAN 90 routines and SHELL scripts, for multiparameter FWI in isotropic-, anisotropic- and visco-elastic media. This package is designed with a modular structure based upon several existing open-source packages including upgraded SPECFEM2D, seisDD and Seismic Unix (SU), which provide the solutions for the forward and adjoint wave equations, flexible inversion workflows and efficient data processing tools. Relative to other open-source packages, this package focuses on providing an ability to vary model parameterizations for isotropic- and anisotropic-elastic FWI, and to vary misfit functions for reducing cycle-skipping and source-structure trade-offs, efficient strategy for attenuation estimation, etc.
Subsurface elastic rock properties (e.g., P-wave velocity to S-wave velocity ratio) are key indicators of hydrocarbon resources (Brown and Korringa, 1975, Castagna et al., 1985, Russell et al., 2003, Gao and Huang, 2019, Fang et al., 2020). An appropriate model parameterization, alongside an optimized inversion strategy, is critical to obtain such elastic properties reliably. In SeisElastic2D, various isotropic-elastic model parameterizations such as velocity-density (DV) and modulus-density (DM) are provided. These parameterizations have been used to optimize isotropic-elastic FWI using both walk-away vertical seismic profile (W-VSP) data (Pan et al., 2018b) and surface reflection seismic data (Pan et al., 2019). Beyond this, anisotropy produces influences on the traveltimes and amplitudes of the waveforms as the seismic waves propagate in real Earth media (Thomsen, 1986, Tsvankin, 1996). Inclusion of anisotropy in elastic FWI can produce more accurate velocity models (Barnes et al., 2008, Lee et al., 2010, Prieux et al., 2011). High-resolution anisotropy models are also important to characterize subsurface lithology distribution (Li et al., 2016). Determining the optimal model parameterization is essential for anisotropic-elastic FWI to estimate the velocity and anisotropy models reliably (Plessix and Cao, 2011, Alkhalifah, 2016). Various model parameterizations for elastic FWI in vertical transverse isotropic (VTI) and tilted transversely isotropic (TTI) media are included in SeisElastic2D. With SeisElastic2D, Pan et al. (2020) examined the performances of different model parameterizations for VTI-elastic FWI using practical W-VSP data. In this paper, we introduce three model parameterizations for TTI-elastic FWI and give numerical examples produced using SeisElastic2D to illustrate the importance of anisotropy in elastic FWI.
The standard FWI misfit function, measuring waveform-difference (WD) between observed and synthetic seismic data, although powerful, suffers from cycle-skipping. Wave equation traveltime inversion with a cross-correlation (CC) misfit function helps reduce nonlinearity of the inverse problem by minimizing the traveltime differences (Luo and Schuster, 1991, Leeuwen and Mulder, 2010). The envelope-difference (ED) misfit function also holds ability of overcoming the cycle-skipping difficulty (Bozdağ et al., 2011, Wu et al., 2014, Yuan et al., 2015, Borisov et al., 2017). Thus, in early stages of FWI, CC and ED misfit functions can be used to recover long wavelength model structures. Subsurface attenuation models are essential to compensate the amplitude loss and phase delay in seismic imaging (Innanen and Lira, 2010, Margrave et al., 2011, Guo et al., 2016, Chen et al., 2017) and characterize attenuative targets (e.g., fluid-filled/fractured reservoirs) (Innanen, 2011, Chen et al., 2018, Pan and Innanen, 2019, da Silva et al., 2019). Compared to the traditional spectral-ratio and central-frequency shift methods (Bath, 1974, Tonn, 1991, Quan and Harris, 1997), FWI also holds the ability of providing more accurate attenuation models in complex medium (Tarantola, 1998). However, the inverted attenuation models may suffer from significant velocity-attenuation trade-off artifacts (Brossier, 2011, Kamei and Pratt, 2013, Fabien-Ouellet et al., 2017, Keating and Innanen, 2019). SeisElastic2D provides several amplitude-based (e.g., root-mean-square amplitude-ratio) misfit functions for attenuation estimation in visco-elastic FWI, which show stronger sensitivities to attenuation anomalies and are able to reduce the trade-off artifacts caused by velocity errors. With SeisElastic2D, Pan and Innanen (2019) applied visco-elastic FWI with amplitude-based misfit functions to practical W-VSP data for attenuative reservoir characterization. However, seismic amplitudes are also affected by other factors such as source mechanism and reflection/transmission loss, which may produce unexpected uncertainties for attenuation inversion. As seismic attenuation reduces high-frequencies more than low-frequencies leading to central-frequency shift of the seismic data, in this paper, a new central-frequency shift (CS) misfit function in SeisElastic2D is designed to estimate the attenuation models in visco-elastic media. Compared to the WD and amplitude-based misfit functions, the CS misfit function is expected to invert for the attenuation models more reliably in the presence of velocity errors and amplitude fluctuations.
Newton-based optimization methods are promising as the inverse multiparameter Hessian accounts for and corrects a significant fraction of interparameter trade-off issues (Operto et al., 2013, Innanen, 2014). However, iteratively solving the Newton equation with linear conjugate gradient algorithms at each nonlinear iteration causes a large jump in computational expense. In this package, the standard steepest descent, nonlinear conjugate gradient, and -BFGS optimization methods (Nocedal and Wright, 2006) are developed. The diagonal Hessian can be applied for preconditioning, which works to compensate for imperfect illumination, accelerates the convergence rate, and balances the updates across different physical parameters (Modrak and Tromp, 2016, Chen and Sacchi, 2017). SeisElastic2D can be used on both large HPC cluster and local workstation supporting massively parallel interface (MPI), to accommodate the high computational requirements of large-scale inverse problem. These tools of SeisElastic2D help overcome the difficulties in multiparameter elastic FWI and bridge the gap between academic studies and industrial applications.
This paper is organized as follows. We first briefly review the theory of elastic FWI with the adjoint-state method. Different model parameterizations for isotropic-, VTI- and TTI-elastic FWI, optimization methods, misfit functions, and the two-stage inversion approach for visco-elastic FWI are introduced. In the numerical modeling section, we first give synthetic examples to illustrate the importance of anisotropy for elastic FWI. Finally, we apply visco-elastic FWI with the two-stage inversion strategy and new CS misfit function to construct near-surface velocity and attenuation models by inversion of surface and body waves.
Section snippets
Multiparameter elastic FWI with the adjoint-state method
The standard FWI algorithm estimates the subsurface properties iteratively by minimizing waveform-difference (WD) between observed data and synthetic seismic data (Tarantola, 1984, Pratt et al., 1998, Tromp et al., 2005, Virieux and Operto, 2009). The inverse problem is solved subject to the elastodynamic wave equation, which gives the augmented Lagrangian function (Liu and Tromp, 2006, Plessix, 2006, de Ridder and Maddison, 2018, Pan et al., 2019):
Structure of the package and inversion workflow
The package is designed with a modular structure, as shown in Fig. 1. The numerical experiments, including forward modeling, kernel calculation, and inversion, can be started by running the job submission script (submit.sh) on HPC cluster or a local workstation. A set of key parameters need to be set in the parameter file. For examples:
The upgraded SPECFEM2D package provides the forward and adjoint simulation solver for calculating the sensitivity kernels and diagonal pseudo-Hessian
Numerical experiments
In this section, we first give synthetic examples to illustrate the importance of anisotropy for estimating velocity models in elastic FWI. Then, we carry out visco-elastic FWI with the two-stage inversion and new CS misfit function for near-surface velocity and attenuation estimation using surface and body waves.
Discussions and conclusions
We present a flexible open-source package SeisElastic2D for isotropic-, anisotropic- and visco-elastic FWI. This package provides various model parameterizations for isotropic-, VTI- and TTI-elastic FWI. Different misfit functions are provided to overcome the cycle-skipping problem, reduce source-structure trade-offs, and invert for attenuation models reliably. With flexible inversion workflows and parallel computing, we have applied isotropic-, VTI- and visco-elastic FWI to the practical
Code availability
The software package, user manual, test example and output results are available at https://github.com/PanIGGCAS/SeisElastic2D_1.1.
CRediT authorship contribution statement
Wenyong Pan: Design and implement the algorithms, Finish the numerical experiments, Prepare the manuscript. Kristopher A. Innanen: Funding support, Ideas, Suggestions for improving the manuscript. Yanfei Wang: Funding support, Parallel computing facilities, Suggestions for organizing the code package.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This research was partially supported by the Consortium for Research in Elastic Wave Exploration Seismology (CREWES) and National Science and Engineering Research Council of Canada (NSERC, CRDPJ 461179-13), and in part from the Canada First Research Excellence Fund. Wenyong Pan is also supported by the IGGCAS Research Start-up Founds E0515402, IGGCAS grant 2019031 and CAS innovation program ZDBS-LY-DQC003. Thanks greatly to Dr. Yanhua O. Yuan (ExxonMobil) for contributing initial codes of
References (89)
Two-dimensional frequency-domain visco-elastic full waveform inversion: Parallel algorithms, optimization and performance
Comput. Geosci.
(2011)- et al.
Time-domain seismic modeling in viscoelastic media for full waveform inversion on heterogeneous computing platforms with OpenCL
Comput. Geosci.
(2017) - et al.
SeisFlows-Flexible waveform inversion software
Comput. Geosci.
(2018) - et al.
FWT2D: A massively parallel program for frequency-domain full-waveform tomography of wide-aperture seismic data—Part 1: Algorithm
Comput. Geosci.
(2009) - et al.
FWT2D: A massively parallel program for frequency-domain full-waveform tomography of wide-aperture seismic data—Part 2: Numerical examples and scalability analysis
Comput. Geosci.
(2009) - et al.
A new full waveform inversion method based on shifted correlation of the envelope and its implementation based on OPENCL
Comput. Geosci.
(2019) Research note: Insight into the data dependency on anisotropy: An inversion prospective
Geophys. Prospect.
(2016)- et al.
Feasibility study for an anisotropic full waveform inversion of crosswell seismic data
Geophys. Prospect.
(2008) Spectral Analysis in Geophysics: Developments in Solid Earth Geophysics
(1974)- et al.
Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion
Geophysics
(2014)