（固体地球物理学专业论文）起伏地表条件下地震波传播数值模拟.pdf

摘要 摘要 山地地震勘探是当前勘探地球物理界面临的一个世界性难题，其中信噪比 低、静校正困难是最突出的两个问题，而它们均与地形起伏有直接的关系，复杂 的地形条件给地震数据采集、资料处理和解释提出了严峻挑战。要解决这些问题， 必须首先认识地形起伏条件下地震波传播的有关规律，它是地震数据合理采集、 正确处理的前提和基础。本论文目的就是用数值模拟方法模拟弹性波在含有起伏 地表模型中的传播过程，分析起伏地表对地震波的影响，为野外数据采集、资料 处理及解释提供模拟工具及理论指导。 起伏地表模拟最大的困难在于自由边界条件的实现方法。众所周知，有限差 分模拟方法已经广泛应用于地震波传播数值模拟中，但该方法的一个重要缺陷是 处理复杂地形比较困难，如何处理自由边界条件是关键所在。在本论文中，通过 在自由边界之上增加一层虚网格点来实现自由边界条件，对不同地形起伏情况下 各自由边界点和各虚网格点不同的离散方式的具体分析，将整个二维空间离散点 分成2 2 类。除内部点外，对另外2 1 类网格点，通过选择合理的自由边界条件实 现方式，分别提出了相应的具有针对性的差分格式，从而为利用有限差分方法进 行起伏地表条件下弹性波传播数值模拟铺平了道路。 通过数值试验发现，地表起伏使得地震波在近地表的传播变得异常复杂，引 起面波、体波等地震波型之间的相互转化，产生了大量的地表散射，从而引起山 地地震勘探中严重的低信噪比问题。通过弹性波数值模拟可以对各种地表干扰加 以识别，来指导实际地震资料的去噪处理。 关键词：起伏地表；数值模拟；有限差分；网格：自由边界 a b s t r a c t a bs t r a c t s e i s m i ce x p l o r a t i o ni nh i l l ya r e ai saw o r l d w i d ep r o b l e m ，i nw h i c hl o ws nr a t i o a n ds t a t i cc o r r e c t i o n ，d e d u c e dt os u r f a c et o p o g r a p h y , a r et h em o s tn o t o r i o u sp r o b l e m s c o m p l e xt o p o g r a p h yg i v e sc h a l l e n g et og e o p h y s i c i s t s t h i st h e s i sa i m st o t h e m o d e l i n go ft h ep r o p a g a t i o no fe l a s t i cw a v e si n t h ep r e s e n c eo ft o p o g r a p h y , p r e s e n t i n gt h ee f f e c t so ft h ei r r e g u l a rs u r f a c eo nt h ep r o p a g a t i o no ft h ew a v e s ，a n d g i v i n gat h e o r e t i c a lg u i d et os e i s m i ce x p l o r a t i o ni nh i h ya r e a t h em o s td i f f i c u l tp o i n tf o rs i m u l a t i o no fs u r f a c ew i t ht o p o g r a p h yi st h ef r e e s u r f a c ec o n d i t i o n sa n dt h ew a yo fc a r r y i n gt h e mo u t t h ef i n i t ed i f f e r e n c e ( f d ) h a s b e e nw i d e l yu s e di nt h es i m u l a t i o no fs e i s m i cw a v e , s ，e s p e c i a l l yi nt h es i m u l a t i o no f s e i s m i cw a v di nc o m p l e xm o d e l ，a n dt h ef i n i t ed i f f e r e n c e ( f d ) i sv e i ye f f i c i e n c y b u t t h em o s ti n f l u e n c ed i s a d v a n t a g ei st h a ti ti sh a r df o rt h ef i n i t ed i f f e r e n c e ( f d ) t od e a l w i t ht h es u r f a c ew i t ht o p o g r a p h y h o wt ou s et h ef i n i t ed i f f e r e n c e ( f d ) a p p r o a c ht o s i m u l a t et h ee l a s t i cw a v ep r o p a g a t i o nw i t ht o p o g r a p h yb e c o m e st h eb i g g e s tp r o b l e m i nt h i sp a p e r , al i n eo fp s e u d o p o i n t si sa d d e da b o v et h et o p o g r a p h yt oi m p l e m e n tt h e f r e e - s u r f a c eb o u n d a r yc o n d i t i o n s a f t e rt h ea n a l y s i sa b o u ta l lk i n d so ft h e s en o d e s ，t h e w h o l es p a c ei sd i v i d e dt o2 2t y p e s ，w h i c hh a v ed i f f e r e n tf i n i t ed i f f e r e n c e ( d ) s c h e m e s b yu s i n gt h i sa p p r o a c h ，f i n i t ed i f f e r e n c e ( f d ) m e t h o dc a nb eu s e dt o s i m u l a t et h ee l a s t i cw a v ep r o p a g a t i o ni nt h ep r e s e n c eo ft o p o g r a p h y n u m e r i c a l m o d e l i n g r e s u l t ss h o wt h a tt h e t o p o g r a p h yc o m p l i c a t e s t h e p r o p a g a t i o no fs e i s m i cw a v e si nn e a rs u r f a c e i tl e a d st ot h ec o n v e r s i o no fd i f f e r e n t w a v e ss u c ha ss u r f a c ew a v ea n db o d yw a v ea sw e l la sg e n e r a t i o no fs c a t t e r i n ga l o n g s u r f a c e t h i ss e v e r e l yd e t e r i o r a t e st h es no fs e i s m i ce x p l o r a t i o ni nm o u n t a i n o u s a r e a s s i m u l a t i o no fe l a s t i cw a v ef a c i l i t a t e st h er e c o g n i t i o na n de l i m i n a t i o no f d i f f e r e n ts u r f a c ei n t e r f e r e n c e s k e yw o r d s ：s u r f a c e w i t h t o p o g r a p h y , n u m e r i c a ls i m u l a t i o n ，f i n i t e d i f f e r e n c e ，m e s h ，f r e e s u r f a c e h 声明 本人郑重声明：本人在导师的指导下，独立进行研究工作所取得的成果， 撰写成硕士学位论文“起迭地麦玺鲑工丝震遮佳搔麴鱼槿塑：。除论文中已经注 明引用的内容外，对本文的研究做出重要贡献的个人和集体，均已在文中以明确 方式标明。本论文中不包含任何未加明确注明的其他个人或集体已经公开发表和 未公开发表的成果。 本声明的法律责任由本人承担。 易t 月 瓤 年 第一章前言 第一章前言弟一早刖i 本章首先介绍起伏地表条件下的地震波数值模拟的意义；然后介绍国内外在 该领域的研究现状以及新的动态；最后阐述本论文的研究内容及所取得的成果。 1 1 起伏地表条件下地震波数值模拟概述 二十一世纪，能源危机仍将威胁着世界，特别是中国，随着国民经济的飞速 发展，能源消耗量与日俱增，能源的供需矛盾愈来愈突出，能否解决这一问题直 接关系到国家下一步的发展。在新的替代能源还没有找到之前，石油、天然气仍 将是重要能源。随着我国东部油区的慢慢枯竭，我国油气勘探战略重点转移到了 西部。在我国西部地区，地形起伏十分剧烈，这对地震勘探工作提出了新的挑战。 不光野外数据采集难，资料的质量也普遍比较差，给处理工作带来了许多难题， 如非常突出的信噪比低、静校正困难等。 由于一些国内外成功的处理方法对该地区资料的处理都得不到令人满意的结 果，即使有，也是针对某一个小范围的，难以广泛使用。而且，山地地震勘探问 题除了南美北部以及北美西部个别地区外，其它地区比较少见，国外在这方面的 研究比较少，没有可以直接引进的方法。要解决这些问题，必须靠国内地球物理 工作者的努力，寻找新方法、新理论和新技术。 另外，在山地地震勘探中，地震资料的低信噪比问题比较突出，在原始单炮 记录上有效反射往往非常弱，资料信噪比一般低于沙漠地区地震资料。这其中除 了近地表岩性横向变化剧烈外，地表起伏是一个重要的原因。地形的变化引起地 震波在地表附近呈现出复杂的传播现象，产生的不同性质的地表干扰淹没了本来 就比较弱的深部反射信号。对于起伏地表产生的干扰波，实际地震资料处理中一 般采用f - k 滤波以及各类线性去噪方法进行处理，以提高资料的信噪比。但是山 地地震记录上各种震相比较复杂，这种不加区分的滤波，可能会将有效反射也滤 掉，也会将某些干扰误认为是有效波而加以保留。因此，在山地地震勘探中，为 了更有效地提高信噪比处理效果，同时也为了更好地设计面向地质目标的野外地 震观测系统，进行起伏地表条件下的地震波传播数值模拟，具体分析地表产生的 各种地震波传播现象和具体研究各类地表干扰波的形成机制是非常必要的。 地震波数值正演模拟是在建立合理的地质物理模型的基础上，用数学模型代 替真实介质，通过数学计算模拟地震记录形成过程，从而获得理论地震记录的手 段。 地震波数值正演模拟有许多方面的应用。它形象再现了地震波场的运动学、 动力学特征，有助于认识、理解地震波的传播规律，它能够帮助我们认识、理解 同济大学申请硕上学位论文起伏地表条件下地震波传播数值模拟吴晓丰 起伏地表对地震波的动力学、运动学特性的影响，为野外数据采集提供指导，为 新的去噪、静校方法提供思路和理论依据，并为这些处理方法研究提供实验数据。 此外，这些正演模拟技术可用于其它用途，如验证地质解释结果；油藏跟踪模拟， 可以预测油藏以后的演变，为生产单位提供决策依据。 1 2 国内外研究现状 随着能源工业对地震勘探要求的提高，弹性波勘探越来越受到人们的重视， 弹性波正演模拟技术也得到了很大发展。前人在对弹性介质性质描述及数值模拟 方法方面做了不少工作，有人甚至考虑了粘弹性及各向异性等复杂介质情况。这 里只介绍一些主要的模拟方法。 k o s l o f f ，d ，r e s h e f ，m ( 1 9 8 4 ) 采用差分法与f o u r i e r 方法的组合( 时间上厢 差分，空间上用f o u r i e r 方法) 。f o u r i e r 方法有精度高的优点，但是它要求速度变 化不能太剧烈。v l d e u x ，j ( 1 9 8 4 、1 9 8 6 ) 利用交错网格差分法求解弹性波速度一 应力方程。这种交错网格较规则网格具有明显优点，如计算精度高、可以直接模 拟液体一固体界面等。l e v a n d e r ( 1 9 8 8 ) 用交错网格有限差分法，在空间上取得了 四阶精度；c r a s e ( 1 9 9 0 ) 用高阶差分法解位移应力方程，将高阶差分的概念引入 了弹性波模拟，取得了时问、空间导数的任意高阶精度，从而提高了模拟效果： r o b e r t s s o n ( 1 9 9 5 、1 9 9 6 ) 用交错网格差分法模拟粘弹性介质中弹性波传播，等等。 在国内，地球物理勘探局研究院开发了有限元法数值模拟粘滞弹性介质中地震波 传播的软件( 1 9 8 8 ) ，并对弹性介质、粘弹性介质以及双相等效介质等进行了研究， 这是国内比较早的。董良国，马在田等( 2 0 0 0 ) 用交错网格高阶差分法解速度一应 力方程，并对该方法的稳定性进行了分析，该方法同时具有交错网格和高阶差分 高精度的优点，模拟效果很好。 地震波在地表传播是弹性波模拟的一个组成部分，但它又具有特殊性，尤其 是地表起伏情况下，地震波在地表传播的模拟研究对地震资料处理中提高信噪比 具有重要意义。 起伏地表条件下的地震波传播数值模拟在八十年代末才开始。在模拟方法问 题上，有限元( f e ) 方法计算效率低，但是处理不规则地表比较方便，成为模拟起 伏地表情况下地震波传播的有效方法。为此，发展了一些有限元和其他方法相结 合的混合方法，例如，m o c z o 等( 1 9 9 7 ) 用离散波数方法模拟震源激发和下部介质 中地震波的传播，继而通过有限元方法来模拟沿起伏地表波的传播( d w - f e ) ；黄 自萍等( 2 0 0 3 ) 用有限元和有限差分( f e e d ) 相结合的方法来模拟起伏地表地震波 传播。这类方法的主要问题是有限元( f e ) 和其他方法计算区域交界处容易产生人 为的反射，吸收边界也不容易处理。边界积分或边界元( b e ) 方法将散射波场通过 2 第一章前言 地表的一个半解析的积分来表示，其中积分项中g r e e n 函数一般在频率波数域中 计算。这种方法在研究起伏地表地震波传播时使用比较多，例如， b o u c h o n ( 1 9 9 6 ) p e d e r s c n ( 1 9 9 4 ) d u r a n d ( 1 9 9 9 ) 等就是用这种方法研究了地形的区域 变化对天然地震波场的影响，分析了地表散射以及地表振动随地形的变化，而 r u - s h a hw u 等( 1 9 9 8 ) 将边界元和屏传播算子相结合提出了一种混合模拟方法，在 处理起伏地表的同时提高了模拟效率。边界元( b e ) 方法的半解析性质决定了该方 法不能适用于地表速度变化较大情况，而实际情况是，由于后期地质作用造成浅 部地层速度变化更为剧烈，因此限制了边界元( b e ) 方法的实际使用。 有限差分( f d ) 方法计算效率高，在模拟复杂模型中地震波传播时应用最为广 泛，但是该方法的一个重要缺陷时处理复杂地形比较困难。t e s s e r , k o s l o f fa n d b e h l e ( 1 9 9 2 ) j 1 恿过坐标变换将具有起伏的地表的模型及弹性波方程( 速度应力方 程) 变换到新的坐标系中，在新坐标系中模型具有水平地表，然后将新的弹性波 方程离散化进行计算模拟，时间上用差分法，空间上横向用f o u d e r 法，纵向采 用c h e b y s h e v 方法，在上表面用自由边界条件处理。s t i gh e s t h o l ma n db e n tr u u d ( 1 9 9 4 ) 借鉴了t e s s e r 等的思路，通过坐标变换将起伏的地表转换成水平地表后， 用交错网格差分方法解弹性波的速度应力方程进行计算模拟。使用交错网格差 分法时，不是在每一个网格点上都同时计算出速度、应力各个分量，而是在某一 些特定的格点上计算速度的一个分量，某些另外的格点上计算另一个速度分量， 对应力亦如此。交错网格在自由界面水平是直接可以用，但是，地表起伏时，坐 标变换后速度应力方程变得更加复杂，在一些点上无法直接差分，必须通过 插值的方法才能将方程离散化，这无疑会降低精度。l e v a n d e r ( 1 9 8 8 ) 提出了一种 镜像法的新思路，他假定应力关于自由地表( 水平情况) 有某种对称关系，将应力 向表面以上延拓两个样点，这样就可以解自由边界条件求出表面处的速度分量， 他采用的是四阶交错网格差分，自由边界与理论非常相符合。c r a s e ( 1 9 9 2 ) 借鉴 了l e v a n d e r 的思路，提出了一种镜像法，假设质点位移关于自由表面是偶函数， 这样速度分量也应是偶对称的；j o l mo ar o b e r t s s o n ( 1 9 9 6 ) 又提出了一种镜像 法，并将它用到了起伏地表情况。他认为在自由表面之上质点的速度应为零。这 些想法来源于直觉，虽然，有些已取得了很好的效果，但目前仍缺乏足够的理论 依据。 j m ，r ，s ，m c l a u g h l i nk la n dd e rz a ( 1 9 8 8 ) 使用有限差分法( f d ) 模拟起伏 地表情况下地震波的传播。文章将地表起伏变化的各种情况归成了几类，如上倾、 下倾、陡坡、缓坡以及多种拐点等，并给出了每一类的差分方程格式，内点采用 常规中心差分。m a s a h i k of u y u k ia n dy o s h i r om a t s u m o t o ( 1 9 8 0 ) 用有限差分方法分 析面波在地堑构造起伏地表中面波的散射，为用有限差分( f d ) 模拟起伏地表的弹 3 同济人学申请硕士学位论文起伏地表条件下地震波传播数值模拟吴晓丰 性波传播提供了一种切实可行的方法。 1 3 论文研究内容及成果 在地形起伏情况下的弹性波传播数值模拟方面，两个关键问题是如何选择适 合不规则地表的数值计算方法以及如何更好地实现自由边界条件，实际模拟中这 两个问题也是有机的联系在一起的。 在利用有限差分( f d ) 方法模拟地表起伏情况下地震波传播方面，m a s a h i k o f u y u k ia n dy o s h i r om a t s u m o t o 在1 9 8 0 年模拟了地堑构造引起的面波散射。由于 假设的地堑构造在地表起伏处的垂直落差比较大，文中对不同起伏地表没有分别 进行处理。本文在其工作的基础上，通过在自由边界之上增加一层虚网格点来实 现自由边界条件，通过对不同地形起伏情况下各自由边界点和虚网格点不同离散 方式的具体分析，将整个二维空间离散点分成2 2 类。除内部点外，对另外2 1 类 网格点，通过选择合理的自由边界条件实现方式，分别提出了相应的差分格式， 从而为利用有限差分方法进行起伏地表条件下弹性波传播数值模拟铺平了道路。 论文分为四章，将详细介绍弹性波传播数值模拟中遇到的各种问题，第二章 主要介绍本文所采用的算法的具体原理，主要包括网格点的分类和每一类网格点 的差分格式及其构造过程。对网格点进行分类是本文所采用的算法的关键部分之 一，分类正确是后面计算正确的前提条件。在本章中将详细介绍网格点分类的依 据和意义。第三章对本文的方法进行模型试验，通过对一些有针对性的模型进行 试验来验证本文所采用的方法的可行性。在理论算法分析和模型试验的基础上， 第四章给出结论。 4 第二章方法原理 第二章方法原理 目前，鉴于有限差分( f d ) 方法具有计算精度高、模拟效率高、使用比较方便 的优点，它已成为地震波传播数值模拟中最常用的方法。但众所周知，有限差分 ( f d ) 方法处理复杂地形比较困难，如何利用有限差分( f d ) 方法来模拟地形起伏条 件下的地震波传播，一直是天然地震学和勘探地震学领域所努力的方向。 本章将系统介绍有限差分( f d ) 方法模拟起伏地表条件下的二维弹性波传播 的方法原理，重点介绍在处理自由边界时所采用的方法，包括计算策略、自由边 界的实施以及不同条件下差分格式的构造思路。同时介绍吸收边界处理方法以及 旨在提高模拟精度的高阶有限差分方法原理。 2 1 二维弹性波位移方程 声波方程是流体介质假设条件下推导出来的，不能用来研究地表起伏条件下 产生的地震波传播规律，需要利用弹性波方程模拟地表起伏时的各种地震波传播 现象。其中两个关键问题时如何选择不规则地表的数值计算方法以及如何更好地 实现自由边界条件。 为了研究问题的方便，假定介质为各向同性的完全弹性半空间，二维弹性波 位移方程为：( e w i n ge ta 1 ，1 9 5 7 ) u f f = 彳1 u 属+ 彳2 u 昭+ a 3 u 荔 ( 2 1 1 ) 其中位移向量u 等于： 矧 u 和w 分别为向量u 的水平分量和垂直分量， 彳一( ；善) ，彳2 。( v ；! y ；y ；：v ；) ， 其中y p 和y s 分别为p 波和s 波速度， 即： 并且系数矩阵为： 小仨 。10 、 害吖2 萨0 2 u + 眈q 蛲一譬 害一軎+ ；t 甓+ 嵋譬 ( 2 1 2 ) ( 2 1 3 ) ( 2 1 4 ) ( 2 1 5 ) 将以上方程离散，可以得到二阶精度差分格式( a t 、a r 和a z 分别是 5 同济大学申请硕士学位论文起伏地表条件下地震波传播数值模拟吴晓丰 时间和空间x 和z 方向的离散间隔) ： v ( i ，加粤妙( f + 1 ，小彬( f ，_ ) + u ( i 一1 ，j ) ) 缸 + 之妙( f ，_ + 1 ) 一撕( f ，j f ) + u ( f ，_ | 一1 ) ) a z + 志移( f + 1 ，_ + 1 ) 一v ( i + 1 ，一1 ) 一( ，( f 一1 ，+ 1 ) + u ( f 一1 ，j f 一1 ) ) ( 2 1 6 ) n ，+ l 嘞n - 和y ；2 恻一) 蚶叫扎2 咖v 孑2 慨，刮，嵋，) 2 2 自由边界条件 ( 2 1 7 ) 2 、 一y ；j ( 2 1 8 ) 起伏地表模型地震波传播数值模拟中的另一个主要问题就是自由边界条件及 其实现方法。 如图2 1 ，设0 点局部地形倾角为矽，原点定义在0 点的两个坐标系分 别为：x轴沿地表切向的xo z局部坐标系和x 轴沿水平方向的 x o z坐标系。根据张量左边变化关系，两个坐标系中应力分量之间满足 关系o r 玎= a 豫a 口o r 七，其中旋转矩阵为： 叫笼】- 【= 瑚 6 ，妒 盖， 1 4 4 + 、i 司a 一 v 川 w 嵋 + + m 厅b ，一犹吐 一 一 川 卜 厅b 吐比饥 气w 、l_l-， - 匐圹 s 。nb 生姚 卜0 一吐 嵋 町 办 鼍b 一 叫 。w 少峨封 嵋 第二章方法原理 lz ，气 ) 、 、 r 扫、 r 、 、 、 、 由此可以得到局部坐标系x o z 中0 点处的法向应力分量，分别为： 吒= t = s i n 2 帆一s i n ) 吒+ c o s 2 驴吒 一吒- 2 f 吒+ ( 厂。) 2 吒 鲁百而广( 2 2 - 1 ) = 吒= 一s i n c , c o s # a 材+ c o s ) 吒+ s i n # c o s t a ：也譬高灶卿) 1 + ( 厂。) 2 卜“叫 由a n n 皇o 和= 0 自由边界条件，可以得到x o z坐标系下自由边 界条件： 吒= ( 厂) 2 吒 ( 2 2 3 ) 0 口= f jo a q 。2 - 4 ) 对于水平地表，上述自由边界条件( 2 2 - 3 ) 和( 2 2 - 4 ) 变为： u z + b 1 u 工一0 ( 2 2 5 ) 对于垂直地表，上述自由边界条件( 2 2 - 3 ) 和( 2 2 - 4 ) 变为： u 工+ b 2 u z = 0 ( 2 2 6 ) 其中： 即b 2 。1 ) 啪 同时： b 22 曰r 2 ( 呈1 2 。毛7 v p ) 2 ) c 2 2 8 ， 7 同济大学申请硕士学位论文起伏地表条件下地震波传播数值模拟吴晓丰 2 3 自由边界点的分类及实现方法 由于地表是不规则的，本文在处理不规则地表时，将其离散后视为地堑结构， 由水平、垂直、9 0 0 拐角和2 7 0 0 拐角等情况组成，不同的类别分别采用不同的差 分格式，下面将详细介绍不规则地表网格点的分类以及分别采用的差分格式公式 推导。 2 3 1 水平条件下自由边界有限差分计算方法 如图2 2 ，点a 位于水平条件下的自由边界，有限差分计算n + l 时刻点a 的 波场值u 蛳) 时，需要知道n 时刻a 点上面一层虚拟网格点中的u ( i , j 一1 ) 、 u ( i 一1 , j 1 ) 、u ( i + l , j - 1 ) 的值，下面介绍如何计算这种条件下的虚拟网格点的波场 值。 l i , j + 1 ) 图2 2 水平条件f 自由地表有限差分网格 以点b 为例，在点a 利用水平条件下的自由边界条件( 2 2 5 ) 可以得到： u z + b i u x = 0 由中心差分得到： 三二垒- 掣+ 毋。吝u ( f ，- ) ：。 即： u ( f ，j 一1 ) = u ( f ，j + 1 ) + 2 a z b l d ；u ( i ，_ ) 于是得到计算b 点波场值的差分格式： 8 第二章方法原理 u ( f ，一1 ) f f iu ( i ，_ + 1 ) + 坳1 妙g + 1 j ) 一v ( i l _ ) ) x ( 2 3 1 ) 由( 2 3 1 ) 得到了n 时刻计算水平条件下虚拟网格点的差分格式，于是在计算 n + l 时刻水平条件下自由边界网格点时，可以利用和内部点相同的差分格式计 算，也就是说差分格式为公式( 2 1 6 ) 。 2 3 2 垂直条件下自由边界有限差分计算方法 垂直条件下分为：介质在左边和介质在右边两种情况。 当介质在左边时，如图2 3 ，网格点a 位于垂直条件下( 左) 的自由边界，有限 差分计算n + l 时刻点a 的波场值u ( i ，j ) 时，需要知道n 时刻a 点右边一层虚拟网 格点中的u ( i + 1 j ) 、u ( i + l , j 一1 ) 、u ( i + 1 j + 1 ) 的值，下面介绍如何计算这种条件下 的虚拟网格点的波场值。 i j - 1 ) 。( i 一1 j )i i ，j ) a 0 + , 1 5 b i j + 1 ) 图2 3 垂直条件下( 左) 自由地表有限差分网格 以点b 为例，在点a 利用垂直条件下的自由边界条件( 2 2 - 6 ) n - i 以得到： u 工+ 口2 u z ；0 由中心差分得到： 业掣掣+ b 2 d ；v ( i ，护o z x 即： v ( i + 1 ，j ) = v ( i 一1 j ) 一2 a x b 2 d 舌u ( f ，j f ) 于是得到计算b 点波场值的差分格式： v ( i + 1 ，j ) - v ( i 一1 j ) 一h x b 2 妙( f ，j + 1 ) 一v ( i ，j 一1 ) ) a z ( 2 3 2 ) 9 同济大学申请硕士学位论文 起伏地表条件下地震波传播数值模拟吴晓丰 由( 2 3 - 2 ) 得到了n 时刻计算垂直条件下( 左) 虚拟网格点的差分格式，于是在计 算n + l 时刻垂直条件下( 左) 自由边界网格点时，可以利用和内部点相同的差分格 式计算，即公式( 2 1 - 6 ) ，这里不再重复。 当介质在右边时，如图2 4 ， ( j - 1 ，j ) b i i - 1 ) i d ) ( i t l 加 a i d + 1 ) 图2 4 垂直条件下( 右) 自由地表有限差分网格 f l 了( 2 3 - 2 ) 可以得到： u ( f 一1 , j ) - u ( i + i , j ) + s x b 2 ( f ，j + 1 ) 一u ( f ，j 一1 ) ) a z ( 2 3 3 ) 同理，1 主i ( 2 3 3 ) 得到了n 时刻计算垂直条件下( 右) 虚拟网格点的差分格式，那 么计算n + l 时刻点a 同样和内部点一样计算，差分格式即公式( 2 1 6 ) 。 2 3 3 拐角为9 0 u 条件下自由边界有限差分计算方法 拐角为9 0 0 条件下分为：介质在左边和介质在右边两种情况。 当介质在左边时，如图2 5 ，网格点p 位于这种条件下自由边界的9 0 0 拐角处， 有限差分计算n + l 时刻点p 的波场值u ( i j ) 时，需要知道1 1 时刻虚拟网格点：点 a ：u ( i , j 1 ) 和点c ：u ( i + 1 j ) 以及点b 的值，下面介绍如何计算点a 和点c 这两 种类型虚拟网格点的波场值。 1 0 第二章方法原理 e ( i - 2 , j - 1 ) e ( i - l n ) b 图2 59 0 0 拐角条件下( 左) 自由地表有限差分网格 首先介绍如何计算点a ， 计算点a 的波场值时，在点d 利用自由边界条件可以得到： d i 妙( f ，j ) + u ( f 一2 ，j ) ) 2 + d ；b 1 u ( f 一】，_ ) 一0 ( 2 3 4 ) 即： u ( f ，j 一1 ) 一u ( f ，_ + 1 ) + u g 一2 ，j + 1 ) 一u ( f 一2 ，一1 ) + 4 z 既d 孑u ( f 一1 , j ) ( 2 3 5 ) 于是得到计算点a 波场值的差分格式： u ( f ，一1 ) 一u ( f ，j + 1 ) + u ( f 一2 ，j + 1 ) 一u ( f 一2 ，j 一1 ) + 2 a z b l ( f ，) 一u ( f 一2 ，j ) ) a x ( 2 3 6 ) 由( 2 3 - 6 ) 得到了n 时刻计算拐角为9 0 0 条件下( 左) 虚拟网格点a 的差分格式。 计算点c 的波场值时，在点e 利用自由边界条件可以得到： d 占妙( f ，d + v ( i ，j + 2 ) ) 2 + d i b 2 u ( f ，j + 1 ) 一0( 2 3 7 ) 即： u g + 1 ，) 一u ( f 一1 ) 一u ( i + l , j + 2 ) + u ( f l j + 2 ) 一4 a x b 2 0 , u ( i ，j + 1 ) ( 2 3 8 ) 得到计算点c 波场值的差分格式： u g + 1 ，j ) 一u ( f 一1 ，j ) 一u ( f + 1 ，j + 2 ) + u ( f lj + 2 ) 一2 a x b 2 妙( f ，j + 2 ) 一u ( f ，j ) ) l a z ( 2 3 9 ) 得到了n 时刻计算拐角为9 0 0 条件下( 左) 虚拟网格点c 的差分格式。 通过上面的计算得到了匕述条件下1 1 时刻的虚拟网格点的值，然后就应该计 t j l ” 萄 掣驴 妒 坐坠兰型塑塑型兰丝垒篁生一墼垡丝耋釜堡! 垫墨鎏箜堡鏊篁垡型 墨堕圭 算n + l 时刻的点p ，此时，如果使用和内部点相同的差分格式计算点p ，那么将 要用到点b 的值，但是点b 是无法计算的，为了避免计算点p 时使用到点b 的 值，应用如下有限差分方法： 将u 麓一1 2 ( 、d ，x d l z + d 二d 二b o ，)( 2 3 1 0 ) 代入( 2 1 1 ) ，即用l ，f d 、，x d ，z + d x d z ，l i 。，j ) i 弋替u x z ， 即： c ，。毒妙m 小刃+ u o 刮” + 毒1 ) 一嬲+ u ( f ，) ) + 五耄塞妙( f + 1 ，_ + 1 ) 一u ( f ，_ + 1 ) 一u o + 1 ，) + 2 u ( f ，) 一u ( f 一1 ，j ) 一u ( f ，一1 ) + u ( f 一1 ，_ 一1 ) ) 于是得到计算拐角为9 0 。条件下( 左) 自由边界网格点p 的差分格式为1 2 3 。1 1 “：；1 - 咄：j 1 + 勉：玎尝) 一知一a ，) + y ；( 尝) 。c l ! m 一劾：，+ “! 一) + 乏1 缸t & 2 1 、，p 2 一v 孑) ( 2 3 1 2 ) 坛u + 1 一嵋川一吆u + 2 吆，一眩u w 抽+ 吆一1 ) 蚶一嵋l + 2 嵋；( 尝) 2 慨，一2 嵋，+ ) + v ；( 尝) 一2 w ；：，+ 吆h ) + 三2 生5 x a z6 ，“)嘲3 ) b 鼻l ，+ l 一甜，+ 1 一i n + l j + 2 “：，一“三1 ，一口，- 1 + 翻厶，一1 ) 当介质在右边时，如图2 6 ，网格点p 位于这种条件下自由边界的9 0 。拐角处， 有限差分计算n + l 时刻点p 的波场值u ( i j ) 时，需要知道n 时刻虚拟网格点：点 a ( u ( i , j _ 1 ) ) 和点c ( u ( i - 1 j ”的值，下面介绍如何计算点a 和点c 这两种类型虚拟 网格点的波场值。 1 2 第二章方法原理 b ( i - l ，i , j ) l ( i - l , j + 1 ) ( 0 ，j + 2 ) f )_ ”“列1 ) 天 。 。 图2 69 0 0 拐角条件下( 右) 自由地表有限差分网格 首先介绍如何计算点a ， 计算点a 的波场值时，利用点d 的自由边界条件可以得到： d 舌移( f ，j ) + u ( f + 2 ，j ) ) 2 + d 吝b 1 u ( f + 1 j ) = 0 那么n 时刻计算拐角为9 0 0 条件下( 右) 虚拟网格点a 的差分格式： u ( f ，j 一1 ) = u ( f + 2 ，j + 1 ) + u ( f ，j + 1 ) 一u ( f + 2 ，j 一1 ) + 2 a z ( u ( i + 2 ，i ) 一u ( f ，j ”缸 ( 2 3 - 1 4 ) ( 2 3 1 5 ) 计算点c 的波场值时，利用点e 的自由边界条件可以得到： d 舌妙( f ，j f ) + u ( f ，j + 2 ) ) 2 + d i b 2 u ( f ，+ 1 ) 一0 ( 2 3 1 6 ) 那么n 时刻计算拐角为9 0 0 条件下( 右) 虚拟网格点c 的差分格式： u ( f 一1 ，) = u ( f + 1 ，) + u ( f + 1 ，j + 2 ) 一u ( 一1 , j + 2 ) + 2 a x ( u ( i ，歹+ 2 ) 一u ( ，j ) ) & ( 2 3 1 7 ) 同理，得到了上述条件下n 时刻的虚拟网格点的值之后，接下来计算n + l 时 刻的点p ，同样，为了避免计算点p 时使用到点b 的值，应用如下有限差分方法： 将u 起；丢：d 三+ d 二d ；砂( f ，歹)( 2 3 1 8 ) 代入( 2 - 1 - 1 ) ，即用丢( d ：d 三+ d 二d ；b ( f ，j ) 代替， 即： 1 3 同济大学申请硕士学位论文起伏地表条件下地震波传播数值模拟吴晓丰 u 一专妙( f 州) 一彬+ u ( f 一1 ，歹” + 鲁妙( f ，j + 1 ) 一彬( f ，小u ( f ，卜1 ) ) + 孟妙( f ，+ 1 ) 一u ( 一1 ，i f + 1 ) 一嬲+ u ( 一1 ，小【，( f + 1 小u ( f + 1 ，| 一1 ) + u ( f ，| 一1 ” ( 2 3 - 1 9 ) 于是得到计算拐角为9 0 0 条件下( 右) 自由边界网格点p 的差分格式为： n ，+ l 一_ u 叫n 7 1 嘞n + y ；( 尝) g n _ 2 u n 一) + y ；斟u n m 一知：一汕吾生a x a z 6 ，“) ( 2 3 2 0 ) 似一逸u + 。一2 嵋+ 畦。吆u 一吆u 以+ 嵋- 1 ) 磁卜嵋l + 2 咖y ；( 尝) 2 k ，叫，帕，) + v ；( 甜一2 w ：朋w n 小丢生a x a z ；) 卿1 ) n ，+ 1 一u n l ，+ l 一2 u i “, j + 口二1 ，+ u 。n + 1 ，一“磊1 卜1 + u i n , j lj 2 3 4 拐角为2 7 0 u 条件下自由边界有限差分计算方法 拐角为2 7 0 0 条件下分为：介质在左边和介质在右边两种情况。 当介质在左边时，如图2 7 ，网格点q 位于这种条件下自由边界的2 7 0 0 拐角 处，网格点e 、f 、h 和i 位于2 7 0 0 拐角的两边上，有限差分计算n + l 时刻点e 、 f 、h 和i 的波场时，需要知道1 1 时刻虚拟网格点g 的值，下面介绍如何计算点 g 这种类型虚拟网格点的波场值。 1 4 第二章方法原理 ( i - l , j 1 ) + ( + 1 州) ( i + 咖1 ) 图2 72 7 0 0 拐角条件下( 左) 自由地表有限差分网格 由图2 7 可以看到，计算e 、f 两点的波场和h 、i 两点的波场时都要用到点 g 的波场值，计算e 、f 两点时用到的点g 值通过在点f 应用垂直条件下的自由 边界条件确定，而计算h 、i 两点时用到的点g 值通过在点h 应用水平条件下的 自由边界条件确定，也就是说点g 需要有两个值分别为：g v 和g h 。 在点f 应用垂直条件下自由边界条件( 2 2 6 ) ，1 主1 ( 2 3 2 ) 可以得到： u p ( f + 1 ，j 一1 ) = u ( f 一1 ，j 一1 ) 一x b 2 妙( f ，) 一u ( f ，_ 一2 ) ) a z ( 2 3 2 2 ) 在点h 应用水平条件下自由边界条件( 2 2 5 ) ，f l t ( 2 3 1 ) 可以得到： u 死( f + 1 ，j 一1 ) = u ( f + l j + 1 ) + a z b l 妙( f + 2 ，j ) 一u ( f ，j ”x ( 2 3 2 3 ) 由( 2 3 - 2 2 ) 和( 2 3 - 2 3 ) n - - j 以得到2 7 0 0 拐角条件下( 左) 虚网格点的值 然后可以计算自由地表点q 、e 、f 、h 和i 的值： 计算网格点q 时，为了避免使用到点g 的值，使用( 2 3 1 1 ) 确定的差分格式； 计算网格点e 时，使用内部点差分格式( 2 1 6 ) ，但是使用的点g 的值由( 2 3 2 2 ) 确定，可以得到： u o ，) 一粤p ( f + 1 ，) 一刃g ，j ) + u ( f 一1 _ ) ) + 鲁妙( f ，_ + 1 ) 一嬲g ，小u ( f ，一1 ) ) & + 志帆( f + 1 j f + 1 ) 一【，( i + 1 ，卜1 ) 一【，( f 一1 ，) + 【，( f 一1 ，) ) ( 2 3 - 2 4 ) 计算网格点f 时，为了避免使用到点h 的值，使用( 2 3 1 9 ) 确定的差分格式， 1 5 同济大学申请硕士学位论文起伏地表条件下地震波传播数值模拟 吴晓丰 同时，使用的点g 的值f l j ( 2 3 2 2 ) 确定，可以得到： u ( i ，_ ) 。之( f + l ) 一彬( f ，) + u ( f 一1 ，_ ” x + 之妙( f ，j + 1 ) 一嬲( f ，小v ( i ，卜1 ” z + 去移( f ，+ 1 ) 一v ( i 一1 ，+ 1 ) 一2 v ( i ，小【，( j l 小v v ( i + 1 小v ( i + 1 ，| 一1 ) + v ( i ，一1 ) ) ( 2 3 - 2 5 ) 计算网格点i 时，使用内部点差分格式( 2 1 - 6 ) ，但是使用的点g 的值由( 2 3 2 3 ) 确定，可以得到： v ( i ，- ) 粤妙g + 1 _ ) 一彬( f ，) + v ( i l _ ” 缸。 + 丢妙( f ，- + 1 ) 一2 v ( i ，小v ( i ，卜1 ) ) z + 去妙( f + 1 ，m ) 一u ( f +