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

辛维， 闫子超, 梁文全， 陈雨红， 杨长春

1 中国科学院地质与地球物理所 中国科学院油气资源研究重点实验室， 北京 100029 2 中国科学院大学， 北京 100049



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

辛维1,2， 闫子超1,2, 梁文全1， 陈雨红1， 杨长春1

1 中国科学院地质与地球物理所 中国科学院油气资源研究重点实验室， 北京 100029 2 中国科学院大学， 北京 100049

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

弹性波方程； 数值模拟； 线性化方法； 模拟退火算法

1 引言

采用有限差分法对弹性波进行数值模拟是一种常用方法(王秀明等，2004；王东等，2006；周竹生等，2007；高静怀等，2012).交错网格有限差分(FD)方法因为稳定性好(Igel et al., 1992; Dong et al., 2000)、计算效率高、所需内存小、实现简单，而广泛应用于地震正演研究(Luo and Schuster, 1990; Graves, 1996; Saenger and Bohlen, 2004)，同时是地震逆时偏移成像技术得以快速发展的基础(Kelly et al., 1976；Dablain, 1986；Levander, 1988；Liu et al., 2012；Yan and Liu, 2013).

网格频散是有限差分中至关重要的问题，直接影响着有限差分法在波动方程中的应用.网格频散是由对时间和空间偏导数的网格化造成的，相速度变成了网格间距的函数，这会导致地震波的数值相速度不等于地球介质的真实相速度，使得波场模拟的精度降低.一般来说，如果存在时间频散，则高频的相速度增大；如果存在空间频散，则高频的相速度减小(Dablain, 1986).

为了缓解或者压制网格数值频散效应，一般采用两种措施：其一是采用低阶差分格式，要求时间步长和空间间距非常小，这会造成所需内存和计算量过大而无法应用于三维的情况；其二是采用高阶差分格式(裴正林，2005；李振春等，2007)，一般情况下，时间方向采用二阶差分格式，空间方向采用高阶差分格式.高阶差分格式主要是利用等波纹准则或最小误差准则优化设计差分系数实现对空间偏导数的近似.伪谱法(朱多林和白超英，2012；唐小平等，2012)则是通过正、反傅里叶变换来实现空间偏导数的精确求解.两种方法相比，高阶空间差分格式因为具有精度适中、计算效率高的特征而得到广泛应用，特别是用于地震逆时偏移成像.Zhang和Yao(2013a，b)通过模拟退火方法确定声波方程有限差分系数，并且给出了万分之一的频散误差容限；Liu(2013)以及Yang等(2013)通过最小二乘方法确定波动方程的有限差分系数.这些方法都减小了数值频散，提高了数值模拟的计算效率.在利用有限差分法求解声波方程时，可以使用线性方法确定交错网格有限差分系数以得到较好模拟效果(Liang et al., 2013; Liang et al., 2014a,b).本文主要把声波方程的线性化方法推广到弹性波方程用于求解有限差分系数，同时提出在既定波数域内以最大最小准则为判定依据并用模拟退火算法对弹性波方程交错网格有限差分系数进行求解.然后将这两种方法与泰勒展开法、最小二乘方法进行了频散分析比较和一阶速度应力交错网格弹性波对比数值模拟，数值模拟对比表明三种方法都较泰勒展开法能有效地压制数值频散.同时试验表明：利用线性化方法直接求解有限差分系数，其计算效率要远远高于最小二乘方法和模拟退火算法.

2 空间域交错网格有限差分系数的计算方法

2.1 基本公式及已有算法

非均质的速度-应力弹性波动方程可以为(Virieux, 1986)

(1)

其中，t为时间，x和z为空间变量，(vx,vz)是速度向量，(τxx,τzz,τxz)是应力向量，λ(x,z)和μ(x,z)是Lamé系数，b=b(x,z)是密度的倒数.求解方程(1)，通常时间方向采用二阶差分格式，空间方向采用高阶差分格式.

2.1.1 泰勒展开法

一般来说，空间导数的精度通过使用高阶有限差分系数来提高，一阶导数的2M阶有限差分交错网格定义为

-u(x-mh+0.5h)]，

(2)

其中am是有限差分系数，h是空间网格间距.把平面波u(x)=u0eikx带入上面公式(2)，其中u0为常数，i为虚数单位，k为波数，并令β=kh/2得到(Yangetal., 2013)：

(3)

公式(3)两边对波数k进行泰勒展开，然后匹配k的系数，就可以得到传统的交错网格有限差分的系数，但是这种方法限定了波数(频率)范围，如果需要覆盖较广的波数域需要增加算子长度增加计算量.

2.1.2 最小二乘算法

计算有限差分系数的目的是用差分方程逼近微分方程，在更广的频率域范围之内使得公式(3)得以满足.为此，建立目标函数为 (Yang et al.,2013)

(4)

其中I是有限差分系数的非线性函数，并使I→min，fm(β)=sin((2m-1)β)，b为积分上限且小于等于π/2.目标函数的极值点在偏导数为零处，由此可以确定方程得到最小二乘交错网格有限差分系数.但是这种方法不能精确的控制误差在既定波数域内的分布，此外，由于需要计算积分，此方法求得差分算子耗时较多.

2.2 线性化方法

从公式(3)出发，可以导出求取弹性波方程交错网格有限差分系数的线性化方程.假设有M个均匀分布的波数点k(i)满足公式(3)所表示的频散关系，则由(3)立即得到线性方程组Ax=b，其中x为有限差分系数，系数矩阵A和右端项b分别为

(5)

(6)

其中β(i)(i=1,2,…,M)均匀分布于rπ/2M与rπ/2之间，其中r是用来确定波数范围的比值，使得在特定的波数范围内频散误差小于万分之一.本文定义r为

(7)

其中ku为波数的上限，f为最大频率，v为波速.解线性方程组Ax=b，就可以得到弹性波方程交错网格有限差分系数.

2.3 模拟退火算法

求解有限差分系数是函数逼近问题，这里可以选择最大最小准则为逼近准则：通过优化系数使得在同样波数覆盖范围内，数值频散带来的最大误差E最小，即可以通过方程(3)建立目标函数(5).在此准则可以使得误差可控性增强，在既定的波数域范围内最大频散误差最小，实现最大风险最小化，使整个波数域范围内的频散误差都在可以接收的范围内，其他准则达不到同样效果.该最大最小准则下的问题为非线性问题，可以利用全局寻优算法在一定准则下对系数进行遍历求解：

(8)

并使得E→min.

模拟退火算法的程序框图如图1所示.

3 频散分析

为了对上述各种方法进行比较，本节首先分别利用泰勒展开算法(TE)、模拟退火算法(SA)、最小二乘算法(LS)和线性化方法(Linearmethod(LM))求解弹性波方程的有限差分系数，然后对其进行频散分析.所求解的有限差分算子长度在M=4、6、8时的差分系数如表1所示，其中最小二乘系数来自文献(Yang et al.,2014).

表1 M=4、6、8时不同方法求得的有限差分系数Table 1 Finite difference coefficients in various methods with M=4, 6, 8

图1 模拟退火算法程序流程Fig.1 Flowchart of simulated annealing algorithm

对于有限差分算子，可以使用公式(9)衡量频散误差，表达式为

(9)

其中E(β)越接近1，误差越小，数值频散对比见图2，从图中可得到：

(1)随着算子长度的增加，各种方法都能够在更大的波数范围内保持频散关系.

(2)在一定的频散误差范围内，泰勒展开方法覆盖的波数范围最小，模拟退火和最小二乘方法覆盖的波数范围最大，线性化方法覆盖的波数范围远大于泰勒展开所覆盖的波数范围而接近于优化方法覆盖的波数范围.

(3)虽然线性化方法覆盖的波数范围略小于优化方法，但是该方法在低波数(低频)时的频散误差小于优化方法.同时考虑到震源的频谱，线性化方法得到的有限差分系数能够有效减少数值频散.

(4)利用线性化方法直接求解有限差分系数，其计算效率要远远高于最小二乘方法和模拟退火算法.以算子长度M=8为例，线性化方法的CPU计算时间在1秒钟以内，而模拟退火算法则需要2个小时左右.当算子长度长时，线性化方法的计算优势更为明显.

4 数值模拟

在本次正演模拟实验中，分别选取了简单的均匀介质模型和较为复杂的Marmousi模型进行对比计算.

4.1 均匀介质模型

在简单的均匀介质模型中，设定纵波速度为2000 m·s-1， 横波速度为1000 m·s-1，密度为1800 kg·m-3.网格大小为350×350，网格间距h=10 m，有限差分系数使用与图1b相同的有限差分系数，时间步长为1 ms.使用主频为30 Hz的雷克子波，震源放在模型的中心位置.震源加载方式不同会影响波场生成，本文震源加在正应力上，因而仅激发纵波.

不同方法得到的弹性波有限差分系数在800 ms时的波场快照如图3所示，其中x为纵波水平分量，z为纵波垂直分量.由图3 可以看出，泰勒展开方法有限差分系数的数值频散最大，模拟退火方法、最小二乘方法和线性方法的有限差分系数数值频散效应都得到了有效的压制.

为了进一步比较频散，抽取图3(a—d)z=1900 m时的切片，如图4所示.由图4可以看到，泰勒展开方法的频散最为严重，其他3种方法都有效减小了数值频散.

4.2 Marmousi模型

图 5显示了Marmousi II速度和密度模型，其纵波速度从1500 m·s-3变化到4800 m·s-3，横波速度从305 m·s-3变化到2800 m·s-3.震源函数使用主频为30 Hz的雷克子波，震源在(10 m, 3000 m)处，检波器记录深度为10 m.模拟试验中，使用了分裂的PML吸收边界(Berenger, 1994).有限差分系数使用与图2c相同的有限差分系数，空间步长为10 m，时间步长为1 ms.

图6显示了泰勒展开法(TE)、模拟退火方法(SA)、最小二乘方法(LS)和线性化方法(LM)计算有限差分系数的地震记录.图6a是泰勒展开法差分格式的地震记录，数值频散明显；图6(b—d)分别是用优化方法和线性化方法有限差分系数得到的x方向的地震记录.可以看出，数值频散得到了明显的压制.

图2 不同方法求得的有限差分系数的频散误差 (a)M=4；(b)M=6；(c)M=8.Fig.2 Dispersion errors of different FD coefficients determination methods

图7抽取了x=3000 m处的地震记录以进一步比较地震记录的数值频散.其中泰勒展开方法有限差分系数的数值频散最为明显，其他3种方法的数值频散都得到了压制，同时从图中可以看出线性方法和模拟退火方法有限差分系数的地震记录几乎重合.

5 结论

本文研究了确定弹性波方程有限差分系数的模拟退火算法和线性化方法.通过对泰勒展开方法、最小二乘方法、模拟退火方法和线性化方法交错网格有限差分系数数值频散的对比分析，可以得到以下几点认识：

(1)线性化方法、模拟退火方法和最小二乘等三种方法的数值频散均小于泰勒展开方法.

(2)在低频范围，线性化方法的数值频散要小于其他方法；同时注意到震源频谱的主能量集中在中低频，因此线性化方法能够有效去除数值频散，其数值效果与最小二乘方法以及模拟退火方法相当.

(3)利用线性化方法直接求解有限差分系数，其计算效率要远远高于最小二乘方法和模拟退火算法.

通过均匀介质模型和复杂构造模型的数值模拟及其结果的对比分析认为，线性化方法、最小二乘方法和模拟退火方法都显著地压制了数值频散(它们的地震记录几乎重合)，因而线性化方法确定交错网格有限差分系数可以用于速度-应力交错网格弹性波数值模拟和逆时偏移成像，以提高计算的精度和效率.致谢 感谢审稿专家提出的宝贵修改意见和编辑部的大力支持！

Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves.JournalofComputationalPhysics, 114(2): 185-200.

Dablain M A. 1986. The application of high-order differencing to the scalar wave equation.Geophysics, 51(1): 54-66.

图3 不同方法有限差分系数的波场快照图 (a)泰勒展开方法(TE); (b)模拟退火算法(SA); (c)最小二乘方法(LS); (d)线性方法(LM)，左为x分量，右为z分量.Fig.3 Wave field snapshots of FD coefficients determined by different methods (a) Taylor method (TE); (b) Simulated annealing algorithm (SA); (c) Least squares method (LS); (d) The linear method (LM); left for the x component, right for the z component.

图4 波场快照的切片图(x分量)Fig.4 Slices of snapshots of the wave field (x component)

图5 Marmousi速度模型 (a) P波；(b) S波；(c) 密度.Fig.5 Marmousi velocity model (a) P waves;(b) S waves;(c) Density.

图6 不同有限差分系数算法模拟得到的Marmousi模型的地震记录 (a)泰勒展开方法(TE);(b)模拟退火算法(SA);(c)最小二乘方法(LS); (d)线性方法(LM).Fig.6 Seismic records of the Marmousi model using different FD coefficient determination methods (a) Taylor method (TE); (b) Simulated annealing algorithm (SA); (c) Least squares method (LS); (d) The linear method (LM).

图7 不同方法有限差分系数地震道的对比Fig.7 Comparison of seismic traces obtained using different FD coefficient determination methods

Dong L G, Ma Z T, Cao J Z. 2000. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation.ChineseJournalofGeophysics, 43(6): 856-864.

Gao J H, He Y Y, Ma Y C. 2012. Comparison of the Rayleigh wave in elastic and viscoelastic media.ChineseJournalGeophysics(in Chinese), 55(1): 207-218, doi: 10.6038/j.issn.0001-5733.2012.01.020. Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences.BulletinoftheSeismologicalSocietyofAmerica, 86(4): 1091-1106.

Igel H, Riollet B, Mora P. 1992. Accuracy of staggered 3-D finite-difference grids for anisotropic wave propagation. //62nd Annual International Meeting, SEG, Expanded Abstracts, 1244-1246. Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms: a finite-difference approach.Geophysics, 41(1): 2-27. Levander A R. 1988. Fourth-order finite-difference P-SV seismograms.Geophysics, 53(11): 1425-1436. Li Z C, Zhang H, Liu Q M, et al. 2007. Numeric simulation of elastic wavefield separation by staggering grid high-order finite-difference algorithm.OilGeophysicalProspecting(in Chinese), 42(5): 510-515. Liang W Q, Wang Y F, Yang C C, et al. 2013. Acoustic wave equation modeling with new time-space domain finite difference operators.ChineseJournalofGeophysics, 56(6): 840-850, doi: 10.1002/cjg2.20076.

Liang W Q, Wang Y F, Yang C C. 2014a. Comparison of numerical dispersion in acoustic finite-difference algorithms.ExplorationGeophysics, doi: 10.1071/EG13072.

Liang W Q, Wang Y F, Yang C C. 2014b. Determining finite difference weights for the acoustic wave equation by a new dispersion-relationship-preserving method.GeophysicalProspecting, 63(1): 11-22. Liu H W, Li B, Liu H, et al. 2012. The issues of prestack reverse time migration and solutions with Graphic Processing Unit implementation.GeophysicalProspecting, 60(5): 906-918.

Liu Y. 2013. Globally optimal finite-difference schemes based on least squares.Geophysics, 78(4): T113-T132.

Luo Y, Schuster G. 1990. Parsimonious staggered grid finite-differencing of the wave equation.GeophysicalResearchLetters, 17(2): 155-158.Pei Z L. 2005. Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method.GeophysicalProspectingforPetroleum(in Chinese), 44(4): 308-315.

Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid.Geophysics, 69(2): 583-591.

Tang X P, Bai C Y, Liu K H. 2012. Elastic wavefield simulation using separated equations through pseudo-spectral method.OilGeophysicalProspecting(in Chinese), 47(1): 19-26.

Virieux J. 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method.Geophysics, 51(4): 889-901.

Wang D, Zhang H L, Wang X M. 2006. A numerical study of acoustic wave propagation in partially saturated poroelastic rock.ChineseJournalofGeophysics(in Chinese), 49(2): 524-532, doi: 10.3321/j.issn:0001-5733.2006.02.027.

Wang X M, Zhang H L, Wang D. 2003. Modelling of seismic wave propagation in heterogeneous poroelastic media using a high order staggered finite difference method.ChineseJournalofGeophysics(in Chinese), 46(6): 842-849.

Yan H Y, Liu Y. 2013. Acoustic VTI modeling and pre-stack reverse-time migration based on the time-space domain staggered-grid finite-difference method.JournalofAppliedGeophysics, 90: 41-52. Yang L, Yan H Y, Liu H. 2013. Least squares staggered-grid finite-difference for elastic wave modelling.ExplorationGeophysics, doi: 10.1071/EG13087. Zhang J H, Yao Z X. 2013a. Optimized finite-difference operator for

broadband seismic wave modeling.Geophysics, 78(1): A13-A18.

Zhang J H, Yao Z X. 2013b. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm.JournalofComputationalPhysics, 250: 511-526.

Zhou Z S, Liu X L, Xiong X Y. 2007. Finite-difference modelling of Rayleigh surface wave in elastic media.ChineseJournalofGeophysics(in Chinese), 50(2): 567-573, doi: 10.3321/j.issn:0001-5733.2007.02.030.

Zhu D L, Bai C Y. 2011. Review on the seismic wavefield forward modelling.ProgressinGeophysics(in Chinese), 26(5): 1588-1599, doi: 10.3969/j.issn.1004-2903.2011.05.011.

附中文参考文献

高静怀, 何洋洋, 马逸尘. 2012. 黏弹性与弹性介质中Rayleigh面波特性对比研究. 地球物理学报, 55(1): 207-218, doi: 10.6038/j.issn.0001-5733.2012.01.020.

李振春, 张华, 刘庆敏等. 2007. 弹性波交错网格高阶有限差分法波场分离数值模拟. 石油地球物理勘探, 42(5): 510-515.

裴正林. 2005. 三维各向同性介质弹性波方程交错网格高阶有限差分法模拟. 石油物探, 44(4): 308-315.

唐小平, 白超英, 刘宽厚. 2012. 伪谱法分离波动方程弹性波模拟. 石油地球物理勘探, 47(1): 19-26.

王东, 张海澜, 王秀明. 2006. 部分饱和孔隙岩石中声波传播数值研究. 地球物理学报, 49(2): 524-532, doi: 10.3321/j.issn:0001-5733.2006.02.027.

王秀明, 张海澜, 王东. 2004. 利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播. 地球物理学报, 46(6): 842-849.

周竹生, 刘喜亮, 熊孝雨. 2007. 弹性介质中瑞雷面波有限差分法正演模拟. 地球物理学报, 50(2): 567-573, doi: 10.3321/j.issn:0001-5733.2007.02.030.

朱多林, 白超英. 2011. 基于波动方程理论的地震波场数值模拟方法综述. 地球物理学进展, 26(5): 1588-1599, doi: 10.3969/j.issn.1004-2903.2011.05.011.

(本文编辑 张正峰)

Methods to determine the finite difference coefficients for elastic wave equation modeling

XIN Wei1,2, YAN Zi-Chao1,2, LIANG Wen-Quan1,CHEN Yu-Hong1, YANG Chang-Chun1

1KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China

Numerical simulation of the elastic wave equation with staggered-grid finite-difference algorithms is widely used to synthesize seismograms theoretically, and is also the basis of the reverse time migration. With some stability conditions, grid dispersion often appears because of the discretization of the time and the spatial derivatives in the wave equation. How to suppress the grid dispersion is therefore a key problem for finite-difference approaches. Different methods have been proposed to address this issue. The most commonly used methods are the high order Taylor expansion (TE) methods. In this paper, we extend the linear method for solving the acoustic wave equation to the staggered grid finite difference method for solving the elastic wave equation. We also apply the maximum-minimum criterion to measure the dispersion error when performing simulated annealing (SA) algorithm. Dispersion analysis and numerical simulation demonstrate that a linear method without iteration is nearly equal to the SA method and the least squares (LS) method in the space domain, and is better than the TE methods.For the finite difference coefficients obtained by the two methods, using homogeneous isotropic and complex structural model, we performed a numerical forward modeling and numerical dispersion analysis firstly, then compared it with the traditional Taylor expansion (TE) method and least squares(LS) method.Dispersion analysis and numerical simulation demonstrate the following conclusions: (1) With the increase of the length of the operator, various methods are able to maintain the dispersion relation in a larger wave number range. (2) The coefficients obtained by the TE method covers the minimal wave number range, coefficients from SA and LS method cover the maximal wave number range, the wave number range of linearization method is much larger than that of TE method, and is very close to that of the optimization method. (3)Although the wave number range of the linearization method is slightly less than the optimization method, but the dispersion error of this method is smaller in lower wave number (low frequency). Taking the spectrum of seismic source into account, the coefficient of linearization finite difference method is able to effectively reduce the numerical dispersion. (4) The linearization method can be used to solve finite difference coefficients directly. Its computational efficiency is much higher than the LS method and SA method.Comparison of the numerical simulation results from the uniform medium model and the complex structure model indicates that the linearization method, LS method and SA method can significantly suppress the numerical dispersion (seismograms almost coincide), thus the coefficients confirmed by the linearization finite difference method can be used to speed-stress staggered grid elastic wave modeling and time-reverse migration to improve the accuracy and efficiency of the calculation.

Elastic wave equation; Numerical simulation; Linear method; Simulated annealing

10.6038/cjg20150724.

国家自然科学基金项目(41325016和11271349)资助．

辛维，男，1986年生，中国科学院地质与地球物理研究所博士毕业，主要从事地震波正演研究和地球物理仪器研发． E-mail：332630765@qq.com

10.6038/cjg20150724

P631

2014-05-14，2015-04-23收修定稿

辛维，闫子超, 梁文全等. 2015. 用于弹性波方程数值模拟的有限差分系数确定方法.地球物理学报，58(7):2486-2495,

Xin W, Yan Z C,Liang W Q, et al. 2015. Methods to determine the finite difference coefficients for elastic wave equation modeling.ChineseJ.Geophys. (in Chinese),58(7):2486-2495,doi:10.6038/cjg20150724.