skip to main content

Title: Accelerating full waveform inversion via source stacking and cross-correlations

Accurate synthetic seismic wavefields can now be computed in 3-D earth models using the spectral element method (SEM), which helps improve resolution in full waveform global tomography. However, computational costs are still a challenge. These costs can be reduced by implementing a source stacking method, in which multiple earthquake sources are simultaneously triggered in only one teleseismic SEM simulation. One drawback of this approach is the perceived loss of resolution at depth, in particular because high-amplitude fundamental mode surface waves dominate the summed waveforms, without the possibility of windowing and weighting as in conventional waveform tomography.

This can be addressed by redefining the cost-function and computing the cross-correlation wavefield between pairs of stations before each inversion iteration. While the Green’s function between the two stations is not reconstructed as well as in the case of ambient noise tomography, where sources are distributed more uniformly around the globe, this is not a drawback, since the same processing is applied to the 3-D synthetics and to the data, and the source parameters are known to a good approximation. By doing so, we can separate time windows with large energy arrivals corresponding to fundamental mode surface waves. This opens the possibility of more » designing a weighting scheme to bring out the contribution of overtones and body waves. It also makes it possible to balance the contributions of frequently sampled paths versus rarely sampled ones, as in more conventional tomography.

Here we present the results of proof of concept testing of such an approach for a synthetic 3-component long period waveform data set (periods longer than 60 s), computed for 273 globally distributed events in a simple toy 3-D radially anisotropic upper mantle model which contains shear wave anomalies at different scales. We compare the results of inversion of 10 000 s long stacked time-series, starting from a 1-D model, using source stacked waveforms and station-pair cross-correlations of these stacked waveforms in the definition of the cost function. We compute the gradient and the Hessian using normal mode perturbation theory, which avoids the problem of cross-talk encountered when forming the gradient using an adjoint approach. We perform inversions with and without realistic noise added and show that the model can be recovered equally well using one or the other cost function.

The proposed approach is computationally very efficient. While application to more realistic synthetic data sets is beyond the scope of this paper, as well as to real data, since that requires additional steps to account for such issues as missing data, we illustrate how this methodology can help inform first order questions such as model resolution in the presence of noise, and trade-offs between different physical parameters (anisotropy, attenuation, crustal structure, etc.) that would be computationally very costly to address adequately, when using conventional full waveform tomography based on single-event wavefield computations.

« less
 ;  ;  
Award ID(s):
Publication Date:
Journal Name:
Geophysical Journal International
Page Range or eLocation-ID:
p. 308-322
Oxford University Press
Sponsoring Org:
National Science Foundation
More Like this

    In previous publications, we presented a general framework, which we called ‘box tomography’, that allows the coupling of any two different numerical seismic wave propagation solvers, respectively outside and inside a target region, or ‘box’. The goal of such hybrid wavefield computations is to reduce the cost of computations in the context of full-waveform inversion for structure within the target region, when sources and/or receivers are located at large distances from the box. Previously, we had demonstrated this approach with sources and receivers outside the target region in a 2-D acoustic spherical earth model, and demonstrated and applied this methodology in the 3-D spherical elastic Earth in a continental scale inversion in which all stations were inside the target region. Here we extend the implementation of the approach to the case of a 3-D global elastic earth model in the case where both sources and stations are outside the box. We couple a global 3-D solver, SPECFEM3D_GLOBE, for the computation of the wavefield and Green’s functions in a reference 3-D model, with a regional 3-D solver, RegSEM, for the computation of the wavefield within the box, by means of time-reversal mirrors. We briefly review key theoretical aspects, showing inmore »particular how only the displacement is needed to be stored at the boundary of the box. We provide details of the practical implementation, including the geometrical design of the mirrors, how we deal with different sizes of meshes in the two solvers, and how we address memory-saving through the use of B-spline compression of the recorded wavefield on the mirror. The proposed approach is numerically efficient but also versatile, since adapting it to other solvers is straightforward and does not require any changes in the solver codes themselves, as long as the displacement can be recovered at any point in time and space. We present benchmarks of the hybrid computations against direct computations of the wavefield between a source and an array of stations in a realistic geometry centred in the Yellowstone region, with and without a hypothetical plume within the ‘box’, and with a 1-D or a 3-D background model, down to a period of 20 s. The ultimate goal of this development is for applications in the context of imaging of remote target regions in the deep mantle, such as, for example, Ultra Low Velocity Zones.

    « less

    Crustal seismic velocity models provide essential information for many applications including earthquake source properties, simulations of ground motion and related derivative products. We present a systematic workflow for assessing the accuracy of velocity models with full-waveform simulations. The framework is applied to four regional seismic velocity models for southern California: CVM-H15.11, CVM-S4.26, CVM-S4.26.M01 that includes a shallow geotechnical layer, and the model of Berg et al. For each model, we perform 3-D viscoelastic wave propagation simulations for 48 virtual seismic noise sources (down to 2 s) and 44 moderate-magnitude earthquakes (down to 2 s generally and 0.5 s for some cases) assuming a minimum shear wave velocity of 200 m s–1. The synthetic waveforms are compared with observations associated with both earthquake records and noise cross-correlation data sets. We measure, at multiple period bands for well-isolated seismic phases, traveltime delays and normalized zero-lag cross-correlation coefficients between the synthetic and observed data. The obtained measurements are summarized using the mean absolute derivation of time delay and the mean correlation coefficient. These two metrics provide reliable statistical representations of model quality with consistent results in all data sets. In addition to assessing the overall (average) performance of different models in the entire study area, wemore »examine spatial variations of the models’ quality. All examined models show good phase and waveform agreements for surface waves at periods longer than 5 s, and discrepancies at shorter periods reflecting small-scale heterogeneities and near-surface structures. The model performing best overall is CVM-S4.26.M01. The largest misfits for both body and surface waves are in basin structures and around large fault zones. Inaccuracies generated in these areas may affect tomography and model simulation results at other regions. The seismic velocity models for southern California can be improved by adding better resolved structural representations of the shallow crust and volumes around the main faults.

    « less

    Precisely constraining the source parameters of large earthquakes is one of the primary objectives of seismology. However, the quality of the results relies on the quality of synthetic earth response. Although earth structure is laterally heterogeneous, particularly at shallow depth, most earthquake source studies at the global scale rely on the Green's functions calculated with radially symmetric (1-D) earth structure. To avoid the impact of inaccurate Green's functions, these conventional source studies use a limited set of seismic phases, such as long-period seismic waves, broad-band P and S waves in teleseismic distances (30° < ∆ < 90°), and strong ground motion records at close-fault stations. The enriched information embedded in the broad-band seismograms recorded by global and regional networks is largely ignored, limiting the spatiotemporal resolution. Here we calculate 3-D strain Green's functions at 30 GSN stations for source regions of 9 selected global earthquakes and one earthquake-prone area (California), with frequency up to 67 mHz (15 s), using SPECFEM3D_GLOBE and the reciprocity theorem. The 3-D SEM mesh model is composed of mantle model S40RTS, crustal model CRUST2.0 and surface topography ETOPO2. We surround each target event with grids in horizontal spacing of 5 km and vertical spacing of 2.0–3.0 km, allowing usmore »to investigate not only the main shock but also the background seismicity. In total, the response at over 210 000 source points is calculated in simulation. The number of earthquakes, including different focal mechanisms, centroid depth range and tectonic background, could further increase without additional computational cost if they were properly selected to avoid overloading individual CPUs. The storage requirement can be reduced by two orders of magnitude if the output strain Green's functions are stored for periods over 15 s. We quantitatively evaluate the quality of these 3-D synthetic seismograms, which are frequency and phase dependent, for each source region using nearby aftershocks, before using them to constrain the focal mechanisms and slip distribution. Case studies show that using a 3-D earth model significantly improves the waveform similarity, agreement in amplitude and arrival time of seismic phases with the observations. The limitations of current 3-D models are still notable, dependent on seismic phases and frequency range. The 3-D synthetic seismograms cannot well match the high frequency (>40 mHz) S wave and (>20 mHz) Rayleigh wave yet. Though the mean time-shifts are close to zero, the standard deviations are notable. Careful calibration using the records of nearby better located earthquakes is still recommended to take full advantage of better waveform similarity due to the use of 3-D models. Our results indicate that it is now feasible to systematically study global large earthquakes using full 3-D earth response in a global scale.

    « less

    For over 40 yr, the global centroid-moment tensor (GCMT) project has determined location and source parameters for globally recorded earthquakes larger than magnitude 5.0. The GCMT database remains a trusted staple for the geophysical community. Its point-source moment-tensor solutions are the result of inversions that model long-period observed seismic waveforms via normal-mode summation for a 1-D reference earth model, augmented by path corrections to capture 3-D variations in surface wave phase speeds, and to account for crustal structure. While this methodology remains essentially unchanged for the ongoing GCMT catalogue, source inversions based on waveform modelling in low-resolution 3-D earth models have revealed small but persistent biases in the standard modelling approach. Keeping pace with the increased capacity and demands of global tomography requires a revised catalogue of centroid-moment tensors (CMT), automatically and reproducibly computed using Green's functions from a state-of-the-art 3-D earth model. In this paper, we modify the current procedure for the full-waveform inversion of seismic traces for the six moment-tensor parameters, centroid latitude, longitude, depth and centroid time of global earthquakes. We take the GCMT solutions as a point of departure but update them to account for the effects of a heterogeneous earth, using the global 3-Dmore »wave speed model GLAD-M25. We generate synthetic seismograms from Green's functions computed by the spectral-element method in the 3-D model, select observed seismic data and remove their instrument response, process synthetic and observed data, select segments of observed and synthetic data based on similarity, and invert for new model parameters of the earthquake’s centroid location, time and moment tensor. The events in our new, preliminary database containing 9382 global event solutions, called CMT3D for ‘3-D centroid-moment tensors’, are on average 4 km shallower, about 1 s earlier, about 5 per cent larger in scalar moment, and more double-couple in nature than in the GCMT catalogue. We discuss in detail the geographical and statistical distributions of the updated solutions, and place them in the context of earlier work. We plan to disseminate our CMT3D solutions via the online ShakeMovie platform.

    « less
  5. The origins of tectonic structures in East Antarctica, such as the Gamburtsev Subglacial Mountains (GSMs), the Wilkes Subglacial Basin (WSB), the Aurora Subglacial Basin (ASB), and the Transantarctic Mountains (TAMs), are not clearly understood. Previous investigations have proposed multiple origin models to explain the formation of these structures; however, existing tomographic images lack resolution and consistency given the sparse seismic coverage in East Antarctica. We use full-waveform ambient noise tomography to model the shear-wave velocity structure beneath East Antarctica to further investigate these features. We extract Rayleigh-wave Empirical Green’s Func-tions (EGFs) between periods of 15 and 340 secs from ambient seismic noise using a frequency time normalization technique. Synthetic waveforms are simulated through a 3-D heterogenous Earth model with a lateral grid spacing of 0.025º (~2.25 km) using a finite-difference wave propagation method. The synthetic seismograms are cross-correlated with the EGFs to measure the phase delays. The fi-nite-frequency sensitivity kernels are calculated using the scattering-integral approach and the shear-wave velocity model is computed by inverting the phase delays using a sparse damped least-square inversion method. Preliminary results show fast seismic velocities beneath the WSB, which may be associated with thick and stable lithosphere, and slow velocities beneath the ASB, possiblymore »reflecting old rift systems or other inherited tectonic structures. Slow upper mantle velocities are also observed beneath the TAMs, possibly associated with a thermal load that contributes to the uplift of the moun-tain range. Slow shear-wave velocities in the vicinity of the GSMs may be associated with rifting along the extended Lambert Rift System. Our final tomographic model and associated tectonic inter-pretations will be shared.« less