（地球探测与信息技术专业论文）高阶差分数值模拟方法研究.pdf

高阶差分数值模拟方法研究 刘庆敏( 地球探测与信息技术) 指争老师：李振春教授 龋要 目前，石油勘探】f 向复杂构造、起伏地表等地区发展。地质构造的复 杂性，要求有更好更精确的方法束模拟真实的地下波场的传播规律。为了 提高地震波数值模拟的精度，本文在前人的基础上，采用了高阶差分方法， 进行了数值模拟和波场特征分析。 文中首先简单的推导r 地震波场导数的高阶差分近似，对声波方程进 行规剐闻格和交错网格的高阶有限差分数值模拟，结果表明，交错网捺高 阶差分具有很高的模拟精度，计算效率很高。另外，起伏地表地震波数值 模拟是波场模拟中的一个难点，采用了声波方程三角网格有限差分法，可 较好的描述起伏地表和弯曲界面。在此基础上，实现了弹性波以及粘声波， 粘弹性波的交错网格高阶差分数值模拟，推出一阶速度应力弹性波波场 分离方程，成功地分离出p 波波场和s 波波场，从模拟中得到的波场快照 和炮记录中分析在各种复杂介质内部的反射、透射、绕射、敞射以及能量 | 勺衰减等运动学和动力学的各种细节特征。 稳定性，数值频散和边界是数值模拟中的关键性问题。通过比较各种 边界条件的优劣，采用绍合边界条件束消除边界反射的影i 响，用通量校正 传输方法来抑制数值频散，提高计算精度。对声波方程和弹性波方程交错 网格高阶差分法稳定性问题进彳亍了简单地分析。 通过对各种介质模錾i 试算表明，本文所述方法稳定性好，模拟精度高。 关键词：高阶差分，交锚网格，声波方程，弹性波方程，粘弹性波方程 m e t h o ds t u d yo fh i g h - o r d e rf i n i t ed i f f e r e n c e n u m e r i c a lm o d e l i n g l i uq i n g m i 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 yp r o f e s s o fl iz h e n - c h u n a b s t r a c t a tp r e s e n t , o i la n d - g a s 函i 豌商6 n i st u r n i n gt ot h ea r c 弱w i t hi r r e g u l a r t o p o g r a p h y a n d c o m p l e xg e o l o g i cs t r u c t u r e s t h ec o m p l e xg e o l o g i c a l s t r u c t u r e ，r e q u i r eb e t t e ra n dm o r ea c c u r a t em e t h o dt os i m u l a t et h er e a l u n d e r g r o u n dw a v ef i e l dp r o p a g a t i o n ，t oi n c r e a s es e i s m i cw a v en u m e r i c a l s i m u l a t i o na c c u r a c y , t h e h i g h o r d e rf i n i t e d i f f e r e n c em e t h o dw h i c hi st o p r o c e s sn u m e r i c a lm o d e l i n ga n dw a v ef i e l da n a l y s i si su s e d f i r s t l y ，t h i sp a p e rd e r i v es i m p l yh i g ho r d e rd i f f e r e n t i a la p p r o x i m a t i o ni nt h e s p a c ea n dt i m ed e r i v a t i v e so fs e i s m i cw a v ef i e l d t h er e g u l a rg r i dh i g h o r d e r f i n i t ed i f f e r e n c ea n ds t a g g e r e dg r i dh i g h o r d e rd i f f e r e n c en u m e r i c a ls i m u l a t i o n i su s e dt ot h ea c o u s t i cw a v ee q u a t i o n ，t h em o d e l i n gr e s u l t si n d i c a t et h a t h i g h o r d e rs t a g g e r e dg r df i n i t eb u ta l s oh a sg o o ds i m u l a t i o nr e s u l t sa n dh i i g h e f f i c i e n c y i na d d i t i o n ，a na n n o y i n gp r o b l e mi nt h es e i s m i cw a v em o d e l i n gi s t h es e i s m i ct o p o g r a p h y t h ef i n i t e d i f f e r e n c es e i s m i cm o d e l i n gm e t h o db a s e d o nt r i a n g u l a rg r i d sm e t h o dc a nb eu s e dt od e s c r i b em o r ea c c u r a t e l yi r r e g u l a r s u r f a c ea n dc u r v e db o u n d a r i e s o nt h i sb a s i s ，w ea c h i e v et h en u m e r i c a l m o d e l i n go f2 一d e l a s t i cw a v ee q u a t i o n ，v i s c o a c o u s t i cw a v ee q u a t i o na n d v i s c o e l a s t i cw a v ee q u a t i o n ( k e l v i nm o d e l ) u s i n gt h eh i g h o r d e rs t a g g e r e dg r i d f i n i t e d i f f e r e n c e ，d e r i v e t h ef i r s t o r d e r v e l o c i t y s t r e s se l a s t i c w a v e f i e l d d e c o m p o s i t i o ne q u a t i o n , a n ds u c c e s s f u l l yo b t a i nt o t a l l ys e p a r a t e dw a v e f i e l d o fp u r ep - w a v ea n ds - w a v ef r o mt h eh y b r i dw a v e f i e l d ，f r o mr e c e i v e dt h e w a v e - f i e l d ss n a pa n ds h o tr e c o r d si nn u m e r i c a lm o d e l i n g ，t oa n a l y z er e f l e c t i o n ， t r a n s m i s s i o n ，d i f f r a c t i o n ，s c a t t e r i n ga n da t t e n u a t i o no fe n e r g yk i n e m a t i ca n d d y n a m i c c h a r a c t e r i s t i c si nc o m p l e xm e d i u m ， s t a b i l i t y , n u m e r i c a ld i s p e r s i o na n db o u n d a r yc o n d i t i o na r et h es e v e r a lk e y p r o b l e m si nt h es e i s m i cm o d e l i n g b yc o m p a r i n gt h em e r i t sa n dd e m e r i t so f b o u n d a r yc o n d i t i o n s ，t h ec o m b i n i n gb o u n d a r yc o n d i t i o ni su s e dt oe l i m i n a t e t h eb o u n d a r yr e f l e c t i o nw e l l ，t h ef l u x - c o r r e c t e dt r a n s p o r tm e t h o di su s e dt o a v o i db o t h s p a c e - g r i dd i s p e r s i o n a n dt i m e g r i dd i s p e r s i o na n dh i g h - o r d e r d i f f e r e n c ec a nt o , s o m ee x t e n tr e d u c et h en u m e r i c a ld i s p e r s i o na n di m p r o v e a c c u r a c y w ea n a l y z es i m p l ys t a b i l i t yo fs t a g g e r e dg r i dh i g h o r d e r f i n i t e d i f f e r e n c em e t h o di nt h ea c o u s t i ca n de l a s t i cw a v ee q u a t i o n t h et e s tr e s u l t so fv a r i o u sm o d e lp r o v e dt h a tt h i sm e t h o da p p l i e di s c h a r a c t e r i s t i co fg o o ds t a b i l i t ya n dh i g hs i m u l a t i o np r e c i s i o n k e yw o r d s ：h i g h - o r d e rd i f f e r e n c e ，s t a g g e r e d 一鲥d ，a c o u s t i cw a v ee q u a t i o n ， e l a s t i cw a v ee q u a t i o n ，v i s c o e l a s t i cw a v ee q u a t i o n 独创性声明 本人声明所呈交的论文是我个人在导师的指导下进行的研究工作及 取得的研究成果。尽我所知，除了文中特别加以杯注和致谢的地方外，论 文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得石油大 学或其它教育机构的学位或证书而使用过的材料。与我一同工作的同志对 本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。 签名：壶l 丞丛知叼年 ；月z 矿日 关于论文使用授权的说明 本人完全了解石油大学有关保留、使用学位论文的规定，即：学校有 权保留送交论文的复印件及电子版，允许论文被查阅和借阅；学校可以公 布论文的全部或部分内容，可以采用影印、缩印或其它复制手段保存论文。 ( 保密论文在解密后应遵守此规定) 学生签名： 导师签名： 到灸故 至盗蠡 如u 7 年 z ，海 岁月 j ，- 月 29 日 l 夕只 中国石油人学( 华尔) 硕士论文第1 苹前言 第1 章前言 地震数值模拟是地震勘探和地震学的重要基础。地震数值模拟就是在 假定地下介质结构模型和相应物理参数己知的情况下，模拟研究地震波在 地下各种介质中的传播规律，并计算在地面或地下各观测点所观测到的数 值地震记录的一种地震模拟方法。这种地震数值模拟方法己在地震勘探和 天然地震领域中得到广泛的应用，它不但在石油、天然气、重金属和非金 属等矿产资源及工程和环境地球物理中得到普遍的应用，而且在地震灾害 预测、地震区带划分以及地壳构造和地球内部结构研究中，也得到相当广 泛的应用。 地震数值模拟在地震勘探和地震学各工作阶段中都有重要的作用。在 地震数据采集设计中，地震数值模拟可用于野外观测系统的设计和评估， 并进行地震观测系统的优化。在地震数掘处理中，地震数值模拟可以检验 各种反演方法的正确性。在地震数掘处理结果的解释中，地震数值模拟又 可以对地震解释结果的于确性进行检验。 由于实际工作中所模拟的介质不同，所用的模拟方程也不一样。根据 模拟方程的不同，波动方程数值模拟主要有：声波模拟、弹性波模拟、粘 弹性波模拟以及裂隙和孔隙弹性模拟等。由于可以用射线理论、积分方程、 微分方程来描述地震波的传播，模拟方法也相应地有射线追踪法、积分方 程数值求解方法以及微分方程数值求解方法。 射线追踪方法通过求解程函方程计算地震波旅行时，通过求解传播方 程计算地震波振幅。该方法以高频近似为前提，适合于物性缓变模型中地 震波传播模拟。模型简单时该方法具有计算速度快的突出优点，正因为如 此，它在地震成像、旅行时层析等方面得到广泛应用。也正是高频近似， 该方法不适合物性参数变化较大模型中地震波的传播模拟。 积分方程数值求解地震波数值模拟方法是基于惠更斯原理而得到的 一种波场计算方法，它又可以分为体积分方法和边界积分方法。该方法的 半解析特征，使其在成像，反演理论研究和公式推导方面具有得天独厚的 优势。由于涉及g r e e n 函数的计算，该方法一般适合于模拟具有特定边界 中国石油人学( 华东) 硕士论文第l 章前言 地质体产生的地震波，而要求该地质体周围为均匀介质。因此，该方法的 适应范围受到严格限制。 微分方程方法使对计算区域网格化，通过数值求解描述地震波传播的 微分方程来模拟波的传播。就目前看来，该方法对模型没有任何限制，在 地震波模拟中使用最为广泛，主要问题是计算量比较大，对计算机内存要 求较高；其中，有限差分法( f d ) 、有限元法( f e ) 以及傅立叶变换法( p s ) 是这类模拟方法中使用较多的方法。近年来还出现界于有限差分法和有限 元法之间的有限体方法( f v ) ，在理论上应该具有有限元法网格剖分的灵活 性，又具有有限差分计算快速的特点，但在简单的矩形网格情况下，该方 法完全退化为有限差分法。上述方法具有各自的优缺点： 有限元法( f e ) 是基于变分原理和剖分插值，考虑的是分段近似，比较 适宜于模拟任意地质体形态，可以任意三角形逼近地层界面，保证复杂地 层形态模拟的逼真性：但算法复杂，计算速度慢，一般对网格要求三角剖 分，基函数是分段线性函数，不具有正交性。算予也是空间局部的，空间 分辨率高，但是频率域中分辨率低，占用内存和运算量均较大。傅立叶变 换法( k o s l o f f a n db a y s a l ，1 9 8 2 ：r e s h e f e ta 1 ，1 9 8 8 ) 是一种逼近空间微 分的方法【1 2 i ，基于空间域中的求导，相当于频率域中的乘积运算，利用 傅里叶变换将波场函数表示为傅里叶级数的展开形式，将波动方程在时间 一波数域或频率域中求解，精度高，占用内存小；但是由于傅里叶变换是 基于整个时间域或空间域的，改变空间中的某一点的值，就会改变频率域 中的所有值，因此每一点的微分结果都要受到计算域中其它点的影响并且 存在众所周知的g i b b s 效应。故不适合复杂模型中地震波模拟，边界处理 也比较困难。有限差分法( k e l l ye ta 1 ，1 9 7 6 ) 使用方便灵活【”，应用比 较广，目前国际勘探地球物理界著名得二维m a r m o u s i 模型以及s a l t 模型 的地震数据合成均是利用f d 方法。 a l t e r m a n 和k a r a l ( 1 9 6 8 ) 首先将有限差分法应用于层状介质弹性波 传播的数值模拟q b 4 1 。此后，b o o r e ( 1 9 7 2 ) 又将有限差分法用于非均匀 介质地震波传播的模拟 5 1 。a l f o r d 等( 1 9 7 4 ) 研究了声波方程有限差分法 2 中国i i 油人学( 华尔j 硕士论文第l 章前言 模拟的精确性1 6j 。d a l a i n ( 1 9 8 6 ) 和m u f t ie ta t ( 1 9 9 0 ) 讨论了用高阶差分 方程解决声波正演模拟问题1 7 4 】。在高阶方程情况下，网格距可以取得很 大，但计算精度并不比二阶差分方程小网格距时低，而且有效的提高了计 算精度。网格距的增大可以大大降低对计算机内存的要求、缩短计算时间。 此后，b a y l i s s ( 1 9 8 6 ) 、l e v a n d e r ( 1 9 8 8 ) 采用了四阶空间有限差分法弹 性波传播的地震记录一4 。此后，人们将高阶差分与交错网格相结合( c r a s e e ，1 9 9 0 ) 1 1 1 。 为了进一步模拟地震波在非完全弹性的实际地层中的传播，c a r e i o n e 等( 1 9 8 8 ) 提出了粘滞声波在地层中传播的模拟方法f 1 2 1 ；t a t e z e r 等( 1 9 9 0 ) 进行了线性粘弹性介质中地震波传播的方法研究3 】；r o b e r t s s o n 等( 1 9 9 4 ) 给出了粘弹性波有限差分模拟方法【i ”。 c a r c i o n e 和h e l l e ( 1 9 9 9 ) 提出了孔隙粘弹性介质中地震波传播的交 错网格有限差分模拟方法；p i t a r k a ( 1 9 9 9 ) 给出了三维各向同性介质中 弹性波的矩形非规则交错网格有限差分模拟方法【l6 】；董良国等( 2 0 0 2 ) 给 出了一阶弹性波方程交错网格高阶差分解法，并且给出了稳定性条件 1 17 - 1 8 ；裴正林( 2 0 0 4 ) 运用交错网格一阶空问导数的任意偶数阶精度展开 式和相应差分系数计算式以及一阶双曲型应力一速度弹性波方程交错网 格任意偶数阶精度差分格式来求解方程【1 9 1 。 交错网格高阶差分法具有很高的模拟效果，计算效率很高。我们用交 错网格高阶差分实现了声波，弹性波以及粘声波，粘弹性波的数值模拟， 从得出的波场快照和炮记录中分析在各种复杂介质内部的反射、透射、绕 射、散射以及能量的衰减等运动学和动力学的各种细节特征。 在地震波传播理论研究中，波动方程用于无限介质空间，通常假设地 球介质为半无限空间介质，微分方程法数值模拟是模拟地震波在这半无限 空间介质中的传播过程，用计算机模拟时，介质的范围必须是有限的，即 人为地限定地球介质的计算区域，由此产生了人工边界，当地震波通过这 种人工边界时，就会产生边界反射，它严重干扰波场，必须消除，因此边 界处理成为了数值模拟的一个关键问题。典型的边界条件按原理来分为两 中国石油入学( 华东) 硕士论文第1 章前言 类：某种单程波构成的吸收边界条件和波动沿波的传播方向逐渐衰减的衰 减边界条件。自边界条件问世以来，许多学者从不同角度提出来多种构造 边界条件的方法：( 1 ) 波动方程分解法，如r e y n o l d s ( 1 9 7 8 ) 的透明边界 条件【2 0 】；( 2 ) 旁轴近似方法，利用不同精度近似的单程波方程作为吸收边 界( c l a y t o n 和e n g q u i s t ，1 9 7 7 ；h i g d o n ，1 9 9 1 ) 1 2 1 也l ；董良国( 1 9 9 9 ) 利用特征分析方法将1 9 7 9 年h e d s t r o m 提出的一维情况下无边界反射概念 推广至三维各向异性介质弹性波的数值模拟中，得到t i 介质中的吸收边 界条件【2 3 1 。( 3 ) 阻尼衰减法，在靠近边界的一定宽度区域设为衰减带，使 得向边界传播的波场在此区域内逐渐得到衰减，降低人为反射，直至没有 明显的反射波回到计算区域( c e i ：i a n ，1 9 8 5 ；k o s l o f f ，1 9 8 6 ；s o c h a c k i ， 1 9 8 7 ) 【2 5 - 2 6 1 ：( 4 ) 最佳匹配层法( p m l - - t h ep e r f e c t l ym a t c h e dl a y e r ) ，一种 比较新的边界构造方法，是在模拟电磁波时被提出的，主要是在边界处加 一个匹配层，在匹配层只能够通过一个阻尼因子来衰减边界反射 ( b e r e n g e r ，1 9 9 4 ) 【2 刀；c o l l i n o 和t s o g k a ( 2 0 0 1 ) 把这种方法成功地运 用于弹性波的波场模拟1 2 。 针对各种边界条件的优缺点，把组合边界( 透射边界和匹配层边界， 特征分析法边界条件和匹配层相结合) 应用于高阶差分数值模拟中，很好 的处理了边界反射问题。 数值频散是有限差分数值模拟的又一个关键问题。地震波在传播过程 中存在物理耗散与物理频散现象：物理耗散是指波的振幅因物理阻尼作用 而衰减的现象；物理频散是由物理介质的原因，波传播相速度随波数发生 变化的现象。用差分方程逼近微分方程时引入了误差项，有时这些误差项 使计算结果振幅值衰减和相速度发生变化，其作用相当于物理耗散和频 散，这种虚假的物理效应称作数值耗散和数值频散。数值频散实质上是一 种因离散化求解波动方程而产生的伪波动，这种频散不同于波动方程本身 引起的物理频散，而是差分方程所固有的本质特征。 为了消除波场模拟中的数值频散问题，许多学者在这方面作了大量的 研究工作。从不同的角度对有限差分方程的数值频散进行了分析，并给出 4 中国石油人学( 华尔，硕十论文第l 章前言 了相应的解决办法。 a l f o r d ( 1 9 7 4 ) 和d a b l a i n ( 1 9 8 6 ) 对声波二阶空间差分的数值频散进 行了分析，指出网格大小和地震波传播方向是影响频散的两个因素。蔡其 新( 2 0 0 3 ) 等人关于有限差分数值模拟的最小频散算法及其应用中，提出 优化算法的主要内容包括高阶有限差分、优化差分参数和f c t 技术( 通量 校j 下传输方法) 2 9 - 3 0 l 。董良国、李培明( 2 0 0 4 ) 地震波传播数值模拟中的 频散问题中分析了影响地震波数值计算中网格频散的各种因素【3 ”，从理论 上以及模拟实例上证明了高阶差分( 特别是交错网格高阶差分) 是提高波 动方程数值计算精度、降低数值频散的有效方法。吴国忱、王华忠( 2 0 0 5 ) 也详细地讨论了波场模拟中的数值频散分析与校正策略【j 。 对高阶差分法声波模拟和交错网格弹性破模拟而言，影响数值频散的 三个因素是地震波传播方向、差分精度和一个波长内离散点数，对交错网 格弹性波模拟而言还包括介质的泊松比。我们用通量校币传输方法( f c t ) 有效的压制了在粗网格情况下差分计算产生的数值频散。其过程主要包括 3 步：有限差分计算、平滑方程的解和反漫射处理。 实现了不同波动方程交错网格高阶差分数值模拟的程序设计，对简单 模型和复杂模型进行了试算，验证方法和算法的正确性和有效性，结果较 为理想。 中国石油大学( 华东) 硕士论文第2 章地震波场的空间导数和时间导数的数值逼近 “( x o m 妒蛾) + 羔j = l 半) 十d ( 舻，( 2 - 1 ) 叫肌缸号争+ 善l - i 盟2 t + 孝l 罱斗2 t + l ( 2 i 川瓯) + o ( 扩卅1 ) 】 乏 缸 智+ 1 ) ! 。” 血掣l2 缸毗棚妒峨一) 】 ( 2 - ，) ，2 甓1 麓委l 善l - i m = l 酱矿州陪4 ， o p l ，2 i、1 ，。、 + 壹鲤娑喾a i n u ( 2 l ：i ) ( 护d ( a x 2 ( i , + i ) + i ) 一 智( 2 l + 1 ) ! ”7 、 2 1 ( ) 2 3 ( 上) 3 6 e ( 2 5 ) 中国tj 油人学( 华尔j 硕七论文第2 章地皮波场的空间导数和时l 司导数的数值逼近 解系数方程( 2 5 ) 百j - 得 2 - 1 ) 。 口m = ( 一1 ) ”兀i 2 一l i = 1 l m 2 m 兀( 聊2 一i 2 ) 兀( f 2 一埘2 ) - l i f f i m + i 由此公式我们可以得到一阶导数不同差分精度的差分权系数( 见表 表2 - 1 一阶导数对应于各阶精度的权系数值 2 l q 口2吒a 4 吼口6 200 050 0 0 0 0 e l oo o66 6 6 6 7 e 一1- 83 3 3 3 3 e 2 00 075 0 0 e i15 0 0 0 0e _ i 8o0 080 0 0 0 0 e 一1o0 0 0 0 0 e 138 0 9 5 2 e - 235 7 1 4 3 e 一3 i o 00 08 ”3 3 3 e l一23 8 0 9 5 e 159 5 2 3 8e | 2，99 2 0 6 3 e 379 3 6 5 1 e 4 l200 085 7 1 4 3 e 1- 26 7 8 5 7 e 179 3 6 5 1e | 2一i7 8 5 7 le - 225 9 7 4 0 e 一3一】8 0 3 7 5 e 一4 中心差分近似的截断误差系数为 铲赤萎舻“k m 中心差分近似的极限，即l _ c o 时，有 ：盟芝，铲0 于是，有 掣i ：a , u ( x o + m a x ) - u ( x o + m a x ) ( 2 - 6 ) 舐l 。m = l m 其中，一阶导数的中心差分算子长度为2 三 二阶导数的各阶精度差分公式【3 4 】 在此仅对关于x 的空间微商进行讨论，并假设差商具有的截断误差为 o f 血2 1 ，上是大于1 的数。 “o + 缸) = “( z ) + 办o _ ! u a x + 五1 万6 7 2 u 【- 血j 2 + 夏1 丽0 3 u 1 觚j - 3 + + 面1 豇瓦d z 万u ( 缸) 儿+ 7 中国石油大学( 华东) 硕士论文第2 章地震波场的空间导数和时间导数的数值逼近 “x - 缸) = “( x ) 一宝血+ 刍雾( 血) 2 一圭雾( 血) 3 两1 萨o z u ( 血) “+ 坐型蒯= 去雾( 缸卜去雾( 缸) 4 + + 面1 萨扩2 l ul 缸) 2 l + 。( 血2 。) ( 2 7 1 ) ”( x + 2 a x ) = “( x ) + 宝( 2 血) + 去雾( 2 缸) 2 + 上31 旦o p 兰f t 2 缸) 3 + ” + 上2 l ! 丝o x 2 t ( 2 缸) “+ “( x - 2 a x ) = “( x ) 一赛( 2 缸) + 刍雾( 2 血) 2 一上31 丛8 x 3 r 、2 缸) 3 + ”- + 。3 v u ( 2 a x ) “ 尘型型业型= 三! 丝d x 2 ( 2 缸h 1c74u。、(2 2o x 4o x 2 缸卜 ( 2 - 7 2 、 、7 1 ”7 f 2 - 2 、 +点2l磐(2缸)“+d(舻)8 x “、 7 、， “o + 上血) = “( x ) + 象( 三血) + 去雾( 缸) 2 1 31 巩c g _ 兰u 。( 三缸) 3 + + 面l 萨扩2 l ul # 缸) “+ “卜一血) = “( x ) 一象( l 缸) + 去雾( 上缸) 2 一圭雾( 三缸) + 去窘( 三缸) “+ 一- u(x+lax)-2u广(x)+一u(x-lax)=圭璺(缸)2+击鲁(缸)4(2-7-l)2,3x 4 , y x21 “7 1 “7 + 一2 上l _ _ ! 出w ：u ，、t a x ) 2 + o ( a x 2 ) 令 8 中国石油大学( 华东) 硕士论文第2 章地震波场的空间导数和时间导数的数值逼近 z ：巫些止攀丝亟型 厶= 亟型爿趔 五：坐兰堕垒2 垃巫型竺! q = 否d 。2 u ( a x , 1 2。d 4 u 4 ，吼= 萨0 2 l u ( 缸) 2 l 由( 2 7 一1 ) ，( 2 - 7 2 ) ，( 2 - 7 一l ) ，结合上述规定，可得如下方程组： 111 ， 五q + 雨0 2 十。+ 面吼2 l 2 22 42 “ 面q + 酉口：+ “+ 面吼。z ( 2 8 ) 譬”譬吒+ f 1 2 + a f ) | 表2 - 2 规则网格二阶导数对应于各阶精度的权系数值 2l 口0 d l口2q d d 口5d 6 22 0 0 0 0 0lo 0 0 0 0 425 0 0 0 0l3 3 3 3 383 3 3 ”e - 2 627 2 2 2 2l5 0 0 0 015 0 0 0 0 一1 l l 】1 1 1 e 一2 8 28 4 7 2 2 i6 0 0 0 0- 2o 0 0 0 0 e 一1 25 3 9 6 8 e - 2- 17 8 5 7 ie - 3 1 029 2 7 2 2l6 6 6 6 723 8 0 9 5 b l39 6 8 2 5 e 一249 6 0 3 2e - 331 7 4 6 0 e - 4 1 2- 29 8 2 7 8l7 】4 2 9- 26 7 8 5 7 e - l52 9 1 0 le - 289 2 8 5 7 e - 310 3 8 9 6 e 一3- 6 0 1 2 5 l e - 5 求解此线性方程组，即可得到口。，a 2 ，a l ，我们仅需q ，即得鲁缸：， 因此存在下式 ：雾确“( x ) + 妻 “x + m a x ) + “( x - m a x ) + 。( 舻) ( 2 9 ) 由此，我们得到规则网格二阶导数各阶精度的权系数值( 见表2 2 ) 。 混合偏导数的各阶差分格式 9 中国自油大学( 华尔) 硕士论文第2 章地震波场的空间导数和时间导数的数值逼近 果沿另一个方向( 如z ) 求取偏导数得到。若函数“( j ，z ) 的某阶混合偏导 数连续，则该导数的结果与求导顺序无关。以二阶混合偏导数为例，可以 写成 磊0 2 , = 耋喜叩烈h 础计尬) ) ( 2 - l 。) 这里，d ，吼均为一阶导数对应的权系数值( 同表2 - 1 ) ，显然满足 q 吼= o 2 1 2 交错网格任意2 l 阶精度有限差分系数计算公式 在交错网格技术中，变量的导数是在相应的变量网格点之间的半程上 计算的【3 5 】。为此，我们用下式计算一阶空间导数。 设“o ) ：a 2 l + 1 阶导数，则“( x ) 在x = x 0 塑岩血处2 上+ l 阶泰勒展 开式为 毗t 2 m - 1 咖等垡娑。俨冲， 由于交错网格一阶导数2 l 阶精度差分近似式可表示为 血掣= 圭m = 1 沁+ 等血卜卜莩缸 ) + 。c ， ( 2 一1 2 ) 将上述个方程带入、化简，有 舻c 护妻c z 训慨协妻雾气焉孚矿c ， 其中，待定系数由以下方程确定 1 0 中国t i 油人学( 华东) 硕十论文第2 章地宸玻场的空间导数和时间导数的数值逼近 13 13 3 2 三1 ( 2 l 1 ) 3 口l 口2 l 0 ( 2 l 一1 ) 2 一1jl , a ，ji o 表2 3 交错网格一阶导数对应丁各阶精度的权系数值 ( 2 1 3 ) 2 d la 2 吒d dd 5a 6 210 0 0 0 0 4i 1 2 5 0 0- 4 1 6 6 6 7 e - 2 61 1 7 1 8 76 5 1 0 4 2 e - 24 6 8 7 5 0 聃 81 1 9 6 2 9 79 7 5 2 6 e 295 7 0 3 le 369 7 5 4 5 e 4 1 0l2 1 l2 489 7 2 17 e 一2 13 8 4 2 8 e 217 6 5 6 6 e 一31l8 6 8 0 e 4 1 2l2 2 1 3 496 9 3 l5 e 2l7 4 4 7 7 e 2。29 6 7 2 9 e 一335 9 0 0 5 e 4- 21 8 4 7 8 e 5 由系数方程( 2 1 3 ) 计算得 l = 1 时，蟊= 1 ； l = 2 时，d = 1 12 5 ，a 2 = 一0 0 4 1 6 6 6 6 7 l 2 时，有 ( 2 m 一1 ) 兀e ( 2 m 由此我们可得到交错网格不同差分精度的差分权系数值( 见表2 3 ) 。 2 22 m 阶精度时间导数中心差分算法 2 2 1 声波方程中时间导数差分精度的近似 a 2 。 二阶双曲型声波方程中波场的二阶时间导数的高阶有限差分近 扰 似是通过波动方程将时间导数转换成对空间导数( d a b l a i n ，1 9 8 6 ) 。利用 泰勒展开，可以得到2 m 阶精度的时间差分近似。 , q 2 ， 首先讨论时间导数芸的四阶差商的推导。 中国4 i 油大学( 华东) 硕士论文第2 章地震波场的空间导数雨i 时间导数的数值逼近 争一a 1 t 2fu ( t 他) 砌( f ) + 啦f ) 】- ，邪a t 删2 m0 扩2 u ( t ) ( 2 - 1 4 ) 在进行正演模拟或偏移成像过程中，要避免差分方程所涉及的时间层 越多，所需的内存也就越大的问题出现。局部速度变化很缓或不变时，利 用声波方程把对时间的高阶微分转嫁到空间微分上去。因此，有 雾= 雾( 挚 = 雾妒c 碧+ 窘， s ， a 2 a 2 卅ia 2j a 2 叫v 、反21 岔2 、- 叫 以四阶时间差分精度为例 如胡= 知d 一嘶刊+ 陋，2 ( 窘+ 爹 + 三4 ( 窘+ 爹) 6 ， + 三4 嚣+ 出4 ) ( 2 - 1 6 ) 式是推导用于正演模拟的高阶差分方程的起始方程。它在时间 方向上截断误差为d ( f 4 ) ，而空间微商的差商的截断误差根据需要而定： 另外，可以根据需要来组合不同阶次的差分格式。 2 2 2 一阶弹性波方程2 m 阶精度的时间导数差分算法 二维各向同性介质一阶应力一速度弹性波方程可表示为： 票：a o ( 2 - 1 7 ) 其中，q = ( v 。，v ，仃。，盯：，盯。) a = o o ( a + 2 u ) l ， 他 肚， 0 土工，0上： pp 00 1 ：上， pp 肛000 ( 旯+ 2 0 t 000 u l ：0 00 1 2 中国j i 油入学( 华尔) 硕十论文第2 章地震波场的空间导数和时间导数的数值逼近 式中，l 。，l ：分别表不对x ，z 的一阶微分算于，q 表不远厦、压力i 司 量，a 表示算子矩阵，五、u 、p 表示各向同性介质l a m e 弹性常数和密 度。( 2 一1 7 ) 式时i q 上2 m 阶精度的显式差分格式为 q o + 了a t ) = q ( f 一筹) + 2 刍面i 1 面了a t ) “孑a 2 m 石- iq + 。( 出2 ”) ( 2 1 8 ) 式中，q = ( v ，v ：，盯。，盯。，盯。) 7 式( 2 一1 8 ) 中的高阶时间导数可以通过方程( 2 一1 3 ) 转移为对空间的导数 来实现，即 0 2 m 。- ：_ 1 q ：a2 m - , q 式中，当m = 1 时，为时问的二阶差分精度。当m = 2 时，则由( 2 1 8 ) 式可得 “，+ 争咄，一争，鲁+ 等等 陪 式中， 盟：一1f 监+ 监1 o t p 出 o z 矿c 。3 v x = 吉知坪，昙c 吉等+ 万ii o c r x = ， +丑毫(土一00xz+_1一ocr,p o x po z ) 】+ 土p 轧o z 鲁( 三p 等 o zo zo x + 土譬) + i 0 【一1 了0 0 x z + - 1 孚0 0 - ) 】 po z o xo呶po z 姒，+ 争味r 一争r 鲁+ 等等 弘2 。， 式中， o v ：一l ( o ( t x ：一十监1 生里塑塑盔堂! 兰垄! 堡主堡苎丝! 兰丝壅壅堑塑至塑量塑塑堕塑量墼盟塑笪丝望 万0 3 v z = 石1 瓦3w 瓦0l 万1 百0 0 = + 7 w 吒- z 帅去( 古百0 0 r z + 1 p0 出0 = ) 1 十土毫睢去01aofxx十一10a习,ppp o z 出出出 + ( 五+ 2 ) 呈( 土孥+ 一1 孚0 c r ) 】 o z p o x p o z 盯。( t + a t ，= 啪r 等+ 等等 式中， 百0 0 x 娟俐警+ 哮 等= ( 2 + 2 , u ) 否0 。石1 瓦0 ( ( 枷肋警+ a + 去毫警+ 声警) 】+ 且瓦0l - 1 瓦0 峥i o v + + = l 出a ( 兄誓o xm 伽) 口出 c 】z 哪圳咆+ 垃等+ 等等 式中， 等= 壕+ c + 2 a ) 老- 争= a 瓦0t 石1 夏0c c 五+ 2 ，警+ 五警，+ 古鲁c 警+ 誓) 】 + ( a + 2 力瓦0 。j l 瓦o 叫_ o v , o x + 警)成口锻 +丢晏(z誓m坪iovzozo x ) 】 d 哪圳砜础等+ 等等 1 4 ( 2 2 1 ) ( 2 2 2 ) ( 2 2 3 ) 土兰坚旦坠笠三! 兰堡竺鉴堡兰墅! j ! 型生塾塑塑兰囹量墼型堕塑量垫堑垫堕堕丝 式中， 等= 豢+ 哮 争2 瓦0t _ 万1 夏3 ( 警+ 警) + 万1 瓦0 豢堋+ z 肋善) 】 + c 3 万1 瓦oc c 五+ 劲，警+ 五誓，+ 吉尝c 豢+ 等刀 a 4 i 琏3 ( 2 - 1 9 卜( 2 2 3 ) 可得弹性波方程的四阶时间差分精度近似： k o + 争咄一争詈c 等十争+ 等z _ - pb 等o x 讹，+ m 。，舞坪一，急急帆0 耐3 0 x z ，但之4 匕c “等，2 匕c 卜等，+ 詈c 等+ 警，+ 丢r 。妥u a - 十c q ， 蚝乜，急地川急：袅坞，等，q 乏5 啦删础k 豢+ c i ，挈+ 筹【c l ，2 紊( 3 3 v 龟+ 凹也 2 4 m 一 蝎川意心，2 蝎一萨。) 急m 鳓 啪础) 吲吵豢蚝争+ 筹軎诋z + 戚出 2 4 ”m 3 铴坞悉地向蚝”训意坞2 万c 8 3 v 】 ”出) _ 啪) + 帆鼍+ c 4 t 挚+ 丽a t 3 2 等+ ( c 4 。2 + g 沈2 4 疗疗一 ，急托舳c i 3 。3 3 一。- c j 嘉坞。：羁3 2 3 q 乏8 5 中国打油人学( 华东) 硕士论文 第3 章一维各向同性介质卢波方程数值模拟 第3 章二维各向同性介质声波方程数值模拟 3 1 均匀介质声波方程规则网格高阶有限差分数值解 均匀各向同性介质二维声波方程可表示为 了1 萨a 2 u = 窘+ 虿9 2 z ( 3 1 ) v 2 甜2叙2 出2 、。1 式中，“( x ，z ，r ) 为声波波场，v 。( x ，z ) 为声波速度场。时间导数采用二 阶中心差分、空间导数为2 n 阶差分精度的二维声波方程的高阶差分格式 为 归u k 杈2 t 坐a x jl + 辩吖m ”) p ：， + 吉( 斟卜+ 喜碱一ko 3 2 非均匀介质中声波方程交错网格高阶有限差分数值解 非均匀各向同性介质中二维声波方程得一阶应力速度方程形式 詈一c 等+ 誓：! 罢( 3 - 3 ) 西p 缸 盟：三塑 a td 貌 其中，v ，v ：是质点速度，“法线应力，p 是密度，v 。是纵波速度。 我们采用交错网格进行削分，相应的声波波场分量和弹性参数的空间 位置见图3 - i 。下面给出2 n 阶空间差分精度、二阶时间差分精度交错网 格高阶有限差分格式，设u ! ；“2 ，p 毛，q 0 。：分别是应力“，速度u ， v ：的离散值。则方程( 3 - 3 ) 的差分格式为： 1 6 中国干油人学( 华尔) 硕士论又 第3 章一维各向同性介质卢波方稃数值模拟 ，。1 ，2 = ：u ! j 1 ，2 a t p 。v r _ _ 五： 薹( ：” 尸：c z 一一，：，f ! c ：一一t ，：，】) ( 3 4 ) 一i a t p v p 2t 善nc r p 。k 。m 】 卷 善掣嘲k - 1 1 ：2 + ( 2 n - 1 ) 1 2 , 锻k - 。1 1 2 七。m j ) 卷 善n 掣k - 1 1 。2m 。飞k - 。1 1 2 2 - ( 2 n - i ) 2 】 7 2 ，p v 。2 v ，p 一 v ：，p 。1 ( 3 5 ) ( 3 - 6 ) 图3 -