Misfit Gradient Research Articles

The connection of velocity and impedance sensitivity kernels with scattering-angle filtering and its application in full waveform inversion

Multi-scale strategies such as starting from the low-frequency and early-arrival part of recorded data are commonly used in full waveform inversion (FWI) to maneuver complex nonlinearity. An alternative way is to apply appropriate filtering and conditioning to the misfit gradient in the model domain. In acoustic constant-density media, we prove that velocity and impedance sensitivity kernels are equivalent to applying a high-pass and a low-pass scattering-angle filter to a conventional single-parameter velocity (CSV) kernel. The high-pass scattering-angle filter allows the velocity kernel to include low-wavenumber updates (tomography component). In contrast, the low-pass scattering-angle filter helps the impedance kernel to yield high-wavenumber updates (migration component). The velocity model can be updated using a hybrid gradient of two components combined with appropriate weights. This FWI scheme is able to overcome the potential nonlinearity and partially mitigate the cycle-skipping problem. Numerical examples for the SEG/EAGE overthrust model and the Marmousi model demonstrate that the hybrid gradient facilitates FWI to converge faster to the true model even in cases when conventional CSV-based FWI fails.

3-D multiparameter full-waveform inversion for ocean-bottom seismic data using an efficient fluid–solid coupled spectral-element solver

SUMMARYOcean-bottom seismic acquisition systems deployed on the seabed give access to three-component geophone data and hydrophone data. Compared with conventional streamer acquisitions, the separation of sources and receivers makes it possible to increase the maximum offset and azimuth coverage for improving the illumination at depth. Furthermore, the three-component geophones naturally capture elastic wave propagation effects. While this information is mostly overlooked up to now, reconstructing jointly P- and S-wave velocities would significantly improve the subsurface characterization. To achieve a 3-D high-resolution multiparameter reconstruction, we design an efficient 3-D fluid–solid coupled full waveform modelling and inversion engine. In this engine, fluid and solid domains are divided explicitly and handled with the acoustic and elastic wave equations, respectively. The numerical implementation is based on a time-domain spectral-element method (SEM) with a flexible 3-D Cartesian-based hexahedral mesh, which contributes to an accurate coupling of the acoustic and elastic wave equations and high computational efficiency through domain-decomposition based parallelization. We select the best acoustic–elastic coupled formulations among 4 possibilities with criteria based on numerical accuracy and implementation efficiency. Moreover, we propose a specific hybrid approach for the misfit gradient building so as to use a similar modelling solver for both forward and adjoint simulations. Synthetic case studies on a 3-D extended Marmousi-II model and a 3-D deep-water crustal-scale model illustrate how our modelling and inversion engine can efficiently extract information from ocean-bottom seismic data to simultaneously reconstruct both P- and S-wave velocities within a full waveform inversion framework.

Approximating the Gauss–Newton Hessian Using a Space-Wavenumber Filter and its Applications in Least-Squares Seismic Imaging

The acquisition footprint, finite-frequency source, and unbalanced subsurface illumination make it difficult for traditional adjoint-based migration to produce a high-quality image for complex structures. By fitting reflection events with linearized simulation data, least-squares migration (LSM) can iteratively incorporate the effects of the Gauss–Newton Hessian (GNH) to produce high-quality depth profiles. However, high computational costs of forward and adjoint simulations limit the LSM applications in production. In this study, we present an efficient approximation approach for the GNH and utilize it as a preconditioner for the misfit gradient of the LSM to accelerate its convergence. The analytic solution of the GNH in homogeneous media reveals that the columns of the GNH are local spatial functions. Based on this observation, we design a space-wavenumber filter to approximate the GNH for heterogeneous media, which can be efficiently computed with the S-transform and spectral division. The mixed-domain property of this filter allows it to automatically take the wavenumber dependence of the GNH into account and, therefore, helps to improve spatial resolution. Numerical examples demonstrate that the approximated GNH can considerably improve the image quality and speed up the convergence of the LSM.

Diffraction-angle filtering of gradient for acoustic full-waveform inversion

Seismic full-waveform inversion (FWI) estimates subsurface velocity structures by reducing data misfit between observed and modeled data. Simultaneous matching of transmitted and reflected waves in seismic FWI causes different updates of different wavenumber components of a given model depending on the diffraction angle between incident and diffracted rays. Motivated by the inverse-scattering imaging condition and elastic FWI, we have applied a diffraction-angle filtering technique in acoustic FWI, which enables us to separate transmission and reflection energy in the partial derivative wavefields. Diffraction-angle filtering is applied to the virtual source, which is the model parameter perturbation acting as a source for the partial derivative wavefields. Diffraction-angle filtering consists of two diffraction-angle filters (DAFs), DAF-I and DAF-II. The former is derived from the particle acceleration of the incident wavefields and suppresses energies at either small or large diffraction angles by changing the sign of the weighting factor. DAF-I is exactly identical to the conventional inverse scattering approach. The DAF-II is derived from the artificial shear strain of the incident P-wave and additionally suppresses energies at intermediate diffraction angles. With this mechanism, we can design various types of diffraction-angle filtering to control the updates of wavenumber components of the misfit gradient with respect to the P-wave velocity. Using the synthetic Marmousi-II data and real ocean-bottom seismic data from the North Sea, we determined that diffraction-angle filtering enables us to control low-, intermediate-, and high-wavenumber components of the gradient direction.

True-amplitude linearized waveform inversion with the quasi-elastic wave equation

We present a quasi-elastic wave equation as a function of the pressure variable, which can accurately model PP reflections with elastic amplitude variation with offset effects under the first-order Born approximation. The kinematic part of the quasi-elastic wave equation accurately models the propagation of P waves, whereas the virtual-source part, which models the amplitudes of reflections, is a function of the perturbations of density and Lamé parameters [Formula: see text] and [Formula: see text]. The quasi-elastic wave equation generates a scattering radiation pattern that is exactly the same as that for the elastic wave equation, and only requires the solution of two acoustic wave equations for each shot gather. This means that the quasi-elastic wave equation can be used for true-amplitude linearized waveform inversion (also known as least-squares reverse time migration) of elastic PP reflections, in which the corresponding misfit gradients are with respect to the perturbations of density and the P- and S-wave impedances. The perturbations of elastic parameters are iteratively updated by minimizing the [Formula: see text]-norm of the difference between the recorded PP reflections and the predicted pressure data modeled from the quasi-elastic wave equation. Numerical tests on synthetic and field data indicate that true-amplitude linearized waveform inversion using the quasi-elastic wave equation can account for the elastic PP amplitudes and provide a robust estimate of the perturbations of P- and S-wave impedances and, in some cases, the density. In addition, true-amplitude linearized waveform inversion provides images with a wider bandwidth and fewer artifacts because the PP amplitudes are accurately explained. We also determine the 2D scalar quasi-elastic wave equation for P-SV reflections and the 3D vector equation for PS reflections.

Elastic least-squares reverse time migration in vertical transverse isotropic media

Using adjoint-based elastic reverse time migration, it is difficult to produce high-quality reflectivity images due to the limited acquisition apertures, band-limited source time function, and irregular subsurface illumination. Through iteratively computing the Hessian inverse, least-squares migration enables us to reduce the point-spread-function effects and improve the image resolution and amplitude fidelity. By incorporating anisotropy in the 2D elastic wave equation, we have developed an elastic least-squares reverse time migration (LSRTM) method for multicomponent data from the vertically transversely isotropic (VTI) media. Using the perturbed stiffness parameters [Formula: see text] and [Formula: see text] as PP and PS reflectivities, we linearize the elastic VTI wave equation and obtain a Born modeling (demigration) operator. Then, we use the Lagrange multiplier method to derive the corresponding adjoint wave equation and reflectivity kernels. With linearized forward modeling and adjoint migration operators, we solve a linear inverse problem to estimate the subsurface reflectivity models for [Formula: see text] and [Formula: see text]. To reduce the artifacts caused by data over-fitting, we introduce total-variation regularization into the reflectivity inversion, which promotes a sparse solution in terms of the model derivatives. To accelerate the convergence of LSRTM, we use source illumination to approximate the diagonal Hessian and use it as a preconditioner for the misfit gradient. Numerical examples help us determine that our elastic VTI LSRTM method can improve the spatial resolution and amplitude fidelity in comparison to adjoint migration.

Efficient time-domain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible Cartesian-based mesh

Viscoelastic full-waveform inversion is recognized as a challenging task for current acquisition deployment at the crustal scale. We have developed an efficient formulation based on a time-domain spectral-element method on a flexible Cartesian-based mesh. We consider anisotropic elastic coefficients and isotropic attenuation. Complete gradient expressions including the attenuation contribution spread into those of elastic components are given in a consistent way. The influence of attenuation on the P-wave velocity reconstruction is illustrated through a toy configuration. The numerical implementation of the forward problem includes efficient matrix-vector products for solving second-order elastodynamic equations for 3D geometries: An original high-order integration for topography effects is performed at nearly no extra cost. Combined adjoint and forward field recomputation from the final state and previously saved boundary values allows the estimation of misfit gradients for density, elastic parameters, and attenuation factors with no I/O efforts. Two-level parallelism is implemented over the sources and domain decomposition, which is necessary for a realistic 3D configuration. The misfit gradient preconditioning is performed by a so-called Bessel filter using an efficient differential implementation based on finite-element ingredients on the forward mesh instead of the often-used, costly convolution approach. A 3D synthetic illustration is provided on a subset ([Formula: see text]) of the SEG Advanced Modeling (SEAM) Phase II Foothills model with 4 lines of 20 sources. The structurally based Bessel filter and a simple data hierarchy strategy considering early body waves before all waves including surface waves allow a precise reconstruction of the P- and S-wavespeeds while keeping a smooth density description.

Imaging disturbance zones ahead of a tunnel by elastic full-waveform inversion: Adjoint gradient based inversion vs. parameter space reduction using a level-set method

We present and compare two flexible and effective methodologies to predict disturbance zones ahead of underground tunnels by using elastic full-waveform inversion. One methodology uses a linearized, iterative approach based on misfit gradients computed with the adjoint method while the other uses iterative, gradient-free unscented Kalman filtering in conjunction with a level-set representation. Whereas the former does not involve a priori assumptions on the distribution of elastic properties ahead of the tunnel, the latter introduces a massive reduction in the number of explicit model parameters to be inverted for by focusing on the geometric form of potential disturbances and their average elastic properties. Both imaging methodologies are validated through successful reconstructions of simple disturbances. As an application, we consider an elastic multiple disturbance scenario. By using identical synthetic time-domain seismograms as test data, we obtain satisfactory, albeit different, reconstruction results from the two inversion methodologies. The computational costs of both approaches are of the same order of magnitude, with the gradient-based approach showing a slight advantage. The model parameter space reduction approach compensates for this by additionally providing a posteriori estimates of model parameter uncertainty.

Element-by-element parallel spectral-element methods for 3-D teleseismic wave modeling

Abstract. The development of an efficient algorithm for teleseismic wave field modeling is valuable for calculating the gradients of the misfit function (termed misfit gradients) or Fréchet derivatives when the teleseismic waveform is used for adjoint tomography. Here, we introduce an element-by-element parallel spectral-element method (EBE-SEM) for the efficient modeling of teleseismic wave field propagation in a reduced geology model. Under the plane-wave assumption, the frequency–wavenumber (FK) technique is implemented to compute the boundary wave field used to construct the boundary condition of the teleseismic wave incidence. To reduce the memory required for the storage of the boundary wave field for the incidence boundary condition, a strategy is introduced to efficiently store the boundary wave field on the model boundary. The perfectly matched layers absorbing boundary condition (PML ABC) is formulated using the EBE-SEM to absorb the scattered wave field from the model interior. The misfit gradient can easily be constructed in each time step during the calculation of the adjoint wave field. Three synthetic examples demonstrate the validity of the EBE-SEM for use in teleseismic wave field modeling and the misfit gradient calculation.

Fast calculations of the gradient and the Hessian in the time-domain global electromagnetic induction inverse problem

The deterministic approach to the inverse problem of the time-domain electromagnetic induction in a spherical Earth requires the calculation of the first derivative of a misfit function in every step of the minimization process. In addition, an a-posteriori error analysis can benefit from the knowledge of the Hessian, the matrix of the second derivatives. We present the derivation of the formulas for the fast calculation of the misfit gradient, the Hessian, and the Hessian-vector product, based on the solution of an adjoint problem. We validate our results on a synthetic model against a slow finite-difference scheme.

Wave equation dispersion inversion using a difference approximation to the dispersion-curve misfit gradient

We present a surface-wave inversion method that inverts for the S-wave velocity from the Rayleigh wave dispersion curve using a difference approximation to the gradient of the misfit function. We call this wave equation inversion of skeletonized surface waves because the skeletonized dispersion curve for the fundamental-mode Rayleigh wave is inverted using finite-difference solutions to the multi-dimensional elastic wave equation. The best match between the predicted and observed dispersion curves provides the optimal S-wave velocity model. Our method can invert for lateral velocity variations and also can mitigate the local minimum problem in full waveform inversion with a reasonable computation cost for simple models. Results with synthetic and field data illustrate the benefits and limitations of this method.

Efficient 3D elastic full-waveform inversion using wavefield reconstruction methods

In reverse time migration (RTM) or full-waveform inversion (FWI), forward and reverse time propagating wavefields are crosscorrelated in time to form either the image condition in RTM or the misfit gradient in FWI. The crosscorrelation condition requires both fields to be available at the same time instants. For large-scale 3D problems, it is not possible, in practice, to store snapshots of the wavefields during forward modeling due to extreme storage requirements. We have developed an approximate wavefield reconstruction method that uses particle velocity field recordings on the boundaries to reconstruct the forward wavefields during the computation of the reverse time wavefields. The method is computationally effective and requires less storage than similar methods. We have compared the reconstruction method to a boundary reconstruction method that uses particle velocity and stress fields at the boundaries and to the optimal checkpointing method. We have tested the methods on a 2D vertical transversely isotropic model and a large-scale 3D elastic FWI problem. Our results revealed that there are small differences in the results for the three methods.

Choice of regularization in adjoint tomography based on two-dimensional synthetic tests

We present synthetic tests of 2-D adjoint tomography of surface wave traveltimes obtained by the ambient noise cross-correlation analysis across the Czech Republic. The data coverage may be considered perfect for tomography due to the density of the station distribution. Nevertheless, artefacts in the inferred velocity models arising from the data noise may be still observed when weak regularization (Gaussian smoothing of the misfit gradient) or too many iterations are considered. To examine the effect of the regularization and iteration number on the performance of the tomography in more detail we performed extensive synthetic tests. Instead of the typically used (although criticized) checkerboard test, we propose to carry out the tests with two different target models—simple smooth and complex realistic models. The first test reveals the sensitivity of the result on the data noise, while the second helps to analyse the resolving power of the data set. For various noise and Gaussian smoothing levels, we analysed the convergence towards (or divergence from) the target model with increasing number of iterations. Based on the tests we identified the optimal regularization, which we then employed in the inversion of 16 and 20 s Love-wave group traveltimes.

Inversion of multicomponent seismic time shifts for reservoir pressure and length: a feasibility study

Pressure drops associated with reservoir production generate excess stress and strain that cause travel-time shifts of reflected waves. Here,we invert time shifts of P-, S-, and PS-waves measured between baseline and monitor surveys for pressure reduction and reservoir length. The inversion results can be used to estimate compaction-induced stress and strain changes around the reservoir. We implement a hybrid inversion algorithm that incorporates elements of gradient, global/genetic, and nearest neighbour methods and permits exploration of the parameter space while simultaneously following local misfit gradients. Our synthetic examples indicate that optimal estimates of reservoir pressure from P-wave data can be obtained using the reflections from the reservoir top. For S-waves, time shifts from the top of the reservoir can be accurately inverted for pressure if the noise level is low. However, if noise contamination is significant, it is preferable to use S-wave data (or combined shifts of all three modes) from reflectors beneath the reservoir. Joint wave type inversions demonstrate improvements over any single pure mode. Reservoir length can be estimated using the time shifts of any mode from the reservoir top or deeper reflectors. We also evaluate the differences between the actual strain field and those corresponding to the best-case inversion results obtained using P- and S-wave data. Another series of tests addresses the inversion of the time shifts for the pressure drops in two-compartment reservoirs, as well as for the associated strain field. Numerical testing shows that a potentially serious source of error in the inversion is a distortion in the strain-sensitivity coefficients, which govern the magnitude of stiffness changes. This feasibility study suggests which wave types and reflector locations may provide the most accurate estimates of reservoir parameters from compaction-induced time shifts.

FULL WAVEFORM INVERSION FOR TIME-DISTANCE HELIOSEISMOLOGY

Inferring interior properties of the Sun from photospheric measurements of the seismic wavefield constitutes the helioseismic inverse problem. Deviations in seismic measurements (such as wave travel times) from their fiducial values estimated for a given model of the solar interior imply that the model is inaccurate. Contemporary inversions in local helioseismology assume that properties of the solar interior are linearly related to measured travel-time deviations. It is widely known, however, that this assumption is invalid for sunspots and active regions, and likely for supergranular flows as well. Here, we introduce nonlinear optimization, executed iteratively, as a means of inverting for the sub-surface structure of large-amplitude perturbations. Defining the penalty functional as the $L_2$ norm of wave travel-time deviations, we compute the the total misfit gradient of this functional with respect to the relevant model parameters %(only sound speed in this case) at each iteration around the corresponding model. The model is successively improved using either steepest descent, conjugate gradient, or quasi-Newton limited-memory BFGS. Performing nonlinear iterations requires privileging pixels (such as those in the near-field of the scatterer), a practice not compliant with the standard assumption of translational invariance. Measurements for these inversions, although similar in principle to those used in time-distance helioseismology, require some retooling. For the sake of simplicity in illustrating the method, we consider a 2-D inverse problem with only a sound-speed perturbation.

Application of weighted early-arrival waveform inversion to shallow land data

Abstract Seismic imaging of deep land targets is usually difficult since the near-surface velocities are not accurately estimated. Recent studies have shown that inverting traces weighted by the energy of the early-arrivals can improve the accuracy of estimating shallow velocities. In this work, it is explained by showing that the associated misfit gradient function tends to be sensitive to the kinetics of wave propagation and insensitive to the dynamics. A synthetic example verifies the theoretical predictions and shows that the effects of noise and unpredicted amplitude variations in the inversion are reduced using this weighted early arrival waveform inversion (WEWI). We also apply this method to a 2D land data set for estimating the near-surface velocity distribution. The reverse time migration images suggest that, compared to the tomogram inverted directly from the early arrival waveforms, the WEWI tomogram provides a more convincing velocity model and more focused reflections in the deeper part of the image.

Determination of the 3-D distribution of electrical conductivity in Earth’s mantle from Swarm satellite data: Frequency domain approach based on inversion of induced coefficients

Mapping the three-dimensional (3-D) electrical conductivity of Earth’s mantle has been identified as one of the primary scientific objectives for the Swarm satellite mission. We present a 3-D frequency domain inversion scheme to recover mantle conductivity from satellite magnetic data. The scheme is based on an inversion of time spectra of internal (induced) spherical harmonic coefficients of the magnetic potential due to magnetospheric sources. Time series of internal and external (inducing) coefficients, whose determination is a prerequisite for this formulation, will be available as a Swarm Level-2 data product. An iterative gradient-type (quasi-Newton) optimization method is chosen to solve our 3-D non-linear inverse problem. In order to make the inversion tractable, we elaborate an adjoint approach for a fast and robust calculation of the data misfit gradient. We verify our approach with synthetic, but realistic time spectra of internal coefficients, obtained by simulating induction due to a realistic magnetospheric source in a 3-D conductivity model of the Earth. In these model studies, both shape and conductivity of a large-scale conductivity anomaly in the mid-mantle are recovered very well. The inversion scheme also shows to be robust with respect to noise and is therefore ready to process Swarm data.

Global 3-D imaging of mantle electrical conductivity based on inversion of observatory C-responses-I. An approach and its verification

SUMMARY We present a novel frequency-domain inverse solution to recover the 3-D electrical conductivity distribution in the mantle. The solution is based on analysis of local C-responses. It exploits an iterative gradient-type method—limited-memory quasi-Newton method—for minimizing the penalty function consisting of data misfit and regularization terms. The integral equation code is used as a forward engine to calculate responses and data misfit gradients during inversion. An adjoint approach is implemented to compute misfit gradients efficiently. Further improvements in computational load come from parallelizing the scheme with respect to frequencies, and from setting the most time-consuming part of the forward calculations—calculation of Green’s tensors—apart from the inversion loop. Convergence, performance, and accuracy of our 3-D inverse solution are demonstrated with a synthetic numerical example. A companion paper applies the strategy set forth here to real data.

Erratum to: ”Fast Calculation of the Sensitivity Matrix for Responses to the Earth’s Conductivity: General Strategy and Examples”

Modern electromagnetic studies of the Earth deal with three-dimensional data sets and their inversion in the class of three-dimensional models. Due to the large scale of the 3D inverse problem, it is usually solved using gradient-type iteration techniques. These techniques require repeated calculations of the gradient of the penalty function (i.e., the sum of the misfit and the regularization term) with respect to the model parameters that describe the spatial distribution of conductivity. Since there are generally no analytical tools for calculating the misfit gradient, its computation requires the application of some numerical techniques. Unfortunately, it is impossible to perform mass calculations of the misfit gradient even using modern computing if one applies conventional numerical differentiation. In reality, in this case the number of calls of the forward problem is proportional to the number of the sought parameters, which is very large in case of a 3D inversion. However, there is a much more efficient method for calculation of the misfit gradient, namely, the adjoint-field technique. At present, this technique is widely used in the numerical schemes of the 3D inversion of electromagnetic and other data. With this technique we can calculate the misfit gradient with the same computational burden as that required for the solution of a single additional forward problem. Although this technique is widely used, we failed to find any comprehensive description of this technique in the literature that would allow one to apply this approach for any particular scheme of sounding, profiling, or transillumination as easily as a general formula is applied for a given case. In the present paper, we obtain general formulas for the quick calculation of the derivatives of the frequency-domain responses and the derivatives of the misfit function with respect to the variations in conductivity. We also show how the derived general formulas are transformed into particular expressions when we consider the response functions of the conventional and generalized magnetotelluric methods. In addition, we demonstrate how to apply our general methodology to the inversion of magnetovariational induction vectors (tippers), the phase tensor, and multipoint transfer functions, namely, the horizontal magnetic and telluric tensors.

General formalism for the efficient calculation of derivatives of EM frequency-domain responses and derivatives of the misfit

Electromagnetic (EM) studies of the Earth have advanced significantly over the past few years. This progress was driven, in particular, by new developments in the methods of 3-D inversion of EM data. Due to the large scale of the 3-D EM inverse problems, iterative gradient-type methods have mostly been employed. In these methods one has to calculate multiple times the gradient of the penalty function—a sum of misfit and regularization terms—with respect to the model parameters. However, even with modern computational capabilities the straightforward calculation of the misfit gradients based on numerical differentiation is extremely time consuming. Much more efficient and elegant way to calculate the gradient of the misfit is provided by the so-called ‘adjoint’ approach. This is now widely used in many 3-D numerical schemes for inverting EM data of different types and origin. It allows the calculation of the misfit gradient for the price of only a few additional forward calculations. In spite of its popularity we did not find in the literature any general description of the approach, which would allow researchers to apply this methodology in a straightforward manner to their scenario of interest. In the paper, we present formalism for the efficient calculation of the derivatives of EM frequency-domain responses and the derivatives of the misfit with respect to variations of 3-D isotropic/anisotropic conductivity. The approach is rather general; it works with single-site responses, multisite responses and responses that include spatial derivatives of EM field. The formalism also allows for various types of parametrization of the 3-D conductivity distribution. Using this methodology one can readily obtain appropriate formulae for the specific sounding methods. To illustrate the concept we provide such formulae for a number of EM techniques: geomagnetic depth sounding (GDS), conventional and generalized magnetotellurics, the magnetovariational method, horizontal gradient sounding (HGS) and a method that combines HGS with GDS. We also show how the developed formalism can be adapted for the inversion of multisite responses—horizontal magnetic and electric tensors.