任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟

基于交错网格有限差分弹性波正演模拟及波场特征分析【摘要】为研究和认识多种储层中弹性波的波场特征，以利于多波地震资料解释，高精度数值模拟是有效的方法之一。

本文在弹性波方程基础上，采用高阶交错网格有限差分技术模拟地震波在各向同性介质和各向异性介质中的传播，可得到不同类型介质的弹性波场。

同时，文中也分析了各向异性系数对多波波场特征的影响。

通过对高精度数值模拟得到的波场快照对比研究表明，该方法可有效获得高精度弹性波正演结果，为研究各种复杂介质中弹性波的波场特征和传播规律奠定了基础。

【关键词】多波多分量波场特征各向异性弹性波正演1 引言随着油气田勘探技术的不断发展[1][2]，人们对地震资料的认识也不断加深，纵波地震资料在含油气的显示上存在一些不确定性，单一纵波资料解释的多解性问题尤为突出。

在地震勘探领域中，过去一直把各向同性弹性体理论作为研究地下介质的前提，但是在实际地层中普遍存在各向异性，地下介质的各向异性（如周期薄互层引起的各向异性、以及裂隙引起的各向异性）产生的弹性波场与各向同性介质产生的弹性波场存在着不可忽略的差异。

由此，多波地震勘探作为油储地球物理的主要方法之一应运而生。

在多波资料解释过程中，要求搞清楚储层的岩性与多波的波场特征之间的关系，因此，多波波场数值模拟技术显得非常重要。

高精度数值模拟技术是联系地震、地质、测井以及油藏工程的纽带，其作用主要体现在提高人们对各种复杂介质中地震波传播规律的认知，并可为新技术、新方法提供试验数据，以满足方法技术研究的需要，同时也可以检验解释结果的正确性。

弹性波波动方程高精度数值模拟可以得到全波场信息，包含了地震波的动力学和运动学特点，为准确描述地震波场特征和波的传播规律奠定基础，本文在弹性波方程基础上，采用高阶交错网格有限差分技术模拟地震波在各向同性介质和各向异性介质中的传播，比较地震波在各向同性介质和各向异性介质中的波场响应异同，并分析了各向异性系数对多波波场特征的影响，这对研究各种复杂介质中弹性波的波场特征和传播规律有着重要的意义。

地震波场的高阶交错网格有限差分模拟霍凤斌;李振鹏;徐发;张涛【摘要】This paper analyzes the stability and convergence of the seismic wave ifeld by using the high-order staggered-grid limited differential method of joining the absorbing boundary condition and attenuating zone to simulate the elastic wave equation. The results of the simulation of both isotropic-and anisotropic-medium models show that the grid frequency dispersion of the high-order differential wave equation simulation is smaller and more accurate. Therefore, this method should improve the efifciency of seismic prospecting and of the associated data interpretation.%应用高阶交错网格有限差分算法，并加入吸收边界条件和衰减带，对弹性波方程进行模拟，分析了其稳定性和收敛性。

通过对各向同性和各向异性介质模型的模拟表明，高阶差分波动方程模拟的网格频散较小、精度较高、效果较好，可为地震勘探及其资料解译提供技术手段。

【期刊名称】《上海国土资源》【年(卷),期】2014(000)001【总页数】4页(P97-100)【关键词】地震波场;波动方程;有限差分;边界条件;交错网格【作者】霍凤斌;李振鹏;徐发;张涛【作者单位】中海石油中国有限公司上海分公司，上海200030;中海石油中国有限公司上海分公司，上海200030;中海石油中国有限公司上海分公司，上海200030;中海石油中国有限公司上海分公司，上海200030【正文语种】中文【中图分类】P315.01随着地震波动理论在天然地震和油气地震中的应用，以及计算机技术的飞速发展，在现代地震数值模拟领域逐渐形成了有限差分法、有限元法、虚谱法和积分方程法等求解波动方程的方法。

二阶弹性波动方程高精度交错网格波场分离数值模拟陈可洋;杨微;刘洪林;吴清岭【摘要】给出了一种等价的二阶弹性波动方程,以解决弹性波场中完全弹性波动方程不能完全分离耦合的纵、横波波场问题.应用高阶交错网格有限差分法求解该波动方程,并使用通量校正技术(FCT)进一步压制频散,采用均匀介质模型和层状介质模型进行波场分离数值试验,精确得到了混合波场、完全分离的纯纵波及纯横波波场.数值结果分析表明,本文方法在均匀介质情况下准确可靠,在分离后的纯纵、横波波场中可观察到较为丰富的能量转换信息,这对认识复杂弹性波的传播规律及弹性波理论具有重要意义.【期刊名称】《物探与化探》【年(卷),期】2009(033)006【总页数】4页(P700-703)【关键词】地震波场分离;高阶交错网格;等价二阶弹性波动方程;数值模拟;通量校正技术【作者】陈可洋;杨微;刘洪林;吴清岭【作者单位】大庆石油学院,地球科学学院,黑龙江,大庆,163318;大庆石油学院,地球科学学院,黑龙江,大庆,163318;大庆石油学院,地球科学学院,黑龙江,大庆,163318;大庆油田有限责任公司,勘探开发研究院,黑龙江,大庆,163712【正文语种】中文【中图分类】P631.4地震波场数值模拟技术一直是勘探地球物理领域内最为活跃的研究内容之一。

利用精确的波动方程数值解模拟复杂地下波场，为研究地震波传播机理、地震资料的特殊处理方法以及复杂地层的解释等许多方面提供更为科学的数学物理依据［1－3］。

在实际勘探中仍然存在许多复杂问题值得深入探讨，这对研究波动理论，指导地震资料的特殊处理与复杂构造解释具有重要意义。

例如，多波多分量地震记录的每一个分量均包含不同的波型，简单地把垂直分量看成P波，水平分量看成S波是很不合理的。

常规地震资料处理总是希望处理单一分量的标量波场，因此从混合波场中解析出纯P波和纯S波场是进行速度分析、偏移成像等地震资料数据处理的前提。

马德堂等［4］提出满足P波为无旋场、S波为无散场的等价方程思路，采用虚谱法来实现波场分离数值模拟［5－6］，得到较好的数值模拟结果，但是虚谱法很难处理吸收边界条件和自由表面边界条件，同时计算时间和存储量的代价很大，计算网格节点数必须满足2N，且可能存在着Gibbs效应。

基于有限差分法的地震波数值模拟研究综述姚铭;高刚;周游;蔡伟祥;周永娇【摘要】针对有限差分算法在进行地震波数值模拟过程中的模拟精度问题,在进行广泛文献调研后,分析了国内外关于提高有限差分法模拟精度的进展研究,主要包括提高有限差分的阶数及改进离散方式2方面。

论述了有限差分法与其他数值模拟方法的优缺点,指出了有限差分法具有占用内存小、计算量小且易于实现等诸多优点,并总结了有限差分法在进行数值模拟过程中存在的主要问题,包括震源函数、边界条件、频散问题以及稳定性分析。

同时指出了有限差分法在地震波场模拟应用中的发展趋势。

【期刊名称】《能源与环保》【年(卷),期】2017(039)010【总页数】6页(P75-79,85)【关键词】有限差分;数值模拟;震源函数;边界条件;频散问题;稳定性分析【作者】姚铭;高刚;周游;蔡伟祥;周永娇【作者单位】[1]长江大学油气资源与勘探技术教育部重点实验室,湖北武汉430100;[2]长江大学地球物理与石油资源学院,湖北武汉430100;[3]中石化重庆涪陵页岩气勘探开发有限公司,重庆408014【正文语种】中文【中图分类】P631.4地震波数值模拟技术作为波动方程正演、逆时偏移和全波形反演等问题的关键，一直都是理论研究的重要内容［1］。

其理论基础主要包括波动理论和射线理论［2］。

其中射线理论的基本原理是直接求出质点运动方程的近似解，这样虽然能够大大提高计算的速度，但其精度较低且缺少波动力学信息。

而波动理论则是通过数值计算方法求出波动方程的解，因此其模拟结果中包含了丰富的波动信息，便于了解地下地质情况，模拟结果也较精确。

基于波动理论的数值模拟方法主要包括边界元法、有限元法、有限差分法以及伪谱法等，这些方法在数值模拟过程中有其各自的优点，但问题也较突出［3-6］。

如边界元法虽然易于处理边界问题，但其计算量较大;有限元法虽然网格剖分灵活，且能够精确模拟复杂地质形态，但存在占用内存较大的问题，无法普及;伪谱法能够精确模拟地震波的传播规律，但其计算速度较慢，在进行并行计算时也会出现诸多问题，因此适用性不广。

第23卷 第1期地 球 物 理 学 进 展V ol.23 N o.12008年2月(页码:40~48)PRO GR ESS IN GEO PH YSICSF eb. 2008地震波有限差分数值模拟框架下的起伏地表处理方法综述王雪秋1, 孙建国2(1.吉林大学地球探测科学与技术学院,长春130026; 2.国土资源部应用地球物理综合解释理论开放实验室,长春130026)摘 要 本文查阅了近30年来国内外关于有限差分方法模拟近地表复杂介质方面的文献,分析了在有限差分法框架下复杂地表地形起伏问题的处理及实现方法.总结了地震波有限差分数值模拟方面的成果.复杂地表类型很多,对地震资料的影响程度也不尽相同,要想完全描述近地表影响因素是不可能的,而关键的问题就是自由边界条件方程的处理.所有处理及实现方法中不外两种:直接法(局部)和间接(全局)法.文中列出的方法,各有特点.我们的目的是寻找合适的方法解决适合的问题.关键词 有限差分,自由边界,直接法,间接法,映射,地表几何形状中图分类号 P631 文献标识码 A 文章编号 1004 2903(2008)02 0040 09The state of the art in numerical modeling includingsurface topography with finite difference methodWANG Xue qiu 1, SUN Jian guo 2(1.Colleg e of Ge o ex p loration sc ienc e an d te chnology ,J ilin Univ er sity ,Ch an gch un 130026,China;2.L abor atory of w av e theory and imag ing tech nology ope n r esearc h labor atory f or integr ate d geop hysicalinter pr eta tion theory of the M inistr y of lan d resourc es.Chang chun 130026,China)Abstract A s a too l for under standing the reg ulatio n f wave pro pag ation in complex media and ex amining the effects o f different kinds of metho ds,seismic numerical modeling plays an impo rtant ro le in modern g eophysical prospect ing.T his pa per briefly intro duced sev eral mo deling metho ds including surface to po gr aphy w ith finite difference in r ecent 30years.T w o implementatio n methods are proposed:der ect met ho d w hich includes r otated coor dinate method,v acuum metho d and general imaging method and indirect method.O ur aim is to find the suit able metho d to so lve the ada pted quest ions.Keywords finit e differ ence,fr ee sur face,direct method,indirect method,mapping surface topog raphy收稿日期 2007 03 10; 修回日期 2007 06 20.基金项目 国家自然科学基金项目(40574052)和教育部骨干教师资助计划项目联合资助.作者简介 王雪秋,女,1974年生,吉林大学地球探测科学与技术学院讲师,在读博士,主要从事地震波理论及波场正演模拟等方面的研究.(E mail:w angxueqiu@)0 引 言经典反射地震勘探的理论基础有3个假设前提:地表是水平的;地下介质是水平层状的;地震波是近垂直反射的.这种假设简化了理论模型和参数,使得到的地震数据便于解释,并且诸如CDP 技术,NMO 和DM O 以及偏移成像技术等都是以这些假设前提为基础建立起来的,这些技术事实上也在很多近似符合条件的地区得到了非常好的结果.但是当在勘探领域向复杂地区深入,勘探技术不断进步,大尺度长排列的应用等情况下前面的假设前提失去了意义,尤其是在我国西部地区复杂的地质条件进行的反射地震勘探,实际情况不允许这些假设条件的存在[1~4].另外,在研究地球介质的性质,地震灾害性质,金属矿勘探及环境工程勘察等方面,更加促使我们抛弃这些假设前提,而注重于更加实际的波场的精细刻画.因此研究复杂地区的波场特征及传播规律成为发展的必然.尤其是复杂地表及近地表1期王雪秋,等:地震波有限差分数值模拟框架下的起伏地表处理方法综述情况对波场的影响规律的研究更是具有重要意义.这是由于我们多数时候都是在地表附近接收到的地震数据,所以地表的形状及介质物理特性对波场的改造相对于一些深层的信息来说是剧烈的.而关于近地表的地震波场 散射理论则成为研究的热点[5~10].一般情况下,静校正可以消除这种影响,但是多数情况下尤其是近地表地质地形非常复杂的地区,静校正的作用是有限的[11~16],而地震波的数值模拟则是完成这种认识的一种直接有效的手段.对复杂地表问题的地震波 散射场的解,一般有两种形式;数值解和解析解.解析解能够解决的地表模型一般比较简单[16~26],所以对复杂的地表情况,采用数值解是必需的.复杂地表一般指两种情况:复杂的几何形状,例如山地、隆起、河谷、凹陷及高差剧烈,地形多变等;近地表介质复杂的物理性质,例如弹性、粘弹性、孔隙性、非均匀性和随机性等.前一种情况就是自由边界条件,也就是来自地下的应力、应变或速度消失的界面.而在一般给出的自由边界条件方程中,都是给出沿界面的法向分量和垂向分量为零的方程.在水平界面情况下,坐标轴与界面重合,这个方程很好处理.但是,当界面与与坐标轴不重合的时候,也就是地表有几何形状起伏时,这个方程的处理就很困难.后一种情况实际上就是介质模型的确定,也即方程的确定.根据地区的地质状况可以采用很多种介质模型的近似,例如弹性波方程、粘弹性介质方程、非线性方程等.而在任何一种方程中都存在自由边界问题,可以说,介质方程的确定是自由边界条件存在的载体.另外一个需要特别指出的问题是,在经典的地震勘探理论中,由于近垂直反射假设条件的前提,人们一般选用声学介质模型来近似地下介质,由于参数简单,易于处理,许多成熟的,在工业界广泛使用的数据处理技术都是以声波模型为基础的.但是声波方程是在流体假设条件推导出来的,不存在起伏形状自由边界问题,不符合实际近地表地震波场研究模型条件[14].所以复杂地表问题最核心的问题就是自由地表几何形状处理问题.关于地表几何形状对波场产生的影响问题,早在上个世纪40年代,M.B.Widess[28]就提出含地形的自由地表响应是影响地震填图的主要因素之一,并且是地震构造解释的主要误差之一.只是由于当时的计算技术落后,无法开展这项工作.20世纪60年代后期,随着计算技术软硬件的发展,地震波的数值模拟技术也迅猛发展起来[29~35].在所有的模拟方法中网格法一直由于其对地形和介质的灵活划分而占主导地位,因此也成为几何状地形起伏地表地震波数值模拟的首选方法.在进行近地表地形模型的网格剖分时,最重要的是考虑如何对地表进行有益而准确的网格剖分,这种剖分既要符合所采用的方法原理,也要符合波动传播的物理机制,最好不要产生不必要的由计算方法而引起的波动现象.在网格法中,针对任意起伏的地表使用有限元方法来模拟近地表的波的传播,能够灵活处理边界的几何形状问题,但是计算比较慢,较其他方法的成本要高.另外一个成本不高,又能够处理自由地表几何形状的方法是边界元方法.已经有人将这个方法应用在了弹性介质的波场模拟.对于内部结构相对简单的介质(均匀或线性变化),这个方法在节省计算内存时相对于其他方法占有优势,原因在于计算是沿着边界来完成的.边界可以是复杂的形状,只要曲线闭合即可.但是仅适用于地表几何形状的模拟,对于介质内部的复杂的岩性和构造则很难处理.而一些谱方法则要求自由边界函数具有光滑的特性,并且会出现一些譬如折叠效应、吉布斯效应等不期望出现的效果[36~44].在网格法中,有限差分方法是较早开展起来的模拟方法,也是一种比较简单明确的模拟方法.在上个世纪的60年代美国的Alterman及其研究小组就应用这种方法来解层状介质的弹性波方程.同时Claerbo ut教授在研究地震成像中也提出了有限差分模拟方法,并给出了公开的算法[4].具有里程碑意义的是R.M.Alfo rd等人对有限差分模拟精度问题的探讨,在他们的研究成果中详细的讨论了网格点取值与差分算子的关系,同时还使用了高阶差分(四阶)[33].到了上个世纪70年代,K.R.Kelly等人给出了2维弹性波方程的完整的有限差分格式,并且给出了当时比较复杂的模型(分块的层,断层等)的时间切片和模拟记录图[32].进入80年代,又发展出了一种时间和空间上可达任意阶的高阶交错网格有限差分法[29,45~46],这种方法的优势是[46]:(1)震源以粒子振动速度或应力的形式表达,可以直接代入方程;(2)能够精确的处理自由边界条件;(3)差分算子是局部的,计算模型不必驻留内存,节省内存空间,节约成本,就算速度快;(4)方便使用高阶算子,提高计算精度;(5)易于与其他方法联合使用;(6)算法容易实现.因此越来越受到重视.而研究者一直也在探讨其在含几何形状起伏的自由地表条件的地震波场数值模拟中的应用.41地 球 物 理 学 进 展23卷在早期的常规有限差分地震波数值模拟中,一般要求自由地表是水平的或者是45度倾斜界面.但是有关含地表形状的自由边界处理方法也出现的很早.MIT的Boore在其博士论文中除了应用有限差分方法来解弹性波方程外,还首次应用有限差分方法研究了倾斜界面中LOVE波的传播机制[46].Ilan等[47~49]就用有限差分方法模拟了P SV波在起伏自由边界弹性介质中的传播,但是在线性分段之间的过度点没有进行处理.并且这种方法要求用不规则网格引起一些不必要的复杂性和在迭代过程中计算精度不够问题.1981年,Boore等应用Ilan的方法模拟了垂直入射的体波在45度倾斜界面的碰撞反射,但是没有处理角点.这些处理手段都是由低阶有限差分方法来实现的,多数情况下受差分阶数及显式实现算法的影响,模拟结果的精度和收敛性都不是很高[50].而现在在使用有限差分方法来进行近地表数值模拟时,人们多数采用高阶交错网格来离散方程[46,50~67].在所有关于含起伏地表几何形状自由边界条件的地震波数值模拟方法中,我们可以依据方法对自由边界的处理方法将其分为两类.一类是直接处理自由边界条件,而与所用介质的波动方程无关,称之为局部法,也可以叫做直接法[46,50~51,54~57,59~62].另一类是对介质模型做变换,即将实际上曲线(曲面)坐标系下的模型方程,变换到规则网格的坐标系下,同时自由边界条件方程也如此变换,我们称这种方式为全局法,也可以叫做间接法[11~14,53~53,64].1 自由边界条件自由边界是实际上就是固体地球与空气接触的表面,这是一个地球介质的不连续界面.这个界面也是来自介质内部的应力应变消失的界面,同时这个界面也强烈地影响着下面的固体介质内的波场的一些物理现象,因为地表的任何起伏都能引起在介质中传播的波场参数的剧烈变化.在数学表达上,自由边界一般都是给出的都是按照牛顿力学原理的应力或位移微分方程,这个边界方程一般是一阶的[66],例如:U z+ Wx=01-2 22Ux+Wz=0,(1)这是一个2D的自由边界条件,式中U,W是水平和垂直方向的位移, , 分别是纵横波速度.这个自由边界条件方程与应力为零是等价的.3D的边界方程与其类似.对这个方程应用有限差分方法时,由于在边界上要降低差分网格的阶数,往往造成了计算结果精度不够或计算结果的不稳定,尤其是在近地表介质物理性质横向剧烈变化或地形起伏剧烈的情况下.因此,在地震波数值模拟中自由边界是对波场影响最大的边界.同时这种影响也能够以最小误差来描述,因为其它层位的参数和位置的不确定因素要远大于自由边界,而自由边界是不需要定位的,它是事实上存在的.因此一旦关于自由地表的准确描述的数值模拟技术可以应用,其前景将十分可观.一般情况下,自由边界的数值模拟(包括地球的不均匀性)产生的影响都能自动在理论合成的结果中提供,并且所提供结果比不加自由边界的更加真实.因此在有限差分法地震波数值模拟中自由边界的处理技巧不但涉及到计算的精度,还涉及到方法的稳定及结果的准确性[17~18,29,32~33,47~49].在有限差分法的常规模拟中,自由边界条件有三种处理方式,一种是自由边界节点位于差分网格线上,这是一种常用的方法[32~33,47~49],如图1所示.另一种是自由边界节点位于差分网格线之间[67],见图2所示.节点位于网格上,便于实现,尤其是在显式差分实现方式中;但是稳定性差,精度不高,因为一般只能采用低阶差分格式.而节点在网格之间可以采用隐式差分来实现,是无条件稳定的.前面两种方法都可以在起伏地表处理中使用,而第三种则是由Levander[29,68]提出的 虚象法,即假设应力分量关于z=0对称,并向上沿拓两个网格,使其自然满足自由边界条件,类似于电磁场中的虚象原理,但是仅能应用在水平的自由地表条件中.421期王雪秋,等:地震波有限差分数值模拟框架下的起伏地表处理方法综述2 起伏地表直接处理方法起伏地表的直接处理方法也称为局部处理方法.能够直接处理起伏地表自由边界条件而不必考虑介质方程的方法主要有三种:局部旋转坐标法,真空法,广义虚象法.2.1 局部旋转坐标法局部旋转坐标法是一种出现比较早的处理自由边界起伏的方法,1975年,Ilan [47~49]等用有限差分方法模拟P SV 波在起伏自由边界弹性介质中的传播时就采用这种方法来实现自由边界条件的处理,1988年,Jih [50]等人发表的成果将这个方法进一步完善.其原理如图3所示,把不规则的地表地形分段线性化,也就是将自由边界按不同的斜率划分成线性段.然后使用局部旋转坐标系的方式(平行于分段的倾斜界面)对应力边界条件进行插值.在两个线性段的衔接处在平分拐角的基础之上使用一阶近似的边界条件.图3 分段起伏地表示意F ig.3Cartesian coo rdinate system f or the 2D finitediffer ence scheme and miscellaneous co rner points on the t opog raphy对于给定如图3这样曲折的地表,进行数值模拟需要进行一些特殊的处理.目前的处理方式是利用边界条件(方程2.1)处理图中所示的角点.这种问题的处理早在1970年的Alterman 和Lo ew enthal 就提出了两种计算90度和270度角点处位移的方法:第一种处理方式是设此时水平和垂直正应力为零,并且有U x = W z =0;第二种处理方式是对角点进行轻微的平滑,使之保持在45度以内,见图4.图4 倾斜界面上的坐标旋转Fig.4 Cartesian coo rdinat e sy st em rotat edfo r an inclined sur face这里 z 是固定的,而 x 则是个变量,他们之间具有这样的关系 x = z cot !,此处的!是X 坐标与自由界面的夹角.设!!45∀,图3中给了A 到F 六个不同的角点,因此有六种计算方式来处理他们.首先要给出随着x 增加而增加的方程,因为要计算角点上的值,就要计算其临近点上的位移值,为了求边界上的值,必须在每一个时间步长上给定其计算顺序,首先用离散的边界差分方程来计算其水平和垂直方向上自由边界段,而后是圈闭点数,最后才计算角点值.其几何形状控制在在每一段上至少要包含3个点,这样每一个角点都被圈闭点分开,空间步长的选取要根据介质的速度、频率等确定.在这个倾斜的自由界面上边界条件方程离散如下:U #i,k =U #0-sin2!(W #i+1,k -W #i,k+1)W #i,k =W #0-sin2!1-222(U #i +1,k -U #i,k+1),(2)这里U #0=∀0U i+1,k+1+(1-∀0)U i,k+2,W #0=∀0W i+1,k +1+(1-∀0)W i,k+2,∀0=2sin 2!.(3)当!∃45∀时,边界方程成为如下格式:U 0=∀0U i+1,k+1+(1-∀0)U i+2,k ,W 0=∀0W i+1,k+1+(1-∀0)W i+2,k ,∀0=2cos 2!.(4)这个只是比较理想的边界处理方式,具体到图3中这样的边界,及到六个节点在自由界面的情况,每一种情况都有相应的节点边界和坐标的处理方式,当然,这与节点的位置及角度有关,但是处理过程与上述类似.2.2 真空法这是一种处理起来比较容易的方法.方法将自由地表上的压缩波和剪切波速度设为零,而密度则43地 球 物 理 学 进 展23卷不为零,避免了由于被零除引起的计算震荡.这个方法发展的很早,使用的人也很多[46,51,55,57~62].这个方法不增加边界上的特殊计算量,尤其在水平的自由地表状态时效果很好,在起伏地表使用时,由于与自由地表条件方程无关,所以要小心使用.通常这个方法的精度和稳定性不高[62].例如,给出的自由边界条件是零应力条件:T z z =0, T xz =0.(5)在这种情况下,边界上只能使用一阶或二阶差分近似来离散,因为应力导数的计算需要边界上的零值区.2.3 广义虚象法该方法是在Levander [29]的虚象法的基础上发展起来能够处理起伏自由边界条件的技术.虚象法是在假设应力分量关于z =0对称,并向上沿拓两个网格基础上处理边界上的差分方程.这是在水平地表的情况下,而在处理起伏地表条件时要按照地表附近差分网格上变量及参数的分布确定其在网格的角内还是角外,然后按照水平和垂直量个方向向外拓展网格,称之为广义虚象法[46,51,55,57~62].可以看出这个方法增加了边界上的计算量,但是方法本身符合 虚象原理 ,经过证明,计算精度及稳定性要好于真空法.同样根据自由边界上的应力条件给出的网格节点分布如图5所示.图5 起伏自由边界上应力节点分布F ig.5 Lo cation of t he str ess in the vicinityo f an irr eg ular free sur face3 起伏地表间接处理方法起伏地表的间接处理方法,也可以称为全局法,则要改变波动方程,将一个在自然坐标系下含有地形起伏自由边界的弹性介质波动方程映射成为一个在规则坐标系下水平界面的方程,这个方程的形式已经改变了,然后用有限差分方法来离散方程求解.全局映射技术由Ekkehart T essmer 等[53]发表.该技术首先假设起伏的地表的函数是光滑的,并且其一阶导数是存在的,取一个曲线坐标系,这个坐标系的上部边界与地表函数一致,将这个坐标系映射于一个规则网格的直角坐标下,在空间导数上,横向上采用Chebychev 方法,垂向上Fo urier 法.尽管Tessmer 的自由界面上含地形变化的处理方法不能直在有限差分方法中使用,但是对于处理边界条件具有指导意义.目前的工作是应用这种变换,然后应用高阶、频散有界及最短差分算子来离散空间导数.在界面上的每一个点都要做相应的坐标变换,以便于自由边界上的法向矢量能够与垂向坐标轴方向重合.在边界上一般使用低阶的,中心交错网格差分就可以.H esthlom [11~13]及董良国[14],王祥春[64]等都是在这个技术的基础上略做改进完成了此种情况下的有限差分数值模拟.这个方法的核心思想就是首先是将起伏的地表映射到一个规则的长方形网格(#,∃)坐标系中,即x =x (#,∃),z =z (#,∃),并且定义#=0,#=#max ;∃=0,∃=∃max ,如图6所示.对于这种映射,采用1D 的线性延伸,即x (#,∃)=#,z (#,∃)=z 0(#)+∃∃max[z max -z 0(#)].(6)此处,z 0(#)定义了界面的形状,并且根据公式这个界面的形状变化是随着深度线性延迟变化的,直到z =z max 达到最大深度,底界面为平界面为止.这里要强调的是映射函数的条件是:单调的,充分光滑的.然后利用复合函数求导法则处理就得到变换后的方程.在这种空间中,x 是#的函数,z 是#和∃的函数,根据所给的模型的实际情况确定函数z 0(#),并估计给出映射函数z (#,∃),并且这个函数是整个网格的全局函数.这里给出的映射函数实现起来精度不高,并且难以选择.而H esthlom 在1994年和1998年分别给出了2D 和3D 的改进映射函数的有限差分数值模拟方法.他在2D 模型中给出的映射函数是:x (#,∃)=#,z (#,∃)=∃∃maxz 0(#).(7)这个映射函数的给出简化了映射的处理过程,441期王雪秋,等:地震波有限差分数值模拟框架下的起伏地表处理方法综述使复合函数的求导更加易于实现.同时自由边界条件方程也要采取同样的变换,得到一个新的坐标系下水平界面的自由边界方程.H esthlo m方法的优势是可以使用有限差分方法,处理灵活.缺陷是:地形函数必须是光滑的,具有一阶导数.如果函数在局部变化的过快,会引起吉布斯效应,甚至会导致模拟结果的震荡.图6 曲网格的拉伸映射示意Fig.6 M aping a curv ed g rid o nto a rectangular gr id with linear v ertical stretching4 评述与分析就有限差分方法解决复杂地表起伏问题来说,近20年来的主要成就就是高阶交错网格的广泛使用,这种方法一般用来解一阶速度应力方程,因为不必对弹性参数求导,从而简化了问题的难度并提高了计算速度.本文所讨论的就是在这种方法的基础上,对自由边界处理技术的实现方式.当然同时也要匹配相应的频散及吸收边界处理研究,对于这方面采用已有的技术就可以[69~70]直接法(局部法)是能够想到的一种简便的处理方法.在所有的直接法中,坐标旋转法处理符合边界上的应力情况的物理原理,但是需要对地表进行折线近似,增加计算量,并且处理起来不是很灵活并且在3D情况下很难处理.真空法处理快捷,灵活,但是在地表起伏时需要对地表进行规则的阶梯状近似,精度不高,稳定性差.广义虚象法借鉴了电磁场论中的 虚象法原理,经过证明是可行的,其精度和收敛性也要好于真空法.但是这些方法普遍存在一个收敛条件问题,也就是都不是无条件收敛的.间接法(全局法)是一种改变所有方程的处理方式,相对直接法而言,其计算不用考虑自由边界上特殊处理,按照常规水平自由边界处理即可.但是需要对自由地表进行连续光滑函数的近似,在实际中,这种条件很难达到,或者说我们想实现的就是那些地表剧烈起伏而非光滑的情况的波场模拟.其实所有问题中最重要的是自由边界条件的特殊处理,用纯计算技巧来处理则忽略了其对波场的物理意义,而考虑到物理意义则地表形状不可能是 任意的.如果考虑把某些方法联合起来使用,对问题的解决则不无裨益.例如将真空法和广义虚象法联合使用.总之,复杂自由地表的处理是数值模拟是研究近地表地震波场传播规律的有效手段,目前的所有处理自由地表方法都不能完全达到地表为 任意形状的条件,需要更深入的研究.参 考 文 献(References):[1] 刘光鼎,张丽莉,祝靓谊.试论复杂地质体的油气地震勘探[J].地球物理学进展,2006,21(3):683~686.Liu G D,Zh ang L L,Zhu L Y.S eism ic pros pecting for oil an d gas on the com plex geological bodies[J].Progress in Geo physics(in Ch ines e),2006,21(3):683~686.[2] 张向林,陶果,刘新茹.油气地球物理勘探技术进展[J].地球物理学进展,2006,21(1):143~151.Zhang X L,Tao G,Liu X R.Progress in oil geophysical ex ploration[J].Progress in Geophys ics(in C hinese),2006,21(1):143~151.[3] 佘德平,严建文,吴继敏.全波场地震勘探技术[J].地球物理学进展,2006,21(2):472~477.Sh e D P,Yan J W,W u J M.Th e full w ave field seismic ex ploration[J].Progres s in Geophy sics(in Chin ese),2006,21(2):472~477.[4] J.F.克莱鲍特.地震成像理论及方法[M].石油工业出版社,1991.45。

高阶交错网格有限差分弹性波场模拟的精度分析
岳晓鹏;白超英;岳崇旺
【期刊名称】《煤田地质与勘探》
【年(卷),期】2017(045)001
【摘要】交错网格波场数值模拟是目前地震正演中广泛使用的方法,为对比分析不同阶数的差分格式下产生的计算效率和精度差异,重新推导了弹性波方程的4种时间4阶、空间2N阶的差分公式及系数,并计算了他们的稳定性条件.利用这4种差分格式进行弹性波场数值模拟,对比分析了波场快照、合成地震记录及CPU时间.结果表明:时间4阶、空间6+6阶精度的交错网格有限差分方法在进行地震波场数值模拟时具有较高的计算精度和计算效率.
【总页数】6页(P125-130)
【作者】岳晓鹏;白超英;岳崇旺
【作者单位】长安大学地质工程与测绘学院地球物理系,陕西西安 710064;许昌学院数学与统计学院,河南许昌 461000;长安大学地质工程与测绘学院地球物理系,陕西西安 710064;长安大学地质工程与测绘学院地球物理系,陕西西安 710064【正文语种】中文
【中图分类】P315.69
【相关文献】
1.各向异性介质弹性波高阶交错网格有限差分模拟 [J], 霍凤斌;李振鹏;徐发;张涛;
2.高阶交错网格有限差分法弹性波叠前逆时深度偏移 [J], 陈可洋
3.二阶各向异性弹性波动方程高阶交错网格有限差分法 [J], 姜宇飞
4.弹性波交错网格高阶有限差分法波场分离数值模拟 [J], 李振春;张华;刘庆敏;韩文功
5.双相TI介质中弹性波交错网格高阶有限差分法数值模拟 [J], 尹学爱;邱光辉因版权原因，仅展示原文概要，查看原文内容请购买。

起伏地表下地震波场数值模拟计算摘要:本文简要总结了地震波场数值模拟的各种方法的基本原理并对比其优缺点，并对最近一些新出现的方法进行简要讨论。

模拟计算方法有很多种，本文重点讨论二维有限差分法，其能够有效地解决起伏地表的情况，是求解波动方程计算效率最高、应用最为广泛的方法之一。

关键词：地震波场数值模拟；二维有限差分法；起伏地表2D Numerical Simulation of Seismic Wave Field with Relief TopographyAbstract:This paper reviews the principles and characteristics of various numerical simulations of seismic wave field，and compares the merits and defects of the simulations. Some newly emerged methods and results are briefly discussed. There are so many computing methods for this, but the keynote of this text is talking about the 2D finite difference method. As one of the most efficient methods，the finite difference technique is widely used in wave equation solution and it is useful to resolve the problem which in surface topography.Key words:seismic wave field numerical simulation; 2D the finite difference method; surface topography地震波场数值模拟是研究地震波在介质中的传播规律、了解地震波在地下介质中的传播特点和协助对观测数据解释的有效手段，而提高计算精度和运算效率一向是所有波场数值模拟计算方法中所追求的目标。

一频率域弹性波场模拟中高阶有限差分法的精度对比与改进策略马超;沈金松;李曦宁【摘要】2D finite-difference operators in second-order as well as fourth-order accuracy based on staggered grid in frequen-cy-domain are derived. By combining optimal difference coefficients with lumped mass methods with a weighted average, nu-merical anisotropy can be suppressed, which enables accurate dispersion comparisons among these difference operators. U-sing a homogeneous isotropic medium, simulation precisions of different operators are compared to contrast the accuracy and error magnitude between second-order and fourth-order operators. It appears that the resulting error is smaller than 2% if the number of grid points required per smallest shear wavelength increases to 3. Meanwhile, in contrast to those based on con-ventional grid, difference operators based on staggered grid can be conducted in model simulations in fluids.%推导基于二维交错网格的二阶和四阶频率域弹性波有限差分算子，并结合最优化差分系数和质量加权平均方法压制数值各向异性。