一种时空域优化的高精度交错网格差分算子与正演模拟

如何有效压制数值频散是有限差分正演模拟研究中的关键问题之一.近年来,许多学者对二阶声波方程的差分算子开展了大量的优化工作,在压制频散方面取得不错的效果.一阶压强-速度方程广泛用于研究地震波在地下变密度模型中传播规律,目前针对一阶方程的优化工作大多只是在空间差分算子上展开.本文在前人研究的基础上,推导出一阶声波方程中压强场与偏振速度场之间的解析关系,据此在传统交错网格基础上给出一种高精度的显式时间递推格式,该递推格式将时间差分与空间差分算子结合在一起,并采用共轭梯度法得到精确时间递推匹配系数,实现时空差分算子的同时优化.在编程实现算法的基础上,通过频散分析与三个典型模型测试表明:本文方法能够较为有效地压制时间频散与空间频散,提高数值计算精度;同时对复杂模型也有很好适用性.     

最优化辛格式广义离散奇异核褶积微分算子地震波场模拟

将波动方程变换至Hamilton体系,构造了一种新的保结构算法,即最优化辛格式广义褶积微分算子(OSGCD). 在时间离散上,首先引入了Lie算子设计二级二阶辛格式,基于最小误差原理得到了优化的辛格式. 在空间离散上,引入广义离散奇异核褶积微分算子计算空间微分,提出了一种有效方法优化GCD并得到了稳定的算子系数. 针对本文发展的新方法,给出了OSGCD稳定性条件. 在数值实验中,将OSGCD与多种方法比较,从精度和计算效率两方面分析了OSGCD的计算优势,计算结果也表明OSGCD长时程以及非均匀介质中地震波模拟亦具有较强能力.     

基于Remez迭代算法求解差分系数的双相介质自适应有限差分正演模拟

利用传统有限差分方法对基于Biot理论的双相介质波动方程进行数值求解时,由于慢纵波的存在,数值频散效应较为明显,影响模拟精度.相对于声学近似方程及普通弹性波方程,Biot双相介质波动方程在同等数值求解算法和精度要求条件下,其地震波场正演模拟需要更多的计算时间.本文针对Biot一阶速度-应力方程组发展了一种变阶数优化有限差分数值模拟方法,旨在同时提高其正演模拟的精度和效率.首先结合交错网格差分格式推导Biot方程的数值频散关系式.然后基于Remez迭代算法求取一阶空间偏导数的优化差分系数,并用于Biot方程的交错网格有限差分数值模拟.在此基础上把三类波的平均频散误差参数限制在给定的频散误差阈值和频率范围内,此时优化有限差分算子的长度就能自适应非均匀双相介质模型中的不同速度区间.数值频散曲线分析表明:基于Remez迭代算法的优化有限差分方法相较传统泰勒级数展开方法在大波数范围对频散误差的压制效果更明显;可变阶数的优化有限差分方法能取得与固定阶数优化有限差分方法相近的模拟精度.在均匀介质和河道模型的数值模拟实验中将本文变阶数优化有限差分算法与传统泰勒展开算法、最小二乘优化算法进行比较,进一步证明其在复杂地下介质中的有效性和适用性.     

基于新的差分结构的时-空域高阶有限差分波动方程数值模拟方法

传统的高阶有限差分波动方程数值模拟方法采用高阶差分算子近似空间偏导数,能有效抑制空间频散.然而,传统的有限差分法仅采用二阶差分算子近似时间偏导数,这使得地震波场沿时间外推的精度较低.当采用较大的时间采样间隔,传统的有限差分法模拟波场会出现明显的时间频散,甚至不稳定.本文基于新的差分结构和中心网格剖分,发展了一种空间任意偶数阶精度、时间四阶和六阶精度的时空域有限差分方法.基于对离散后的频散关系进行泰勒展开,本文推导了时空域高阶有限差分算子的差分系数.相速度分析表明时间四阶、六阶精度的差分方法能显著地减小传统时间二阶精度差分方法的时间频散.在相同的精度下与传统差分法比较,本文发展的时间四阶、六阶有限差分方法的计算效率比传统方法高.均匀和非匀均介质中的波场数值模拟实验进一步证实本文研究的时空高阶有限差分方法的优越性.     

Optimized staggered-grid finite-difference operators using window functions

The staggered-grid finite-difference (SGFD) method has been widely used in seismic forward modeling. The precision of the forward modeling results directly affects the results of the subsequent seismic inversion and migration. Numerical dispersion is one of the problems in this method. The window function method can reduce dispersion by replacing the finite-difference operators with window operators, obtained by truncating the spatial convolution series of the pseudospectral method. Although the window operators have high precision in the low-wavenumber domain, their precision decreases rapidly in the high-wavenumber domain. We develop a least squares optimization method to enhance the precision of operators obtained by the window function method. We transform the SGFD problem into a least squares problem and find the best solution iteratively. The window operator is chosen as the initial value and the optimized domain is set by the error threshold. The conjugate gradient method is also adopted to increase the stability of the solution. Approximation error analysis and numerical simulation results suggest that the proposed method increases the precision of the window function operators and decreases the numerical dispersion.     

模拟地震波场的伪谱和高阶有限差分混合方法

伪谱法是一种高效、高精度计算非均匀介质地震波传播的数值方法,但是由于它的微分算子的全局性,使得该方法不适用于分散内存的并行计算.本文将有限差分算子的局部性和伪谱法算子的高效、高精度相结合,发展基于两种方法的伪谱/有限差分混合方法.该方法在一个空间坐标方向上利用交错网格高阶有限差分算子,在另外的空间坐标方向上利用交错网格伪谱法算子,既保留了后者的高效、高精度优势,又便于在PC集群上实现并行计算.对二维模型的计算显示,混合方法能有效处理介质不连续面,在保证伪谱法计算精度的情况下,提供了一种并行计算的可能途径.     

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

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

一种优化的频率域三维声波有限差分模拟方法

为提高频率域有限差分（FD,finite-difference）正演模拟技术的计算精度和效率,基于旋转坐标系统的优化差分格式被广泛应用,但是只应用于正方形网格的情况.基于平均导数法（ADM）的优化差分格式,应用于正方形和长方形网格模拟.这些频率域有限差分算子,各自具有不同的差分格式和对应的优化系数求解表达式.本文基于三维声波方程发展了一种新的优化方法,只要给定FD模板形式,可直接构造频散方程,求取FD模板上各节点的优化系数.此方法的优点在于频率域FD算子的优化系数对应各个节点,可扩展优化其他格式.运用此优化方法,计算得到了不同空间采样间距比情况下27点和7点格式的优化系数.数值实验表明,优化27点格式与ADM 27点格式具有相同的精度,优化7点格式比经典的7点格式具有更小的数值频散.     

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

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

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

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

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.     

Acoustic VTI modeling and pre-stack reverse-time migration based on the time–space domain staggered-grid finite-difference method

Reverse-time migration (RTM) is based on seismic numerical modeling algorithms, and the accuracy and efficiency of RTM strongly depend on the algorithm used for numerical solution of wave equations. Finite-difference (FD) methods have been widely used to solve the wave equation in seismic numerical modeling and RTM. In this paper, we derive a series of time–space domain staggered-grid FD coefficients for acoustic vertical transversely isotropic (VTI) equations, and adopt these difference coefficients to solve the equations, then analyze the numerical dispersion and stability, and compare the time–space domain staggered-grid FD method with the conventional method. The numerical analysis results demonstrate that the time–space domain staggered-grid FD method has greater accuracy and better stability than the conventional method under the same discretizations. Moreover, we implement the pre-stack acoustic VTI RTM by the conventional and time–space domain high-order staggered-grid FD methods, respectively. The migration results reveal that the time–space domain staggered-grid FD method can provide clearer and more accurate image with little influence on computational efficiency, and the new FD method can adopt a larger time step to reduce the computation time and preserve the imaging accuracy as well in RTM. Meanwhile, when considering the anisotropy in RTM for the VTI model, the imaging quality of the acoustic VTI RTM is better than that of the acoustic isotropic RTM.     

Dispersion-relationship-preserving seismic modelling using the cross-rhombus stencil with the finite-difference coefficients solved by an over-determined linear system

Finite-difference modeling with a cross-rhombus stencil with high-order accuracy in both spatial and temporal derivatives is a potential method for efficient seismic simulation. The finite-difference coefficients determined by Taylor-series expansion usually preserve the dispersion property in a limited wavenumber range and fixed angles of propagation. To construct the dispersion-relationship-preserving scheme for satisfying high-wavenumber components and multiple angles, we expand the dispersion relation of the cross-rhombus stencil to an over-determined system and apply a regularization method to obtain the stable least-squares solution of the finite-difference coefficients. The new dispersion-relationship-preserving based scheme not only satisfies several designated wavenumbers but also has high-order accuracy in temporal discretization. The numerical analysis demonstrates that the new scheme possesses a better dispersion characteristic and more relaxed stability conditions compared with the Taylor-series expansion based methods. Seismic wave simulations for the homogeneous model and the Sigsbee model demonstrate that the new scheme yields small dispersion error and improves the accuracy of the forward modelling.     

Acoustic Reverse-time Migration using Optimal Staggered-grid Finite-difference Operator Based on Least Squares

Reverse-time migration (RTM) directly solves the two-way wave equation for wavefield propagation; therefore, how to solve the wave equation accurately and quickly is very important for RTM. The conventional staggered-grid finite-difference (SFD) operators are usually based on the Taylor-series expansion theory. If they are used to solve wave equation on a larger frequency content, a strong dispersion will occur, which directly affects the seismic image quality. In this paper, we propose an optimal SFD operator based on least squares to solve acoustic wave equation for prestack RTM, and obtain a new antidispersion RTM algorithm that can use short spatial difference operators. The synthetic and real data tests demonstrate that the least squares SFD (LSSFD) operator can mitigate the numerical dispersion, and the acoustic RTM using the LSSFD operator can effectively improve image quality comparing with that using the Taylor-series expansion SFD (TESFD) operator. Moreover, the LSSFD method can adopt a shorter spatial difference operator to reduce the computing cost.     

地震波传播的三维伪谱和高阶有限差分混合方法并行模拟

基于交错网格伪谱法和高阶精度有限差分方法,发展了模拟非均匀介质地震波传播的三维伪谱和有限差分混合算法.该方法在两个水平方向利用交错网格伪谱算子计算空间微分,保留了该方法高效、高精度的优势,在垂直方向采用交错网格高阶精度有限差分算子实现空间微分计算．利用有限差分方法的局部性特征,将三维计算区域在垂直方向上划分为一系列子区域,并分配给不同的处理器,实现了在并行计算机集群上的三维并行计算．通过模拟算例,与离散波数法比较,检验了该算法的精度．为了检验该方法的实用性,在64个处理器上,对三维沉积盆地模型进行了67108864个网格点的并行计算,模拟的波场主频率为1.25Hz,讨论了沉积盆地深度对三维沉积盆地地面运动的影响．      

基于加权平均导数的频率-空间域正演模拟及GPU实现

传统基于旋转坐标系的频率-空间域正演模拟方法仅适用于方形网格,而实际生产中矩形网格广泛存在,本文提出一种适用性广的正演差分算子,不仅适用于方形网格而且适用于矩形网格.通过综合运用平均导数法、加速项加权平均、模拟退火法压制频散和减少单个波长所需网格点数,从而提高算法精度和减少计算量.在该方法的基础上采用不完全LU分解作为求解Helmholtz方程的预条件,并利用图形处理器加速计算速度,很大程度上提高了频率域正演的效率.     

三维复杂构造中地震波模拟的单程波方法

复杂构造中单程波与双程波方法模拟结果的比较表明,就地震勘探中主要关心的一次反射波而言,单程波算法已具有足够的精度. 使用单程波方程将极大地减少数值计算的计算量,同时对介质的几何和物理参数建模也降低了要求. 单程波算法可视为深度偏移的“逆运算”,这样可以很好地借用已知的深度偏移方法及其程序系统. 基于计算效率和计算精度的双重考虑,本文在介质速度结构较复杂时采用显式短算子波场延拓方法,而在介质速度结构相对简单时采用分裂步相移法. 反射系数的计算中考虑了其随入射角的变化.     

基于辛-谱元-FK混合方法的保结构远震波场模拟

远震全波形层析成像能获得研究区域下方岩石圈乃至地幔过渡带高分辨率速度结构,是研究地球深部构造与动力学过程的有效工具.该类方法需以高精度及长时程远震波场正演模拟为基础,这为设计高精度长时程稳定的正演算法带来了挑战.在此背景之下,本文提出了一种适用于远震波场模拟的保结构算法.该方法采用谱元法（SEM）对研究区域进行空间离散,在不考虑耗散项情况下,将空间离散后的常微分方程变换为哈密顿系统形式,采用保辛分部龙格-库塔方法数值求解.在三级保辛分部龙格-库塔算法基础上添加额外空间离散项,得到修正辛算法.本文将该时间-空间全离散形式称为修正辛-谱元法（SSEM）,并将SSEM算法与频率波数域（FK）方法结合,发展了可模拟高频远震波场在局域模型内传播的SSEM-FK混合方法.该方法结合了FK方法模拟层状介质中平面波传播的高效性和SSEM计算复杂介质中弹性波传播的精确性.数值实验表明,SSEM-FK能够准确模拟高频远震波场在研究区域内的传播,结合该方法在计算效率上的优势,可为高效、高精度的远震全波形层析成像打下基础.     

最小二乘傅立叶有限差分偏移

一般偏移算法是用反演算子通过解析方法求解.最小二乘偏移方法采用另一种思路,即采用数值方法,通过解一个线性离散反问题来索求解.这样我们试着寻找一个模型匹配地震数据并能表现出其某些特点来得到偏移图像.最小二乘法能减少偏移赝像,得到更精确的偏移效果.Kirchhoff算子在最小二乘偏移方法中应用较广,但需要较多的迭代次数,而且具有Kirchhoff偏移的缺点.本文把最小二乘方法运用到基于波长延拓的波动方程偏移方法中,为提高最小二乘偏移的效率,可采用效率较高的正传播算子和反传播算子.我们利用效率较高,能适应剧烈横向变速的傅立叶有限差分正传播和反传播算子来做叠后最小二乘偏移.数值实例表明,通过少数的共轭梯度法迭代,就能得到与真实模型差别很微小的偏移效果.对于傅式变换我们采用了数值软件FFTW,其变换速度比常规FFT算法一般要快六倍以上,进一步提高了效率.本文算法很容易在并行机上实现,这些特点在处理大型数据时大有裨益.     

A new kind of optimal second-order symplectic scheme for seismic wave simulations

Here we introduce generalized momentum and coordinate to transform seismic wave displacement equations into Hamiltonian system. We define the Lie operators associated with kinetic and potential energy, and construct a new kind of second order symplectic scheme, which is extremely suitable for high efficient and long-term seismic wave simulations. Three sets of optimal coefficients are obtained based on the principle of minimum truncation error. We investigate the stability conditions for elastic wave simulation in homogeneous media. These newly developed symplectic schemes are compared with common symplectic schemes to verify the high precision and efficiency in theory and numerical experiments. One of the schemes presented here is compared with the classical Newmark algorithm and third order symplectic scheme to test the long-term computational ability. The scheme gets the same synthetic surface seismic records and single channel record as third order symplectic scheme in the seismic modeling in the heterogeneous model.