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.

/pdf/an-overview-of-full-waveform-inversion-in-exploration-59sei7fzvh.pdf

An overview of full-waveform inversion in exploration geophysics

Seismic waveforms contain much information that is ignored under standard processing schemes; seismic waveform inversion seeks to use the full information content of the recorded wavefield. In this paper I present, apply, and evaluate a frequency-space domain approach to waveform inversion. The method is a local descent algorithm that proceeds from a starting model to refine the model in order to reduce the waveform misfit between observed and model data. The model data are computed using a full-wave equation, viscoacoustic, frequency-domain, finite-difference method. Ray asymptotics are avoided, and higher-order effects such as diffractions and multiple scattering are accounted for automatically. The theory of frequency-domain waveform/wavefield inversion can be expressed compactly using a matrix formalism that uses finite-difference/finite-element frequency-domain modeling equations. Expressions for fast, local descent inversion using back-propagation techniques then follow naturally. Implementation of these methods depends on efficient frequency-domain forward-modeling solutions; these are provided by recent developments in numerical forward modeling. The inversion approach resembles prestack, reverse-time migration but differs in that the problem is formulated in terms of velocity (not reflectivity), and the method is fully iterative. I illustrate the practical application of the frequency-domain waveform inversion approach using tomographic seismic data from a physical scale model. This allows a full evaluation and verification of the method; results with field data are presented in an accompanying paper. Several critical processes contribute to the success of the method: the estimation of a source signature, the matching of amplitudes between real and synthetic data, the selection of a time window, and the selection of suitable sequence of frequencies in the inversion. An initial model for the inversion of the scale model data is provided using standard traveltime tomographic methods, which provide a robust but low-resolution image. Twenty-five iterations of wavefield inversion are applied, using five discrete frequencies at each iteration, moving from low to high frequencies. The final results exhibit the features of the true model at subwavelength scale and account for many of the details of the observed arrivals in the data.

Seismic waveform inversion in the frequency domain; Part 1, Theory and verification in a physical scale model

SUMMARY By specifying a discrete matrix formulation for the frequency^space modelling problem for linear partial diierential equations (‘FDM’ methods), it is possible to derive a matrix formalism for standard iterative non-linear inverse methods, such as the gradient (steepest descent) method, the Gauss^Newton method and the full Newton method We obtain expressions for each of these methods directly from the discrete FDM method, and we refer to this approach as frequency-domain inversion (FDI)The FDI methods are based on simple notions of matrix algebra, but are nevertheless very general The FDI methods only require that the original partial diierential equations can be expressed as a discrete boundary-value problem (that is as a matrix problem) Simple algebraic manipulation of the FDI expressions allows us to compute the gradient of the mis¢t function using only three forward modelling steps (one to compute the residuals, one to backpropagate the residuals, and a ¢nal computation to compute a step length) This result is exactly analogous to earlier backpropagation methods derived using methods of functional analysis for continuous problems Following from the simplicity of this result, we give FDI expressions for the approximate Hessian matrix used in the Gauss^Newton method, and the full Hessian matrix used in the full Newton method In a new development, we show that the additional term in the exact Hessian, ignored in the Gauss^Newton method, can be e⁄ciently computed using a backpropagation approach similar to that used to compute the gradient vector The additional term in the Hessian predicts the degradation of linearized inversions due to the presence of ¢rst-order multiples (such as free-surface multiples in seismic data) Another interpretation is that this term predicts changes in the gradient vector due to second-order non-linear eiects In a numerical test, the Gauss^Newton and full Newton methods prove eiective in helping to solve the di⁄cult non-linear problem of extracting a smooth background velocity model from surface seismic-re£ection data

Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion

Magnetotelluric (MT) data are inverted for smooth 2-D models using an extension of the existing 1-D algorithm, Occam’s inversion. Since an MT data set consists of a finite number of imprecise data, an infinity of solutions to the inverse problem exists. Fitting field or synthetic electromagnetic data as closely as possible results in theoretical models with a maximum amount of roughness, or structure. However, by relaxing the misfit criterion only a small amount, models which are maximally smooth may be generated. Smooth models are less likely to result in overinterpretation of the data and reflect the true resolving power of the MT method. The models are composed of a large number of rectangular prisms, each having a constant conductivity. Apriori information, in the form of boundary locations only or both boundary locations and conductivity, may be included, providing a powerful tool for improving the resolving power of the data. Joint inversion of TE and TM synthetic data generated from known models al...

/pdf/occam-s-inversion-to-generate-smooth-two-dimensional-models-2xddoz7i9r.pdf

Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data

Prestack migration and/or inversion may be implemented in either the time or the frequency domain. In the frequency domain, it is possible to discretize the frequencies with a much larger sampling interval than that dictated by the sampling theorem and still obtain an imaging result that does not suffer from aliasing (wrap around) in the depth domain. The selection of input frequencies can be reduced when a range of offsets is available; this creates a redundancy of information in the wavenumber coverage of the target. In order to optimize the use of this information, we define a new discretization strategy that depends on the maximum effective offset present in the surface seismic survey: the larger the range of offsets, the fewer frequencies are required. The strategy, exact in a homogeneous 1D earth, selects frequencies by making use of the well-known effect of image stretch in normal-moveout (NMO) correction and in migration (usually considered detrimental for the imaging). The strategy is also useful in more general earth models: we apply it to the 2D Marmousi model and recover a continuous range of wavenumbers using only three input frequencies. The Marmousi inversion result accurately predicts all other data frequencies, demonstrating the redundancy of the data.

Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies

A new approach to seismic migration formalizes the classical diffraction (or common-tangent) stack by relating it to linearized seismic inversion and the generalized Radon transform. This approach recasts migration as the problem of reconstructing the earth’s acoustic scattering potential from its integrals over isochron surfaces. The theory rests on a solution of the wave equation with the geometrical-optics Green function and an approximate inversion formula for the generalized Radon transform. The method can handle both complex velocity models and (nearly) arbitrary configurations of sources and receivers. In this general case, the method can be implemented as a weighted diffraction stack, with the weights determined by tracing rays from image points to the experiment’s sources and receivers. When tested on a finite-difference simulation of a deviated-well vertical seismic profile (a hybrid experiment which is difficult to treat with conventional wave-equation methods), the algorithm accurately reconstructed faulted-earth models. Analytical reconstruction formulas are derived from the general formula for zero-offset and fixed-offset surface experiments in which the background velocity is constant. The zero-offset inversion formula resembles standard Kirchhoff migration. Our analysis provides a direct connection between the experimental setup (source and receiver positions, source wavelet, background velocity) and the spatial resolution of the reconstruction. Synthetic examples illustrate that the lateral resolution in seismic images is described well by the theory and is improved greatly by combining surface data and borehole data. The best resolution is obtained from a zero-offset experiment that surrounds the region to be imaged.

/pdf/a-new-slant-on-seismic-imaging-migration-and-integral-26neyv6pqn.pdf

A new slant on seismic imaging: Migration and integral geometry

A three-dimensional (3D) time-domain numerical scheme for simulation of ground penetrating radar (GPR) on dispersive and inhomogeneous soils with conductive loss is described. The finite-difference time-domain (FDTD) method is used to discretize the partial differential equations for time stepping of the electromagnetic fields. The soil dispersion is modeled by multiterm Lorentz and/or Debye models and incorporated into the FDTD scheme by using the piecewise-linear recursive convolution (PLRC) technique. The dispersive soil parameters are obtained by fitting the model to reported experimental data. The perfectly matched layer (PML) is extended to match dispersive media and used as an absorbing boundary condition to simulate an open space. Examples are given to verify the numerical solution and demonstrate its applications. The 3D PML-PLRC-FDTD formulation facilitates the parallelization of the code. A version of the code is written for a 32-processor system, and an almost linear speedup is observed.

/pdf/finite-difference-time-domain-simulation-of-ground-3mktyd7a25.pdf

Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils

We describe a numerical method for time‐stepping Maxwell’s equations in the two‐dimensional (2-D) TE‐mode, which in a conductive earth reduces to the diffusion equation. The method is based on the classical DuFort‐Frankel finite‐difference scheme, which is both explicit and stable for any size of the time step. With this method, small time steps can be used at early times to track the rapid variations of the field, and large steps can be used at late times, when the field becomes smooth and its rates of diffusion and decay slow down. The boundary condition at the earth‐air interface is handled explicitly by calculating the field in the air from its values at the earth’s surface with an upward continuation based on Laplace’s equation. Boundary conditions in the earth are imposed by using a large, graded grid and setting the values at the sides and bottom to those for a haft‐space. We use the 2-D model to simulate transient electromagnetic (TE) surveys over a thin vertical conductor embedded in a half‐space...

Diffusion of electromagnetic fields into a two-dimensional earth; a finite-difference approach

A new inversion algorithm for the simultaneous reconstruction of permittivity and conductivity recasts the nonlinear inversion as the solution of a coupled set of linear equations. The algorithm is iterative and proceeds through the minimization of two cost functions. At the initial step the data are matched through the reconstruction of the radiating or minimum norm scattering currents; subsequent steps refine the nonradiating scattering currents and the material properties inside the scatterer. Two types of basis functions are constructed for the nonradiating currents: “invisible” (global) basis functions, which are appropriate for discrete measurements and nonradiating (local) basis functions, which are useful in studying the limit of continuous measurements. Reconstructions of square cylinders from multiple source receiver measurements at a single frequency show that the method can handle large contrasts in material properties.

Simultaneous nonlinear reconstruction of two-dimensional permittivity and conductivity

An inverse scattering formula is derived for an ensemble of experiments in which a three-dimensional scattering object is illuminated by an impulsive point source which successively occupies all positions on a surface surrounding the object. For each source position, the scattered field is recorded by receivers at all positions on the surface. The data from such a complete set of experiments contains redundant information about the scattering potential. The formula provides a way of integrating over the redundancy to give a perfect reconstruction of the scattering potential within the Born approximation. The formula involves three steps: first, the data as a function of receiver position are backpropagated (or focused) to an arbitrary image point inside the surface using the Kirchhoff integral; second, the source array is focused on the same point by a second Kirchhoff integral over source positions; finally, the image field created by the double focusing is filtered as a function of frequency and integrated over all frequencies. The formula relies upon an unusual representation of the spatial delta function which follows from an integral relation between the causal Green function, the anti-causal Green function, and their difference. The method can be extended to the case of an arbitrary background medium if ray-asymptotic approximations are used for the Green functions.

https://iopscience.iop.org/article/10.1088/0266-5611/5/6/015/pdf

Michael L. Oristaglio

Papers

A new slant on seismic imaging: Migration and integral geometry

Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils

Diffusion of electromagnetic fields into a two-dimensional earth; a finite-difference approach

Simultaneous nonlinear reconstruction of two-dimensional permittivity and conductivity

An inverse scattering formula that uses all the data