（地球探测与信息技术专业论文）二维地震正演模拟方法技术研究.pdf

摘要 地震正演模拟是地震数据采集、处理、解释三大环节的分析基础，可为地震数据采 集、处理、解释提供理论依据以及科学的评估方法。通过地震正演模拟，可以检验采集 设计的合理性，处理和解释成果的可靠性以及反演方法和结果的正确性。 正演模拟的主要研究方法分为射线追踪法和波动方程法。射线追踪法主要理论基础 是几何光学，在高频近似条件下，地震波的主能量沿射线轨迹传播，主要优点是计算速 度快，所得地震波的传播时间比较准确。 采用波动方程理论模拟地震波的传播，不仅能保持地震波的运动学特征，还能保持 动力学特征。有限差分法是基于差分原理的一种数值计算法。有限差分法是一种求解波 动方程的快速有效方法，可以很好地适应剧烈变化的地下介质情况，而且地质模型的复 杂程度并不影响运算速度。 本文在详细分析了地震数值模拟基本理论和方法技术的基础上，针对二维正演模拟 进行了研究，详细论述了射线追踪方法和波动方程方法的基本原理，并分别论述了其基 本算法。其中在二维波动方程的正演模拟中，对不同阶精度有限差分算子进行了推导， 对差分格式的稳定性及边界条件进行了讨论。最终在w i n d o w s 平台上，利用c + + 和 f o r t r a n 语言，并在前人工作的基础上，完成了两点射线追踪试射法和时间方向四阶精度 有限差分正演算法功能模块，经多种模型的测试，计算结果正确，算法稳定，适应性较 好。 文中结合能源勘探和工程勘查中可能存在的地质问题，设计了多种地质模型，包括 断层模型、透镜体模型、薄层模型以及空洞模型等，分别进行了射线和有限差分法正演 模拟，并对其效果进行了分析和比较。结果表明，两种算法均可进行一般地质问题的正 演模拟，其中有限差分算法对复杂模型的适应性更好，模拟结果的精度和可靠性较高， 波场动力学特征变化可信，可直接用于辅助地震勘探( 查) 设计、资料处理和解释。 关键词：正演模拟、射线追踪、波动方程、有限差分、透镜体、空洞 a b s t r a c t s e i s m i cf o r w a r dm o d e l i n gi st h ea n a l y s i sb a s i so fs e i s m i cd a t aa c q u i s i t i o n ，p r o c e s s i n g a n di n t e r p r e t a t i o n ，a n dc a np r o v i d et h e o r e t i c a lf o u n d a t i o na n ds c i e n t i f i ca s s e s s m e n tf o r s e i s m i cd a t aa c q u i s i t i o n ，p r o c e s s i n ga n di n t e r p r e t a t i o n b ys e i s m i cf o r w a r dm o d e l i n g ，w ec a n c h e c kt h e r a t i o n a l i t yo fd e s i g n ，a n d t h er e l i a b i l i t yo ft h er e s u l t so fp r o c e s s i n ga n d i n t e r p r e t a t i o n ，a n dt h ev a l i d i t yo fi n v e r s i o nm e t h o da n dr e s u l t s t h em o s t l yr e s e a r c hm e t h o d so ff o r w a r dm o d e l i n ga r er a yt r a c i n gm e t h o da n dw a v e e q u a t i o nm e t h o d r a yt r a c i n gm e t h o di sb a s e do ng e o m e t r i c a lo p t i c s u n d e rt h ec o n d i t i o no f h i g h 行e q u e n c ya p p r o x i m a t i o n ，t h em a i ne n e r g yo fs e i s m i cw a v es p r e a d sa l o n gt h er a yt r a c k , r a yt r a c i n gm e t h o dh a sv e l o c i t ye x c e l l e n c y ，a n dt h et r a v e lt i m ei sn i c e t y w a v ee q u a t i o nt h e o r yi su s e dt os i m u l a t et h es p r e a do fs e i s m i cw a v e ，c a nk e e pt h e k i n e m a t i c sa n dd y n a m i cc h a r a c t e ro fs e i s m i cw a v e ，a n da l s ob ea b l et ok e e pt h ed y n a m i c c h a r a c t e r f i n i t ed i f f e r e n c em e t h o di sas o r to fn u m e r i c a lm e t h o d ，b a s e do nd i f f e r e n c et h e o r y f i n i t ed i f f e r e n c em e t h o di saf a s ta n da v a i l a b l ea p p r o a c ht os o l v ew a v ee q u a t i o n ，c a na d a p t t h ed r a s t i cc h a n g e so ft h eu n d e r g r o u n dm e d i a ，m o r e o v e rt h ei n t r i c a c yo fm o d e ld o e sn o t i m p a c tt h eo p e r a t i o ns p e e d i nt h i sa r t i c l e ，b a s i n go na n a l y z i n gt h eb a s i ct h e o r ya n dm e t h o dt e c h n o l o g yo fs e i s m i c n u m e r i c a ls i m u l a t i o n ，a n dw eh a v er e s e a r c h e dt h e2 一df o r w a r dm o d e l i n g ，a n dd i s c u s s e di n d e t a i lr a yt r a c i n gm e t h o da n dw a v ee q u a t i o nm e t h o dt h e o r y ，a n dd i s c u s s e dt h e i rb a s i c a l g o r i t h m s i n2 - dw a v ee q u a t i o nf o r w a r dm o d e l i n g ，w eh a v ed e d u c t e dt h ef i n i t ed i f f e r e n c e o p e r a t o rw i t hd i f f e r e n to r d e rp r e c i s i o n ，a n dd i s c u s s e ds t a b i l i t yo fd i f f e r e n c es c h e m ea n d b o u n d a r yc o n d i t i o n o nt h ew i n d o w sp l a t f o r m ，u s i n gc + + l a n g u a g ea n df o r t r a n ，a n d b a s e do nt h ef o r e r u n n e r s r e s e a r c h ，w eh a v ef i n i s h e dt h ef u n c t i o n a lm o d u l e so ft w op o i n tr a y t r a c i n gs h o o t i n gm e t h o da n df o u ro r d e rt e m p o r a lp r e c i s i o nf i n i t ed i f f e r e n c e ，a n dc a r r i e do u t m u l t i m o d e l st e s t s ，a n dp r o v e dt h a tt h er e s u l tb er i g h t ，a n dt h es t a b i l i z a t i o na n d a d a p t a b i l i t yo f a l g o r i t h m i nt h i sa r t i c l e ，b yg e o l o g i c a lp r o b l e m st h a tm a yb ee x i s t e di ne n e r g yp r o s p e c t i n g a n de n g i n e e r i n gs u r v e ye x i s t e n t ，w ed e s i g n e dg e o l o g i c a lm o d e l s ，i n c l u d ef a u l tm o d e l ，l e n s m o d e l ，t h i nl a y e rm o d e la n dc a v i t ym o d e la n ds oo n ，a n dp r o c e e d e dw i t hr a yt r a c i n ga n d f i n i t ed i f f e r e n c em e t h o df o r w a r dm o d e l i n g ，a n da n a l y z e da n dc o m p a r e dt h er e s u l t s t h er e s u l t i n d i c a t et h a tt h e s et w om e t h o d sc a nb ep r o c e s s e di nt h ef o r w a r d m o d e l i n go fg e n e r a l g e o l o g i c a lp r o b l e m s ，a n dt h ef i n i t ed i f f e r e n c em e t h o dh a v eb e t t e ra d a p t a b i l i t yw i t hi n t r i c a c y m o d e l ，h i g h e rp r e c i s i o na n dr e l i a b i l i t yo fr e s u l t s ，c r e d i b i l i t yo fd y n a m i cc h a r a c t e ro fs e i s m i c w a v e ，a n dc a l lb ea p p r o v e dd i r e c t l yt os e i s m i cs u r v e yd e s i g n ，d a t ap r o c e s s i n ga n d i n t e r p r e t a t i o n k e yw o r d s ：f o r w a r dm o d e l i n g ；r a yt r a c i n g ；w a v ee q u a t i o n ；f i n i t ed i f f e r e n c e ；l e n s ；c a v i t y 论文独创性声明 本人声明：本人所呈交的学位论文是在导师的指导下，独立进行研究工 作所取得的成果。除论文中已经注明引用的内容外，对论文的研究做出重 要贡献的个人和集体，均已在文中以明确方式标明。本论文中不包含任何 未加明确注明的其他个人或集体已经公开发表的成果。 本声明的法律责任由本人承担。 论文作者签名：弗孚 2 卯7 年乒月秒日 论文知识产权权属声明 本人在导师指导下所完成的论文及相关的职务作品，知识产权归属学 校。学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权 利。本人离校后发表或使用学位论文或与该论文直接相关的学术论文或成 果时，署名单位仍然为长安大学。 ( 保密的论文在解密后应遵守此规定) ：南_ 年笋其夸e l 列年名月芗日 步 岳叮、川 书砖 名 孔 别 名 者 签 作 节 划 师 长安大学硕十学位论文 1 1 研究意义及目的 第一章引言 正演模拟是地震勘探方法中非常重要的一项技术。地震正演模拟就是求取地震波在 已知的地下地质模型中的传播规律，包括传播时间、路径、能量等。通过模型正演，可 以正确认识地震波在复杂介质中传播的运动学和动力学特征，准确地分析地下地质构造 所产生的反射地震波场特征。在人工给定的地质模型中进行正演得到的规律能够增进人 们对未知模型的认识，从而有助于地质问题的解决。 地震模拟技术可分为地震物理模拟与数值模拟两种，其中地震数值模拟是地震勘探 和地震学的重要基础。所谓地震数值模拟就是在假定地下介质结构模型和相应物理参数 已知的情况下，模拟研究地震波在地下各种介质中的传播规律，并计算在地面或地下各 观测点所应观测到的数值地震记录的一种地震模拟方法，其理论基础就是表征地震波在 地下各种介质中传播的地震波传播理论。地震数值模拟方法己在地震勘探和天然地震领 域中得到广泛的应用。它不但在石油、天然气、煤、金属和非金属等矿产资源及工程和 环境地球物理中得到普遍的应用，而且在地震灾害预测、地震区带划分以及地壳构造和 地球内部结构研究中，也得到相当广泛的应用。地震数值模拟在地震勘探和地震学各工 作阶段中都有重要的作用。在地震数据采集设计中，地震数值模拟可用于野外地震观测 系统的设计和评估，并进行地震观测系统的优化。在地震数据处理中，地震数值模拟可 以检验各种反演方法的正确性。在地震数据处理结果的解释中，地震数值模拟又可以对 地震解释结果的正确性进行检验。 总的来说，地震正演模拟有以下几个方面的用处： ( 1 ) 对波场特征进行理论研究； ( 2 ) 指导野外数据采集工作设计； ( 3 ) 作为反演解释的理论基础； ( 4 ) 帮助资料解释。 本论文在详细分析射线理论和波动理论的原理和方法基础上，针对目前比较流行的 地震正演模拟方法进行研究和对比，并结合能源勘探和工程勘查领域中比较典型的地质 和灾害问题，设计模型，进行正演模拟，比较不同正演方法的优缺点，为更好地指导野 外数据采集、资料处理、解释和反演工作奠定基础。 第一章引言 1 2 正演模拟的基本方法及研究现状 1 2 1 正演模拟基本理论 地球物理正问题是在给出地球模型，包括全球或局部的地质结构，在假定的自然源 或人工源的作用下，用相应的地球物理场方程求解地面或地下某处的地球物理场的响应 函数【9 1 。 地震波的正演过程数学上可描写为 d = 彳( 所) ( 1 1 ) 其中，d 一合成地震数据向量；a 一正演算子；m 一模型向量。d 的精度受m 的离 散化精度和正演算子彳计算精度的影响。 地震波传播方程正演问题的研究内容： ( 1 ) 地震波数值模拟原理； ( 2 ) 地震波数值模拟算法； ( 3 ) 计算程序的质量。它受三个主要因素影响，即计算精度，计算效率和计算的稳 定性。 图1 1 为地震数值模拟方法示意图。 图1 1 地震数值模拟方法示意图 地震数值模拟的实现方法 1 ) 地质一地球物理模型离散化 不同的研究目的，构成地质一地球物理模型的地球物理量也不同对于均匀各向同性 介质声波方程，地质一地球物理模型m 可表示为聊= m ( p ，) ，而对均匀各向同性介质 2 长安大学硕士学位论文 弹性波方程而言，其地质一地球物理模型则表示为肌= m ( p ，v p ，匕) 即模型m 离散化为包 含有限个模型物理参数的向量m 。 地质一地球物理模型的离散化是通过对模型的空间剖分来实现的。目前地质一地球 物理模型的空间剖分方法可分为两种，即正交网格剖分和非正交网格剖分。正交网格在 平面上是矩形网格，而非正交网格剖分在平面上是三角形网格和不规则四边形网格。对 于地下介质进行非正交网格剖分可以充分考虑到地下介质分布的几何形状，并且不受边 界几何形状的限制。从这一点来说，基于非正交网格的数值模拟方法要优于基于正交网 格的数值模拟方法。 为了准确地刻画地下介质物理性质空间变化，模型网格剖分必须要求足够精细，只 有这样才能以较小的误差逼近地下介质的实际物性分布。模型网格剖分的越细，则空间 网格点的数目越多，这必然会占用大量的内存量，加大计算的工作量，降低计算的效率， 增加计算的成本。因此，模型离散化时既要考虑数值模拟的分辨率( 或网格大小) 又要考 虑计算成本。 2 ) 地震波方程的离散化 模型空间的网格化必然带来波场分量的网格化。由于这种网格化把一个连续的地震 波动问题转化成一个离散的地震波动问题，因此必然涉及到波场逼近问题，并在空间网 格化以后，尽可能以较小的逼近误差表示离散波场的空间微分。图1 2 表示地震数值模 拟各种概念之间的关系【1 8 】。 ( h ，a t 0 ) 图1 2 数值模拟各种概念之间的关系 1 2 2 研究现状 地震波场数值模拟是研究复杂地区地震资料采集、处理和解释的有效辅助手段，地 震波场数值模拟的主要方法包括两大类，即几何射线法和波动方程法。波动方程数值模 3 第一章引言 拟方法实质上是求解地震波波动方程，因此模拟的地震波场包含了地震波传播的所有信 息，但其计算速度相对于几何射线法要慢。几何射线法也就是射线追踪法，属于几何地 震学方法，由于它将地震波波动理论简化为射线理论，主要考虑的是地震波传播的运动 学特征，缺少地震波的动力学信息，因此该方法计算速度快。因为波动方程模拟包含了 丰富的波动信息，为研究地震波的传播机理和复杂地层的解释提供了更多的佐证，所以 波动方程数值模拟方法一直在地震模拟中占有重要地位【1 1 。 1 、几何射线法 几何射线法又称为射线追踪法，其主要理论基础是，在高频近似条件下，地震波的 主能量沿射线线路径问题。基于这种认识，运用惠更斯原理和费马原理来重建射线路径， 并利用程函方程来计算射线的旅行时，以期获得快速的解。 射线理论的优点是计算速度快，所得地震波的传播时间比较准确，但是计算结果很 难保持动力学特征，对复杂的地质构造会出现盲区，对于复杂构造进行两点三维射线追 踪往往比较麻烦。 在地震学中，有两类地震射线追踪问题：一类是1 点射线追踪，即已知射线初始点 位置和初始出射方向求地震波的传播路径问题；另一类是两点射线追踪，即已知射线初 始点和另一观察点接收点的位置，不知射线初始出射方向，求两点之间的射线路径问题。 1 ) 有关初值问题的试射法。这种方法根据由震源出发的一束射线到达接收点的情 况对射线出射角及其密度进行调整，最后由最靠近接收点的两条射线走时内插求出接收 点走时。 2 ) 有关边值问题的弯曲法。这种方法从震源与接收点之间的一条假想初始路径开 始，根据最小时准则对路径进行扰动，从而求出接收点处的走时及射线路径。费马时间 稳定原理是弯曲法依据的基本原理，即在固定的两个端点之间，波实际传播的路径是在 真实射线附近变动的路径中，能使波的旅行时稳定的路径，也就是说，时间取极值，通 常是极小值的路径。 以上两种方法的局限是难于处理介质中较强的速度变化，难于求出多值走时中的全 局最小走时，计算效率低，阴影区内射线覆盖密度不足。 2 0 世纪8 0 年代末以来，随着k i r c h h o f f 积分叠前深度偏移在解决复杂构造成像中 获得一系列成功，作为其算法基础之一的射线追踪方法也得到了很大的促进和发展，出 现了大量不同于传统算法的新型算法。这些算法的主要特点在于不再局限于地震波的路 径描述，而是直接从惠更斯原理或费马原理出发，采用等价的波前描述地震波场的特征。 4 长安人学硕士学位论文 近年来随着三维勘探的发展，为了适应三维数据处理和复杂地质体研究，新的射线追踪 算法不断涌现。这些改进主要有：在传统的试射法和弯曲法的基础之上产生了很多改进 的射线追踪方法，如波前重建法；最小走时算法的改进，使之可适应多值走时计算，如 慢度匹配法。下面对这些方法作一简要的介绍。 1 ) v i d a l e 方法。与传统方法不同，v i d a l e 方法计算的是波阵面而不是射线路径， 以二维情况为例，用正方形的网格对速度模型进行离散化。 2 ) 改进的v i d a l e 方法。q i n ( 1 9 9 2 ) 在v i d a l e 扩展方阵的基础上实现了扩展波前的 递推方法，但计算量增加很大；p o d v i n ( 1 9 9 1 ) 也是按扩展方阵的方式求取走时，对每 一个网格点，系统地比较来自各个方向的透射波、衍射波和首波。 3 ) v a nt r i e r 法。v a nt r i e r ( 1 9 9 1 ) 首先将程函方程化为守恒型程函方程，然后用有 限差分上风法直接求解变换后的方程，进而求出地震波场的最小走时。 4 ) w f r t 法。w f r t 方法的基本出发点为惠更斯原理。根据介质的非均匀程度，将 所要研究的介质分割成大小相等的矩形网格，每个矩形网格单元内的速度可视为均匀 的，称为第一次分割；然后再根据计算精度的要求将每一个矩形网格进一步分成等份的 小矩形网格，称为第二次分割。 5 ) 最短路径法。最短路径法的基础是费马原理及图论中的最短路径理论。利用最 短路径的思路求解程函方程，很多人都作了大量工作，这些方法的基本思路相同，只不 过具体实现步骤上存在差异。 6 ) s a v a 和f o m e l ( 1 9 9 8 ) 结合传统射线追踪方法和程函方程地有限差分法，提出了 h w t ( h u y g e n sw a v e f r o n tt r a c i n g ) 法。 7 ) 慢度匹配法。慢度匹配法的目的仍是求解多值走时问题，s y m e s ( 1 9 9 8 ) 假设地震 波的射线都是下行的，该算法的计算量很大，因此在目前还不是一种可实用化的方法。 8 ) 三维多值走时地震波场重建及格林函数计算。这种方法根据三维复杂介质地震 波传播数值计算中出现的多值走时情况，以及地震波场重建和格林函数计算中的困难， 提出在相位空间拉格朗日流形上的波场重建及格林函数计算方法。该方法基于波前重建 思路，通过在相位空间中将位形空间中具有多值且不可微的射线参数曲面展为单值可微 曲面，以克服地震波场射线解的复杂性态。 辛几何算法。在非均匀介质中，常常需要用m a s l o v 渐进理论解决求地震波路径 时出现的焦散问题。m a s l o v 渐进理论需要对空间坐标进行傅氏变换，即将部分空间坐标 变为慢度坐标。因此，在数值计算中，尽量保持位置和慢向量之间的关系是十分重要的。 5 第一章引言 辛几何算法正是保持他们之间关系的算法【。 2 、波动方程法 采用波动方程理论模拟地震波的传播，不仅能保持地震波的运动学特征，而且还能 保持地震波的动力学特征。求解波动方程的方法有多种，有限差分法是一种求解波动方 程的快速有效方法，而且地质模型的复杂程度并不影响运算速度；虚谱法是对波动方程 的波场函数进行傅氏变换，其优点是运算速度快；克希霍夫积分法是基于均匀介质模型， 利用格林公式计算曲面积分，求出空间波场值，但这种方法不能适应变速介质；有限元 法能对复杂地质构造的地质模型进行地震波场数值模拟，但计算速度较慢。 ( 1 ) 有限差分方法 有限差分法通过有限差分算子将波动方程离散化，以差分代替微分，将微分方程问 题化为代数问题，然后求解相关的线性代数方程组以获得微分方程问题的数值解。有限 差分法只考虑有限个离散点上的函数值，而不考虑点的领域函数值如何变化，差分算子 是个空间局部的算子，在空间域具有较高的分辨率，可以很好地适应剧烈变化的地下介 质情况。但是在频率域中，有限差分算子的分辨率却很低。算法的稳定性和收敛性受空 间采样率和时间采样率的影响，但算法的速度较快。 有限差分数值模拟技术是声波或弹性波场数值模拟中最为流行的方法之一，现在正 向提高精度的方向发展。a l t e r m a n 和k a r a l ( 1 9 6 8 ) 首先将有限差分法应用于层状介质 弹性波传播的数值模拟中，之后，a l f o r d 等研究了声波方程有限差分法的精度，v i r i e u x ， 提出了稳定的二阶空间和时间弹性波有限差分格式，它适用于任何泊松比的介质， l e v a n d e r 将v i r i e u x 的方法推广到空间四阶、时间二阶的情况。c r a s e 则发展了精度可 达任意阶的高阶交错网格法，提高了地震模拟的精度和稳定性，但其计算量和内存要求 比低阶有限差分法大幅度增加。g r a v e s ( 1 9 9 6 ) 给出了三维速度一应力方程交错网格有 限差分法弹性波传播的模拟方法。三维速度一应力方程交错网格有限差分弹性波数值模 拟方法使得地震数值模拟对地震波在弹性介质中传播模拟能力有了明显提高。 为了进一步模拟地震波在非完全弹性的实际地层中的传播，c a r c i o n e 等( 1 9 8 8 ) 提 出了粘滞声波在地层中传播的模拟方法。r o b e r t s s o n 等( 1 9 9 4 ) 给出了粘弹性波有限差 分模拟方法。而为了模拟地震波在实际地层中的各向异性性质。i g e l 等( 1 9 9 5 ) 给出了 各向异性介质中地震波传播的有限差分模拟方法，此后，为了模拟地震波在双相或多相 孔隙弹性介质中的传播，d a i 等( 1 9 9 5 ) 提出了非均匀孔隙介质中速度应力方程有限 差分地震模拟方法。o z d e n v a r 和m c m e c h a n ( 1 9 9 7 ) 给出了孔隙弹性介质中地震波的交 6 长安大学硕士学位论文 错网格有限差分模拟方法。c a r c i o n e 和h e l l e ( 1 9 9 9 ) 提出了孔隙粘弹性介质中地震波传 播的交错网格有限差分模拟方法。董良国等给出了一阶弹性波方程交错网格高阶差分解 法。裴正林等给出了三维任意各向异性介质中弹性波传播的的交错网格高阶有限差分模 拟方法。 一般的有限差分地震模拟方法基于笛卡尔坐标系中的规则网格，在模拟复杂地质构 造和复杂地质体的复杂界面时，为了减弱必然会出现阶梯状的边界，在这种边界上必然 引起人为的虚假绕射波，这种虚假绕射波，则必须采用精细网格，而这将导致数值模拟 计算储存量的增加和计算量的增大为此，则必须发展基于可变网格和不规则网格的地震 数值模拟方法。 j a s t r a m 和t e s s m e r ( 1 9 9 4 ) 提出了垂直间距可变网格的弹性波模拟方法。张剑锋， o p r s a l 和z a h r a d n i k ( 1 9 9 9 ) 提出了非均匀介质中弹性波的矩形不规则网格有限差分法模 拟方法。p i t a r k a 给出了三维各向同性介质中弹性波的矩形非规则交错网格有限差分模拟 方法。n o r d s t r o m 和c a r p e n t e r 提出了曲线坐标下变形网格高阶有限差分法地震数值模拟 方法。 随着山地、沙漠等地表复杂区地震勘探的发展，起伏地表地震数值模拟技术受到了 地球物理勘探学家的广泛关注和重视。h e s t h o l m 和r u u d ( 1 9 9 4 ) 对起伏地表下弹性波 传播有限差分法数值模拟进行了研究。r o b e r t s s o n 给出了含起伏地表弹性和粘弹性介质 中地震波传播的有限差分模拟方法，裴正林等对任意起伏地表下2 d 弹性波传播交错网 格高阶有限差分法数值模拟进行了研究【1 】【1 刚。 ( 2 ) 虚谱法 有限差分法的算法简单快速，但难以克服频散效应，而要解决频散问题，须加密数 值计算的网格这势必会导致计算量增加，效率下降。因此，选择一种既能精确计算，又 有较高计算效率的方法就显得非常必要，而虚谱法正好符合这种要求。虚谱法是一种有 效的数值模拟方法，权衡精度和效率，有其他方法不可替代的优点。 虚谱法利用f o u r i e r 变换将波场函数表示为f o u r i e r 级数的展开形式，将波动方程在 时间一波数域或频率中求解，虚谱法对微分算子的逼近程度可以达到n y q u i s t 频率或波 数，并且收敛速度快。但由于f o u r i e r 变换是基于整个时间域或空间域的，改变空间中 的某一点的值，就会改变频率域中的所有值，因此每一点的微分结果都要受到计算域中 其它点的影响并且存在g i b b s 效应，虚谱法在频率域分辨率高，在时间域分辨率低。虚 谱法无法处理空间物性剧烈变化的情形。 7 第一章引言 g a z d a g ( 1 9 8 1 ) 提出了声波传播的虚谱法地震模拟方法。k o d o f f 等( 1 9 8 2 ) 给出了 虚谱法地震波模拟方法。f o m b e r g ( 1 9 8 7 ) ，对弹性波传播模拟的虚谱法和有限差分法进 行了对比研究。r e s h e f 等( 1 9 8 2 ) 给出了三维声波模拟的虚谱法。n i e l s e n 等将伪谱法 用到了曲线网格的二维声波方程的模拟中。t e s s m e n 等则用伪谱法在存在表面地形的情 况下模拟了弹性波的传播。 ( 3 ) 有限元法 有限元法也是正演模拟的有效手段。由于剖分的任意性及它所依据的变分原理，对 含有多种介质和自然边界条件的处理，非常方便有效，已成为解决地震波传播数值模拟 的一种重要方法。它是目前为止最精确的一种正演模拟方法，但它对计算机内存要求很 高，计算量大。同时为满足不同研究的需要，应当设置不同的边界条件。 有限元法的主要优点是适宜于模拟任意地质体形态，可以任意三角形逼近地层界 面，保证复杂地层形态模拟的逼真性。l y s m e r 和d r a k e ( 1 9 7 2 ) 最早将有限元法应用于 地震数值模拟。m a r f u r t ( 1 9 8 4 ) 研究对比了模拟弹性波传播的有限差分法和有限元法的 精度。s e r o n 等( 1 9 9 0 ，1 9 9 6 ) 给出了弹性波传播有限元模拟方法。p a d o v a n i 等( 1 9 9 4 ) 研究了地震波模拟的低阶和高阶有限元法。s a r m a 等( 1 9 8 8 ) 给出了弹性波传播有限元 模拟中的无反射边界条件。王尚旭利用有限元法研究了双向各向同性介质中弹性波传播 规律。周辉等给出了各向异性介质中波动方程的有限元数值解法。杨顶辉对孔隙各向异 性介质中基于b i s q 模型的弹性波传播进行了有限元方法模拟。 近年来，国内外许多著名学者在复杂介质中地震波传播问题的研究方面已作了很多 杰出的工作。以吴如山等为代表发展了广义屏方法，基于d e w o l f 近似推导了双域方程 的表达式，将复杂介质用背景速度和扰动速度近似，相应的波场用背景场和扰动场近 似，背景场在波数域中延拓，扰动场在空间域延拓，同时处理边界等。这种方法对于地 下复杂的介质己取得明显的效果。 3 、常用的地震数值模拟方法分类 表1 1 给出了常用的地震数值模拟方法分类比较。 8 长安大学硕上学位论文 表1 1 数值模拟方法分类 分类理论方法优点缺点 基于惠更斯原理和斯试射法 能避开射线盲区， 复杂模型收敛慢， 射线追踪数值奈尔定律 实现精确追踪 计算效率低 模拟模拟地震波场的运动 对于简单模型，计 复杂模型，有时会 学特征 弯曲法陷入局部收敛，得 算效率高 到非全局最优解 运算速度快，实现频散严重，难应付 普通有限差分法 简单 起伏地表问题 精度高，运算速度难应付起伏地表 高阶有限差分法 快，实现简单问题 基于弹性或粘弹性理 有限元法 处理地表起伏问 计算量大，计算效 波动方程数值 论和牛顿力学题强率低 模拟 傅氏变换法运算速度快难应付横向变速 模拟地震波场的动力 计算量大，计算效 学特征 谱元法处理边界灵活 率低 表征问题准确，适 构造的方程系非 边界元法 应各种边界条件常大，计算营大 介质速度变化小， 广义屏法运算速度快 入射角小 1 3 主要研究内容及研究成果 1 3 1 主要研究内容 在详细分析地震数值模拟的基本原理和方法技术的基础上，借助前人部分工作成 果，针对2 d 地质模型开展地震模拟的研究工作，主要内容包括： ( 1 ) 射线追踪和有限差分方法理论分析和数学模型的建立； ( 2 ) 根据能源勘探和环境与工程勘查领域中的典型地质和灾害问题，总结并建立 相关的地质模型； ( 3 ) 利用可视化编程环境，完成上述相关算法的正演功能模块； ( 4 ) 结合方法原理及理论，以及各种模型的正演结果，分析各方法之间的效果优 劣及适用性，同时对有限差分的边界条件进行讨论。 1 3 2 主要研究成果 经过近一年的论文工作，基本完成了论文开题所涉及的内容，主要获得了下述成果： ( 1 ) 利用可视化编程环境，完成了两点射线追踪和有限差分法正演功能模块，经 过各种模型的计算验证，算法实现正确，效果良好； 9 第一章引言 ( 2 ) 采用插值加圆滑的方法完成了速度建模模块，为有限差分正演的实现奠定了 基础； ( 3 ) 在相同时间精度条件下，经过对不同空间精度的比较，分析认为4 阶有限差 分方法在不考虑运算速度的情况下，效果较好； ( 4 ) 在有限差分的边界处理上，经过试验和对比分析认为采用阻尼吸收边界在同 时计算速度和精度的要求上比c l a y t o n 吸收边界效果好。 ( 5 ) 分析并建立了多种典型地质模型，包括断层模型、透镜体模型、薄层模型以 及空洞模型等，获得了较丰富的地震模拟结果，为指导工程实践奠定了基础。 1 0 长安人学硕上学位论文 第二章地震波传播的基本理论 2 1 地震勘探地层假设 地震波是在地下岩层中传播的弹性波，其传播规律满足波动方程。由于实际岩层的 非均匀性，使得波动方程的边界条件非常复杂，以至于无法解出该波动方程。在一定条 件下，用射线理论来近似波动理论，使数学问题大大简化，从而能方便的研究地震波在 地下介质中的传播路径、传播速度、旅行时间等运动学特征。 地震波在岩层中的传播速度与岩层的性质，诸如弹性参数、岩石成分、密度、埋藏 深度、地质年代、孔隙度等因素有关，这导致了地震速度分布规律及波传播规律的复杂 性，给研究问题造成了困难。因此，地震勘探理论一般把地下介质作以下几种近似： 均匀介质：这种介质假设地震波在岩层各处的传播速度都相同。 层状介质：假设岩层是成层分布的，地震波在每一层中的传播速度都是相同的，而 各层之间的地震速度不同。 连续介质：假设地震波在岩层中的传播速度是随空间位置连续变化的。 连续介质的物性随空间坐标连续变化，显然是非均匀介质。层状介质实际上也是一 种非均匀介质，只不过物性的不均匀只出现在分界面处，是一种比较特殊的非均匀介质。 岩层介质结构不同，地震波的传播规律也会不同。对地震波在不同介质中传播特征 的研究，是地震勘探的重要内容，也是进行定性和定量解释的依据。层状介质模型虽然 只是实际介质的一种近似，但它仍然是地震勘探中最常用的物理模型。 2 2 惠更斯原理 惠更斯原理可做如下描述：在弹性介质中，已知t 时刻的同一波前面上的各点，可 以把这些点看作从该时刻产生的新的点震源，经过任何一个t 时刻后，这些子波的包 络面就是原波到达的t + a t 时刻新的波前面。根据这个原理，只要知道某一时刻的波前 位置，就能够确定出地震波在各种不同时间的波前位置。 惠更斯原理只给出了波传播的空间几何位置，而没有涉及到波到达该位置的物理状 态。菲涅尔补充了惠更斯原理，他指出，从同一波阵面上的各点所发出的子波，经传播 而在空间相遇时，可以相互叠加而产生干涉现象，从而在该点观测到的是总扰动。当波 遇到另一种介质时，这个原理可以形象地说明波的反射、透射、折射和绕射等现象。 第二章地震波传播的基本理论 2 3 费马原理 费马原理是从射线角度描述波传播特点的。它指出，地震波沿射线传播的时间和沿 其它任何路径传播的时间比较为最小，或者说波沿传播时间为最小的路径传播。根据这 个原理可以确定地震波在已知传播速度的介质中的射线形状，如在均匀介质中射线路径 为直线，而在非均匀介质中，射线路径不再是直线，是与波前面垂直的曲线。 2 4 地震波在层状介质中的传播规律 2 4 1 斯奈尔定理 当地震波传至两种不同介质的分界面时，一部分能量在原介质形成反射波，另一部 分能量透过分界面，在第二介质中传播形成透射波。s n e l l 定理指出： ( 1 ) 入射线、透射线位于界面法线的两侧，入射线、透射线和法线在同一射线平面内； ( 2 ) 入射角o l 的正弦和透射角岛的正弦之比，等于入射波速度k 与透射波速度k 之 比，即 掣。掣 ( 2 1 ) k 、 2 4 2 反射与透射 当地震波入射到两种介质的分界面时，通常会分成两部分。一部分回到原介质中形 成反射波；另一部分透入第二种介质中传播，在物理学中叫做折射波，而在地震学中习 惯叫做透射波。反射波和透射波的情况如图2 1 所示。设有两种介质的分界面，用n 和 p ：分别代表第一种介质和第二种介质的质量密度，用h 和 ，：分别代表波在两种介质中的 传播速度。第一种介质的波阻抗是z 1 = m ，第二种介质的波阻抗是z ：= 胁屹。研究证 明，只有在z l - z ：的条件下，地震波才会发生反射；z 1 和z ：的差别越大，反射波越强。 长安大学硕上学位论文 反射i 2 介1 介质特 8 l 夕 p 、 介质2 令弋 运 图2 1 入射波、反射波、透射波 2 4 3 地震波在多层介质中的传播规律 地震波在多层介质中传播时，当波到达上层界面时，一部分能量反射，另一部分能 量透射。设入射振幅为单位1 的波，则反射振幅为r - 的波返回地面，透射振幅为( 1 一墨) 的波到达第二界面。当振幅为( 1 一r ) 的波到达第二界面时，又会产生反射和透射，反射 波的振幅为( 1 一r ) r ，透射波的振幅为( 1 一r ) ( 1 一是) 。如此等等，在每一个分界面上 都能形成反射波和透射波，如图2 2 所示。 r 1 r 2 r os 图2 2 多层介质中波的反射 在具有多界面的介质中，由于各层介质的速度不同，波的传播不再是以直线形式传 播，而是以折线形式传播。一般来讲，上、下界面的反射波是独立的、彼此互不干扰、 依次的向上传播。若下层介质的波速大于上层介质的波速，还会在分界面产生滑行波， 1 3 第二章地震波传播的摹本理论 并在上层介质中产生折射波。由震源出发向外传播，没有遇到分界面直接到达接收点的 波叫做直达波。因此，地震波在多层介质中传播时，按照波在传播过程中传播路径的特 点，可以产生直达波、反射波、透射波、折射波等波形，如图2 3 所示。 激发点 地面 透射波 图2 3 多层介质中的各种地震波 1 4 长安大学硕- 上学位论文 3 1 射线追踪基本原理 第三章射线追踪正演 3 1 1 射线追踪理论 射线追踪是常用的一种正演模拟方法，射线追踪法主要理论基础是几何光学，在高 频近似条件下，地震波的主能量沿射线轨迹传播，它对模型的几何尺寸要求是大大超过 子波的波长。基于这种认识，运用惠更斯原理和费马原理来重建射线路径，并利用程函 方程来计算射线的旅行时。 射线追踪就是把在介质中连续传播的波离散成一条射线来研究，从震源出发将射线 分成若干段，逐段求出射线的轨迹和走时，依次累加，最后得出地下介质中射线分布和 地面上的时距曲线。 众所周知，在介质速度变化不太大，波长比较短的情况下，波的传播可以用波前方 程来描述： ( 詈) 2 + ( 詈) 2 + ( 警) 2 = 咒2 。m 式中n 为折射率，w 为波阵面( 波前) 波前方程( 3 1 ) 是一阶偏微分方程，它的解只是空间坐标的函数。具有非常简单的形 式， w ( x ，y ，z ) = c ( c 为常量) ( 3 2 ) 它表示了三维空间曲面簇，也就是波的传播可以视为从震源散开的一个个空间曲 面。每一时刻对应一个曲面，这个曲面就是该时刻的波阵面，波阵面的法线就是所谓的 射线，射线的方向就是波传播的方向。波强度的各种变化可以用射线的收缩与扩张来描 述，更重要的是到达波的时间可以由沿射线路经的积分得到。 从波前方程很容易看出波阵面的形状与介质的速度分布有关。在均匀介质中速度是 常数，相同时间间隔的波阵面之间对应点的距离相等，即对应点的曲率不变，因此法线 的方向也不会变，也就是说均匀介质中射线是直线。在非均匀介质中，速度分布是空间 坐标函数，在相同时间间隔内波阵面前进的距离不一样，波阵面的曲率发生变化，波阵 面的法线方向也跟着发生变化，也就是说非均匀介质中的射线是曲线。 如果已知速度分布和边界条件，可以求出波前方程的解析解或数值解。也可以把非 1 5 第三章射线追踪正演 均匀介质中弯曲的射线，用若干个圆弧来逼近。即把非均匀介质分成若干薄层，把每个 薄层的走时视为是线性变化的。所以，速度线性变化的介质，地震射线的轨迹是圆弧。 因此很容易分别求出射线在各个薄层的轨迹和走时，依次累加，最后得到射线在非均匀 介质中的射线分布和走时【4 1 。 3 1 2 程函方程和输运方程 对波动方程 窖g c 2 v ：p ( 3 3 ) o t 、 其中的p 为质点位移，c 为速度，我们作关于t 的f o u d e r 变换，则得到 一2 p 。c 2 v p 其中为角频率