（地球探测与信息技术专业论文）基于高阶交错网格的有限差分地震波场数值模拟.pdf

摘要 地震波场数值模拟是研究复杂地区地震资料采集处理和解释的有效辅助手 段，地震波场数值模拟的主要方法包括两大类，即几何射线法和波动方程法。有 限差分法是波动方程中很重要的一种方法，它的主要优点是计算速度快，占用内 存小；缺点是精度低，仅适用于相对较简单的地质模型。 本文在前人的基础上，从波动方程的建立出发，推导了位移方程的中心差分 格式和一阶速度应力方程的任意阶空间和时间精度的交错网格差分格式。用有 限差分法模拟地震波在半无限介质中的传播问题，由于受计算机存储空间、计算 速度等的限制，需要用有界区域代替无界区域，即引入人工边界。人工边界的反 射使波场受到严重的破坏，因此一种有效的吸收边界的方法，成为地震波正演中 一个非常重要的研究问题。解决边界问题的方法很多，论文主要讨论了透射边界 和完全匹配层边界，同时还进行了稳定性条件、频散问题的分析。用f o r t r a n 语 言实现了对两种方程的各向同性介质的模拟，并且对于一阶速度一应力方程还实 现了横向各向同性介质的模拟。 在以上的理论基础上，通过模型试算及效果验证可以看出：基于一阶速度一 应力方程的交错网格差分格式精度高，可以适用于较为复杂介质模型的正演模 拟，得到的正演模拟结果可以较清晰的反映出地震波场中波传播的运动学和动力 学特征。由此可以证明，高阶交错网格有限差分具有高效、精确、实用的特点， 具有广阔的应用前景和实际意义。 关键词：弹性波动方程，有限差分，交错网格，吸收边界条件 a bs t r a c t s e i s m i cw a v ef i l e dn u m e r i c a lm o d e lc a l c u l a t i o ni sa l la v a i l a b l em e t h o df o r s e i s m i ci n f o r m a t i o na 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 no fac o m p l e xa r e a i t h a st w om a j o rs p e c i e s ，o n ei sg e o m e t r yr a y sm e t h o da n dt h eo t h e ri sw a v ee q u a t i o n m e t h o d f i n i t ed i f f e r e n t i a lm e t h o di sa ni m p o r t a n to n eo ft h ew a v ee q u a t i o ns o l v i n g m e t h o d s i t sa d v a n t a g e sa r ef a s ts p e e di nc a l c u l a t i o na n dl e s sm e m o r yo c c u p a t i o n i t s d i s a d v a n t a g e s i sl o wp r e c i s i o n i to n l yf i t ss i m p l eg e o p h y s i c s b a s e do nt h ep r e v i o u sr e s e a r c h ，s t a r t i n gf r o mt h ee s t a b l i s h m e n to fw a v ee q u a t i o n , id e r i v et h ec e n t e rf i n i t e - d i f f e r e n c eo fd i s p l a c e m e n te q u a t i o na n dt h ea n yo r d e rs p a c e a n da c c u r a t et i m ef o r m u l a t i o no fo n e - o r d e rs t r e s s - v e l o c i t ye q u a t i o n i no r d e rt od e a l w i t hl a r g ec o m p u t e rs t o r a g ea n dl o n gt i m eo fc a l c u l a t i o ni nt h ef i n i t e d i f f e r e n c e m o d e l i n go fs e i s m i cw a v e ，w eu s ef i n i t ed o m a i n st os u b s t i t u t ef o ri n f i n i t ed o m a i n s ， w h i c hi st os a ya na r t i f i c i a lb o u n d a r yn e e dt ob ei n t r o d u c e d r e f l e c t i o nw a v ef r o m a r t i f i c i a lb o u n d a r yi ss os t r o n gt h a ti tc a nd e s t r o ys e i s m i cw a v ef i e l d ，t h e r e f o r et h e w a yo fe f f e c t i v e l ya b s o r b i n gr e f l e c t i o nw a v ef r o ma r t i f i c i a lb o u n d a r yi s av e r y i m p o r t a n tr e s e a r c hs u b j e c ti ns e i s m i cw a v en u m e r i c a ls i m u l a t i o n t h e r ea r em a n y k i n d so fs o l v e n t sf o rt h ea r t i f i c i a lc o n d i t i o n ，w ef o c u s e do nt h ef o l l o w i n gt w om e t h o d s i nt h ep a p e r , t h et r a n s m i t t i n gb o u n d a r yc o n d i t i o na n dt h ep e r f e c t l ym a t c h e dl a y e r a t t h es a m e ，t h ep a p e ra n a l y z ei t ss t a b i l i t ya n dd i s p e r s i o np r o p e r t i e s t h e ni m p l e m e n t t w od i f f e r e n te q u a t i o n sm o d e l i n gi ni s t o t r o p i cm e d i u ma n dm o d e l i n go n e o r d e r s t r e s s v e l o c i t ye q u a t i o ni nt r a n s v e r s e l yi s t o t r o p i cm e d i u mw i t h f o r t r a nl a n g u a g e b a s e do nt h et h e o r ya b o v e ，f r o mt h ec o m p u t a t i o no fm o d e la n dt h ee f f e c t v a l i d a t i o nw ec a nk n o wt h a ts t a g g e r e d g r i df i n i t e d i f f e r e n c es c h e m eb a s e do nt h e o n e o r d e rs t r e s s v e l o c i t ye q u a t i o nh a sah i g ha c c u r a c gi tc a nb eu s e dt of o r w a r d s i m u l a t i o nf o rm o r ec o m p l e xm e d i u mm o d e l t h et e s tr e s u l tc a l lr e f l e c tt h ec h a r a c t e r s o fd y n a m i c sa n dk i n e m a t i c so fw a v e si ns e i s m i cw a v e - f i e l d t h et e s t r e s u l t s d e m o n s t r a t et h a tt h eh i g h - o r d e rs t a g g e r e d 一班df i n i t e d i f f e r e n c es c h e m ec a nb e e f f i c i e n ta n ds u f f i c i e n t l ya c c u r a t e i tm a ya l s ob ea ne f f e c t i v ea n dr o b u s tt o o li nw i d e a r e a si n v o l v i n gt h en u m e r i c a lc o m p u t a t i o no fp a r t i a ld i f f e r e n t i a le q u a t i o n s k e y w o r d s ：e l a s t i cw a v ee q u a t i o n ；f i n i t e d i f f e r e n c e ；s t a g g e r e d - g r i d ； a b s o r b i n gb o u n d a r yc o n d i t i o n 论文独创性声明 本人声明：本人所呈交的学位论文是在导师的指导下，独立进行研究工 作所取得的成果。除论文中已经注明引用的内容外，对论文的研究做出重 要贡献的个人和集体，均已在文中以明确方式标明。本论文中不包含任何 未加明确注明的其他个人或集体已经公开发表的成果。 本声明的法律责任由本人承担。 论文作者签名：正致 砑年朔万日 论文知识产权权属声明 本人在导师指导下所完成的论文及相关的职务作品，知识产权归属学 校。学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权 利。本人离校后发表或使用学位论文或与该论文直接相关的学术论文或成 果时，署名单位仍然为长安大学。 ( 保密的论文在解密后应遵守此规定) 论文作者签名：丘茎 导师签名：多缮 硼年y , 9 巧日 加年s 其鞠e l 长安大学硕士学位论文 1 1 研究的目的和意义 第一章绪论弟一早殖下匕 地震波场数值模拟是地震勘探和地震学的重要基础，进行数值模拟和偏移处 理不但可以检验和验证地震资料处理的效果，而且还对地震资料的解释有着非常 重要的指导意义。地震波场数值模拟就是对特定的地质、地球物理问题作适当的 简化，从而形成一个简单的数学模型，模拟研究地震波在地下各种介质中的传播 规律，并计算在地面或地下各观测点所观测到的数值地震记录的一种地震模拟方 法。它是研究各种地震地质条件下构造、物性和岩性等各种地质因素与地震波相 应特征之间的一门技术，它对人们认识地震波的传播规律，解释实际地震资料等 均具有重要的理论和实际意义。 地震数值模拟方法已在地震勘探和天然地震领域中得到广泛的应用。它不但 在石油、天然气、煤、金属和非金属等矿产资源及工程和环境地球物理中得到普 遍的应用，而且在地震灾害预测、地震区带划分以及地壳构造和地球内部结构研 究中，也得到相当广泛的应用。地震数值模拟在地震勘探和地震学各个工作阶段 中都有非常重要的作用。在地震数据采集设计中，地震数值模拟能够获得成像构 造关键部位的反射记录，同时能够保证采集到高质量的地震数据，从而达到缩短 勘探周期，提高经济效益的目的。在地震资料处理中，地震数值模拟可以检验各 种反演方法的正确性。在地震数据处理结果的解释中，地震波数值模拟又可以对 地震解释结果的正确性进行检验。 地震波场的数值模拟作为联结以反射波为主要课题的反射地震学和以岩层 为主要研究课题的钻井、测井、地质和油藏工程等学科的纽带，数值模拟技术的 应用几乎贯穿于地震资料采集、处理、解释、和油气藏开发工程的各个环节。随 着石油天然气工业的发展对地震勘探提出越来越复杂的勘探目标和越来越高的 勘探精度要求，我们必须在提高对复杂勘探目标中地震波传播规律认识的基础上 采取一些特殊的、有针对性的勘探方法和技术，从而发展新的勘探技术方法和技 术。要达到这一目标，我们必须深入对地震波场计算方法的研究，提高对地震波 传播规律的认识，从而对于解决现代油气勘探、开发工作中所面临的各种棘手问 题等具有非常重要的意义。 第一章绪论 1 2 有限差分的发展历史与现状 地震波场的数值模拟技术是在已知地下介质结构和参数的情况下，利用理论 计算的方法研究地震波在地下介质中的传播规律，合成地震记录的一种技术。随 着地震勘探技术的发展，数值模拟成为贯穿地震数据采集、处理和解释全过程的 一种重要方法，在确定观测系统的合理性，检验处理和解释的正确性等方面有着 越来越广泛的应用。地震波动方程的数值模拟方法主要有【l l = 有限差分法、有限 元法、反射率法、傅立叶伪谱法等，但这些方法都各具优缺点。傅立叶变换法由 于利用空间的全部信息对波场函数进行三角函数插值，所以能更加精确的模拟地 震波的传播规律。同时，利用快速傅立叶变换( f f t ) 进行计算，还可以提高运 算效率。该方法的主要优点是精度高，占用内存小，缺点是计算速度较慢。波动 方程有限元法的主要优点是适宜于模拟任意地质体形态，可以任意三角形逼近地 层界面，保证复杂地层形态模拟的逼真性，缺点是占用内存和运算量均较大。有 限差分法是对波动方程中的微分用差分近似代替，把波动理论上的微分形式变成 适于数值计算的差分形式，得到差分方程。它的主要优点是计算速度快，占用内 存小，缺点是精度相对比较低，适合相对较简单的地质模型。 有限差分法数值模拟技术开始于上世纪7 0 年代初，a l t e r m a n 等人( 1 9 6 8 ) 作了开创性的工作【2 】，使用显式有限差分格式获得了层状介质二阶弹性波方程的 离散数值解。a l t e r m a n 等人实际上得到的是均匀介质弹性波数值解，只在内界 面运用了应力和位移连续的内边界条件，使得波能通过弹性界面传播，对于结 构复杂和不规则的岩性层面，必须使用适应非均匀介质模型的方法，即自动满足 内界面处应力和位移连续的有限差分格式，由于其开创性的工作，差分方法很快 被用于解决各种勘探地震学的实际问题，并在应用中不断得到发展。b o o r e ( 1 9 7 2 ) 提出了非均匀介质二阶弹性波有限差分方法【3 】，a l f o r d 等( 1 9 7 4 ) 研究了声波方程 有限差分法模拟的精确性【4 1 ，k e l l y 等( 1 9 7 6 ) 改进发展了这一方法【5 】o 除了规 则网格外，一种较为先进的交错网格最早由m a d a r i a g a ，r ( 1 9 7 6 ) 提出【6 1 。v i r i e u x ( 1 9 8 4 ) 发展的交错网格方法使用一阶的速度应力方程 7 1 ，在不同的时间步上使 用不同的网格节点，分别计算应力和速度的传播，其差分精度为o ( t 2 + a x 2 ) ， 在不增加计算工作量和存储空间的前提下，和常规网格相比局部精度提高了4 倍，收敛速度也较快。由于仅使用一阶差分，解的光滑性较好，但它的稳定性条 2 长安大学硕士学位论文 件是最严格的。l e v a n d e r ( 1 9 8 8 ) 这种差分网格的精度提高到d ( 垃2 + a x 4 ) ，但精 度仍有待于提高【引。在解决计算精度和运算效率这一矛盾方面，d a b l a i n ( 1 9 8 6 ) 在求解声波方程时提出了高阶差分法【9 】。c r a s e ( 1 9 9 0 ) 将这一方法运用到求解二 阶弹性波方程中【1 0 】，这一方法的最大的优点是在保证计算精度的基础上，采用合 适的时间和空间差分阶数，从而采用较大的空间步长，提高计算效率。董良国 ( 2 0 0 0 ) 采用将速度( 应力) 对时间的奇数阶高阶导数转化为应力( 速度) 对空间的 导数【1 1 】，在不增加所需要内存量的前提下，将交错网格和高阶差分法有机结合， 求解横向同性( t i ) 介质一阶速度应力弹性波方程。另外还有很多学者如王秀明 ( 2 0 0 3 ) 1 2 】，杨旭明、裴正林( 2 0 0 7 ) 1 3 】，等将这一格式运用到波场模拟中，揭示了 波在地下传播的一些特性。 在地震波传播理论的研究中，波动方程用于无限介质空间，我们通常假设地 球介质为半无限空间介质，微分方程法数值模拟就是模拟地震波在这半无限空间 介质中的传播过程，但是用计算机进行模拟时，介质的范围必须是有限的，即人 为地限定地球介质的计算区域，由此产生了人工边界，当地震波通过这种人工边 界时，就会产生边界反射，它严重干扰波场，得不到很好的正演模拟结果，必须 消除，因此边界处理成为了数值模拟的一个关键问题。典型的边界条件按原理来 分为两类：某种单程波构成的吸收边界条件和波动沿波的传播方向逐渐衰减的衰 减边界条件。自边界条件问世以来，许多学者从不同角度提出来多种构造边界条 件的方法：( 1 ) 波动方程分解法，如r e y n o l d s ( 1 9 7 8 ) 的透明边界条件【1 4 】。( 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 5 1 ；董良匡1 ( 1 9 9 9 ) 利用特征分析方法将1 9 7 9 年h e d s t r o m 提出的维 情况下无边界反射概念推广至三维各向异性介质弹性波的数值模拟中，得到t i 介质中的吸收边界条件【l6 1 。( 3 ) 阻尼衰减法，在靠近边界的一定宽度区域设为衰 减带，使得向边界传播的波场在此区域内逐渐得到衰减，降低人为反射，直至没 有明显的反射波回到计算区域( c e r j a n ，1 9 8 5 ) b 7 1 。( 4 ) 多次透射边界法，廖振鹏等 模拟波传播的物理过程，提出的一种离散的吸收边界法( 2 0 0 ) b s 】。( 5 ) 最佳匹配层 法( p m l - t h ep e r f e c t l ym a t c h e dl a y e r ) ，一种比较新的边界构造方法，是在模拟电磁 波时被提出的，主要是在边界处加一个匹配层，在匹配层只能够通过一个阻尼因 子来衰减边界反射。c o l l i n o 和t s o g k a ( 2 0 0 1 ) 把这种方法成功地运用于弹性波的 第一章绪论 波场模拟【1 9 】。 数值频散是有限差分数值模拟的又一个非常关键问题。地震波在传播过程 中存在物理耗散和与物理频散现象：物理耗散是指波的振幅因物理阻尼作用而衰 减的现象；物理频散是由物理介质的原因，波传播相速度随波数发生变化的现象。 用差分方程逼近微分方程时引入了误差项，有时这些误差项使计算结果振幅值衰 减和相速度发生变化，其作用相当于物理耗散和频散，这种虚假的物理效应称作 数值耗散和数值频散。数值频散实质上是一种因离散化求解波动方程而产生的伪 波动，这种频散不同于波动方程本身引起的物理频散，而是差分方程所固有的本 质特征。 为了消除波场模拟中的数值频散问题，许多学者在这方面作了大量的研究 工作，从不同的角度对有限差分方程的数值频散进行了分析，并给出了相应的解 决办法。a f f o r d ( 1 9 7 4 ) 1 4 1 和d a b l a i n ( 1 9 8 6 ) 9 1 对声波二阶空间差分的数值频散进行 了详细的分析，指出网格大小和地震波传播方向是影响频散的两个重要因素。董 良国、李培i j ) 1 ( 2 0 0 4 ) 在地震波传播数值模拟中的频散问题中分析了影响地震波数 值计算中网格频散的各种因素【2 0 1 ，从理论上以及模拟实例上证明了高阶差分( 特 别是交错网格高阶差分) 是提高波动方程数值计算精度、降低数值频散的有效方 法。吴国忱、王华忠( 2 0 0 5 ) 也详细的讨论了波场模拟中的数值频散分析与校正策 略【2 1 】。 稳定性的问题是有限差分数值模拟方法的另一个重要的问题。所谓稳定性 就是指所有时间步的计算中，数值解的有界性。有限差分数值模拟的稳定性是由 多个因素决定的。由于计算过程中数值参数选择不合理，可能产生无物理意义的 按指数增大的数值计算结果，造成模拟结果网格频散严重，影响对问题的分析， 严重时会造成溢出而使计算无法进行。对于不同类型的波动方程稳定性条件是不 同的，而且随着模拟空间维数的提高稳定性条件变得更加严格。另外边界条件也 可引起有限差分数值模拟的不稳定。董良国( 2 0 0 0 ) 对一阶弹性波方程交错网格高 阶差分解法的稳定性做了详细的讨论【2 2 1 ，这种交错网格高阶差分解法和传统的精 度为d f 2 + 缸4 ) 的解法相比，随着差分精度的提高，其稳定性要求也有所提高， 即没有必要为提高差分精度而加密差分网格影响计算效率，因此高阶交错网格差 分法是一种即精确又高效的弹性波方程数值解法。 4 长安大学硕士学位论文 1 3 本文所做的工作 本文从地震波动方程出发，根据前人对位移波动方程的二阶中心差分方法 格式和对一阶速度一应力方程的高阶交错网格任意阶时间和空间差分格式的推 导，用f o r t r a n 语言进行编程，实现了对位移方程用中心差分格式对三层水平模 型介质的模拟，得到了每个时刻的波场快照和单炮记录，对一阶速度一应力方程 采用高阶交错网格法，在各向同性介质中进行了三种模型的模拟：( 1 ) 断层模型， 分别采用纵波震源和全波震源进行模拟，并且进行了对比。另外，还用平面波入 射法，对其进行了模拟，可以得到不同时刻的波场快照和单炮记录。( 2 ) 在三层 水平模型中采用纵波震源与用位移方程的中心差分法得到的结果相对比，不但得 到了波场快照和单炮记录，而且还可以得到零偏和非零偏的v s p 记录。( 3 ) 在复 杂模型中采用纵波震源进行了模拟，得到不同时刻的波场快照和单炮记录。在横 向各向同性介质中在应用纵波震源对均匀介质采用不同的系数进行模拟，比较了 这些不同系数产生的区别。 1 3 1 正演数值模拟 在这部分内容里主要研究和讨论了地震波场计算中需要应用的弹性波波动 方程，在不同类型的介质中弹性参数不同，则对应的波动方程有不同的表达形式， 这里重点讨论了位移波动方程二阶精度的中心差分格式和一阶速度一应力方程 的高阶交错网格差分格式，并且在前人的基础上对位移波动方程推导了均匀和非 均匀介质中的差分格式，对一阶速度一应力波动方程推导了其任意阶时间和空间 精度的交错网格差分格式，同时还研究讨论了交错网格差分方法的稳定性、频散 现象和边界条件，研究和讨论了不同震源正演模拟的区别。通过理论试算，验证 了使用吸收边界条件后的波场计算效果，并且针对一阶速度一应力方程进行了不 同震源不同介质的波场模拟计算。 1 3 2 模型试算及效果分析 对于位移波动方程设计了各种不同的地质模型，对于一阶速度应力方程则 分别设计了各向同性介质和各向异性介质模型，对一阶速度应力方程用高阶交 错网格有限差分方法进行了弹性波波场正演模拟计算。分析结果表明，一阶速度 第一章绪论 一应力方程的高阶交错网格法的效果要比位移波动方程的二阶精度的中心差分 格式的效果好，该方法的正演模拟结果比较清晰，能够较好地反映出各种弹性波 的传播特征，这就说明用高阶交错网格法的一阶速度一应力方程要比用中心差分 法的位移方程精度高，而且在做各向异性模型模型时也比较容易，只需要稍微改 变参数的形式，不用重新推导方程。通过正演模拟的结果可以帮助我们深入研究 和认识地震波传播规律和地下介质的分布情况。 6 长安大学硕士学位论文 第二章波动方程有限差分格式的建立 2 1 波动方程的建立 机械振动在弹性介质中的传播过程称为弹性波地震波则是机械振动在岩石 中的传播过程【2 3 1 。在一定的条件下，如温度、压力和作用力大小适当，岩石具有 弹性性质，可视为弹性体，从而地震波可视为岩石中传播的弹性波。建立波动方 程要利用弹性力学中的基本问题，即在已知物体本身的形状位置、弹性参数以及 外力分布情况下求物体的位移、应变和应力的分布状态。 由弹性理论可知，弹性体由应力表示的运动平衡微分方程式为1 2 4 j ： 等+ 等+ 警+ p f , 一p 粤o t = 。a x铆 a z 1 等+ 等+ 等+ 矾一夕窑o t = 。一十一r 一- r 删一u 一u瓠如 包 “ 。 等+ 等+ 丝o z + 阢一p 窑o t = 。叙却 2 2 ( 2 1 ) 式( 2 1 ) 中，t 为时间分量，u 、v 、w 分别表示弹性体中质点在( x ，y ，z ) - - 个方 向上的位移分量，他们是关于( x ，y ，z ) 及时间t 的函数，p 为弹性密度，为应 力，z 为外力分量( i ，j - - - x ，y ，z ) 。 在弹性力学中，表示弹性体在外力作用下体内各点的应变白与位移u 、v 、w 之间关系的几何方程为【2 4 】： o u g 属2 _ d 譬 加 2 万 o w 2 瓦 i = , , 1 , 0 + 2 1 z s 。= 肛砂 o x x = 彳, 8 + 2 m s 。仃秽= 肛叫 l o r 搿= 兄0 + 2 l z s 麒 仃砂2 f 叫 7 ( 2 3 ) 刁仁 锄一砂加一钯伽一知却一缸挑一砂抛一钯 i i = 1 1 秒 弦 口 鼬 印 勿 第二章波动方程有限差分格式的建立 其d p t r ( i ：x ，y ，z ) 为正应力，( f ，j = x ，y ，z ) 为切应力，岛( 汪x ，y ，z ) 为正 应变，毛( f ，y = x ，y ，z ) 为切应变，五为拉梅系数，秒为体变。 将( 2 3 ) 式带入( 2 2 ) 式，可得下面新的关系式： ( 2 4 ) 平衡方程( 2 1 ) 式是在物体处于平衡状态时的关系式，当有外力作用而使物体 处于不平衡时，物体就要作加速运动，此时要用牛顿第二运动定理来建立新的关 系式： a 2 “ 2 p 萨 a 2 ， 2 p 萨 a 2 w 。p 哥。 ( 2 5 ) 现在我们只考虑二维情况下x 和z 方向上的位移分量而忽略y 分量。则对于 ( 2 5 ) 式中的x 分量( 第一个方程) 和z 分量( 第三个方程) ，利用( 2 4 ) 式中的对应 关系，我们可以把它换算为： 当外作用力厂( 六+ z ) 为零时： a 2 ” p 哥 a 2 w p 哥。 ( 2 6 ) ( 2 6 ) 式就是二维非均匀介质中的弹性波方程。在均匀介质中，p 、见和不 随空间位置而变化，所以( 2 6 ) 式进一步简化为： 8 ) 、， ， ， 玑一硝饥一挑一锄锄一砂伽一砂锄昆 h 坝 坝 = = = 咖 咖 彩 织 织 破 织 碱 碱 + + + 抛一缸却一砂伽一瑟 “ 心 “ 矽 正 z + + + 眈i 盟七眈i 堕砂峨一砂魄一砂 眈i峨i魄卜毽 正 正 + + 刀川 锄一瑟 锄一玉 + + 加一锄 伽一苏 “ i1一il a 一瑟 a 一缸 + + 1j r_j 锄一锄伽一玉 l f + + 、，j、j跏一瑟伽一瑟 + + 锄一叙锄一魂 ，、，i 兄 允 a一缸a一瑟 川州训 望奶盟瑟 业锄垒厶nr引一 z 卜，r 吖斟旦缸 1 i r j 1 叫1 _ 1 z 0 、，、- 一挑一昆伽一七 + + 锄一苏锄一缸 ，l， 兄 兄 i_l i_。，l a一缸a一瑟 i l = 謦k 叼 长安大学硕士学位论文 f p 窘= c 名+ 2 ( 雾+ 塞 + ( 窘一是 1p 可a 2 w 地+ 2 p 芒蛊+ 善 + 捂一囊 仁7 9 第二章波动方程有限差分格式的建立 得到。这种方法是随着计算机的诞生和应用而发展起来的。其计算格式和程序的 设计都比较直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的 重要组成部分。 有限差分法是求解波动方程的一种近似数值解法，所求的近似解能否收敛于 真解取决于空间差分网格与时间延拓步长的大小，也即两者要满足一定的稳定性 问题，所以有限差分法的关键问题是网格剖分，尤其是在倾斜地层和复杂构造部 位。网格剖分直接影响着近似解的精度，一般来说，网格剖分的越细，精度越高， 但计算量也随之剧增；另外，所用近似波动方程的阶数越高，所求得的解精度也 越高，但用有限差分法解高阶偏微分方程存在着不少的问题和实际困难；同一偏 微分方程用不同的差分格式来近似也影响其解的稳定性，一般来说，在差分的过 程中所用到的网格节点数越多，其解的精度也相应要高，但同样也带来大的计算 量的问题。与其它数值模拟方法相比，有限差分法在理论上简单，实际应用也相 对简单，而且理论和实际的应用都比较的成熟。 2 2 2 有限差分方程的建立过程 有限差分法的具体操作可以分为两个部分：( 1 ) 用差分代替微分方程中的微 分，将连续变化的变量离散化，从而得到差分方程组的数学形式；( 2 ) 求解差分 方程组。 在第一步中，我们通过网格剖分法，将函数定义域分成大量相邻而不重合的 子区域。通常采用的是规则的剖分方式，最常用的是正方形网格。这样可以便于 计算机自动实现和减少计算的复杂性。网格线划分的交点称为节点。若与某个节 点p 相邻的节点都是定义在场域内的节点，则p 点称为正则节点；反之，若节点 p 有处在定义域外的相邻节点，则p 点称为非正则节点。 在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在 所有这些节点上的离散近似值。目前最常用的两种有限差分方法包括：基于位移 波动方程的二阶中心差分法和基于一阶速度一应力波动方程的高阶交错网格法， 前者算法简单，易于实现，但差分精度具有局限性，最后得到的是节点上鼽z 分 量的位移离散近似值，后者算法稍复杂，但可以提高差分精度，最终得到的是节 点上的位移速度离散近似值。本文将两种差分格式都进行了详细的推导，经过对 比可得一阶速度一应力波动方程的高阶交错网格法计算效率高，速度快，适于模 l o 长安大学硕士学位论文 拟复杂模型。 2 2 3 有限差分的差分格式 一个函数在x 点上的一阶和二阶微商，可以近似的用它所临近的两点上的函 数值的差分来表示。如对一个单变量函数厂( x ) ，x 为定义在区间【口，b 】的连续变 量。以步长办- - a x e , , ，b 】区间离散化，我们得到一系列节点五 - a ，恐= 五+ 五， x 3 = 屯+ 办= a + 2 a x ，毛+ l = 毛+ = 6 ，然后求出f ( x ) 在这些节点上的近似值。 显然步长j j l 越小，近似解的精度就越高。与节点而相邻的节点有毛一h ，而+ ， 因此在而点可以构造如下形式的差值： f ( x t + j j l ) 一厂( 墨) ，节点而的一阶向前差分 f ( x i ) 一厂( 而- h ) ， 节点薯的一阶向后差分 ( 而+ 五) 一厂( 而- h ) ，节点鼍的一阶中心差分 与毛相邻两点的泰勒展开式可以写为： f ( x ，- h ) = 厂( 玉) 一矿( 而) + 了h 2 _ ( 薯) 一酉h 3 厂( 再) + 可h 4 厂( 而) 一 ( 2 8 ) 厂( 而+ ) = 厂( 玉) + 矽( t ) + 等厂( 薯) + 鲁厂( 薯) + 等厂( 而) + ( 2 9 ) 利用( 2 8 ) 、( 2 9 ) 式，并忽略h 的平方和更高阶的项得到一阶微分的中心差分 表示： 厂( 咖丝掣掣 ( 2 1 0 ) z ，l 利用( 2 1 0 ) 和( 2 1 1 ) 式我们还可以得到一阶微分的向前，向后差分表示： 厂( 而) 丝攀攀盟 ( 2 1 1 ) 厂( 薯) 丝皇华生 ( 2 1 2 ) 将( 2 1 0 ) 和( 2 1 1 ) 式相加，忽略h 的立方及更高阶得到二阶微分的中心差分： 厂( x 3 盟业等巡 利用( 2 1 0 ) ( 2 1 3 ) 式，我们就可以构造出微分方程的差分格式。这里要指出 的是：在构造差分格式时，究竟应该选择向前，向后还是中心差分来代替微分方 程中的微分，应当根据由此得到的差分方程解的稳定性和收敛性来考虑。同时还 应该兼顾到差分格式的简单和求解的方便。 第二章波动方程有限差分格式的建立 2 3 二维弹性波动方程有限差分法 2 3 1 均匀介质中弹性波动方程的有限差分的建- 式r t 2 5 1 x 、z 为水平及垂直的直角坐标轴系，z 轴正向朝下，于是对二维情形可以用 两个耦合的二阶偏微分方程来描述p 波及s v 波，其运动方程为： 式中u 、w 分别为水平和垂直位移分量，口、为纵横波速度： 式中p 为介质的密度，叁为拉梅系数。 ( 2 1 4 ) 为建立方程( 2 1 4 ) 对应的差分方程，设x = m a x = m h ，z = n a z = n h ，t = l a t ，这 里h 是x 、z 方向的网格间距，出为步长。为方便起见，记： u ( m ，刀，) = 甜( x ，z ，) ，u ( m l ，刀l ，1 ) = “( x 办，z h ，t + a t ) 对于( 2 1 4 ) 式中的二阶偏导数，采用二阶中心差分格式来近似： 窘= 矿1 “( m n l ) 也( 聊朋+ “m , n , l - 1 ) ， 巫地2 业h 型 1 w ( m + l ，刀，) 一w ( m + l ，n - 1 ，) w ( m - 1 ，n + l ，1 ) - w ( m - 1 ，n - 1 ，) ( 2 1 5 ) = 去 w ( 所+ l 栉，旷w ( 坍+ 1 ，川，旷w ( 所_ 1 川，卅w ( m l ，) 对其它的二阶偏微分也采用类似的近似，对于方程( 2 1 4 ) 中的第一式可写成： 1 2 舄去堕劳謦 瞪 堆 叫 汁 为患如一掰伽一智 v 文 呐 略 ”一2 一0 挑一铲挑一酽 rh_j rl a 一瑟 。一舶 = i i 塑跳 长安大学硕士学位论文 【- “( 聊，刀，+ 1 ) 一三”( 聊，刀，) + 甜( 朋，刀，一1 ) = i a t = 等 “( 朋+ 1 ，叫) 也( 聊删+ 甜( m - l , n , 1 ) + 等肿u ) - 2 + u ( m , n - l , 1 ) + c 。_ 2 霸- - 乒p 一2 w ( 聊+ 1 ，疗+ 1 ，) 一w ( m - l , n + l , 1 ) ( 2 1 6 ) 若记f ：- a a - t ，：星，则上式可简写成： ”( 所，”，l + 1 ) = 2 u ( m ，胛，1 ) - u ( m ，刀，i - 1 ) + f 2u ( m 1 ，玎，1 ) - 2 u ( m ，刀，1 ) + u ( m - 1 明 + f 2 ，2 “( 聊，n + l ，1 ) - 2 u ( m ，刀，) + 甜( m ，n - 1 ，伽 ( 2 1 7 ) + ! 兰玉毛 垒【w ( 肌+ 1 ，”+ 1 ，) 一w ( 聊一1 ，刀+ 1 ，) - w ( m + l ，n - 1 ，1 ) + w ( m - 1 ，n - 1 ，咖 同理方程( 2 1 4 ) 中第二式可写成： w ( m ，刀，l + 1 ) = 2 w ( m ，”，1 ) - w ( m ，刀，1 - 1 ) + f 2 w ( 肌，n + l ，1 ) - 2 w ( m ，甩，咖w ( 朋，n - 1 ，伽 + ，2 ，2 w ( m + l ，刀，1 ) - 2 w ( m ，刀，1 ) + w ( m - 1 ，”，明 ( 2 1 8 ) + 掣【“( ， + 1 ，刀+ 1 ，) 一“( 朋一1 ，刀+ 1 ，) 一u ( m + l ，刀一1 ，) + “( 肌一1 ，刀一1 ，) 】 由( 2 1 7 ) 、( 2 1 8 ) 两个方程构成的显示格式，可利用前两个时刻，、，一1 的运 动来算出，+ 1 时刻的位移u 、w ，从而可以求出各个时刻各个网格节点上的值。 2 3 2 非均匀介质中弹性波动方程有限差分的建立【2 5 】 均匀介质弹性波动方程对介质为水平的情况是适用的，但当界面为倾斜或 更复杂时，使用上面的波动方程来处理时就存在很多困难，尤其是在边界条件的 处理上，要克服这些困难就得用非均匀介质情况下的弹性波动方程。本文为了更 好的模拟三层水平模型的波场记录，对于位移方程采用的是非均匀介质中弹性波 动方程的差分格式。 在二维非均匀介质的情况下，p 波及s v 波的运动方程为： 第二章波动方程有限差分格式的建立 a 耄2 u z 雩0 20 妻u + 口20 宝w 一2 p 20 宓w ) + 昙c 2 尝+ 2 量x 。2。9，0 窘= 昙( 口2 警+ 口2 罢一2 2 隶+ 言c j _ 1 丢o w o u+ 2 参 。 引入纵横波速度舐，上面两式变为： 当p 、五和为常数时，上面的式子就同均匀介质弹性波方程一致。 为了能差分出( 2 2 0 ) 式，前人做了大量研究工作，目前认为比较理想的方法 如下： 先考虑方程右端的一个典型项： 舡似z ) 瓦o u 若口2 ( x ，z ) 用节点( 所，厅) 处的离散值代替，并定义 口2 ( 朋割= 垫掣 则典型项的差分近似式为： 口2 ( 册+ 三，刀) “( 聊+ t ，”，) 一“( 聊，以，) 一口2 ( ，刀一互1 ，力) “( 研，刀，) 一甜m - l , n , 1 ) 缸2 再考虑一混合偏导项： 善 口2 ( x ，z ) 昙“( x ，z ，) 兰岳 c ( x ，z ，r ) 这里的c ( x ，z ，f ) 是为了简化公式而引入的。上面的式子就可以用一阶中心差分来 近似了： 奇c ( 圳) - c ( m , n + 1 1 , 1 ) - f c ( m 一, n - l , 1 ) ，c ( 朋删崭( m ，) u ( m + l , n , l z ) 似- u ( m - l , n , 1 ) 于是有： 1 4 02 亿 、，、， 锄一钯锄一瑟 2 2 口 廖 伽一缸跏一缸 口 口 ( ，l a一瑟a一缸 l i l ) + + 跏一勃础一苏 2 也 肛 格 乙 2 挑一昆锄一苏 矿 口 “ 斗 钆毽挑一如 2 2 ， 砬、 、= ：， ( a一苏a一勿 ”一2 w 一0 如一酽挑一扩 长安大学硕士学位论文 丢 口2 ( x ，z ) 丢 ( x ，z ，r ) =4 a x a z 钕2m ，力+ 1 ) u ( m + l ，n + l ，1 ) - u ( m - 1 ，n + l ，明一口2m ，疗一1 ) “( 肌+ 1 ，n - 1 ，1 ) - u ( m - 1 ，n - i ，明 用类似的方法就可得( 2 2 0 ) 式的差分方程： 去c 出2 竺：! 竺! ：! ! 竺：! 竺：尘 2 竺坠：! ! 竺：( 竺二! ：吐 2 u ( m + l ，刀，1 ) - u ( m ，刀，) u ( m , n , 1 ) - u ( m - l , n , 1 ) a x + 瓦1 口2m + l , n ) r v ( m + 1 , n + 1 , 1 1 ) - f w ( m + 1 , n 一- l , 1 ) 口2 m - l , n ) w ( m - l , n + 1 , l z ) 业- w ( m - l , n - l , 1 ) 一石1 2 ( m + 1 ，刀) 2 ( m - i ，刀) w ( m + l ，n + l ，) 一w ( m + l ，r 一1 ，) 2 z w ( m - l , n + 1 , 1 z ) 心- w ( m - l , n - l , 1 ) z 心 + 去【2 m , n + 1 ) w ( m + 1 , n + _ 1 , 1 ) 石- w ( m 一- l , n + 1 , 1 ) 2m ，n - 1 ) + 击r w ( m + 1 , n - 1 , l z ) 出- w ( m - 1 , n - 1 , 1 ) z 出 壁：( ! ：! ! ! 壁：( 竺：竺2 2 壁：( 竺：竺! 壁：( 竺：竺二1 2 以及： 2 u ( m ，n + l ，1 ) - u ( m ，2 ，) u ( m , n , 1 ) - u ( m , n - l , 1 ) a z ( 2 2 1 ) 第二章波动方程有限差分格式的建立 ! ! 竺：! ：! ! ! 二三! ! 竺：! ：! ! ! ! ! ：竺：! 二! ! ： 址2 1 a 2 ( m , n + l ! + a 2 ( m , n ) w ( m , n