（地球探测与信息技术专业论文）基于波动方程的瑞利面波有限差分法正演模拟.pdf

摘要 利用瑞利波信息反演层状半空间介质的性质和状态，是地震勘探、 岩土工程领域关注的研究课题。本文根据弹性波方程，采用交错网格有 限差分数值求解算法，对浅层各向同性弹性介质进行了包括瑞雷面波和 体波在内的全波场模拟。在数值模拟中，采用地震勘探中常用的爆炸点 源激发产生的瑞利面波，讨论了震源埋深对瑞利面波激发的影响。 为了更好的处理模拟过程中的边界条件，在正演模拟使用吸收边界 条件，并给出了角点的处理方法。针对工程勘察中常见的连续和层状介 质模型进行了模拟，获得了更加接近实际情况的地震记录。结合模拟记 录，探讨了瑞雷面波的形成条件。 此外，本文最后分析了现存的瑞利面波资料处理及反演方法。对在 资料处理过程中各环节不同方法的优缺点做了详细的讨论。对反演方法 的未来发展方向作了展望，对当前极有潜力的时频分析方法进行了详细 的了解，对其在反演中的应用前景及其本身的局限性做了论述。事实表 明时频分析方法获得的频散曲线不及多道频率波数分析得到的瑞利波频 散曲线准确。不过由于时频分析仅需要一个接收道的数据就能给出结果， 比多道频率波数分析方法要简便，如果能对此方法加以改进，还是一种 具有前景的分析方法。 总的来说，本文系统地给出了用有限差分进行瑞利面波在分层半空 间介质的正演模拟方法，并写出了相应的程序。为用瑞利面波探测地层 结构奠定了可靠的基础。在此基础上，可进一步发展相应的反演解释方 法，可望能提供给勘探地层层状结构以指导。 关键词：瑞利面波，全波场模拟，有限差分，交错网格，边界条件，反演 a b s t r a c t ni sat o p i ci nt h ef i e l d so ft h es e i s m i ce x p l o r a t i o na n dg e o - t e c h n i c a l e n g i n e e rt oi n v e r s et h ep r o p e r t i e so f t h es t r a t i f i e dm e d i ai nh a l f s p a c eb yu s i n g t h er a y l e i g hw a v e s t h i sp a p e ra p p l i e st h ef i n i t ed i f f e r e n c em e t h o dw i t h s t a g g e rg r i d st os i m u l a t ef u l l w a v e f i e l d ，i n c l u d e sb o t hr a y l e i g hs u r f a c ew a v e a n db o d yw a v e s ，f o rs h a l l o wh o m o g e n e o u sm e d i ab a s e do ne l a s t i cw a v e e q u a t i o n s i no u rn u m e r i c a ls i m u l a t i o n ，t h er a y l e i g hw a v e sw e r ee x c i t e db y t h ee x p l o s i o no fap o i n ts o u r c e n i sm e t h o di sw i d e l yu s e di nt h es e i s m i c e x p l o r a t i o n t h ei n f l u e n c eo ft h ed e p t ho fe p i c e n t e rt ot h ee n e r g yo fs u r f a c e w a v ei sd i s c u s s e da l s o m e a n w h i l e ，aa b s o r p t i o nb o u n d a r yc o n d i t i o nw h i c hm a k e st h ed i s p o s a l o fb o u n d a r yc o n d i t i o nm o r es i m p l ea n dm o r ee f f e c t i v ei sp u tf o r w a r d f u r t h e r m o r e ，am e t h o dt op r o c e s sa n g l ep o i n t si sa l s op r o p o s e d s ot h e b o u n d a r yc o n d i t i o ni sp r o c e s s e dv e r yw e l l s t i m u l a t i o n sh a v eb e e nd o n ef o r c o n t i n u o u sa n ds t r a t i f i e dm e d i aw h i c ha r ec o m m o n l ys a wi nt h ef i e l do f e n g i n e e r i n gs u r v e y , t h e r e a f t e r , s e i s m o g r a m sw h i c ha r ev i v i dt or e a l i s t i co n e s a r eo b t a i n e d a c c o r d i n gt os t i m u l a t e dr e c o r d s ，t h ep r o d u c e dc o n d i t i o no f r a y l e i g hs u r f a c ew a v ei sp r o b e d i na d d i t i o n ，m a n yp r o c e s sm e t h o da n di n v e r s i o nm e t h o db es t u d i e di nt h i s p a p e r t h ea d v a n t a g ea n dd i s a d v a n t a g eo f t h ep r o c e s sm e t h o dt h a tb eu s e di n d i f f e r e n tp r o c e d u r ea n ds t e pb ed i s c u s s e dv e r yd e t a i l e d t h en - e n do fi n v e r s i o n b ep r o s p e c t e d t h et i m e f r e q u e n c ya n a l y s i s ( t f a lm e t h o dw a sc o m p r e h e n d e d a n dd i s c u s s e di t d e v e l o pa n dd i s a d v a n t a g ei n f u t u r e t h er e s u l to ft h e d i s p e r s i o nc u r v e so b t a i n e df r o mt f aa r en o ta sg o o da st h o s ef r o mt h ef - k a n a l y s i sm e t h o d b t i tt h ea d v a n t a g eo ft h et f ai st h a ti to n l yn e e d st h ed a t a f r o mm u l t i c h a n n e l sa r en e e d e d i f s o m ei m p r o v e m e n ta d d e d t h et f ac o u l db e e x p e c t e dt ob eap r o m i s i n gm e t h o d i naw o r d ，t h i sp a p e rs h o wt h em e t h o dt h a ta p p l i e st h ef i n i t ed i f f e r e n c e m e t h o dt os i m u l a t ep r o m u l g a t i n go f r a y l e i g hs u r f a c ew a v ei ns t r a t i f i e dm e d i a a n dt h e c o r r e s p o n d i n gi m p l e m e n t a t i o ns o f t w a r e s i m u l a t i o n r e c o r dw i l l b e c o m et h eb a s eo ft h es t r a t i f i e dm e d i ad e t e c t i o n m o r e ，t h ei n v e r s i o nm e t h o d c a nb ed e v e l o p e da n df i n dt h ed i r e c t i o no f p r o s p e c tf o rs t r a t u m k e y w o r d s ：r a y i e i g hs u r f a c ew a v e f u i 卜w a v e f i e i dm o d e ii n g ，f i n i t e d i f f e r e n c e s t a g g e rg r i d b o u n d a r yc o n d i t i o n ，i n v e r s i o n l 硕卜学位论文 1 1 瑞利面波的发展概况 第一章绪论 瑞利面波足沿地表传播的一种弹性波。早在1 8 8 7 年，英国数学物理学家r a y l e i g h 在求解自由表面半空间中的平面弹性波场时，就预占了一种沿纵向指数衰减的面波的 存在，此面波的速度小于横波的速度，质点运动轨迹为一逆向椭圆。后来人们从天然 地震记录中证实了这种匾波的存在，并称其为r a y l e i g h 波，也叫瑞利面波。 在天然地震中，瑞利面波的能量最强、破坏性最大，早期人们曾将瑞利面波视为 干扰波。5 0 年代仞，人们发现了瑞利面波在层状介质中具有频散特性，并根据这一 特性利用天然地震记录中的瑞利面波信号来研究地球内部结构“1 。6 0 年代初， h o y k a l l e n 等提出面波勘探的半波长解释方法，利用微震器产生正弦波信号，实测道 路断面的速度分布，将稳态瑞利面波首先应用于地基勘察。8 0 年代初，日本v i c 株 式会社推出了g r 8 l o 佐滕式全自动地下勘探系统，并应用于工程勘察中。3 。1 9 8 9 年，杨成林”1 也自行研制了稳态瑞利面波勘探系统，并将其应用于第四系地层分层和 地基处理效果评价。 1 9 8 2 年。s t o k e ”1 等采用锤击震源，通过两个检波器之间的互谱相位信息求取面 波的相速度，成为最初的瞬态面波勘探试验。1 9 8 6 年，n a z r i a n “1 等用表面波谱分析 方法( s a s w ) 对高速公路路面及路基进行了探测，为瞬态瑞利面波法在工程中的广泛 应用奠定了基础。1 9 9 6 年，刘云祯3 等研制了s w s 瞬态面波多道数据采集处理系统， 并将之应用于机场工程勘探、浅层煤田勘探、地下煤巷探测等方面的工作中，取得了 好的效果。 自从r a y l e i g h 发现瑞利面波之后，国内外大苣的学者对瑞利面波在各种介质中 的传播理论进行了广泛而深入的研究。瑞利面波勘探技术近几年来在工程领域中的应 用得到了很大发展，应用的范围也逐渐扩大，如在地质分层、地层物理力学参数计算 和密实情况判断、寻找地下空洞、探测地下目标物体、检测灰土及强夯地基密实情况、 公路路基及强化路面检测、滑坡调查等方面，部取得了较好的应用效果。目前，瑞利 面波勘探己成为浅层或超浅层地球物理勘探和工程岩体及施工质量检测的主要勘探 手段之一。 碗十学位论文 1 2 瑞利面波的正演研究现状 1 2 1 标准层状介质中瑞利面波传播的正演研究 界面互相平行，且界面处连续的均匀各项同性自由表面层状半空间弹性介质称之 为标准层状介质。对于在这种介质中传播的瑞利面波的频散曲线的正演问题，大量学 者进行了研究，并提出了各种方法，取得了较好的效果，也为其他复杂的层状介质的 正演问题奠定了基础。其主要的方法有t h r o m s o n - - h a s k e l l 法、s c h w a b - - - k n o p o f f 法、 8 矩阵法、a b o - z e n a 法、r t 矩阵法等。 1 9 5 3 年，h a s k e l l “1 在t h r o m s o n 的研究基础上，通过相邻两界面的传递矩阵公式 以及自由表面边界条件和无穷远处的辐射条件导出了层状介质中平面瑞利面波的频 散方程，人们称之为t h r o m s o n - - h a s k e l l 法。 为了避免高频数值溢出及精度丢失，k n o p o f f ”1 在h a s k e l l 的矩阵方法的基础上， 提出了一种新的求取频散曲线的方法，s c h w a b 州通过一系列的行列式变换将矩阵变为 块对角阵，减小了计算量，这种方法称为快速s c h w a b - - - k n o p o f f 法。 1 9 6 3 年，p e s t e l 和l e c k i e “”介绍了6 矩阵法。1 9 6 5 年，t h r o w e r 1 ”和d u n k i n “” 首次将6 矩阵法用于瑞利面波频散方程的推导，有效避免了数值精度的丢失。1 9 9 6 年，b u c h e n 和b e n h a d o n 提出了快速6 矩阵法。1 9 9 8 年，i v a n s s o n “”用更简便的方 法得到了快速6 矩阵法。 1 9 7 9 年，a b o - z e n a “”用一系列4 4 阶反对称矩阵的循环计算得到了瑞利面波的 频散方程，解决了高频数值不稳定问题。m e n k e “5 1 解释了a b o z e n a 解释了反对称矩阵 的物理含义，发现此矩阵与p - s v 波的基础解密切相关，并引入了d 矩阵和e 矢量， 导出了新的瑞利面波频散方程，此方法既有效的避免了高频数值不稳定性问题，又加 快了计算速度，并适用于有低速夹层存在的层状介质。1 9 6 8 年，b e n - m e n a h e n “”等提 出b 、p 、c 坐标系，大大简化了层状介质中三维声波问题的求解过程。1 9 8 2 年，李 幼铭等“”在b 、p 、c 坐标下将a b o - z e n a 的算法采用h a s k e l l 矩阵本征值分解的6 矩 阵运算用五个形式简单的矩阵连乘形式重新组织了a b o - z e n a 算法，使运算步骤大量 简化。1 9 9 6 年，张碧星“”在a b o - z e n a 和m e n k e 研究的基础上用b 、p 、c 坐标系进一 步研究了层状介质中的频散问题，将m e n k e 的e 矢量的传递矩阵f 分解为简单的三个 矩阵的乘积形式，其中两个与频率无关，另一个为对角阵，从而简化了计算，且避免 了数值不稳定性i 口j 题，并提出了关于e 矢量的一个新结论，定义了两个物理量61 、 b2 ，发现6l 、b2 与介质的层状结构密切相关，还可以用来辨别接收到的信号是否 为瑞雷面波。 1 9 7 4 年及1 9 7 9 年，k e n n e t t “”圳，提出并发展了r t 矩阵法，此方法使用了反射 2 硕十学位论文 和透射子矩阵技术，巧妙地避开了数值不稳定l 口j 题。但这种方法要求所有的计算均使 用复数运算形式，故计算量较大。 以上介绍了几种求职层状介质中瑞利面波频散方程的经典方法。 1 2 2 非标准层状介质中瑞利面波传播的正演研究 对 e 杯准层状介质，如果面不连续、界面不平坦，倾斜界面、非弹性层状介质等 情况，瑞利面波的频散曲线的求取以及瑞利面波波场的正演研究叫层状介质更为复 杂，有些研究以保证层状介质的正演理论为基础，有些则利用到了其他方法。 1 9 7 1 年，s c h w a b 等”基于b i s w a s k n o p o f f 变换研究了非弹性层状介质中面波的 频散衰减舰律，并应用于径向椎均匀| 乍弹性介质中的乐夫面波波场的频散及衰减计 算。1 9 7 8 年，t i e r s t e n 和s n i e d e r 啪1 利用线型压电弹性体方程的格林函数公式得到 了压电弹性体的扰动理论，由特征值的第一扰动方程的反复使用得到了想要的精度的 面波衰减，以及大波长在覆盖层介质表面传播的相速度频敬。1 9 8 6 年，s n i e d e r 1 用线型化散射理论对不规则表面的层状介质中的面波散射进行了研究，并计算了周期 为1 5 秒到3 0 秒之间的面波基础模相速度的扰动。1 9 8 6 年，b h a t t a c h a r y a “1 在6 矩 阵法及a b o z e n a 的方法的基础上，提出了球形层状介质中瑞利面波的频散方程的求 解方法，此方法对高频情况不会出现数值精度丢失问题。1 9 8 8 年，m a l 对层状介质中 含非完好连接界面时的面波场进行了研究，用矩阵方法分别对两种j f 完好连接界面模 型进行了面波频散方程的求解，此方法也避免了高频数值不稳定问题。1 9 9 1 年， r i c h a r d s 等瞄1 提出了弹性波在有倾斜界面存在的层状介质中传播的射线理论，对面 波有倾斜界面存在的层状介质中传播规律有一定的指导作用。1 9 9 6 年，牛滨华、何 继善渊从虚功变分原理出发，引入等效面波震源概念，用有限元方法导出了含多个 倾斜界面的层状介质以及垂向各向异性和粘滞性多层介质中瑞利面波的频散方程，并 进行了数值计算。1 9 9 7 年a r o r a 等。1 利用t h o m s o n - - h a s k e l l 矩阵也导出了球形层状 介质中瑞利面波的频散方程，将各层的传递函数表达为一种简单明确的形式。1 9 9 9 年，b e n h a d o r 和8 u c h e n ”1 基于w h i t r a m 的一阶扰动理论导出了存在倾斜界面的层状 介质中乐夫面波和瑞利面波频散方程，并数值计算了二层及三层模型下面波的频率、 波数、幅度及相位的扰动。 非标准层状介质中的含倾斜界面的层状介质中瑞利面波研究对于实际的地球物 理勘探具有较大的现实意义，在工程中它们比层状介质更贴近于实际地层情况，在层 状介质的瑞利面波法日益完善之后，人们无疑应大力度开展这种介质中的瑞利面波法 勘探研究。 相对于非标准层状介质，瑞利面波在存在裂缝、缺陷、拐角、阶梯、障碍物、空 洞、沉积盆地等复杂结构的介质中的传播规律更为复杂，几乎再也不能以层状介质中 硕十学位论文 的瑞利面波理论为基础进行研究了，目前研究这类问题的方法形形色色，而求解的问 题己远离简单的频散和衰减，大多数研究均以求解面波散射场为主要目标，求解的方 法大多用到了数值近似方法，如边界元法、有限元法、有限差分法、b o r n 近似法、 波函数展开法、扰动法等，这些研究对充实瑞利面波理论及弹性理论有不可估量的作 用。 1 3 本文的主要工作 利用瑞利波信息反演层状半空日j 介质的性质和状态，是地震勘探、岩土工程领域 关注的研究课题。本文由层状半空间瑞利面波的实验数据，模拟了瑞利面波在层状半 空间中瑞利波的传播。本文根据弹性波方程，采用交错网格有限差分数值求解算法， 对浅层各向同性弹性介质进行了包括瑞雷面波和体波在内的全波场模拟。在数值模拟 中，采用地震勘探中常用的爆炸点源激发产生的瑞利面波，讨论了震源埋深对瑞利面 波激发的影响。针对工程勘察中常见的连续和层状介质模型进行了模拟，获得了更加 接近实际情况的地震记录。结合模拟记录，探讨了瑞雷面波的形成条件。 本文各章节内容如下： 第一章简述了瑞利面波勘探的发展概况及正演模拟的现状。 第二章主要阐述了瑞利面波在均匀半空日j 介质中以及层状半空日j 介质中传播的 理论基础。 第三章对波动方程的有限差分解法作了详尽的描述。 第四章为对各种模型的正演模拟。 第五章对瑞利面波的资料处理及反演作了详细的讨论，力争找到改进现有方法缺 点的突破方向。 第六章足全文的总结及展望。 4 硕 学付论文 第二章瑞利面波的基本理论 2 1 弹性介质中的波动方程 介质因受力而产生形变，当外力撤销时，能完全恢复原来状态的介质为理想弹性 介质。介质的弹性性质不仅同空间方向无关，而且同空间分布无关的理想弹性介质为 均匀各向同性介质。 2 1 1 应力 弹性介质受夕l - ) j 后发生弹性形变，由于质点间具有弹性联系，介质内部产生反抗 形变的内力与外力抗衡，称为内应力。作用在单位面积上的内应力称为应力。 应力在直角坐杯系中是二阶张星，共有九个分量，称为应力张垦。对于无旋的均 衡应力体，在这九个分量中应有f 。= 7 因而只有六个独立分量，故应力为二阶对 称张量。应力通常排成一个矩阵 h 勺k 1 l 2 1 ) l 吃勺乞j 2 1 2 应变 所谓应变是指当物体受到应力的作用时，它的大小和形状的相对改变量。应变在 隧刭 汜z ， 2 1 3 应力与应变的关系 胡克定律表示了应力和应变之间的线性关系。而表达应力和应变关系的应力与应 变方程也称为广义胡克方程，即本构方程。利用本构方程可完全地确定各向异性介 质的应力和应变系统，其数学表达式如下 硕p 学位论丈 c l ， c 2 ， g 5 c 4 5 g 5 c 6 5 各向同性介质中，弹性常数仅有f 再个，把它表示成拉梅常数，即 c i 2 = c 1 3 = c 2 i = = c 3 2 = g l = 旯 c 4 4 = g 5 = c k = 2 z c 1 i = c = c 3 3 = 3 , + 2 k t 其余的2 4 个系数都为0 ，应力应变关系可简化为 k = ( 五+ 2 ) 气+ a s w + a 屯 f 盯= 珧h + a + 2 8 眄+ 勉。 t 。= a 。+ 娩w + o r + 2 z ) g = f q = 2 坪口 f 口= 2 坪 f 4m m 2 ，昭“ 2 1 4 运动方程 ( 2 3 ) ( 2 4 ) ( 2 5 ) 假设质点是在没有体力作用下作自由振动，由牛顿定理可得到各向同性介质中的 矢量弹性波波动方程： 户箬= 以 4 - 2 【f ) v v 一胛v ” ( 2 6 ) 叫 其中“为位移矢量，可以用以下的势函数来表示 甜= v 伊+ v x ( z e ：) + v x v x ( r e ：) ( 2 7 ) 代入波动方程( 2 6 ) 式可得到 其中的是妒， - f ，z 分别对应着时间域中p 波，s v 波和s h 波的三个势函数。 巧= i z 忑刃万，k = 刀万分别为纵波和横波的波速，且 告= 厩- j 掐2 0 - a1 五+ 劲1 ) 6 ( 2 9 ) 和和彩缸和跏函易锄艮厶厶 o o o 已。靠 函白氏如如靠 厶凸凸d以& 缸啊勿缸伽撕 互 矿 z 甲| n 脚一扩嘞一扩押一酽 上哆上曙上曙 = = = 矿 z r 吻 吮 咖 硕十学位论文 这黾6 为泊松比，盯= a 2 以+ ) ，它表示物体的纵向拉伸变化与横向收缩变化之比，其值在 o _ o 5 之间变化，对于1 常坚硬、刚性很强的岩石，泊松比可达0 0 5 左右；随岩石刚度的降低， 泊松比增大。对f 软的、胶结度羊的岩石，泊松比可达0 4 5 。 由( 2 8 ) 式可知，矢量波动方程( 2 6 ) 式被转换成了形如( 2 8 ) 式所示的三个标量 波动方程。假设场量与y 无关，仅是x 和z 的函数，则( 2 8 ) 式具有如下形式的通解 i 妒= q ，o ( z ) e “训 甲= ( z ) p “训 ( 2 1 0 ) lz = x o ( z ) e “1 将( 2 1 0 ) 式代入( 2 8 ) 式，可得 其中a 2 = 睇一k 2 ， 睁口2 卜= 。 睁6 2 卜= 。 眩 ( 导彬 腓瑚 6 2 _ 搿叫= 号，霹= 等该式的通解为 j 仍( 力= 彳p 2 + b e 。2 ( z ) = “+ 眈“ ( 2 1 2 ) i x o ( z ) = e e “+ f e 。妇 2 2 均匀半空间中的瑞利面波 2 2 1 均匀半空间中瑞利面波的形成 假设再自由界面以下的无限办空间充满了弹性常数为和以及密度为的均匀弹性 介质，x 、l ，轴取在自由表面上，z 轴指向弹性介质内部，z 0 为真空区域，在z = 0 平面上，波场满足自由边界条件，即应力为零。 由于z 轴为体系对称轴，而且，在这罩只考虑s v 波，对于平行于自由表面运动 的剪切波( s h 波) 将不作考虑，故引入位势移 甜= v r p + v x v x ( e ，) ( 2 1 3 ) 将( 2 1 3 ) 式代入( 2 8 ) 式中，可得 巾：胪证p 铡+ v x v x 防一p 务 - o 汜 从而 7 硕十学位论文 ：三壹 砖西2 1a 2 甲 2 虿万 ( 2 1 5 ) 设瑞利面波的速度为珞，z o x 平面内沿z 轴方向传播，在j ，轴方向的振幅和 相位完全相同，即讨论平面二维情况。令其位移位为 艇八咖= (216)g(z)e【t f ，= 俩“ 静争功删 汜 l 窘卅瑚删 “ 其中2 昙，b 2 詈，k2 罢。式( 2 1 7 ) 的解为 1 缈：一+ & 。 ( 2 1 8 ) f 靠= 挥虿 【6 - 舸 上面系数c 、d 的项当z 趋于无穷大时， 应为零。于是有 ( 2 1 9 ) 不满足远处收敛的条件，因此，c 和d i 伊= a e 一。e “k r x - 。t i = b e “一( k r z “) 在自由界面的边界条件是正应力和切应力为零， 即 弹性常数五和与介质密度及纵波、横波波速的关系分别为 8 ( 2 2 0 ) ( 2 2 1 ) ( 2 2 2 ) 硕 学位论艾 j 五= a v ；一2 )( 2 2 3 ) 垆讥2 噍塞犁。 汜。， ， i z 塞+ 害一害，l = 。 。 j 西( 口2 一) 妒+ 2 蝌2 唧2 一i k k b 妒 ) = o( 2 2 5 ) p 一2 = 兰2 2 2 i k r 撅一霹拈o ( 2 2 6 ) 将式( 6 2 0 ) 式代入式( 6 2 2 ) 式，并利用式( 6 2 6 ) 消去b 可得 取其实部 式( 6 2 9 ) 为瑞利面波的位移表达式，可见，当z 一一时，u 。一o ，u 。一o ，即在x 和z 方向的位移都为零，说明瑞利面波的分布深度是有限的。当介质为泊松固体( i l = 0 2 5 ) 时 9 一中霉 一 p | 旦k p 一 r 卸l r 刎l 醐 吣 k 口 旦k垫砜 硕十学位论丈 “ ( 2 3 0 ) l b = 0 3 9 3 3 k 。 将式( 2 3 0 ) 和式( 2 2 7 ) 代入式( 2 2 9 ) ，得 j 虬2 4 k e 。0 ”一o 5 7 7 3 p ”。) s i n ( 耐，x ) ( 2 3 1 ) 【“z = a k r ( - o 8 4 7 5 e _ 0 8 4 7 5 h7 + 1 4 6 7 9 e - 0 3 9 3 3 k z ) c o s ( 吐蟮一k 工) 当z = o 时，即在自由表面上， m 扩0 4 2 d 8 瓤卅) ( 2 3 2 ) l 也l 0 6 2 d c o s ( 耐一k 月x ) 其中，d = a k 。将式( 2 3 2 ) 的两个等式平方后相加并整理，得 志j _ + 志j = s i n 2 ( c o t k r s 2 ( t o t k r ) - 1 ( 2 3 3 ) 式( 3 3 3 ) 为椭圆方程，这表明在自由表面附近沿波传播方向得垂直平面内，瑞 利面波质点运动的轨迹时椭圆，椭圆的水平轴与垂直轴之比约为2 ：3 ，且质点的垂直 位移比水平位移相位超前 。 介质的泊松比为0 2 5 时，根据式( 6 3 1 ) 可以计算水平位移u ，和垂直位移u 。 的振幅随深度及泊松比o r 的变化，如图( 2 1 ) 所示。从图中可以看出，当z 2 o 1 9 3 ，两者的符号不同，质点振动轨迹为顺时针转动的椭圆。质 点振动轨迹和振幅随单位波长深度的变化如图( 2 2 ) 所示。 图2 1 均匀半空间中瑞雷波位移分量与深度x 的关系 图2 2 瑞雷波中质点的运动轨迹示意图 1 0 硕十学位论丈 2 3 层状半空间中的瑞利面波 在层状介质条件下，可以寻求一个面波传播的解析解。如大地为多层弹性介质， 如均匀半空日j 条件一样。瑞利面波仍在z o x 平面内传播，在这样的条件下，在自由表 面上，仍有2 个边界条件，垂直应力和水平应力为零；有两种介质的分界面处的边界 条件为：垂直位移和水平位移连续，因此有4 个边界条件。对于n 层介质，计算面波 的传播问题的解共有4 n - 2 个边界条件，即有4 n - 2 个齐次联立方程。为了了解层状介 质对面波速度的影响，取一种简单的2 层大地情况，对位函数的取法见式( 2 2 0 ) ， 假如介质的自由表面之上有一个非弹性覆盖层，设该层的厚度可以忽略，坐标的取法 如图( 2 3 ) 所示。在这样的条件下，由于覆盖层的影响，该分界面上法向应力不再 等于零，而是等于p 兽，切向应力仍然等于零( 因为覆盖层是椎弹性介质) 。这时的 边界条件为 或者写成 卜恪豺2 t ( 0 甜2 0 + 斟p 鲁t 警割叫。， 眵z 塞一铡。= 。 将式( 2 2 0 ) 代入式( 2 3 5 ) ，并化简得到 l k ( 2 一) 一2 碟一七；j 4 + f p 肚一t ；一砖一p r 0 2 p = o ，。m i 一2 i k 。一一， + 毕h 2 一k ；一霹归- - 0 令这一方程组的系数行列式等于零，就得到 ( 2 碟一牡( 2 睇一) 一p r 0 2 、f z ；= 虿j _ 2 k 虿= 巧( 2 肚。0 虿二虿一2 ) = o ( 2 3 7 ) 由上式可以看出，此时面波速度的解与频率有关，即存在画波的速度频散。所谓 速度频散是指谐波传播的速度随频率( 或波长) 变化。在瑞利波位移位的表达式中， 有因子p “k x - “t ) = e b r ( 嘞”，这假定瑞利面波足渚波，而v 。是谐波同相面的传播速度， 称之为相速度，用v 表示。实际上地震波是有许多不同频率、不同振幅的谐波叠加而 成，各谐波按各自的相速度传播。地震波的频谱一般是连续的，一些频率的谐波波峰 遇到一起会相互叠加使振幅增大，反之会互相抵消使振幅减小，这种叠加而成的大振 幅传播的速度称为群速度，用u 表示。因此波的群速度即是波的能量传播的速度。在 物理学中，相速度和群速度的关系为 一 睾卜 等蒜 硕t 学位论文 硕十学位论文 第三章波动方程的有限差分解法 3 1 有限差分法简介 有限差分法是数值求解常微分方程或偏微分方程的一种方法。地球物理学的许多 问题在被分析研究之后，往往可以归结为偏微分方程的求解问题。一般说来，处理一 个特定的地球物理学问题，除了需要知道它满足的数学方程外，还应当同时知道这个 问题的定解条件，然后才能设计出行之有效的计算方法来求解。 有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续 取值。在有限差分方法中，放弃了微分方程中独立变量可以取连续值的特征，而关注 独立变量离散取值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满 意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间隔，或者 通过离散点上的函数值插值计算来近似得到。这种方法是随着计算机的诞生和应用而 发展起来的。其计算格式和程序的设计都比较直观和简单。因而，它的实际应用已经 构成了计算数学的重要组成部分。 3 1 1 有限差分法的差分格式 一个函数在x 点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值 的差分来表示。如对一个单变量函数f ( x ) ，x 为定义在区日j a ，b 的连续变量。以步 长h = x 将 a ，b 区间离散化，我们得到一系列节点x l _ a ，x z = x ，+ h ，x 3 - - - - - x ：+ h = a + 2 a x ，】【n 1 - x 。+ h = b ，然后求出f ( x ) 在这些点上的近似值。显然步长h 越小，近似 解的精度就越i 茼。与节点x 。相邻的节点有x l + h 和x t - h ，因此在点x 可以构造如下形 式的差值： 厂( + ) 一( t ) 节点x 。的一阶向前差分 ( 葺) 一八一h ) 节点x 。的一阶向后差分 f ( x j + - f ( x , - h ) 节点x i 的一阶中心差分 与点相邻两点的泰勒展开式可以写为 厂( x t - - ) ：，( ) 一h f 。“) + 等厂( ) 一譬厂瓴) + 等厂r ) ( 3 1 ) j二 ，( 为+ 厅) ：厂( t ) + 妒( t ) + 等，( t ) + h _ 3 厂( t ) + 等，( t ) + ( 3 2 ) 忽略( 3 1 ) 式和( 3 2 ) 式中h 的平方和及更高阶的项，得到一阶微分的向前、向后和 硕学位论文 中心差商表示： ，( ) ；丝掣 ( 3 3 ) ，( 薯) 。f ( x t ) - 瓦f 一( x , - h ) ( 3 4 ) 厂( _ ) 。f ( x , + _ h ) i - f 一( x , - h ) ( 3 5 ) 将( 3 1 ) 和( 3 2 ) 式相加，忽略h 的立方及更高阶的项得到二阶微分的中心差商表示 ，。( 一) 。丝坐哩等也盟生 ( 3 6 ) 月 利用( 3 3 ) 一( 3 6 ) 式，就可以构造出微分方程的差分格式。在构造差分格式时， 究竟应该选择向i i ，向后还是中间差分或差商来代替微分方程中的微分或微商，应当 根据由此得到的差分方程解的稳定性和收敛性来考虑，同时兼顾到差分格式的简单和 求解的方便。 3 1 2 有限差分法的计算步骤 有限差分法的具体操作分为两个部分： ( 1 ) 用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方 程组的数学形式； ( 2 ) 求解差分方程组。 在第一步中，通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子 区域。通常采用的是规则的分割方式。这样盯以便于计算机自动实现和减少计算的复 杂性。网络线划分的交点称为节点。在第二步中，数值求解的关键就是要应用适当的 计算方法，求得特定问题在所有这些节点上的离散近似值。 3 1 3 有限差分法的边值条件 利用微分的差分表示，可很容易地将微分方程离散化为差分方程的形式。但足由 差分方程所得的解完全取决于待求微分方程的特性。边界条件的情况变化将会引起差 分方程组的不同。在求解微分方程中会遇到两类问题：一类足初值问题；另一类是边 值f j 题。在地球物理学的初值问题中，函数值和函数偏导值是给定的。而在边值问题 中，由于模型的有限，必须作大量的工作消除边界引起的各种问题。本文中重点考虑 边值问题。 假定某方程形式上可以写为： 工= g ( 3 7 ) 其中l 为含偏微商的算符。 1 4 硕 学位论文 它的边界条件一般可写为： t 4 b a + g l 等i c _ 9 2 g 表示场域d 的边界，g t ( s ) ，9 2 ( s ) 为边界上s 点的逐点函数。 三类边界条件： ( 1 ) 第一类边界条件，或称为狄利克莱( d i r i c h l e t ) 问题( g l = 0 ，9 2 o ) 。 硎g = g ( j ) ( 2 ) 第二类边界条件，或称诺伊曼( n e u m a n n ) 口j 题( g l o ，9 2 - o ) 。 娑i 。：g ( 5 ) 一i g 一5 、4 ， ( 3 ) 第三类边界条件，或称混合问题( g l o ，9 2 0 ) 。 酬川心) 鼍i o = 9 2 3 1 4 区域d 的离散化 ( 3 8 ) ( 3 9 ) ( 3 1 0 ) ( 3 1 1 ) 即通过任意的网格划分方法把区域离散为许许多多的小单元。原则上讲这种网格 分割是可以任意的，但是在实际应用中，常常足根据边界g 的形状，采用最简单，最 有规律，和边界的拟合程度最佳的方法来分割。常用的有正方形分割法和矩形分割法 ( 如图3 1 ) 。有时也用三角形分割法( 见图3 2 ) 。对圆形区域，应用图( 3 3 ) 所示 的极网络格式也许更方便些。这些网格单元通常称为元素，网格点称为节点。 3 1 求解区域的矩形分割3 2 求解区域的二角形分割3 3 求解区域的极网络分割 硕 学位论文 3 2 波动方程的有限差分格式 二维均匀各向j 司性搏性动力学方程的标量形式如f p 婆：坠+ 堕 ( 3 1 2 ) p 蒂2 蔷+ 菅 堡：丝+ 丝(313)p= j o + ol j o t o xd z 誓= 以脚i 3 u , + 五警 ( 3 1 4 ) 等= 以伽i 3 u , + 五警 ( 3 1 5 ) 等= 4 警+ 等) 其中虬 和z 表示质点在x 和z 方向的位移，k 、吃和f 。是应力张量，表示质 点在x z 平面内的剪应力。 将位移对时间的二阶偏导用质点振动速度( 位移的一阶偏导) 的一阶偏导代替， 可将式( 3 1 2 ) 降为一阶，即得 口盟：- a - 7 k - - + 一c o r = ( 3 1 7 ) o t d x d z p 盟：錾+ 一c o t = ( 3 1 8 ) - 升g k - ( 2 + 2 p ) 警+ a 鲁班韶 誓= 以脚i 3 v + 五警凹c2珊 拿o t = 4 誓+ 等) 仍， 其中屹= 3 u ，a t 和匕= 3 u ：a 表示质点在x 和z 方向的振动速度， 在x 和z 方向的正应力，吒表示质点在x z 平面内的剪应力。 上述标量方程的矢量形式为 翌：爿竺+ c a w 3 t 融知 ( 3 1 9 ) ( 3 2 0 ) ( 3 2 1 ) k 和f 。表示质点 ( 3 2 2 ) 式中矽= ( v x ，也，氏，屯，屯) 7 ，其中心和v z 为x 和z 方向的速度，皖为面力( 体应力) 分量：a 和c 分别为仅与介质密度和弹性系数有关的系数矩阵。 用交错网格法( 见图3 4 ) 求解式( 3 1 3 ) ，这样v ，和v ，不在源点上计算，避免 了振幅出现无限大的问题，但在同一个节点上，不能同时得到两个分量，与矩形网格 在同一个节点上得到所有分量不同。但正因如此，在交错网格法中使用中心差分时， 避免了矩形网格法中使用中心差分时需要插值而降低精度的缺点。 1 6 硕p 学位论文 幽3 4 夏错嘲霹 不葸圈 如图3 4 所示，设工= i a x 或x = ( f 玲；z = j a z 或z = o g a z ；t = k a t 或 t - - ( k 9 出。其中，血和z 是空日j 采样日】隔，t 是时间采样间隔，则式( 3 1 4 ) 对应的数值计算的差分公式为： ，、 u “5 = u + b v 尝咬r _ + b 。塞忆r 也 ) ( 3 z s ) v 冀一= v 麓一州一， 尝也u r 眦+ ) 州畦+ 尝e 抄。一t 扣) ( 3 2 4 ) r 啪k + l = r ： j + 以+ 2 l + j 尝( u 笔一u ) + + 兰嘛r 嘣， ) ( 3 2 5 ) t 麓j = t j + o + 2 ) f 书兰峨r v 葛 ) + + 尝( i 茂一u 纠 ( 3 2 6 ) h 麓= h ：+ ；+ 心+ 尝( u 譬一u 5 ) + 。+ 老( 、穗+ 一v , q t + ) ( 3 2 7 ) 其中o 和v 分别为u 和叱的离散量；r 、t 和h 分别为f 。、吃和r 。的离散量5 上标k 、下标i 和j 分别为时间t 及空间x 和z 的离散索引；a t 、a x 和z 分别为t 、 x 和z 的离散日j 隔：b 为密摩的倒麴。 3 3 波动方程有限差分解的初始值及震源选择 进行模拟时各节点的初始值为 、 v ，l ，0 = 0 v ：i ，= 0 k b = 0 ( 3 2 8 ) f 。l ，。= 0 f 。1 s o = 0 在弹性波传播的有限差分模拟中，须在网格节点上对应的炮点位置加入震源影 响，零相位子波比其他类型的子波能取得更好的效果。本文采用主频为2 8 h z 的 r i c k e r 子波作为激发震源。其表达式为 1 7 硕十学位论文 f ( o = e 2 c o s ( 2 矿。，) 其中f - 为子波主频。 3 4 波动方程有限差分解的边界条件 3 4 1 自由边界条件的处理 在水平自由衷面上，应力分量必须满足”： f 。= 0 = 0 由上式可得到： 堕：一蔓 a za x ( 3 2 9 ) ( 3 3 0 ) ( 3 3 1 ) ( 3 3 2 ) 盟：- f 上1 盟 ( 3 3 3 ) 出 l a + 2 ，o x 将( 3 3 3 ) 式