Theory of equivalent staggered‐grid schemes: application to rotated and standard grids in anisotropic media

The previous finite‐difference numerical schemes designed for direct application to second‐order elastic wave equations in terms of displacement components are strongly dependent on Poisson's ratio. This fact makes theses schemes useless for modelling in offshore regions or even in onshore regions where there is a high Poisson's ratio material. As is well known, the use of staggered‐grid formulations solves this drawback. The most common staggered‐grid algorithms apply central‐difference operators to the first‐order velocity–stress wave equations. They have been one of the most successfully applied numerical algorithms for seismic modelling, although these schemes require more computational memory than those mentioned based on second‐order wave equations. The goal of the present paper is to develop a general theory that enables one to formulate equivalent staggered‐grid schemes for direct application to hyperbolic second‐order wave equations. All the theory necessary to formulate these schemes is presented in detail, including issues regarding source application, providing a general method to construct staggered‐grid formulations to a wide range of cases. Afterwards, the equivalent staggered‐grid theory is applied to anisotropic elastic wave equations in terms of only velocity components (or similar displacements) for two important cases: general anisotropic media and vertical transverse isotropy media using, respectively, the rotated and the standard staggered‐grid configurations. For sake of simplicity, we present the schemes in terms of velocities in the second‐ and fourth‐order spatial approximations, with second‐order approximation in time for 2D media. However, the theory developed is general and can be applied to any set of second‐order equations (in terms of only displacement, velocity, or even stress components), using any staggered‐grid configuration with any spatial approximation order in 2D or 3D cases. Some of these equivalent staggered‐grid schemes require less computer memory than the corresponding standard staggered‐grid formulation, although the programming is more evolved. As will be shown in theory and practice, with numerical examples, the equivalent staggered‐grid schemes produce results equivalent to corresponding standard staggered‐grid schemes with computational advantages. Finally, it is important to emphasize that the equivalent staggered‐grid theory is general and can be applied to other modelling contexts, e.g., in electrodynamical and poroelastic wave propagation problems in a systematic and simple way.     

A hybrid absorbing boundary condition for elastic staggered‐grid modelling

We recently proposed an efficient hybrid scheme to absorb boundary reflections for acoustic wave modelling that could attain nearly perfect absorptions. This scheme uses weighted averaging of wavefields in a transition area, between the inner area and the model boundaries. In this paper we report on the extension of this scheme to 2D elastic wave modelling with displacement‐stress formulations on staggered grids using explicit finite‐difference, pseudo‐implicit finite‐difference and pseudo‐spectral methods. Numerical modelling results of elastic wave equations with hybrid absorbing boundary conditions show great improvement for modelling stability and significant absorption for boundary reflections, compared with the conventional Higdon absorbing boundary conditions, demonstrating the effectiveness of this scheme for elastic wave modelling. The modelling results also show that the hybrid scheme works well in 2D rotated staggered‐grid modelling for isotropic medium, 2D staggered‐grid modelling for vertically transversely isotropic medium and 2D rotated staggered‐grid modelling for tilted transversely isotropic medium.     

Optimal implicit staggered‐grid finite‐difference schemes based on the sampling approximation method for seismic modelling

We propose new implicit staggered‐grid finite‐difference schemes with optimal coefficients based on the sampling approximation method to improve the numerical solution accuracy for seismic modelling. We first derive the optimized implicit staggered‐grid finite‐difference coefficients of arbitrary even‐order accuracy for the first‐order spatial derivatives using the plane‐wave theory and the direct sampling approximation method. Then, the implicit staggered‐grid finite‐difference coefficients based on sampling approximation, which can widen the range of wavenumber with great accuracy, are used to solve the first‐order spatial derivatives. By comparing the numerical dispersion of the implicit staggered‐grid finite‐difference schemes based on sampling approximation, Taylor series expansion, and least squares, we find that the optimal implicit staggered‐grid finite‐difference scheme based on sampling approximation achieves greater precision than that based on Taylor series expansion over a wider range of wavenumbers, although it has similar accuracy to that based on least squares. Finally, we apply the implicit staggered‐grid finite difference based on sampling approximation to numerical modelling. The modelling results demonstrate that the new optimal method can efficiently suppress numerical dispersion and lead to greater accuracy compared with the implicit staggered‐grid finite difference based on Taylor series expansion. In addition, the results also indicate the computational cost of the implicit staggered‐grid finite difference based on sampling approximation is almost the same as the implicit staggered‐grid finite difference based on Taylor series expansion.     

Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity‡

This paper presents a Lebedev finite difference scheme on staggered grids for the numerical simulation of wave propagation in an arbitrary 3D anisotropic elastic media. The main concept of the scheme is the definition of all the components of each tensor (vector) appearing in the elastic wave equation at the corresponding grid points, i.e., all of the stresses are stored in one set of nodes while all of the velocity components are stored in another. Meanwhile, the derivatives with respect to the spatial directions are approximated to the second order on two‐point stencils. The second‐order scheme is presented for the sake of simplicity and it is easy to expand to a higher order. Another approach, widely‐known as the rotated staggered grid scheme, is based on the same concept; therefore, this paper contains a detailed comparative analysis of the two schemes. It is shown that the dispersion condition of the Lebedev scheme is less restrictive than that of the rotated staggered grid scheme, while the stability criteria lead to approximately equal time stepping for the two approaches. The main advantage of the proposed scheme is its reduced computational memory requirements. Due to a less restrictive dispersion condition and the way the media parameters are stored, the Lebedev scheme requires only one‐third to two‐thirds of the computer memory required by the rotated staggered grid scheme. At the same time, the number of floating point operations performed by the Lebedev scheme is higher than that for the rotated staggered grid scheme.     

Pseudo‐spectral method using rotated staggered grid for elastic wave propagation in 3D arbitrary anisotropic media

Staggering grid is a very effective way to reduce the Nyquist errors and to suppress the non‐causal ringing artefacts in the pseudo‐spectral solution of first‐order elastic wave equations. However, the straightforward use of a staggered‐grid pseudo‐spectral method is problematic for simulating wave propagation when the anisotropy level is greater than orthorhombic or when the anisotropic symmetries are not aligned with the computational grids. Inspired by the idea of rotated staggered‐grid finite‐difference method, we propose a modified pseudo‐spectral method for wave propagation in arbitrary anisotropic media. Compared with an existing remedy of staggered‐grid pseudo‐spectral method based on stiffness matrix decomposition and a possible alternative using the Lebedev grids, the rotated staggered‐grid‐based pseudo‐spectral method possesses the best balance between the mitigation of artefacts and efficiency. A 2D example on a transversely isotropic model with tilted symmetry axis verifies its effectiveness to suppress the ringing artefacts. Two 3D examples of increasing anisotropy levels demonstrate that the rotated staggered‐grid‐based pseudo‐spectral method can successfully simulate complex wavefields in such anisotropic formations.     

二维各向异性介质中地震波场的高阶同位网格有限差分模拟

本文将DRP/opt MacCormack有限差分格式用于模拟二维各向异性介质中的地震波传播.DRP/opt MacCormack是一种同位网格下的差分格式,避免了传统的交错网格在计算各向异性问题时由于变量插值而导致的误差.而且相对于低阶同位网格差分格式,它具有低色散、低耗散的优点.此格式将中心差分算子分成前向和后向两个空间单边差分,然后在4-6步Runge-Kutta时间积分中使用单边差分组合.在具有垂直对称轴的横向各向同性（VTI）模型下,通过对比DRP/opt MacCormack有限差分和谱元方法的模拟结果,验证了前者具有很高的精度和稳定性.由于实际地质条件下TI介质的对称轴通常是倾斜的(TTI),本文在二维三分量框架下模拟TTI介质中的地震波场.结果显示横波分裂和切平面/反平面运动耦合的特征.数值实验表明DRP/opt MacCormack是一种有效的研究各向异性介质中地震波传播规律的差分格式.     

一种基于全局优化的交错网格有限差分法

在地震波场数值模拟中, 交错网格有限差分技术得到了广泛的应用, 但是在弹性模量变化较大时, 通常会因插值而导致模拟误差增大. 旋转交错网格可以很好地克服这个缺点, 因而适合于各向异性介质正演模拟. 但是对于同样大小的网格单元, 旋转交错网格需要的步长比常规交错网格要大, 这会使梯度和散度算子的误差增大因而更易产生空间数值频散. 针对这些问题, 本文提出了旋转交错网格与紧致有限差分相结合的方法, 并基于模拟退火算法进行全局优化, 压制数值频散, 拓宽波数范围. 数值模拟结果表明, 此方法可以有效地压制数值频散, 且具有较高的模拟精度.      

Finite-difference elastic wave propagation in 2D heterogeneous transversely isotropic media1

The velocity-stress formulation for propagation of elastic seismic waves through 2D heterogeneous transversely isotropic media of arbitrary orientation is presented. The equations are recast into a finite-difference scheme and solved numerically using fourth-order spatial operators and a second-order temporal operator on a staggered grid. Absorbing, free-surface and symmetry boundary conditions have been implemented. Test cases compare well with other published solutions. Synthetic seismograms are calculated over two idealized models: (i) vertical fractures in granite with a dolerite sill reflector and (ii) a dipping anisotropic shale. Comparisons with the isotropic counterparts show significant differences which may have to be accounted for in seismic processing in the future.     

裂隙型单斜介质中多方位地面三分量记录模拟

针对裂隙型储集层中更具代表性的各向异性介质模型,即在各向同性背景介质中含有两组斜交的垂直裂隙所构成的单斜各向异性介质模型,利用时间和空间上可达任意阶的高阶交错网格有限差分技术,对具有不同裂隙填充物性质的单斜介质中波的传播快照进行了模拟.结果证实各向异性介质中波的传播速度随传播方向的不同而产生明显的差异;裂隙填充物的性质对于速度各向异性具有很大的影响.另外,利用坐标旋转法,对水平层状各向异性介质中多方位地面三分量记录进行了模拟,结果表明了方位各向异性介质中,波的传播速度不仅随入射角的变化而变化,同时也随观测方位的不同而产生差异.数值模拟结果为进一步利用地面多方位地震属性进行各向异性参数的反演及裂隙参数的描述提供理论基础.     

ELASTIC WAVE PROPAGATION IN TRANSVERSELY ISOTROPIC MEDIA USING FINITE DIFFERENCES1

When treating the forward full waveform case, a fast and accurate algorithm for modelling seismic wave propagation in anisotropic inhomogeneous media is of considerable value in current exploration seismology. Synthetic seismograms were computed for P-SV wave propagation in transversely isotropic media. Among the various techniques available for seismic modelling, the finite-difference method possesses both the power and flexibility to model wave propagation accurately in anisotropic inhomogeneous media bounded by irregular interfaces. We have developed a fast high-order vectorized finite-difference algorithm adapted for the vector supercomputer. The algorithm is based on the fourth-order accurate MacCormack-type splitting scheme. Solving the equivalent first-order hyperbolic system of equations, instead of the second-order wave equation, avoids computation of the spatial derivatives of the medium's anisotropic elastic parameters. Examples indicate that anisotropy plays an important role in modelling the kinematic and the dynamic properties of the wave propagation and should be taken into account when necessary.     

Modified image algorithm to simulate seismic channel waves in 3D tunnel model with rugged free surfaces

The existence of rugged free‐surface three‐dimensional tunnel conditions in the coal seams, caused either by geological or mining processes, will inevitably influence wave propagation characteristics when the seismic waves go through the coal mines. Thus, a modified image algorithm has been developed to account for seismic channel waves propagating through this complicated topography with irregular free surfaces. Moreover, the seismic channel waves commonly exhibit damped and dispersive signatures, which is not only because of their own unique sandwich geometry of rock–coal–rock but also because of the viscoelastic behavior of coal. Considering the complexity of programming in three‐dimensional tunnel models with rugged free surfaces, an optimized vacuum grid search algorithm, enabling to model highly irregular topography and to compute efficiently, is also proposed when using high‐order staggered finite‐difference scheme to simulate seismic channel wave propagations in viscoelastic media. The numerical simulations are implemented to investigate the accuracy and stability of the method and the impact of coal's viscoelastic behavior on seismic channel wave propagation characteristics. The results indicate that the automatic vacuum grid search algorithm can be easily merged into high‐order staggered finite‐difference scheme, which can efficiently be applied to calculate three‐dimensional tunnel models with rugged free surfaces in the viscoelastic media. The simulation also suggests that the occurrence of a three‐dimensional tunnel with free surfaces has a remarkable influence on the seismic channel wave propagation characteristics and elastic energy distribution.     

Effective velocities in fractured media: a numerical study using the rotated staggered finite‐difference grid

The modelling of elastic waves in fractured media with an explicit finite‐difference scheme causes instability problems on a staggered grid when the medium possesses high‐contrast discontinuities (strong heterogeneities). For the present study we apply the rotated staggered grid. Using this modified grid it is possible to simulate the propagation of elastic waves in a 2D or 3D medium containing cracks, pores or free surfaces without hard‐coded boundary conditions. Therefore it allows an efficient and precise numerical study of effective velocities in fractured structures. We model the propagation of plane waves through a set of different, randomly cracked media. In these numerical experiments we vary the wavelength of the plane waves, the crack porosity and the crack density. The synthetic results are compared with several static theories that predict the effective P‐ and S‐wave velocities in fractured materials in the long wavelength limit. For randomly distributed and randomly orientated, rectilinear, non‐intersecting, thin, dry cracks, the numerical simulations of velocities of P‐, SV‐ and SH‐waves are in excellent agreement with the results of the modified (or differential) self‐consistent theory. On the other hand for intersecting cracks, the critical crack‐density (porosity) concept must be taken into account. To describe the wave velocities in media with intersecting cracks, we propose introducing the critical crack‐density concept into the modified self‐consistent theory. Numerical simulations show that this new formulation predicts effective elastic properties accurately for such a case.     

选择性滤波同位网格有限差分法在地震波数值模拟中的应用.

引入计算空气声学领域的选择性滤波同位网格有限差分算法（SFFD法）用于二维地震波数值模拟．SFFD法使用经过优化的11点DRP同位网格差分格式,对空间一阶导数进行离散近似,同时采用选择性滤波方法来消除同位网格差分所产生的格点高频振荡,它既提高了数值模拟的精度, 又保证了求解过程的稳定性．数值实验结果表明,SFFD法能够达到O(Delta;x8, Delta;t4)阶交错网格算法同样的精度,同时该方法还具有很强的适应性,能够应用于存在着强泊松比差异的介质模型中,完整地模拟地震波传播过程中各类型的波场,并且对复杂非均匀介质的适应能力也很好．此外,由于避免了交错网格算法在曲线坐标系和一般各向异性介质的数值模拟时所需进行的复杂的插值运算, SFFD法在这些问题上也有着很好的应用前景．      

裂缝介质旋转交错网格正演模拟

油气勘探实践表明, 裂缝是油气的储存空间或运移通道, 裂缝介质地震波场的研究越来越受到关注. 实际地层中裂缝形成受多种因素控制, 物理属性比较复杂, 表现出强烈的各向异性. 由于上覆地层的压力使得水平或低陡倾角裂缝存在较少, 大多数为高陡倾角裂缝, 利用线性滑动裂缝介质的等效理论将高陡倾角裂缝介质等效为横向各向同性介质, 便于实际应用. 本文采用各向异性弹性波旋转交错网格模拟方法对含裂缝介质单炮记录进行模拟与分析. 结果表明: 裂缝的存在相当于人为增加反射界面; 裂缝密度越大, 裂缝纵横比越小, 裂缝充填物与背景介质弹性性质差别越大, 引起的反射波能量变化越大. 本文模拟结果为利用地震数据进行裂缝介质参数反演与储层识别及油气预测提供了依据.      

Seismic wavefield modeling in media with fluid-filled fractures and surface topography

We present a finite difference (FD) method for the simulation of seismic wave fields in fractured medium with an irregular (non-flat) free surface which is beneficial for interpreting exploration data acquired in mountainous regions. Fractures are introduced through the Coates-Schoenberg approach into the FD scheme which leads to local anisotropic properties of the media where fractures are embedded. To implement surface topography, we take advantage of the boundary-conforming grid and map a rectangular grid onto a curved one. We use a stable and explicit second-order accurate finite difference scheme to discretize the elastic wave equations (in a curvilinear coordinate system) in a 2D heterogeneous transversely isotropic medium with a horizontal axis of symmetry (HTI). Efficiency tests performed by different numerical experiments clearly illustrate the influence of an irregular free surface on seismic wave propagation in fractured media which may be significant to mountain seismic exploration. The tests also illustrate that the scattered waves induced by the tips of the fracture are re-scattered by the features of the free surface topography. The scattered waves provoked by the topography are re-scattered by the fractures, especially Rayleigh wave scattering whose amplitudes are much larger than others and making it very difficult to identify effective information from the fractures.     

Multi‐block finite‐difference method for 3D elastodynamic simulations in anisotropic subhorizontally layered media

Prediction of elastic full wavefields is required for reverse time migration, full waveform inversion, borehole seismology, seismic modelling, etc. We propose a novel algorithm to solve the Navier wave equation, which is based on multi‐block methodology for high‐order finite‐difference schemes on curvilinear grids. In the current implementation, the blocks are subhorizontal layers. Smooth anisotropic heterogeneous media in each layer can have strong discontinuities at the interfaces. A curvilinear adaptive hexahedral grid in blocks is generated by mapping the original 3D physical domain onto a parametric cube with horizontal layers and interfaces. These interfaces correspond to the main curvilinear physical contrast interfaces of a subhorizontally layered formation. The top boundary of the parametric cube handles the land surface with smooth topography. Free‐surface and solid–solid transmission boundary conditions at interfaces are approximated with the second‐order accuracy. Smooth media in the layers are approximated up to sixth‐order spatial schemes. All expected properties of the developed algorithm are demonstrated in numerical tests using corresponding parallel message passing interface code.     

波场模拟中的数值频散分析与校正策略

波动方程有限差分法正演模拟,对认识地震波传播规律、进行地震属性研究、地震资料地质解释、储层评价等,均具有重要的理论和实际意义.但有限差分法本身固有存在着数值频散问题,数值频散在正演模拟中是一种严重的干扰,会降低波场模拟的精度与分辨率.针对TI介质波场模拟的交错网格有限差分方法,本文从空间网格离散、时间网格离散和算子近似等三个方面对其产生的数值频散进行了分析,并结合其他学者的研究成果给出了TI介质波场模拟中压制数值频散的方法与策略：在已知介质频散关系时,对差分算子可实施算子校正;通过提高差分方程的阶数来提高波场模拟精度;采用流体力学中守恒式方程的通量校正传输方法来压制波场模拟中的数值频散;在实际正演模拟时,采用交错网格高阶有限差分方程,不仅在空间上采用高阶差分,而且在时间上也要采用高阶差分,否则只在单一方向上（空间或时间）提高方程的阶数对压制数值频散也不会取得理想的效果.     

VTI介质起伏地表地震波场模拟

起伏地表下地震波场模拟有助于解释主动源和被动源地震探测中穿过山脉和盆地的测线所获得的资料.然而传统的有限差分法处理起伏的自由边界比较困难,为了克服这一困难,我们将笛卡尔坐标系的各向异性介质弹性波方程和自由边界条件变换到曲线坐标系中,采用一种稳定的、显式的二阶精度有限差分方法离散(曲线坐标系)VTI介质中的弹性波方程;对...     

Characterization of dipping fractures in a transverselyisotropic background

Although it is believed that natural fracture sets predominantly have near‐vertical orientation, oblique stresses and some other mechanisms may tilt fractures away from the vertical. Here, we examine an effective medium produced by a single system of obliquely dipping rotationally invariant fractures embedded in a transversely isotropic with a vertical symmetry axis (VTI) background rock. This model is monoclinic with a vertical symmetry plane that coincides with the dip plane of the fractures. Multicomponent seismic data acquired over such a medium possess several distinct features that make it possible to estimate the fracture orientation. For example, the vertically propagating fast shear wave (and the fast converted PS‐wave) is typically polarized in the direction of the fracture strike. The normal‐moveout (NMO) ellipses of horizontal reflection events are co‐orientated with the dip and strike directions of the fractures, which provides an independent estimate of the fracture azimuth. However, the polarization vector of the slow shear wave at vertical incidence does not lie in the horizontal plane – an unusual phenomenon that can be used to evaluate fracture dip. Also, for oblique fractures the shear‐wave splitting coefficient at vertical incidence becomes dependent on fracture infill (saturation). A complete medium‐characterization procedure includes estimating the fracture compliances and orientation (dip and azimuth), as well as the Thomsen parameters of the VTI background. We demonstrate that both the fracture and background parameters can be obtained from multicomponent wide‐azimuth data using the vertical velocities and NMO ellipses of PP‐waves and two split SS‐waves (or the traveltimes of PS‐waves) reflected from horizontal interfaces. Numerical tests corroborate the accuracy and stability of the inversion algorithm based on the exact expressions for the vertical and NMO velocities.     

Seismic reflection coefficients of faults at low frequencies: a model study

We use linear slip theory to evaluate seismic reflections at non‐welded interfaces, such as faults or fractures, sandwiched between general anisotropic media and show that at low frequencies the real parts of the reflection coefficients can be approximated by the responses of equivalent welded interfaces, whereas the imaginary parts can be related directly to the interface compliances. The imaginary parts of low frequency seismic reflection coefficients at fault zones can be used to estimate the interface compliances, which can be related to fault properties upon using a fault model. At normal incidence the expressions uncouple and the complex‐valued P‐wave reflection coefficient can be related linearly to the normal compliance. As the normal compliance is highly sensitive to the infill of the interface, it can be used for gas/fluid identification in the fault plane. Alternatively, the tangential compliance of a fault can be estimated from the complex‐valued S‐wave reflection coefficient. The tangential compliance can provide information on the crack density in a fault zone. Coupling compliances can be identified and quantified by the observation of PS conversion at normal incidence, with a comparable linear relationship.