地球物理学进展  2015, Vol. 30 Issue (2): 565-570   PDF    
基于辛算法及FCT的地震波场模拟
黄明曦, 薛霆虓 , 王有学    
桂林理工大学地球科学学院, 桂林 541004
摘要:为了解决在长时间和复杂结构的地震波场数值模拟情况时的数值频散问题, 现在辛算法的基础上, 主要结合通量校正传输(FCT)、褶积微分算子、完全匹配层(PML)等数值模拟技术, 寻找一种更为优秀的地震波场数值模拟方法.地震波场的数值模拟结果表明, 辛算法不仅具有保持体系结构的特性, 并且具有长时间跟踪能力, 具有很强的数值模拟稳定性;FCT方法基于通量守恒原理, 压制网格频散效果明显;褶积微分算子突出了空间微分的局部属性.因此, 通过合理应用各种技术, 可以对地震波场特征进行更精确的数值模拟.
关键词辛算法     通量校正传输     褶积微分算子     完全匹配层     数值频散    
Seismic wave field simulation on the symplectic and FCT algorithm
HUANG Ming-xi, XUE Ting-xiao , WANG You-xue    
College of Earth Sciences, Guilin University of Technology, Guilin 541004, China
Abstract: To find out a solution of numerical dispersion in long-time tracing and complicated structure modeling, based upon the symplectical algorithm, this paper introduced a 2D seismic wavefield modeling technique which combins with the Flux-corrected Transport method (FCT) and convolutional differentiator, and uses the Perfectly Matched Layer (PML) as a absorbing boundaries. The numerical modeling shows that the Symplectical algorithm not only has a structure-preserving feature but also has a long-time tracing ability, it has a good stability during numerical modeling. Due to adopting convolutional differentiator and FCT techniques, the spatial feature is enhanced, and the numerical dispersion caused by the gridding is suppressed to produce a more accurate seismic wavefield. So we can get a better result with combins these technique reasonable.
Key words: Symplectical algorithm     FCT     convolutional differentiator     PML     numerical dispersion    
0 引 言

对于复杂介质中的地震波场来说,直接求其解析解是几乎不可能的事情,因此一般是通过地震波场的数值模拟来实现.在地震波场数值模拟中,弹性波动方程模拟较之运动学的射线追踪方法具有更丰富的动力学信息,因而,地震波场数值模拟是地震学研究中的重要手段之一.在弹性波动方程的数值模拟方面前人已经提出了很多种方法,1968年Alterman和Karal(1986)提出了弹性波动方程的有限差分法;1970年Aki和Larner(1970)提出了谱分析方法;1975年Smith(1975)提出了有限元法等.

数值模拟中所使用的算法应当尽可能的保持系统本身的性质,辛算法正是如此理念的算法.1984年,冯康和秦孟兆(1984)先生首次系统的提出了哈密尔顿系统的辛几何算法,“哈密尔顿体系是动力系统的重要体系,一切真实的、耗散可忽略不计的物理过程都可表示成哈密尔顿体系”.哈密尔顿系统的辛几何算法较之传统方法有保持体系结构稳定的优点,使其在计算时具有长时间的跟踪能力.经过对辛算法的基本理论的建立(Yoshida,1990;Suzuki,1992),其在在科学计算领域有着广泛的应用,在地震波场在高频近似条件下(李世雄,2001),辛算法可以用于研究地震波场的射线走向和波阵面形态,这对地震波场传播的渐近解研究、模型研究及波场成像的运用具有广泛的实用价值.因此本文使用辛算法对波动方程的时间微分项进行处理,而在空间微分项上将使用褶积微分算子.褶积微分算子是Zhou和Greenhalgh(1992)推导出的用于求解微分方程的方法,算子通过取不同的窗系数来灵活控制衰减速度,具有计算精度高,稳定性好的特点,可作为有限差分算子的一种替代.本文使用的是Shannon奇异核褶积微分算子(龙桂华等,2009),该算子被证明优于一般的高阶有限差分算子,且在与辛算法的结合使用时取得了较好的效果(李小凡等,2011李小凡等,2012).在地震波场的数字模拟过程中存在频散问题和边界吸收处理问题,这些问题的存在会严重制约地震波场的数字模拟计算速度与效果.通量校正传输(FCT-Flux-corrected Transport)方法,是由Boris和Book(1973)等在对流体动力学连续方程的求解过程中提出的方法,而后,有效地压制差分计算产生的数值频散.完全匹配层(PML-Perfectly Matched Layer)方法在处理边界吸收方面是一种广为所使用的有效方法,在理论上能够在边界上完全吸收任意频率和任意角度入射的弹性波,不产生任何反射.在本文中将使用辛算法以及其他数值计算技术对二维空间中的地震波场进行数值模拟和分析. 1 地震波场数值模拟方法

在二维各向同性介质中,假定外力为零的情况下,用速度-应力表示的弹性波方程为

其中vx、vz分别表示弹性介质中质点的在水平和垂直两个方向上的振动速度,τxx、τzz和τxz分别为主应力和切应力,λ和μ为拉梅系数,ρ为介质密度.

1.1 辛差分格式

在波动方程(1)中,令 v =(vx,vz),τ =(τxx,τzz,τxz),则(v,τ)=(vx,vz,τxx,τzz,τxz),根据辛算法理论,k级n阶显式辛格式可以写为

其中Δt表示时间步长,矩阵I为与 P、Q 同阶的单位矩阵.

根据(2)式可以得到波动方程(1)式的单步三级三阶显式辛格式解为

其中的参数ai、biIwatsu(2009)给出

1.2 通量校正传输

通量校正传输(FCT)是应用于流体动力学连续方程求解过程中的一种方法,该方法是基于守恒型方程的差分格式,要求交界面两侧流入和流出的通量守恒,这符合算法系统中能量守恒的原则.该方法的原理是计算每个网格边界地震波场的扩散,对数值频散进行压制,其步骤称之为通量校正.

在通量校正的第一阶段,首先计算第n个时间点的波场扩散通量为

其中η1可以为一个常数或者一个函数,其值越大代表平滑能力越强,如果过大则真实波场将被抑制的模糊不清,过小则失去压制频散的效果,通常其取值范围为0.01≤η1≤0.1.然后,利用下式计算n+1时间点的扩散通量为

其中η2的取值与η1类似,但一般应比η1大10%至15%,这是因为要补偿差分运算带来的损失.η1与η2的数值大小的应按实际情况酌情选择.最后再对第n+1个时间点进行扩散校正得

在通量校正的第二阶段,首先计算反扩散通量为

然后,对第n+1个时间点进行反扩散校正得

1.3 褶积微分算子

根据连续信号的离散原理,一个连续信号u(x)的离散信号可以表示为连续信号u(x)与Delta函数的褶积,而Shannon奇异核函数可以很好地近似Delta函数,将其归一化及离散后,Shannon奇异核函数、离散形式及其一阶微分算子分别为

其中Δ是x步长,σ为高斯函数包络线宽度.在使用微分算子(12)计算时,通常利用Hanning窗对其进行截断处理,以便有效抑制Gibbs效应,并且能够让能量更加集中在主瓣附近.Hanning窗表达式为

其中α和β为度量Hanning窗函数的常数,这里选择为α=0.54,β=8.最终进过截断处理后的微分算子可表示为

其中W为信号长度之半.

1.4 完全匹配层

在处理吸收边界时,本文采用完全匹配层(PML)吸收边界.PML吸收边界的实质是在匹配层内设置一种各向异性耗散介质,仅在垂直于边界方向上引入阻尼因子,通过选取适当的参数,可使自由介质与人工边界的分界面上仅产生很微弱的反射.其方法实现为对(1)式的时间域变量进行分裂,得到带有衰减因子的吸收边界方程组

其中dx、dz是与x、z方向导数有关的衰减因子

这里L为PML边界总厚度,lx、lz为由PML边界内到自由介质分界面的距离,d0为吸收初值,R=10-6为理论反射系数,β为控制PML边界的吸收效果的参数,可以根据情况给定,通常为匹配层内的最大纵波速度vpml. 2 数值模拟

2.1 均匀各项同性介质

对于均匀各项同性的弹性介质模型,假设其纵波速度2000 m/s,横波速度1154.7 m/s,密度2073.1 kg/m3,模型长度和深度均为600 m,模型的网格大小为2 m×2 m.震源使用80 Hz的Ricker子波,设在模型的正中央(300 m,300 m)处,时间步长为0.2 ms.

利用高阶交错网格算法计算得到的80 ms(400个时间步长)时垂直分量的波场快照如图 1a所示;而结合以上几种算法(以下简称辛算法)得到的80 ms时垂直分量的波场快照如图 1b所示.从波场快照来看,高阶交错网格算法(图 1a)和辛算法(图 1b)都能准确的模拟出地震波场,说明两种方法在短时间内都是可靠的数值模拟的手段.当地震波场经过1 s(5000个时间步长)时,利用高阶交错网格算法和辛算法计算得到的波场快照如图 2所示.此时,两种算法的差异明显的体现了出来:由高阶交错网格算法得到的波场(图 2a)数字频散非常严重,完全无法分辨出波前;然而,根据辛算法得到的波场(图 2b)却能清晰地分辨出波前,体现出辛算法长时间跟踪能力的优越性.

图 1 80 ms时的垂直分量波场快照
(a)高阶交错网格算法;(b)辛算法
Fig. 1 The vertical component of seismic wavefield at 80 ms
(a)Staggered-grid;(b)Symplectical.

图 2 1 s时的垂直分量波场快照
(a)高阶交错网格算法;(b)辛算法
Fig. 2 The vertical component of seismic wavefield at 1 s
(a)Staggered-grid;(b)Symplectical.
2.2 Marmousi模型

为了进一步验证辛算法对复杂介质的适用性,将使用的Marmousi模型(图 3)做正演模拟.本文使用的Marmousi模型长度为3072 m,深度为976 m,网格大小为4 m×4 m.震源使用20 Hz的Ricker子波,设在模型的左上角(0 m,0 m),时间步长为0.3 ms.在模型网格数较少,地层的速度变化快的情况下,其对算法的考验更高.通过对比地表(z=0 m)地震记录可以看出高阶交错网格算法(图 4a)的数值频散足以造成虚假地层信息,而辛算法(图 4b)仍然能模拟出这些细小的地层速度变化.

图 3 Marmousi模型Fig. 3 Marmousi model

图 4 Marmousi模型的水平分量地震记录Fig. 4 The horizontal component seismograms of Marmousi
(a)Staggered-grid;(b)Symplectical.
2.3 各项同性复杂介质

对一般的复杂结构的地震波场数值模拟,我们可以观察理论上的模型应有的震相信息是否准确.如图 5含有空洞的两层介质速度模型,模型长度为800 m,深度为600 m,介质分界面的埋深为160 m,空洞顶面的埋深为400 m,空洞的左边界位于140 m处,宽为160 m,模型的网格大小为2 m×2 m.在进行波场模拟时,在模型的四周各增加90 m(45个网格宽度)的PML吸收边界,震源仍使用80 Hz的Ricker子波,设在模型的原点(0,0)处,时间步长为0.2 ms.由数值模拟得到的地表处的地震记录(图 6)可以看出,各种震相都有清晰楚的波至,数值频散得到了很好的压制.

图 5 各项同性复杂介质Fig. 5 The complicated isotropic structure

图 6 各项同性复杂介质中的地震记录
(a)垂直分量;(b)水平分量;(c)模型
Fig. 6 The complicated isotropic structure and seismograms
(a)vertical;(b)horizontal;(c)model.

为了方便接下来的分析,本文规定地震波在第一层传播的纵、横波编码分别为P、S,在第二层传播的纵、横波编码分别为p、s,左侧的绕射波用下标dl表示,左侧的绕射波用下标dr表示,直达波为g,折射波为rf.

震源首先激发了直达波Pg和介质分界面产生的纵波的折射波PPrf和转换横波的折射波PSrf,在介质分界面不仅产生了反射波PP、PS一次反射波震相,而且在介质分界面与空洞顶面还形成多次反射纵波PppppP;同时在空洞顶部两个角点产生了绕射波纵波PppPdl、PppPdr、PppSdl、PpsPdl、PpsPdr、PspPdr、PpsSdl、PpsSdl、PpsSdr、PpsSdl、PssSdl、PssSdr.由于震源为胀缩源,地震波在各项同性均匀介质中球面波前法向方向上的振幅随着远离震源而迅速衰减,在地表接收到的垂直分量振幅急剧减小(图 6a),Pg震相尤为明显,而水平分量在地表恰好与球面波波前的法向方向一致,因而其能力较大(图 6b),但是仍然呈现为逐渐衰减的趋势.在介质分界面处的一次反射波PP在垂直分量的能量较强(图 6a),而在水平分量的能量较弱(图 6b),反射转换波PS则正好相反;由于震源为胀缩源,随着波前远离震源,分界面处的入射波的入射角度逐渐增大,PP震相的垂直分量逐渐衰减,而PP震相的水平分量和PS震相的垂直分量与水平分量均呈现逐渐增大继而逐渐减小的趋势;对于折射纵波PPrf、折射转换横波PSrf,他们分别在300m、200m左右处与反射波PP、PS相切,而且能量变弱.

在空洞的顶面也产生了一系列的反射波,其中在垂直分量的地震波场中PppP的具有较强的能量,而且可以观测到清晰的多次反射波PppppP;而在水平分量上PppP、PppS、PpsS、PssS各震相具有大小相似的能量分布,但是多次反射波PppppP能量十分微弱以至于无法被观察到.同时,在空洞顶面的左、右两个角点还有一系列的绕射波产生,虽然能量较弱,但由于旅行时间相近从而波形相互叠加增强因此能够被观察到.如在水平分量的波场中,绕射波PspSdr与PpsSdr的旅行时间极小点清晰地指示了绕射点的水平位置. 3 结 论

本文通过对辛算法的研究,结合FCT等数值模拟技术,对各向同性弹性介质中的地震波场进行了数值模拟和结果分析,并且通过对比,可以得到如下结论:

(1)辛算法不仅具有保持体系结构的特性,并且具有长时间跟踪能力,虽然在计算时增加了运算量,但具有很强的数值模拟稳定性;

(2)采用了FCT及褶积微分算子数值模拟技术,突出了空间微分的局部属性,压制网格频散,可以对地震波场特征进行更精确的数值模拟;

(3)对多种技术进行适当整合,合理利用,可以使地震波场模拟结果更加准确,适用于数值模拟中长时间及模型结构复杂的需要.

致 谢 感谢我的导师对我的指导,感谢专家、编辑部老师对论文的中肯意见.
参考文献
[1] Aki K, Larner K L. 1970. Surface motion of a layered medium having an irregular interface due to incident plane SH waves[J]. J. Geophys. Res, 75(5): 933-954.
[2] Alterman Z, Karal F C Jr. 1968. Propagation of elastic waves in layered media by finite difference methods[J]. Bull. Seism. Soc. Am., 58(1): 367-398.
[3] Boris J P, Book D L. 1973. Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works[J]. J. Comp. Phys., 11(1): 38-69.
[4] Feng K, Qin M Z. 2003. Symplectlc Geometric Algorlthm for Hamiltonian Systems (in Chinese)[M]. Hangzhou: Zhejiang Science and Technology Press.
[5] Huang C, Dong L G. 2009. Staggered-grid high-order finite-difference method in elastic wave simulation with variable grids and local time-steps[J]. Chinese Journal of Geophysics (in Chinese), 52(11): 2870-2878, doi: 10.3969/j.issn.0001-5733.2009.11.022.
[6] Iwatsu R. 2009. Two new solutions to the third-order symplectic integration method[J]. Phys. Lett. A, 373(34): 3056-3060.
[7] Li S X. 2001.The High-Frequency Approximation of Wave Equations and Symplectic Algorithm (in Chinese)[M]. Beijing: Science Press.
[8] Li Y Q, Li X F, Zhu T. 2011. The seismic scalar wave field modeling by symplectic scheme and singular kernel convolutional differentiator[J]. Chinese J. Geophys. (in Chinese), 54(7): 1827-1834, doi: 10.3969/j.issn.0001-5733.2011.07.016.
[9] Li Y Q, Li X F, Zhang M G. 2012.The elastic wave fields modeling by symplectic discrete singular convolution differentiator method[J]. Chinese J. Geophys. (in Chinese), 55(5): 1725-1731. doi: 10.6038/j.issn.0001-5733.2012.05.029.
[10] Long G H, Li X F, Zhang M G. 2009. The application of convolutional differentiator in seismic modeling based on Shannon singular kernel theory[J]. Chinese J. Geophys. (in Chinese), 52(4): 1014-1024.
[11] Luo M Q, Liu H, Li Y M. 2001. Hamiltonian description and symplectic method of seismic wave propagation[J]. Chinese Journal of Geophysics (in Chinese), 44(1): 120-128.
[12] Qin M Z, Chen J B. 2000. Maslov asymptotic theory and symplectic algorithm[J]. Chinese Journal of Geophysics (in Chinese), 43(4): 522-533.
[13] Smith W D. 1975. The Application of finiteelement analysis to body wave propagation problems[J]. Geophysics J. Roy Astr. Soc., 42(2): 747-768.
[14] Sun L J, Yin X Y. 2011. A finite-difference scheme based on PML boundary condition with high power grid step variation[J]. Chinese Journal of Geophysics (in Chinese), 54(6): 1614-1623, doi: 10.3969/j.issn.0001-5733.2011.06.021.
[15] Suzuki M. 1992. General theory of higher-order decomposition of exponential operators and symplectic integrators[J]. Phys. Lett. A, 165(5-6): 387-395.
[16] Sun Y H, Zhou Y C, Li S G, et al. 2006. A windowed Fourier pseudospectral method for hyperbolic conservation laws[J]. J. Comp. Appl. Math., 214(2): 466-490.
[17] Wang W S, Li X F, Lu M W, et al. 2012. Structure-preserving modeling for seismic wavefields based upon a multisymplectic spectral element method[J]. Chinese Journal of Geophysics (in Chinese), 55(10): 3427-3439, doi: 10.6038/j.issn.0001-5733.2012.10.026.
[18] Wei G W. 1999. Discrete singular convolution for the solution of the Fokker-Planck equation[J]. J. Chem. Phys., 110(18): 8930-8942.
[19] Yang K D, Song G J, Li J S. 2011. FCT compact difference simulation of wave propagation based on the Biot and the squirt-flow coupling interaction[J]. Chinese J. Geophys. (in Chinese), 54(5): 1348-1357, doi: 10.3969/j.issn.0001-5733.2011.05.024.
[20] Yoshida H. 1990. Construction of higher order symplectic integrators[J]. Phys. Lett.A, 150(5-7): 262-269.
[21] Zhou B, Greenhalgh S A. 1992. Seismic scalar wave equation modeling by a convolutional differentiator[J]. Bull. Seism. Soc. Am., 82(1): 289-303.
[22] 冯康, 秦孟兆. 2003. 哈密尔顿系统的辛几何算法[M]. 杭州: 浙江科学技术出版社.
[23] 黄超, 董良国. 2009. 可变网格与局部时间步长的交错网格高阶差分弹性波模拟[J]. 地球物理学报, 52(11): 2870-2878, doi: 10.3969/j.issn.0001-5733.2009.11.022.
[24] 李世雄. 2001. 波动方程的高频近似与辛几何[M]. 北京: 科学出版社.
[25] 李一琼, 李小凡, 朱童. 2011. 基于辛格式奇异核褶积微分算子的地震标量波场模拟[J]. 地球物理学报, 54(7): 1827-1834, doi: 10.3969/j.issn.0001-5733.2011.07.016.
[26] 李一琼, 李小凡, 张美根. 2012. 基于辛格式离散奇异褶积微分算子的弹性波场模拟[J]. 地球物理学报, 55(5): 1725-1731, doi: 10.6038/j.issn.0001-5733.2012.05.029.
[27] 龙桂华, 李小凡, 张美根. 2009. 基于Shannon奇异核理论的褶积微分算子在地震波场模拟中的应用[J]. 地球物理学报, 52(4): 1014-1024, doi: 10.3969/j.issn.0001-5733.2009.04.018.
[28] 罗明秋, 刘洪, 李幼铭. 2001. 地震波传播的哈密顿表述及辛几何算法[J]. 地球物理学报, 44(1): 120-128.
[29] 秦孟兆, 陈景波. 2000. Maslov渐近理论与辛几何算法[J]. 地球物理学报, 43(4): 522-533.
[30] 孙林洁, 印兴耀. 2011. 基于PML边界条件的高倍可变网格有限差分数值模拟方法[J]. 地球物理学报, 54(6): 1614-1623, doi: 10.3969/j.issn.0001-5733.2011.06.021.
[31] 汪文帅, 李小凡, 鲁明文,等. 2012. 基于多辛结构谱元法的保结构地震波场模拟[J]. 地球物理学报, 55(10): 3427-3439, doi: 10.6038/j.issn.0001-5733.2012.10.026.
[32] 杨宽德,宋国杰,李静爽. 2011. Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟[J]. 地球物理学报, 54(5): 1348-1357, doi: 10.3969/j.issn.0001-5733.2011.05.024.