有限差分法地震波场数值模拟中的吸收边界条件

交错网格有限差分正演模拟的联合吸收边界胡建林;宋维琪;张建坤;邢文军;徐文会【摘要】三维声波方程交错网格有限差分正演模拟中的边界问题一直是热点问题.完全匹配层吸收边界(PML)具有较强且稳定的吸收效果,但必须具有一定的边界厚度才能吸收干净,这就增大了三维正演模拟的模型空间,即增加了运算量;Higdon边界能消除任意角度入射波的边界反射,也具有较强稳定性,但该高阶吸收边界离散化后过于复杂,而低阶时吸收效果不如PML边界.因此,基于对PML吸收层中的平面波传播规律的研究,重新推导PML最外层的Higdon吸收边界条件,得到含PML吸收系数的新的Higdon吸收边界条件.联合吸收边界不仅可使用较小厚度(相对于单纯PML边界)的PML层对分量进行衰减,而且在PML边界外层,能应用新推导的Higdon吸收边界条件对反射波进行匹配吸收.在相同吸收效果下,联合吸收边界大幅度降低了PML厚度,减小了运算量,得到精确的模拟结果.【期刊名称】《石油地球物理勘探》【年(卷),期】2018(053)005【总页数】7页(P914-920)【关键词】三维声波方程;交错网格有限差分;正演;PML边界;Higdon边界;联合吸收边界【作者】胡建林;宋维琪;张建坤;邢文军;徐文会【作者单位】中国石油大学(华东)地球科学与技术学院,山东青岛266555;中国石油大学(华东)地球科学与技术学院,山东青岛266555;中国石油冀东油田公司勘探开发研究院,河北唐山063004;中国石油冀东油田公司勘探开发研究院,河北唐山063004;中国石油冀东油田公司勘探开发研究院,河北唐山063004【正文语种】中文【中图分类】P6311 引言复杂地下介质中，地震波的传播过程繁冗，难以得到解析解，因此，一般是通过正演模拟探究地下地震波的传播。

在地震波正演模拟中，利用波动方程的正演模拟比用运动学的射线追踪法可获得更丰富的动力学信息，因此地震波场的数值模拟是地震波场传播研究中的重要手段之一[1-8]。

吸收边界条件的研究及其应用吸收边界条件的研究及其应用摘要：边界条件是数值计算领域中一个重要的概念，在传统的数值计算过程中，常常需要定义边界条件来模拟现实问题中的边界情况。

本文主要研究并探讨了吸收边界条件的概念、数值模拟方法以及在不同领域的应用。

通过对吸收边界条件的深入研究与应用分析，可以为数值计算领域中的边界问题提供解决思路与方法。

一、引言边界条件是数值计算中模拟实际问题时必不可少的一部分，因为在实际问题中物体通常都是有界限的，通过定义边界条件可以模拟出物体在边界处的行为。

而吸收边界条件则是在计算模拟过程中为了有效地模拟无穷远处场景时使用的一种边界条件。

二、吸收边界条件的概念1. 吸收边界条件的定义吸收边界条件是一种能够吸收远场无穷大距离处的波前能量的边界条件。

在数值计算模拟中，通常通过设定一定的参数来实现吸收边界条件，以模拟出远场无穷远处能量的衰减过程。

2. 吸收边界条件的分类吸收边界条件可以按照其应用领域和具体实现方式进行分类。

按照应用领域可分为电磁、声波、流体力学等各个领域的吸收边界条件；按照实现方式可分为频域方法、时间域方法等。

三、吸收边界条件的数值模拟方法1. 建立数值模型在吸收边界条件的数值模拟中，首先需要建立数值模型。

根据具体问题的不同，可以选择合适的数值计算方法，如有限差分法、有限元法等。

2. 设定吸收边界条件参数在数值模拟中，需要通过设定一定的吸收边界条件参数来实现吸收边界条件。

参数的设定需要考虑问题的具体情况，如波长、波速等。

3. 选取吸收边界条件算法根据具体模拟问题的特点，可以选择合适的吸收边界条件算法。

常用的方法有海绵层吸收、复杂是吸收边界条件法等。

四、吸收边界条件的应用1. 电磁学领域的应用在电磁学领域中，吸收边界条件可以用于模拟天线辐射、电磁波传播等问题。

通过合理地选择和设定吸收边界条件，可以精确地模拟电磁场在边界处的行为。

2. 声波学领域的应用在声波学领域中，吸收边界条件可以用于模拟声波传播、声源辐射等问题。

第22卷　第2期地　球　物　理　学　进　展Vol.22　No.22007年4月(页码:487～491)PRO GRESS IN GEOP H YSICSApr.　2007地震波有限差分模拟综述冯英杰,　杨长春,　吴　萍(中国科学院地质与地球物理研究所,北京100029)摘　要　本文从有限差分法数值模拟技术的各个方面对地震波有限差分模拟的发展和现状进行了论述.波场的数值模拟技术是认识地震波传播规律,检验各种处理方法正确性的重要工具,地震波的数值模拟是地震波传播规律研究的必要手段,贯穿于地震资料的采集、处理、解释的整个过程中.有限差分法数值模拟技术相对于射线方法具有更高的精度,同时比有限元方法计算量小,因此在实际应用中占很重要的地位.关键词　有限差分,差分格式,震源,边界条件,数值频散中图分类号　P631 文献标识码　A 文章编号　100422903(2007)022*******The review of the f inite 2difference elastic w ave motion modelingFEN G Y ing 2jie ,　YAN G Chang 2chun ,　WU Ping(I nstit ute of Geology and Geophysics ,Chinese A cadem y of sciences ,B ei j ing 100029,Chi na )Abstract The numerical seismic wave propagation modeling is a powerf ul tool in the oil exploration ,such as the date collection ,the processing and the interpretation and so on .It can not only find out the properties of the media ,but also check the validity of processing methods ,recognize the law of the wave propagation.In all the numerical meth 2ods ,the finite 2difference method is more usef ul with its advantages ,such as high precision ,flexibility ,costless.In this paper ,several parts of the finite 2difference method are discussed ,such as the finite -scheme ,the source prob 2lem ,the boundary condition and the numerical dispersion dumbness.K eyw ords finite 2difference ,source ,boundary condition ,wave propagation ,numerical dispersion收稿日期　2006210208;　修回日期　2006212220.基金项目　国家973项目(2005CB422104)和中国科学院知识创新工程重大项目资助(KZCX12SW 218204)联合资助.作者简介　冯英杰,女,1980年生,山东昌邑人,硕士,中国科学院地质与地球物理研究所,主要从事油储地球物理方面的研究.(E 2mail :fyj@ )0　引　言地震波场的数值模拟技术是在已知地下介质结构和参数的情况下,利用理论计算的方法研究地震波在地下介质中的传播规律,合成地震记录的一种技术.随着地震勘探技术的发展,数值模拟成为贯穿地震数据采集、处理和解释全过程的一种重要方法,在确定观测系统的合理性,检验处理和解释的正确性等方面有着越来越广泛的应用.地震勘探中的数值模拟方法主要以射线理论和波动方程理论为基础,有射线追踪法、柯希霍夫积分法、有限元法、有限差分法和伪谱法,还有将有限元和有限差分结合到一起的区域分裂法等.有限差分法是最常用的一种正演模拟方法,它将波动方程中波场函数的空间导数和时间导数用相应空间和时间的差分来代替.有限差分法虽然计算精度较有限元低一些,但是它的计算速度较有限元要快.1　有限差分模拟的历史有限差分法数值模拟技术开始于上世纪70年代初,Alterman 等人(1968)作了开创性的工作,使用显式有限差分格式获得了层状介质二阶弹性波方程的离散数值解.Alterman 等人实际上得到的是均匀介质弹性波数值解,只在内界面运用了应力和位移连续的内边界条件,使得波能通过弹性界面传播,对于结构复杂和不规则的岩性层面,必须使用适应非均匀介质模型的方法,即自动满足内界面处应力和位移连续的有限差分格式.Boore (1972)提出了非均匀介质二阶弹性波有限差分方法,Kelly 等(1976)改进发展了这一方法.Madariaga (1976)提地　球　物　理　学　进　展22卷出了非均匀介质速度-应力弹性波方程组交错网格有限差分方法,Virieux(1984,1986)利用这一格式完成了对弹性介质的P-SV和SH波的速度-应力方程组的正演计算,成为弹性波数值模拟的经典之作.Igel等人(1995)实现了各向异性介质交错网格有限差分波传播模拟,1996年他又在柱坐标和球坐标下实现了有限差分模拟.国内也有很多学者(王秀明,2003,王德利,2005)将这一格式运用到波场模拟中,揭示了波在地下传播的一些特性.为了适应地下介质多尺度非均匀性和不规则自由边界,避免局部采样过疏或过密的问题,后来又发展了一系列不规则网格的有限差分模拟(J ast ram,1992,1994; Falk,1996,1998;张剑锋,1998,2000;Tessmer, 2000;杨顶辉,1996).Carcione(2001)一直致力于粘弹性、各向异性、孔隙多相流体介质地震波传播的研究和数值模拟,他在2002年发表在Geop hysics上的文章是对数值模拟技术现状很好的总结.在有限差分正演中,通常有以下几方面的问题需要考虑:差分格式、震源函数、边界条件、数值稳定性和频散效应,以下将这几个方面来论述其发展现状.2　差分格式有限差分数值模拟与其他数值分析方法一样,必须把连续问题离散化.因此首先要对求解区域也就是弹性介质模型进行网格剖分,然后用有限差分算子近似微分算子,得到差分方程.因此高精度有限差分算子的求取和误差估计可以说是有限差分模拟的核心.目前数值模拟中常用的有限差分数值模拟可以分为二阶波动方程(Dablain,1986;Kneib,1993)和速度2应力一阶方程组(virieux,1984,1986)两种. Levander(1988)发展了交错网格格式的四阶差分格式,使得模拟精度有了很大提高.但是经典交错网格格式存在本身固有的缺点,如图1所示,拉梅常数定义在所有的半网格点和整网格点上,但是实际中通常只定义在半网格点上;对于切应力的计算,需要对拉梅常数进行插值或者用周围的值来近似,如果变化很大时,就会出现计算的不稳定.在自由边界处,由于固体和空气性质的强对比性,就需要引入专门处理边界的问题(Graves,1996;Hest holm&ruud, 1998;Opral&zahradhik,1999),也带来许多不便. Igel(1995)分析了交错网格格式的缺点,此后又有一系列文章指出交错网格存在的不足之处(K oma2 tit sch,2002;Carcio ne,2002).Saenger(2000)年研究了该问题出现的原因,提出了旋转交错网格格式(RSG(Rotated Staggered grid))并将这种格式应用于各种不同的模型.相对于交错网格格式,RSG可以得到更稳定、更可信的解,在自由界面处使用与内部相同的差分格式来处理不会引起数值不稳定.如图2所示,RSG只要沿着坐标轴方向作差分来求波函数的微分值即可得到稳定的解.3　震源函数子波是震源的时间函数,描述震源的时间延续特征.对于地震子波而言,子波延续时间越短,频带越宽,地震子波的垂直分辨率就越高.但是有限差分模拟很大的问题就是数值频散,子波中的高频成分对网格间距很敏感,当空间采样不足时,高频成分频散很严重.因此要根据模型的速度参数和网格间距选取子波主频.常用的地震子波有Ricker子波、Gauss子波及其导数.Ricker子波是零相位的,零相位子波可以达到分辨率的极限.任义庆(1998)模拟了从爆炸到地震子波形成的过程,对于研究地震子波的频率变化有一定的意义.震源函数给定通常有两种方法,一种是用理论结果作为初始值来给出,即初值法;另一种是以力源8842期冯英杰,等:地震波有限差分模拟综述的方式给出,即力源法.这两种方法各有优势.初值法避开了震源位置的奇异性,可以定义在模型任意处,但是震源却不能放在自由表面或内界面附近.力源法,震源虽然可以定义在自由表面附近,但是必须在网格点处.在速度-应力方程组中,是将震源赋在两个法向应力处来模拟点爆炸震源,而不是赋在速度处,这样就很好的避开了震源处的无穷速度问题(Virieux,1986).董清华(2000)介绍了胀缩力源、剪切力源和方向力源的给定方法.Graham.J.Hick (2002)论述了震源函数的模拟,给出了震源函数的最佳窗函数的形式,最优的逼近了实际震源的效果.4　边界条件在有限差分数值模拟中,计算区域是有限的,不可能模拟无限区域的情况,因此有限差分数值模拟的一个重要问题就是人工边界处理.如果在模型边界直接采用刚性边界即位移为0,或者自由边界即应力为0,两种边界都是完全反射边界,即反射系数绝对值都是1,都会导致严重的边界反射,破坏有效区域的数值解.目前主要有5种方法用于消除模型边界效应, (1)运动边界条件,即计算区域随计算时间的推移而扩大,在计算时间内波不能传到介质的边界.可以想象该方法一定可以很好的模拟无限边界,但是其对于内存和机时的需求也是可观的;(2)Smit h边界条件(Smit h,1974),即综合Neumann边界条件和Dirichlet边界条件,因为在Neumann边界,介质的反射系数是+1,Dirichlet边界上介质的反射系数是-1,将这两种边界上的反射结果相加,则得到无反射的波场.这种边界条件对于消除一次波的效果比较理想,对于多次波效果很差,而且随着边界数目的增加,计算量也迅速增大;(3)吸收边界条件(Clayton&Engquist,1977;Engquist&Majda, 1977;Reynolds,1978;Keys,1985;Hidgon,1987; Long,1990;Hagst rom,1997),即在边界处,运用单程波方程来计算波场,由于单程波方程的导出有其自身的假设条件,所以这种方法对于垂直入射波吸收效果较好,而对于大角度入射波吸收效果则不理想;(4)加吸收层技术(Cerjan et al,1985;K osloff. R&K o sloff.D,1986;Sochacki,1987;),也称吸收边界,即在模型以外,增加多层网格,对波函数值进行衰减,目前最佳吸收层技术(Berenger,1994; peng,1994,1995;Hasting,1996;)堪称是该类方法中的首选,但是这类方法的缺点是计算量和存储空间增加;(5)波场外推法,这种方法最先是由Jianlin zhu(1999)在Geop hysics上提出的,他把它称为透明边界条件.该方法是利用模型内部的数值计算结果,根据同一波前面上的质点具有相同的振动相位和波传播过程中的振幅变化规律,计算得到边界上的波函数值.罗大清(2000)将该方法用于消除模型的角点反射,田小波(2004)改进了这一技术,在理论计算中都取得很好的效果.对于起伏自由表面的处理是目前处理的关键, Erik H.Saenger(2000,2004)提出从自由表面开始按一定的函数形式把介质划分为不规则网格,通过数学变换,将不规则网格变换为规则网格,在规则网格上计算波场.但是这种方法只能处理一阶可导的光滑自由表面.陈伟(2005)用渐变的速度模型进行了起伏地表的模拟.5　数值稳定性和频散消除数值稳定条件是显式有限差分格式必须要分析的问题,波动方程有限差分格式一般都是按时间逐层推进的,这样前一时间波函数值的舍入误差必然影响到后一时间的波场.这就有必要分析误差传播和积累情况,使误差不至于随时间的推进而迅速增长,破坏整个数值解,甚至导致计算溢出.根据Lax 等价定理,稳定性也保证了差分格式的收敛性.稳定性分析方法一般是利用Von Neumann提出的Fou2 rier谱分析方法,影响稳定性的关键参数就是网格比p=Δt/Δx.董良国(2000)进行了交错网格高阶差分的稳定性研究.在实际介质中,地震体波的频散并不明显(谢里夫等,1999).波动方程有限差分数值解可以理解为波在离散化的网格上以差分格式传播,这种离散直接导致各个频率成分传播速度不同,一般是高频成分相速度明显下降,因此可以说网格频散是有限差分的固有数值问题,当网格大小不合适时,会表现出严重的频散现象,在合成记录上可以看到主要震相之后有很长的拖尾,降低了分辨率,主波长上的网格点数以及差分格式精度是影响合成记录的关键因素.压制频散最简单的办法就是减小网格步长.蔡其新等(2003)曾经研究了优化差分参数的一种公式,用来确定空间步长.其他的还有高阶差分格式(Fornberg,1987;吴国忱,2005),通量校正传输法(FC T)(Fei,1996).Fornberg对比高阶有限差分和伪谱法后指出,当有限差分算子的阶数逼近无穷时,等价于伪谱法,逼近阶数越高,模拟的数值频散越984地　球　物　理　学　进　展22卷小.FC T是Boris(1973)研究流体运移问题提出的方法,Fei将其用于消除数值频散,其基本原理是假设所有的极值点都是由数值频散引起的,传统的FCT方法对所有网格点进行通量校正处理,再对局部极值点进行补偿的逆通量校正.Tong Fei(1995)提出了优化的FC T,通量校正只用在局部极大值上,节省了大约40%的计算量.同时FCT方法可以放大时间和空间步长,从而抵消计算FC T带来的计算量的增加.6　结　论有限差分法数值模拟是数值模拟中一种很重要的方法,该方法在理论研究和实际应用中发挥着越来越重要的作用.但是数值模拟作为一门博大精深的学问,无论在理论上还是实际应用中需要突破的地方还很多.本文作者仅就自己的研究领域所涉及的范围做了一些论述.参　考　文　献(References):[1]　王秀明,张海澜,王东.利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播[J].地球物理学报,2003, 46(6):842～849.[2]　王德利,何樵登,韩立国.裂隙型单斜介质中多方位地面三分量记录模拟[J].地球物理学报,2005,48(2):386～393.[3]　张剑锋.弹性波数值模拟的非规则网格差分法[J].地球物理学报,1998,41(增刊):357～365.[4]　张剑锋.各向异性介质中弹性波的数值模拟[J].固体力学学报,2000,21(3):234～242.[5]　张剑锋,刘铁林.三维非均匀介质中弹性波传播的数值模拟[J].固体力学学报,2001,22(4):356～360.[6]　杨顶辉,滕吉文,张中杰.三分量地震波场的近似解析离散模拟技术[J].地球物理学报.1996,39,(增刊):283～291.[7]　杨顶辉.各向异性介质弹性波方程的正反演方法研究[D].北京:中国科学院地质与地球物理所,1996.[8]　任义庆,李勤学,马在田.地震波爆炸震源模拟[J].石油物探,1998,37(3):15～21.[9]　董清华.震源数值模拟[J].世界地质工程,2000,16(3):27～32.[10]　罗大清,宋炜,吴律.一种有效的处理模型角点反射的方法[J].石油物探,2000,39(4):26～31.[11]　田小波,吴庆举,曾融生.弹性波数值模拟的延迟边界方法[J].地球物理学报,2004,47(2):268～273.[12]　陈伟.起伏地表条件下二维地震波场的数值模拟[J].勘探地球物理进展,2005,28(1),25～31.[13]　董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报,2000,43(6):856～864.[14]　谢里夫,吉尔达特著,初英,等译.勘探地震学(第二版)[M].北京:石油工业出版社,1999.[15]　蔡其新,何佩军,秦广胜,等.有限差分数值模拟的最小频散算法及其应用[J].石油地球物理勘探,2003,38(3).247～251.[16]　吴国忱,王华忠.波场模拟中的数值频散分析与校正策略.[J]地球物理学进展,2005,20(1):58～65.[17]　常旭,刘伊克.地震正反演与成像[M].北京:华文出版社,2001.[18]　牛滨华,孙春岩.地震波理论研究进展———介质模型与地震波传播[J].地球物理学进展,2004,19(2):255～163. [19]　王红落.地震波传播与成像若干问题的研究[D].北京:中国科学院地质与地球物理所,2004.[20]　孙若昧.地震波传播有限差分模拟的人工边界条件[J].地球物理学进展,1996,11(3):53～58.[21]　何兵寿,魏修成,刘洋.三维波动方程的数值频散关系及其叠前和叠后数值模拟[J].石油大学学报(自然科学版),2001,25(1):67～71.[22]　黄自萍,张铭,吴文清,等.弹性波传播数值模拟的区域分裂法[J].地球物理学报,2004,47(6):1094～1100..[23]　裴正林,牟永光.地震波传播数值模拟[J].地球物理学进展,2004,19(4):933～941.[24]　Alterman Z,Karal F C.Propagation of seismic wave in lay2ered media by finite difference met hods[J].BSSA,1968,58(1):367～398.[25]　Boore D M.Finite2difference met hods for seismic wave prop2agation in heterogeneous materials in Met hods in computa2tional physics[J].1972,11: B. A.Bolt,ed.,AcademicPress,inc.[26]　Kelly K R,Ward R W,Treitel S,Alford R M.Synt hetic seis2mograms2a finite2difference approach[J].Geophysics,1976,41:2～27.[27]　Madariaga R.Dynamics of an expanding circular fault[J].Bull Seism Soc Am,1976,65:163～18.[28]　Virieux J.S H wave propagation in heterogeneous media:ve2locity2st ress finite difference met hod[J].Geophysics,1984,49(11):1933～1957.[29]　Virieux J.P2SV wave propagation in heterogeneous media:velocity2stress finite difference met hod[J].Geophysics,1986,51(4):889～901.[30]　Igel H,Mora P,et al.Anisotropic wave propagation t hroughfinite2difference grids[J].Geophysics,1995,60(4):1203～1216.[31]　Igel H,weber M.P2SV wave propagation in t he Eart h’smantle using finite differences:application to heterogeneouslowermost mantle structure[J].Geophys,Res.Lett,1996,23:415～418.[32]　J astram C,Behle A.Acoustic modeling on a grid of velocityvarying spacing[J].Geophysical prospecting,1992,40:157～169.[33]　J astram C,Tessemer E.Elastic modeling on a grid wit h ver2tically varying spacing[J].Geophysical prospecting,1994,42:357～370.[34]　Falk J,Tessmer E,Gajewski D.Tube wave modeling by t hefinite2difference met hod wit h varying grid spacing[J].Pa20942期冯英杰,等:地震波有限差分模拟综述geoph,1996,148:77～93.[35]　Falk J,Tessmer E,Gajewski D.Efficient finite2differencemodeling of seismic waves using locally adjustable time steps[J].Geophysical Prospecting,1998,46:603～616.[36]　Tessmer E.seismic finite difference modeling wit h spatiallyvarying time step[J].Geophysics,2000,65(4):1290～1293.[37]　Carcione J M.wave fields in real media:wave propagation inanisot ropic an elastic and porous media[J].U K:Elsevier Sci2ence L TD,2001.[38]　Dablain M A.The application of high2order differencing tot he scalar wave equation[J].Geophysics,1986,51(1):54～66.[39]　Kneib G,Kerner C.Accurate and efficient seismic modelingin random media[J].Geophysics1993,58(4):576～588. [40]　Levander A.R Fourt h2order finite2difference P2SV seismo2grams[J].Geophysics,1988,53(11):1425～1436.[41]　Graves R W.Simulating seismic wave propagation in3D elas2tic media using staggered2grid finite2difference[J].BSSA,1996,86:1091～1106.[42]　Hest holm S O,Ruud B O.32D finite2difference elastic wavemodeling[J].Geophysics,1998,63:613～622.[43]　Oprsal I A,Zahradnik J.Elastic finite2difference met hod forirregular grids[J].Geophysics1999,64:240～250.[44]　K omatit sch D,Barnes C,Tromp J.Simulation of t he diffrac2tion by single cracks:An accuracy st udy.72nd Annual Inter2national meeting,SEG,Abstract s,2002,2007～2010. [45]　Carcione J M,et al.Seismic Modeling[J].Geophysics,2002,67(4):1304～1325.[46]　Saenger E H,G old N,Shapiro S A.Modeling t he propaga2tion of elastic waves using a modified finite2difference grid[J].Wave Motion,2000,31:77～92.[47]　Hicks G J.Arbitrary source and receiver positioning in finite2difference scheme using Kaiser windowed sinc functions[J].Geophysics,2002,67(1):156～166.[48]　Smit h W D.A nonreflecting plane boundary for wave propa2gation problems[J].J Comp Phys,1974,15:492～503. [49]　Majda E B.A Absorbing boundary conditions for t he numeri2cal simulation of wave[J].Mat h Comp,1977,629～651. [50]　Clayton R,Engquist B.Absorbing boundary conditions for a2coustic and elastic wave equation[J].BSSA,1977,67:1529～1540.[51]　Reynolds A C.Boundary conditions for t he numerical solu2tion of wave propagation problem[J].Geophysics,1978,43:1099～1110.[52]　Keys R G.Absorbing boundary conditions for acoustic media[J].Geophysics,1985,50:892～902.[53]　Higdon R L.Numerical absorbing boundary conditions fort he wave equation[J].Mat h Comp,1987,49:65～90. [54]　Long L T,Liow J S.A transparent boundary for finite2difference wave simulation[J].Geophysics,1990,55:201～208.[55]　Hagstrom T.On high2order radiation boundary condition,inEngquist B,Kriegsmann G A Eds[J].Computational WavePropagation:Springer2Verlag New Y ork,Inc,1997,86:1～21.[56]　Cerjan C,K osloff D,K osloff R,et al.A nonreflectingboundary condition for discrete acoustic and elastic wave e2quations[J].Geophysics,50:705～708.[57]　K osloff R,K osloff D.Absorbing boundary conditions forwave propagation problems[J].J Comp Phys,1986,63:363～376.[58]　Sochacki J,Kubichek R,et al.Absorbing boundary condi2tions and surface wave[J].Geophysics,1987,52:60～71. [59]　Berenger J.A perfectly matched layer for t he absorption ofsorption of electromagnetic wave[J].J Comput Phys,1994,114:185～200.[60]　Peng C,Toksoz M N.Optimal absorbing boundary condi2tions for finite2difference modeling of acoustic and elasticwave propagation[J].J Acoust Soc Am,1994,95:733～745.[61]　Peng C,Toksoz M N.An optical absorbing boundary condi2tion for elastic wave modeling[J].Geophysics,1995,60:296～301.[62]　Hasting F,Schneider J B,et al.Application of t he perfectlymatched layer(PML)absorbing boundary condition to elasticwave propagation[J].J Acoust Soc Am,1996,100:3061～3069.[63]　Jianlin Zhu.A transparent boundary technique for numericalmodeling of elastic waves[J].Geophysics,1999,64(3):963～966.[64]　Saenger E H,Shapiro S.A Effective velocities in fracturedmedia:A numerical study using t he rotated staggered finite2difference grid[J].Geophysical Prospecting,2002,50:183～194.[65]　Saenger E H,Thomas Bohlen.Finite2difference modeling ofviscoelastic and anisotropic wave propagation using t he rota2ted staggered grid[J].Geophysics,2004,69(2):583～591. [66]　Fornberg B.The pseudo2spectral met hod:comparisons wit2hfinite differences for t he elastic wave equation[J].Geophys2 ics,1987,52(4):483～501.[67]　Fei T,larner K.Elimination of numerical dispersion in finite2difference modeling and migration by flux2corrected transport[J].Geophysics,1995,60(6):1830～1842.[68]　Boris J,Book D.Flux2corrected transport.I.SHASTA,Afluid transport algorit hm t hat works[J].J Comput.Phys,1973,11:38～69.[69]　Muller T M,Shapiro S A.Most probable seismic pulses insingle realizations of two2and t hree2dimensional random media[J].Geophysical Journal International,2001,144:83～95.[70]　Muller T M,Shapiro S A,Sick C M.A most probable ballis2tic waves in random media:A weak2fluctuation approximationand numerical result s[J].Waves in Random media,2002,12:223～245.194。

地震声波数值模拟中的透射边界条件
邢丽
【期刊名称】《科学技术与工程》
【年(卷),期】2009(009)013
【摘要】地震波传播的数值模拟中,在有限区域建立吸收边界条件是一个很重要的问题.在不同类型的吸收边界条件中,比较经典地是Clayton和Engquist提出得Clayton-Engquist吸收边界条件,透射边界条件(MTF)是通过直接模拟各种单向波动的共同运动学特征来建立的一种边界条件,算法简单容易.重点讨论了MTF边界条件并通过数值计算与被传统应用的CE吸收边界条件作比较.
【总页数】5页(P3607-3611)
【作者】邢丽
【作者单位】上海第二工业大学数学系,上海,201209
【正文语种】中文
【中图分类】P315.31
【相关文献】
1.地震声波数值模拟中的吸收边界条件 [J], 邢丽
2.声波数值模拟中改进的非分裂式PML边界条件 [J], 熊章强;毛承英
3.基于优化有限差分和混合吸收边界条件的三维VTI介质声波和弹性波数值模拟[J], 徐世刚;刘洋
4.透射边界条件在点声源海底地震波数值计算中的应用 [J], 卢再华;张志宏;顾建农
5.地震声波数值模拟中人工边界条件的差别与组合 [J], 韩复兴;王若雯;孙章庆;高正辉
因版权原因，仅展示原文概要，查看原文内容请购买。

二维地震波场的组合型紧致有限差分数值模拟汪勇;石好果;周成尧;桂志先【摘要】地震波场数值模拟在地球物理勘探和地震学中具有重要的支撑作用.本文将组合型紧致差分格式用于声波和弹性波方程的数值模拟中.根据泰勒级数展开和声波方程,建立了位移场时间四阶离散格式,并将组合型紧致差分格式用于位移场空间导数的求取,然后对该差分格式进行了精度分析、误差分析、频散分析和稳定性分析.理论研究结果表明:① 该差分格式为时间四阶、空间六阶精度,与常规七点六阶中心差分和五点六阶紧致差分相比,具有更小的截断误差和更高的模拟精度;② 每个波长仅需要5.6个采样点,且满足稳定性条件的库郎数为0.792,可以使用粗网格和较大时间步长进行计算.所以该方法具有占用内存少、计算效率高和低数值频散等优势.最后,本文进行了二维各向同性完全弹性介质的声波和弹性波方程的数值模拟,实验结果表明本文提出的方法具有更高的计算精度,能够大幅度的节约计算量和内存需求,对于三维大尺度模型问题具有更好的适应性.【期刊名称】《地球物理学报》【年(卷),期】2018(061)011【总页数】16页(P4568-4583)【关键词】组合型紧致差分;声波方程;弹性波方程;数值模拟;数值频散;稳定性条件【作者】汪勇;石好果;周成尧;桂志先【作者单位】油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100;胜利油田勘探开发研究院西部分院,山东东营 257001;油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100;油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100【正文语种】中文【中图分类】P6310 引言地震波场数值模拟是应用数值方法模拟地震波在地下介质中的传播过程，计算在地面或地下各观测点地震记录的一种数值模拟方法，是地球物理勘探和地震学的重要基础.目前，常用的数值模拟方法主要有射线追踪法和波动方程法，其中波动方程法有伪谱法、有限元法、边界元法、谱元法和有限差分法(Alterman and Karal，1968；Graves，1996；董良国等，2000a；Saenger et al.，2000，2004).随着数值模拟技术的发展和生产实践的要求，围绕着提高有限差分计算效率(黄超和董良国，2009；唐佳等，2016)、模拟精度(Liu，2013；李振春等，2016)、算法稳定性(Virieux, 1986；董良国等，2000；杜启振等，2015)、压制数值频散、处理复杂介质(Fletcher et al.，2008；王璐琛等，2016；程玖兵等，2013)和吸收边界条件(Berenger，1994；Pan et al.，2012)等方面，研究人员已经发展了多种方法.特别是在如何压制数值频散和提高计算效率方面，杨顶辉等在这方面做了大量相关工作(Yang et al.，2009，2012a,b，2014；Yang and Wang, 2010，Ma and Yang，2017；Tong et al.，2013；Zhou et al.，2015；He et al.，2015；Liu et al.，2017；Jing et al.，2017)，取得了许多有意义研究成果.提高差分格式数值计算精度最直接的方法就是在差分计算时增加网格节点数，但这也增加了计算量和所需的存储空间.紧致差分方法恰好能够较好地解决这个矛盾，同时紧致差分是一种隐式差分格式，具有较好的稳定性，这些优势也使得它成为目前研究较多的有限差分方法之一.1989年，Dennis和Hudson(1989)针对Navier-Stokes方程提出了空间四阶的紧致格式，能够使用粗网格获得较高的精度.1992年，Lele(1992)构造了求解一阶导数和二阶导数的紧致差分格式.Adams 和Shariff(1996)提出了紧致ENO格式，用于求解激波湍流相互作用问题.Chu和Fan(1998)提出了三点组合型紧致差分格式，并将其用于求解对流扩散方程.随后人们针对不同的问题和方程，提出了许多种不同的紧致差分格式，广泛用于空气动力学、流体力学和电磁波方程等方面.在地震波场数值模拟方面，王书强等(2002)和Zhou和Greenhalgh(2003)先后将Lele(1992)提出的紧致差分格式用于弹性波方程的数值模拟中.Yang等(2003)年给出各向异性情况下的紧致差分方法以及声波和弹性波数值模拟结果，并给出相应的吸收边界条件.Du等(2009a,b)利用紧致交错网格差分方法对横向各向同性介质地震波场进行了数值模拟，Liu和Sen(2009)研究了任意阶空间导数的紧致格式和一阶空间导数的交错紧致格式，分析了它们的精度和计算效率，并将其用于声波和弹性波数值模拟中.Kosloff等(2010)提出可以利用当前点两边各任意多个点计算导数值的一般紧致格式，并基于Remes方法求取差分系数.杨宽德等(2011)研究了Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟方法.此外，Chu和Stoffa(2010)、Liu(2014)还研究了频率域地震波数值模拟中的紧致差分方法.而组合型紧致差分格式在地震波场模拟中的应用，尚未见到文献报道.本文在声波和弹性波方程时间四阶离散格式的基础上，将组合型紧致格式应用到位移场空间偏导数的求取，实现了二维各向同性介质地震波场的数值模拟.1 组合型紧致有限差分方法原理早期的紧致差分格式是基于Hermite多项式构造而来的，在Hermite多项式中，相邻三个节点的函数值、一阶和二阶导数值关系可以表示为其中fi表示节点函数值，和分别表示一阶和二阶导数值.1992年，Lele(1992)对Hermite公式进行了扩展，构造了求解一阶导数和二阶导数的紧致差分格式，表示为(2)其中a,b,c,α,β为差分系数，h表示网格间距.对上述Lele差分格式进行泰勒公式展开和求解差分系数，可以得到一阶导数五点中心六阶精度差分格式：(3)二阶导数五点中心六阶精度差分格式：(4)由公式(3)和(4)可以看出，紧致差分(Compact Difference，以下简称CD)格式的特点就是用相邻节点的函数值求解若干节点上的导数值.而常规中心差分(Central Finite Difference，以下简称CFD)仅求解中心节点的导数值.另一方面，也可以认为常规中心差分是紧致差分的特例，即公式(2)中的α=β=0，例如一阶和二阶导数的六阶精度差分格式表示为(5)其中b1=0.75, b2=-0.15, b3=0.0167, c0=-2.7222, c1=1.5, c2=-0.15,c3=0.0111.常规中心差分如果要得到六阶精度的差分格式，需要用到七个节点，而紧致差分只需要五个点，这使得紧致差分格式比常规差分处理边界节点上更为方便.且一般情况下，在相同网格间距时，紧致差分格式比常规差分具有更高的精度，具有更小的数值频散(王书强等，2002).Chu和Fan(1998)等构造了精度更高的三点六阶组合型紧致差分格式(Combined Compact Difference，以下简称CCD)：(6)与公式(3)和(4)对比可以看出，上述组合型紧致差分格式只需要相邻的三个节点就可以同时求得一阶和二阶导数的六阶精度近似值，比常规紧致差分的节点数更少.公式(6)中的一阶导数和二阶导数是耦合的，既可以同时求出也有利于波形的保真性.2 声波方程组合型紧致有限差分方法分析二维声波方程可以表示为(7)式中u(x,z)为地震波位移，v(x,z)为地震波速度.利用截断的泰勒公式表示n+1和n-1时刻的位移场可以得到：(9)两式相加，略去高次项，并代入声波方程，得到位移场时间四阶精度的差分格式：(10)公式(10)为位移场的三层显式差分格式，利用它就可以进行地震波场时间层的推进.在公式中含有位移对x和z的二阶和四阶导数，对这些导数将利用组合型紧致差分格式(6)进行求取.假设模型纵向和横向节点数为m和n，h为空间网格大小.利用公式(6)求公式(10)中的空间偏导数∂2u/∂z2和∂2u/∂x2，矩阵表示为AE=BU， FC=UD,(11)其中A和C为公式(6)左端的差分系数矩阵，大小分别为2m×2m和2n×2n；E和F为待求位移场空间一阶和二阶导数矩阵，大小分别为2m×n和m×2n；B和D为公式(6)右端的差分系数矩阵，大小分别为2m×m和n×2n；U为位移场矩阵，大小为m×n.这些矩阵分别表示为(12)(13)(15)(16)(17)(18)公式(11)可以用追赶法进行求解.由公式得出位移场ui,j可同时求得它的空间一阶和二阶导数∂ui,j/∂x、∂ui,j/∂z、∂2ui,j/∂x2和∂2ui,j/∂z2，表示为E=A-1BU， F=UDC-1,(19)差分公式(10)中的空间四阶导数可以在求得二阶导数后，将二阶导数值当作U，再次使用公式(19)进行求取.2.1 精度分析不论是利用CCD格式还是利用常规的七点六阶CFD或五点六阶CD格式，进行声波方程数值模拟时，它们在时间层推进方式上是相同的，即都是利用差分格式(10)，时间差分为四阶精度.CD、CFD和CCD不同之处在于它们分别利用公式(4)、(5)和(6)计算位移空间高阶导数.虽然在声波方程的差分格式中并没有空间一阶偏导数，但在下文中的第4节进行弹性波场数值模拟中，存在空间的混合偏导数，需要通过一阶导数的迭代求取，所以这里对这三种格式计算的一阶和二阶导数的截断误差进行对比，结果见表1，α、β、a、b和c为公式(2)中的差分系数.从表1可以看出，虽然三种方法都能达到空间六阶精度，但截断误差有较大的差别.CFD和CD差分格式计算一阶偏导数的截断误差约是CCD的40倍和4.4倍，计算二阶偏导数的截断误差约是CCD的36倍和8.5倍，这说明CCD方法的具有更高的差分精度.表1 一阶和二阶导数不同差分方法的截断误差比较Table 1 Truncation errors in various difference schemes for the first and second derivative calculationsαβabc截断误差uxCFD003/2-3/51/1036/(7!)×f(7)h6CD1/3014/91/904/(7!)×f(7)h6CCD/≈-0.9/(7!)×f(7)h62ux2CFD003/2-3/51/1072/(8!)×f(8)h6CD2/11012/113/110≈-17/(8!)×f(8)h6CCD/-2/(8!)×f(8)h62.2 误差分析通过模拟二维平面谐波初值问题并计算模拟误差，定量分析CCD方法的数值模拟精度.二维平面波初值问题可以表示为(20)其中v是平面波的波速，θ 0是初始时刻波阵面法线方向(即传播方向)与x轴的夹角，f0是平面简谐波的频率.上述初值问题的解析解为(21)在二维波场数值模拟中，假设f0=20 Hz，θ0=π/4，均匀介质的波速为3600 m·s-1，模型长度和深度均为2000 m，纵横向网格长度相同，采样时间为1 s.在不同的空间网格长度和时间步长条件下，计算数值模拟的相对误差.相对误差定义为(22)式中表示数值解，u(tn,xi,zj)表示解析解.研究表明(Yang et al.，2014)，优化的近似解析离散化方法(Optimal Nearly-Analytic Discrete，以下简称ONAD)也具有较高的数值模拟精度和计算效率，为了比较，在此计算CCD、CD、ONAD和CFD四种方法在不同模拟参数情况下的相对误差曲线见图1.从图1可以看出，随着空间网格长度和时间步长的增加，四种方法的相对误差均会逐渐增加，并且随模拟时间的增加而增加.便于比较，将实验结果列于表2.从表中可以看出，四种情况下，误差相对关系基本一致，例如当空间步长为20 m，时间步长1 ms时，CCD算法的相对误差仅为0.066%，CD次之，为0.089%，CFD相对误差为0.943%，这与精度分析的结论是一致的.需要说明的是，本次实验的ONAD算法相对误差比CFD方法大，这是由于ONAD只有空间四阶精度，而本文的CFD空间精度为六阶.当采用较细网格时，如空间步长10 m，时间步长0.5 ms时，四阶ONAD方法的误差要小于六阶CFD方法，这也说明了ONAD方法在提高模拟精度方面也有积极的意义.用较小的空间步长(10 m)和时间步长(0.5 ms)时，CCD算法的精度显著提高，相对误差仅有0.0008%，并且相对误差增长缓慢，这说明了采用CCD格式的声波模拟结果精度较高，能进行较长时间的地震波场数值模拟.图1 四种数值模拟算法的相对误差曲线(a) CCD; (b) CD; (c) ONAD; (d) CFD.Fig.1 Relative error-time curves of the numerical simulation using four different methods and models(a) CCD; (b) CD; (c) ONAD; (d) CFD.表2 四种方法在不同情况下的最大相对误差(%)比较Table 2 Relative errors in four different schemes and parameters模拟参数CCDCDONADCFDdx=10 m,Δt=0.5ms0.00080.00160.1190.596dx=15 m,Δt =0.5 ms0.0120.0150.690.62dx =15m,Δt =1 ms0.0080.0190.5560.629dx =20 m,Δt =1 ms0.0660.0892.160.943 2.3 频散分析在数值模拟时，如果空间网格长度使用过大，会造成较大的求解误差，并会产生所谓的数值频散现象，所以频散关系分析既是判断数值模拟方法优劣的重要方法，也是确定空间网格大小的重要依据.通过对声波方程的组合型紧致差分格式进行数值频散分析，进一步分析该方法的适用条件.首先考虑x方向，令简谐波解=Aexp×[I(ihkx+jhkz-kv′nΔt)]，且：(23)如果数值模拟时不存在误差，则A,B,C应该满足以下关系：B=(Ikx)A, C=-(kx)2A,(24)但在差分的实际计算中，所产生的数值波数(又称修正波数)与解析波数不同.令：(25)将公式(23)和(25)代入差分格式(6)，可以得到方程组：(26)其中：解上述方程组可得：(27)同理z方向可得：(28)式中φx=φcosθ，φz=φsinθ，且φ=kh.θ是波的传播方向与x轴的夹角，k为波数，h为空间步长.将公式(23)、(27)和(28)代入声波方程的差分格式(10)中，并令x=kv′Δt，α=vΔt/h.其中v′为数值波速，v是实际波速，Δt是时间步长，α称为Courant数，利用欧拉公式，化简可得CCD频散关系：]]/12.(29)通过上述频散关系确定φ后可以解得对应的x.定义数值波速与真速度的比值为(30)在压制数值频散方面，优化的近似解析离散化方法(ONAD)方法也具有较好的效果，所以选取CFD、CD和ONAD 三种方法与之进行比较.在理想情况下，如果不存在数值频散则速度比γ恒等于1.γ偏离1越大，则说明该方法的数值频散越严重，反之则说明该方法能更好地压制数值频散.图2显示了以上四种方法在不同的α和θ的时候的速度比曲线.取φ∈[0,π]作为横坐标，它是波数与空间步长的乘积，单位波长内采样点数N=2π/φ，所以横坐标也可以看作N由逐渐减小至2.从图中的速度比曲线可以看出：①随着空间采样点数的减小，四种方法的频散现象逐步加剧.CCD、CD和ONAD方法的数值频散比CFD方法要小，其频散曲线更趋近于1，这也说明了紧致差分(CCD和CD)和近似解析离散类(ONAD)方法在压制数值频散方面的都具有较大的优势，都适合使用较大的空间步长，进而提高计算效率；②假设数值速度在理论速度的99.9%～100.1%以内表示不存在频散，则CCD、CD、ONAD和CFD 对应的φ的极小值分别为1.12、0.98、1.06和0.91，所以它们在最小主波长内需要使用的网格点数分别为5.6、6.4、5.9和6.9个，即CCD和ONAD方法对于空间网格长度的要求最低，其次是CD方法，CFD方法最为严格；③CCD方法的空间网格长度过大时，速度比曲线上翘，数值速度大于真速度，这与另外三种方法相反.这一结论也能从图3c中可以看出，当φ=π，即一个周期内取2个样点，蓝色的速度比曲线大部分要大于1.图2 不同差分方法的速度比曲线(a) θ=0,α=0.25; (b) θ=π/8,α=0.25; (c)θ=π/4,α=0.25; (d) θ=0,α=0.45; (e) θ=π/8,α=0.45; (f) θ=π/4,α=0.45.Fig.2 Velocity ratio curves of the numerical simulation in different methods图3 不同差分方法的方位角-速度比曲线(a) CFD; (b) CD; (c) CCD; (d)ONAD.Fig.3 Azimuth-speed ratio curves of the numerical simulation in different methods图4 空间网格32 m，时间步长1 ms，420 ms时刻波场快照(a) CFD ; (b) CD; (c) CCD; (d) ONAD.Fig.4 Snapshots at 440ms: Δx=32 m，Δt=1 ms为研究该差分方法的数值各向异性特征，即速度比与传播方向的关系，绘制极坐标速度比曲线，如图3所示.图中的极角表示地震波传播方向与x轴的夹角，极径表示速度比.四种方法在空间采样(φ=π/3)满足频散要求时，都没有各向异性特征.在空间采样不足时(φ=7π/5，φ=9π/5)均会出现速度各向异性，0°(或90°)与45°方向上的速度差别最大.利用二维模型来说明CCD方法在压制数值频散方面的优势.模型长度和深度均为4 km，纵波速度3600 m·s-1，震源子波为f(t)=sin(2πf0t)exp×主频20 Hz.根据频散分析结论，CCD方法在最小主波长内只需要使用约5.6个点，所以设置空间网格间距32 m，时间步长1 ms.图4为四种方法数值模拟的波场快照，可以看出CCD和ONAD方法的结果清晰，没有频散现象，CD方法在0°和90°方向上存在微弱的频散现象，而CFD的频散非常严重，与前面理论分析结果一致.2.4 稳定性分析稳定性条件也是有限差分数值模拟中一个非常重要的问题，是影响差分方法计算效率的重要因素.采用Fourier方法(张文生，2006)对声波方程的组合型紧致差分格式进行稳定性分析.定义：(31)其中为空间网格大小，kx和kz为x和z方向的视波数，i、j和n分别为空间和时间网格坐标，代入公式(6)中可得：(32)令θ1=kxh，-π/2≤θ1≤π/2，并利用欧拉公式，将上式化简可得：(33)解方程可得：=Aξn/h2,(34)其中(35)同理，定义：=eIh(kxi+kzj), =eIh[kxi+kz(j-1)],=eIh[kxi+kz(j+1)], =eIh(kxi+kzj),=eIh[kxi+kz(j-1)],=eIh[kxi+kz(j+1)],且令θ2=kzh， -π/2≤θ2≤π/2，代入公式(6)可以得到：=B ξn/h2,(37)其中(38)公式(10)是一个三层显式差分格式，为了分析其增长矩阵，将其改写为：(39)式中=φneIh(kxi+kzj).将公式(34、35、37和38)代入(∂2u/∂x2=exp[Ih(kxi+kzj)]和(∂2u/∂z2=exp[Ih(kxi+kzj)]，并代入声波差分格式(39)中，化简并写成矩阵格式为(40)其中(41)为增长矩阵，其中α=vmaxΔt/h为该差分格式的Courant数.为使差分格式(40)满足稳定性条件，需要其增长矩阵的谱半径小于等于1.直接求解该增长矩阵的特征值较为复杂，用二分法求该问题的解(宋国杰，2008)，计算得到的近似稳定性条件为α=vmaxΔt/h<0.792.2.5 计算效率分析采用2.3节所用模型进行数值模拟，采样时间1 s，分析CCD、CFD、CD和ONAD四种方法的计算效率，各方法使用的数组大小、个数和算法时间复杂度等列于表3.表3 不同方法计算效率比较Table 3 Comparison of computational efficiency between different methodsCCDCFDCDONAD网格半径3753每个波长采样点5.66.96.45.9Δx=Δz/(m)32262831数组/个1281225Courant0.790.980.930.527Δt/ms7774时间复杂度O(n3)O(n3)O(n2)O(n2)时间/s0.6670.8681.041.79从表3中可以看出：CCD和ONAD差分方法的网格半径最少，仅需要3个点，这对于边界处理最有利，但ONAD空间精度仅为四阶，要低于CCD的六阶精度.CFD和CD差分方法虽然空间精度能达到六阶，但网格半径也随之增加.CCD和ONAD均为低数值频散的差分方法，满足无数值频散要求的空间网格长度大于CFD和CD方法.但是ONAD方法需要计算位移场的梯度，这就增加了数组的个数和计算时间.同样的声波方程，CCD仅需要12个数组，而ONAD需要25个，在网格长度相近的情况下，CCD方法占用内存最小.从稳定性来看，要达到同样的时间四阶精度，CFD和CD方法的Courant数最大，CCD次之，ONAD最小.考虑空间步长的不同，CCD、CFD和CD能允许较大的时间步长(7 ms).从算法的时间复杂度来看，CCD和CD算法的复杂度要大于CFD和ONAD，其中n表示网格点数.值得注意的是，在模型大小相同的时候，由于网格长度的差异，会造成四种算法的空间网格节点数n不同，也就造成了实际所需运算时间的不同，表中运算时间是在同一台计算机和相同编程环境下模拟得到的.综合模拟精度、空间网格长度、时间步长和占用内存四个因素来看，CCD方法进行声波方程数值模拟时，它具有高精度(时间四阶、空间六阶)、高计算效率(较大Courant数、占用内存少)和低数值频散(粗网格)的优势.需要说明的是，表3中四种方法各参数的选取时，均不考虑计算精度，空间网格和时间步长仅满足无数值频散和稳定性的要求.但在实际计算时，为了精度考虑，二者均需要按要求适当减小，这势必会进一步增加计算量和存储量.3 弹性波方程组合型紧致有限差分方法各向同性完全弹性介质中的二维波动方程可以表示为(43)其中，u(x,z)和w(x,z)表示x和z方向的质点位移，VP(x,z)和VS(x,z)表示纵横波波速.同声波方程一样，利用公式(8)和(9)，并代入波动方程，将u对时间的导数转换为u和w对空间的导数，可以得到时间四阶精度差分格式：(44)同理可以得到w的时间四阶精度差分格式为(45)弹性波差分格式与声波方程不同，公式(44)和(45)中除了位移u的空间二阶导数∂2ui,j/∂x2和∂2ui,j/∂z2外，还有二阶混合偏导数∂2ui,j/∂x∂z和其他四阶偏导数.假设偏导数∂2/∂x2，∂2/∂z2，∂/∂x和∂/∂z的差分算子为和则二阶混合偏导数∂2/∂x∂z的差分算子为四阶偏导数∂4/∂x4的差分算子表示为即对于二阶混合偏导数，可以在求出∂ui,j/∂x和∂ui,j/∂z之后，再次使用公式(19)，对另一个方向求导，从而求出空间二阶混合偏导数∂2ui,j/∂x∂z.同样，对差分公式中出现的u和w的其他高阶偏导数，也可以按照上述同样的方法进行迭代求取，这里不再赘述.4 模型试算为了验证方法的数值性能，我们采用组合型紧致方法对公式(43)和(44)表示的二维各向同性完全弹性介质的弹性波方程进行数值模拟.4.1 均匀介质模型模型长度和深度均为10 km，介质为泊松体，纵波速度3460 m·s-1，横波速度2000 m·s-1，在模型中心处激发10 Hz主频的雷克子波.由于最小速度为2000 m·s-1，按频散分析的结果，每个波长内需采样5.6个点，所以设置网格长度为35 m.同时，由于最大速度为3460 m·s-1，按稳定性条件，取时间步长为8 ms.长时程数值模拟是数值模拟中的一个重要内容(Chen,2006; Li et al.,2012; Gao et al.,2016).对于大尺度模型和长时程传播，即便是十分微弱的误差也会随着传播时间的增加而不断累积，最终可能导致不容忽视的数值假象.为了分析CCD方法的长时程数值模拟结果，数值模拟时没有使用吸收边界条件，图5(a—d)显示的是1 s、2 s、5 s和10 s时刻的位移水平分量波场快照，地震波场特征清晰，无可见的数值频散现象.为了比较，对该模型采用相同的参数，使用时间四阶空间六阶精度CFD格式进行数值模拟，图5(e—h)为相同时刻的波场快照.可以看出，在满足频散要求和稳定性条件的情况下，CCD方法具有计算效率较高、长时程稳定的特点，能适应较长时间的数值模拟，计算结果较为理想.在相同参数下，CFD结果存在很明显的频散现象，并且随时间的增加，频散更加明显，严重地影响了数值模拟的效率和精度.图5 均匀介质模型弹性波数值模拟波场快照(a) CCD(1 s); (b) CCD(2 s); (c)CCD(5 s); (d) CCD(10 s); (e) CFD(1 s); (f) CFD(2 s); (g) CFD(5 s); (h) CFD(10 s).Fig.5 Wave field snapshots using elastic wave equation in homogeneousmedium(a—d) are the snapshots at 1,2,5,10 s for the displacement component in X direction using CCD; (e—h) are the snapshots at 1,2,5,10 s for the displacement component in X direction using CFD.4.2 水平层状介质模型设置水平层状介质模型，共20层，每层厚度400 m，模型总厚度8 km，模型长度同样是8 km.第一层纵波速度为3500 m·s-1，往下各层纵波速度比上一层增加100 m·s-1.模型为泊松体，即纵横波波速比为1.732，所以模型中最小和最大波速为2020 m·s-1和5400 m·s-1.在地面中间处(X=4 km，Z=0 m)激发10 Hz雷克子波震源，人工边界采用PML边界条件(刘有山等，2012)处理.按前文分析结果，设置纵横向网格长度为36 m，时间步长为5 ms.图6显示的是地面接收到的地震记录，记录长度8 s，由于直达波能量明显强于反射波，所以图中对地震记录进行了增益处理.多层模型中地面模拟接收到的地震记录非常清晰，没有数值频散和不稳定数值结果，地震记录中的直达波、反射波和转换波显示清楚，说明组合型紧致有限差分算法可以有效地模拟弹性波在多层各向同性介质模型中的传播.为了检验本文算法长时程数值模拟的结果，图7给出在较粗网格条件下(网格长度36 m，时间步长5 ms)，CCD和CFD方法分别得到的检波点(X=4 km，Z=0 m)处的波形图，并与利用CFD方法在较细的网格(网格长度10 m，时间步长1 ms)情况下得到的近似理论解的参考波形进行对比.图7a显示三种情况下得到的反射波到达时间一致，从图7(b，c)的局部放大结果来看，在粗网格和长时程模拟条件下，CCD方法得到的地震记录平滑无数值频散，与之前的分析结果一致，且CCD方法和近似解析解的波形曲线基本重合，这说明了CCD方法的正确性.而CFD方法在相同条件下存在较严重的数值频散，且模拟时间越长，模拟结果与近似解偏差越大.这些结果和对比说明，CCD方法能够长时间计算稳定性，而且也能很好地压制数值频散.5 结论与建议有限差分法是地震波场数值模拟的重要方法之一，在目前常规有限差分数值模拟中为了获得更高的模拟精度，需要使用更多的节点，这样会增加计算量和所需的存储空间，还会增加人工边界处理的难度.本文从组合型紧致差分格式出发，建立了时间四阶、空间六阶的声波方程的差分格式，对该差分格式进行了模拟精度、频散曲线和稳定性分析，并与常用的几种格式进行了综合比较，该格式可同时计算空间一阶和二阶导数，且具有使用节点少(3个)，低数值频散(每个波长采样5.6个点)，较高稳定性(Courant数为0.792)和模拟精度(空间6阶)的特点，适用于大尺度和较长时程地震波场数值模拟.文中还建立了各向同性介质弹性波方程的差分格式，并进行了数值模拟，模拟结果显示地震波场特征清晰，说明了该方法的适用性.该方法不仅可以用于弹性介质的波场模拟，还可以进一步推广到二维或三维的各向异性介质、黏弹介质和双相介质等复杂介质的数值模拟中，为今后研究地震波传播规律、逆时偏移、全波形反演等工作提供了一种有效的方法.组合型紧致差分格式仅需使用三个节点，就能使一阶和二阶空间偏导数的差分精度达到传统差分方法的六阶精度，因而能够大幅度节约内存需求和计算量，对于三维大尺度模型的正演和反演等问题均具有重要意义.图6 水平层状介质模型弹性波方程数值模拟地面地震记录(a) X方向位移分量； (b) Z方向位移分量.Fig.6 Seismic records of the model′s surface using elastic wave equation in horizontal layers media(a) Displacement component in X direction； (b) Displacement component in Z direction.图7 接收点(4 km，0 km)处接收到的反射波水平位移波形图(a) 整道反射波记录对比； (b) 4000～4500 ms局部放大对比； (c) 6000～6500 ms局部放大对比.其中，红色加点实线是粗网格条件下的CCD数值结果，绿色实线是粗网格条件下。

井间地震正演模拟【摘要】井间地震以其高精度和高分辨率的显著特点在石油和天然气的勘探和开发中发挥着越来越重要的作用。

本文根据各向同性介质中的一阶速度—应力弹性波方程，构造出了此方程在交错网格中的高阶有限差分格式以及它的稳定性条件，研究了该方程的完全匹配层（PML）吸收边界条件，在此基础上实现井间地震正演模拟，得到了高精度的地震合成记录。

结果表明：井间地震可以观测到地震勘探中可能遇到的多种波场，利用数值模拟可以获得令人满意的合成地震记录。

【关键词】井间地震有限差分正演模拟完全匹配层1 引言随着海内外含油田勘探程度的日益增高，可能勘探到的新油田越来越少，勘探花费也在不断增加，寻找新油田变得越来越困难。

为了增加油气采收率，降低成本，勘探时对地震分辨率的要求也越来越高。

在这种情况下，凭借其高精度高分辨率的特点，井间地震技术越来越受到油气勘探开发行业的重视。

井间地震是将震源与检波器放置在相邻的两口井中，在目的层内部进行地震波的激发与接收。

井间地震资料包含丰富的波场信息，可以提供超高频率的纵、横波资料，将这些资料与其他资料综合研究，可以解决薄互层、储层连通性、流体分布、注气效果、压裂效果等复杂的油藏地质问题。

2 各向同性介质弹性波方程数值模拟各向同性介质中质点速度和应力表示的一阶弹性波方程为：（式3）其余参量的差分格式同理可得。

通过波数域特征值分析，可以得出三维弹性波方程交错网格高阶有限差分数值解法的稳定性条件为：5 数值模拟算例5.1 理论模型模型长高各为2000m分A，B，C三层。

A层深800m，横波速度为1500m/s，纵波速度为2500m/s，密度为2000kg/㎡；B层深1000m㎡，横波速度为2000m/s，纵波速度为4000m/s，密度为2500kg/㎡；C层深200m，横波速度为3000m/s，纵波速度为5000m/s，密度为3000kg/㎡。

网格步长△x=△z=10m，△t=1ms，记录长度为1000ms，。

地震波传播模拟中的数值方法一、引言对地球上发生的自然灾害进行研究和预测一直是人类所探究的课题之一。

其中，地震是一种造成极大灾害的自然现象，它的预测和探测对减轻地震对社会影响，提高人类对灾害的应对能力，具有重要意义。

地震波传播模拟是地震研究领域的重要课题，为了更好地预测地震和应对地震灾害，需要对地震波传播的数值模拟方法进行深入研究。

二、地震波传播数值模拟的方法1. 有限差分法（FDTD）有限差分法，英文全称为Finite Difference Time Domain，是一种常用的求解电磁场和声场传播问题的数值方法。

FDTD方法利用有限差分逼近微分算符，将偏微分方程离散化，然后通过差分方程组求解离散化问题。

FDTD方法的优点是较为简便和直观，对于一些基础场问题可以精确求解，但是FDTD方法在离散化问题域时会导致误差，对于具有复杂形状、边界不规则和含有多个介质的问题，其求解需要繁琐的预处理工作和较为复杂的网格划分，求解过程也较为复杂。

2. 有限元法（FEM）有限元法，英文全称为Finite Element Method，是一种广泛应用于工程和科学计算领域的数值方法。

它是通过将一个复杂的问题域分解成多个小问题域，用简单的数学公式在每个小问题域内求解，通过对这些小问题域的求解累加得到整个问题域的解。

FEM方法的特点是能够对不规则的计算域进行处理，求解过程较为直观和简单，对于多介质、弹性、非线性等问题也有很好的处理能力。

但FEM方法对于较为复杂的问题各向异性和自由面的处理比较困难。

3. 间接边界积分法（BEM）边界积分法，英文全称为Boundary Element Method，是近年来发展起来的一种求解偏微分方程的数值方法。

BEM方法将待求解的域分为界面和域外两部分，通过界面上的边界积分求解内部问题。

BEM方法对于不规则和异形问题的边界条件求解有很好的处理能力，并且具有较高的精度和较低的计算量。

但是对于非线性问题处理不够准确，对纯内部问题的求解效果不如其他方法。

Kelvin-Voigt黏弹性介质地震波场数值模拟与衰减特征严红勇;刘洋【摘要】利用高阶交错网格有限差分模拟Kelvin-Voigt黏弹性介质中传播的地震波,同时将完全匹配层吸收边界条件引入到其边界处理中.数值模拟结果表明,完全匹配层吸收边界效果好,高阶有限差分能模拟得到的黏弹性介质波场精度较高.对模拟的黏弹性波场进行分析,表明介质的粘滞性使地震反射波的能量变弱,高频衰减明显,并比低频衰减得快,主频向低频方向移动,有效频带变窄,即降低了地震波的分辨率；并且反射转换波比反射纵波要衰减得快；而且还随着传播距离的增加,其峰值频率也逐渐降低.通过数值模拟分析具有不同的粘滞系数介质对地震波的吸收和衰减,结果表明随着粘滞系数的增大,地下介质对地震波的吸收衰减更明显.%This paper uses finite difference algorithm of high-order staggered-grid simulate Kelvin-Voigt viscoelaslic media of seismic waves and meanwhile introduces the perfectly matched layer(PML) absorbing boundary condition into its boundary. Numerical simulation demonstrates that the effect of this algorithm of absorbing boundary is very good and the wavefield of viscoelastic media obtained from high-order finite difference is relatively accurate. An analysis of viscoelastic wavefield simulation shows that the energy of the reflected wave becomes weaker,the attenuation of the high frequency wave is much more apparent in comparison with that of the low frequency wave,the main frequency becomes closer to the low frequency,and the effective bandwidth is narrower,which all induce low resolution of seismic wave according to the simulation of viscoelastic wavefields. Besides,the attenuation of PS-wave is much more rapid incomparison with that of PP-wave and the peak frequency becomes lower gradually with the increase of the propagating distance. It is also shown that the absorption and attenuation are more apparent with the increaseof viscosity coefficient by analysis of the absorption and attenuation of seismic wave in different viscosity coefficient media.【期刊名称】《物探与化探》【年(卷),期】2012(036)005【总页数】7页(P806-812)【关键词】黏弹性;交错网格;有限差分;完全匹配层;衰减【作者】严红勇;刘洋【作者单位】中国石油大学油气资源与探测国家重点实验室,北京102249;中国石油大学CNPC物探重点实验室,北京102249;中国石油大学油气资源与探测国家重点实验室,北京102249;中国石油大学CNPC物探重点实验室,北京102249【正文语种】中文【中图分类】P631.4目前，在讨论地震波传播理论时，绝大部分情况下是把地震波看作弹性波。

FDTD方法中的吸收边界条件有几种吸收边界条件可以用于时域有限差分时间-域（FDTD）方法中，以模拟开放空间的无穷边界。

这些吸收边界条件的目标是吸收或反射尽可能少的能量，并且尽量减小边界对计算域中场分布的影响。

第一种常见的吸收边界条件是完美匹配层（Perfectly Matched Layer，简称PML）。

PML是一种人工吸收层，模拟了在模拟区域边界的理想吸收层。

PML通过嵌入在实际计算区域之外的等效吸收层来实现，使得电磁波在边界处被吸收。

PML通过引入Lossy介质来吸收能量，其吸收效果由材料的吸收率、传播常数和吸收边界厚度等因素决定。

PML在FDTD方法中已经被广泛应用，并且在模拟计算区域边界时表现出了较好的吸收效果。

第二种常见的吸收边界条件是吸收边界条件（Absorbing Boundary Condition，简称ABC）。

ABC是一类在边界处引入的边界条件，通过使边界处的反射系数减少来实现。

常见的ABC包括Mur ABC、Berenger ABC和磁流体（Magnetic Fluid）吸收边界等。

ABC通常需要通过增加节点的数目或周围嵌入吸收材料来实现。

这些方法通常需要对计算域进行扩展或使用专门的吸收元件。

第三种常见的吸收边界条件是开放边界条件（Open Boundary Condition，简称OBC）。

OBC试图模拟一个无穷大的开放空间，以模拟场的自由传播。

OBC通常在计算域边界上应用自然边界条件或辐射条件。

自然边界条件要求边界上电场、磁场的法向分量为零，而切向分量与场的梯度成正比。

辐射条件则要求边界上场的入射和反射能量相对较小。

OBC通常在计算域的边界上要求场的梯度不受限制，并且在边界处不反射能量。

这些吸收边界条件在FDTD方法中有各自的优缺点，适用于不同的应用和场景。

PML是一种广泛使用的有效吸收边界条件，但需要较复杂的计算和实现。

ABC通过减小边界反射系数来实现吸收，可以有效减小计算域边界对场的影响，但需要对计算域进行扩展或使用专门的吸收元件。

二维地震波场的五点八阶超紧致有限差分数值模拟周诚尧;汪勇;蔡伟祥;桂志先【摘要】首先将迎风机制引入五点八阶超紧致有限差分(CCD8)格式,得到五点七阶迎风超紧致(UCCD7)格式,并对两种格式进行数值频散分析和精度分析;其次建立了二阶声波方程的位移场时间四阶离散格式,将五点CCD8格式和五点UCCD7格式分别应用于位移场空间导数的求取,并分析这两种格式的稳定性条件;最后基于PML 边界条件,将上述两种格式分别应用于声波方程的均匀介质、水平层状介质及Marmousi模型的数值模拟和波场特征分析及对比.研究结果表明:相较于普通紧致差分,五点CCD8格式具有小截断误差、高模拟精度、低数值频散、高稳定性、所需网格点数少的优点;引入迎风机制后,声波方程的五点UCCD7格式稳定性得到进一步提高;模型试算的结果验证了五点CCD8格式适用于复杂介质的数值模拟,模拟精度和计算效率都高.【期刊名称】《石油物探》【年(卷),期】2019(058)002【总页数】11页(P176-186)【关键词】五点八阶超紧致差分;五点七阶迎风超紧致差分;数值模拟;数值频散;稳定性条件【作者】周诚尧;汪勇;蔡伟祥;桂志先【作者单位】油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100;油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100;中石化重庆涪陵页岩气勘探开发有限公司,重庆408014;油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100【正文语种】中文【中图分类】P631地震正演模拟是探索地震波在介质中传播规律的重要手段,基于波动方程的数值模拟方法是一种地震正演模拟方法[1]。

目前波动方程法主要包括伪谱法[2]、有限元法[3]、边界元法[4]、谱元法[5]以及有限差分法。

有限元法和有限差分方法相比,二者精度相同时,有限元法耗时长[6]。

有限差分法因其具有简单灵活、计算效率高以及占用内存小等优势被广泛应用于地震波场数值模拟[7]。

第42卷第2期2003年6月石　油　物　探GE OPHY SIC A L PROSPECTI NG FOR PETRO LE UMV ol.42,N o.2Jun.,2003文章编号:100021441(2003)022*******地震波场数值模拟方法张永刚(中国石油化工股份有限公司科技发展部,北京100029)摘要:简要总结了地震波场数值模拟的各种方法的基本原理及其主要特点,对最近在该领域出现的一些方法和研究结果做了简要的阐述,并对比了各种方法的优缺点。

在此基础上提出了运用波动方程数值模拟作为基础,结合射线方法辅助识别波场类型,用于分析异常波的产生机理和出现特点的基本思想,这对复杂条件下的地震勘探具有指导和借鉴意义。

关键词:地震波场;数值模拟;射线追踪;有限元;伪谱法;正演模拟中图分类号:P63114+1 文献标识码:AOn numerical simulations of seismic w avefieldZhang Y onggang(Department of Science and T echnology Development,SI NOPEC,Beijing100029,China)Abstract:This paper reviews the principles and characteristics of various numerical simulations of seismic wavefield,and com2 pares the merits and defects of the simulations.S ome newly emerged methods and results are briefly discussed.The author pro2 poses to study the generation mechanism and characteristics of abnormal waves based on wave equation numerical simulation supplemented by ray tracing.K ey w ords:seismic wavefield;numerical simulation;ray tracing;finite element;pseudo2spectrum;forward m odeling 地震波场数值模拟是研究复杂地区地震资料采集、处理和解释的有效辅助手段,地震波场数值模拟的主要方法包括2大类,即波动方程法和几何射线法。

地震波交错网格高阶差分数值模拟研究摘要: 地震波数值模拟技术是勘探地球物理学中的重要组成部分，研究通过弹性波一阶速度——应力方程，采用交错网格高阶有限差分法实现了地震波在各向同性介质中的高精度的数值模拟，并采用完全匹配层( PML) 吸收边界来消除边界反射，可取得较好的效果。

通过模型的正演计算和复杂模型的处理结果表明，交错网格高阶有限差分法数值模拟是一种快速有效的地震波数值模拟方法。

关键词: 地震勘探; 交错网格; 有限差分; 数值模拟引言地震数值模拟是模拟地震波在介质中传播的一种数值模拟技术，随着地震波理论在天然地震和地震勘探中的应用，地震模拟技术便应运而生，并随着地震波理论和计算机技术的发展，地震数值模拟技术自20世纪60年代以来也得到了飞速发展，形成了目前具有有限差分法、有限元法、虚谱法和积分方程法等各种数值模拟方法的现代地震数值模拟技术。

有限差分法是偏微分方程的主要数值解法之一。

在各种地震数值模拟方法中，最早出现的数值模拟方法是有限差分法。

Ａlterman和Karal(1968)首先将有限差分法应用于层状介质弹性波传播的数值模拟中。

此后，Boore(1972)又将有限差分法用于非均匀介质地震波传播的模拟。

Alford等（1974）研究了声波方程有限差分法模拟的精确性。

Kelly等（1976）研究了用有限差分法制作人工合成地震记录的方法。

Virieux(1986)提出了应用速度——应力一阶方程交错网格有限差分法模拟P——SV波在非均匀介质中的传播。

交错网格方法提高了地震模拟的精度和稳定性，并消除了部分假想。

有限元法也是偏微分方程的数值解法之一。

Lysmer和Drake（1972）最早将有限元法应用于地震数值模拟。

Marfurt(1984)研究对比了模拟弹性波传播的有限差分法和有限元法的精度。

Seron等（1990,1996）给出了弹性波传播有限元模拟方法。

Padovani等（1994）研究了地震波模拟的低阶和高阶有限元法。