1、地震波有限差分数值模拟框架下的起伏地表处理方法综述第 23 卷第 1 期2008 年 2 月(页码:4O48)地球物理学进展PR0GRESSINGE0PHYSICSVo1.23No.1Feb.2008地震波有限差分数值模拟框架下的起伏地表处理方法综述王雪秋,孙建国.(1.吉林大学地球探测科学与技术学院,长春 130026;2.国土资源部应用地球物理综合解释理论开放实验室,长春 130026)摘要本文查阅了近 3O 年来国内外关于有限差分方法模拟近地表复杂介质方面的文献,分析了在有限差分法框架下复杂地表地形起伏问题的处理及实现方法.总结了地震渡有限差分数值模拟方面的成果.复杂地表类型很多,对地震
2、资料的影响程度也不尽相同,要想完全描述近地表影响因素是不可能的,而关键的问题就是自由边界条件方程的处理.所有处理及实现方法中不外两种:直接法(局部 )和间接(全局)法.文中列出的方法,各有特点.我们的目的是寻找合适的方法解决适合的问题.关键词有限差分,自由边界,直接法,间接法,映射,地表几何形状中图分类号 P631 文献标识码 A 文章编号 10042903(2008)02004009Thestate-of-the-artinnumericalmodelingincludingsurfacetopographywithfinite?differencemethodWANGXueqiu,SUNJ
3、ianguo(1.CollegeofGeoexplorationsciencbandtechnology,JilinUniversity,Changchun130026,China;2.LaboratoryofwavetheoryandimagingtechnologyopenresearchlaboratoryforintegratedgP0pscfn PrprP 0ntheoryoftheMinistryoflandresources,Changchun130026,China)AbstractAsatoolforunderstandingtheregulationfwavepropaga
4、tionincomplexmediaandexaminingtheeffectsofdifferentkindsofmethods,seismicnumericalmodelingplaysanimportantroleinmoderngeophysicalprospecting.Thispaperbrieflyintroducedseveralmodelingmethodsincludingsurfacetopographywithfinite-differenceinrecent30years.Twoimplementationmethodsareproposed:derectmethod
5、whichincludesrotatedcoordinatemethod,vacuummethodandgeneralimagingmethodandindirectmethod.Ouraimistofindthesuitablemethodtosolvetheadaptedquestions.Keywordsfinite-difference,freesurface,directmethod,indirectmethod,mappingsurfacetopography0 引言经典反射地震勘探的理论基础有 3 个假设前提:地表是水平的;地下介质是水平层状的;地震波是近垂直反射的.这种假设简化
6、了理论模型和参数,使得到的地震数据便于解释,并且诸如 CDP 技术,NMO 和 DMO 以及偏移成像技术等都是以这些假设前提为基础建立起来的,这些技术事实上也在很多近似符合条件的地区得到了非常好的结果.但是当在勘探领域向复杂地区深入,勘探技术不断进步,大尺度长排列的应用等情况下前面的假设前提失去了意义,尤其是在我国西部地区复杂的地质条件进行的反射地震勘探,实际情况不允许这些假设条件的存在1.另外 ,在研究地球介质的性质,地震灾害性质,金属矿勘探及环境工程勘察等方面,更加促使我们抛弃这些假设前提,而注重于更加实际的波场的精细刻画.因此研究复杂地区的波场特征及传播规律成为发展的必然.尤其是复杂地表
7、及近地表收稿日期 20070310;修回日期 20070620.基金项目国家自然科学基金项目(40574052)和教育部骨干教师资助计划项目联合资助.作者简介王雪秋,女,1974 年生,吉林大学地球探测科学与技术学院讲师,在读博士,主要从事地震波理论及波场正演模拟等方面的研究(Emall:)1 期王雪秋,等:地震波有限差分数值模拟框架下的起伏地表处理方法综述 41情况对波场的影响规律的研究更是具有重要意义.这是由于我们多数时候都是在地表附近接收到的地震数据,所以地表的形状及介质物理特性对波场的改造相对于一些深层的信息来说是剧烈的.而关于近地表的地震波场“散射“理论则成为研究的热点_5.一般情况
8、下,静校正可以消除这种影响,但是多数情况下尤其是近地表地质地形非常复杂的地区,静校正的作用是有限的 l_l 卜,而地震波的数值模拟则是完成这种认识的一种直接有效的手段.对复杂地表问题的地震波“散射“ 场的解,一般有两种形式;数值解和解析解.解析解能够解决的地表模型一般比较简单_lq 引,所以对复杂的地表情况,采用数值解是必需的.复杂地表一般指两种情况:复杂的几何形状,例如山地,隆起,河谷,凹陷及高差剧烈,地形多变等;近地表介质复杂的物理性质,例如弹性 ,粘弹性,孔隙性,非均匀性和随机性等.前一种情况就是自由边界条件,也就是来自地下的应力,应变或速度消失的界面.而在一般给出的自由边界条件方程中,
9、都是给出沿界面的法向分量和垂向分量为零的方程.在水平界面情况下,坐标轴与界面重合,这个方程很好处理.但是,当界面与与坐标轴不重合的时候,也就是地表有几何形状起伏时,这个方程的处理就很困难.后一种情况实际上就是介质模型的确定,也即方程的确定.根据地区的地质状况可以采用很多种介质模型的近似,例如弹性波方程,粘弹性介质方程,非线性方程等.而在任何一种方程中都存在自由边界问题,可以说,介质方程的确定是自由边界条件存在的载体.另外一个需要特别指出的问题是,在经典的地震勘探理论中,由于近垂直反射假设条件的前提,人们一般选用声学介质模型来近似地下介质,由于参数简单,易于处理,许多成熟的,在工业界广泛使用的数
10、据处理技术都是以声波模型为基础的.但是声波方程是在流体假设条件推导出来的,不存在起伏形状自由边界问题,不符合实际近地表地震波场研究模型条件.所以复杂地表问题最核心的问题就是自由地表几何形状处理问题.关于地表几何形状对波场产生的影响问题,早在上个世纪 40 年代,M.B.Widess_28就提出含地形的自由地表响应是影响地震填图的主要因素之一,并且是地震构造解释的主要误差之一.只是由于当时的计算技术落后,无法开展这项工作.20 世纪 60 年代后期,随着计算技术软硬件的发展,地震波的数值模拟技术也迅猛发展起来 35.在所有的模拟方法中网格法一直由于其对地形和介质的灵活划分而占主导地位,因此也成为
11、几何状地形起伏地表地震波数值模拟的首选方法.在进行近地表地形模型的网格剖分时,最重要的是考虑如何对地表进行有益而准确的网格剖分,这种剖分既要符合所采用的方法原理,也要符合波动传播的物理机制,最好不要产生不必要的由计算方法而引起的波动现象.在网格法中,针对任意起伏的地表使用有限元方法来模拟近地表的波的传播,能够灵活处理边界的几何形状问题,但是计算比较慢,较其他方法的成本要高.另外一个成本不高,又能够处理自由地表几何形状的方法是边界元方法.已经有人将这个方法应用在了弹性介质的波场模拟.对于内部结构相对简单的介质(均匀或线性变化),这个方法在节省计算内存时相对于其他方法占有优势,原因在于计算是沿着边
12、界来完成的.边界可以是复杂的形状,只要曲线闭合即可.但是仅适用于地表几何形状的模拟,对于介质内部的复杂的岩性和构造则很难处理.而一些谱方法则要求自由边界函数具有光滑的特性,并且会出现一些譬如折叠效应,吉布斯效应等不期望出现的效果_3.在网格法中,有限差分方法是较早开展起来的模拟方法,也是一种比较简单明确的模拟方法.在上个世纪的 60 年代美国的 Alterman 及其研究小组就应用这种方法来解层状介质的弹性波方程.同时Claerbout 教授在研究地震成像中也提出了有限差分模拟方法,并给出了公开的算法.具有里程碑意义的是 R.M.Alford 等人对有限差分模拟精度问题的探讨,在他们的研究成果
13、中详细的讨论了网格点取值与差分算子的关系,同时还使用了高阶差分(四阶)333.到了上个世纪 70 年代,K.R.Kelly 等人给出了 2 维弹性波方程的完整的有限差分格式,并且给出了当时比较复杂的模型(分块的层,断层等)的时间切片和模拟记录图323.进入 80 年代,又发展出了一种时间和空间上可达任意阶的高阶交错网格有限差分法,这种方法的优势是 l_4:(1)震源以粒子振动速度或应力的形式表达,可以直接代入方程;(2)能够精确的处理自由边界条件;(3)差分算子是局部的,计算模型不必驻留内存,节省内存空间,节约成本,就算速度快;(4)方便使用高阶算子,提高计算精度;(5)易于与其他方法联合使用
14、;(6)算法容易实现.因此越来越受到重视.而研究者一直也在探讨其在含几何形状起伏的自由地表条件的地震波场数值模拟中的应用.,42 地球物理学进展在早期的常规有限差分地震波数值模拟中,一般要求自由地表是水平的或者是 45 度倾斜界面.但是有关含地表形状的自由边界处理方法也出现的很早.MIT 的 Boore 在其博士论文中除了应用有限差分方法来解弹性波方程外,还首次应用有限差分方法研究了倾斜界面中 LOVE 波的传播机制_4.Ilan 等_49就用有限差分方法模拟了 P-SV 波在起伏自由边界弹性介质中的传播,但是在线性分段之间的过度点没有进行处理.并且这种方法要求用不规则网格引起一些不必要的复杂
15、性和在迭代过程中计算精度不够问题.1981年,Boore 等应用 Ilan 的方法模拟了垂直入射的体波在45 度倾斜界面的碰撞反射,但是没有处理角点.这些处理手段都是由低阶有限差分方法来实现的,多数情况下受差分阶数及显式实现算法的影响,模拟结果的精度和收敛性都不是很高_5.而现在在使用有限差分方法来进行近地表数值模拟时,人们多数采用高阶交错网格来离散方程_4.在所有关于含起伏地表几何形状自由边界条件的地震波数值模拟方法中,我们可以依据方法对自由边界的处理方法将其分为两类.一类是直接处理自由边界条件,而与所用介质的波动方程无关,称之为局部法,也可以叫做直接法_4.另一类是对介质模型做变换,即将实
16、际上曲线(曲面)坐标系下的模型方程,变换到规则网格的坐标系下,同时自由边界条件方程也如此变换,我们称这种方式为全局法,也可以叫做间接法_1 卜.1 自由边界条件自由边界是实际上就是固体地球与空气接触的表面,这是一个地球介质的不连续界面.这个界面也是来自介质内部的应力应变消失的界面,同时这个界面也强烈地影响着下面的固体介质内的波场的一些物理现象,因为地表的任何起伏都能引起在介质中传播的波场参数的剧烈变化.在数学表达上,自由边界一般都是给出的都是按照牛顿力学原理的应力或位移微分方程,这个边界方程一般是一阶的 l_6,例如:+一 of121+一 o,(1).,dZd这是一个 2D 的自由边界条件,式
17、中【,W 是水平和垂直方向的位移,a, 分别是纵横波速度 .这个自由边界条件方程与应力为零是等价的.3D 的边界方程与其类似.对这个方程应用有限差分方法时,由于在边界上要降低差分网格的阶数,往往造成了计算结果精度不够或计算结果的不稳定,尤其是在近地表介质物理性质横向剧烈变化或地形起伏剧烈的情况下.因此,在地震波数值模拟中自由边界是对波场影响最大的边界.同时这种影响也能够以最小误差来描述,因为其它层位的参数和位置的不确定因素要远大于自由边界,而自由边界是不需要定位的,它是事实上存在的.因此一旦关于自由地表的准确描述的数值模拟技术可以应用,其前景将十分可观.一般情况下,自由边界的数值模拟(包括地球
18、的不均匀性)产生的影响都能自动在理论合成的结果中提供,并且所提供结果比不加自由边界的更加真实.因此在有限差分法地震波数值模拟中自由边界的处理技巧不但涉及到计算的精度,还涉及到方法的稳定及结果的准确性_1.“.在有限差分法的常规模拟中,自由边界条件有三种处理方式,一种是自由边界节点位于差分网格线上,这是一种常用的方法_3.“,如图 1 所示 .另一种是自由边界节点位于差分网格线之间_6,见图 2 所示.节点位于网格上,便于实现,尤其是在显式差分实现方式中;但是稳定性差,精度不高,因为一般只能采用低阶差分格式.而节点在网格之间可以采用隐式差分来实现,是无条件稳定的.前面两种方法都可以在起伏地表处理
19、中使用,而第三种则是由Levander_2 胡提出的“虚象法“,即假设应力分量关于 z 一 0 对称 ,并向上沿拓两个网格 ,使其自然满足自由边界条件,类似于电磁场中的虚象原理,但是仅能应用在水平的自由地表条件中.杏由 1 凸 r1 凸 r1 凸 r 自由边界图 1 自由边界节点在网格上Fig.1Thenodesoffreesurfaceboundaryonthegrids一自由 j 直界图 2 自由边界节点在网格问Fig.2Thenodesoffreesurfaceboundaryinthegrids1 期王雪秋,等:地震波有限差分数值模拟框架下的起伏地表处理方法综述 432 起伏地表直接处
20、理方法起伏地表的直接处理方法也称为局部处理方法.能够直接处理起伏地表自由边界条件而不必考虑介质方程的方法主要有三种:局部旋转坐标法,真空法,广义虚象法.2.1 局部旋转坐标法局部旋转坐标法是一种出现比较早的处理自由边界起伏的方法,1975 年,Ilan4 叼等用有限差分方法模拟 PSV 波在起伏自由边界弹性介质中的传播时就采用这种方法来实现自由边界条件的处理,1988 年,Jih5 叼等人发表的成果将这个方法进一步完善.其原理如图 3 所示,把不规则的地表地形分段线性化,也就是将自由边界按不同的斜率划分成线性段.然后使用局部旋转坐标系的方式(平行于分段的倾斜界面)对应力边界条件进行插值.在两个
21、线性段的衔接处在平分拐角的基础之上使用一阶近似的边界条件.图 3 分段起伏地表示意Fig.3Cartesiancoordinatesystemforthe2Dfinitedifferenceschemeandmiscellaneouscornerpointsonthetopography对于给定如图 3 这样曲折的地表,进行数值模拟需要进行一些特殊的处理.目前的处理方式是利用边界条件(方程 2.1)处理图中所示的角点.这种问题的处理早在 1970 年的 Alterman 和 Loewenthal 就提出了两种计算 90 度和 270 度角点处位移的方法:第一种处理方式是设此时水平和垂直正应力为零,Jt:ffOU 一一 0;第二种处理方式aZaZ是对角点进行轻微的平滑,使之保持在 45 度以内,见图 4./,/l,:I 二lI/IlI
