地震波数值模拟的褶积微分算子法与伪谱法的对比分析①

刘红艳， 陈宇坤， 刘 芳

(天津市地震局，天津 300201)

地震波数值模拟的褶积微分算子法与伪谱法的对比分析①

刘红艳， 陈宇坤， 刘 芳

(天津市地震局，天津 300201)

将基于Forsyte广义正交多项式的褶积微分算子法运用于复杂非均匀介质地震波场模拟中，并将计算结果与伪谱法计算结果进行分析比较。通过二者的计算时间对比发现：在同样的计算条件下，褶积微分算子法的采样时间始终小于伪谱法，这是其进行地震波数值模拟的一个明显优势。通过波场快照的对比，褶积微分算子法的模拟结果与伪谱法数值模拟结果的频散效应相当，可为地震波场的值计算提供一种新的选择。

褶积微分算子； 伪谱法； 地震波传播； 数值模拟

0 引言

地震波场的数值算法和正演模拟历来都是地震学研究的热门领域。迄今为止，地震波场正演模拟的数值算法已有很多，如射线追踪法[1-3]，积分方程法[1-6]，或基于波动方程直接解法的有限差分法[7-12]、伪谱法[12-17]、有限元法[18-19]和谱元法[20-22]，但大家很容易忽略如细胞自动机法[23]、辛几何算法[24]、褶积微分算子法[25-31]、辛格式奇异核褶积微分算子[32-34]和Shannon奇异核褶积微分算子[35],这些也非常重要的数值计算方法。实际上它们都是吸收计算数学的一些思想，并成功地将其应用于地震波场的数值模拟。

现有的各种地震波场正演数值模拟方法都存在自身的优越性与局限性。本文的主要目标就是研究一种已推出的、快速的、高精度的地震波场模拟方法(基于Forsyte广义正交多项式的迭积微分算子法)，并将该算子所构造的正演模拟算法运用于复杂非均匀介质模型的地震波场模拟中，同时将计算结果与伪谱法的计算结果做比较，以分析该方法的优越性。

1 基本方程

二维均匀介质中，一阶速度-应力弹性波动方程为

(1)

(2)

在二维各向异性介质中，一阶速度-应力弹性波动方程是：

(3)

2 Forsyte广义正交多项式褶积微分算子求导的基本原理

Forsyte 多项式是一个广义正交多项式，其插值函数可写为[29]：

(4)

其中，

p0=1

f(xi)为被插值函数f(x)在点xi处的值。式(4)中的p0(x)，…，pj+1(x)定义为Forsyte多项式系统。

对式(4)中的x求导，可得：

(5)

其中

Forsyte 多项式微分算子可写为

(6)

将式(6)离散化可得

(7)

式中，i为采样指标；Δx为沿着x轴的采样间隔。在实际应用中须将微分算子截成短算子，势必引起Gibbs现象。另外，多项式的引入还将引起Runge现象。为了消除这些现象，必须采用窗函数以截断长微分算子。本文采用的Gaussian窗函数为

(8)

其中，mx为单边截断长度的采样数；c为常数；a(0.1≤a≤0.75)为衰减因子。将式(5)用式(6)截断并锯齿化后，可得实用的一阶迭积微分算子：

(9)

通过对算子长度的调节及系数的优化，可同时兼顾波场解的全局信息与局部信息。本文运用算子长度为9 点的一阶褶积微分算子求解弹性波场的一阶速度-应力方程。9 点迭积算子的最优权系数为：-0.002 34，0.008 74，-0.027 4，0.085，0.0，-0.085，0.027 4，-0.008 74及0.002 34。可以看出算子的系数是反对称的。图1是9点微分算子振幅随采样点的变化曲线,从图中可以看出算子振幅随采样点很快地衰减。

图1 9点褶积微分算子振幅随采样点的变化曲线Fig.1 Amplitude versus sampling point for nine-point convolutional differentiator

3 伪谱法求导的基本原理

伪谱法是地震波场数值模拟诸多方法中使用较多的一种方法，其基本原理就是在空间域中通过FFT将空间域变换到频率-波数域中，然后再通过傅里叶求导方式来进行求导，从而避免了空间域中的直接求导。

对于二维波动方程，水平分量u的傅里叶变换为

(10)

由傅里叶变换的微分性质

(11)

对式(11)进行反傅里叶变换得

(12)

传统伪谱法在求解时间导数的时候采用二阶中心差分的形式，其求解方式为：

(13)

(14)

其中，Δt是时间采样率。

4 数值模拟实例及分析

利用三个不同特性和几何结构的模型对褶积微分算子法和经典伪谱法作比较，检验两种方法在刻画波场精细细节和抑制频散方面的效果。

4.1 大角度倾斜界面模型

图2所示为大角度倾斜界面模型。网格大小为256×256，网格空间步长为10 m，时间步长为1 ms，震源位于格点 (140,130) 处，介质其他参数见表1。

表1 大角度倾斜界面模型参数

图2 大角度倾斜界面模型示意图Fig.2 Diagram of large-angle dipping interface model

①直达波P1；②直达波S1；③同类反射波P1P1；④转换反射波P1S1；⑤转换反射波S1P1；⑥同类反射波S1S1； ⑦同类透射波P1P2；⑧转换透射波P1S2；⑨同类透射波S1S2

图3为波场快照。可以看到：两层介质波速差异很大，使其波阻抗差别也很大，导致两层界面之间出现了很强且十分明显的反射和透射。同样，在介质中由于波传播时发生分裂(即分裂为横波和纵波)，且纵波波速较横波大，可以判断外层圆为纵波快照，内层圆为横波快照。横波和纵波在与波传播方向垂直与平行方向上的差异同单一均匀介质相似。但当波通过界面时，在其左侧可以看到四个圈层，由左向右分别是透射P波和转换S波，以及转换P波和透射S波；在右侧亦对应地出现反射P波和S波以及转换P波和透射S波。由它们的反射点，可以很清晰地辨别两层几面所在。比较二者的波场快照可以发现，二者的频散效应基本相当。但是从图4所示的时间曲线上可以发现，在计算相同空间步数的时候，褶积微分算子法的用时要明显少于伪谱法，这是由于伪谱法是一种全局算法，而褶积微分算子法是一种局部算法。

图4 褶积微分算子法与伪谱法计算耗时比较Fig.4 Time-consuming comparison of calculating by the convolution differentiator method and pseudo-spectral method

4.2 Marmousi经典模型

图5(a)为原始Marmousi速度模型示意图，图5(b)、(c)分别是利用褶积微分算子法和伪谱法得到的波场快照。从波场快照来看，震源发出的地震波在介质界面和突变点处产生的反射和绕射等现象非常明显，并且两种方法得到的波场快照的频散程度基本相当。说明基于Forsyte广义正交多项式褶积微分算子法的地震数值模拟方法可以很好地模拟地震波在速度变化比较剧烈的介质中的传播。

图5 Marmousi速度模型及某时刻的波场快照Fig.5 The Marmousi velocity model and wave field snapshot at a given time

图6所示为褶积微分算子法与经典伪谱法原始Marmousi模型同一点地震图的比较。从图中可以看出，除了振幅有微小的差别外，相位基本上一致，这也从另一个侧面说明了褶积微分算子法的有效性。

4.3 各向异性介质模型

采用Thomsen[30]选出的实际页岩岩样作为检验模型。模型网格大小为128×128，网格空间步长为10 m，时间步长为1 ms，震源为主频20 Hz的Ricker子波，位置位于介质的中心位置。介质的弹性参数见表2。

表2 各向异性介质模型参数

图6 褶积微分算子法与经典伪谱法原始 Marmousi 模型同一点振动图的比较Fig.6 Comparison of vibration at the same point of original Marmousi model using the convolutional differentiator method and psendo-spectra method

图7 各向异性页岩介质中的波场快照Fig.7 Snapshot of wave field in anisotropic shale media

比较二者的波场快照(图7)可以发现，二者的频散效应基本相当。两种方法得到的波场快照都清晰地刻画了波前面的形态。我们知道，各向异性介质模型中的波场传播特征对确定地层构造特征和类型具有重要的理论意义与应用前景。

5 与伪谱法的计算时间理论分析及与有限差分法及伪谱法的理论精度分析

5.1 与伪谱法的计算时间理论分析

假设二维空间计算区域的大小为N×N，褶积微分算子的长度为M。对声波方程进行数值计算时，使用伪谱法需要(N2log2N)次复数乘法；使用褶积微分算子法只需要M×N2次实数乘法。显然两者理论计算时间比为M/(4log2N)，只要M<(4log2N)，则褶积微分算子法将比伪谱法的计算速度快[25]。

5.2 与有限差分法及伪谱法的理论精度分析

在波动方程基础上，检验数值模拟方法准确性的最精确实验是均匀介质中脉冲响应的计算，因为在数值计算中单脉冲能够产生最宽的频段。Zhou等[25，31]分别计算了不同方案所对应的零炮检距脉冲响应，计算结果表明，伪谱法可以产生很好的效果，子波延续时间很短，并且没有二次扰动；四阶有限差分法的准确性要相对差一点，子波开始时是正确的，但由于网格频散，子波延续时间长，并且有二次扰动产生。九点褶积微分算子法计算结果接近于伪谱法，比四阶有限差分法要好得多[28]。

6 结论

本文所使用的方法是将褶积微分算子应用于地震波场数值模拟的又一次成功尝试。通过引入高斯窗函数，较大程度地压制了算子截断引起的吉普斯效应。三种不同地质模型理论计算结果表明，该方法计算速度快，精度高，是继戴志阳[27]发展的褶积微分算子法之后的又一种弹性波数值计算的有效工具。将本方法与伪谱法进行比较，发现本文的褶积微分算子法对数值频散的抑制能力与伪谱法相当，但计算时间明显占据优势，可以更高效地模拟地震波在复杂介质中的传播过程。

)

[1]ChapmanCH，ShearerPM.RayTracinginAzimuthallyAnisotropicMedia-Ⅱ ——Quasi-shearWaveCoupling[J].GeophysicalJournalInternational，1989，96(1)：65-83.

[2] Li Y G，Leary P C，Aki K.Seismic Ray Tracing for VSP Observations in Homogeneous Fractured Rock at Oroviue[J].EOS Transaction America Geophysical Union，1986，67(44)：1117.

[3] Cerveny V，Firbas P.Numerical Modeling and Inversion of Travel Times of Seismic Body Waves in Inhomogeneous Anisotropic Media[J].Science & Mathematic Geophysical Journal International，1984，76(1)：41-51.

[4] Pao Y H，Varatharajulu V.Huygen’s Principle Radiation Conditions and Integral Formulas for the Scattering of Elastic Waves[J].The Journal of the Acoustical Society of America，1976，59(6)：1361-1371.

[5] Bouchon M.Diffraction of Elastic Waves by Cracks or Cavities Using the Discrete Wave Number Method[J].The Journal of the Acoustical Society of America，1987，81：1671-1676.

[6] Fu L Y，Bouchon M.Discrete Wavenumber Solutions to Numerical Wave Propagation in Piecewise Heterogeneous Media-Ⅰ——Theory of Two-dimensional SH Case[J].Geophysical Journal International，2004，157：481-498..

[7] Alterman Z，Karal C.Propagation of Seismic Waves in Layered Media by Finite Difference Methods[J].Bulletin of the Seismological Society of America，1968，58(1)：367-398..

[8] Virieux J.S H Wave Propagation in Heterogeneous Media：Velocity-stress Finite Difference Method[J].Geophysics，1984，49(11)：1933-1957.

[9] Virieux J，P-SV Wave Propagation in Heterogeneous Media：Velocity-stress Finite-difference Method：Shear Waves[J].Geophysics，1986，51(4)： 889-901.

[10] Oprsal I，Zahradnik J.Elastic Finite Difference Method for Irregular Grids[J].Geophysics，1999，64(1)：240-250.

[11] Carcione J M，Herman G C，Kroode A P E.Seismic Modeling[J].Geophysics，2002，67(4)：1304-1325.

[12] 严武建，王彦宾，石玉成.基于伪谱和有限差分混合方法的兰州盆地强地面运动二维数值模拟[J].地震工程学报，35(2)：302-310.YAN Wu-jian，WANG Yan-bin，SHI Yu-cheng.2D Simulation of the Strong Ground Motion in Lanzhou Basin with Hybrid PSM/FDM Method[J].Northwestern Seismological Journal，2013，35(2)：302-310.(in Chinese)

[13] Gazdag J.Modeling of the Acoustic Wave Equation with Transform Method[J].Geophysics，1981，46：854.

[14] Kosloff D，Baysal E.Forward Modeling by a Fourier Method[J].Geophysics，1982，47(10)：1402-1412.

[15] Kosloff D，Reshef M，Loewenthal D.Elastic Wave Calculation by the Fourier Method[J].Bulletin of the Seismological Society of America，1984，74(3)：875-891.

[16] Fornberg B.The Pseudo-spectral Method：Comparisons with Finite Differences for the Elastic Wave Equation[J].Geophysics，1987，52(4)：483-501.

[17] Reshef M，Kosloff D.Applications of Elastic Forward Modeling to Seismic Interpretation[J].Geophysics，1985，50：1266-1272.

[18] 周辉，徐世浙，刘斌，等.各向异性介质中波动方程有限元法模拟及其稳定性[J].地球物理学报，1997，40(6)：833-841.ZHOU Hui，XI Shi-zhe，LIU Bin，et a1.Modeling of Wave Equations in Anisotropic in Anisotropic Media Using Finite Element Method and Its Stability Condition[J].Chinese Journal of Geophysics，1997，40(6)：833-841.(in Chinese)

[19] 张美根，王妙月，李小凡，等.各向异性弹性波场的有限元数值模拟[J].地球物理学进展，2002，17(3)：384-389，413.ZHANG Mei-gen，WANG Miao-yue，LI Xiao-fan，et a1.Finite Element Forward Modeling of Anisotropic Elastic Waves[J].Progress in Geophysics，2002，17(3)：384-389，413.(in Chinese)

[20] Komatitsc D，Tromp J.Introduction to the Spectral-element Method for 3-D Seismic Wave Propagation[J].Geophysical Journal International，1999，139(3)：806-822.

[21] Komatitsch D，Ritsema J，Tromp J.The Spectral-element Method，Beowulf Computing and Global Seismology [J].Science，2002，298(29)：1737-1742.

[22] Sinclair C，Greenhalgh S，Zhou B.2.5-D Modeling of Elastic Waves in Anisotropic Media Using the Spectral-element Method——a Preliminary Investigation[C]//AESC，Melbourne，Australia，2006.

[23] 李幼铭，胡健行.细胞自动机在地震波传播中的应用[J].地球物理学报，1995，38(5)：651-661.LI You-ming，HU JIAN-XING.Application of Cellijlar Automata Approach to the Study of Seismic Wave Propagation[J].Chinese Journal Geophysics，1995，38(5)：651-661.(in Chinese)

[24] 罗明秋，刘洪，李幼铭.地震波传播的哈密顿表述及辛几何算法[J].地球物理学报，2001，44(1)：120-128.LUO Ming-qiu，LIU Hong，LI You-ming.Hamiltonian Description and Symplectic Method of Seismic Wave Propagation[J].Chinese Journal Geophysics，2001，44(1)：120-128.(in Chinese).

[25] Zhou B，Greenhalgh S A，Zhe J.Numerical Seismogram Computations for Inhomegeneous Media Using a Short，Variable Length Convolutional Differentiator[J].Geophysical Prospecting，1993，41(6)：751-766..

[26] 张中杰，滕吉文，杨顶辉，等,声波与弹性波场数值模拟中的褶积微分算子法[J].地震学报，1996，18(1)：60-69.ZHANG Zhong-jie，TENG Ji-wen，YANG Dang-hui.Acoustic and Elastic Wave Modeling by a Convolutional Differentiator[J].Acta Seismoligica Sinica，1996，18(1)：63-69.(in Chinese)

[27] 戴志阳，孙建国，查显杰.地震波场模拟中的褶积微分算子法[J].吉林大学学报：地球科学版，2005，35(4)： 520-524.DAI Zhi-yang，SUN Jian-guo，CHA Xian-jie.Seismic Wave Field Modeling with Convolutional Differentiator Algorithm[J].Journal of Jilin Unviersity：Earth Science Edition，2005，35(4)：520-524.(in Chinese)

[28] 戴志阳，孙建国，查显杰.地震波混合阶褶积算法模拟[J].物探化探计算技术，2005，27(2)：111-114.DAI Zhi-yang，SUN Jian-guo，CHA Xian-jie.Seismic wave Field Modeling With Hybrid Order Convolution Differentiator[J].Computing Techniques for Geophysics and Geochem ical Exploration，2005，27(2)：111-114.(in Chinese)

[29] 谢靖.物探数据处理的数学方法[M].北京：地质出版社，1981，97-104.XIE Jing.Numerical Method for Geophysical Data-process[M].Beijing：Geological Publishing House，1981：97-104.

[30] Thomsen L.Weak-elastic Anisotropy[M].Geophysics，1986，51(10)：1954-1996.

[31] Zhou B，Greenhalgh S A.Seismic Scalar wave Equation Modeling by a Convolutional Differentiator[J].Bulletin of the Seismological Society of America，1992，82(1)：289-303.

[32] 贺同江，刘红艳，李小凡.粘弹性介质地震波传播的褶积微分算子法数值模拟研究[J].西北地震学报，2010，32(4)318-324.HE Tong-jiang，LIU Hong-yan，LI Xiao-fan.Numerical Simulation of Seismic Wave Propagation in Viscous-elastic Media by the Convolutional Differentiator Method[J].Northwestern Seismological Journal，2010，32(4)：318-324.(in Chinese)

[33] 刘少林，李小凡，汪文帅，等.最优化辛格式广义离散奇异核褶积微分算子地震波场模拟[J].地球物理学报，2013，56(7)：2452-2462.LIU Shao-lin，LI Xiao-fan，WANG Wen-shuai，et al.Optimal Symplectic Scheme and Generalized Discrete Convolutional Differentiator for Seismic Wave Modeling[J].Chinese Journal Geophysics，2013，56(7)：2452-2462.(in Chinese)

[34] 李一琼，李小凡，朱童.基于辛格式奇异核褶积微分算子的地震标量波场模拟[J].地球物理学报，2011，54(7)：1827-1834.LI Yi-qiong，LI Xiao-fan，ZHU Tong.The Seismic Scalar wave Field Modeling by Symplectic Scheme and Singular Kernel Convolutional Differentiator[J].Chinese Journal Geophysics，2011，54(7)：1827-1834.(in Chinese)

[35] 龙桂华，李小凡，张美根.基于Shannon奇异核理论的褶积微分算子在地震波场模拟中的应用[J].地球物理学报，2009，52(4)：1014-1024.LONG Gui-hua，LI Xiao-fan，ZHANG Mei-gen.The Application of Convolutional Differentiator in Seismic Modeling Based on Shannon Singular Kernel Theory[J].Chinese Journal Geophysics，2009，52(4)：1014-1024.(in Chinese)

Comparison of Convolutional Differentiator Method and the Pseudo-spectral Method Used in Seismic Wave Simulation

LIU Hong-yan， CHEN Yu-kun， LIU Fang

(EarthquakeAdministrationofTianjinMunicipality,Tianjin300201，China)

This study employed the convolutional differentiator seismic simulation method constructed using Forsyte polynomials to simulate the seismic wave propagation in a complex heterogeneous media and compared the simulation results to those obtained from the pseudo-spectral method.By comparing the computational time of the pseudo-spectral and convolutional differentiator methods，it was found that the computing time of the convolutional differentiator method is always less than that of pseudo-spectral method when the same computing environment is employed.Consequently，the convolutional differentiator seismic simulation method has an advantage in this regard.By comparing the snapshots obtained using the above mentioned methods，it was found that the dispersion of the snapshots by the convolutional differentiator is similar to that obtained by the pseudo-spectral method，which constituted another advantage for the convolutional differentiator method for seismic wave numerical simulation.Overall，the convolutional differentiator method is a viable alternative for seismic wave propagation simulation.

convolutional differentiator； pseudo-spectral method； seismic wave propagation； numerical simulation

2015-04-01

金基项目：地震科技星火计划项目(XH13002)．

刘红艳(1983-)，女，工程师，主要从事地震活动性、工程地震、复杂介质的地震波场数值模拟研究等工作.E-mail：liuhongyan_012@163.com

P315； P631.4

A

1000-0844(2015)02-0594-07

10.3969/j.issn.1000-0844.2015.02.0594