A Nearly Analytical Discrete Method for Wave-Field Simulations in 2D Porous Media

The nearly analytic discrete method (NADM) is a perturbation method originally proposed by Yang et al. (2003) [26] for acoustic and elastic waves in elastic media. This method is based on a truncated Taylor series expansion and interpolation approximations and it can suppress effectively numerical dispersions caused by the discretizating the wave equations when too-coarse grids are used. In the present work, we apply the NADM to simulating acoustic and elastic wave propagations in 2D porous media. Our method enables wave propagation to be simulated in 2D porous isotropic and anisotropic media. Numerical experiments show that the error of the NADM for the porous case is less than those of the conventional finite-difference method (FDM) and the so-called Lax-Wendroff correction (LWC) schemes. The three-component seismic wave fields in the 2D porous isotropic medium are simulated and compared with those obtained by using the LWC method and exact solutions. Several characteristics of wave propagating in porous anisotropic media, computed by the NADM, are also reported in this study. Promising numerical results illustrate that the NADM provides a useful tool for large-scale porous problems and it can suppress effectively numerical dispersions.     

Multi-Frequency Contrast Source Inversion for Reflection Seismic Data

In the contrast source inversion (CSI) method, the contrast sources (equivalent scattering sources) and the contrast (parameter perturbation) are iteratively reconstructed by an alternating optimization scheme. Traditionally integral equation CSI method is formulated for transmission tomography using analytic Green's function in homogeneous background. To extend the method to the case of reflection seismology, in this paper, we use WKBJ method to compute the Green's function of depth dependent background media and the solving method of equations to initialize the contrast source of different frequencies, resulting in an efficient method to invert multi-frequency reflection seismic data – multi-frequency contrast source inversion method (MFCSI). Numerical results for the Marmousi model show that MFCSI method can obtain good results even when low frequency data are missing, in which case the conventional FWI fails.     

Subspace Methods in Multi-Parameter Seismic Full Waveform Inversion

In full waveform inversion (FWI) high-resolution subsurface model parameters are sought. FWI is normally treated as a nonlinear least-squares inverse problem, in which the minimum of the corresponding misfit function is found by updating the model parameters. When multiple elastic or acoustic properties are solved for, simple gradient methods tend to confuse parameter classes. This is referred to as parameter cross-talk; it leads to incorrect model solutions, poor convergence and strong dependence on the scaling of the different parameter types. Determining step lengths in a subspace domain, rather than directly in terms of gradients of different parameters, is a potentially valuable approach to address this problem. The particular subspace used can be defined over a span of different sets of data or different parameter classes, provided it involves a small number of vectors compared to those contained in the whole model space. In a subspace method, the basis vectors are defined first, and a local minimum is found in the space spanned by these. We examine the application of the subspace method within acoustic FWI in determining simultaneously updates for velocity and density. We first discuss the choice of basis vectors to construct the spanned space, from linear updates by distinguishing only the contributions of different parameter classes towards nonlinear updates by adding the contributions of higher-order perturbations of each parameter class. The numerical character of FWI solutions generated via subspace methods involving different basis vectors is then analyzed and compared with traditional FWI methods. The subspace methods can provide better reconstructions of the model, especially for the velocity, as well as improved convergence rates, while the computational costs are still comparable with the traditional FWI methods.     

A Weighted Runge-Kutta Discontinuous Galerkin Method for 3D Acoustic and Elastic Wave-Field Modeling

Numerically solving 3D seismic wave equations is a key requirement for forward modeling and inversion. Here, we propose a weighted Runge-Kutta discontinuous Galerkin (WRKDG) method for 3D acoustic and elastic wave-field modeling. For this method, the second-order seismic wave equations in 3D heterogeneous anisotropic media are transformed into a first-order hyperbolic system, and then we use a discontinuous Galerkin (DG) solver based on numerical-flux formulations for spatial discretization. The time discretization is based on an implicit diagonal Runge-Kutta (RK) method and an explicit iterative technique, which avoids solving a large-scale system of linear equations. In the iterative process, we introduce a weighting factor. We investigate the numerical stability criteria of the 3D method in detail for linear and quadratic spatial basis functions. We also present a 3D analysis of numerical dispersion for the full discrete approximation of acoustic equation, which demonstrates that the WRKDG method can efficiently suppress numerical dispersion on coarse grids. Numerical results for several different 3D models including homogeneous and heterogeneous media with isotropic and anisotropic cases show that the 3D WRKDG method can effectively suppress numerical dispersion and provide accurate wave-field information on coarse mesh.     

Numerical Modeling of Anisotropic Elastic-Wave Sensitivity Propagation for Optimal Design of Time-Lapse Seismic Surveys

Reliable subsurface time-lapse seismic monitoring is crucial for many geophysical applications, such as enhanced geothermal system characterization, geologic carbon utilization and storage, and conventional and unconventional oil/gas reservoir characterization, etc. We develop an elastic-wave sensitivity propagation method for optimal design of cost-effective time-lapse seismic surveys considering the fact that most of subsurface geologic layers and fractured reservoirs are anisotropic instead of isotropic. For anisotropic media, we define monitoring criteria using qP- and qS-wave sensitivity energies after decomposing qP- and qS-wave components from the total elastic-wave sensitivity wavefield using a hybrid time- and frequency-domain approach. Geophones should therefore be placed at locations with significant qP- and qS-wave sensitivity energies for cost-effective time-lapse seismic monitoring in an anisotropic geology setting. Our numerical modeling results for a modified anisotropic Hess model demonstrate that, compared with the isotropic case, subsurface anisotropy changes the spatial distributions of elastic-wave sensitivity energies. Consequently, it is necessary to consider subsurface anisotropies when designing the spatial distribution of geophones for cost-effective time-lapse seismic monitoring. This finding suggests that it is essential to use our new anisotropic elastic-wave sensitivity modeling method for optimal design of time-lapse seismic surveys to reliably monitor the changes in subsurface reservoirs, fracture zones or target monitoring regions.     

Concepts in the Direct Waveform Inversion (DWI) Using Explicit Time-Space Causality

The Direct Waveform Inversion (DWI) is a recently proposed waveform inversion idea that has the potential to simultaneously address several existing challenges in many full waveform inversion (FWI) schemes. A key ingredient in DWI is the explicit use of the time-space causality property of the wavefield in the inversion which allows us to convert the global nonlinear optimization problem in FWI, without information loss, into local linear inversions that can be readily solved. DWI is a recursive scheme which sequentially inverts for the subsurface model in a shallow-to-deep fashion. Therefore, there is no need for a global initial velocity model to implement DWI. DWI is unconditionally convergent when the reflection traveltime from the boundary of inverted model is beyond the finite recording time in seismic data. In order for DWI to work, DWI must use the full seismic wavefield including interbed and free surface multiples and it combines seismic migration and velocity model inversion into one process. We illustrate the concepts in DWI using 1D and 2D models.     

Target-Oriented Inversion of Time-Lapse Seismic Waveform Data

Full waveform inversion of time-lapse seismic data can be used as a means of estimating the reservoir changes due to the production. Since the repeated computations for the monitor surveys lead to a large computational cost, time-lapse full waveform inversion is still considered to be a challenging task. To address this problem, we present an efficient target-oriented inversion scheme for time-lapse seismic data using an integral equation formulation with Gaussian beam based Green's function approach. The proposed time-lapse approach allows one to perform a local inversion within a small region of interest (e.g. a reservoir under production) for the monitor survey. We have verified that the T-matrix approach is indeed naturally target-oriented, which was mentioned by Jakobsen and Ursin [24] and allows one to reduce the computational cost of time-lapse inversion by focusing the inversion on the target-area only. This method is based on a new version of the distorted Born iterative T-matrix inverse scattering method. The Gaussian beam and T-matrix are used in this approach to perform the wavefield computation for the time-lapse inversion in the baseline model from the survey surface to the target region. We have provided target-oriented inversion results of the synthetic time-lapse waveform data, which shows that the proposed scheme reduces the computational cost significantly.     

Stability Conditions for Wave Simulation in 3-D Anisotropic Media with the Pseudospectral Method

Simulation of elastic wave propagation has important applications in many area such as inverse problem and geophysical exploration. In this paper, stability conditions for wave simulation in 3-D anisotropic media with the pseudospectral method are investigated. They can be expressed explicitly by elasticity constants which are easy to be applied in computations. The3-D wave simulation for two typical anisotropic media, transversely isotropic media and orthorhombic media, are carried out. The results demonstrate some satisfactory behaviors of the pseudospectral method.     

Towards a Theoretical Background for Strong-Scattering Inversion – Direct Envelope Inversion and Gel'fand-Levitan-Marchenko Theory

Strong-scattering inversion or the inverse problem for strong scattering has different physical-mathematical foundations from the weak-scattering case. Seismic inversion based on wave equation for strong scattering cannot be directly solved by Newton's local optimization method which is based on weak-nonlinear assumption. Here I try to illustrate the connection between the Schrödinger inverse scattering (inverse problem for Schrödinger equation) by GLM (Gel'fand-Levitan-Marchenko) theory and the direct envelope inversion (DEI) using reflection data. The difference between wave equation and Schrödinger equation is that the latter has a potential independent of frequency while the former has a frequency-square dependency in the potential. I also point out that the traditional GLM equation for potential inversion can only recover the high-wavenumber components of impedance profile. I propose to use the Schrödinger impedance equation for direct impedance inversion and introduce a singular impedance function which also corresponds to a singular potential for the reconstruction of impedance profile, including discontinuities and long-wavelength velocity structure. I will review the GLM theory and its application to impedance inversion including some numerical examples. Then I analyze the recently developed multiscale direct envelope inversion (MS-DEI) and its connection to the inverse Schrödinger scattering. It is conceivable that the combination of strong-scattering inversion (inverse Schrödinger scattering) and weak-scattering inversion (local optimization based inversion) may create some inversion methods working for a whole range of inversion problems in geophysical exploration.     

A Spectral Iterative Method for the Computation of Effective Properties of Elastically Inhomogeneous Polycrystals

We report an efficient phase field formalism to compute the stress distribution in polycrystalline materials with arbitrary elastic inhomogeneity and anisotropy. The dependence of elastic stiffness tensor on grain orientation is taken into account, and the elastic equilibrium equation is solved using a spectral iterative perturbation method. We discuss its applications to computing residual stress distribution in systems containing arbitrarily shaped cavities and cracks (with zero elastic modulus) and to determining the effective elastic properties of polycrystals and multilayered composites.     

Forward Scattering and Volterra Renormalization for Acoustic Wavefield Propagation in Vertically Varying Media

We extend the full wavefield modeling with forward scattering theory and Volterra Renormalization to a vertically varying two-parameter (velocity and density) acoustic medium. The forward scattering series, derived by applying Born-Neumann iterative procedure to the Lippmann-Schwinger equation (LSE), is a well known tool for modeling and imaging. However, it has limited convergence properties depending on the strength of contrast between the actual and reference medium or the angle of incidence of a plane wave component. Here, we introduce the Volterra renormalization technique to the LSE. The renormalized LSE and related Neumann series are absolutely convergent for any strength of perturbation and any incidence angle. The renormalized LSE can further be separated into two sub-Volterra type integral equations, which are then solved non-iteratively. We apply the approach to velocity-only, density-only, and both velocity and density perturbations. We demonstrate that this Volterra Renormalization modeling is a promising and efficient method. In addition, it can also provide insight for developing a scattering theory-based direct inversion method.     

Numerical Regularized Moment Method for High Mach Number Flow

This paper is a continuation of our earlier work [SIAM J. Sci. Comput., 32(2010), pp. 2875–2907] in which a numerical moment method with arbitrary order of moments was presented. However, the computation may break down during the calculation of the structure of a shock wave with Mach number M₀≥3. In this paper, we concentrate on the regularization of the moment systems. First, we apply the Maxwell iteration to the infinite moment system and determine the magnitude of each moment with respect to the Knudsen number. After that, we obtain the approximation of high order moments and close the moment systems by dropping some high-order terms. Linearization is then performed to obtain a very simple regularization term, thus it is very convenient for numerical implementation. To validate the new regularization, the shock structures of low order systems are computed with different shock Mach numbers.     

Direct Simulation Methods for Scalar-Wave Envelopes in Two-Dimensional Layered Random Media Based on the Small-Angle Scattering Approximation

This study presents stochastic methods to simulate wave envelopes in layered random media. High-frequency Seismograms of small earthquakes are so complex due to lithospheric inhomogeneity that seismologists often analyze wave envelopes rather than wave traces to quantify the subsurface inhomogeneity. Since the statistical properties of the inhomogeneity vary regionally, it is important to develop and examine direct envelope simulation methods for non-uniform random media. As a simple example, this study supposes plane wave propagation through two-layer random media in 2-D composed of weak and strong inhomogeneity zones. The characteristic spatial-scale of the inhomogeneity is supposed to be larger than the wavelength, where small-angle scattering around the forward direction dominates large-angle scattering. Two envelope simulation methods based on the small-angle scattering approximation are examined. One method is to solve a differential equation for the two-frequency mutual coherence function with the Markov approximation. The other is to solve the stochastic ray bending process by using the Monte Carlo method based on the Markov approximation for the mutual coherence function. The resultant wave envelopes of the two methods showed excellent coincidence both for uniform and for two-layer random media. Furthermore, we confirmed the validity of the two methods comparing with the envelopes made from the finite difference simulations of waves. The two direct envelope simulation methods presented in this study can be a mathematical base for the study of high-frequency wave propagation through randomly inhomogeneous lithosphere in seismology.     

A Level Set Method for the Inverse Problem of Wave Equation in the Fluid-Saturated Porous Media

In this paper, a level set method is applied to the inverse problem of 2-D wave equation in the fluid-saturated media. We only consider the situation that the parameter to be recovered takes two different values, which leads to a shape reconstruction problem. A level set function is used to present the discontinuous parameter, and a regularization functional is applied to the level set function for the ill-posed problem. Then the resulting inverse problem with respect to the level set function is solved by using the damped Gauss-Newton method. Numerical experiments show that the method can recover parameter with complicated geometry and the noise in the observation data.     

An Edge-Based Smoothed Finite Element Method with TBC for the Elastic Wave Scattering by an Obstacle

Elastic wave scattering has received ever-increasing attention in military and medical fields due to its high-precision solution. In this paper, an edge-based smoothed finite element method (ES-FEM) combined with the transparent boundary condition (TBC) is proposed to solve the elastic wave scattering problem by a rigid obstacle with smooth surface, which is embedded in an isotropic and homogeneous elastic medium in two dimensions. The elastic wave scattering problem satisfies Helmholtz equations with coupled boundary conditions obtained by Helmholtz decomposition. Firstly, the TBC of the elastic wave scattering is constructed by using the analytical solution to Helmholtz equations, which can truncate the boundary value problem (BVP) in an unbounded domain into the BVP in a bounded domain. Then the formulations of ES-FEM with the TBC are derived for Helmholtz equations with coupled boundary conditions. Finally, several numerical examples illustrate that the proposed ES-FEM with the TBC (ES-FEM-TBC) can work effectively and obtain more stable and accurate solution than the standard FEM with the TBC (FEM-TBC) for the elastic wave scattering problem.     

Linear Wavefield Optimization Using a Modified Source

Recorded seismic data are sensitive to the Earth's elastic properties, and thus, they carry information of such properties in their waveforms. The sensitivity of such waveforms to the properties is nonlinear causing all kinds of difficulties to the inversion of such properties. Inverting directly for the components forming the wave equation, which includes the wave equation operator (or its perturbation), and the wavefield, as independent parameters enhances the convexity of the inverse problem. The optimization in this case is provided by an objective function that maximizes the data fitting and the wave equation fidelity, simultaneously. To enhance the practicality and efficiency of the optimization, I recast the velocity perturbations as secondary sources in a modified source function, and invert for the wavefield and the modified source function, as independent parameters. The optimization in this case corresponds to a linear problem. The inverted functions can be used directly to extract the velocity perturbation. Unlike gradient methods, this optimization problem is free of the Born approximation limitations in the update, including single scattering and cross talk that may arise for example in the case of multi sources. These specific features are shown for a simple synthetic example, as well as the Marmousi model.     

Performance Analysis of a High-Order Discontinuous Galerkin Method Application to the Reverse Time Migration

This work pertains to numerical aspects of a finite element method based discontinuous functions. Our study focuses on the Interior Penalty Discontinuous Galerkin method (IPDGM) because of its high-level of flexibility for solving the full wave equation in heterogeneous media. We assess the performance of IPDGM through a comparison study with a spectral element method (SEM). We show that IPDGM is as accurate as SEM. In addition, we illustrate the efficiency of IPDGM when employed in a seismic imaging process by considering two-dimensional problems involving the Reverse Time Migration.     

Amplitude Compensation for One-Way Wave Propagators in Inhomogeneous Media and Its Application to Seismic Imaging

The WKBJ solution for the one-way wave equations in media with smoothly varying velocity variation with depth, c(z), is reformulated from the principle of energy flux conservation for acoustic media. The formulation is then extended to general heterogeneous media with local angle domain methods by introducing the concepts of Transparent Boundary Condition (TBC) and Transparent Propagator (TP). The influence of the WKBJ correction on image amplitudes in seismic imaging, such as depth migration in exploration seismology, is investigated in both smoothly varying c(z) and general heterogeneous media. We also compare the effect of the propagator amplitude compensation with the effect of the acquisition aperture correction on the image amplitude. Numerical results in a smoothly varying c(z) medium demonstrate that the WKBJ correction significantly improves the one-way wave propagator amplitudes, which, after compensation, agree very well with those from the full wave equation method. Images for a point scatterer in a smoothly varying c(z) medium show that the WKBJ correction has some improvement on the image amplitude, though it is not very significant. The results in a general heterogeneous medium (2D SEG/EAGE salt model) show similar phenomena. When the acquisition aperture correction is applied, the image improves significantly in both the smoothly varying c(z) medium and the 2D SEG/EAGE salt model. The comparisons indicate that although the WKBJ compensation for propagator amplitude may be important for forward modeling (especially for wide-angle waves), its effect on the image amplitude in seismic imaging is much less noticeable compared with the acquisition aperture correction for migration with limited acquisition aperture in general heterogeneous media.     

A Stable Q Compensated Reverse Time Migration Method Based on Excitation Amplitude Imaging Condition

The stability and efficiency, especially the stability, are generally concerned issues in Q compensated reverse time migration (Q-RTM). The instability occurs because of the exponentially boosted high frequency ambient noise during the forward or backward seismic wavefield propagation. The regularization and low-pass filtering methods are two effective strategies to control the instability of the wave propagation in Q-RTM. However, the regularization parameters are determined experimentally, and the wavefield cannot be recovered accurately. The low-pass filtering method cannot balance the selection of cutoff frequency for varying Q values, and may damage the effective signals, especially when the signal-to-noise ratio (SNR) of the seismic data is low, the Q-RTM will be a highly unstable process. In order to achieve the purpose of stability, the selection of cutoff frequency will be small enough, which can cause great damage to the effective high frequency signals. In this paper, we present a stable Q-RTM algorithm based on the excitation amplitude imaging condition, which can compensate both the amplitude attenuation and phase dispersion. Unlike the existing Q-RTM algorithms enlarging the amplitude, the exponentially attenuated seismic wavefield will be used during both the forward and backward wavefield propagation of Q-RTM. Therefore, the new Q-RTM algorithm is relative stable, even for the low SNR seismic data. In order to show the accuracy and stability of our stable Q-RTM algorithm clearly, an example based on Graben model will be illustrated. Then, a realistic BP gas chimney model further demonstrates that the proposed method enjoys good stability and anti-noise performance compared with the traditional Q-RTM with amplitude amplification. Compare the Q-RTM images of these two models to the reference images obtained by the acoustic RTM with acoustic seismic data, the new Q-RTM results match the reference images quite well. The proposed method is also tested using a field seismic data, the result shows the effectiveness of our proposed method.