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.     

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

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

Application of the nearly perfectly matched layer for seismic wave propagation in 2D homogeneous isotropic media

Numerical modelling plays an important role in helping us understand the characteristics of seismic wave propagation. The presence of spurious reflections from the boundaries of the truncated computational domain is a prominent problem in finite difference computations. The nearly perfectly matched layer has been proven to be a very effective boundary condition to absorb outgoing waves in both electromagnetic and acoustic media. In this paper, the nearly perfectly matched layer technique is applied to elastic isotropic media to further test the method's absorbing ability. The staggered‐grid finite‐difference method (fourth‐order accuracy in space and second‐order accuracy in time) is used in the numerical simulation of seismic wave propagation in 2D Cartesian coordinates. In the numerical tests, numerical comparisons between the nearly perfectly matched layer and the convolutional perfectly matched layer, which is considered the best absorbing layer boundary condition, is also provided. Three numerical experiments demonstrate that the nearly perfectly matched layer has a similar performance to the convolutional perfectly matched layer and can be a valuable alternative to other absorbing layer boundary conditions.     

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.     

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.     

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.     

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

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

用于弹性波方程数值模拟的有限差分系数确定方法

有限差分方法是波场数值模拟的一个重要方法,交错网格差分格式比规则网格差分格式稳定性更好,但方法本身都存在因网格化而形成的数值频散效应,这会降低波场模拟的精度与分辨率.为了缓解有限差分算子的数值频散效应,精确求解空间偏导数,本文把求解波动方程的线性化方法推广到用于求解弹性波方程交错网格有限差分系数;同时应用最大最小准则作为模拟退火(SA)优化算法求解差分系数的数值频散误差判定标准来求解有限差分系数.通过上述两种方法,分别利用均匀各向同性介质和复杂构造模型进行了数值正演模拟和数值频散分析,并与传统泰勒展开算法、最小二乘算法进行比较,验证了线性化方法和模拟退火方法都能有效压制数值频散,并比较了各个算法的特点.     

Elastic wave modelling method based on the displacement–velocity fields: an improving nearly analytic discrete approximation

In this paper, we present an improvement to our previously published nearly analytic discrete method (NADM) to solve acoustic and elastic wave equations. We compare the numerical errors of the improved NADM with the original NADM and also the fourth-order Lax–Wendroff correction and present examples of three-component wave fields in 2D transversely isotropic media with strong velocity contrasts. Comparing with the original NADM, we find that the improved method requires significantly less storage space and can increase the time accuracy from second order of the original NADM to fourth order, while the space accuracy remains the same as that of the original one. Theoretical analyses and numerical results suggest that our improved NADM is suitable for large-scale numerical modelling as it can effectively suppress numerical dispersion and source-generated noises caused by discretising wave equations when too-coarse grids are used.     

TTI介质的交错网格伪P波正演方法

研究了三维弱各向异性近似下,利用伪P波(伪纵波)模拟弹性波场P分量在倾斜对称轴的横向各向同性(TTI)介质中的传播过程,并对比了分别基于弹性Hooke定律、弹性波投影和运动学色散方程所建立的三种二阶差分伪P波方程的正演特点.目前这些伪P波方程数值计算主要采用规则网格差分,但是规则网格在TTI模拟中有低效率、低精度以及不稳定的缺点.为了提高计算的精度,本文构建出相应方程的交错网格有限差分格式.通过对比伪P波方程在三维TTI介质中不同的数值模拟的表达形式,本文认为基于色散方程所建立的伪P波方程在模拟弹性波中P波传播的过程中具有最小的噪声.本文分析不同的各向同性对称轴空间角度的频散特征,并引入适当的横波速度维持计算的稳定.二维模型算例表明,本文提出的交错网格正演算法可以得到稳定光滑的伪P波正演波场.使用本文交错网格算法对二维BP TTI模型的逆时偏移也具有较稳定的偏移结果.     

Finite-difference modelling of S-wave splitting in anisotropic media

We have implemented a 3D finite‐difference scheme to simulate wave propagation in arbitrary anisotropic media. The anisotropic media up to orthorhombic symmetry were modelled using a standard staggered grid scheme and beyond (monoclinic and triclinic) using a rotated staggered grid scheme. The rationale of not using rotated staggered grid for all types of anisotropic media is that the rotated staggered grid schemes are more expensive than standard staggered grid schemes. For a 1D azimuthally anistropic medium, we show a comparison between the seismic data generated by our finite‐difference code and by the reflectivity algorithm; they are in excellent agreement. We conducted a study on zero‐offset shear‐wave splitting using the finite‐difference modelling algorithm using the rotated staggered grid scheme. Our S‐wave splitting study is mainly focused on fractured media. On the scale of seismic wavelenghts, small aligned fractures behave as an equivalent anisotropic medium. We computed the equivalent elastic properties of the fractures and the background in which the fractures were embedded, using low‐frequency equivalent media theories. Wave propagation was simulated for both rotationally invariant and corrugated fractures embedded in an isotropic background for one, or more than one, set of fluid‐filled and dry fractures. S‐wave splitting was studied for dipping fractures, two vertical non‐orthogonal fractures and corrugated fractures. Our modelling results confirm that S‐wave splitting can reveal the fracture infill in the case of dipping fractures. S‐wave splitting has the potential to reveal the angle between the two vertical fractures. We also notice that in the case of vertical corrugated fractures, S‐wave splitting is sensitive to the fracture infill.     

Determining finite difference weights for the acoustic wave equation by a new dispersion‐relationship‐preserving method

Numerical simulation of the acoustic wave equation is widely used to theoretically synthesize seismograms and constitutes the basis of reverse‐time migration. With finite‐difference methods, the discretization of temporal and spatial derivatives in wave equations introduces numerical grid dispersion. To reduce the grid dispersion effect, we propose to satisfy the dispersion relation for a number of uniformly distributed wavenumber points within a wavenumber range with the upper limit determined by the maximum source frequency, the grid spacing and the wave velocity. This new dispersion‐relationship‐preserving method relatively uniformly reduces the numerical dispersion over a large‐frequency range. Dispersion analysis and seismic numerical simulations demonstrate the effectiveness of the proposed method.     

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

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

A Hamiltonian Particle Method with a Staggered Particle Technique for Simulating Seismic Wave Propagation

We present a Hamiltonian particle method (HPM) with a staggered particle technique for simulating seismic wave propagation. In the conventional HPM, physical variables, such as particle displacement and stress, are defined at the center, i.e., at the same position, of each particle. As most seismic simulations using finite difference methods (FDM) are practiced with staggered grid techniques, we know the staggered alignment of space variables could improve the numerical accuracy. In the present study, we hypothesized that staggered technique could improve the numerical accuracy also in the HPM and tested the hypothesis. First, we conducted a plane wave analysis for the HPM with the staggered particles in order to verify the validity of our strategy. The comparison of grid dispersion in our strategy with that in the conventional one suggests that the accuracy would be improved dramatically by use of the staggered technique. It is also observed that the dispersion of waves is dependent on the propagation direction due to the difference in the average spacing of the neighboring two particles for the same parameters, as is usually observed in FDM with a rotated staggered grid. Next, we compared the results from the conventional Lamb’s problem using our HPM with those from an analytical approach in order to demonstrate the effectiveness of the staggered particle technique. Our results showed better agreement with the analytical solutions than those from HPM without the staggered particles. We conclude that the staggered particle technique would be a method to improve the calculation accuracy in the simulation of seismic wave propagation.     

三角网格有限元法声波与弹性波模拟频散分析

本文对声波与弹性波方程进行有限元法离散,构造有限元法频散关系的一般特征值问题,分析了时间离散格式为中心差分的三角网格有限元法声波与弹性波模拟的频散特性. 比较了三种质量矩阵即分布式质量矩阵、集中质量矩阵和混合质量矩阵对有限元法频散的影响;选取四种典型三角网格,分析了混合质量矩阵有限元（MFEM）频散的方向各向异性;数值频散、方向各向异性随插值阶数的增加逐渐减弱,当空间为三阶插值时,频散主要表现为随采样率的变化而几乎无明显方向各向异性, 其频散幅值也较小. 控制其他影响因素不变的情况下,研究了不同波速比介质中弹性波的数值频散. 最后给出了三角网格MFEM的数值耗散性.     

Interactions between topographic irregularities and seismic ground motion investigated using a hybrid FD-FE method

A hybrid method combining finite element and 4th-order finite difference techniques is developed to model SH and P-SV seismic wave propagation in a 2D elastic medium with irregular surface topography. Both the classic staggered grid finite difference scheme and the partially staggered grid scheme are tested. The accuracy of the hybrid method is studied by comparison with a semi-analytical and another numerical method. Subsequently, to study the amplification, numerical simulations of seismic wave propagation in a series of hills are carried out and compared with the single-hill case. Depending on the position of the source in relation to the topography, the ratio between the heights and lengths of the hills or the ratio between the lengths of the hills and the wavelength, the presence of several hills as opposed to a single one can increase the amplification effect due to topography. This study highlights the fact that, when evaluating topographic site effects, surrounding topography must be taken into account in addition to local topography.     

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.     

波动方程的应力-速度有限元解耦交叠格式

本基于有限差分交叠格式和解耦有限元方法的基本概念,以应力－速度为变量,提出了求解波动的应力－速度有限元解耦交叠格式,这一格式不仅时空解耦,而且为显式,它适合于线性及非线性波动问题的数值模拟,已有的应力－速度有限元交叠格式（即格子法）为本的特例。通过解析解数值检验表明,本建议的方法具有较高的精度,而格子法计算精度较低。     

用交错网格有限差分法计算三维频率域电磁响应

用交错网格有限差分法(SFD)，实现了三维频率域电磁场响应 的数值模拟. 该方法适用于任何方向的磁偶极子源. 经与解析方法、积分方程等 其他方法的计算结果对比表明，交错网格有限差分法结合散度校正和不完全乔累斯基分解预 处理的双共轭梯度迭代方法进行正演计算，速度快、精度高、结果稳定，能适应三维复杂介 质的数值模拟，为三维电磁反演奠定了基础.     

基于保辛算法的声波叠前逆时偏移

叠前逆时偏移是目前成像精度最高的地震偏移方法之一,其实现过程中的一个重要步骤是数值求解全波方程,所以快速有效求解全波方程的数值算法对逆时偏移至关重要. 四阶近似解析辛可分Runge-Kutta (NSPRK) 方法是近年发展的一种具有高效率、高精度的数值求解波动方程的保辛差分方法, 能在粗网格条件下有效压制数值频散, 从而提高计算效率, 节省计算机内存需求量. 本文利用四阶NSPRK方法构造的基本思想,发展了具有六阶空间精度的NSPRK方法,并对新的六阶NSPRK方法进行了详细的稳定性和数值频散分析,以及计算效率比较和波场模拟. 同时将该方法用于声波叠前逆时偏移中, 得到一种时间上保辛、空间具有六阶精度、低数值频散、可应用大步长进行波场延拓并能长时计算的叠前逆时偏移方法,对Sigsbee2B模型进行了偏移成像, 并和四阶NSPRK方法、传统的六阶差分方法、四阶Lax-Wendroff correction (LWC) 方法进行了对比. 数值结果表明, 基于六阶NSPRK方法的叠前逆时偏移能得到更好的成像结果, 是一种优于四阶NSPRK方法、传统的六阶差分方法、四阶LWC叠前逆时偏移的方法, 尤其是在粗网格情况下具有更明显的优越性.