（地球探测与信息技术专业论文）波动方程变网格步长有限差分数值模拟.pdf

波动方程变网格步长有限差分数值模拟 李胜军( 地球探测与信息技术) 指导教师：孙成禹副教授 摘要 有限差分算法是常用的正演模拟方法之一，其优点是包含的地震信 息丰富，且实现简单。传统的常网格差分算法在模拟含有低速地表或低 速高速夹层模型时，常存在不稳定或频散问题，为了得到较高精度的模 拟结果，常采用减小网格步长的方法。这样即浪费了计算机内存资源又 降低了计算效率。 为解决上述问题，本文在前人研究的基础上，发展了纵向变网格步 长差分算法，该算法不需要进行波场插值，实现简单，计算效率较高， 且相对常网格算法节省了大量的计算机内存资源。为验证算法的有效性， 本文对不同模型进行了试算，并从模拟精度，计算效率，内存需求三个 方面与常网格算法做了对比分析，结果表明纵向变网格步长算法优于常 网格算法，该方法可视为对传统的有限差分算法的有效改进。 对比表明，纵向变网格步长差分算法与完全匹配层吸收边界算法有 着类似的实现思想，两者之间可以建立良好的过渡关系，这为将来实现 变网格完全匹配层吸收边界提供了可能。文中对处理边界反射的透明边 界，旁轴近似吸收边界，完全匹配层吸收边界进行了定量的对比，并结 合理论模型进行了试算。结果表明完全匹配层边界吸收效果较好，实现 也比较简单。 关键词：数值频散，变网格，计算时间，边界条件 i i w a v e e q u a t i o nn u m e r i c a lm o d e l i n go n am e s ho f v a r y i n g s p a c i n g l is h e n g d u n ( g e o p h y s i c a lp r o s p e c t i n ga n di n f o r m a t i o nt e c h n o l o g y ) d i r e c t e db ya s s o - p r o f e s s o rs u nc h e n g - y u a b s t r a c t f i n i t ed i f f e r e n c ei so n eo ft h ep o p u l a rm e t h o & f o rs e i s m i cw a v e f i d d m o d e l i n g ，s i n c ei t s r e c o r dm a yc o n t a i na b u n d a n ti n f o r m a t i o na n dt h e a l g o r i t h mi se a s yt ob er e a l i z e d t h et r a d i t i o n a lf i n i t ed i f f e r e n c ea l g o r i t h m h a st h ep r o b l e mo fi n s t a b i l i t ya n d o rn u m e r i c a ld i s p e r s i o nw h e nd e a l 吣w i t h am o d e lw i t hl o w - s p e e ds u r f a c em e d i u mo rl o w h i g hs p e e di n t e r l a y e r t og e t s i m u l a t i o nr e s u l t sw i t h 】1 i g h e rp 陀c i s i o l l ，am e t h o do f d e c r e a s i n gg r i ds t e p si s a d o p t e du s u a l l y , b u tw h i c hm a k e s t h ec o m p u t i n gt i m el o n g e ra n dc o m p u t e r m e m o r yb ew a s t e d t os o l v et h ea b o v ep r o b l e m ，b a s e do nt h ef o u n d a t i o no ft h ep r e c e d i n g r e l a t e dr e s e a r c h , av a r y i n gg r i da l g o r i t h mi nt h el o n g i t u d i n a ld i r e c t i o ni s i n t r o d u c e da n di m p r o v e di nt h i sp a p e r , w h i c hi se a s yt ob er e a l i z e dw i t h o u t w a v e - f i e l di n t e r p o l a t i o n , c a ni m p r o v ec o m p u t a t i o ne f f i c i e n c y , a n ds a v e m e m o r yr e s o u r c e se f f i c i e n t l yc o m p a r e dw i t hi m d i t i o n a lo n e f o rv e i l 刖n g t h ev a l i d i t yo fa l g o r i t h m ，s e v e r a lp r a c t i c a lm o d e l sa r et e s t e d c o m p a r e dw i t h t h et r a d i t i o n a lf i n i t ed i f f e r e n c e a l g o r i t h mf r o mt h ea s p e c to fs i m u l a t i o n p r e c i s i o n , c o m p u t a t i o ne f f i c i e n c ya n dm e m o r yr e q u i r e m e n t ，t h er e s u l t s i i l d | c a t ct h a tt h en e wm e t h o di sm u c hb e t t e rt h a nt h ec o n v e n t i o n a lo n e a n d c a nd e c r e a s en u m e r i c a ld i s p e r s i o n e f f e c t i v e l y , w h i c hi sr e g a r d e da sa n 1 i i e f f e c t i v ei m p r o v e m e n to v e rt h et r a d i t i o n a lo n e c o m p a r e dw i t ht h ep e r f e c t l ym a t c h e dl a y e r ( p m l ) a b s o r b i n gb o u n d a r y c o n d i t i o na l g o r i t h m ，i ti sp r o v e dt h a tt h ev a r y i n gg r i do n ei nt h ev e r t i c a l l y d i r e c t i o nh a st h es i m i l a rr e a l i z a t i o nt h o u g h t s t h e r ei s g o o d t r a n s i t i o n r e l a t i o n s h i pb e t w e e nt h e m , w h i c ho 舒e r sap o s s i b i l i t y 协r e a l i z ev a r y i n gg r i d p e r f e c t l ym a t c h e dl a y e r ( p m l ) a b s o r b i n gb o u n d a r yc o n d i t i o ni nt h ef u t u r e t h r o u g hq u a n t i t a t i v ec o m p a r i s o na m o n gt h et z a n s p a r e n tb o u n d a r y , p a r a x i a l a p p r o x i m a t i o n sa b s o r b i n gb o u n d a r ya n dp m lo n e ，a n dc a l c u l a t i o nb y e x a m p l e s ，t h er e s u l ts h o w st h a tt h ep e r f e c t l ym a t c h e dl a y e ra b s o r b i n g l x , u n a a r yc o n d i t i o ni sb e t t e rt h a no t h e ra b s o r b i n gb o u n a , a r yc o n d i t i o n sa n di s e a s yt ob er e a l i z e d k e yw o r d s ：n u m e r i c a ld i s p e r s i o n ，v a r y i n gg d d ，c o m p u t i n gt i m e ，8 0 u n a a r y c o n d i t i o n i v 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论 文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得中国 石油大学或其它教育机构的学位或证书而使用过的材料。与我一同工作 的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了 谢意。 签名： 2 0 0 7 年岁月谚日 关于论文使用授权的说明 本人完全了解中国石油大学有关保留、使用学位论文的规定，即： 学校有权保留送交论文的复印件及电子版，允许论文被查阅和借阅；学 校可以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手 段保存论文。 ( 保密论文在解密后应遵守此规定) 学生签名：兰嫱 导师签名：= 三丕毖 2 0 0 7 焦 2 0 0 7 焦 譬只 ， 只 订日 二j ，日 中国石油大学( 华东) 硕士论文第1 章前言 第1 章前言 地震波场正演贯穿于地震学发展的始终，其中预测给定地质模型中 特定观测位置上地震响应的地震波传播模拟是正演的一个重要方面。在 地震采集、处理与解释过程中，正演模拟都是一个不可缺少的内容，许 多科研单位都将该部分内容列入目前的科研规划中，对其进行重点研究。 地震勘探的理论基础是波动理论，研究波动现象对于地震勘探有着 重要的意义。对于波动现象的深入认识，不仅可以指导地震勘探的实施， 而且还为地震数据处理提供了理论依据，同时地震模拟技术也是验证地 震数据处理效果的重要方法【l 】。 同其它物理现象一样，波动现象可以用偏微分方程来描述。一般而 言，这些偏微分方程很难用解析方法求解，尤其在较为复杂的非均匀介 质中，因此必须借助数值计算方法或其它渐近解。直接用数值计算法求 解波动方程的方法是对波动理论方程的渐近近似，一般精度较低。波动 理论方法基于最基本的波动理论方程，因此精度较高，波场全面。本文 的研究也就是基于波动理论。 1 1 有限差分算法的研究现状 有限差分是波动方程正演模拟和偏移处理使用最广泛的数值计算方 法。用有限差分法解决勘探地震学波动问题最早可追溯到2 0 世纪6 0 年 代。自从a l t e m 柚闭( 1 9 6 8 ) 将有限差分( f d ) 法引入到地震波模拟中以来， 差分方法就很快被用于解决各种勘探地震学的实际问题，并在应用中不 断发展。7 0 年代，c l a e r b o u t 对有限差分方法在勘探地震学领域的应用做 了比较大的推进，到7 0 年代中后期，差分算法已经牢固建立了它在勘探 地震学中不可替代的地位，发展也比较成熟，与差分相关的理论和技术 ( 如频散现象、吸收边界条件的处理) 也逐步发展、完善。如何进一步提 高模拟精度，减小数值频散，减小计算量，改善边界的吸收效果一直都 是有限差分法研究的核心。a l f o r d 3 1 等研究了有限差分的精度，指出了差 分网格必须明显小于最短波长的一半，如果网格足够精细，就能得出精 中国石油大学( 华东) 硕士论文第1 章前言 确的结果。v i d e u x 4 提出了稳定的二阶( 空间和时间) 弹性波有限差分格 式，它适用于任何泊松比的介质，l e v a n d e r 5 1 将v i r i e u x 的方法推广到空 间四阶、时间二阶的情况。c r a s e 6 】则发展了精度可达任意阶的高阶交错 网格法，但其计算量和内存需求与低阶有限差分法算法相比大幅度增加。 m a g l l i c r 【7 】等提出了最小网格有限差分法，它能压制非最小网格的人为现 象。在国内，周家纪和贺振华【8 】用大网格快速算法模拟地震波的传播， 其空间网格可以取得很大，达到每个最短波长只需3 个网格点，大大缩 短了计算时间，等等。然而，在许多情况下，这一问题并没有得到根本 解决，尤其是对低速介质中的弹性波模拟时，数值频散仍很严重。 地震勘探是一个复杂的动力学系统，在这个系统中发生着能量和信 息的转换。由于各种条件的限制，地震勘探所得到的地震记录( 数据) 只能是有限的。数据越多，对地下诸因素的反推就越准确。由于地震波 的传播是满足波动方程的，但波动方程只在特定的条件( 初、边值条件) 下解才唯一。因此，需给定初、边值条件才能确定出波场的情况。值得 注意的是人为边界条件的引入，由于边界条件的不同，对波场是有影响 的。c l a y t o n 和e n g q u i s t 9 1 ( 1 9 7 7 ) 根据旁轴近似理论导出了声波与弹性波 方程的吸收边界条件，r e y n o l d s 【l 川( 1 9 7 8 ) 发展了声波方程的透明边界条 件；王华忠、马在田【i i 】在最小二乘意义下，用共轭梯度法导出了优化求 解吸收边界条件方程中的系数，使得吸收边界条件方程尽可能与内部方 程相逼近，从而使得边界反射波被全部或部分吸收掉，得到最佳吸收系 数的吸收边界条件方程；近些年p m l ( p e r f e c t l ym a t c h e dl a y e r ) 边界条件在 地震勘探领域也得到推广，且效果不错。 传统的有限差分算法都是基于固定步长的矩形网格实现的。在很多 实例中，由于采用固定网格步长使得算法缺乏稳定性。如果这些方法被 用到垂向速度相对较高的介质中，模型的大部分区域可能被过采样；如 果将这种方法应用到有低速夹层或低速地表的介质中，网格不够密的话 会产生数值频散，分辨率降低等问题，为得到精确的模拟结果就必须对 整个模拟区域用较密的小网格步长进行采样模拟，这就造成了内存需求 的增大，同时也降低了计算效率。 2 中国石油大学( 华东) 硕士论文第1 章前言 针对传统有限差分算法的上述限制，可通过改变网格步长的方法加 以改善。变网格步长有限差分算法的基本思想是根据计算区域的速度将 模型划分为不同的区块，对速度高的区域用大网格步长，对速度低的区 域用小网格步长，对需要重点模拟的区域也用小网格步长，从而可实现 提高分辨率，节约计算机内存，提高计算效率的目的。因为网格步长可 根据物质的特性和模型的形状而改变，所以，与常网格差分算法相比， 这种方法对声波、弹性波模拟的灵活性、有效性都能有很大的提高。 1 2 变网格步长差分算法的提出 在地震波模拟中，有限差分法按模型的基本方程形式不同可分为位 移差分、位移应力差分以及速度应力差分等方法；以物理的定义点不 同而分普通网格差分和交错网格差分( 交错网格通过将物理量定义在同 一网格不同的半格点和格点以及参数平均，避免了在普通网格中采用中 心差分时所产生的误差) ；以网格的形状还可以分为均匀网格差分和不均 匀网格差分，其中不均匀网格差分又分为不均匀连续网格和不均匀不连 续网格【1 2 1 。 由于传统的有限差分算法采用了固定的网格步长，使得该方法在一 些实例中的有效性受到限制。从空间采样的角度考虑，为避免对小网格 步长对高速区域的过采样，很多学者做了大量的研究。m o c z o 【1 3 】( 1 9 8 9 ) ， 用网格连续变化的方法来来减小计算量，j a s t r a m 和b e h l e 1 4 1 ( 1 9 9 2 ) 提出 了二维声波方程针对某一深度变网格步长的算法；j a s t r a m 和 t e s s e m c r l l s ( 1 9 9 4 ) , 1 羽了一种类似的方法针对二维弹性波模型，提出了在 垂向上网格步长逐渐变化的算法；w a n g 1 6 】( 1 9 9 5 ) 也发展了二维声波、弹 性波变网格有限差分算法，该方法在网格边界用了双倍的步长。w a n g 和s c h u s t e r l l 7 1 ( 1 9 9 6 ) 又将j a s t r a m s 算法扩展到对不同模型时三维和轴对 称粘弹性波的有限差分模拟。f a l k i l 8 】等( 1 9 9 8 ) 把j a s t r a m 和b e h e l 的二维 方法发展到三维井间地震模拟。 在国内，张剑割旧1 ( 1 9 9 8 ) 基于应力、速度混合变量的弹性波动方程 及任意四边形网格差分算子，给出了交错计算应力及速度的非规则网格 3 中国石油大学( 华东) 硕士论文第1 章前言 弹性波应力速度差分法。该方法有较高的计算精度，所需的计算存储空 间也较少，计算效率相对较高。但是变网格步长算法在地震模拟上研究 较少，在热学领域，周建兴【2 0 r 2 1 ( 1 9 9 9 ，2 0 0 0 ) 针对有限差分问题中的均 匀网格技术存在的弊端，提出了变网格技术，给出了区分租、细网格的 假设与原则，并把它应用到了温度场的数值模拟。 前人的上述算法无论从内存需求，还是计算效率上都有了很大的改 进，同时使得有限差分算法对不同的模型的适应性也更加灵活。但是调 整网格步长时，存在一个中间过渡带，有的算法必须进行插值，如j a s t r a m 等【1 4 1 ( 1 9 9 2 ) 给出的变网格算法，在某深度处要进行波场插值，这样在过 渡带区域边界也就导致了稳定性问题的出现：有的在过渡带边界用平滑 函数进行平滑，使得差分实现起来变得困难。 本论文中，发展了一种新的对二维声波、弹性波模拟的变网格步长 有限差分算法。该方法与前面介绍的方法的不同是，步长在特定的边界 发生变化，但不需要进行插值计算。也不需要平滑函数进行平滑，使得 算法简单，且易于实现；同时对某些模型，由于网格的变化可节省大量 的计算时间和内存。 1 3 本文的研究内容 在导师的指导之下，本论文研究了均匀各向同性均匀介质条件下的 波动方程数值模拟的变网格算法。通过变网格差分算法解决了常网格差 分算法模拟含低速高速夹层、低速地表等特殊模型时存在的过采样、计 算效率低等问题，实现了保证较高模拟精度同时有较高计算效率，内存 需求不是很大的目的。本文就以下几个方面进行了研究： ( 1 ) 基础理论研究。对波动方程的不均匀化差分方法进行了介绍； 减小数值频散，提高模拟精度是差分算法面临的问题，本文对减小数值 频散的高阶差分算法进行了研究，并对空间、时间上不同精度对频散的 影响进行了分析。 ( 2 ) 差分算法对计算机内存要求高，且计算时间比较长，为解决上 述问题，文中完成了对声波方程、弹性波方程变网格差分算法的推导， 4 中国石油大学( 华东) 硕士论文第l 章前言 并从一维方程入手给出了详细的计算方法和离散公式。 ( 3 ) 设计低速地表、起伏地表、低速夹层等模型，分别用变网格差 分算法与常网格差分算法进行了试算。从模拟精度，内存资源的需求， 以及计算效率上对结果进行了详细的分析对比，结果显示，变网格差分 算法模拟精度较高，且节省计算机内存资源，计算效率得到显著提高。 “) 数值模拟时，加入了人为边界，所以必须对边界条件进行处理， 文中对完全匹配层( p m l ) 边界条件进行了分析研究，并定量的描述了反 射系数随匹配层的厚度的变化，该边界条件与入射角度没有关系，模型 模拟结果表明，完全匹配层边界条件效果较好，符合数值模拟的要求。 ( 5 ) 对完全匹配层吸收边界与变网格差分算法的关系进行了分析， 为下一步的研究工作点定了理论基础。 ( 6 ) 文中分别对透明边界，基于旁轴近似的吸收边界进行了详细的 分析，从定量的角度分析了反射系数随入射角的变化情况，并与p m l 边界条件进行对比。 ( 7 ) 通过算法验证和模型试算证实，文中提到的变网格差分算法不 会在变网格边界差生人为的反射；对速度变化剧烈的模型该算法模拟精 度较高，与相同精度的常网格差分算法相比，计算效率有显著提高，同 时节省了大量的计算机内存资源。 5 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 第2 章波动方程的有限差分数值解法 2 1 波动方程及其有限差分解法 在地震波传播的实际介质中，弹性参数可以看成是位置坐标的函数， 而密度一般变化不大，可视为一个常数。这样波的传播速度也将是位置 坐标的函数。在这种条件下，我们可以由弹性体的运动微分方程式导出 位移方程瞄l 。假设其中的体力项为零，可有 p 警= 誓+ 等+ 誓 p 雾= - d - + 誓+ 誓 ( 2 - 1 ) p 争= 誓+ 等+ 誓 铲吖2 l 堕o u j + 甜下骱搠惭一矧三维全波波动 睾替+ 考黜2 辩 曲 + 导l + 甜昙l 隹+ 剀 一叫 6 瑚b 塑砂丝砂 一一+ - 斗引垒渖。飞 u儿a一七 剀心 +、；叫 f s 一砂西苏 + _ 钰一融，础一渺、=，ril 廿 引两旦缸 一 十 塑酽 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 孝嚣絮隅。+ 丢 ( 老+ 警) + 昙 声( 笔十考) 峥叫 刁、失一股性，1 殳半面波的波前面平行于y 轴，此时波动问题与y 轴 无关，变为二维问题。即有= o ，此时，纵波( p 波) 和垂直极化波 ( s v 波) 耦合在一起，描述这两个波的两个互相关联的二阶偏微分方 程可以写成 p 害= 去 旯罡+ 尝) 伽割+ 鲁 瞎+ 老) 劬 p 字= 永罢+ 警) 伽卦甜瞎+ 老) 协 当密度p 和拉梅系数五和为常数时，方程变为均匀介质中的方程 ( 2 - 4 ) 。 p 等川+ 刎昙降+ 毫卜陪+ 是) q 4 砂 p 窘出伽，瓦a 商ra 2 u + 斜+ 降嚎 如果利用公式( 2 4 ) 所示的均匀化公式进行地震模拟，需要对速度变 化的地下各点再补充边界性条件，这样做起来很麻烦。通常情况下，在 进行地震模拟的时，都使用非均匀介质中的公式( 2 3 ) ，同时假设整个介 质的密度p 是一个常数，可以把上述方程( 2 3 ) 写成纵波和横波速度随空 间变化的函数，即九g ，z ) 和p g ，z ) 可以用v ，g ，z ) 和v ，x ，z ) 来代替。 现在，用有限差分去逼近( 2 3 ) 式。等号左边的时间微分可以用简单 7 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 的微分项的差分。在方程( 2 - 3 ) 中出现的有两类项。一类仅仅与一个空间 坐标变量工或z 有关的微分商，y j - - 类是对两个空间变量x 和z 有关的偏 微商。对于第一类的偏微商，倒如 昙限z ) 割 ( 2 - s ) 在任意网格点( m ，厅) ，用离散值v ；d ，行) 代替v ；g ，z ) 。规定v ；，甩) 是 中心在如缸，玎& ) ，边长为x 和z 的矩形上v ；g ，z ) 的平均值。式( 2 5 ) 项的差分逼近可以写成 生堂坠竺塑掣：逊：竺坐：! ! ! q d k ) 2 其中平均值v ；白+ l ，2 ) 和v ；0 一l 2 ) 用下式确定 v ；) = 塑萼型 对于第二类微商，例如 昙lv ；g ，：) 昙“g ，z ，r ) i = 丢k g ，z ，r ( 2 8 ) 其中为了方便，引入函数c g ，z ，f ) 表示( 2 8 ) 等号左边方括号中的项。这样 整个右边的项可用中心一阶差分逼近 知而必血业去警趔 ( 2 9 ) 瞄z z 对c g ，z ，f ) = v ；云“g ，z ，f ) 用一阶中心差分进行逼近 c ) 一v 胁) 亟型毫等趔 ( 2 1 0 ) 8 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 孙) 去岫，) l “志 v ； ，栉+ l h b + 1 ，”+ 1 ，) 一“一l ， + 1 ，) 1 ( 2 - 1 1 ) 一v ；b ，以一1 ) 【“+ 1 ，月一1 ，) 一“如一1 ，丹一1 ，为) 式( 2 - 3 ) t 9 的其它各项也可以用类似的方法进行处理，从而得到具有 二阶精度的弹性波有限差分格式。对声波方程可做类似的推导。 2 。2 波动方程差分解中的几个问题 2 2 1 震源函数 在波动方程的正演模拟过程中，震源的处理对模拟结果有重要的影 响。在实际地震波动过程中，震源机制是十分复杂的，一般在地震模拟 过程中，通常用点震源或线震源来表示。 文中采用指数衰减的正弦子波作为震源函数，主频为2 5 h z 。式( 2 1 2 ) 是所用的子波的函数表达式，图2 - 1 为子波波形。 w ( t ) = e x p ( - a t 2 ) s i l l ( 2 石厶f ) ( 2 1 2 ) 式中厶为子波的主频，t 2 为衰减系数。 1 0 0 o 舶 0 m 粤 鞲咖 _ 0 舶 - 0 瑚 0 1 k i加加加6 0 卸珊 时间1 1 1 s 图2 - 1 子波波形图 9 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 2 2 2 波动方程的高阶精度差分 这里所说的商阶差分，主要是指空间偏导数的高阶精度。因此要取 得高阶精度的有限差分，就要用到波函数在空间多个点上的值。基本思 想就是： 设函数l l ( x ) 连续，且具有2 k ( k + 1 ) 阶导数，将其展开成泰勒级数的表 达式。经整理后可以得到，具有2 n 阶精度的奇数m 阶导数的近似表达 式为 害：杰w。u(x十心)(2-13)dx ”幺 其中，砒为权因子，满足 w - ，= w i ，m = 0 ( 2 1 4 ) 般情况下，去1 2 阶精度就足以满足导数计算的要求。表2 1 给出 了当m = l ，n = l , - 6 及0 0 时，一阶导数对应与各阶精度情况下各采样点的 权系数值。其中，当精度为0 0 时，系数未能全部写下。 表2 - l 一阶导数对应与各阶精度的权系数值 同理，也可以得到函数l l ( x ) 具有2 n 阶精度的偶数m 阶导数的近似 表达式，具有与奇数阶导数同样的表达式 1 0 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 其中，为权因子，满足 岩= 薹甜g + 心) w _ l = ，= 0 ( 2 - 1 6 ) l 表2 2 给出了当m = l ，n - 1 “及0 0 时，二阶导数对应与各阶精度情 况下各采样点的权系数值。同上，当精度为o 。时，系数未能全部写下 表2 - 2 二阶导数对应与各阶精度的权系数值 2 2 3 波动方程有限差分解的稳定性 一个有实际意义的数值计算要求有限差分算法是稳定的，即有限差 分数值解与解析解之差必须是受限制的。稳定性是一个实用差分格式必 备的条件。差分格式有两类，一是隐式差分格式，二是显式差分格式。 一般来讲，隐式差分格式为无条件稳定，而显式差分格式是有条件稳定 的，并随差分精度的不同而有所不同。当差分精度的阶数提高时，差分 格式的稳定性变坏。研究得出，对于声波方程而言，2 疗阶差分精度的稳 定性条件为 佗一1 7 ) 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 其中谚”为表2 - 2 中精度为2 刀阶时第f 个系数，a x ，血分别为横向和 纵向的采样步长，f 为时间采样步长。当缸= a z = h 时，变为 f ( 2 - 1 8 ) 可见，随着差分精度2 珂的增大，时间采样要求越来越密，这样才能 满足稳定性的需要。也就是说，高阶的空间差分精度能够在不提高空间 采样率的前提下提高差分精度，其代价是必须提高时间采样率。 2 3 有限差分中频散问题的理论分析 在地震波的正演模拟中，所追求的目标是：模拟精度要高，计算速 度要快。微分方程数值模拟方法在求解波动方程时，或多或少都会产生 不期望的数值频散或称网格频散，导致了数值模拟结果分辨率降低】。 究其根源，是因为微分方程法是先对介质、波场函数进行了网格剖分， 利用离散化的微分方程去逼近微分波动方程，从而使得波动方程的系数 发生变化，使波传播的相速度变成了离散空间间隔的函数。因此，当每 一波长内空间采样点太少( 即空间网格太粗) 时，就会产生数值频散。这种 频散即不同于波动方程本身引起的物理频散，也不同于因波传播的速度、 频率和角度而引起的频散，它是微分方程法求解波动方程时所固有的本 质特征，无法避免的。为了消除这种数值频散，如何提高模拟精度是微 分方程方法所必须面对的关键问题之一。为了达到这一目标，许多人都 作了努力：在时间空间域有限差分方面，a l f o r d t 3 l 等( 1 9 7 6 ) 的二阶差分方 法到d a h l a i l l 剀( 1 9 8 6 ) 的高阶差分方法；从规则网格到交错网格( v i r i e u x ， 1 9 8 6 ；0 z d e n v a r & m c m e c h a n ，1 9 9 7 ；董良国【2 5 | 2 6 1 1 2 7 1 ，2 0 0 0 ：k e t i lh o k s t a d ， 2 0 0 3 ) ；在频率空间域有限差分方面，c h u l r l h y u n ( 1 9 9 6 ) 提出了旋转九点 差分格式对各向同性的p 波模拟：p r a r ( 1 9 9 5 ) 采用了2 5 点差分格式和优 化方法对弹性波正演模拟，减少了数值频散提高了地震波场的模拟精度。 还有的学者采用了求解守恒方程的通量校正传输方法( f c t ) 【2 s 1 【2 9 1 。 1 2 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 网格频散是相速度成为频率和网格大小的函数，从理论上讲，网格 频散并不成问题，十分精细的网格总能使频散变的很弱，但过分精细的 网格会使计算所需的存储量和计算时间非常大，而过粗的网格又会使计 算结果面目全非。 2 3 1 空间数值频散和时间数值频散 ( 1 ) 空间数值频散 为方便起见，剿扮析二维波动方程害+ 窘= 专窘的空间2 n 阶差分精度近似时的频散特征。 设平面波传播方向与x 轴的夹角为0 ，将平面谐波 u ( x ，z ，f ) = o x p f f ( , o t k x c o s o 一乜s i n o ) ( 2 1 9 ) 代入2 n 阶空间差分公式( 这里假定缸= a z ) ， 嚣= 麓瞧篙：篆篙，z ) ， 可以得到： 茜= 厅坛卅【c o s c o s o ) + c o s b i i l o ) _ z 】 ( 2 2 t 2 ) 其中，矿= 兰为地震波相速度，妒= 七缸= 2 万等为相位角。 图2 - 2 、2 3 、2 - 4 分别显示了0 = 0 。、0 = 2 2 5 。和0 = 4 5 。，空间差分 精度分别为2 、4 、6 、8 和1 0 时的空间数值频散随酬九的变化曲线。从 图中可以看出，空间差分引起的数值频散由三个因素决定，一是地震波 传播方向，二是空间差分精度，三是一个波长内离散点数目，它们和频 散之间的关系是： l 、随着地震波传播方向和离散坐标轴之问夹角的增大，频散降低。 传播方向与x 轴的夹角为0 的平面波，其空间数值频散与传播方向与x 轴 1 3 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 的夹角为9 0 。一0 的平面波相同。也就是说0 = 4 5 。时离散数值频散最小。 2 、无论地震波传播方向如何，随着空间差分精度2 n 的提高，上述 高阶差分解法产生的数值频散将会逐渐减小。因此，数值频散和差分精 度存在直接的关系，可以通过提高差分精度的办法来减小数值频散。一 般情况下，采用8 阶空间差分精度即可满足抑制数值频散的要求。 3 、对任意阶空间差分精度，随着一个波长内离散点数的增加，离散 网格频散越严重。 另外，从图2 - 3 和图2 - 4 可以看出，y z o ，也就是说，由于 时间离散造成的数值频散在波形上在正常到达时之前出现；( 2 ) 随着差 分精度2 m 的提高，上述高阶差分解法产生的数值频散会逐渐减小。 1 15 0 0 1 4 0 x 13 0 12 0 11 0 1 09 0 ， ， ， ，。 么夕 ，? 多 ， 盏兰 一：= 2 阶 4 阶 6 阶 8 阶 1 0 阶 0 0 0 50 1 00 1 5 0 2 0 0 2 5n 3 003 50 4 00 4 60 5 0 址i t 图2 - 7 不同时间差分精度数值频散随a t l t 变化曲线 从式( 2 2 2 ) 可以看出，时间差分引起的数值频散由两个因素决定，一 是差分精度，二是一个时间周期内离散点数目，如图2 7 所示。一般情 况下，采用四阶时间差分精度即可以满足抑制数值频散的要求，更高的 时间差分精度对解决数值频散问题效果不明现，反而会影响模拟效率。 通过对有限差分方法空间和时间离散数值频散的理论分析，说明了 1 8 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 影响频散的主要原因是空间和时间网格大小、差分格式和差分精度、介 质速度以及子波频率，证明了利用高阶差分法可以减小数值频散的结论。 2 3 2 网格步- k 比对数值频散的影响 类似上面的推导，设平面波传播方向与x 轴的夹角为0 ，将平面谐 波( 如式( 2 1 9 ) ) 代入2 n 阶空间差分公式( 这里假定缸= 地= 胁) ， 专窘= 糟窿篙二乏篙二4 删l p 2 s ， 可以得到： 嚣= 蠢倭w 枷训) m 如e ) 一2 ( 2 + 1 ) 】) y f 2 2 5 ) 其中，v ：旱为地震波相速度，k ：孕。令 膏 r = 谚) b s o 枷c o s o ) + c o s 0 肼2 h s i n o ) 一2 ( f 2 + 1 ) 】 ( 2 2 6 ) 从式( 2 2 5 ) 可以看出，当式( 2 2 6 ) 6 p 的r 取极小值时频散最小。由2 3 1 节的结论可知，对不同精度时频散随采用点的变化是一致的，我们以二 阶精度为例，做出不同横向与纵向采样步长比与频散最小时入射角的变 化曲线，如图( 2 8 ) 所示。 5 0 4 0 之一3 0 篓z o 1 0 0 t f 、_ 。咔一 051 01 52 02 5 d x ，d z 图2 - 8 采样步长比与频散最小时入射角的交化曲线 1 9 中国石油大学( 华东) 硕士论文第2 章波动方程的有限差分数值解法 从图2 8 可以看出，横向与纵向网格步长比越大，频散程度最小的 入射角越小。如当使用正方形网格( 横纵网格步长比为1 ：1 ) 时，波沿与垂 直方向呈4 5 。夹角传播时的频散程度最轻，而当网格的纵横比为1 ：2 时， 波沿与垂直方向呈2 6 5 6 5 。夹角传播时的频散程度最轻。考虑到在实际 地震勘探中，入射角一般不大于3 0 。，即地震波的传播路径偏向铅垂方 向。同时，地层速度在横向上的变化也远没有在纵向上那么剧烈。因此， 对待频散可以在两个方向区别对待，在纵向上加密网格更有利于压制数 值频散接近实际情况。 中国石油大学( 华东) 硕士论文第3 章波动方程变网格有限差分算法 第3 章波动方程变网格有限差分算法 3 1 引言 有限差分算法已经被广泛的应用于地震波在介质中传播的模拟。相 对其它模拟方法，有限差分算法可以模拟任意密度和速度模型的波场。 传统的有限差分算法用固定的网格步长对整个模型区域进行离散，为得 到较精确的模拟结果，整个区域都必须使用较小的采样间隔( a l f o r de ta 1 1 9 7 4 ) 3 1 ，这样降低了计算效率。每个波长上太少的采样又会导致数值 假象，又叫数值频散。a l f o r d 等( 1 9 7 4 ) 和k e l l y l 3 0 1 等( 1 9 7 6 ) 总结出为削弱 数值频散，在时间二阶精度，空间上二阶精度情况下每个波长上最少需 要1 1 个网格点，空间上四阶精度时每个波长上最少需要6 个采样点。 虽然常网格步长差分算法比较容易实现，但是他们对大部分模型而 言都不必要的增大了计算量。如对有浅层低速层沉积盆地模型地面运动 响应的模拟，由于这种地层能增大地面运动响应的振幅和延续时间，所 以这种模型显得比较重要( o l s e n e t a l ，1 9 9 5 ) 3 t 】。不像深层的高速层， 这些浅层低速层需要用小网格。当然，常网格算法必须用小网格离散整 个模型，这就导致增加不必要的代价，如内存需求的增大、计算效率的 降低 3 2 1 1 3 3 1 。 近年来，超小规模的非规则体非均质的洞、礁、孔、缝储集体的模 拟也提到议事日程上。常规的均匀网格的算法难以适应这些需求，为了 模拟小尺度( 米级、厘米级，或更小) 的地质介质，需要大幅度地减小 计算网格。如果对设计的全模型都按照小尺度介质确定网格，其所需的 存储空间和运算时间的增加十分巨大，即使大规模并行计算机也难以承 受。如果只在所研究的目的介质上加密网格，需要增加的存储空间和运 算时间就不会太多。 通过调整网格步长，也就是在任意层内每个波长上的采样点不再是 固定的，可以提高计算效率，减小内存需求。例如，有一两层地质模型， 2 l 中国石油大学( 华东) 硕士论文第3 章波动方程变网格有限差分算法 顶层速度是底层速度的一半，如果顶层的采样步长是底层采样网格步长 的一倍，这样的网格分布是最为理想的。这对大区域模型的模拟来说， 将节省大量的计算时间和内存需求。 很多学者对变网格差分算法做了大量的研究，其中，j a s t r a m 和 b e h l e 1 4 1 ( 1 9 9 2 ) 提, q 4 , 了二维声波方程针对某一深度变网格步长的算法，受 到了广泛关注。j a s t r a m 和t e s s e m e r is l ( 1 9 9 4 ) 用了一种类似的方法针对二 维弹性波模型，提出了在垂向上网格步长逐渐变化的算法。上述算法无 论从内存需求，还是计算时间上都有了很大的改进，同时使得有限差分 算法对不同的模型的适应性也更加灵活。但是在调整网格步长时，存在 一个中间过渡带，有的算法必须进行插值，这样在过渡带区域边界将会 导致了稳定性问题的出现；有的算法对过渡带边界用平滑函数进行平滑， 使得差分实现起来变得困难1 3 4 ”l 。 本文，针对上述限制，发展了一种新的对二维声波、弹性波模拟的 变网格步长有限差分算法。该方法与前面介绍的方法的不同是，步长在 特定的边界发生变化，但不需要进行插值计算，使得算法简单，且易于 实现；同时对某些模型，实用性较高，且由于网格步长的变化，可大幅 度减少计算时间，节省内存。 3 2 声波方程变网格有限差分算法 变网格步长差分算法在步长由大到小( 或由小到大) 调整的时候，在 数值计算时会出现一个特殊的区域，我们定义该区域为过渡带( 四阶精度 时的过渡带如图3 - 1 中阴影部分所示，下侧过渡带与上侧类似) 。为计算 过渡带的波场值，在过渡带内要用不同的网格步长进行重叠覆盖，直到 全部过渡到目标网格。j a s t r a m 和b e h l e ( 1 9 9 2 ) 提出的纵向变网格步长差 分算法( 网格分布如图3 - 1 所示) ，在网格由大到小过渡时要通过一种简 单的插值，实现网格步长的调整。即，在深度z 0 上侧和z i 下侧有一部 分没有波场值，均需要进行插值，这在算法的实现就增加了难度。另外， 据w a n g 和s c h u s t e r ( 1 9 9 6 ) 1 拘i 算例验证及分析，通过这种插值来计算过渡 中国石油大学( 华东) 硕士论文第3 章波动方程变网格有限差分算法 带的波场值，在网格大小变化梯度不大，模拟精度不高时，也就是过渡 带比较窄时，这种变网格算法是比较好的，但是，如果网格调整较大， 模拟精度较高，将会在过渡边界处产生一个弱的人为边界反射，同时插 值边界也可能会存在稳定性问题。 辽 菠 撩- f 一 yl ，v 1 1 i 1 i i yi f i i i i fvi 月l1ii b书书 j lq，l，1l lilll y yi i ff li1 l 半 i l l i i l i l 一 i ，i1 lli1i门 llli l ili iiiii ll ii ii ll 1 过 渡 带 图3 - 1 变网格步长插值算法示意图 为了解决上述问题，避免在过渡带边界插值时产生人为的边界反射 以及不稳定问题，从易于实现的角度出发，我们总结前人的算法，对上 述变网格算法做了进一步简化，取消了在过渡带的插值。为简单起见，