（地质工程专业论文）任意差分精细积分逆时偏移在复杂构造成像中的研究与应用.pdf

摘要 本文在高精度任意差分精细积分地震模型正演研究的基础上，将任意差分精细积 分应用于逆时偏移中，提出了任意差分精细积分逆时偏移的方法。推导了任意差分精 细积分逆时偏移的公式，给出了差分格式的稳定性条件。该方法可以适用于不规则网 格划分，实现局部点的加密。由于逆时偏移采用全波动方程，没有做近似，并且对方 程的求解采用的是高精度任意差分精细积分，因此从理论上此方法是一种精度较高的 差分偏移方法。由于没有对速度的限制，且得到的是深度剖面，因此属于深度偏移的 范围。数值模拟对各种复杂地质构造的偏移都得到了较满意的结果，表明了精细积分 逆时偏移应用于复杂构造的可行性。并且表明此方法具有空间和时间精度高，能压制 多次波和速度扰动，不受复杂构造地层的限制，受随机噪音干扰小的特点。与其他方 法的对比表明任意差分精细积分逆时偏移成像精度优于其他方法。但任意差分精细积 分逆时偏移计算比较费时，且偏移的效果与所给定的速度模型的精确程度有关。最后 本文对任意差分精细积分的不足做了分析，对其研究方向给出了建议。 关键词：逆时偏移，任意差分，精细积分，复杂构造，波动方程 s t u d ya n da p p l i c a t i o no na r b i t r a r yd i f f e r e n c ep r e c i s ei n t e g r a t i o n r e v e r s e t i m em i g r a t i o ni nc o m p l e xs t r u c t u r ei m a g i n g l i ub a o m i ( g e o l o g i c a le n g i n e e r i n g ) d i r e c t e db yp r o f s u nc h e n g y u a b s t r a c t b a s e do nh ig h - a c c u r a c yf o r w a r dm o d e l i n gm e t h o do fa r b i t r a r yd i f f e r e n c ep r e c i s e i n t e g r a t i o nm e t h o d ( a b b r e v i a t ea d p i ) ，t h i sp a p e ra d a p t sa d p it ot w o w a yr e v e r s e t i m e m i g r a t i o na n db r i n g sf o r w a r dt h ea r b i t r a r yd i f f e r e n c ep r e c i s ei n t e g r a t i o nr e v e r s e t i m e m i g r a t i o nm e t h o d t h ec o r r e s p o n d i n gd i f f e r e n c ef o r m u l ai sd e r i v e da n dt h es t a b l e c o n d i t i o ni sp r o v i d e d 眦sm i g r a t i o nm e t h o dc a nf l e x i b l yd e a lw 曲u n e q u a l 鲥dp l o t ， t h e r e b yi tc a nt h i c k e nl o c a lp o i n t s b e c a u s et h er e v e r s e t i m em i g r a t i o ni st w o w a yw a v e e q u a t i o nm i g r a t i o na n dt h ea d p ii st h eh i g h - a c c u r a c ys o l u t i o nt ot h ep a r t i a ld i f f e r e n c e e q u a t i o n ，t h i sm e t h o di s o n eo ft h eh i g h a c c u r a c ym i g r a t i o nm e t h o d st h e o r e t i c a l l y f u r t h e r m o r e ，i th a sn ol i m i t a t i o nt ov e l o c i t ya n dt h eo u t p u ti sd e p t hp r o f i l e ，t h ea d p i r e v e r s e t i m em i g r a t i o ni so n eo fd e p t h m i g r a t i o nm e t h o d s a l lk i n d so fg e o l o g i c a l s t r u c t u r e s s u c c e s s f u li m a g i n gp r o v e dt h ef e a s i b i l i t yo fa d p ir e v e r s e - t i m em i g r a t i o na n di t h a st h ea d v a n t a g e so fh i 曲一a c c u r a c y , 谢t l la b i l i t yo fs u p p r e s s i n gm u l t i p l ea n dv e l o c i t y d i s t u r b ，a n di ss u t i a b l ef o ra n yc o m p l e xs t r u c t u r ea n dn o ts e n s i t i v et or a n d o mn o i s e c o m p a r i n gw i n lo t h e rm i g r a t i o nm e t h o d si n d i c a t et h a ta d p ir e v e r s e - t i m em i g r a t i o ni sa m i g r a t i o nm e t h o d 、析t l la d v a n t a g eo fh i g hi m a g ea c c u r a c y i nt h es a m et i m eb e c a u s eo fi t s 址g h - o r d e rs p a c ea p p r o x i m a t i o na n da l lw a v ef i e l dc o m p u t a t i o n ，t h i sm e t h o di st i m e c o n s u m i n g ，a n dn u m e r i c a lp r a c t i c eh a si n d i c a t ei t si m a g i n gi sh i g h l yr e l a t e d 、析t l lt h eg i v e n v e l o c i t y l a s t l yt h ed e f i c i e n c yo fa d p ir e v e r s e t i m em i g r a t i o ni sa n a l y z e da n dt h e s t u d y i n gd i r e c t i o ni ss u g g e s t e d k e yw o r d s ：r e v e r s e t i m em i g r a t i o n ，a r b i t r a r yd i f f e r e n c e ，p r e c i s ei n t e g r a t i o n ， c o m p l e xs t r u c t u r e ，w a v ee q u a t i o n 关于学位论文的独创性声明 本人郑重声明：所呈交的论文是本人在指导教师指导下独立进行研究工作所取得 的成果，论文中有关资料和数据是实事求是的。尽我所知，除文中已经加以标注和致 谢外，本论文不包含其他人已经发表或撰写的研究成果，也不包含本人或他人为获得 中国石油大学( 华东) 或其它教育机构的学位或学历证书而使用过的材料。与我一同 工作的同志对研究所做的任何贡献均已在论文中作出了明确的说明。 若有不实之处，本人愿意承担相关法律责任。 学位论文作者签名：菱蜢 日期：z 口矿夕年多月圹日 学位论文使用授权书 本人完全同意中国石油大学( 华东) 有权使用本学位论文( 包括但不限于其 印刷版和电子版) ，使用方式包括但不限于：保留学位论文，按规定向国家有关 部门( 机构) 送交学位论文，以学术交流为目的赠送和交换学位论文，允许学位 论文被查阅、借阅和复印，将学位论文的全部或部分内容编入有关数据库进行检 索，采用影印、缩印或其他复制手段保存学位论文。 保密学位论文在解密后的使用授权同上。 翥等篡嘉篓霎芋一 指导教师签名：二三复之型型 日期：知哆年多月尹日 日期：勿，罗年石月卢日 l 中国石油大学( 华东) 工程硕士学位论文 第一章绪论弟一早珀t 匕 地震偏移技术，无论是过去，现在和将来都是地震勘探的重要内容【6 1 ，因为它是 与地震勘探的最终目的确定地下构造展布，地层接触关系，从而确定有利的油气聚 集带密切相关的。它是现代地震勘探数据处理的三大主要技术( 反褶积、叠加、偏移) 之一【3 】。现在偏移已经从研究简单的探测目标的几何形态进而发展成为研究反射界面 空间的波场特征，振幅和反射率等问题的有效方法。 然而随着勘探程度的进一步提高，地表、地下地质条件较好的地区，大部分已查 明，剩下的多是沙漠、戈壁、黄土塬、悬崖峭壁的山地、水陆交替的海滩、以及老油 田的边缘、深部、低渗透、小断块、陡倾角、特殊岩体等复杂构造。这些探区地表 地震条件恶劣，地下地质构造复杂，地震资料普遍存在信噪比低，成像困难的问题。 而目前我国地震勘探的主战场西部探区则大多属于这种高难地区【1 2 】。在这种情况下， 信噪比低，成像困难成为了首要问题。因为没有准确的成像，就不能正确认识地下真 实的地质构造。 偏移概念的发展经历了三个阶段，第一，把偏移过程看作是数据的空间映射；第 二，把它看作是非实时的合成聚焦成像；第三，偏移是一个逆传播过程。这三种观点 也代表了偏移技术发展的三个阶段。现在对偏移和偏移处理可定义为：地震记录上的 反射同相轴因为受波的传播特性的影响和记录方式的限制往往与其相应的反射地质 体在形态和位置上不一致，这种不一致性称为偏移，设法消除偏移影响的方法叫偏移 处理。在实际工作中，为了简单，常将偏移处理称为偏移【2 9 1 。 偏移技术由射线理论到波动理论代表着偏移技术的一大飞跃。目前常用的基于波 动方程的方法有k i r c h h o f f 积分方法，有限差分法和频率波数域方法。但各种方法都 有其适用范围和限制条件【l9 1 。如：k i r c h h o f f 积分能适应任意倾角，但不易确定偏移 孔径和格林函数，难以处理横向变速问题，且计算量大；有限差分偏移用近似波动方 程，有偏移倾角限制，而要进行大倾角偏移就要用高阶方程，势必会增大计算量。频 率波数域偏移用的是常速模型，当速度有变化是会出现明显的精度问题。早在二十 世纪9 0 年代初，马在田院士就指出了波动方程偏移的发展方r 甸t 3 1 ，即： l 、必须使用更精确的或很接近准确的波动方程。 2 、方法必须适应地下复杂的地质结构和速度模型。 3 、叠前偏移和深度偏移是研究重点。 第一章绪论 在前两种发展方向的指引下，本文提出了适合复杂构造成像的任意差分精细积分 逆时偏移方法。由于此方法使用全波动方程的同时，对偏微分方程的求解采用高精度 的任意差分精细积分，这势必会提高方法的精度。 实践证明逆时偏移该方法适用于二维和三维任意复杂的地质构造成像。由于没有 对方程的近似，避免了对方程的近似，同时也没有对速度的限制，因此可偏移任意倾 角的界面，适用于层间参数强烈间断等复杂情况，且由于逆时偏移采用全波动方程， 因此有较高的精度。同时这种方法能压制多次波和速度的扰动，成像精度较高【1 9 1 。 此外此方法既可以用于叠前偏移成像，也可用于叠后成像。将单炮集记录作为地 面边界条件，进行波场逆时延拓，采用激发时间成像的成像原理，就能实现叠前偏移。 将叠加剖面作为地面边界条件，采用爆炸反射面成像原理，就能实现叠后偏移。 鉴于其上述优势，对其进行深入研究有很大的现实意义。 本文将高精度精细积分应用于逆时偏移中。推导了任意差分精细积分逆时偏移的 公式，分析了差分格式的稳定性，从理论上证明了此方法是现有的精度较高的差分偏 移方法。各种复杂构造地质模型的成像表明了此方法的可行性，与其他偏移方法的对 比表明任意差分精细积分逆时偏移具有成像精度高，能灵活处理不规则网格划分，能 压制多次波和速度扰动，适合陡倾角地层构造，受噪音干扰小，成像精度高等优点。 由于避免了对方程的近似，同时也没有对速度的限制，而且最终得到的是偏移后的深 度剖面，因此该偏移方法属于深度偏移方法的范围。最后通过与现有的k i r e h h o f f 深 度偏移，相移法，显式有限差分深度偏移的对比，证明任意差分精细积分在速度模型 精确的情况下成像精度最高。但由于任意差分精细积分的高阶空间展开和逆时偏移的 计算使得此方法比较费时，并且逆时偏移的效果与所给定的速度模型的精确程度有 关，速度越接近真实速度成像精度越高。 2 中国石油大学( 华东) 工程硕士学位论文 第二章有限差分正演模拟 在介绍任意差分精细积分逆时偏移时，首先讨论有限差分正演问题。大家知道地 震波数值模拟是波动方程的正演问题，是在地下介质结构和参数已知的情况下，利用 数值计算的方法来研究地震波在地下介质中的传播规律，从而获得理论上地震记录的 一种技术。数值模拟技术是勘探地球物理学一个重要组成部分，其作用主要有两点【2 0 】： 一是作为一种认知工具，帮助我们提高对复杂介质中地震波传播规律的认知水平，从 而指导地震野外观测系统方案的设计以及对野外采集数据进行实时质量控制，帮助进 行地球物理资料解释与地质解释；二是作为一种检验工具，检验各种方法技术的应用 效果及适用范围。 地震波场数值模拟的主要方法包括两大类【2 1 1 ，即波动方程法( 克希霍夫积分法、有 限差分法、有限单元法、边界单元法和虚谱法) 和几何射线法( 射线追踪法) 。波动方程 数值模拟方法实质上是求解地震波波动方程，因此模拟的地震波场包含了地震波传播 的所有信息，但其计算速度相对于几何射线法要慢。几何射线法也就是射线追踪法， 属于几何地震学方法，由于它将地震波波动理论简化为射线理论，主要考虑的是地震 波传播的运动学特征，缺少地震波的动力学信息，因此该方法计算速度快。因为波动 方程模拟包含了丰富的波动信息，为研究地震波的传播机理和复杂地层的解释提供了 更多的佐证，所以波动方程数值模拟方法一直在地震模拟中占有重要地位。 本章首先介绍了有限差分法正演模拟的概念和任意差分精细积分数值模拟的实 现，继而介绍了数值模拟中经常用到的爆炸反射面模拟的概念，最后详细阐述了差分 法数值模拟中需注意的差分格式稳定性，数值频散和边界的处理。 2 1 有限差分法数值模拟 有限差分法是一种最常用的正演模拟方法。1 9 7 0 年，a l t e r m a n 等首先将有限差 分法应用于地震波动方程模拟中。之后，a l f o r d 等研究了有限差分法的精度，v i r i e u x 提出了稳定的二阶( 空间和时间) 弹性波有限差分格式，它适用于任何泊松比的介质， l e v a n d e r 将v i r i e u x 的方法推广到空间四阶、时间二阶的情况。c r a s e 则发展了精度可 达任意阶的高阶交错网格法，但其计算量和内存要求比低阶有限差分法大幅度增加。 m a g n i e r 等提出了最小网格有限差分法，它能压制非最小网格的人为现象。周家纪和 贺振华用大网格快速差分算法模拟地震波的传播，其空间网格可以取得很大，达到每 最短波长只需3 个网格点，大大缩短了计算时间；刘洋在“任意偶数阶精度有限差分 第二章有限差分数值模拟 法数值模拟”文苹中将差分方法的空间精度推广至任意偶数阶，但时间精度未见提高 【8 】o 虽然当前的地震波数值模拟有许多方法，但在振幅保真、数值频散、模拟精度、 边界吸收、计算速度方面还有待提高【6 】。 以二维连续介质为例，在二维连续介质中，变速声波方程可描述为 a d 2 f p 。_ c 2 ( x , z ) 0 d x 2 p _ - 0 a z 2 p ： i 1 - 厂( 2 1 ) 式中：p = p ( x ，z ，t ) ，为质点运动的位移位；f = - f f x ，z ，t ) ，为震源分布函数：e ( x ，z ) 为c r 质的速度分布。对上式离散化，时间导数项采用二阶中心差分来逼近，空间导数项则 根据t a y l o r 级数展开式来构造，得到高阶差分逼近式，而且我们知道差分的阶数愈高， 精度愈高。 一个函数p ( x ) 的二阶空间导数等0 2 的高阶差分可表示为l 】 p ( 斛池m ( 卜池) - 2 p - 2 喜鲁) 上式中：a x 为网格步长；i 为偏离x 位置的网格点数；p ( 2 n 为p ( x ) 的2 n 阶导数。 馨的2 m 阶差分逼近式可写作： o x 。 窆o x i x f f i x o 黑学删矿， 2 乌 血2 p 。， 其中： c 腑：上掣f i ，z f 矗( 磊：一一卉( ，z 一朋z ) m i 主鑫m s 鑫嘉1) 所= 1 ，2 ，m 据此得到式( 2 - 1 ) 的2 m 阶显式有限差分时间递推格式 p 舢_ 2 旷p 似+ 警，新p ( 矿池帅“酬像2 ) 这是一个显格式的时间递推差分格式，其中c - m = c m ，则由t - a t 和t 时刻的波场值外 推可得到t + a t 时刻的波场值。这是传统的有限差分波动方程解法。在此方法的基础 上发展了高精度正演模拟的任意差分精细积分方法。 2 2 任意差分精细积分正演的原理 与有限差分相比，任意差分精细积分以其空间离散的任意性和时间积分的精细性 4 中国石油大学( 华东) 工程硕士学位论文 在提高正演模拟精度方面具有很大优势。 这里仍旧以二维连续介质声波方程为例： 吉窘= 窘+ 害- 万c x x o ，z - - z 0 - i - 小，( 2 引 一一= 一十一， i7i f -i ，一， 1 ，2a 2 苏2 瑟2 。、 7 八7。 其中：8 ( x x 0z z o ) 厂( f ) 为震源项。 首先，对式( 2 - 3 ) 进行二维空间坐标离散，即在任一离散点j ( x ，z ) 处，有 1 可a 2 u j = 碧+ c 0 2 u ：jb拂删xj-xo,zj-zoaz z j ) ，2西2苏2 2 。( 。，：) = ( 、7 八7 其次，将u 在j 点的x 方向邻域做t a y l o r 展开， = ，+ 善去( 她昙广”k _ + o 【( c 卅】 式中，缸，= x 广x ，i 是j 邻域内的任意点，这正是a d p i 法中“任意差分”的意 义，这使得a d p i 法可以方便地处理非均匀介质与边界条件( 并且可以保证差分近似 的精度) ，以及实现局部点的加密【4 吲。对上式移项并加权得 善n ( u i - - u j ) = 喜( 善击( 血，昙广乩：勺+ 。【( 缸厂l 】) 其中：口；是加权系数，m 原则上可任意选取，但m n 使下式的方程组成为矛盾 方程组，m n 时下式的方程组为欠定的，因此一般取m - - - n ( 以m = n 为例来展开讨论) ， 则由上式展开，有 这是一个关于口，的m 维线性方程组，解之可求出加权系数。同理，可将u 在j 点的z 方向邻域做t a y l o r 展开( 以均匀网格划分为例) 。最后，声波方程右端的二阶导 数项可用相邻点的任意差分来表示( 暂不考虑震源) ，那么式( 2 1 ) 可表示为( 暂不考虑 震源项) 2 = 2 、_ 、 x ，- 一 口 瑚 0 = 4 、 x ，- ， 口 所瑚 0 = 3 、， x _ ， 口 册m 0 1 1 圻 、- 、 xa l 伢 卅渊 o i i o所 、， x ，l 口 所 第二章有限差分数值模拟 专等=塞咖，1)ot1 ，2 2 鲁v “ 在一般情况下( 有限差分) ，该式的右端可近似为f 。时刻的值，对左端作二阶时 间差分近似，得到与( 2 - 2 ) 式类似的递推关系式，但在任意差分精细积分法中，认 为右端的铭，没有必要做近似，因此上式可写为： 专争f u j 塞：艺q 比 丁蒂缶2 备q 比 若右端采用时间的零阶近似，则此时的半解析解( 此即精细积分的意义) 为 “，= 叩x p ( i a v ( f - f 。) ) + c 2e x p v ( h 棚+ 争 上式中，口2 = ，b n = 口，甜矿 在i t i ，“l 】区间求解，消去c ic ：可以得到如下的递推公式： 肿l - 2 。c o s 心) 叫肋一l + 2 1 - c o s 心) 】( 2 - 4 ) 若将右端项展开到，一t 。的二阶，即： 坼= “枷+ “枷。一乙) + 去“枷。一乙) 2 这时的半解析解为： 吩邓x p g 川一m c ：e x p ( - i m ( t 训) 寺+ 掣+ 尘垡二刍! 一堕 口2d 4 v 2 2 m2 * m2 * m 12 * m 其中：口2 = ，吃= 甜枷，气= l j i , n 吒= 去“枷 1 = 1_ li = l- ，= l 同理，选择积分子域为 乙，乙+ 。】，则利用t n - ! , r 。消去c 。，c ：，进而可以求出，州时 刻的值： 咋。槲= 弛。s ( a v e ) 一吩川+ 2 夸一番) 一c 。s 心) 】+ 等( 的2 ( 2 - 5 ) 在上式中加入震源项，则可以据此公式进行地震波的正演传播模拟。由于此方法 在对空间高阶离散的同时对时间采用精细积分，求出解析解，从而得到递推的差分格 式，因此其精度与传统的有限差分相比精度有所提高【8 1 。 6 中国石油大学( 华东) 工程硕士学位论文 2 3 爆炸反射界面模拟 我们知道水平叠加剖面可以看作是沿测线单点自激自收的观测剖面。在水平均匀 介质情况下，水平叠加剖面和沿测线单点自激自收的观测剖面完全等价，在非均匀介 质非水平层界面的情况下水平叠加剖面和沿测线单点自激自收的观测剖面只是近似 相等。由于经济和效率的原因，地震勘探从来也没有使用过自激自收观测方式。所以 常常把水平叠加地震剖面看作自激自收的地震剖面【3 】。 自激自收剖面的特点就是各个点上的入射波和反射波沿着同一条法线射线传播 的【l 】。因此在地面上观测到的反射波可以认为是反射界面上同时激发地震波在地面上 观测的结果。这就是l o e w e n t h a l e t a l 等人首先提出来的爆炸反射面的原理。该原理把 地下反射界面想象成具有爆炸性的物质或爆炸源，爆炸源的形状，位置，与反射界面 的形状和位置一致。它所产生的波为脉冲波，其强度，极性与界面反射系数的大小和 正负一致。并且假定在t = o 时刻所有爆炸反射界面同时起爆，发射上行波到达地面各 观测点，波的传播速度为真实地层速度的一半【2 9 1 。如图2 1 a 、2 1 b 所示。 0 z 图2 1 a x 0 z 图2 1 b x 图2 1 a 自激自收反射界面模型，上行波与下行波路径相同，波的传播速度为v ， 图2 1 b 爆炸反射界面模型，t = o 时，波自界面发出( 只有上行波) ，波的传播速度为 v e = v 2 ，由于v e 比v 小一倍，因此模型a 和b 中波到达地面的时间相同。 由于爆炸界面法正演模拟可以减少计算量，并且对于叠后偏移方法的检验与单炮 资料处理后得到的叠加剖面是一致的。因此这里采用的是由爆炸反射模拟的地震记录 剖面来代替叠j 3 n 静 面。所谓爆炸界面法是指将激发源置于地下各个层面，通过接收自 下而上传来的波来求取地震记录。值得注意的是在计算时要将真实介质的速度除以 2 ，即所用于正演模拟的声波方程为 7 第二章有限差分数值模拟 粤o t = 地4 i 磐0 x + 磐0 zl + 厂 2 i 2 。2 i 2 4 有限差分正演模拟中稳定性，频散和边界的处理 在用有限差分法对偏微分方程做离散近似的时候，不可避免的会带来数值计算稳 定性和频散问题，同时由于正演模拟的空间的有限性，对边界的处理也是一个不容忽 视的问题。 2 4 1 稳定性研究 对于二维弹性介质中声波波动方程的2 n 阶精度差分方程【1 1 】 p 川n + l 却”，n - i 蛋( p n j - m , i 4 - p n j + m j ) + c o 办 ，2 f 2缸2 c ，( p “川一m + p j ，f + m ) + c 。p ” 型三! a z 2 令缸= z = d ，a 2 = v 2 a t 2 a d 2 ， 则上式可化为 + j p 肘1 = 一g ”_ ，+ a 2 c 。( p ”一。，+ p ”+ m ，+ p ”- ，一脚+ p ”，+ m ) + 2 c o p ”+ 2 ig 肿1 = p ”j , t 令u nj 。| = l p n q nj i 、) t = wn e 谩l j 厶de 如| d 则上式可简化为 缈肘1 = g 矿”= 言二1 ” 其中g 为过渡矩阵，g 的表达式为： g = 2 a 2 c 。 c o s ( m k ，a d ) + c o s ( m k ：a d ) + 2 a 2 c o + 2 m = l 我们知道当过渡矩阵所有特征值的模小于1 时，则差分格式稳定，可以证明过渡 矩阵所有特征值的模小于1 的充要条件为l g l 2 。 考虑最大空间频率( 即尼奎斯特频率) ，则k x = 屯= x a d ， 将此式代入上面的表达式，并注意到c o = - 2 c 。，因此上式可化简为 8 中国石油大学( 华东) 工程硕士学位论文 p k l 2 m = l f 一8 口2 - l + 2 l ll 其中n i 为不超过n 的最大奇数，进一步化简，可得稳定性条件 面v a t 一瓯。 对比其与刘志敏、包吉山在高阶有限差分格式频散中得出的稳定性条件可知其结 论是一致的f 2 3 1 。 为了与有限差分法对比，在任意差分精细积分方法的稳定性研究中，我们假设 a x l = a x 2 = a d 。 对于任意差分精细积分零阶时间展开的差分格式： + l - 2 。c o s ( 舭沪l + 2 争 1 _ c o s ( 砒f ) 】 采用v o nn e u m a n n 方法来分析其稳定性，注意到齐次线性差分方程的误差传播 方程与原方程的形式一样，及j 的邻点取在与j 的距离相等的两侧，我们将误差的付 氏级数的一个分量表示为： 面加= 圪e x p ( i k x ) 面，卜l = 圪一le x p ( i o g ) 8 u 伽+ l = le x p ( i k r ) 面，i 。= 圪e x p ( i k ( x 一缸) ) 面+ 1 。= 圪e x p ( i k ( x + a x ) ) 式中：f = 一l ，得到误差传播方程为： 8 u 肼l _ 2 t 勋j nc o s 。v a t ) 一( 蜀u j n _ i j r2 斋 1 - c o s ( a v a t ) 将其代入上式，并消去公因子得到： 圪+ l = 2 圪e o s ( a v a t ) - v 一l + 2 圪c o s ( c a x ) ( 1 一e o s ( a v a t ) ) 我们设误差传播因子为g ，则v n + l - g 母v n ，代入上式并消去公因子得到关于g 的 方程： g 2 2 ( c o s 辈) + c o s ( ) - - c o s f a v a t ) c o s ( k a x ) ) g + 1 ：0 xa x 上式中，g 的模小于等于1 的充要条件是：c o s ( 竺垒! ) 0 ，也就是： x 9 第二章有限差分数值模拟 2 n t r - t r 2 _ a w a t 2 n z r + ，r 2 ，n ：o ，l ，2 ，。 出 实际算例说明：只有当n = 0 时的区间是有效的，其他区间的解虽然稳定，但误差 太大，已不可用【引。因此对于二阶差分实际的稳定性条件为： 一v a ts 万( 2 压) ：1 11 0 8 9 a x 对比有限差分法，_ v a t l ，可知a d p i 方法吏为稳定，对比八阶有限差分格式的 工 稳定性条件_ v a t 0 5 5 4 6 ，任意差分精细积分的稳定性条件为堂0 6 5 8 3 ，可知任 xx 意差分精细积分格式更加稳定。 2 4 2 频散的控制 限差分法是声波或弹性波场数值模拟中最为流行的方法之一，在求解波动方程 时，会产生不期望的数值频散或称网格频散，导致了数值模拟结果分辨率的降低。究 其根源，是因为基于波动方程的有限差分求解过程，通常是利用离散化的有限差分方 程去逼近波动方程，从而使得相速度变成了离散空间间隔的函数。所谓数值频散实质 上是一种因离散化求解波动方程而产生的伪波动，这种频散既不同于波动方程本身引 起的频散，也不同于因波传播的速度、频率和角度而引起的频散【2 1 1 ，它是有限差分方 法求解波动方程时所固有的本质特征，无法避免，但可以通过一定的参数的选取来控 制使其得到削弱。 我们知道，当差分格式精度较低时，频散现象严重，随着阶数的增加，频散现象 逐渐降低，模拟波场的精度逐渐提高【l l 】。图1 2 a 和1 2 b 给出了分别为二阶和四阶差 分的对比结果。 显然，随着差分精度的增加，参与计算差分的点数也会增加，计算量随之增大。 所以，要获得较高的精度，必须花费较多的计算时间。要想快速得到正演模拟结果， 就要以降低精度为代价。总之，根据所要求的计算速度和精度，找到一个阶数最佳的 差分格式来进行计算。 为了消除这种数值频散，人们进行了大量的研，他们的结论基本一致，即为了消 除数值频散，在使用二阶有限差分方法时，每个1 2 功率对应的波长至少必须使用1 1 个网格点【2 3 】，a l f o r d 的结论为不少于l o 个【2 2 】，而四阶有限差分则可用二阶差分网格 点数的一半。 1 0 中国石油上学( 牛东) 工程硕士学位论立 ；， 图2 - 2 a2 阶差分的时闯快照图2 2 b4 阶差分的时间快照 值得注意的是通常人们认为有跟差分的精度与差分阶数和空间网格步长有关。所 以一般做法是在介质速度、差分阶数、子波频率和空问时恪步长选定的情况下，根据 有限差分的稳定和收敛条件，选取时问递推步长。在刘志敏、包吉山高阶有限差分格 式频敞分析的基础上分析可知时间递推步长与频散也存在很大关系1 2 ”，选用优化时间 递推步长可以在很大程度上控制频散。 确定优化时间递推步长的具体做法是：通过分析波的传播方程，把平面波形式解 “= a e ( 。0 5 “+ 2s i n “一删 代入式( 2 2 1 所示的2 m 阶显示有限差分格式， ! 型：。二! ! ：! ：! ! 型：三：! ! ! ：! ：! ! ：! ：a ! ：+ 一出2 a x 2 c ，( p ”j m + p ”，m ) + c o p “ ! ! i a z 。 可得出有限差分数值相速度与介质速度的关系，即 生：三。i 。塑 ci 抽、2 其中：h = a x ，= c a r a x ，k = 2 n ，c p 为数值相速度：c 为给定介质速度。c 以越接 近于1 ，差分格式的频散越小，差分精度就越高。因此一旦给定介质参数与子波频率， 待确定的参数只有空间网格步长x 及时f h j 递推步长t 。依据a l f o r d 等的研究成果 根据子波频率和介质速度来设置空b j 采样率，保证每个二分之一功率对应波长的采样 率足够高就能有效抑制频散，但对于时间递推步长t ，由上式不能得出显式求解表 第章有限差分数值模拟 达成，但可用迭代算法求取t 。 速度c 取介质的平均速度值，给定空间网格步长x ，未知量只有t ，反复修f l 值使c 以尽可能地接近于1 ，最终得到适宣的时间递推步长。这里给出了接近最佳 步长与任意时间步长数值模拟效果对比图( 图2 - 3 a ，2 - 3 b ) 。该理论模型为一均匀介质， 速度为2 0 0 0 m s ，点震源位于模型中央。j 下演计算的子波为雷克子波，其主频为3 5 h z ， 恻格步蚝b x = a z - 1 0 m 。图2 3 a 和图2 3 b 分别是采用任意时间步长l m s 和最佳时间 步长2 m s 存某时刻的波场快照图。 一 图2 - 3 a4 阶差分l m s 采样的时间快照圈2 - 3 b4 阶差分2 m s 采样的时间快照 通过对比可以看出，采用最佳的时间步长2 m s 得到的时间切片( 图2 - 3 b ) 频散小， 波形稳定，能量集中。而i m s 时间步长对应的时间切片( 图2 - 3 a ) 则达不到如此效 果。崮此采用虽佳的时间步长可以对差分频散起到一定的压制作用。 对于任意差分精细积分方法无法得出显式的差分参数( 空间步长，时间步长) 的 表达式但是从其差分原理和实际计算表明当差分格式精度越高，频散越小。唰时当 相同阶数的有限差分格式满足要求时a d p i 方法也能达到很好的频散控制效果。 2 43 边界条件 在有限差分法波动方樗j 下演及偏移成像过程中，对正演过程而占，模型计算的区 域不可能为无限。对偏移成像而盲接收信号的检波器分布范围也是有限的。这就给 上述两种计算引人了人为的边界。这些人为的边界会反射由计算区域中传来的波，使 得正演结果及偏移成像结果受到干扰。为了消除这些边界反射波，人们提出了各种方 法，比较常用的有以下几种：一是扩展边界法，即在给定的计算区域外，扩展若干道， 中国石油大学( 华东) - 1 - 程硕士学位论文 使人为边界远离给定的计算区域，这样边界反射对计算区域的干扰自然会小一些；二 是衰减边界法，即在计算区域的边界部分设置若干衰减道，使得向计算区域中传播的 边界反射波得到衰减；三是吸收边界法，即导出吸收边界条件方程，使之与计算区域 内的波动方程结合在一起求解，使得从计算区域传到边界上的波被全部或部分吸收掉 【2 6 】 o c l a y t o n 和e n g q u i s t 导出吸收边界条件方程，是通过把波动方程按幂级数展开分 解成各阶左行波和右行波方程。其吸收机理是：在左边界上吸收右行波，即相当于左 行波到达左边界后，大部分能量继续向前传送，反射回来的能量只占很小一部分；在 右边界上吸收左行波，即右行波到达右边界后，同样是大部分能量仍然一直向前传送， 很少能量被反射回来。这自然便可以达到吸收边界反射波的目的。c l a y t o n ( 1 9 8 0 ) 提出 如下三种形式的右边界条件： 零阶边界条件 一阶边界条件 二阶边界条件 堕：6 ( 1 9 诹：一 。k ，t - 二= 一a + d 2 国缈 堕：一( d b v k r ) ( 1 一c 堕) 上述二柙彤瓦边界条仟网堤出，均妥涿征瓦甲小h 岜出士见k 一网局于一次明坝，如 果出现了这些高次项后，在边界条件微分方程中就会出现波场关于x 的二阶以上的微 分项，在进行差分计算时便会涉及到边界外的值。因此在计算时采用的是一阶近似的 边界条件。将其通过傅立叶反变换变换到时空域，得到如下的边界方程： 娶一! 娶+ v 擎：0左边界 a x 8 tva 2 ta 2z 7 昙善+ ! 宴一v 娶：0右边界 a x a t va 2 t a 2z 喜一三擎+ v 拿：0下边界 a z a tv8 2 t a 2 x 娶+ 三擎一v 宰：0上边界 8 z a tva z t8 zx 。j 近年来b e r e n g e r 针对电磁波传播情况，给出了一种高效的完全匹配层( p m l ) 嚼婀沩界条件算存殚诊卜证明该方法可以完伞吸收来自各个方向、各种频率的电磁 第二章有限差分数值模拟 波，而不产生任何反射1 2 6 j 。 对于二维声波方程 - a 2 p = v 2 0 , 2 + 0 如2 _ 卫p 明 引入中间变量u l ，u 2 ，a l ，a 2 则此方程可以写作： 甜= “l + 甜2 鲁= 叭扎z ) 堕o x 等：y ：( 即) 警 at 、j。 az 堕：堕+ 堕 0t0 x8x 至生：堕+ 垒生 0t0zaz 其相应的完全匹配层控制方程为： 等“( 咖。= 暇础百a a l 鲁峨驴暇础百a a 2 鲁“c 圳。= 等十等。 鲁“小等+ 百o n 2 上式中d l ( x ) 和d 2 ( z ) 分别为x 方向和z 方向的衰减系数，也就是说d l ( x ) 起到衰减 x 方向传播的波、d 2 ( z ) 起到衰减z 方向传播的波的作用。对于任意方向传播的波：可 以通过矢量分解，分解成x 方向和z 方向传播的波，分别进行衰减。波场是按传播距 离的指数规律衰减，衰减速度很快。当衰减系数d l ( x ) 、d 2 ( z ) 随空间位置变化时，不 会在介质中产生任何反射。这些性质使得这种介质特别适合用作声波方程的边界吸收 介质。 图2 - 4 p m l 边界划分 1 4 中国石油大学( 华东) 工程硕士学位论文 利用完全匹配层作为吸收边界的基本做法是在所研究区域的四周引入完全匹配 层。如图2 - 4 所示，区域a b c d 为所要研究的区域，即我们要在此区域中研究波的传 播问题。在区域的周围加上完全匹配层，在区域1 ( 四个角区域) 中，令d l ( x ) 0 ，d 2 ( z ) 0 ， 速度v 都等于角点的速度。在区域2 ( 上下区域) 中，令d l ( x ) = 0 ，d 2 ( z ) o ，速度v 在 z 方向为常数，在x 方向和边界的速度相等。在区域3 ( 左右区域) 中，令d l ( x ) o ； d 2 ( z ) = 0 ，速度v 在x 方向为常数，在z 方向和边界的速度相等。这样在计算边界的 周围都有完全匹配层吸收介质，波由区域内通过边界传播到完全匹配层时，不会产生 任何反射。波在完全匹配层中传播时，不会产生反射并且按传播距离的指数规律衰减。 当波传播到完全匹配层的边界时，波场近似为零也不会产生反射。数值模拟( 平界面 单炮记录，图2 5 ) 的结果表明， 2 0 0 2 0 0 2 0 0 0 m s 3 0 0 0 m s 图2 - 5 平界面模型 用最佳匹配层吸收边界合成的地震记录，在振幅增强的显示下，也看不到明显的 边界反射，而c l a y t o n 吸收边界就不能达到相同的效果，因此最佳匹配层吸收边界比 c l a y t o n 吸收边界在压制边界反射方面效果要好。如图2 6 a 和2 6 b 所示。基于上述原 因，在做波场的数值模拟时使用的是p m l 边界。 第：章冉限差分数值模拟 图2 - 6 ac l a y t o n 一阶近似吸收边界效果 图2 6 bp m l 边界效果 中国石油大学( 华东) t 程硕士学位论文 第三章任意差分精细积分逆时偏移 逆时偏移方法最早是1 9 8 2 年w h i t e m o r e 在美国的d a l l a s 召开的第5 2 届s e g 年 会上提出的，后经过了多位学者的发展和完善。1 9 8 3 年b a y s a l ，e ，d d k o s l o f f , a n dj w c s h e r w o o d 提出了不同的逆时偏移的概念【1 7 1 ，同年，l o e w e n t h a ld a n di r m u m 将其应用在频率一空间域偏移【1 8 】。1 9 8 4 年，l e v i ns a 概括了逆时偏移的基本原理 和实现方法【16 1 。1 9 8 7 年h i l d e b r a n ds t 将其应用于波阻抗成像，取得了很好的效果。 1 9 8 8 年e s m e r s o yc a n dm o r i s t a g l i o 研究了逆时波场的外推，成像和反演【1 4 】。同年 z h uj a n dl l i n e s 比较了逆时偏移与克希霍夫积分，得出前者对m a r m o u s i 模型成像 精度更高，但费时也多【1 2 】。s u nr a n dga m c m e c h a n 于2 0 0 1 进行了标量波动方程 的逆时深度偏移，对纵波和横波成像表明比单纯声波效果要好。逆时偏移方法适用于 二维和三维任意复杂的地质构造，由于没有对方程的近似，避免了对方程的近似，同 时也没有对速度的限制，因此可偏移任意倾角的界面，适用于层间参数强烈间断的情 况，且由于逆时偏移采用全波动方程，因此有较高的精度。同时这种方法能够压制多 次波和速度的扰动，成像精度较耐1 9 1 。 本章首先介绍了传统的有限差分深度方向延拓偏移方法，继而对比介绍了单程和 双程波动方程逆时偏移的原理，着重讲述了任意差分精细积分全声波方程逆时深度偏 移的应用，最后通过一系列数值例子的检验验证了任意差分精细积分逆时偏移，针对 复杂构造成像的效果。 3 1 逆时偏移原理 本节首先介绍单程波方程向下延拓有限差分法波动方程偏移的实现过程，继而对 比介绍单程和双程波动方程逆时偏移的原理。 3 1 1 向下延拓有限差分法波动方程偏移 延拓又