When reliable a priori information is not available, it is difficult to correctly predict near-surface S-wave velocity models from Rayleigh waves through existing techniques, especially in the case of complex geology. To tackle this issue, we have developed a new method: two-grid genetic-algorithm Rayleigh-wave full-waveform inversion (FWI). Adopting a two-grid parameterization of the model, the genetic algorithm inverts for unknown velocities and densities at the nodes of a coarse grid, whereas the forward modeling is performed on a fine grid to avoid numerical dispersion. A bilinear interpolation brings the coarse-grid results into the fine-grid models. The coarse inversion grid allows for a significant reduction in the computing time required by the genetic algorithm to converge. With a coarser grid, there are fewer unknowns and less required computing time, at the expense of the model resolution. To further increase efficiency, our inversion code can perform the optimization using an offset-marching strategy and/or a frequency-marching strategy that can make use of different kinds of objective functions and allows for parallel computing. We illustrate the effect of our inversion method using three synthetic examples with rather complex near-surface models. Although no a priori information was used in all three tests, the long-wavelength structures of the reference models were fairly predicted, and satisfactory matches between “observed” and predicted data were achieved. The fair predictions of the reference models suggest that the final models estimated by our genetic-algorithm FWI, which we call macromodels, would be suitable inputs to gradient-based Rayleigh-wave FWI for further refinement. We also explored other issues related to the practical use of the method in different work and explored applications of the method to field data.

The treatment of multioffset seismic data as an acoustic wave field is becoming increasingly disturbing to many geophysicists who see a multitude of wave phenomena, such as amplitude‐offset variations and shearwave events, which can only be explained by using the more correct elastic wave equation. Not only are such phenomena ignored by acoustic theory, but they are also treated as undesirable noise when they should be used to provide extra information, such as S‐wave velocity, about the subsurface. The problems of using the conventional acoustic wave equation approach can be eliminated via an elastic approach. In this paper, equations have been derived to perform an inversion for P‐wave velocity, S‐wave velocity, and density as well as the P‐wave impedance, S‐wave impedance, and density. These are better resolved than the Lamé parameters. The inversion is based on nonlinear least squares and proceeds by iteratively updating the earth parameters until a good fit is achieved between the observed data and the modeled data corresponding to these earth parameters. The iterations are based on the preconditioned conjugate gradient algorithm. The fundamental requirement of such a least‐squares algorithm is the gradient direction which tells how to update the model parameters. The gradient direction can be derived directly from the wave equation and it may be computed by several wave propagations. Although in principle any scheme could be chosen to perform the wave propagations, the elastic finite‐ difference method is used because it directly simulates the elastic wave equation and can handle complex, and thus realistic, distributions of elastic parameters. This method of inversion is costly since it is similar to an iterative prestack shot‐profile migration. However, it has greater power than any migration since it solves for the P‐wave velocity, S‐wave velocity, and density and can handle very general situations including transmission problems. Three main weaknesses of this technique are that it requires fairly accurate a priori knowledge of the low‐ wavenumber velocity model, it assumes Gaussian model statistics, and it is very computer‐intensive. All these problems seem surmountable. The low‐wavenumber information can be obtained either by a prior tomographic step, by the conventional normal‐moveout method, by a priori knowledge and empirical relationships, or by adding an additional inversion step for low wavenumbers to each iteration. The Gaussian statistics can be altered by preconditioning the gradient direction, perhaps to make the solution blocky in appearance like well logs, or by using large model variances in the inversion to reduce the effect of the Gaussian model constraints. Moreover, with some improvements to the algorithm and more parallel computers, it is hoped the technique will soon become routinely feasible.

全波形反演技术是当前勘探地球物理领域的研究热点之一。新一轮的全波形反演研究触及了波场模拟、梯度估计、数据预处理、目标泛函的选择等诸多深层次的内容,逐渐将全波形反演对实际观测系统、地震数据等要求的局限性以及其具备的潜力揭示出来。通过对全波形反演理论和技术研究进展及应用现状的全面调研,介绍了全波形反演算法研究从时间域到频率域、再到混合域和拉普拉斯域的发展进程；阐明了全波形反演技术已实现海上地震资料应用但对陆上地震资料还没有真正成功实例的应用现状；着重分析了陆上资料全波形反演应用的瓶颈主要在于观测系统的限制、低频数据的缺失、数据预处理面临的挑战以及近地表条件和激发-接收的影响等；指出了发展分步骤、分尺度的反演方法和反演策略以及多种手段的有效联合是实现陆上资料全波形反演的有效途径。

<p>&nbsp;Full waveform inversion(FWI)is one of the hot issue in the geophysics exploration. The new round study involvees in many deeper issues,such as wavefield simulation,gradient estimation,data pre-processing,objective function and so on. Those further studies reveal the limtations of geometry and seismic data, as well as the potiential of FWI,through an overall investigation on the FWI theory,technology and application status,we introduce FWI development,from time domain to frequency domain and Laplace domain,especially to the hybrid domain. It is shown that FWI has been already successfully used in the marine data,however it is still a challenge in land data application;bottlenecks of the land data application are analyzed,such as limitation of geometry,lack of low-frequency component of data,the challenge of data pre-processing,the near surface conditions and source-receiver limitation of land acquisition;Finally it is concluded that step-by-step hierarchical multi-scale strategies and combination of different methods would be available to realize the land data FWI application.</p>

We analyze the kinematic properties of offset‐domain common image gathers (CIGs) and angle‐domain CIGs (ADCIGs) computed by wavefield‐continuation migration. Our results are valid regardless of whether the CIGs were obtained by using the correct migration velocity. They thus can be used as a theoretical basis for developing migration velocity analysis (MVA) methods that exploit the velocity information contained in ADCIGs.

We present a general formula for the back projection of traveltime residuals in traveltime tomography. For special choices of an arbitrary weighting factor this formula reduces to the asymptotic back‐projection term in ray‐tracing tomography (RT), the Woodward‐Rocca method, wavepath eikonal traveltime inversion (WET), and wave‐equation traveltime inversion (WT). This unification provides for an understanding of the differences and similarities among these traveltime tomography methods. The special case of the WET formula leads to a computationally efficient inversion scheme in the space‐time domain that is, in principle, almost as effective as WT inversion yet is an order of magnitude faster. It also leads to an analytic formula for the fast computation of wavepaths. Unlike ray‐tracing tomography, WET partially accounts for band‐limited source and shadow effects in the data. Several numerical tests of the WET method are used to illustrate its properties.

本文首先对20世纪80年代发展起来的全波形反演应用及其在勘探地球物理领域的发展进行了分析；其次，面对定量化、精细化的地震勘探要求，提出了将地震勘探全波形反演与其他数据处理环节或处理技术相结合的研究设想，并展望了全波形反演的发展趋势；最后，论述了全波形反演研究中地震波场数值模拟、反演初始速度模型获取、目标函数形式选择、寻优算法启用及各向异性介质中的应用等关键问题，并总结了通过Laplace域的全波形反演获取反演初始速度模型、结合射线追踪并充分发挥并行计算之于波动方程方法来模拟地震波场的巨大优势，及灵活选用反演目标函数形式和寻优算法更新速度模型参数来加快全波形反演方法的实用化进程。

全波形反演已成功应用于海上地震资料，而陆上地震资料全波形反演应用还面临诸多难点。研究给出了一种针对陆上地震资料的全波形反演策略：首先采用相位拟合互相关目标泛函增强全波形反演的适用性；其次构建动态波数域滤波梯度预处理方法压制梯度噪声；最后采用构造倾角信息约束的自适应高斯平滑算子改善反演结果质量。二维陆上地震资料的反演结果表明，该策略可有效实现陆上地震资料高波数全波形反演速度建模。

Full waveform inversion of marine seismic data has been successfully applied,however,the application of this technology for onshore seismic data faces many difficulties.In this study,a robust strategy of full waveform inversion for onshore seismic data is presented.Firstly,the phase-matched cross-correlation cost function is adopted to improve the applicability of full waveform inversion.Secondly,dynamic wavenumber domain gradient preprocessing method is used to suppress gradient noise.Finally,an adaptive Gaussian smoothing operator is applied to improve the quality of inversion results.The inversion case of two-dimensional onshore seismic data shows that this strategy presented in this study can effectively realize high wavenumber onshore seismic data velocity modeling.

The migration, imaging, or inversion of wide‐aperture cross‐hole data depends on the ability to model wave propagation in complex media for multiple source positions. Computational costs can be considerably reduced in frequency‐domain imaging by modeling the frequency‐domain steady‐state equations, rather than the time‐domain equations of motion. I develop a frequency‐domain approach in this note that is competitive with time‐domain modeling when solutions for multiple sources are required or when only a limited number of frequency components of the solution are required.

The nonlinear inverse problem for seismic reflection data is solved in the acoustic approximation. The method is based on the generalized least‐squares criterion, and it can handle errors in the data set and a priori information on the model. Multiply reflected energy is naturally taken into account, as well as refracted energy or surface waves. The inverse problem can be solved using an iterative algorithm which gives, at each iteration, updated values of bulk modulus, density, and time source function. Each step of the iterative algorithm essentially consists of a forward propagation of the actual sources in the current model and a forward propagation (backward in time) of the data residuals. The correlation at each point of the space of the two fields thus obtained yields the corrections of the bulk modulus and density models. This shows, in particular, that the general solution of the inverse problem can be attained by methods strongly related to the methods of migration of unstacked data, and commercially competitive with them.

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 thewave-physics complexity. Different hierarchical multiscale strategies are designed to mitigate the nonlinearity and ill-posedness of FWI by incorporating progressively shorter wavelengths in the parameter space. Synthetic and real-data case studies address reconstructing various parameters, from [Formula: see text] and [Formula: see text] 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 techniques to make 3D elastic FWI feasible.

The problem of interpretation of seismic reflection data can be posed with sufficient generality using the concepts of inverse theory. In its roughest formulation, the inverse problem consists of obtaining the Earth model for which the predicted data best fit the observed data. If an adequate forward model is used, this best model will give the best images of the Earth’s interior. Three parameters are needed for describing a perfectly elastic, isotropic, Earth: the density ρ(x) and the Lamé parameters λ(x) and μ(x), or the density ρ(x) and the P-wave and S-wave velocities α(x) and β(x). The choice of parameters is not neutral, in the sense that although theoretically equivalent, if they are not adequately chosen the numerical algorithms in the inversion can be inefficient. In the long (spatial) wavelengths of the model, adequate parameters are the P-wave and S-wave velocities, while in the short (spatial) wavelengths, P-wave impedance, S-wave impedance, and density are adequate. The problem of inversion of waveforms is highly nonlinear for the long wavelengths of the velocities, while it is reasonably linear for the short wavelengths of the impedances and density. Furthermore, this parameterization defines a highly hierarchical problem: the long wavelengths of the P-wave velocity and short wavelengths of the P-wave impedance are much more important parameters than their counterparts for S-waves (in terms of interpreting observed amplitudes), and the latter are much more important than the density. This suggests solving the general inverse problem (which must involve all the parameters) by first optimizing for the P-wave velocity and impedance, then optimizing for the S-wave velocity and impedance, and finally optimizing for density. The first part of the problem of obtaining the long wavelengths of the P-wave velocity and the short wavelengths of the P-wave impedance is similar to the problem solved by present industrial practice (for accurate data interpretation through velocity analysis and “prestack migration”). In fact, the method proposed here produces (as a byproduct) a generalization to the elastic case of the equations of “prestack acoustic migration.” Once an adequate model of the long wavelengths of the P-wave velocity and of the short wavelengths of the P-wave impedance has been obtained, the data residuals should essentially contain information on S-waves (essentially P-S and S-P converted waves). Once the corresponding model of S-wave velocity (long wavelengths) and S-wave impedance (short wavelengths) has been obtained, and if the remaining residuals still contain information, an optimization for density should be performed (the short wavelengths of impedances do not give independent information on density and velocity independently). Because the problem is nonlinear, the whole process should be iterated to convergence; however, the information from each parameter should be independent enough for an interesting first solution.

全波形反演利用叠前地震波场的运动学和动力学信息重建地下地球物理参数，具有揭示复杂地质背景下构造细节及岩性的潜在能力。介绍了频率域反演联合时间域正演的优化算法的三维声波全波形反演方法及算法实现。该算法的正演实现使用了伪保守形式的一阶速度-压力场声波方程。应用三维SEG/EAGE模型反演实例验证了三维声波全波形反演方法的有效性。

The full waveform inversion (FWI) method makes use of the kinematic and dynamic information of pre-stack seismic data to rebuild underground geophysical parameters.Due to its high resolution,it has the potential of revealing structure details and lithology characteristic under complex geological background.We present some methodological and algorithmic developments of 3D acoustic full waveform inversion by frequency-domain inversion combined with time-domain forward.The FWI algorithm relies on a pseudo-conservative form of the velocity-stress acoustic wave equation.The numerical simulation data of 3D SEG/EAGE target model verified 3D acoustic full waveform inversion algorithm.

A crosshole experiment was carried out in a layered sedimentary environment in which a normal fault is known to cut through the section. Initial traveltime inversions produced stable but low‐resolution images from which the fault could be only vaguely inferred. To image the fault, wavefield inversion was used to produce a velocity model consistent with the detailed phase and amplitude of the data at a number of frequencies. Our wavefield inversion scheme uses a classical, descent‐type algorithm for decreasing the data misfit by iteratively computing the gradient of this misfit by repeated forward and backward propagations. Our propagator is a full‐wave equation, frequency‐domain, acoustic, finite‐difference method. The use of the frequency‐space domain yields computational advantages for multisource data and allows an easy incorporation of viscous effects. By running wavefield inversion on the field data, a quantitative velocity image was produced that yielded a significantly improved image of the fault (when compared with the original traveltime inversions). Because the original field data were noisy and contained a high degree of multiple scattering (from the layering of the sediments), the transmitted arrivals were selectively windowed to enhance the image. The sediments at the site were strongly attenuating; we therefore used a viscoacoustic model during the modeling and the inversion that correctly simulated the observed decrease in amplitude with increasing frequency and source‐receiver offset. Furthermore, since the traveltime inversion indicated a high degree of anisotropy at the site, a fixed, homogeneous level of anisotropy was used during the inversion. Tests at varying levels of anisotropy confirmed the improvement in image quality and in data fit when anisotropy was incorporated. The final image was verified by examining the distribution of the residuals in the frequency domain, by comparing time‐domain modeled wavefields with the observed data, and by direct comparison with borehole logs.

Estimation of reflector depth and seismic velocity from seismic reflection data can be formulated as a general inverse problem. The method used to solve this problem is similar to tomographic techniques in medical diagnosis and we refer to it as seismic reflection tomography. Seismic tomography is formulated as an iterative Gauss‐Newton algorithm that produces a velocity‐depth model which minimizes the difference between traveltimes generated by tracing rays through the model and traveltimes measured from the data. The input to the process consists of traveltimes measured from selected events on unstacked seismic data and a first‐guess velocity‐depth model. Usually this first‐guess model has velocities which are laterally constant and is usually based on nearby well information and/or an analysis of the stacked section. The final model generated by the tomographic method yields traveltimes from ray tracing which differ from the measured values in recorded data by approximately 5 ms root‐mean‐square. The indeterminancy of the inversion and the associated nonuniqueness of the output model are both analyzed theoretically and tested numerically. It is found that certain aspects of the velocity field are poorly determined or undetermined. This technique is applied to an example using real data where the presence of permafrost causes a near‐surface lateral change in velocity. The permafrost is successfully imaged in the model output from tomography. In addition, depth estimates at the intersection of two lines differ by a significantly smaller amount than the corresponding estimates derived from conventional processing.

Over the past 10 years, ray-based postmigration grid tomography has become the standard model-building tool for seismic depth imaging. While the basics of the method have remained unchanged since the late 1990s, the problems it solves have changed dramatically. This evolution has been driven by exploration demands and enabled by computer power. There are three main areas of change. First, standard model resolution has increased from a few thousand meters to a few hundred meters. This order of magnitude improvement may be attributed to both high-quality, complex residual-moveout data picked as densely as [Formula: see text] to [Formula: see text] vertically and horizontally, and to a strategy of working down from long-wavelength to short-wavelength solutions. Second, more and more seismic data sets are being acquired along multiple azimuths, for improved illumination and multiple suppression. High-resolution velocity tomography must solve for all azimuths simultaneously, to prevent short-wavelength velocity heterogeneity from being mistaken for azimuthal anisotropy. Third, there has been a shift from predominantly isotropic to predominantly anisotropic models, both VTI and TTI. With four-component data, anisotropic grid tomography can be used to build models that tie PZ and PS images in depth.

This efficient multiscale method for time-domain waveform tomography incorporates filters that are more efficient than Hamming-window filters. A strategy for choosing optimal frequency bands is proposed to achieve computational efficiency in the time domain. A staggered-grid, explicit finite-difference method with fourth-order accuracy in space and second-order accuracy in time is used for forward modeling and the adjoint calculation. The adjoint method is utilized in inverting for an efficient computation of the gradient directions. In the multiscale approach, multifrequency data and multiple grid sizes are used to overcome somewhat the severe local minima problem of waveform tomography. The method is applied successfully to 1D and 2D heterogeneous models; it can accurately recover low- and high-wavenumber components of the velocity models. The inversion result for the 2D model demonstrates that the multiscale method is computationally efficient and converges faster than a conventional, single-scale method.

The nonlinear problem of inversion of seismic waveforms can be set up using least‐squares methods. The inverse problem is then reduced to the problem of minimizing a lp;nonquadratic) function in a space of many [Formula: see text] variables. Using gradient methods leads to iterative algorithms, each iteration implying a forward propagation generated by the actual sources, a backward propagation generated by the data residuals (acting as if they were sources), and a correlation at each point of the space of the two fields thus obtained, which gives the updated model. The quality of the results of any inverse method depends heavily on the realism of the forward modeling. Finite‐difference schemes are a good choice relative to realism because, although they are time‐consuming, they give excellent results. Numerical tests performed with multioffset synthetic data from a two‐dimensional model prove the feasibility of the approach. If only surface‐recorded reflections are used, the high spatial frequency content of the model (but not the low spatial frequencies) is recovered in few (≃5) iterations. By using transmitted data also (e.g., between two boreholes), all the spatial frequencies are recovered. Since the problem is nonlinear, if the initial guess is far enough from the true solution, the iterative algorithm may converge into a secondary solution. A nonlinear inversion with 8 shots, each shot recorded at 400 receiver locations, with 700 samples in each seismogram, corresponding to a 2-D model described by 40 000 grid points, takes approximately 1 hour in a CRAY 1S supercomputer.

Green's functions for an elastic layered medium can be expressed as a double integral over frequency and horizontal wavenumber. We show that, for any time window, the wavenumber integral can be exactly represented by a discrete summation. This discretization is achieved by adding to the particular point source an infinite set of specified circular sources centered around the point source and distributed at equal radial interval. Choice of this interval is dependent on the length of time desired for the point source response and determines the discretized set of horizontal wavenumbers which contribute to the solution. Comparisons of the results obtained with those derived using the two-dimensional discretization method (Bouchon, 1979) are presented. They show the great accuracy of the two methods.

When finite‐difference methods are used to solve the elastic wave equation in a discontinuous medium, the error has two dominant components. Dispersive errors lead to artificial wave trains. Errors from interfaces lead to circular wavefronts emanating from each location where the interface appears “jagged” to the rectangular grid. The pseudospectral method can be viewed as the limit of finite differences with infinite order of accuracy. With this method, dispersive errors are essentially eliminated. The mappings introduced in this paper also eliminate the other dominant error source. Test calculations confirm that these mappings significantly enhance the already highly competitive pseudospectral method with only a very small additional cost. Although the mapping method is described here in connection with the pseudospectral method, it can also be used with high‐order finite‐difference approximations.

Recent interest in the extraction of fine detail from field seismograms has stimulated the search for numerical modeling procedures which can produce synthetic seismograms for complex subsurface geometries and for arbitrary source‐receiver separations. Among the various techniques available for this purpose, the replacement of the two‐dimensional wave equation by a suitable finite‐difference representation offers distinct advantages. This approach is simple and may be readily implemented. It automatically accounts for the proper relative amplitudes of the various arrivals and includes the contributions of converted waves, Rayleigh waves, diffractions from faulted zones, and head waves. Two computational schemes have been examined. For the so‐called “homogeneous formulation,” the standard boundary conditions between media of different elastic properties must be satisfied explicitly. In the case of the alternate “heterogeneous formulation”, these elastic properties may be specified at each grid point of a finite‐difference mesh, and the boundary conditions are satisfied implicitly. The proper simulation of the source requires special treatment for both cases. Synthetic seismograms computed for several models of exploration interest serve to illustrate how the technique may help the interpreter. The examples also illustrate various implementational aspects of the finite‐difference approach, which involve such phenomena as grid dispersion, artificial reflections from the edge of the model, and choice of spatial and temporal sampling increments.

The treatment of multioffset seismic data as an acoustic wave field is becoming increasingly disturbing to many geophysicists who see a multitude of wave phenomena, such as amplitude‐offset variations and shearwave events, which can only be explained by using the more correct elastic wave equation. Not only are such phenomena ignored by acoustic theory, but they are also treated as undesirable noise when they should be used to provide extra information, such as S‐wave velocity, about the subsurface. The problems of using the conventional acoustic wave equation approach can be eliminated via an elastic approach. In this paper, equations have been derived to perform an inversion for P‐wave velocity, S‐wave velocity, and density as well as the P‐wave impedance, S‐wave impedance, and density. These are better resolved than the Lamé parameters. The inversion is based on nonlinear least squares and proceeds by iteratively updating the earth parameters until a good fit is achieved between the observed data and the modeled data corresponding to these earth parameters. The iterations are based on the preconditioned conjugate gradient algorithm. The fundamental requirement of such a least‐squares algorithm is the gradient direction which tells how to update the model parameters. The gradient direction can be derived directly from the wave equation and it may be computed by several wave propagations. Although in principle any scheme could be chosen to perform the wave propagations, the elastic finite‐ difference method is used because it directly simulates the elastic wave equation and can handle complex, and thus realistic, distributions of elastic parameters. This method of inversion is costly since it is similar to an iterative prestack shot‐profile migration. However, it has greater power than any migration since it solves for the P‐wave velocity, S‐wave velocity, and density and can handle very general situations including transmission problems. Three main weaknesses of this technique are that it requires fairly accurate a priori knowledge of the low‐ wavenumber velocity model, it assumes Gaussian model statistics, and it is very computer‐intensive. All these problems seem surmountable. The low‐wavenumber information can be obtained either by a prior tomographic step, by the conventional normal‐moveout method, by a priori knowledge and empirical relationships, or by adding an additional inversion step for low wavenumbers to each iteration. The Gaussian statistics can be altered by preconditioning the gradient direction, perhaps to make the solution blocky in appearance like well logs, or by using large model variances in the inversion to reduce the effect of the Gaussian model constraints. Moreover, with some improvements to the algorithm and more parallel computers, it is hoped the technique will soon become routinely feasible.

Iterative inversion methods have been unsuccessful at inverting seismic data obtained from complicated earth models (e.g. the Marmousi model), the primary difficulty being the presence of numerous local minima in the objective function. The presence of local minima at all scales in the seismic inversion problem prevent iterative methods of inversion from attaining a reasonable degree of convergence to the neighborhood of the global minimum. The multigrid method is a technique that improves the performance of iterative inversion by decomposing the problem by scale. At long scales there are fewer local minima and those that remain are further apart from each other. Thus, at long scales iterative methods can get closer to the neighborhood of the global minimum. We apply the multigrid method to a subsampled, low‐frequency version of the Marmousi data set. Although issues of source estimation, source bandwidth, and noise are not treated, results show that iterative inversion methods perform much better when employed with a decomposition by scale. Furthermore, the method greatly reduces the computational burden of the inversion that will be of importance for 3-D extensions to the method.

I describe the properties of a fourth‐order accurate space, second‐order accurate time, two‐dimensional P-SV finite‐difference scheme based on the Madariaga‐Virieux staggered‐grid formulation. The numerical scheme is developed from the first‐order system of hyperbolic elastic equations of motion and constitutive laws expressed in particle velocities and stresses. The Madariaga‐Virieux staggered‐grid scheme has the desirable quality that it can correctly model any variation in material properties, including both large and small Poisson’s ratio materials, with minimal numerical dispersion and numerical anisotropy. Dispersion analysis indicates that the shortest wavelengths in the model need to be sampled at 5 gridpoints/wavelength. The scheme can be used to accurately simulate wave propagation in mixed acoustic‐elastic media, making it ideal for modeling marine problems. Explicitly calculating both velocities and stresses makes it relatively simple to initiate a source at the free‐surface or within a layer and to satisfy free‐surface boundary conditions. Benchmark comparisons of finite‐difference and analytical solutions to Lamb’s problem are almost identical, as are comparisons of finite‐difference and reflectivity solutions for elastic‐elastic and acoustic‐elastic layered models.

The nonlinear inverse problem for seismic reflection data is solved in the acoustic approximation. The method is based on the generalized least‐squares criterion, and it can handle errors in the data set and a priori information on the model. Multiply reflected energy is naturally taken into account, as well as refracted energy or surface waves. The inverse problem can be solved using an iterative algorithm which gives, at each iteration, updated values of bulk modulus, density, and time source function. Each step of the iterative algorithm essentially consists of a forward propagation of the actual sources in the current model and a forward propagation (backward in time) of the data residuals. The correlation at each point of the space of the two fields thus obtained yields the corrections of the bulk modulus and density models. This shows, in particular, that the general solution of the inverse problem can be attained by methods strongly related to the methods of migration of unstacked data, and commercially competitive with them.

We recognized that the envelope fluctuation and decay of seismic records carries ultra low-frequency (ULF, i.e., the frequency below the lowest frequency in the source spectrum) signals that can be used to estimate the long-wavelength velocity structure. We then developed envelope inversion for the recovery of low-wavenumber components of media (smooth background), so that the initial model dependence of waveform inversion can be reduced. We derived the misfit function and the corresponding gradient operator for envelope inversion. To understand the long-wavelength recovery by the envelope inversion, we developed a nonlinear seismic signal model, the modulation signal model, as the basis for retrieving the ULF data and studied the nonlinear scale separation by the envelope operator. To separate the envelope data from the wavefield data (envelope extraction), a demodulation operator (envelope operator) was applied to the waveform data. Numerical tests using synthetic data for the Marmousi model proved the validity and feasibility of the proposed approach. The final results of combined [Formula: see text] (envelope-inversion for smooth background plus waveform-inversion for high-resolution velocity structure) indicated that it can deliver much improved results compared with regular full-waveform inversion (FWI) alone. Furthermore, to test the independence of the envelope to the source frequency band, we used a low-cut source wavelet (cut from 5 Hz below) to generate the synthetic data. The envelope inversion and the combined [Formula: see text] showed no appreciable difference from the full-band source results. The proposed envelope inversion is also an efficient method with very little extra work compared with conventional FWI.

Conventional full-waveform inversion (FWI) using the least-squares norm as a misfit function is known to suffer from cycle-skipping issues that increase the risk of computing a local rather than the global minimum of the misfit. The quadratic Wasserstein metric has proven to have many ideal properties with regard to convexity and insensitivity to noise. When the observed and predicted seismic data are considered to be two density functions, the quadratic Wasserstein metric corresponds to the optimal cost of rearranging one density into the other, in which the transportation cost is quadratic in distance. Unlike the least-squares norm, the quadratic Wasserstein metric measures not only amplitude differences but also global phase shifts, which helps to avoid cycle-skipping issues. We have developed a new way of using the quadratic Wasserstein metric trace by trace in FWI and compare it with the global quadratic Wasserstein metric via the solution of the Monge-Ampère equation. We incorporate the quadratic Wasserstein metric technique into the framework of the adjoint-state method and apply it to several 2D examples. With the corresponding adjoint source, the velocity model can be updated using a quasi-Newton method. Numerical results indicate the effectiveness of the quadratic Wasserstein metric in alleviating cycle-skipping issues and sensitivity to noise. The mathematical theory and numerical examples demonstrate that the quadratic Wasserstein metric is a good candidate for a misfit function in seismic 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.

Seismic waveform inversion is one of many geophysical problems which can be identified as a nonlinear multiparameter optimization problem. Methods based on local linearization fail if the starting model is too far from the true model. We have investigated the applicability of “Genetic Algorithms” (GA) to the inversion of plane‐wave seismograms. Like simulated annealing, genetic algorithms use a random walk in model space and a transition probability rule to help guide their search. However, unlike a single simulated annealing run, the genetic algorithms search from a randomly chosen population of models (strings) and work with a binary coding of the model parameter set. Unlike a pure random search, such as in a “Monte Carlo” method, the search used in genetic algorithms is not directionless. Genetic algorithms essentially consist of three operations, selection, crossover, and mutation, which involve random number generation, string copies, and some partial string exchanges. The choice of the initial population, the probabilities of crossover and mutation are crucial for the practical implementation of the algorithm. We investigated the effects of these parameters in the inversion of plane‐wave seismograms in which a normalized crosscorrelation function was used as the objective or fitness function (E). We also introduce the concept of “update” probability to control the influence of past generations. The combination of a low value of mutation probability (∼0.01), a moderate value of the crossover probability (∼0.6) and a high value of update probability (∼0.9) are found to be optimal for the convergence of the algorithm. Further, we show that concepts from simulated annealing can be used effectively for the stretching of the fitness function which helps in the convergence of the algorithm. Thus, we propose to use exp (E/T) rather than E as the fitness function, where T (analogous to temperature in simulated annealing) is a properly chosen parameter which can change slowly with each generation. Also, by repeating the GA optimization procedure several times with different randomly chosen initial model populations, we derive “a very good subset” of models from the entire model space and calculate the a posteriori probability density σ(m) ∝ exp (E(m)/T). The σ(m) ’s are then used to calculate a “mean” model, which is found to be close to the true model.

There is a deep and useful connection between statistical mechanics (the behavior of systems with many degrees of freedom in thermal equilibrium at a finite temperature) and multivariate or combinatorial optimization (finding the minimum of a given function depending on many parameters). A detailed analogy with annealing in solids provides a framework for optimization of the properties of very large and complex systems. This connection to statistical mechanics exposes new information and provides an unfamiliar perspective on traditional optimization problems and methods.

基于常规模拟退火算法的零偏VSP全波形反演面临着计算量大和耗时长的问题。为此提出了一种不同阶段对应不同扰动模型和退火方式的分段快速模拟退火(segmented fast simulated annealing,SFSA)反演策略,以提高零偏VSP资料全波形反演的效率。在反演前期采用大模型扰动空间和较慢温度衰减速度,充分发挥全局搜索能力,而在后期引入限制因子产生扰动模型,在迭代不断增加的时候逐渐减小模型的扰动空间,同时采用较快的温度衰减速度,有效提高反演的速度,使反演快速收敛到最优解。基于相同的初始温度和马尔可夫链长度,分别利用基于SFSA和非常快速模拟退火(very fast simulated annealing,VFSA)方法进行零偏VSP纵波速度全波形反演测试。结果表明,基于SFSA的反演方法的反演效率提高约50%,在迭代次数更少的条件下能获得更好的反演效果。基于SFSA的零偏VSP全波形反演具有高效和高精度的特点,其反演结果为地震地质层位标定、成果解释及油气预测奠定了基础。

<p>&nbsp;The application of the full-waveform inversion technique based on the conventional simulated annealing algorithm to zero-offset VSP data is time consuming and complex.This paper therefore proposes an efficient segmented fast simulated annealing (SFSA) strategy,which uses different disturbance models and annealing methods at different stages to improve the waveform inversion efficiency of zero-offset VSP data.In the early stage of inversion,large model perturbation space and slow decay rate of temperature were adopted to fully perform global search;in late stage,the perturbation model with restriction factor was adopted;the model perturbation space was gradually reduced as the iteration increased,with a faster decay rate of temperature,to improve the efficiency of the inversion and help the algorithm converge to the optimal solution quickly.Finally,SFSA and conventional very fast simulated annealing (VFSA) algorithms were compared for P-wave velocity inversion of zero-offset VSP.The test results showed that the SFSA algorithm can improve the inversion efficiency by 50%,while a better inversion effect could be obtained using less iteration.The proposed method could lay the foundation for the calibration of seismic geological horizons,seismic interpretation,and reservoir prediction.</p>