褶积微分算子法数值模拟的边界问题研究

贺同江,刘红艳,李小凡

(1.天津市地震局,天津 300201;2.中国科学院 地质与地球物理研究所,北京 100029)

褶积微分算子法数值模拟的边界问题研究

贺同江1,刘红艳1,李小凡2

(1.天津市地震局,天津 300201;2.中国科学院 地质与地球物理研究所,北京 100029)

褶积微分算子法是一种全新的数值模拟方法,已被广泛应用于复杂介质的地震波场数值模拟,但是其边界反射问题一直没有解决。这里将最佳匹配层(PML)吸收边界条件引入到褶积微分算子法中,此方法是在研究区域的边界上加入吸收层,使边界上传入吸收层的波,随传播距离按指数规律衰减,不产生任何反射,以达到消除边界反射的目的。构造不同的模型,通过对比分析证明,PML吸收边界条件能比较好地解决褶积微分算子法的边界问题,从而验证了完全匹配层吸收效果的优越性。

褶积微分算子;数值模拟;PML吸收边界

边界吸收的方法很多,可谓举不胜举,但是真正能有效的用于作者在本文所构造的褶积微分算子法数值模拟的却少之又少。通过模型试算分析可知,C layton吸收边界条件适用于有限差分数值模拟,但是对作者所构造的广义正交多项式褶积微分算子法的数值模拟的边界反射无效。Cerjan吸收边界条件适用于伪谱法数值模拟,理论上也能适用于作者在本文所构造的基于Fo rsyte广义正交多项式褶积微分算子法的数值模拟,但是采用Cerjan吸收边界条件并不能完全适用本文方法。因此作者在本文中采用PML吸收边界条件,通过模拟实例,证明了PML吸收边界能比较好地适用于褶积微分算子法。

1 Fo rsyte广义正交多项式褶积微分算子求导原理

作者在本文中利用错格有限差分褶积微分算子法来计算弹性介质中的地震波场,其主导思想是:利用基于Forsyte广义正交多项式褶积微分算子,有效表示计算波场对空间的偏导数,采用错格有限差分法计算对时间的偏导数,其计算思路类似于伪谱法[1]。

首先给出Fo rsyte多项式微分算子。Forsyte多项式是一个广义正交多项式,其插值函数为:

其中 P0=1

f(xi)为被插值函数f(x)在点xi处的值。在式(1)中的P0(x)、…、Pj+1(x)定义为Forsyte多项式系数。

对式(1)中的x求导,可得:

其中

Forsyte多项式微分算子可写为式(3)。

将上述微分算子式(3)离散化,可得式(4)。

其中 i为采样指标;Δx为沿着x轴的采样间隔。

就实际应用而言,须将微分算子截成短算子,这势必引起Gibbs现象。另外,多项式的引入还将引起Runge现象,为了消除这些现象,必须采用窗函数,以截断长微分算子。作者在本文采用的是下列Gaussian窗函数:

其中 m x为单边截断长度的采样数;c为常数;a(0.1≤a≤0.75)为衰减因子。

将微分算子式(4)用式(5)截断并锯齿化后,可得如下实用的一阶褶积微分算子:

通过对算子长度的调节及算子系数的优化,可同时兼顾波场解的全局信息与局部信息。作者运用算子长度为9点的一阶褶积微分算子,求解波动方程。而9点褶积算子的最优权系数为:

可以看出,算子的系数是反对称的。

2 二维各向同性介质一阶速度(应力弹性波动方程及PML吸收边界)

利用完全匹配层作为吸收边界的基本做法,是在所研究区域的四周引入完全匹配层。如图1所示,区域ABCD为要研究的区域,即我们要在此区域中研究波的传播问题。在区域的周围加上完全匹配层,在区域“1”中,令d(x)≠0;d(z)≠0,速度V都等于角点的速度。在区域“2”中,令d(x)=0;d(z)≠0,速度V在z方向为常数,在x方向和边界的速度相等。在区域“3”中,令d(x)≠0;d(z)=0,速度V在x方向为常数,在z方向和边界的速度相等。这样在计算边界的周围都有完全匹配层吸收介质,波由区域内通过边界传播到完全匹配层时,就不会产生任何反射。波在完全匹配层中传播时,也不会产生反射,并且按传播距离的指数规律衰减。当波传播到完全匹配层的边界时,波场近似为零,也不会产生反射。

图1 PML吸收边界示意图Fig.1 Schem atic diagram of the PML abso rbing boundary

对于二维非均匀各向同性介质,一阶速度,即应力弹性波动方程(假定体力为零)为式(7):

式中 vx=∂ux/∂t,vz=∂uz/∂t为速度的x、z分量;σxx、σzz为正应力;σxz为剪应力;λ、μ为拉梅常数;λ=ρ(V2p-2V2S);μ=ρVS2;VP、VS分别为介质纵波与横波速度;ρ为介质密度。

将一阶应力,即速度波动方程的波场分为与x方向和与z方向有关的二个子分量:

在式(8)中,上标x和z代表该项只与相应的空间导数有关。

根据二维非均匀各向同性介质波动方程组,可以分裂得到二维非均匀各向同性介质波动方程带有衰减因子的PML吸收边界系统方程组:

其中

式中 dPML是PML吸收层的厚度,作者取二十个网格点数;x、z是PML内计算点到内部区域和PML层交界面的距离;R为理论反射系数(作者在数值模拟中取为10-6),式(10)的解是衰减的。

d(x)和d(z)分别为x方向和z方向的衰减系数,也就是说d(x)起到衰减x方向传播的波,d(z)起到衰减z方向传播的波的作用。对于任意方向传播的波,可以通过矢量分解,分解成x方向和z方向传播的波,分别进行衰减。波场是按传播距离的指数规律衰减,衰减速度很快。当衰减系数d(x)、d(z)随空间位置变化时,不会在介质中产生任何反射。

3 二维一阶速度(应力弹性波动方程PML吸收边界条件下的数值解格式)

对模型区间离散后,设n、m、k分别是沿着空间x轴、z轴、时间t轴的采样点数,Δx、Δz、Δt分别是沿x轴、z轴、t轴的采样间隔,m x、m z是沿x轴、z轴采样数的半算子长度,则式(10)二维非均匀各向同性介质完全匹配层弹性波一阶速度——应力方程,其离散化的时间错格差分褶积微分算子法格式为(体力为零)式(11)。

4 数值算例

模型网格点数为256×256,震源位于模型正中,震源为45°倾斜集中力源,震源R icker子波主频为20 Hz。采用不同的吸收边界条件进行对比,模型参数如表1所示。

(1)图2(见下页)不加边界条件时,边界反射很强,并且在边界条件处,由于边界反射引起很严重的数值频散现象,使模拟精度降低。

(2)图3(见下页)是加入Cerjon吸收边界条件的模拟结果,相对图2边界反射引起的数值频散已几乎消失,但是边界反射还很严重。从图3可以看出,波遇到边界产生的反射波。

(3)图4(见下页)是加入PML吸收边界条件的模拟结果,看不到任何边界反射。从模拟结果看,边界条件是数值模拟中需解决的关键问题,作者在本文中所构造的方法,Cerjon吸收边界条件不能完全适用,但是PML吸收边界条件却能取得较好的结果,适于本文中所构造的方法。

我们利用一个比较复杂的非均匀介质模型,对本文方法加PML边界和Cerjon吸收边界作比较。复杂介质模型结构如后面图5所示。模型参数如后面表2所示。模型的网格大小为256×256,网格间距为d x=d z=10m,时间步长为1m s,震源放置于介质①中,坐标为(128,50),震源为45°倾斜集中力源,震源R icker子波主频为20 Hz。

表1 均匀各向同性介质模型参数Tab.1 M odelparam etersof homogeneous and isotrop icm edia

图2 均匀各向同性介质模型,未加入吸收边界条件400m s时刻波场快照Fig.2 W ave-field snapshots in homogeneous and isotrop icmodelw ith non-absorbing boundary(t=400m s,left:x component,right:z component)

图3 均匀各向同性介质模型,加入Cerjon吸收边界条件400m s时刻波场快照Fig.3 W ave-field snapsho ts in homogeneous and iso trop icmodelw ith Cerjon abso rbing boundary(t=400m s,left:x component,right:z component)

图4 均匀各向同性介质模型,加入PML吸收边界条件400m s时刻波场快照Fig.4 W ave-field snap shots in homogeneous and isotrop icmodelw ith PML absorbing boundary(t=400m s,left:x component,right:z component)

表2 复杂非均匀介质模型参数Tab.2 M odelparam etersof comp lex and inhomogeneousm edia

图5 复杂非均匀介质模型Fig.5 Com p lex and inhomogeneousmodel

从图6、图7(见下页)可以看出,作者在本文所构造的基于Forsyte广义正交多项式褶积微分算子法,能较好地模拟特别复杂的非均匀介质中的地震波场。采用Cerjon吸收边界条件,不能有效地消除边界反射,边界反射比较明显,影响了模拟的精度,因而不能完全适用于作者在本文所构造的,基于Forsyte广义正交多项式褶积微分算子法的波场数值模拟中。而采用PML吸收边界,不会产生任何反射,可达到消除边界反射的目的。因此采用PML吸收边界条件,能够完全适应本文方法。

5 结论

(1)采用完全匹配层(PML)吸收边界条件,是消除边界反射的理想方法之一。在研究区域的四周加入很薄的匹配层,以很小的计算代价,可换取很好的吸收效果。

(2)在研究区域的边界上加入吸收层,使边界上传入吸收层的波,随传播距离按指数规律衰减,不产生任何反射,可达到消除边界反射的目的。

图6 加入Cerjon吸收边界条件t=500m s时波场快照Fig.6 W ave-field snapshotsw ith Cerjon absorbing boundary of comp lexmodel(t=500m s,left:x component,right:z component)

(3)基于Forsyte广义正交多项式褶积微分算子法,能较好地模拟特别复杂的非均匀介质中的地震波场,PML吸收边界解决了本文方法的边界问题。作者在本文为褶积微分算子法的后续研究做了铺垫工作。可以预期,褶积微分算子法的推出及后续研究的成功开展,将为高精度地震波模拟,地震波偏移,地震反演,地震波成像,以及地震波在复杂非均匀介质中的传播等研究,提供更为广泛的选择。

[1]程冰洁,李小凡,龙桂华.基于广义正交多项式褶积微分算子的地震波场数值模拟方法[J].地球物理学报,2008,51(2):531.

[2]程冰洁,李小凡.2.5维地震波场褶积微分算子法数值模拟[J].地球物理学进展,2008,23(4):1099.

[3]刘红艳,李小凡.非均匀介质中地震波传播的数值模拟[J].物探化探计算技术,2008,30(3):173.

[4]李信富,李小凡.复杂介质地震波传播的褶积微分算子数值模拟[J].地震学报,2008,30(4):377.

[5]CLAYTON R,ENGQU ISTB.Absorbing boundary conditions for acoustic and elastic wave equations[J].Bu ll.Seis.Soc.Am.,1997,67:1529.

[6]CERJAN C,KOSLOFFD.A non-reflection boundary condition for discrete acoustic and elasticw ave equation[J].Geophysics,1985,50:705.

[7]BERENGER JP.Three-dim ensionalperfectlym atched layer for the absorp tion electrom agnetic waves[J].Comp.Phys,1996,127:363.

[8]CHW EW C,L IU Q H.Perfectly m atched layers for elastodynam ics:A new absorbing boundary condition[J].J.Comp.A coust.,1996,4(4):72.

[9]FRANC ISCOLL INO,CHRYSOULA TSOGKA.App lication of the perfectlym atched absorbing lamodel to the linear elastodynam ic p rob lem in anisotrop ic heterogeneousm edia[J].Geophysics,2001,66(1):294.

[10]FESTA G,N IELSEN S.PML Absorbing Boundaries[J].Bulletin of the Seismo logical Society Am erica,2003,93(2):891.

[11]BECACHE E,FAUQUEUX S,JOLY P.Stability of perfectly m atched layers,group velocities and anisotrop ic waves[J].Journal of Com pu tational Physics,2003,188:399.

[12]ZENG Y Q,HE JQ,L IU Q H.The app lication of the perfectlym atched layer in num ericalmodeling ofwave p ropagation in po roelastic m edia[J].Geophysics,2001,66(4):1258.

[13]SONGRUOLONG,MA JUN,WANG KEX IE.The App lication of the nonsp litting perfectlym atched layer in num ericalmodeling of wave p ropagation in poroelastic m edia[J].App lied Geophysics,2005,2(4):216.

[14]SHUO MA,PENGCHENG L IU.M odeling of the PerfectlyM atched Layer absorbing boundary and intrinsic attenuation in exp licit Finite-Elem entm ethods[J].Bu lletin of the Seismo logical Society of Am erica,2006,96(5):1779.

[15]赵海波,王秀明,王东,等.完全匹配层吸收边界在孔隙介质弹性波模拟中的应用[J].地球物理学报,2007,50(2):581.

图7 加入PML吸收边界t=500m s时波场快照Fig.7 W ave-field snapshotsw ith PML abso rbing boundary of comp lexmodel(t=500m s,left:x component,right:z component)

P 631.4

A

1001—1749(2010)06—0571—07

0 前言

国家自然科学基金重点项目资助(40437018)

2009-07-28

贺同江(1979-),男,硕士,现在天津市地震局监测预报中心工作。

褶积微分算子法作为一种全新的数值模拟方法,已被广泛应用于复杂介质的地震波场数值模拟中,该方法同时具有广义正交多项式方法的高精度和短算子有限差分算法的高速度。通过对算子长度的调节及算子系数的优化,可同时兼顾波场解的全局信息与局部信息。该算法的计算速度快,计算效率高,能够直观、高效地反映介质中波场的传播规律[1～4]。但是在进行波动方程数值模拟时,考虑到计算机的有限内存及有限计算时间,要对考虑问题的无限区域进行截取,使进行中的数值模拟在有限的区域内完成,即可以在有限的区域中得到较高精度的无限空间解。为此需要引入人工边界来达到此目的,但这样会在人工边界处产生人为反射,如不消除或者压制这种虚假反射,就会影响数值模拟的结果和精度。为此,边界条件问题是必须要解决的,也是较难解决的问题。地震波数值模拟中,在人工边界上使用吸收边界条件,应尽可能地降低边界上所产生的反射波强度,使得反射波在边界上好像被“吸收”了一样,从而大大减小对计算区域内精度的影响,以达到边界吸收的目的。应用较普遍的是C layton、Engquist等人[5]提出的吸收边界条件(即CE边界),该条件在旁轴近似理论的基础上导出,在特定的入射角和频率范围内,具有较好的吸收效果。Cerjan etal.[6]引入了有损吸收层来衰减外行波,但这种方法需要足够厚的吸收层,才可以得到满意的效果,因此计算量较大。Berenger[7]对海绵吸收边界条件作了进一步的改善,提出一种新的人工边界条件,即最佳匹配层法(PML吸收边界),这种方法几乎达到零反射。目前,PML技术已广泛应用于声波和弹性波的有关问题。Chew和L iu.[8]首先证明了PML技术可应用于弹性波的数值模拟中;Co llino和Tsogka[9]把PML用在各向异性非均匀的介质情况;Festa和N ielson[10]针对各向同性弹性介质PML情况下,离散波动方程稳定性进行了详细讨论,并推导了PML理论反射系数计算公式;Becache et al.[11]在理论上讨论了PML在各向异性情况下稳定性问题。另外PML也应用到其它问题上,如黏弹性介质、双相孔隙弹性介质中[12～15]。