A staggered-grid lowrank finite-difference method for elastic wave extrapolation

Migration is an important step in seismic data processing for oil and gas exploration. The accuracy of migration directly affects the accuracy of subsequent oil and gas reservoir characterization. Reverse-time migration is one of the most accurate migration methods at present. Multi-wave and multicomponent seismic data contain more P- and S-wave information. Making full use of multi-wave and multicomponent seismic data can offer more information about underground structure and lithology, as well as improve the accuracy of seismic exploration. Elastic reverse-time migration (ERTM) has no dip restriction and can be applied to image multi-wave and multicomponent seismic data in complex structural areas and some special lithology structures. However, the surface topography of complex regions has an influence on wavefield and seriously degrades the quality of ERTM’s migration results. We developed a new ERTM method to migrate multi-wave and multicomponent seismic data in the region with complex surface topography. We first fill the layers between the highest and lowest undulating surface with near-surface elastic parameters in a complex topography model to obtain a new model with a horizontal surface. This allows the finite difference (FD) method based on the regular rectangular grid to be used to numerically solve elastic wave equations in the model with complex topography. The decoupled wave equations are used to generate source P- and S-waves and receiver P- and S-waves to reduce crosstalk artefacts in ERTM. A topography-related filter is further used to remove the influence of surface topography on migration results. The scalar imaging condition is also applied to generate PP and PS migration images. Some numerical examples with different complex topographies demonstrate that our proposed ERTM method can remove the influence of complex topography on ERTM’s images and effectively generate high-quality ERTM images.

Vector-based excitation amplitude imaging condition for elastic RTM

The earth [Formula: see text] filtering causes poor illumination of the subsurface. Compensating for the attenuated amplitude and distorted phase is crucial during elastic reverse time migration (ERTM) to improve the imaging quality. Conventional [Formula: see text]-compensated ERTM ([Formula: see text]-ERTM) methods tend to boost the attenuated energy to inverse the [Formula: see text] effects. These methods usually suffer from severe numerical instability because of the unlimited amplification of the high-frequency noise. Low-pass filtering is generally used to stabilize the process, however, at the expense of precision. We have developed a stable compensation approach in this paper. Based on the decoupled fractional Laplacians viscoelastic wave equation, two compensation operators are obtained by extrapolating wavefield in the dispersion-only and viscoelastic media. Because no explicit amplification is included, these two operators are absolutely stable for implementation. To improve the division morbidity for calculating the compensation operators, we use the excitation amplitude criterion and embed the operators into a vector-based [Formula: see text]-compensated excitation amplitude imaging condition. With the derived imaging condition, we could compensate for the absorption accurately without needing to concern the stability issue. The [Formula: see text]-ERTM results for noise-free data are carried out over a simple layered model and a more realistic Marmousi model with an attenuating area to verify the feasibility of the proposed approach. The migration results for noisy data from the Marmousi model further prove that the proposed method performs better fidelity, adaptability, and antinoise performance compared with conventional compensation method with low-pass filtering.

Non-artifact vector P- and S-wave separation for elastic reverse time migration

Summary In this abstract, we discuss some practical issues related to the implementation of reverse time migration (RTM) in a tilted transversely isotropic (TTI) medium. We have proposed a methodology to effectively reduce the cost. Through optimization of the coefficients of high-order difference scheme on the variable grid, larger depth steps at greater depth can be used instead of fixed step size. Synthetic data examples are used to demonstrate the effectiveness of the methodology. Introduction Reverse time migration (RTM) based on solving a two-way wave equation has been one of the most accurate imaging tools available in dealing with large lateral velocity variations and complicated geology such as large dips. RTM for isotropic and vertical transversely isotropic (VTI) medium has been routinely in use in the industry, but the application of RTM in TTI medium is still in the early stage. This lag results from the difficulties in numerical formulations for non-vertical symmetric axes and the subsequent instabilities. The TTI implementation also leads to higher computational costs compared with its isotropic or VTI counterparts. Instead of solving the original anisotropic wave equation, several modified quasi-P wave equations were derived by different authors, that can be used to perform anisotropic RTM for P-waves imaging without cause significant S wave scattering. Alkhalifah (2000) derived a 4th-order partial differential equation (PDE) in time to describe TTI propagation. To avoid solving a high order PDE, Zhou et al. (2006) split Alkhaifah's equation into two coupled 2nd-order PDEs. However it has been observed that their equations cause instability in practice when large dip is at presence. To circumvent this stability issue, Zhang et al. (2009) pursue an equivalent TTI formulation that preserves physical energy with the kinematics retained. Fletcher et al. (2009) proposed adding non-zero shear wave velocity terms to mitigate the instability problem. Unfortunately, their solution adds numerical complexity. In our implementation, to deal with both the stability and accuracy issues in TTI environment, we have to cut the cost dramatically to run it in production environment. To reduce the computational cost of TTI RTM, we use the least-square-root method in the Fourier domain to construct optimized finite difference (FD) schemes of variable grid in depth to approximate the first and second order derivatives. With larger steps in space in deep layers of higher velocity while keeping the same criteria for controlling the numerical dispersion, the total number of grid point is reduced. We analyzed the stability condition for this proposed FD scheme to determine the size of time step, which is related to minimum grid space, maximum velocity as well as the coefficients of FD scheme. With this implementation our RTM has been enhanced greatly in efficiency, and gained success in numerous 3D land imaging projects. Here we will demonstrate the effectiveness of this method using the 2D BP TTI model as a test case.

The purpose of this paper, which builds on previous work (Int. J. Numer. Meth. Engng 2009; 77:1646–1669), is to improve a numerical scheme based on the partition of unity finite element method (PUFEM) for the solution of the time harmonic elastic wave equations. The approach consists to approximate the displacement field by the standard finite element shape functions, enriched locally by superimposing pressure (P) and shear (S) plane waves. The aim is to accurately model two‐dimensional elastic wave problems on relatively coarse mesh grids, capable of containing many wavelengths per nodal spacing, for wide ranges of frequencies. This allows us to relax the traditional requirement of about 10 nodal points per S wavelength. In this work, an exact integration scheme for the linear triangular finite element is developed to evaluate the oscillatory integrals arising from the use of the PUFEM. The main contribution here consists in developing an explicit closed‐form solution for two‐dimensional wave‐based integrals, when the phase variation is linear in the local coordinate element system. The evaluation of the element mass matrix is performed from appropriate edge integrals. All other element matrices, obtained by adequate splitting of the element stress tensor matrix, are simply deduced from the element mass matrix entries.The results show clearly that the proposed integration scheme evaluates accurately the entries of the global matrix with drastic reduction of the computational time. Numerical tests dealing with the scattering of S elastic plane waves by a circular rigid body show that, for the same discretization level, it is possible to improve the accuracy by using large elements associated with high numbers of approximating plane waves rather than using small elements with less plane waves. However, this increases the conditioning and the fill‐in of the global matrix. At high frequency, it is even possible to push the number of degrees of freedom per S wavelength under 2 and still achieve good accuracy. Finally, some remarks on the choice of the numbers of P and S plane waves leading to better accuracy and conditioning are discussed. Copyright © 2010 John Wiley & Sons, Ltd.

Summary We propose a new lowrank finite-difference method for elastic wave propagation in time domain. Firstly, based on elastic wavefield vector decomposition, we derive recursive integral time extrapolation operators for P-and S-wave, respectively. And we employ lowrank decomposition method to approximate these extrapolation operators. Then, we design the corresponding finite-difference scheme. Compared with conventional finite-difference method, the lowrank finite-difference method has higher accuracy and allows us to use large time-step in elastic wave simulation. Numerical experiments confirm that the proposed method can be an effective tool to simulate elastic waves.

Summary Recently a time domain nearly constant Q (NCQ) wave equation derived from Kjartansson's constant Q model has been developed for modeling viscoacoustic wavefield. The wave equation introduces decoupled attenuation and dispersion terms based on two separate fractional Laplacians, which can be easily calculated by spatial Fourier pseudo-spectral method. The fractional orders of the Laplacians are related to Q, and that means the orders are actually spatially varying. However, no desirable approach is presented in the current literatures to handle the varying orders. The fractional Laplacian with a spatially varying order can be exactly represented by a wavenumber-space domain operator. In this abstract we use a lowrank decomposition method to approximate the mixed-domain operator, thus making the NCQ wave equation adapt to large Q contrasts. Additionally, we reformulate the existing velocity-stress-strain NCQ formulation as an equivalent compact velocity-pressure system. The staggered-grid pseudo-spectral (SGPS) method and unsplit convolutional perfectly matched layer (CPML) are adopted in numerical simulations.

Summary Borehole dipole shear reflection imaging has attracted widespread attention in the evaluation of heterogeneous reservoirs in recent years due to its good performance in reflection amplitude, detecting depth and directivity. However, the conventional migration method based on oneway acoustic equation has been proved to be improper for the wave excited by dipole source in this paper. An elastic wave reverse time migration (RTM) is developed for more exactly imaging of dipole reflected wave to find the location of effective reservoir. The problem of storage in RTM is solved by a boundary storage strategy improved for borehole environment in our research, which has been proved to be very effective by the simulated data. The processing of 2D simulated data of formation models with cave, crack and formation interface indicate the good performance of this migration method. The migration results of the field data with this method agree well with the borehole electric imaging logging data. Furthermore, our researches reveal that the effects of the unknown formation parameters, shielding interface of former reflector, would make the interpretation of the migration image difficult.

The conventional reverse time migration utilizes regularly sampled computational grids to simulate wave propagation. Selecting the appropriate grid sampling is important for computational accuracy and efficiency. In general, the uniform-size grid cannot represent the complexity of the geology well. The grid may appear sparse in the low-velocity zone, especially in shallow depths where dispersion may occur. Conversely, it may appear excessively dense in the high-velocity zone, such as at greater depths or within a salt body, which results in higher computational memory and time consumption. To overcome these issues, we developed an efficient and accurate adaptive variable grid discretization method that automatically selects the vertical grid size based on the velocity, depth, and dominant frequency of the wavelet in elastic medium. Then we reformulated the elastic equations based on the adaptive variable grid by introducing a mapping relationship. To test the effectiveness, accuracy, and efficiency of the equation, we implemented it to both the forward propagation and migration of elastic wavefield. Synthetic numerical examples demonstrate that our proposed method can achieve elastic wavefield separation and no significant dispersion phenomenon. The multi-component imaging accuracy of reverse time migration is nearly equivalent to the traditional method, while significantly improving computational efficiency and saving storage space.

P- and S-wave separated elastic reverse time migration for OBC data from fluid-solid coupled media with irregular seabed interfaces

Reverse time migration (RTM) exhibits great advantages over other imaging methods because it is based on computing numerical solutions to a two-way wave equation. It does not suffer from dip limitation like one-way downward continuation techniques do, thus enabling overturned reflections to be imaged. As well as correctly handling multipathing, RTM has the potential to image internal multiples when the boundaries responsible for generating the multiples are present in the model. In isotropic media, one can use a scalar acoustic wave equation for RTM of pressure data. In anisotropic media, P- and SV-waves are coupled together so, formally, elastic wave equations must be used for RTM. A new wave equation for P-waves is proposed in tilted transversely isotropic (TTI) media that can be solved as part of an acoustic anisotropic RTM algorithm, using standard explicit finite differencing. If the shear velocity along the axis of symmetry is set to zero, stable numerical solutions can be computed for media with a vertical axis of symmetry and ε not less than δ. In TTI media with rapid variations in the direction of the axis of symmetry, setting the shear velocity along the axis of symmetry to zero can cause numerical solutions to become unstable. A solution to this problem is proposed that involves using a small amount of nonzero shear velocity. The amount of shear velocity added is chosen to remove triplications from the SV wavefront and to minimize the anisotropic term of the SV reflection coefficient. We show modeling and high-quality RTM results in complex TTI media using this equation.

The parametric coupling between traveling electromagnetic and elastic waves in a ferrimagnet supplied by magnetic pumping is theoretically investigated. The pumping threshold is calculated for steady-state wave propagation in the presence of loss for plane electromagnetic and longitudinal elastic waves in an infinite ferrite medium. This threshold is compared to the threshold for second-order spin-wave instabilities and it is found that these spin wave processes should not prevent the magnetoelastic process from occurring. It is shown that this traveling wave system may be useful for amplifying electromagnetic or elastic waves.

It is important to apply certain intrinsic properties (e.g., the wave propagation directions and particle motion information) in reconstructed wavefields in order to solve the image problems exhibited in elastic reverse time migration (RTM). This paper presents a procedure to decompose the source and receiver wavefields into local plane waves as well as to separate them into pure P and S modes. We generate the partial PP and PS images by crosscorrelation of these plane waves along different directions, and then formulate an imaging condition as a product of an angle-domain operator and the partial images. The new angle-domain imaging condition substantially reduces the artifacts in the PP image and produces the PS image with correct polarizations. However, the imaging procedure involves intensive computations. A numerical algorithm, which can greatly reduce the calculation, is proposed to improve the efficiency.

Elastic reverse time migration with Gaussian beams combines the computational efficiency of ray-based methods with the precision of wave-based approaches, making it a flexible and adaptable imaging tool for multicomponent data. Aiming at the problem of rapid divergence of elastic Gaussian beams width with the increase of propagation distance and the widespread anisotropy in underground media. This paper proposes an anisotropic reverse time migration with Gaussian beams method for multicomponent seismic data. Firstly, an anisotropic ray tracing algorithm based on the phase velocity for elastic wave was derived to correct the impact of anisotropy on imaging. Subsequently, the dynamically focused beams operator was constructed to constrain the energy of the seismic wave within the effective range, which further enhancing the imaging accuracy. The forward and backward continuation wavefields and imaging formulas are constructed by using the Green’s function performed by elastic dynamically focused beams. Finally, the numerical tests on a complex structural model indicate the effectiveness and applicability of the proposed method.