Stéphane Operto

An overview of full-waveform inversion in exploration geophysics

Full-waveform inversion FWI is a challenging data-fitting procedure based on full-wavefield modeling to extract quantitative information from seismograms. High-resolution imaging at half the propagated wavelength is expected. Recent advances in high-performance computing and multifold/multicomponent wide-aperture and wide-azimuth acquisitions make 3D acoustic FWI feasible today. Key ingredients of FWI are an efficient forward-modeling engine and a local differential approach, in which the gradient and the Hessian operators are efficiently estimated. Local optimization does not, however, prevent convergence of the misfit function toward local minima because of the limited accuracy of the starting model, the lack of low frequencies, the presence of noise, and the approximate modeling of the wave-physics complexity. Different hierarchical multiscale strategiesaredesignedtomitigatethenonlinearityandill-posedness of FWI by incorporating progressively shorter wavelengths in the parameter space. Synthetic and real-data case studies address reconstructing various parameters, from VP and VS velocities to density, anisotropy, and attenuation. This review attempts to illuminate the state of the art of FWI. Crucial jumps, however, remain necessary to make it as popular as migration techniques. The challenges can be categorized as 1 building accurate starting models with automatic procedures and/or recording low frequencies, 2 defining new minimization criteria to mitigate the sensitivity of FWI to amplitude errors and increasing the robustness of FWI when multiple parameter classes are estimated, and 3 improving computational efficiency by data-compression techniquestomake3DelasticFWIfeasible.

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 inversionFWIof wide-aperture seismic data for recovering high-resolution P- and S-wave velocity models of complex onshore structures. FWI of land data is challengingbecauseoftheincreasednonlinearityintroducedbyfreesurface 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 modeliflowfrequenciesareunavailableinthedata.Weevaluated 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 that outperforms the more classic preconditioned conjugate-gradient algorithm. The forward problem is based upon a discontinuous Galerkin DG method on triangular mesh, which allows accuratemodelingoffree-surfaceeffects.Sequentialinversionsofincreasingfrequenciesdefinethemostnaturallevelofhierarchyin 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 complexvalued frequencies is necessary and provides convergence of the inversion toward acceptable P- and S-wave velocity models. Amongthepossiblestrategiesforsamplingfrequenciesintheinversion, 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 consideringcomplexwavephenomena.

3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study

We present a finite-difference frequency-domain method for 3D visco-acoustic wave propagation modeling. In the frequency domain, the underlying numerical problem is the resolution of a large sparse system of linear equations whose right-hand side term is the source. This system is solved with a massively parallel direct solver. We first present an optimal 3D finite-difference stencil for frequency-domain modeling. The method is based on a parsimonious staggered-grid method. Differential operators are discretized with second-order accurate staggered-grid stencils on different rotated coordinate systems to mitigate numerical anisotropy. An antilumped mass strategy is implemented to minimize numerical dispersion. The stencil incorporates 27 grid points and spans two grid intervals. Dispersion analysis shows that four grid points per wavelength provide accurate simulations in the 3D domain. To assess the feasibility of the method for frequency-domain full-waveform inversion, we computed simulations in the 3D SEG/EAGE overthrust model for frequencies 5, 7, and 10 Hz. Results confirm the huge memory requirement of the factorization (several hundred Figabytes) but also the CPU efficiency of the resolution phase (few seconds per shot). Heuristic scalability analysis suggests that the memory complexity of the factorization is O(35N(4)) for a N-3 grid. Our method may provide a suitable tool to perform frequency-domain full-waveform inversion using a large distributed-memory platform. Further investigation is still necessary to assess more quantitatively the respective merits and drawbacks of time- and frequency-domain modeling of wave propagation to perform 3D full-waveform inversion.

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,definitionofmultiparameterclasses,andinaccuratemodeling of wavefield amplitudes. We have investigated the performance of different minimization functionals as the least-squares norm 2, the least-absolute-values norm 1, and combinations of both the Huber and so-called hybrid criteria with reference to twonoisyoffshoreValhallmodelandonshoreoverthrustmodel 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 1-norm provided the most reliable models for P- and S-wave velocities VP and VS, even when strongly decimated data sets that correspond to fewfrequencieswereusedintheinversionandwhenoutlierspolluted the data. The 2-norm can provide reliable results in the presence of uniform white noise for VP and VS if the data redundancyisincreasedbyrefiningthefrequencysamplingintervalin the inversion at the expense of computational efficiency. The 1-norm and the Huber and hybrid criteria, unlike the 2-norm, allowforsuccessfulimagingoftheVSmodelfromnoisydataina 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 1-norm provides a robust alternativetothe2-normforinvertingdecimateddatasetsinthe frameworkofefficientfrequency-domainFWI.

Velocity model building by 3D frequency-domain, full-waveform inversion of wide-aperture seismic data

We assessed 3D frequency-domain FD acoustic fullwaveforminversionFWIdataasatooltodevelophigh-resolution velocity models from low-frequency global-offset data. The inverse problem was posed as a classic leastsquares optimization problem solved with a steepest-descent method. Inversion was applied to a few discrete frequencies, allowing management of a limited subset of the 3D data volume. The forward problem was solved with a finite-difference frequency-domain method based on a massively paralleldirectsolver,allowingefficientmultiple-shotsimulations. The inversion code was fully parallelized for distributedmemory platforms, taking advantage of a domain decomposition of the modeled wavefields performed by the direct solver. After validation on simple synthetic tests, FWI was applied to two targets channel and thrust system of the 3D SEG/EAGEoverthrustmodel,correspondingto3Ddomains of 78.752.25 km and 13.513.54.65 km, respectively. The maximum inverted frequencies are 15 and 7 Hz for the two applications.Amaximum of 30 dual-core biprocessor nodes with 8 GB of shared memory per node were used for the second target. The main structures were imaged successfully at a resolution scale consistent with the inverted frequencies. Our study confirms the feasibility of 3D frequency-domain FWI of global-offset data on large distributed-memory platforms to develop high-resolution velocity models. These high-velocity models may provide accurate macromodelsforwave-equationprestackdepthmigration.

An efficient frequency-domain full waveform inversion method using simultaneous encoded sources

Three-dimensional full waveform inversion (FWI) still suffers from prohibitively high computational costs that arise because of the seismic modeling for multiple sources that is performed at each nonlinear iteration of FWI. Building supershots by assembling several sources allows mitigation of the number of simulations per FWI iteration, although it adds crosstalk artifacts because of interference between the individual sources of the supershots. These artifacts themselves can be reduced by encoding each individual source with a random phase shift during assembling of the sources. The source encoding method is applied to an efficient frequency-domain FWI, in which a limited number of discrete frequencies or coarsely sampled frequency groups are inverted successively following a multiscale approach. Random codes can be regenerated at each FWI iteration or for each frequency of a group during each FWI iteration, to favor the destructive summation of crosstalk artifacts over FWI iterations. Either a limited ...

Fast 2-D ray+Born migration/inversion in complex media

In this paper, we evaluate the capacity of a fast 2-D ray+Born migration/inversion algorithm to recover the true amplitude of the model parameters in 2-D complex media. The method is based on a quasi‐Newtonian linearized inversion of the scattered wavefield. Asymptotic Green’s functions are computed in a smooth reference model with a dynamic ray tracing based on the wavefront construction method. The model is described by velocity perturbations associated with diffractor points. Both the first traveltime and the strongest arrivals can be inverted. The algorithm is implemented with several numerical approximations such as interpolations and aperture limitation around common midpoints to speed the algorithm. Both theoritical and numerical aspects of the algorithm are assessed with three synthetic and real data examples including the 2-D Marmousi example. Comparison between logs extracted from the exact Marmousi perturbation model and the computed images shows that the amplitude of the velocity perturbations...

Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures

Full Waveform Inversion (FWI) is a powerful tool for quantitative seismic imaging from wide-azimuth seismic data. The method is based on the minimization of the misfit between observed and simulated data. This amounts to the resolution of a large-scale nonlinear minimization problem. The inverse Hessian operator plays a crucial role in this reconstruction process. Accounting accurately for the effect of this operator within the minimization scheme should correct for illumination deficits, restore the amplitude of the subsurface parameters, and help to remove artifacts generated by energetic multiple reflections. Conventional preconditioned gradient-based minimization methods only roughly approximate the effect of this operator. We are interested in this study to another class of minimization methods, named as truncated Newton methods. These methods are based on the computation of the model update through a matrix-free conjugate gradient resolution of the Newton linear system. The aim of this study is to present a feasible implementation of this method for the FWI problem, based on a second-order adjoint state formulation for the computation of Hessian-vector products. We compare this method with the nonlinear conjugate gradient and the l-BFGS method within the context of 2D acoustic frequency FWI for the reconstruction of P-wave velocity models. Two test cases are investigated. The first is the synthetic BP 2004 model, representative of the Gulf Of Mexico geology with high velocity contrasts associated with the presence of salt structures. The second is a 2D real data-set from the Valhall oil field in North sea. These tests emphasize the interesting properties of the truncated Newton method regarding conventional optimization methods within the context of FWI.

Deep structure of the southern Kerguelen Plateau (southern Indian Ocean) from ocean bottom seismometer wide-angle seismic data

Wide-angle seismic data collected during the Kerguelen ocean bottom seismometer experiment provide the first images of the deep structure of the southern Kerguelen Plateau and support a new interpretation of the origin of the plateau. Velocity models based on travel time inversions and reflectivity synthetic seismograms show a 22-km-thick crust composed of ∼1.6 km of sedimentary cover, ∼5.3 km of upper crust, ∼11.0 km of lower crust, and a 4- to 6-km-thick reflective zone immediately above Moho. Velocities in the upper crust (from 3.8–4.5 km/s at top to 6.0–6.5 km/s at bottom) are consistent with the basaltic nature of this layer, the top of which was sampled during the Ocean Drilling Program. Velocities in the lower crust increase continuously from 6.60 km/s at the top to 6.90 km/s at 19.5 km depth. The reflective zone at the base of the crust identified by wide-angle reflections is observed only along the NNW-SSE direction. It consists of alternating high- and low-velocity layers with an average velocity of 6.70 km/s in the NNW-SSE direction and ≥6.90 km/s in the perpendicular direction. Strong azimuthal anisotropy is also observed in the upper mantle with velocities of 8.60 and 8.00 km/s, in the NNW-SSE and E-W directions, respectively. The absence of high velocities at the base of the crust that characterizes many large-volume mafic provinces, the reflective lower crust, and anisotropy in upper mantle suggest that the southern Kerguelen Plateau represents a stretched continental fragment overlain by basaltic flows isolated from the Antarctic margin during the early opening of the Indian Ocean.

Shallow-structure characterization by 2D elastic full-waveform inversion

Assessing the effectiveness of elastic full-waveform-inversion (FWI) algorithms when applied to shallow 2D structures in the presence of a complex topography is critically important. By using FWI, we overcome inherent limitations of conventional seismic methods used for near-surface prospecting (acoustic tomography and multichannel spectral analysis of surface waves). The elastic forward problem, formulated in the frequency domain, is based on a mixed finite-element P0P1 discontinuous Galerkin method to ensure accurate modeling of complex topography effects at a reasonable computing cost. The inversion problem uses an FWI algorithm to minimize the misfit between observed and calculated data. Based on results from a numerical experiment performed on a realistic landslide model inspired from the morphostructure of the Super-Sauze earthflow, we analyzed the effect of using a hierarchical preconditioning strategy, based on a simultaneous multifrequency inversion of damped data, to mitigate the strong nonlinearities coming from the surface waves. This strategy is a key point in alleviating the strong near-surface effects and avoiding convergence toward a local minimum. Using a limited-memory quasi-Newton method improved the convergence level. These findings are analogous to recent applications on large-scale domains, although limited sourcereceiver offset ranges, low-frequency content of the source, and domination of surface waves on the signal led to some difficulties. Regarding the impact of data decimation on the inversion results, we have learned that an inversion restricted to the vertical data component can be successful without significant loss in terms of parameter imagery resolution. In our investigations of the effect of increased source spacing, we found that a sampling of 4 m (less than three times the theoretical maximum of one half-wavelength) led to severe aliasing.