Full-wavefield seismic inversion (FWI) estimates a subsurface elastic model by iteratively minimizing the difference between observed and simulated data. This process is extremely computationally intensive, with a cost comparable to at least hundreds of prestack reverse-time depth migrations. When FWI is applied using explicit time-domain or frequency-domain iterative-solver-based methods, the seismic simulations are performed for each seismic-source configuration individually. Therefore, the cost of FWI is proportional to the number of sources. We have found that the cost of FWI for fixed-spread data can be significantly reduced by applying it to data formed by encoding and summing data from individual sources. The encoding step forms a single gather from many input source gathers. This gather represents data that would have been acquired from a spatially distributed set of sources operating simultaneously with different source signatures. The computational cost of FWI using encoded simultaneous-source gathers is reduced by a factor roughly equal to the number of sources. Further, this efficiency is gained without significantly reducing the accuracy of the final inverted model. The efficiency gain depends on subsurface complexity and seismic-acquisition parameters. There is potential for even larger improvements of processing speed.

Fast full-wavefield seismic inversion using encoded sources

We present a review of the recent development in finite-difference time-domain modeling of seismic wave propagation and earthquake motion. The finite-difference method is a robust numerical method applicable to structurally complex media. Due to its relative accuracy and computational efficiency it is the dominant method in modeling earthquake motion and it also is becoming increasingly more important in the seismic industry and for structural modeling. We first introduce basic formulations and properties of the finite-difference schemes including promising recent advances. Then we address important topics as material discontinuities, realistic attenuation, anisotropy, the planar free surface boundary condition, free-surface topography, wavefield excitation (including earthquake source dynamics), non-reflecting boundaries, and memory optimization and parallelization.

The finite-difference time-domain method for modeling of seismic wave propagation

SUMMARY
In this study, we propose a new numerical method, named as Traction Image method, to accurately and efficiently implement the traction-free boundary conditions in finite difference simulation in the presence of surface topography. In this algorithm, the computational domain is discretized by boundary-conforming grids, in which the irregular surface is transformed into a ‘flat’ surface in computational space. Thus, the artefact of staircase approximation to arbitrarily irregular surface can be avoided. Such boundary-conforming gridding is equivalent to a curvilinear coordinate system, in which the first-order partial differential velocity-stress equations are numerically updated by an optimized high-order non-staggered finite difference scheme, that is, DRP/opt MacCormack scheme. To satisfy the free surface boundary conditions, we extend the Stress Image method for planar surface to Traction Image method for arbitrarily irregular surface by antisymmetrically setting the values of normal traction on the grid points above the free surface. This Traction Image method can be efficiently implemented. To validate this new method, we perform numerical tests to several complex models by comparing our results with those computed by other independent accurate methods. Although some of the testing examples have extremely sloped topography, all tested results show an excellent agreement between our results and those from the reference solutions, confirming the validity of our method for modelling seismic waves in the heterogeneous media with arbitrary shape topography. Numerical tests also demonstrate the efficiency of this method. We find about 10 grid points per shortest wavelength is enough to maintain the global accuracy of the simulation. Although the current study is for 2-D P-SV problem, it can be easily extended to 3-D problem.

Traction image method for irregular free surface boundaries in finite difference seismic wave simulation

[1] The southernmost San Andreas fault has a high probability of rupturing in a large (greater than magnitude 7.5) earthquake sometime during the next few decades. New simulations show that the chain of sedimentary basins between San Bernardino and downtown Los Angeles form an effective waveguide that channels Love waves along the southern edge of the San Bernardino and San Gabriel Mountains. Earthquake scenarios with northward rupture, in which the guided wave is efficiently excited, produce unusually high long-period ground motions over much of the greater Los Angeles region, including intense, localized amplitude modulations arising from variations in waveguide cross-section.

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2005GL025472

Strong shaking in Los Angeles expected from southern San Andreas earthquake

SUMMARY

Perfectly matched layers (PMLs) provide an exponential decay, independent of the frequency, of any propagating field along an assigned direction without producing spurious reflections at the interface with the elastic volume. For this reason PMLs have been applied as absorbing boundary conditions (ABCs) and their efficiency in attenuating outgoing wave fields on the outskirts of numerical grids is more and more recognized. However, PMLs are designed for first-order differential equations and a natural extension to second-order Partial Differential Equations (PDEs) involves either additional variables in the time evolution scheme or convolutional operations. Both techniques are computationally expensive when implemented in a spectral element (SE) code and other ABCs (e.g. paraxial or standard sponge methods) still remain more attractive than PMLs.



Here, an efficient second-order implementation of PMLs for SE is developed from interpreting the Newmark scheme as a time-staggered velocity–stress algorithm. The discrete equivalence with the standard scheme is based on an L2 approximation of the stress field, with the same polynomial order as the velocity. In this case, PMLs can be introduced as for first-order equations preserving the natural second-order time stepping. Subsequently, a non-classical frequency-dependent perfectly matched layer is introduced by moving the pole of the stretching along the imaginary axis. In this case, the absorption depends on the frequency and the layer switches from a transparent behaviour at low frequencies to a uniform absorption as the frequency goes to infinity. It turns out to be more efficient than classical PMLs in the absorption of the incident waves, as the grazing incidence approaches, at the cost of a memory variable for any split component. Finally, an extension of PMLs to general curvilinear coordinates is proposed.

The Newmark scheme as velocity–stress time-staggering: an efficient PML implementation for spectral element simulations of elastodynamics

[1] Robust absorbing boundary conditions are central to the utility and advancement of 3-D numerical wave propagation methods. It is in general preferred that an absorbing boundary method be capable of broadband absorption, be efficient in terms of memory and computation time, and be widely stable in connection with sophisticated numerical schemes. Here we discuss these issues for a promising absorbing boundary method, perfectly matched layers (PML), as implemented in the widely used fourth-order accurate three-dimensional (3-D) staggered-grid velocity-stress finite difference (FD) scheme. Numerical results for point (explosive and double couple) and extended sources, velocity structures (homogeneous, 1-D and 3-D), and different thickness PML zones are excellent, in general, leaving no observable reflections in PML seismograms compared to the amplitudes of the primary phases. For both homogeneous half-space and 1-D models, typical amplitude reduction factors (with respect to the maximum trace amplitude) range between 1/100 and 1/625 for PML thicknesses of 5–20 nodes. A PML region of thickness 5 outperforms a simple exponential damping region of thickness 20 in a homogeneous half-space model by a factor of 3. We find that PML is effective across the simulation bandwidth. For example, permanent offset artifacts due to particularly poor absorption of long-period energy by the simple exponential damping are effectively absent when PML is used. The computational efficiency and storage requirements of PML, compared to the simple exponential damping, are reduced due to the need for only narrow absorbing regions. We also discuss stability and present the complete PML model for the 3-D velocity-stress system.

https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2002JB002235

On the implementation of perfectly matched layers in a three‐dimensional fourth‐order velocity‐stress finite difference scheme

Method for simultaneous full-wavefield inversion of gathers of source (or receiver) encoded geophysical data to determine a physical properties model for a subsurface region, especially suitable for surveys where fixed receiver geometry conditions were not satisfied in the data acquisition. First, a shallow time window of the data (202) where the fixed receiver condition is satisfied is inverted by simultaneous encoded (203) source inversion (205). Then, the deeper time window of the data (208) is inverted by sparse sequential source inversion (209), using the physical properties model from the shallow time window (206) as a starting model (207). Alternatively, the shallow time window model is used to simulate missing far offset data (211) producing a data set satisfying the stationary receiver assumption, after which this data set is source encoded (212) and inverted by simultaneous source inversion (214).

Hybrid method for full waveform inversion using simultaneous and sequential source method

Method for simultaneous full-wavefield inversion of gathers of source (or receiver) encoded geophysical data to determine a physical properties model (118) for a subsurface region, especially suitable for surveys where fixed receiver geometry conditions were not satisfied in the data acquisition. Simultaneous source separation (104) is performed to lessen any effect of the measured geophysical data's not satisfying the fixed-receiver assumption. A data processing step (106) coming after the simultaneous source separation acts to conform model-simulated data (105) to the measured geophysical data (108) for source and receiver combinations that are missing in the measured geophysical data.

Simultaneous source encoding and source separation as a practical solution for full wavefield inversion

Summary We developed a seismic-attenuation-tomography algorithm for improving image amplitude fidelity and other geophysical applications. This algorithm is based upon the centroid-frequency-shift method. However, compared with the conventional centroid-frequency-shift method, our algorithm has been significantly improved through the

A Robust And Accurate Seismic Attenuation Tomography Algorithm

Method for reconstructing subsurface profiles for seismic velocity or other geophysical properties from recorded seismic data. In one embodiment, a starting model of seismic velocity is assumed ( 10 ). The computational domain is divided into two (or more) subdomains by horizontal planes based on an analysis of velocity model ( 30 ), and the allowed maximum grid size for each subdomain is determined ( 50 ). Auxiliary perfectly matched layers (PML's) are attached to each planar interface between subdomains ( 80 ), e.g. two PML's on each side of the interface between the coarse and fine subdomains. Simulated seismic data are computed using the SG-DO technique ( 100 - 230 ). The simulated seismic data are compared to the recorded seismic data, then the residual is calculated ( 240 ) and used to update the model ( 320 ). The method may be iterated until the model is suitably converged ( 260 ).

Carey Marcinkovich

Papers

On the implementation of perfectly matched layers in a three‐dimensional fourth‐order velocity‐stress finite difference scheme

Hybrid method for full waveform inversion using simultaneous and sequential source method

Simultaneous source encoding and source separation as a practical solution for full wavefield inversion

A Robust And Accurate Seismic Attenuation Tomography Algorithm

Full Waveform Inversion Using Perfectly Reflectionless Subgridding