Frequency-domain Elastic Full Waveform Inversion Research Articles

Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? The Gorgon data FWI case study

Frequency-domain full-waveform inversion (FWI) is potentially amenable to efficient processing of full-azimuth long-offset stationary-recording seabed acquisition carried out with a sparse layout of ocean-bottom nodes (OBNs) and broadband sources because the inversion can be performed with a few discrete frequencies. However, computing the solution of the forward (boundary-value) problem efficiently in the frequency domain with linear algebra solvers remains a challenge for large computational domains involving tens to hundreds of millions of parameters. We illustrate the feasibility of 3D frequency-domain FWI with a subset of the 2015/2016 Gorgon OBN data set in the North West Shelf, Australia. We solve the forward problem with the massively parallel multifrontal direct solver MUMPS, which includes four key features to reach high computational efficiency: an efficient parallelism combining message-passing interface and multithreading, block low-rank compression, mixed-precision arithmetic, and efficient processing of sparse sources. The Gorgon subdata set involves 650 OBNs that are processed as reciprocal sources and 400,000 sources. Monoparameter FWI for vertical wavespeed is performed in the viscoacoustic vertically transverse isotropic approximation with a classical frequency continuation approach proceeding from a starting frequency of 1.7 Hz to a final frequency of 13 Hz. The target covers an area ranging from 260 km2 (frequency ≥ 8.5 Hz) to 705 km2 (frequency ≤ 8.5 Hz) for a maximum depth of 8 km. Compared to the starting model, FWI dramatically improves the reconstruction of the bounding faults of the Gorgon horst at reservoir depths as well as several intrahorst faults and several horizons of the Mungaroo Formation down to a depth of 7 km. Seismic modeling reveals a good kinematic agreement between recorded and simulated data, but amplitude mismatches between the recorded and simulated reflection from the reservoir suggest elastic effects. Therefore, future works involve multiparameter reconstruction for density and attenuation before considering elastic FWI from hydrophone and geophone data.

Frequency-domain elastic full-waveform inversion based on Green functions

In this paper, the elastic full-waveform inversion in the frequency domain based on Green functions is investigated. The method allows to image media velocities of compressional wave and shear wave by using vertical and horizontal components. The forward problem is solved by the finite difference scheme with the perfectly matched layer in the frequency domain. The inversion is an optimization iterative process to minimize the residual between the synthetic data and the observed data on the surface. A gradient method is used to find the optimization direction, preconditioned with the diagonal part of the inverse Hessian matrix. The gradient is constructed based on Green functions, which serves to compute the gradient of the objective function with respect to model parameters effectively. The inversion is completed from low frequency to high frequency so that smaller wavelengths are progressively introduced in inversion. Numerical computations for the benchmark Overthrust model demonstrate the full-waveform inversion method in this paper can inverse media velocities of complex structures accurately.

Elastic full-waveform inversion of steeply dipping structures with prismatic waves

High-resolution reconstruction of steeply dipping structures is an important but challenging subject in seismic exploration. Prismatic reflections that contain information on these structures are helpful for reconstructing steeply dipping structures. Elastic full-waveform inversion (EFWI) is a powerful tool that can accurately estimate subsurface parameters from multicomponent seismic data, which can provide information useful for characterizing oil and gas reservoirs. We have constructed the relationship between the forward and inverse problems related to the prismatic reflections by considering the multiparameter exact Hessian in realistic elastic media. We numerically analyze the characteristics of the multiparameter exact Hessian and determine that, when prismatic reflections are apparent in multicomponent data, the multiparameter delta Hessian has a strong influence. We develop this in more detail through forward analysis and determine that the multiparameter delta Hessian considers not only the prismatic reflections but also compensates for the primary reflections in multicomponent data. To use the prismatic waves, we develop a migration/demigration approach-based truncated Newton (TN) method in frequency-domain EFWI, whose storage requirements and computational costs are the same as those of the truncated Gauss-Newton (TGN) method. Realistic 2D numerical examples demonstrate that, compared with the TGN method based on the first-order Born approximation, the TN method can converge faster and obtain higher accuracy in the reconstruction of steeply dipping structures.

A numerical study of 3D frequency-domain elastic full-waveform inversion

We have developed a 3D elastic full-waveform inversion (FWI) algorithm with forward modeling and inversion performed in the frequency domain. The Helmholtz equation is solved with a second-order finite-difference method using an iterative solver equipped with an efficient complex-shifted incomplete LU-based preconditioner. The inversion is based on the minimization of the data misfit functional and a total variation regularization for the unknown model parameters. We implement the Gauss-Newton method as the optimization engine for the inversions. The codes are parallelized with a message passing interface based on the number of shots and receivers. We examine the performance of this elastic FWI algorithm and workflow on synthetic examples including surface seismic and vertical seismic profile configurations. With various initial models, we manage to obtain high-quality velocity images for 3D earth models.

Depth Scaling Strategy Using a Flexible Damping Factor for Frequency-Domain Elastic Full Waveform Inversion

We introduce a depth scaling strategy to improve the accuracy of frequency-domain elastic full waveform inversion (FWI) using the new pseudo-Hessian matrix for seismic data without low frequency components. The depth scaling strategy is based on the fact that the damping factor in the Levenberg-Marquardt method controls the energy concentration in the gradient. In other words, a large damping factor makes the Levenberg Marquardt method similar to the steepest-descent method, by which shallow structures are mainly recovered. With a small damping factor, the Levenberg-Marquardt method becomes similar to the Gauss-Newton methods by which we can resolve deep structures as well as shallow structures. In our depth scaling strategy, a large damping factor is used in the early stage and then decreases automatically with the trend of error as the iteration goes on. With the depth scaling strategy, we can gradually move the parameter-searching region from shallow to deep parts. This flexible damping factor plays a role in retarding the model parameter update for shallow parts and mainly inverting deeper parts in the later stage of inversion. By doing so, we can improve deep parts in inversion results. The depth scaling strategy is applied to synthetic data without lowfrequency components for a modified version of the SEG/EAGE overthrust model. Numerical examples show that the flexible damping factor yields better results than the constant damping factor when reliable low frequency components are missing.

2D frequency-domain elastic full-waveform inversion using the block-diagonal pseudo-Hessian approximation

Elastic full-waveform inversion (EFWI) of multicomponent seismic data is a powerful tool for estimating the subsurface elastic parameters with high accuracy. However, the trade-offs between multiple parameters increase the nonlinearity of EFWI. Although the conventional diagonal-approximate Hessian matrix describes the illumination and limited bandwidth effects, it ignores the trade-off effects and decreases the convergence rate of EFWI. We have developed a block-diagonal pseudo-Hessian operator for 2D frequency-domain EFWI to take into account the approximate trade-offs among the P-wave (compressional-wave) velocity, S-wave (shear-wave) velocity, and density without extra computational costs on forward simulations. The Hessian matrix tends toward a block-diagonal matrix as the frequency grows to infinity; thus, the proposed block-diagonal pseudo-Hessian matrix is more accurate at higher frequencies. The inverse of the block-diagonal pseudo-Hessian matrix is used as a preconditioner for the nonlinear conjugate-gradient method to simultaneously reconstruct P- and S-wave velocities and density. This approach effectively mitigates the crosstalk artifacts by correcting the gradients from the trade-off effects and produces more rapid inversion convergence, which becomes more significant at higher frequencies. Synthetic experiments on an inclusion model and the elastic Marmousi2 model demonstrate its feasibility and validity in EFWI.

Fast 3D frequency-domain full-waveform inversion with a parallel block low-rank multifrontal direct solver: Application to OBC data from the North Sea

Wide-azimuth long-offset ocean bottom cable (OBC)/ocean bottom node surveys provide a suitable framework to perform computationally efficient frequency-domain full-waveform inversion (FWI) with a few discrete frequencies. Frequency-domain seismic modeling is performed efficiently with moderate computational resources for a large number of sources with a sparse multifrontal direct solver (Gauss-elimination techniques for sparse matrices). Approximate solutions of the time-harmonic wave equation are computed using a block low-rank (BLR) approximation, leading to a significant reduction in the operation count and in the volume of communication during the lower upper (LU) factorization as well as offering great potential for reduction in the memory demand. Moreover, the sparsity of the seismic source vectors is exploited to speed up the forward elimination step during the computation of the monochromatic wavefields. The relevance and the computational efficiency of the frequency-domain FWI performed in the viscoacoustic vertical transverse isotropic (VTI) approximation was tested with a real 3D OBC case study from the North Sea. The FWI subsurface models indicate a dramatic resolution improvement relative to the initial model built by reflection traveltime tomography. The amplitude errors introduced in the modeled wavefields by the BLR approximation for different low-rank thresholds have a negligible footprint in the FWI results. With respect to a standard multifrontal sparse direct factorization, and without compromise of the accuracy of the imaging, the BLR approximation can bring a reduction of the LU factor size by a factor of up to three. This reduction is not yet exploited to reduce the effective memory usage (ongoing work). The flop reduction can be larger than a factor of 10 and can bring a factor of time reduction of around three. Moreover, this reduction factor tends to increase with frequency, namely with the matrix size. Frequency-domain viscoacoustic VTI FWI can be viewed as an efficient tool to build an initial model for elastic FWI of 4C OBC data.

Weighted pseudo-Hessian for frequency-domain elastic full waveform inversion

In full waveform inversion, an appropriate preconditioner for scaling the gradient improves the inversion results and accelerates the convergence speed. The Newton-based method uses the Hessian to scale the gradient and ensures a correctly inverted image with good convergence. However, calculating the full Hessian or approximated Hessian and obtaining their inverse matrices require extensive computation, which has remained an obstacle to the application of Newton-based methods to full waveform inversion. Many attempts have been made to reduce the computational cost of obtaining the Hessian and its inverse; among these alternatives, the pseudo-Hessian method has been widely used. The use of the pseudo-Hessian reduces the computational cost by regarding the zero-lag correlation of the impulse response as a unit matrix, but the pseudo-Hessian cannot properly scale the deep portion of a model. Therefore, we proposed a weighted pseudo-Hessian which can overcome the limitations of the conventional pseudo-Hessian. The weighted pseudo-Hessian was generated by combining a weighting matrix with the conventional pseudo-Hessian. The weighting matrix is an amplitude field and helps the pseudo-Hessian to scale the gradient properly. Therefore, the weighted pseudo-Hessian can effectively scale the gradient from the shallow part to the deep part of the model with great balance. Calculating the weighting function required minimal computation, such that the computational cost of generating the weighted pseudo-Hessian was nearly the same as the computational cost needed to calculate the conventional pseudo-Hessian. To verify the proposed algorithm, Marmousi-2 data were used for the synthetic test. The results indicate that the weighted pseudo-Hessian can effectively scale the gradient from the shallow to deep portions of a model.

2D frequency-domain elastic full-waveform inversion using time-domain modeling and a multistep-length gradient approach

To decouple the parameters in elastic full-waveform inversion (FWI), we evaluated a new multistep-length gradient approach to assign individual weights separately for each parameter gradient and search for an optimal step length along the composite gradient direction. To perform wavefield extrapolations for the inversion, we used parallelized high-precision finite-element (FE) modeling in the time domain. The inversion was implemented in the frequency domain; the data were obtained at every subsurface grid point using the discrete Fourier transform at each time-domain extrapolation step. We also used frequency selection to reduce cycle skipping, time windowing to remove the artifacts associated with different source spatial patterns between the test and predicted data, and source wavelet estimation at the receivers over the full frequency spectrum by using a fast Fourier transform. In the inversion, the velocity and density reconstructions behaved differently; as a low-wavenumber tomography (for velocities) and as a high-wavenumber migration (for density). Because velocities and density were coupled to some extent, variations were usually underestimated (smoothed) for [Formula: see text] and [Formula: see text] and correspondingly overestimated (sharpened) for [Formula: see text]. The impedances [Formula: see text] and [Formula: see text] from the products of the velocity and density results compensated for the under- or overestimations of their variations, so the recovered impedances were closer to the correct ones than [Formula: see text], [Formula: see text], and [Formula: see text] were separately. Simultaneous reconstruction of [Formula: see text], [Formula: see text], and [Formula: see text] was robust on the FE and finite-difference synthetic data (without surface waves) from the elastic Marmousi-2 model; satisfactory results are obtained for [Formula: see text], [Formula: see text], [Formula: see text], and the recovered [Formula: see text] and [Formula: see text] from their products. Convergence is fast, needing only a few tens of iterations, rather than a few hundreds of iterations that are typical in most other elastic FWI algorithms.

Weighting technique using backpropagated wavefields incited by deconvolved residuals for frequency-domain elastic full waveform inversion

Weighting technique using backpropagated wavefields incited by deconvolved residuals for frequency-domain elastic full waveform inversion

Frequency-domain elastic full waveform inversion for VTI media

SUMMARY To describe subsurface structures in anisotropic media properly, particularly in transversely isotropic media with a vertical symmetry axis (VTI), which frequently appear in sedimentary basin environments, we develop a frequency-domain elastic full waveform inversion algorithm for 2-D VTI media. The inversion algorithm is based on the cell-based finite-difference modelling method and the adjoint state of the wave equation. Because the anisotropic inversion for VTI media deals with more elastic constants than the isotropic inversion, it is more prone to obtain local minimum solutions. For this reason, we may not succeed in properly describing the elastic constants of subsurface media if we only apply the standard inversion techniques to anisotropic waveform inversion. To compensate for the ill-posedness of the anisotropic waveform inversion, we couple elastic constants C11 and C33 based on Thomsen's relationship, which is also supported by the sensitivity analysis with respect to the parameters. To enhance the inversion results, we apply the frequency-selection strategy, moving from lower to higher frequencies and we carry out the inversion process over two stages. In both stages, all of the elastic constants are simultaneously optimized, as is done in the conventional waveform inversion. However, we only accept the inversion results for C11, C33 and C44 at the first stage, which will be used as the starting models for the second stage and C13 is reinitialized as a linearly increasing model. We apply our waveform inversion algorithm to a simple 3-layered model and a part of the overthrust model. For the 3-layered model, the first inversion stage is enough to yield reasonable inversion results for all of the elastic constants. For the overthrust model, the second stage is needed to enhance the inversion results for C13.

Which data residual norm for robust elastic frequency-domain full waveform inversion?

Elastic full-waveform inversion is an ill-posed data-fitting procedure that is sensitive to noise, inaccuracies of the starting model, definition of multiparameter classes, and inaccurate modeling of wavefield amplitudes. We have investigated the performance of different minimization functionals as the least-squares norm [Formula: see text], the least-absolute-values norm [Formula: see text], and combinations of both (the Huber and so-called hybrid criteria) with reference to two noisy offshore (Valhall model) and onshore (overthrust model) synthetic data sets. The four minimization functionals were implemented in 2D elastic frequency-domain full-waveform inversion (FWI), where efficient multiscale strategies were designed by successive inversions of a few increasing frequencies. For the offshore and onshore case studies, the [Formula: see text]-norm provided the most reliable models for P- and S-wave velocities ([Formula: see text] and[Formula: see text]), even when strongly decimated data sets that correspond to few frequencies were used in the inversion and when outliers polluted the data. The [Formula: see text]-norm can provide reliable results in the presence of uniform white noise for [Formula: see text] and [Formula: see text] if the data redundancy is increased by refining the frequency sampling interval in the inversion at the expense of computational efficiency. The [Formula: see text]-norm and the Huber and hybrid criteria, unlike the [Formula: see text]-norm, allow for successful imaging of the [Formula: see text] model from noisy data in a soft-seabed environment, where the P-to-S-waves have a small footprint in the data. However, the Huber and hybrid criteria are sensitive to a threshold criterion that controls the transition between the criteria and that requires tedious trial-and-error investigations for reliable estimation. The [Formula: see text]-norm provides a robust alternative to the [Formula: see text]-norm for inverting decimated data sets in the framework of efficient frequency-domain FWI.

Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion

Quantitative imaging of the elastic properties of the subsurface at depth is essential for civil engineering applications and oil- and gas-reservoir characterization. A realistic synthetic example provides for an assessment of the potential and limits of 2D elastic full-waveform inversion (FWI) of wide-aperture seismic data for recovering high-resolution P- and S-wave velocity models of complex onshore structures. FWI of land data is challenging because of the increased nonlinearity introduced by free-surface effects such as the propagation of surface waves in the heterogeneous near-surface. Moreover, the short wavelengths of the shear wavefield require an accurate S-wave velocity starting model if low frequencies are unavailable in the data. We evaluated different multiscale strategies with the aim of mitigating the nonlinearities. Massively parallel full-waveform inversion was implemented in the frequency domain. The numerical optimization relies on a limited-memory quasi-Newton algorithm thatoutperforms the more classic preconditioned conjugate-gradient algorithm. The forward problem is based upon a discontinuous Galerkin (DG) method on triangular mesh, which allows accurate modeling of free-surface effects. Sequential inversions of increasing frequencies define the most natural level of hierarchy in multiscale imaging. In the case of land data involving surface waves, the regularization introduced by hierarchical frequency inversions is not enough for adequate convergence of the inversion. A second level of hierarchy implemented with complex-valued frequencies is necessary and provides convergence of the inversion toward acceptable P- and S-wave velocity models. Among the possible strategies for sampling frequencies in the inversion, successive inversions of slightly overlapping frequency groups is the most reliable when compared to the more standard sequential inversion of single frequencies. This suggests that simultaneous inversion of multiple frequencies is critical when considering complex wave phenomena.

Parsimonious finite-volume frequency-domain method for 2-DP-SV-wave modelling

SUMMARY A new numerical technique for solving 2-D elastodynamic equations based on a finite-volume frequency-domain approach is proposed. This method has been developed as a tool to perform 2-D elastic frequency-domain full-waveform inversion. In this context, the system of linear equations that results from the discretization of the elastodynamic equations is solved with a direct solver, allowing efficient multiple-source simulations at the partial expense of the memory requirement. The discretization of the finite-volume approach is through triangles. Only fluxes with the required quantities are shared between the cells, relaxing the meshing conditions, as compared to finite-element methods. The free surface is described along the edges of the triangles, which can have different slopes. By applying a parsimonious strategy, the stress components are eliminated from the discrete equations and only the velocities are left as unknowns in the triangles. Together with the local support of the P0 finite-volume stencil, the parsimonious approach allows the minimizing of core memory requirements for the simulation. Efficient perfectly matched layer absorbing conditions have been designed for damping the waves around the grid. The numerical dispersion of this FV formulation is similar to that of O(Δx2) staggered-grid finite-difference (FD) formulations when considering structured triangular meshes. The validation has been performed with analytical solutions of several canonical problems and with numerical solutions computed with a well-established FD time-domain method in heterogeneous media. In the presence of a free surface, the finite-volume method requires 10 triangles per wavelength for a flat topography, and fifteen triangles per wavelength for more complex shapes, well below the criteria required by the staircase approximation of O(Δx2) FD methods. Comparisons between the frequency-domain finite-volume and the O(Δx2) rotated FD methods also show that the former is faster and less memory demanding for a given accuracy level, an attractive feature for frequency-domain seismic inversion. We have thus developed an efficient method for 2-D P–SV-wave modelling on structured triangular meshes as a tool for frequency-domain full-waveform inversion. Further work is required to improve the accuracy of the method on unstructured meshes.

Frequency-Domain Elastic Full Waveform Inversion Using the New Pseudo-Hessian Matrix: Experience of Elastic Marmousi-2 Synthetic Data

Abstract A proper scaling method allows us to find better solutions in waveform inversion, and it can also provide better images in true-amplitude migration methods based on a least-squares method. For scaling the gradient of a misfit function, we define a new pseudo-Hessian matrix by combining the conventional pseudo-Hessian matrix with amplitude fields. Because the conventional pseudo-Hessian matrix is assumed to neglect the zero-lag autocorrelation terms of impulse responses in the approximate Hessian matrix of the Gauss–Newton method, it has certain limitations in scaling the gradient of a misfit function relative to the approximate Hessian matrix. To overcome these limitations, we introduce amplitude fields to the conventional pseudo-Hessian matrix, and the new pseudo-Hessian matrix is applied to the frequency-domain elastic full waveform inversion. This waveform inversion algorithm follows the conventional procedures of waveform inversion using the backpropagation algorithm. A conjugate-gradient method is employed to derive an optimized search direction, and a backpropagation algorithm is used to calculate the gradient of the misfit function. The source wavelet is also estimated simultaneously with elastic parameters. The new pseudo-Hessian matrix can be calculated without the extra computational costs required by the conventional pseudo-Hessian matrix, because the amplitude fields can be readily extracted from forward modeling. Synthetic experiments show that the new pseudo-Hessian matrix provides better results than the conventional pseudo-Hessian matrix, and thus, we believe that the new pseudo-Hessian matrix is an alternative to the approximate Hessian matrix of the Gauss–Newton method in waveform inversion.