（地球探测与信息技术专业论文）复杂地表条件下弹性波数值模拟方法研究.pdf

捅晏 目前，石油勘探正向复杂构造、复杂地表等地区发展，而中国西部复杂地表条件 下的地震勘探是一个世界性难题。复杂地表起伏和风化带对地震数据处理结果的影响 非常大，简单地形的影响可以通过静校j 下很好地消除，而复杂地表地形所产生的散射 和地滚波对地震记录的严重污染及畸变无法通过简单的静校正消除掉。因此研究复杂 地表地震波数值模拟方法对于认识复杂地表下地震波传播的规律、消除地表地形对地 震资料的影响，以及研究基于复杂地表的地震叠前深度偏移成像方法都有非常重要 的意义。本文在前人的基础上，采用了有限差分方法对复杂地表条件下的声波、无反 射递推算法、弹性波进行了数值模拟及传播特征分析。 文中首先对复杂地表条件下的地震波数值模拟各种方法进行了概述。为了分析起 伏地表对直达波、反射波的影响，利用坐标变换原理的坐标映射关系，推导出映射后 坐标系下的二阶标量声波方程，分析了起伏地表对直达波、反射波的能量、波形、振 幅、频谱等参数的影响。另外推导出复杂地表条件下的无反射递推算法，利用逆时偏 移原理，实现了复杂地表条件下的叠后数值模拟及逆时偏移。 对于复杂地表条件下的弹性波数值模拟方面，本文利用不规则起伏自由表面介质 的交错网格有限差分新算法处理起伏自由边界，然后与交错变网格高阶有限差分法相 结合对复杂地表条件下的弹性波方程进行数值求解，此算法不必处理现有的有限差分 必须考虑的其他限制，既减少了机时，又提高了效率。通过典型起伏地表和近地表速 度模型的弹性波方程数值模拟结果，来认识、分析地震波波场特征、地震波主要类型 以及地震波传播规律，从而分析复杂地表对地下反射波场的影响。 关键词：复杂地表交错变网格坐标变换数值模拟逆时偏移 m e t h o ds t u d yo fe l a s t i cw a v en u m e r i c a lm o d e l i n g u n d e rc o m p l e xs u r f a c ec o n d i t i o n s z h a n gh u a ( g e o p h y s i c a lp r o s p e c ta 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 rh a nw e n g o n ga n dp r o f e s s o rl 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 dg a se x p l o r a t i o ni st u r n i n gt ot h ea r e a sw i t hi r r e g u l a rt o p o g r a p h y a n dc o m p l e xg e o l o g i c s t r u c t u r e s ，a n dt h ec o n d i t i o n so f c o m p l e x s u r f a c es e i s m i c e x p l o r a t i o n i sa g l o b a lp r o b l e mi nc h i n a sw e s t e r nr e g i o n c o m p l e xs u r f a c ef l u c t u a t i o na n d w e a t h e r i n gb e l th a v eg r e a ti m p a c to nt h er e s u l t so ft h es e i s m i cd a t ap r o c e s s i n g ，t h ei m p a c t o fs i m p l et o p o g r a p h yc a nb ev e r yg o o de l i m i n a t eb ys t a t i cc o r r e c t i o n ，b u tc o m p l e xs u r f a c e t o p o g r a p h yo ft h es c a t t e r i n ga n dt h er o l l i n gw a v e so fs e i s m i cr e c o r d so fs e r i o u sp o l l u t i o n a n dd i s t o r t i o nc a nn o tb ee l i m i n a t e db ys i m p l ys t a t i cc o r r e c t i o n s ot h es t u d yo fs u r f a c e s e i s m i cn u m e r i c a ls i m u l a t i o nm e t h o df o ru n d e r s t a n d i n gs e i s m i cw a v e p r o p a g a t i o n r e g u l a r i t yo nc o m p l e xn e a r s u r f a c ec o n d i t i o n ，t h ee l i m i n a t i o no fs u r f a c et o p o g r a p h yo nt h e i m p a c to fs e i s m i cd a t a ，a n dt h es t u d yo ft h ep r e - s t a c kd e p t hm i g r a t i o ns e i s m i ci m a g i n g m e t h o db a s i n go f ft h er o l l i n gt o p o g r a p h yh a v ey e r yi m p o r t a n ts i g n i f i c a n c e i nt h i sp a p e r ， b a s i n go nt h ep r e v i o u sw o r k ，u s i n gf d em e t h o dt on u m e r i c a ls i m u l a t ea n dc h a r a c t e r i s t i c s a n a l y s i sf o rt h es e i s m i cw a v ee q u a t i o nu n d e rt h ec o n d i t i o n so fc o m p l e xs u r f a c e f i r s t ，t h i sp a p e ro u t l i n e sav a r i e t yo fs 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 nm e t h o di nt h e c o m p l e xs u r f a c ec o n d i t i o n s i no r d e rt oa n a l y z et h er o l l i n gt o p o g r a p h yo ft h ei m p a c to f d i r e c tw a v ea n dr e f l e c t i o nw a v e ，u s i n gc o o r d i n a t et r a n s f o r m a t i o no ft h ep r i n c i p l eo f r e l a t i o n sb e t w e e nt h ec o o r d i n a t em a p p i n g ，d e r i v i n gt h es e c o n ds c a l a ra c o u s t i cw a v e e q u a t i o nu n d e rt h em a p p i n gc o o r d i n a t e s ，a n a l y z i n go ft h er o l l i n gt o p o g r a p h yo ft h ei m p a c t o ft h ee n e r g ya n dw a v ea m p l i t u d e ，f r e q u e n c ys p e c t r u mo ft h ed i r e c tw a v ea n dr e f l e c t i o n w a v e t h i sp a p e ra l s od e r i v e dt h en or e f l e c t i o nw a v ee q u a t i o nu n d e rt h ec o n d i t i o n so f c o m p l e xs u r f a c e ，a n da c h i e v e dap o s t s t a c kn u m e r i c a ls i m u l a t i o na n dr e v e r s em i g r a t i o n f o rt h en u m e r i c a ls i m u l a t i o no fe l a s t i cw a v eu n d e rt h ec o n d i t i o n so fc o m p l e xs u r f a c e ， t h i sp a p e rh a sab a n do fe l a s t i cw a v ee q u a t i o nn u m e r i c a ls o l u t i o nu s i n gt h ev a r i a b l eg r i d a n dt h es t a g g e r e dg r i dw i t hh i g h o r d e rf i n i t ed i f f e r e n c em e t h o d ，t h e nu s i n gt h ei r r e g u l a r u n d u l a t i n gs u r f a c eo ft h em e d i u ms t a g g e r e dg r i df i n i t e d i f f e r e n c ea l g o r i t h md e a l sw i t ht h e f l u c t u a t i o nf r e eb o u n d a r y t h r o u g ht h ee l a s t i cw a v ee q u a t i o nn u m e r i c a ls i m u l a t i o nr e s u l t s o ft h et y p i c a lf l u c t u a t i o ns u r f a c ea n dn e a r - s u r f a c ev e l o c i t ym o d e ，u n d e r s t a n d i n ga n d a n a l y z i n gt h ec h a r a c t e r i s t i c so fs e i s m i cw a v ef i e l d ，m a i nt y p e so fs e i s m i cw a v e sa n d s e i s m i cw a v ep r o p a g a t i o n ，t h u st h en e a r - s u r f a c ea n a l y s i so ft h ei m p a c to fr e f l e c t e dw a v e f i e l d k e yw o r d s ：c o m p l e xs u r f a c e ，v a r i a b l es t a g g e r e dg r i d ，c o o r d i n a t et r a n s f o r m a t i o n ， n u m e r i c a ls i m u l a t i o n ，r e v e r s et i m em i g r a t i o n 关于学位论文的独创性声明 本人郑重声明：所呈交的论文是本人在指导教师指导下独立进行研究工作所取得 的成果，论文中有关资料和数据是实事求是的。尽我所知，除文中己经加以标注和致 谢外，本论文不包含其他人己经发表或撰写的研究成果，也不包含本人或他人为获得 中国石油大学( 华东) 或其它教育机构的学位或学历证书而使用过的材料。与我一同 工作的同志对研究所做的任何贡献均己在论文中作出了明确的说明。 若有不实之处，本人愿意承担相关法律责任。 学位论文作者签名：蔓迭亡仁 日期： 神g 年 月 5 日 学位论文使用授权书 本人完全同意中国石油大学( 华东) 有权使用本学位论文( 包括但不限于其 印刷版和电子版) ，使用方式包括但不限于：保留学位论文，按规定向国家有关 部门( 机构) 送交学位论文，以学术交流为目的赠送和交换学位论文，允许学位 论文被查阅、借阅和复印，将学位论文的全部或部分内容编入有关数据库进行检 索，采用影印、缩印或其他复制手段保存学位论文。 保密学位论文在解密后的使用授权同上。 日期：葛寸形年参月r 日 日期：一辞6 月产日 中国石油大学( 华东) 硕十学位论文 第一章前言 消除地表地形对地震资料的影响一直都是地震勘探中的一个世界性的难以克服 的问题，尤其是在中国。众所周知，中国的地形地势是非常复杂的，但是在当今对油气 资源大力需求的形式下，地震勘探的前缘正在逐步向复杂地表地区转移，尤其是在西 部的沙漠、高原、山地和丘陵地区，地震工作已经大面积的展开了。所以对这些地区 的地表地形影响的研究就显得日益迫切。一般情况下，对地表地形影响的消除主要是 用静校正方法。但是这里有两个问题需要解决：其一是静校正量的确定，对于复杂的地 表条件，这是一个很难确定的问题，需要做大量的正演模拟和野外调查；另外一个致命 的因素就是，目前的静校f 方法是以几何地震学的假设前提( 地表水平，表层介质均匀 为前提的) ，显然复杂地表地形的影响是不在这个假设之内的。 近年来，地震勘探所涉及的复杂地表包括近地表的不均匀性( 引起近地表散射) ， 起伏地表、厚的黄土覆盖层，地形起伏变化复杂，地表高差大( 山地) 。山地的特点是地 形剧烈起伏，岩石出露，山间地表松软，低速层厚度变化大。山前戈壁砾石区近地表速 度变化大，高速岩层直接出露地表( 不存在低速带，速度横向变化大，近地表射线传播 路径不垂直) 复杂地表使速度的提取精度受到影响，另外复杂的地表地形散射也影响 了深层信息的提取。例如世界上公认的地震勘探困难地区黄土高原，主要表现在两个 方面：巨厚的黄土，干燥疏松，地震波吸收衰减严重，激发和接收条件差；表层结构复杂， 地形高差大，黄土层厚度和低速带厚度、速度横向变化剧烈，静校正问题突出，测线布 设以沿沟弯线为主。黄土塬地区，一方面是激发能量弱，另一方面是黄土对地震波的吸 收衰减严重。相对于白垩系砂泥岩而言，黄土层对地震波的吸收是严重的。比较同 检波点，接收在基岩与黄土中激发折射波初至的能量，含膏砂泥岩是黄土的5 l 倍， 相差1 4 ，- - - 2 4 d b ，地震波在黄土层中的衰减是基岩的2 8 倍。结论是在黄土上激发区 无论什么埋置条件，资料总是很差，红土的激发能量比黄土强。 众所周知，在中国西部复杂地区实施地震勘探是一项集高难度、高技术、高投 入和高风险于一身的极富挑战性的工作。要在复杂的地表起伏条件下获得较高品质的 地震剖面，一方面要克服恶劣的地表条件给野外施工作业带来的影响；另一方面又必 须较好地解决和攻克在地震勘探方法和技术上因复杂地表因素而诱发的各种难题，如 地震资料的信噪比很低，干扰类型复杂，产生的散射和地滚波对地震记录的严重污染 和畸变无法通过简单的静校f 去消除，速度分析、叠加成像和偏移归位及构造成像困 第一章前言 难等。随着油气勘探难度的加大( 如复杂地表条件) ，传统的基于平缓构造勘探的地震 数据采集、处理和解释方法技术面临一系列挑战，许多技术在这类复杂地表地区已经 不再适用。 地震数值模拟是地震勘探和地震学的重要基础。地震数值模拟就是在假定地下 介质结构模型和相应物理参数己知的情况下，模拟研究地震波在地下各种介质中的传 播规律，并计算在地面或地下各观测点所应观测到的数值地震记录的一种地震模拟方 法。这种地震数值模拟方法已在地震勘探和天然地震领域中得到广泛的应用。它不但 在石油、天然气、重金属和非金属等矿产资源及工程和环境地球物理中得到普遍的应 用，而且在地震灾害预测、地震区带划分以及地壳构造和地球内部结构研究中，也得 到相当广泛的应用。 ，通过地震波传播数值模拟分析，弄清地震波在复杂地表地形中的传播特征，可以 根据不同地区的实际情况，采用更加合理的野外观测系统( 如测线方向和排列长度) ， 使用更加合理的处理流程和处理参数，进行合理的地震解释，以及研究基于复杂地表 的地震叠前深度偏移成像方法都有非常重要的意义。 在复杂地表条件下的地震波传播数值模拟方面，两个关键问题是如何选择适合不 规则地表的数值计算方法以及如何更好地实现自由边界条件，实际模拟中两个问题也 是有机联系在一起的。 本文在前人研究的基础上，以声波、弹性波波动方程理论为基础，利用坐标变换 原理对起伏地表条件下的声波方程进行数值求解，并进行特征分析；推导出适合做逆 时偏移的无反射递推算法，结合坐标变换原理，对复杂地表条件下的模型进行叠后数 值模拟及逆时偏移；而对于复杂地表条件下的弹性波方程，采用交错网格高阶有限差 分数值计算法求解，利用不规则起伏自由边界新算法处理起伏边界，同时利用交错网 格变步长法对起伏边界进行加密采样处理，对复杂地表条件下的弹性波进行了数值模 拟及传播规律分析。 中国石油大学( 牛东) 硕十论文 第二章复杂地表条件下数值模拟概况 众所周知，在中国西部复杂地区实施地震勘探是一项集高难度、高技术、高投入 和高风险于一身的极富挑战性的工作。要在复杂的地表起伏条件下获得较高品质的地 震剖面，一方面要克服恶劣的地表条件给野外施工作业带来的影响；另一方面又必须 较好地解决和攻克在地震勘探方法和技术上因复杂地表因素而诱发的各种难题，如地 震资料的信噪比很低，干扰类型复杂，产生的散射和地滚波对地震记录的严重污染和 畸变无法通过简单的静校正去消除，速度分析、叠加成像和偏移归位及构造成像困难 等。而传统的基于平缓构造勘探的地震数据采集、处理和解释方法技术面临一系列挑 战，许多技术在这类复杂地区已经不再适用。因此，研究复杂近地表地震波传播问题 数值模拟方法，对于认识复杂地表下地震波传播的规律，弄清地震波在复杂地表条件 下的传播特征，以及研究基于起伏地表的地震资料处理、采集和解释方法都有非常重 要的意义 2 1复杂地表数值模拟定性分析 复杂地表的类型很多，对地震资料的影响程度也不尽相同，要想完全描述近地表 的影响因素是不可能的。为了模拟地表严重起伏和表层强非均匀性对地震波传播影响， 算法的选择要符合以下几点：( 1 ) 不平的自由地表面必须在几何上精确描述并且自由 地表面容易用算法实现；( 2 ) 所有前向、反向散射，包括在局部复杂速度构造层内和强 不均匀介质内的鸣震，必须在模型正演过程中得到正确处理。 对于地球中地震波的传播来说，地震数值模拟是一项重要的技术。其目标就是在 假定地质构造已知的情况下，预测能够在地面接收到的理论地震图。对于地震解释和 地震反演方法来说，这也是一种非常重要的技术。其另一个比较重要的方面是为地震 勘探提供评估和设计。地震波数值模拟方法有很多种，可以分成3 类：直接法、积分方 程法和射线追踪法，各种方法有其优缺点。 直接法就是解波动方程，即将地质模型近似放在一个数值网格上，也就是个有 限数字离散点的模型。这些方法也可以叫做网格法或全波动方程法，后者的叫法是因 为方程的解给出了全波场的特征。直接法对于介质的的变化没有什么限制，并且在足 够小的网格上能够得到非常精确的解。这种方法能够处理含不同流体的问题，而且非 常容易得到波场传播快照图( 辐射图) ，这对结果的解释非常重要。但其缺陷就是相对 于其他方法来说，耗机时较多。 第一：市复杂地表条件下数值模拟概况 积分方程法是基于点源波场波级数的积分表达。该方法起源于惠更斯原理，详细 研究惠更斯原理可以认为：波场一方面可以表示为一系列点源所形成的波场的叠加， 另一方面可以表示成由边界点源所形成的波场的叠加。惠更斯原理现在仍然在使用， 称这两种方式为体积分方程和边界积分方程，并且各有自己的适用范围。这两种方法 在应用上比之于直接方法有更多的限制。然而对于特定的几何形态，例如均匀介质中 包含异常绕射体问题、井l 日j 问题及微裂隙等问题，积分方程法的计算效率和精度都很 高。由于具有良好的解析性，该方法对于建立在b o r n 近似基础之上的成像方法也很用。 渐进法或射线追踪法是地震模拟和成像领域中应用最广泛的技术。根据射线理论 确定射线轨迹和波传播时间，用波动方程理论确定波的振幅。这种方法快速简单，得 到的波的传播时问等运动学特点较准确，但振幅等动力学特点则是不精确的：且由于 没有考虑全波场特征，所以是一种近似方法。但是从另外一个角度看，这种方法也是本 文所讨论的方法中效率最高的方法。 所以在上述方法中，能够适应复杂地表地质条件的模拟方法一般选取直接法( 波 动方程方法) 。这是由实际问题的要求和该方法本身的特点所决定的。 2 2 复杂地表正演模拟定量研究 波动方程的解法有很多，总体上以网格法为主，包括有限差分、有限元、傅立叶变 换等方法及这三类方法的联合方法等。其基本特点是能够提供的波场快照图，可以帮 助研究者细致分析波场在各种复杂介质内部的反射、透射、绕射、散射以及能量的衰 减等运动学和动力学的各种细节特征。 由于实际工作中所模拟的介质不同，所用的模拟方程也不一样。根据模拟方程的 不同，波动方程数值模拟主要有：声波模拟、弹性波模拟、粘弹性波模拟以及裂隙和孔 隙弹性模拟等。波动方程模拟方法的数学实现是数值模拟的一项重要研究内容，它需 要极高的计算机资源，这也是长期制约该方法走向生产实际的最主要因素，并因此产 生了针对各种不同地质问题的数值实现方法。例如：有限元( f e ) 、有限差分( f d ) 、 f o u r i e r 变换、边界元( b e ) 、谱元( s e m ) 、广义屏等方法，但这些方法都各具优缺点。 ( 1 ) 有限元法( f e ) 是基于变分原理和剖分插值，考虑的是分段近似，比较适宜于 模拟任意地质体形态，可以任意三角形逼近地层界面，保证复杂地层形态模拟的逼真 性：但算法复杂，计算速度慢，一般对网格要求三角剖分，基函数是分段线性函数， 不具有f 交性。算子也是空 日j 局部的，空间分辨率高，但是频率域中分辨率低，占用 内存和运算量均较大。但是此方法处理不规则地表比较方便，成为模拟起伏地表情况 4 中国石油大学( 华东) 硕十论文 下地震波传播的有效方法。为此，发展了些有限元和其他方法结合的混合方法，例 如：m o c z o 等( 1 9 9 7 ) 用离散波数方法模拟震源激发和下部介质中地震波的传播，而 通过有限元方法来模拟沿起伏地表波的传播( d w f e ) ：黄自萍等( 2 0 0 4 ) 提出了一种 有限元和有限差分( f e f d ) 耦合的区域分裂方法，它灵活的网格剖分，克服了单纯 用差分方法对区域的依赖性，而同样的精度所需的计算量比有限元法小，利用这种方 法成功地模拟了起伏地表地震波传播。马德堂等( 2 0 0 4 ) 用有限元法和伪普法相结合 的方法模拟起伏地表模型。用有限元法计算起伏地表区域和模型第一区的波场值，而 利用伪谱法计算模型内部区域的波场值，对于第二区采用“过渡区域 技术解决两种 不同算法的衔接问题，这样使得自由边界条件能自动满足，人工边界的反射也容易消 除，既提高了模拟精度也节省了计算量。但是这类混合方法的主要问题是有限元法和 其它方法计算区域交界处易产生人为反射，吸收边界也不易处理。 ( 2 ) 傅立叶变换法( 伪谱法、虚谱法) 是一种逼近空间微分的方法，基于空间域 中的求导，相当于频率域中的乘积运算，利用傅里叶变换将波场函数表示为傅旱叶级 数的展开形式，将波动方程在时间一波数域或频率域中求解，精度高，占用内存小；但 是由于傅里叶变换是基于整个时间域或空间域的，改变空间中的某一点的值，就会改 变频率域中的所有值，因此每一点的微分结果都要受到计算域中其它点的影响并且存 在众所周知的g i b b s 效应。实际上，求导运算应该是一种局部运算。对于空间物性剧 烈变化的情形，这显然是不适合的，即在处理地表结构复杂或地表剧烈起伏以及地下 结构复杂的情况时存在较大的难度。例如：n i e l s e n ( 1 9 9 4 ) 、赵景霞等( 2 0 0 3 ) 将伪 谱法用到了曲线网格的二维声波方程的模拟中，提高了运算效率；t e s s m e r 等( 1 9 9 4 ) 通过坐标变换的思路，利用傅立叶变换法与有限差分法相结合的方法对起伏地表的弹 性波波场进行了模拟。 边界积分或边界元法( b e ) 是一种发展迅速的计算方法，此方法将散射波场通过 地表的一个半解析的积分来表示，其中积分项中g r e e n 函数一般在频率波数域中计 算。这种方法在研究起伏地表地震波传播时使用较多，例如，b o u c h o n ( 1 9 9 6 ) 、p e d e r s e n ( 1 9 9 4 ) 、d u r a n d ( 1 9 9 9 ) 等就是用这种方法研究了地形的区域变化对天然地震波场 的影响，分析了地表散射以及地表振动随地形的变化，而r u - s h a nw u 等( 1 9 9 8 ) 将 边界元和屏传播算子相结合提出了一种混合模拟算法，在处理起伏地表的同时提高了 模拟效率。b e 方法的半解析性质决定了该方法不能适用于地表速度变化较大情况， 而实际情况是，由于后期地质作用造成浅部地层速度变化更为剧烈，因此限制了b e 第二章复杂地表条件下数值模拟概况 方法的实际应用。 ( 3 ) 有限差分法( f d ) 是通过有限差分算子将波动方程离散化，以差分代替微分， 将微分方程问题化为代数问题，然后求解相关的线性代数方程组以获得微分方程问题 的数值解。差分法是点近似，只考虑有限个离散点上的函数值，而不考虑点的领域函 数值如何变化。差分算子是一个空间局部的算子，在空间域具有较高的分辨率，可以 很好地适应剧烈变化的地下介质情况。但是在频率域中，有限差分算子的分辨率却很 低，即精度低，仅适合于相对较简单的地质模型。但算法的计算速度快，占用内存小， 在模拟复杂介质中地震波传播时应用最为广泛，但该方法的一个主要缺陷是处理复杂 地形比较困难。为此，t e s s m e ra n dk o s l o f f ( 1 9 9 2 ) 提出了一种新的思路，即通过 坐标变换将具有起伏地表的模型及弹性波方程变换到新的具有水平地表的坐标系中， 在新坐标系中求解弹性波方程，时间上用差分法，空间上横向用傅立叶法，纵向采用 c h e b y s h e v 方法计算波场对空间的导数，使用的吸收边界和自由边界均由g o t t l i e b ( 1 9 8 2 ) 提出的边界条件。h e s t h o l ma n dr u u d ( 1 9 9 4 ，1 9 9 8 ，1 9 9 9 ) 借鉴了t e s s m e r 等的思路，通过坐标变换将起伏的地表转换成水平地表后，完全通过差分法求解弹性 波方程，在适应地表起伏的同时提高了计算效率。s t i g h e t h o l m 等( 1 9 9 8 ，2 0 0 1 ) 用交错 网格高阶有限差分法解速度一应力方程，研究波场p 波在起伏地表处的散射波和转换 波，关键技术是用仿射坐标变换将直边梯形网格变为直角网格，将起伏地表的自由边 界条件作相应变换，他们已经证明了该方法是条件稳定的。甘文权( 2 0 0 1 ) 利用类似 方法对位移波动方程进行了数值求解。裴正林( 2 0 0 4 ) 运用交错网格一阶空间导数的 任意偶数阶精度展开式和相应差分系数计算式以及一阶双曲型应力一速度弹性波方 程交错网格任意偶数阶精度差分格式求解方程；然后采用将零速度法和广义虚像法相 结合的方法来处理自由边界，并在自由边界上采用四阶精度差分格式，成功地对任意 起伏地表弹性波波场进行了模拟。 ( 4 ) 近年来提出的谱元法( s e m ) 把有限元法和谱展开法相结合，它的基本思想是 选取以正交多项式表示的形函数，在各个单元上通过配置点插值，提高解的级数表示 的收敛速度。这种方法把计算范围分成许多子区域，在每个子区域中把解表示成截断 的c h e b v s h e v 或勒让德多项式级数的乘积，用伽辽金方法求解正交问题的变分格式， 得到离散方程。谱元法是弹性波方程在空间近似的高阶有限元方法，它可以通过选取 适当的形函数，使得离散方程中等效质量矩阵对角化，从而大大降低了计算成本。而 且由于具有有限元法处理不规则地表比较方便的特征。所以在处理起伏地表的数值模 中罔油人学( 华东) 硕十论文 拟方法中谱元法是一种比较有前途的方法。例如：d i m i t r i k o m a t i s c h 等( 1 9 9 8 ，2 0 0 0 ) 用谱元法实现了起伏地表二维和三维弹性波数值模拟，讨论了复杂地表引起的各种干 扰，该方法具有费时少、精度高等优点，并利用该方法模拟了按正弦规律起伏的海底 反射、透射特点。 这些数值解法的基础都是在离散网格上计算基于各种方法实现的波动方程代数 方程或方程组。因此，在解决复杂地表问题的上，基于网格化这个层次上可以分为两 类：一类是从模型出发，改变方程，建立适应这种情况的方程，如映射网格，可调节网格， 坐标变换等；另外一种方法就是建立适当的边界条件。 在处理起伏界面问题时，离散化数值计算只能将连续界面变成阶梯状界面，在大间 距的网格上做各种数值计算必然会产生不必要的人为绕射和由于界面条件不准而引 起的积累误差。为了消除这些绕射和误差，产生了许多方法。 2 2 1 映射网格法 在曲网格上做数值计算的思想最早由f o r n b e r g 提出，但是f o r n b e r g 的映射技术只 能使网格沿着主要的界面变化，n i e l s e n 等采用了计算流体力学中常用的分块映射技 术使网格可以沿着所有界面变化，提高了算法的精确性。使数值网格沿着所有地下界 面变化是n i e l s e n 的主要思想。为实现这一目标，可将实际的几何构造影射成某一曲坐 标系下的规则图形，比如矩形。再利用插值法产生区域内部的各个点，从而生成所需的 曲网格( 如图2 2 一1 ) 。 n 如 z o f 。r i ) i l 墨q 憾 图2 。2 - 1 直网格变曲网格示意图 具体实现如下：将直角坐标系x z 下的任意一个不规则多边形的边分成4 段，每 段对应于四边形中的一条边。将2 条相对的边分别映射成曲坐标系孝一7 7 域下的2 条平 行于7 7 轴的对边，类似地将另2 条相对边映射成益坐标系下的2 条平行于f 轴的对边。 通过这4 条边界，任意一个直角坐标系的不规则区域可映射成了一个矩形，相应地建立 了一个曲坐标系，然后在曲坐标系下建立一个规则网格，就可以在规则网格上进行数 7 第二章复杂地表条件下数值模拟概况 值计算了。一般来说网线常常取成整数，如善= 1 ，2 ，r - - 1 ，2 ，曲坐标系下网点的 坐标为善，7 7 ，而其实际坐标值为j - - x ( # ，r ) ，z = z ( 善，r ) 。图形内部点实际坐标值采 用一种多维插值法一超限插值法来计算，即利用4 条边界上的点产生内部点的一种插 值方法。 为了使网格沿所有的地下界面变化，并防止因界面过于复杂而造成网格过于不规 则或过于倾斜。采用在计算流体力学中应用广泛的分块方法，每一层被分成若干相邻 的块，将每一块单独映射生成网格后再把所有的块连起来形成统一的整体。连接时应 注意网格的连续性。因为数值运算是在曲坐标系下的规则网格上进行的，所以应把直 角坐标系的波动方程变换成曲坐标系下的方程。这样映射后在曲坐标系下的网格是一 个规则网格，波动方程就可以在这样的网格上计算。 2 2 2 非规则网格法 一 非规则网格差分算子与上面的分段映射方法有相似之处，通过在离散混合变量弹 性波方程时引入非规则网格差分算子，应用于应力一速度方程可以适用于任意形状的 边界和内部边界。该方法的本质就是一个精确处任意形状自由表面的边界条件的坐标 变换问题。这个方法从任意四边形网格着手，在非正交网格上定义任意四边形网格差 分算子，由任意四边形的4 个角点的函数值计算函数在域内的空间导数，将( x ，y ) 坐标 系下的任意四边形映射到( f ，7 7 ) 坐标系下的标准四边形，则( x ，y ) 坐标系下任意四边 形网格内的差分( 微分) 计算可以转化为在( 考，7 7 ) 坐系下的正方形网格上进行，定义映 射函数为： i z = 去( 1 均( 1 刊”去( 1 咱( 1 训矿去( 1 一期刊墨+ 丢( 1 均( 1 刊托 i z = 去( 1 + 善) ( 1 + 7 7 ) z 。+ 去( 1 一d ( 1 + 刁) z ：+ 去( 1 一d ( 1 一，7 ) z ，+ 去( 1 + 9 ( 1 一叩) z 。 其中( x j ，z ，) ( i = 1 , 2 ，3 ，4 ) 是( x ，z ) 坐标系下四边n 4 个角点的坐标。该映射函数可以将图 中左侧的任意四边形映射为右侧的边长为2 的标准正方形，若在( 孝，7 ) 坐标系下对函 数尸( x ，z ) 求导，则有： 一a p ：塑鱼+ 塑鱼 望：望鱼+ 望鱼 a 专i 坟a 考8 za弓8，1 8 xa qa z8 研 由上式可以得到： 中国石油火学( 仁东) 硕j j 论文 瓠8 z 0 4a 专 o xo z o r a 7 7 设函数尸( 爿，z ) 在四边形4 个角点处的函数值p 。( 卢1 ，2 ，3 ，4 ，则( 孝，1 1 ) 坐标系下函数 尸( 置z ) 可以用如下双线性插值描述： p = 三( 1 + 4 ) 0 + r 1 ) p 。+ 去( 1 一善) ( 1 + ，7 ) p ：+ 三( 1 一舌) ( 1 7 7 ) p ，+ 三( 1 + f ) ( 1 一巧) p 。 对上式求导，得： 印 鸳 印 o r 1l1 + 7 71 7 7l + 7 71 7 7l 41 + 孝1 一善1 + 善1 一亭l 将这个式子代入上面的矩阵中，即可得到任意四边形4 个角点的函数值，计算函数 在域内空间上导数的差分算式。类似的方法还有尧德中等利用六角形网格的周期采样， 通过线性变换直接利用傅立叶变换方法求解波动方程，也可以得到处理复杂边界条件 的目的。 另外由m a r t i n ，k a s e r 等人提出的不规则三角网格，可以任意剖分不规则边界，通 过一个不规则加权函数来离散波动方程，来求解复杂条件下的波动方程的解，包括复 杂地表条件。 其实改变网格形状或大小的不规则网格还有很多种，例如最小网格、不均匀网格、 力口权网格等很多( d o n g j o omi n ，c h a n g s o os hi n ，b y u n g d o ok w o na n ds e u n g h wa l l c h u n g ，2 0 0 0 ) 。在有些情况下，这样的剖分方式更有助于算法的实现和精度的保证。 2 2 3 坐标变换法 董良国( 2 0 0 5 ) 利用纵向坐标变换思路，将不规则计算区域转变为新坐标系下的 规则计算区域，就可以使用统一的差分进行计算，并在零应力自由边界概念基础上， 导出了新坐标下的自由边界条件，为利用交错网格差分法提高了方便，同时将特征分 析法吸收边界条件推广到起伏地表情况，利用速度一应力方程模拟了地表起伏情况下 弹性波的传播。具体实现思路为： 设地表高程变化z 。= 厂g ) 为一个单值函数，最高点高程为z 。，将x o z 坐标系作 如下线性变换： 第二节复杂地表条件下数值模拟概况 j 善3 乏( o x x 。a x ) 卜焉z p 一一“一7 通过这样的纵向拉伸变换，原来顶界面起伏的不规则计算区域( 图a ) 就变换为水平顶 界面的规则计算区域( 图b ) ，曲网格也相应地变换为矩形网格，为进行有限差分数值 计算提供了方便。 图2 2 - 2 纵向拉伸变换示意 在上述纵向坐标变换。下，x o z 坐标系f 用速度一应力表不的二维弹性波方程： 一o u ：型士一一o u a b = 士 a t瓠a z 在新坐标系乒7 下变换为： 豢：a 豢+ p 皓归尝+ q g ，7 7 p 等 = + p l fl d + 口i f ，刀i 爿 a t a 芒 i 。 a n。“8 玎 式中u = u e , u , “劈，甜盯叩，“勃) 7 为乒7 坐标下速度一应力构成的列向量， p 皓) 2 考蒿和q ( 亭，刁) = 一7 茜厂( 孝) 为举标拉伸系数，f ( f ) 2 善为x = 孝处地形坡度。 上式实质是线性变换，决定了乒7 7 坐标系下物理量u 和x o z 坐标系下物理量u 是一一对应的映射关系。因此，只要在规则坐标系下求得每个网格点上的 u 。= 0 f ，“。，“菇，甜。，“却) 7 ，就可以确定弯曲坐标系x o z 下每个网格点上的物理量 u 。= 0 ，“：，“。，都。，“。) 7 。上述思路实际上是通过复合求导来数值计算x o z 坐标系下 弯曲网格点处波场的空间导数。 2 2 4 混合网格法 黄自萍、董良国等人发表的弹性波传播数值模拟的区域分裂法( 2 0 0 4 ) 针对弹性 波波场传播的数学模型，提出了一种给予有限元和有限差分耦合的区域分裂方法，而 中国杠油人学( 华东) 硕 j 论文 网格剖分也是基于这两种算法上的混合网格，以此来处理起伏地表的介质模型。 将计算区域q 分成部分重叠的两部分q 和q ：，如图2 2 3 所示为了求解方便， 设q 的上边界在靠近左右两边界的地方是水平的，q 。是由地表及边界r l 包围起来的 不规则区域，q ，是山人工边界和边界r ，以及部分水平地表包围起来的凹字型区域。 为了网格剖分的方便，令d is t ( r i ，r ) = 2 k h ( k 是大于等于1 的整数，h 为差分计 算所需的空间网格步长) 。由r 1 ，r ，和部分水平地表包围起来的区域，我们称之为q ， 显然有q nq ：= q ，。 首先，在q ，上作四边形网格剖分：先将i ，等分，假设等分后的点从依次为 只2 ，露，只2 ，焉，那么有d i s t ( 只2 ，p , i ) = 2 h ；然后在q 。，q 3 上以p 2 ( 1 i n ) 为区域q ，q ，部分边界1 1 ，的节点从作四边形剖分；并且将区域q ，划分为以2 h 为边 长的正方形组成的区域这样就得到了在整个q ，区域上的正规四边形剖分。最后，在 q ，上做步长为h 的正方形网格剖分。 网格剖分好后在区域q ，上既有正方形又有四边形的网格剖分，且正方形网格的 节点落在四边形的顶点或中点上，如图2 2 3 所示。 霄强廖锻猕忑熬千 * 灌】一 7 ：i ：一：卜挎 ；r 1 1 t1 ltl iiii ll ；l土! 1； o 工强。 j ；1 毒1 e； i i 图2 2 3 区域分裂示意图和网格剖分示 2 2 5 不同阶差分法 用波动方程有限差分方法做正演模拟时，如果使用二阶差分格式，网格间距必须 取得很小，以保证计算精度和稳定性，而这样做明显地降低了计算效率。为此，d a l a i n 矛 1 m u f t i 提出了高阶有限差分方法，可以适用于复杂介质强横向变速情况下的正演模 拟。在高阶方程的情况下，网格间距可以取得很大，但计算精度并不l l - - 阶差分方程 小网格间距时低。而网格阳怔巨增大可以显著提高计算效率。一般称截断误差高于四阶 的差分方程为高阶差分方程。 第_ 二章复杂地表条件下数值模拟概况 三维声波方程：筹+ 矿0 2 u + 可0 2 l g 2 瓦击萼o t ，其中吣戊z ) 为地表记录的 苏2 却2 瑟2v2 ( x y ，z 1 2 。 一 一一一 压力波场，v ( x ，y ，z ) 为纵横向可变化的介质速度。在局部速度变化很缓或不变时，截 断误差o ( 血m ，a y m ，心m ) 的统一的三维正演模拟高阶差分方程是： “( t + a t ) = 2 砸) 叫卜+ ( v 2 窘+ 窘+ 鲁o z 出咖 + 未( v m 4 【 矿0 4 u + 萨o a 4 + 可q 0 4 u j _ + ( , - a t ) 4 【翥+ 意+ 急】+ o 池4 ， 这是一个推导用于正演模拟高阶差分方程的起始方程。它在时间方向上截断误差 为o ( a t 4 ) ，而空间的截断误差则要根据需要来定( 至少是四阶以上的) 。另外，可以 根据需要来组合不同阶次的差分格式。当局部速度横向变化很大，即v = v ( x ，少，z ) 时， 截断误差为o ( 血m ，a y 盯，a z m ) 的统一的三维正演模拟高阶差分方程为： 肌卅砌”舡纠+ 陋) 2 ( 警+ 雾+ 窘) + 等 ( 芸) 2 + v 窘+ ( 多 2 + v 雾+ ( 翁+ v 象1 ( 窘+ 雾+ 万02 u v3 f 4laf ，3 “矽3 甜 夕3 材1夙f ，汐3 扰矽3 扰d 3 “、 + 了i 瓦l 万+ 丽+ 蕊肟【、丽十矿+ 面j + 氧急+ 而03 1 d + 针1 4 丛o k 2 & ：+ 嘉+ 舄 + 瓦i 蕊+ 丽i + 可j + 石【比叫+ 万万+ 矛虿j + 击( v 叫窘+ 雾+ 雾 + 。p 4 ) 其中十阶差分方程的系数是： 当m = 1 0 时，( i ) o = - 5 8 5 4 4 4 4 5 ，l = 3 3 3 3 3 3 3 ，0 ) 2 = 一0 4 7 6 1 9 0 1 ， o 3 = 0 0 7 9