一种全局优化的隐式交错网格有限差分算法及其在弹性波数值模拟中的应用

在数值模拟中,隐式有限差分具有较高的精度和稳定性.然而,传统隐式有限差分算法大多由于需要求解大型矩阵方程而存在计算效率偏低的局限性.本文针对一阶速度-应力弹性波方程,构建了一种优化隐式交错网格有限差分格式,然后将改进格式由时间-空间域转换为时间-波数域,利用二范数原理建立目标函数,再利用模拟退火法求取优化系数.通过对均匀模型以及复杂介质模型进行一阶速度-应力弹性波方程数值模拟所得单炮记录、波场快照分析表明:这种优化隐式交错网格差分算法与传统的几种显式和隐式交错网格有限差分算法相比不但降低了计算量,而且能有效的压制网格频散,使弹性波数值模拟的精度得到有效的提高.     

可变网格与局部时间步长的交错网格高阶差分弹性波模拟

交错网格高阶差分方法是一种在保持效率的前提下提高弹性波模拟精度的有效方法.本文将可变空间网格与变化的时间步长技术引入到交错网格高阶差分弹性波模拟中,提出一种空间网格可任意奇数倍变化与时间步长任意变化的交错网格高阶差分弹性波模拟方法.一系列数值试验表明,该方法能够在保证模拟精度的同时,通过有效降低空间与时间维度上的过采样来显著提高弹性波模拟的效率.同时,该方法还能够精细刻画含孔缝洞介质以及横向变化剧烈介质的局部细微结构,减小弹性波模拟误差,提高介质细微结构处的弹性波传播模拟精度.     

拟空间域弹性波方程交错网格有限差分逆时偏移

弹性波逆时偏移是一种高精度的复杂构造地震成像方法.然而,在传统的基于矩形网格离散化的逆时偏移中,介质界面通常会产生畸变.另外,因使用双程波动方程进行波场延拓,其产生的反射波会在成像过程中产生偏移假象.为解决这些问题,本文提出了一种拟空间域弹性波方程高阶交错网格有限差分格式,并给出了差分格式的稳定性条件,进而实现了高精度的拟空间域弹性波方程有限差分逆时偏移.模型实验表明,若在计算拟空间域采样间隔时引入速度界面信息,则拟空间域弹性波方程高阶交错网格有限差分逆时偏移能够避免常规弹性波方程逆时偏移中弯曲界面形态畸变问题;此外基于该方法进行波场延拓时可有效压制弯曲界面的假散射现象,并能有效压制层间反射波,因此可以减少剖面上的偏移假象,从而显著提高成像的质量.

     

二维弹性及粘弹性TTI介质中地震波场数值模拟:四种不同网格高阶有限差分算法研究

本文利用交错网格、辅助网格、旋转交错网格、同位网格有限差分方法分别模拟了二维弹性TTI介质和二维黏弹性TTI介质中的地震波传播.在稳定性条件内,选用不同的网格间距及时间间隔,通过波场快照、合成理论地震图较为系统分析对比了这四种不同网格有限差分数值模拟在计算精度、CPU时间、相移、频散、以及保幅方面的优缺点.数值模拟结果表明:1)这四种不同网格有限差分算法都是很好的波场数值模拟算法;2)就CPU计算时间而言,旋转交错网格有限差分算法的计算效率最高;3)从计算精度来看,同位网格有限差分的计算精度最高;4)从振幅保护方面来看,四种网格的保护振幅的能力相当;5)相移方面,当网格间距增大时,交错网格和旋转交错网格有可能出现相移现象;6)频散方面,同位网格的频散现象不明显.     

波场模拟有限差分参数优化选取

有限差分方法(Finite-difference Method, FD)广泛用于地震波场数值模拟, 但其存在固有的数值频散问题, 影响模拟的计算效率和数值精度.本文主要研究了有限差分方法的空间数值频散误差和网格划分精度以及差分算子的关系, 基于计算量最小准则, 提出了最优化有限差分参数选取流程, 为有限差分数值模拟参数选取提供理论指导.本文主要工作包括: (1) 提出了空间数值频散正变换过程(Forward Space Dispersion Transform, FSDT)方法, 该方法可以高效模拟出不同网格划分精度(波长采样点数)的带有空间数值频散的波场; (2) 提出了波场空间数值频散误差衡量准则, 可以定量地判断出数值模拟导致的波形频散程度, 选取合适的频散误差阈值; (3) 研究了给定空间数值频散误差阈值下, 差分算子系数、差分算子阶数、网格划分精度与计算量之间的关系.文中基于雷米兹交换方法(Remez Exchange Method, RE)和泰勒级数展开方法(Taylor-series Expansion Method, TE)的差分系数, 在空间数值频散误差阈值0.01时, 数值模拟了不同差分算子阶数、网格划分精度与计算量的关系, 并给出了有限差分参数选取的参考值.

     

Some unintended features of elastic finite-difference models

While finite-difference methods have been used extensively for many years to model wave propagation in elastic media, some of the more subtle effects observable in such models are very inadequately documented in the geophysical literature, especially in regard to their practical numerical consequences. In addition to the intended travelling waves, and the undesirable exponential instability revealed by the von Neumann test, typical second-order-time finite-difference equations also support drifting linear solutions, as can be verified, both theoretically and by numerical experiment. The necessity of these solutions, and their relationship to the incompleteness of the set of travelling-wave eigenfunctions of the finite-difference operator, can be exposed by a matrix-based analysis, and exact expressions for them can be obtained by using standard algebraic techniques. A further peculiarity of the finite-difference formulation is numerical anisotropy, which emerges in a grid of more than one spatial dimension, even when the modelled medium is intended to be isotropic. This anisotropy can be explained and quantified in terms of the exact eigenfunction solutions to the finite-difference equation, which, it is found, can be obtained in a simple, closed form, for a typical modern 3D staggered scheme.     

Poroelastic finite-difference modeling for ultrasonic waves in digital porous cores

Scattering attenuation in short wavelengths has long been interesting to geophysicists. Ultrasonic coda waves, observed as the tail portion of ultrasonic wavetrains in laboratory ultrasonic measurements, are important for such studies where ultrasonic waves interact with small-scale random heterogeneities on a scale of micrometers, but often ignored as noises because of the contamination of boundary reflections from the side ends of a sample core. Numerical simulations with accurate absorbing boundary can provide insight into the effect of boundary reflections on coda waves in laboratory experiments. The simulation of wave propagation in digital and heterogeneous porous cores really challenges numerical techniques by digital image of poroelastic properties, numerical dispersion at high frequency and strong heterogeneity, and accurate absorbing boundary schemes at grazing incidence. To overcome these difficulties, we present a staggered-grid high-order finite-difference (FD) method of Biot’s poroelastic equations, with an arbitrary even-order (2L) accuracy to simulate ultrasonic wave propagation in digital porous cores with strong heterogeneity. An unsplit convolutional perfectly matched layer (CPML) absorbing boundary, which improves conventional PML methods at grazing incidence with less memory and better computational efficiency, is employed in the simulation to investigate the influence of boundary reflections on ultrasonic coda waves. Numerical experiments with saturated poroelastic media demonstrate that the 2L FD scheme with the CPML for ultrasonic wave propagation significantly improves stability conditions at strong heterogeneity and absorbing performance at grazing incidence. The boundary reflections from the artificial boundary surrounding the digital core decay fast with the increase of CPML thicknesses, almost disappearing at the CPML thickness of 15 grids. Comparisons of the resulting ultrasonic coda Q _sc values between the numerical and experimental ultrasonic S waveforms for a cylindrical rock sample demonstrate that the boundary reflection may contribute around one-third of the ultrasonic coda attenuation observed in laboratory experiments.     

基于自适应优化有限差分方法的全波VSP逆时偏移

与地面地震资料相比,VSP资料具有分辨率高、环境噪声小及能更好地反映井旁信息等优点.常规VSP偏移主要对上行反射波进行成像,存在照明度低、成像范围受限等问题.为了增加照明度、拓宽成像范围、提高成像精度,本文采用直达波除外的所有声波波场数据(全波),包括一次反射波、多次反射波等进行叠前逆时偏移成像.针对逆时偏移中的四个关键问题,即波场延拓、吸收边界条件、成像条件及低频噪声的压制,本文分别采用自适应变空间差分算子长度的优化有限差分方法(自适应优化有限差分方法)求解二维声波波动方程以实现高精度、高效率的波场延拓,采用混合吸收边界条件压制因计算区域有限所引起的人工边界反射,采用震源归一化零延迟互相关成像条件进行成像,采用拉普拉斯滤波方法压制逆时偏移中产生的低频噪声.本文对VSP模型数据的逆时偏移成像进行了分析,结果表明:自适应优化有限差分方法比传统有限差分方法具有更高的模拟精度与计算效率,适用于VSP逆时偏移成像;全波场VSP逆时偏移成像比上行波VSP逆时偏移的成像范围大、成像效果好;相对于反褶积成像条件,震源归一化零延迟互相关成像条件具有稳定性好、计算效率高等优点.将本文方法应用于某实际VSP资料的逆时偏移成像,进一步验证了本文方法的正确性和有效性.     

地震波有限差分模拟综述

本文从有限差分法数值模拟技术的各个方面对地震波有限差分模拟的发展和现状进行了论述.波场的数值模拟技术是认识地震波传播规律,检验各种处理方法正确性的重要工具,地震波的数值模拟是地震波传播规律研究的必要手段,贯穿于地震资料的采集、处理、解释的整个过程中.有限差分法数值模拟技术相对于射线方法具有更高的精度,同时比有限元方法计算量小,因此在实际应用中占很重要的地位.     

可变网格与局部时间步长的高阶差分地震波数值模拟

提高计算精度与效率是所有地震波正演方法所追求的目标.本文通过将变化的空间网格与变化的时间步长技术相结合,提出一种空间网格大小与时间步长均可任意变化的高阶有限差分模拟方法.一系列数值试验表明,该方法在保证模拟精度的同时,显著提高了模拟的效率.这种可变空间网格与局部时间步长的模拟方法,能够精细刻画含孔缝洞介质以及横向变化剧烈介质的微小结构,减小地震波模拟误差,提高介质细微结构情况下的地震波传播模拟精度与效率.     

Using the pseudospectral method on curved grids for 2D elastic forward modelling1

When applying the conventional Fourier pseudospectral method (FSM) on a Cartesian grid that has a sufficient size to propagate a pulse, spurious diffractions from the staircase representation of the curved interfaces appear in the wavefield. It is demonstrated that these non-physical diffractions can be eliminated by using curved grids that conform to all the interfaces of the subsurface. Methods for solving the 2D acoustic wave equation using such curved grids have been published previously by the authors. Here the extensions to the full 2D elastic wave equations are presented. The curved grids are generated by using the so-called multiblock strategy which is a well-known concept in computational fluid dynamics. In principle the sub-surface is divided into a number of contiguous subdomains. A separate grid is generated for each subdomain patching the grid lines across domain boundaries to obtain a globally continuous grid. Using this approach, even configurations with pinch outs can be handled. The curved grid is taken to constitute a generalized curvilinear coordinate system. Thus, the elastic equations have to be written in a curvilinear frame before applying the numerical scheme. The method implies that twice the number of spatial derivatives have to be evaluated compared to the conventional FSM on a Cartesian grid. However, it is demonstrated that the extra terms are more than compensated for by the fewer grid points needed in the curved approach.     

Propagation of Kelvin waves along irregular coastlines in finite-difference models

In this paper, we examine the behavior of internal Kelvin waves on an f-plane in finite-difference models using the Arakawa C-grid. The dependence of Kelvin wave phase speed on offshore grid resolution and propagation direction relative to the numerical grid is illustrated by numerical experiments for three different geometries: (1) Kelvin wave propagating along a straight coastline; (2) Kelvin wave propagating at a 45° angle to the numerical grid along a stairstep coastline with stairstep size equal to the grid spacing; (3) Kelvin wave propagating at a 45° angle to the numerical grid along a coarse resolution stairstep coastline with stairstep size greater than the grid spacing. It can be shown theoretically that the phase speed of a Kelvin wave propagating along a straight coastline on an Arakawa C-grid is equal to the analytical inviscid wave speed and is not dependent on offshore grid resolution. However, we found that finite-difference models considerably underestimate the Kelvin wave phase speed when the wave is propagating at an angle to the grid and the grid spacing is comparable with the Rossby deformation radius. In this case, the phase speed converges toward the correct value only as grid spacing decreases well below the Rossby radius. A grid spacing of one-fifth the Rossby radius was required to produce results for the stairstep boundary case comparable with the straight coast case. This effect does not appear to depend on the resolution of the coastline, but rather on the direction of wave propagation relative to the grid. This behavior is important for modeling internal Kelvin waves in realistic geometries where the Rossby radius is often comparable with the grid spacing, and the waves propagate along irregular coastlines.©1998 Published by Elsevier Science Limited. All rights reserved     

2D finite-difference elastic wave modelling including surface topography1

A 2D numerical finite-difference algorithm accounting for surface topography is presented. Higher-order, dispersion-bounded, cost-optimized finite-difference operators are used in the interior of the numerical grid, while non-reflecting absorbing boundary conditions are used along the edges. Transformation from a curved to a rectangular grid achieves the modelling of the surface topography. We use free-surface boundary conditions along the surface. In order to obtain complete modelling of the effects of wave propagation, it is important to account for the surface topography, otherwise near-surface effects, such as scattering, are not modelled adequately. Even if other properties of the medium, for instance randomization, can improve numerical simulations, inclusion of the surface topography makes them more realistic.     

傅里叶有限差分法三维波动方程正演模拟

傅里叶有限差分(FFD)法兼有相位屏法和隐式有限差分法二者的优势,能够处理复杂地质构造中的波传播问题,但在三维情形下,算子的双向分裂会引起明显的方位各向异性误差.本文用Fourier变换计算双向分裂过程中的高阶交叉项,消除了方位各向异性误差.该方法充分利用了FFD法在双域实现的算法结构,明显减少了由于引入误差校正所带来的计算量.将该方法应用于修改后的三维French模型的地震正演问题,并将得到的叠后记录、单炮记录同全波有限差分法的模拟结果进行对比,结果证实了该方法对一次反射波具有较高的模拟精度,在内存需求和计算效率方面则具有更大的优势.     

A staggered-grid high-order finite-difference modeling for elastic wave field in arbitrary tilt anisotropic media

Introduction The real Earth usually presents anisotropy. Therefore, it is of theoretical and practical sig- nificance for many fields as oil and gas, seismic exploration and production, earthquake prediction, detection of deep structure and so on to study on seismic wave theory, numerical simulation method and its applications in the anisotropic media (Crampin, 1981, 1984; Crampin et al, 1986; Hudson et al, 1996; Liu et al, 1997; Thomsen, 1986, 1995; TENG et al, 1992; HE and ZHANG, 1996)…     

结合边界条件的截断高阶隐式有限差分方法(英文)

有限差分方法广泛应用于求解许多科技领域所涉及的偏微分方程,高阶显式有限差分方法通常用来提高求解精度,已经提出的高阶隐式有限差分方法和截断高阶显式有限差分方法可用来进一步提高模拟精度而不增加计算量。本文首先计算了针对常规网格上的一阶导数和二阶导数、交错网格上的一阶导数的有限差分系数,发现高阶隐式有限差分系数中存在一些小的系数。频散分析结果表明:忽略这些小的差分系数能够近似维持有限差分的精度,但是显著减小了计算量。然后,引入镜像对称边界条件来提高隐式有限差分方法的精度和稳定性,采用混合吸收边界条件来减小来自模型边界所不需要的反射。最后,给出了针对均匀和非均匀介质模型的弹性波模拟例子,表明了本文方法的优点。     

二维频率空间域25点优化系数差分格式弹性波数值模拟(英文)

频率空间域地震波数值模拟具有独特的优势:可以同时模拟多源的波传播、每个频率之间独立并行地计算、计算频带选择灵活、不存在累计误差、容易模拟粘弹性介质中地震波传播.但是该方法的最大瓶颈是对于计算机内存的巨大需求.我们使用压缩存储系数矩阵的方法,极大地减少了计算机内存的需求量.同时为了减少短筹分算子的数值频散,引用了频率空间域25点弹性波波动方程的差分格式,并使用了最小二乘意义下求出的优化差分系数.为了克服边界反射,采用了最佳匹配层吸收边界条件.数值模拟试验证明:用压缩存储系数矩阵及优化差分系数的频率空间域25点差分格式进行弹性波正演模拟,可以减少数值频散,提高计算精度.使用较大的网格间距,降低计算机内存需求,并保持较高的计算效率.该正演方法为后续弹性波偏移和弹性参数反演提供较好的基础.     

Efficient finite-difference modelling of seismic waves using locally adjustable time steps

Conventional finite-difference modelling algorithms for seismic forward modelling are based on a time-stepping scheme with a constant (global) time step. Large contrasts in the velocity model or in the spatial sampling rate cause oversampling in time for some regions of the model. The use of locally adjustable time steps can save large amounts of computation time for certain modelling configurations. The computation of spatial derivatives across the transition zone between regions of the model with different temporal sampling requires the definition of the wavefield at corresponding time levels on both sides of the transition zone. This condition can be obtained by extrapolation in time, which is inaccurate, or by multiple time integration in the transition zone. The error in the latter solution is of the same order as the conventional time-stepping scheme because both methods are based on the same iteration formula. The technique of multiple time integration simply requires the use of different sizes of time step. It is applicable only for certain factors of variation of the time step.     

求解弹性波方程的辛RKN格式

将弹性波方程变换至Hamilton体系, 构造适用于弹性波模拟的高效显式二阶辛Runge-Kutta-Nyström (RKN)格式, 运用根数理论得到此格式的阶条件方程组. 通过给定系数的限定条件, 得到方程的对称解. 为了使时间离散误差达到极小,提出数值频率与真实频率比较,通过Taylor展开,得到关于辛系数的限定方程,求解方程组得到最小频散辛RKN格式. 对比分析时间演进方程的稳定性,得到使库朗数达到极大值的限定方程,求解方程组得到最稳定辛RKN格式. 发现此两种格式为同一格式. 新得到的辛RKN格式不依赖于空间离散方法,为了对比的需要,选取有限差分法进行空间离散. 在频散、稳定性分析中,与常见辛格式对比,从理论上分析了本文提出的格式在数值频散压制、稳定性提升等方面的优势, 数值实验进一步证实了理论分析的正确性.