（应用数学专业论文）地震数据曲网格模拟与叠后深度偏移.pdf

摘要 在直角坐标系下用传统伪谱法模拟地下界面时，弯曲起伏的连续界面 会因离散化而成为阶梯状界面，这种阶梯状界面在模拟波的传播时，将产 生人为绕射，降低合成记录的精度。本论文实现了一种精确有效的方法来 模拟曲界面波场。 该方法结合了多块映射和超限插值技术，这项技术出自计算流体力学， 它可使产生的网格网线沿着地下所有界面变化。这一映射技术将直角坐标 系下的曲界面变换成曲坐标系下的规则界面，在曲坐标系下用伪谱法模拟 波场，最后利用插值得到实际波场的合成记录。这种方法能有效地消除在 直角坐标系下因离散化曲界面产生的不必要的人为绕射。曲网格伪谱法用 较少的网格数可以获得与直网格伪谱法相同的精度，从而节省了计算机内 存，这一点对于三维地震模型的计算具有重要意义。通过声波模型和弹性 波模型在直网格下的合成记录与在曲网格下的合成记录的比较验证了方法 的有效性与精确性。 傅立叶有限差分( 简称f f d ) 偏移算法是一种叠后深度偏移算法，其向 下延拓算子是一种混合算子，包括三项：一项是处理常速的相移算子，一 项是一阶相移修正算子，最后一项是类似4 5 度方程的有限差分算子，用来 处理剧烈横向变速。f f d 算法既保持了相移法的精确性，又可以处理剧烈 横向变速问题，在一定程度上克n t n 移法偏移和有限差分法偏移的缺点， 是一种比较有效的偏移方法。本论文用连分式近似单程波波动方程中的平 方根导出f f d 算法的基本公式，并对f f d 算法中的有限差分算子进行了 系数优化，进一步提高了计算的有效性。 关键词：地震模型；深度偏移；超限插值；伪谱法；合成记录：波场快照 波场延拓 a b s t r a c t w h e n a p p l y i n gt h ec o n v e n t i o n a lf o u r i e rp s e u d o s p e c t r a l m e t h o do nc a r t e s i a n g r i d s ，c u r v e di n t e r f a c e sa r er e p r e s e n t e d i na s t a t r c a s ef a s h i o n c a u s i n gs p u r i o u s d i f f r a c t i o n s i ti s d e m o n s t r a t e dt h a tt h e s e n o n p h y s i c a l d i f f r a c t i o n sc a nb e e l i m i n a t e db yu s i n gc u r v e dg r i d st h a tg e n e r a l l yf o l l o wa l lc u r v e di n t e r f a c e s t h ec u r v e dg r i d sa r eg e n e r a t e db yu s i n gt h et r a n s f i n i t ei n t e r p o l a t i o na n dt h e s o c a l l e dm u l t i b l o c kt e c h n i q u e st h a to r i g i n a l l yd e v e l o p e df o rc o m p u t a t i o nf l u i d d y n a m i c s t h ec u r v e dg r i d i st a k e nt oc o n s t i t u t eag e n e r a l i z e dc u r v i l i n e a r s y s t e m t h u s ，t h ew a v ee q u a t i o n sh a v e t ob ew r i t t e ni nc u r v i l i n e a rc o o r d i n a t e s b e f o r ea p p l y i n gt h en u m e f i c a is c h e m e b e c a u s et h eg r i di sc a r t e s i a ni nt h e c u r v i l i n e a rd o m a i n ，s t a n d a r dp s e u d o s p e c t r a lt e c h n i q u ec a l lb ea p p l i e d a t1 a s t t h es y n t h e t i cr e c o r do fp h y s i c a lc o n f i g u r a t i o nc a nb ee v a l u a t e d c o m p a r e dt o s i m i l a rm e t h o d so nc a r t e s i a ng r i d s ，t h es a m ea c c u r a c yi so b t a i n e dw i t hal o t f e w e tg r i dp o i n t s ，w h i c hm e a n st h a tc o n s i d e r a b l es a v i n g si nc o m p u t e rm e m o r y c a l lb eo b t a i n e d t h i sf a c th a sa ni m p o r t a n ti m p l i c a t i o nf o re x t e n s i o nt o3 d c o n f i g u r a t i o n s c o m p a r i n gt h es y n t h e t i cs e i s m o g r a mb yu s i n gc a r t e s i a ng r i d s w i t ht h es y n t h e t i cs e i s m o g r a mb yu s i n gt h ec u r v e dg r i d i tc a np r o v et h e 一一一 a p p r o a c h i se f f e c t i v e a h y b r i dm i g r a t i o nm e t h o d ，n a m e d “f o u r i e rf i n i t e d i f f e r e n c em i g r a t i o n ” i sap o s t s t a c kd e p t hm i g r a t i o ns c h e m e t h ed o w n w a r de x t r a p o l a t i o no p e r a t o ri s s p l i ti n t o t h r e eo p e r a t o r s ：o n eo p e r a t o ri sap h a s e s h i f to p e r a t o rf o rac h o s e n c o n s t a n tb a c k g r o u n dv e l o c i t y , a n o t h e ro p e r a t o ri sd l ew e l l k n o w nf i r s t o r d e r c o r r e c t i o nt e r m ，a n dt h et h i r do p e r a t o ri saf i n i t e d i f i e r e n c eo p e r a t o rf o rt h e v a r y i n g o ft h e v e l o c i t y f u n c t i o n p h a s e s h i f td o w n w a r de x t r a p o l a t i o na n d f i n i t e d i f i e r e n c ed o w n w a r d e x t r a p o l a t i o np r e s e r v e st h ea d v a n t a g eo fp h a s e ，s h i f t m e t h o da n df i n i t e d i f f e r e n c em e t h o d s of f d m i g r a t i o nm e t h o di sa ne f f e c t i v e s c h e m e i nt h et h e s i s ，t h ef u n d a m e n t a lf o r m u l ao ff f dm e t h o dd e r i v e sf r o mt h e s q u a r er o o t t h a ti s a p p r o x i m a t e db yac o n t i n u e df r a c t i o ne x p a n s i o ni nt h e o n e w a y w a v e e q u a t i o n o p t i m i z a t i o n s o ft h e p a r a m e t e r s o ft h e f i n i t e d i f f e r e n c eo p e r a t o ri m p r o v et h ev a l i d i t yo f t h em e t h o d k e y w o r d s ：s e i s m i cm o d e l i n g ；d e p t h m i g r a t i o n ：t r a n s f i n i t ei n t e r p o l a t i o n ； p s e u d o s p e e t r a l ；s y n t h e t i cr e c o r d ；s n a p s h o t ；w a v e f i e l d e x t r a p o l a t i o n 查堡堡! - 大兰塑丝壅! ! 童 前言 地震勘探是探测石油与天然气的重要方法。世界上绝大多数油气田主 要是通过地震勘探方法找到的。利用地震勘探方法并结合其它地球物理资 料就可以确定地下岩层的几何形态和物理性质，并进一步判断地下含油气 构造的存在，为下一步的钻井提供恰当的井位。 地震勘探方法包括野外数字地震信号采集，地震资料( 即信号) 处理 和地震资料解释三大步骤。地震信号采集采用人工方法激发地震波，在地 表按一定规则布置检波器，接收人工激发的地震波在地下传播时遇到地下 界面反射回的地震波构成所旧地震资料。再使用各种数值技术在计算机 上对这些地震资料进行处理，以得到地下界面的形态( 即地震剖面) 及地 层的物性参数。地震资料解释的任务是通过划地震资料( 即信号) 处理得 到的地震剖面，综合各种地质资料，进行分析解释，确定油气藏的空间位 置，以最后确定井位。本文仅涉及到地震资料处理的一部分：偏移成像和 模型计算。 根据地面接收的地震资料用波动方程进行偏移处理的目的是：( 1 ) 通 过从地面得到的地震资料确定反射点( 或绕射点) 的空间位置，也就是使 地下界面反射轴归位( 偏移) ，绕射收敛( 成像) 得到偏移剖面。( 2 ) 恢复 反射波在地下空间位置上的波形和振幅特性，为岩性分析提供依据。通过 对偏移剖面进行解释，可得到有用的地质构造信息，并最终用于判断油气 藏的存在。地震数据偏移是地震资料处理的核心内容和必不可少的环节， 它与地震勘探的最终目的相连，并直接影响着地震勘探的最终效果。近三 十多年来，由于经济发展的需要和计算机技术的飞速发展，地震数据偏移 技术，特别是基于波动理论的偏移技术从无到有，得到了迅速发展，不仅 形成了一门崭新的学科，而且得到了广泛的工业应用，获得了巨大的经济 效益。在类型上，偏移算法可从二维或三维、叠后或叠前、时间或深度来 分类，如二维叠后时间偏移、二维叠后深度偏移、二维叠前时间偏移、：二 维叠前深度偏移、及三维叠后时间偏移等八种不同类型。每种类型都反映 了偏移算法发展的不同阶段，计算的规模，难易程度和成像的要求。本文 仅涉及到二维叠后深度偏移。 叠后深度偏移有两个基本要素：( 1 ) 波场延拓或波场外推( 2 ) 波场成 像。为解决波场成像问题，可假设把炮点分布在所有反射界面上，并在零 大连理下丈举硕士论文 前言 时刻同时激发，激发产生的子波振幅同该点的反射系数成正比，并令介质 的速度为真实速度的一半，则可获得与自激自收等价的叠加剖面。这就是 d ，l o w e n t h a l 提出的爆炸反射面模型峨根据这一物理模型，并结合上行波 方程，地球物理的偏移问题就化为求解如下的数学问题：己知叠加剖面p ( x ， z = 0 t 1 ，通过上行波方程求偏移剖面p ( x ，z ，t = o ) 。爆炸反射层模型解决了 波场成像问题。不用的波场延拓方法形成了不同的偏移算法。地震偏移技 术中最基本的算法有：缈一耳域的相位移算法哆x t 域的有限差分算法”， f k 域的s t o l t 算法1 4j ，积分法”j 。 有限差分法最大的优点是介质速度可以纵横向变化，但由于使用的方 程是近似上行波方程，有倾角限制，精度略差，由差分带来的频散难以消 除。当速度没有横向变化时，相位移算法是一种精确有效的偏移算法。上 行波方程在频率波数域( 即珊一k 。域) 可精确求解，但在曲一k ，域无法处 理横向变化的速度。为了克服这一缺点，二十多年来先后产生了一系列算 法。如相位移加内插偏移算法( p s p i ) p j ，该算法允许横向变速；裂步法深度 偏移算法 6 1 ，该算法把相位移法与透镜方程( 也称5 度方程) 相结合，可以适 应较大的横向速度变化：傅立叶有限差分深度偏移( f f d ) 算法【7 1 ，f f d 方法 将相移法和一x 域有限差分法相结合，保持了两种算法的优点，通过增加 有限差分项提高裂步法深度偏移近似的精度。f f d 方法可以适应较剧烈的 横向速度变化，是目前较为流行的一种深度偏移算法。其基本做法是：在 每一延拓步将速度场分成两部分：常速速度场和变速速度场；向下延拓算 子也分成两部分：处理不变速度的相移算子和处理速度变化的有限差分算 子。论文第三章介绍了这种算法。 地震的正演模拟技术与偏移相反，它假设地下地质模型己知，在此条 件下，模拟地震波的传播和地震勘探的采集过程，以得到理论的人工合成 地震记录，把人工地震记录与野外实际的地震记录进行比较可以检验地震 解释的结论，修正地震资料地质解释方案。在此基础上进行地震波的最优 拟和约束反演，用以预测含油、气地层的横向分布。为地质钻探与油、气 开发提供科学依据悼j 。地震正演模拟也是研究、评价地震处理方法的基础。 地震正演模型法主要包括波动方程模型法和射线模型法。由于波动方 程从本质上描述了地震波的传播，地震波的旅行时间、振幅、相位、速度、 频率等等动力学信息反映了地下地层的各利，几何学特征及物理性质，所以 2 查堡塞三查堂翌主笙苎! ! 三 经典的波动理论在地震勘探中得到了广泛的应用。求解波动方程是波动理 论模型法的主要手段。波动方程模型算法主要有积分法【9 】、有限差分法”1 、 傅立时变换法郛有限元法l 】2 】。积分法速度快，理论简单，方法直观，但 精度不高，应用并不广泛。有限差分法一经提出便因其快速、有效性得到 广泛的关注。但这种方法为了保证一定的精度需对模型过分采样，从而消 耗大量的计算机时间。采用高阶劳分技术【t3 和交错网格差分技术，可以 提高计算精度，相应地也会增加计算量。有限元法比较适宣于地震波的模 拟，它可以模拟任意地质形态，也能处理波动方程强间断解，但其严重缺 点是运算量大，所需计算机内存大，因此计算成本较高，致使实际应用受 到严重限制。傅立时变换法是用傅氏变换求鳃空间导数，故精度高。其最 大优点是最小波长所需网格少，但当地下界面比较复杂时，大距离采样会 产生人为绕射。为了消除这种绕刺，产生了很多办法，最简单的方法是平 均弹性常数，但结果不那么令人满意。之后又发展了等介质概念，这种技 术不适子解决多层趋界面的情况。f o r n b e r g 提出在曲丽格上做数值计算的 想法1 1 “。但是f o r n b e r g 的映射技术只能使网格沿着主要的界面变化。 n i e l s e ne ta l 采用了计算流体力学中常用的分块映射技术【1 6 使网格可以 沿着所有界面变化，提高了算法的精确性。 本文以n i e l s e n 的曲网格映射技术为基础，对二维声波模型和弹性波 模型进行了模拟，并给出了与传统伪谱法比较的结果。论文第一章介绍了 伪谱法模型的算法，并对算法进行了稳定性分析，推导了稳定性条件并作 了误差估计，同时具体说明了如何生成曲网格，完成直角坐标系到曲坐标 系的映射。第二章在曲网格上运用伪谱算法计算了几个具体的地震模型， 包括声波模型和弹性波模型，并导出曲网格伪谱法的稳定性条件。 综上，本论文主要完成了如下工作： 1 从理论上阐明了伪谱法模型算法，分析了算法的稳定性，捺导了馥 网格伪谱法声波方程稳定性条件，讨论了反映射，改进了f f t 计算，给出 了简单的误差估计，完善了n i e l s e n 的算法。 2 利用超限插值和分块技术生成曲网格。 3 用曲网格伪谱法原理，编制了相应的声波方程和弹性波方程模型程 序，并在计算机上实现了曲刚格伪谱法计算。计算了几个具体的声波和弹 性波模型，并与传统伪谱法的结果进行了比较。 大连理工大学顺士论文 刖舌 4 阐明了二维叠后f f d 偏移算法原理及有限差分项优化系数的求法。 5 利用实际算例比较了四种不同偏移算法，即相位移法，4 5 。方程有 限差分法，裂步法，f f d 法的结果，编制了相应的程序，以此，比较了各 种算法的优缺点及时间偏移与深度偏移之间的差别，进一步验证了f f d 算 法的有效性。 大连理丁大学硕士论文曲删格伪谱法的原理及算法 第一章曲网格伪谱法的原理及算法 第一节傅立叶变换法 傅立叶变换法也称伪谱法( p s e u d 。s p e c t r a l ) ，是用快速傅立叶变换 即f f t 来计算波动方程的空间导数项，用中心时间差分来进行递推。其主 要优点是可以在大空间网格上得到精确的波场值。因为每个最短波长仅需 数量很少的网格数，从而降低了计算存储量和c p u 时间，所以被广泛应用 于波场模拟。 一傅立叶变换法模型的算法 已知声波方程【1 6 】： p 詈一罢怕小，“，。欲 敏一：， p 詈一考 办，“，西出 “， 印，抛。却 a ta xa z ， ( 1 1 1 a ) ( 1 1 1 b ) f 1 12 ) 其中p 2 一盯表示压力( 盯表示应力) ，“和v 表示位移速度分量，p 表示 密度，兄( = p c 2 ，c 表示波速) 表示l a m e 常数，屯表示震源。 ( 1 1 ，1 ) 一( 1 i 2 ) 式的右端的空间导数项用傅立叶变换计算，因为 p ( x , z , t ) 2 荨p 懒p g 2 善比删e x p q k z 杰( 1 1 3 a ) 其中i = 一以， 芦( 屯，z ，f ) 和卢( x ，砭，) 分别为p ( x ，z ，f ) 对变量。和：的 傅立叶变换： 声( t ，z ，) _ z ，p ( 工，z ，懒p ( i kx ) ， ( 1 1 3 b ) 大连理工火学硕士论文 曲网格伪谱法的原理及算法 芦( x ，七) ：yp ( 耶，t ) e x p ( i k p ( i k z z ) 盹，) 5 乙p ( 耶， 散射关系式为： 甜：c ( ；+ 女；) 1 坨， 于是我们有： ， 老= 蛳咖则纠 k 2 ( 1 1 3 c ) ( 1 1 4 ) f 1 1 5 a ) ( 1 ，1 5 b ) 抛a v 同理可得n - ，西。方程( 1 1 1 ) - ( 1 1 2 ) 左端的时间导数可以用时间差分求 买际上我们璺得剑的仅是坡场j 土力p ，所以口 把万程( 112 ) 改写成： 警= 五夏0r 一1 瓦0 p p+ 墨p 卜丑尝( 一上p 塞+ 曼p ) ， ( 1 1 - 6 ) a 积幽由m l 1 1 1 ”j 由差分格式： _p(t+at)-2p(t)+p(t-出)。以晏(害+垒)一丑昙(三害+兰) 矿 m p 出p也p 赢p 、 ( 1 ，1 ，7 ) 得递推公式： p c h z z p 一加一a t ) + & 2 卜丑昙c 一古老+ 詈，一丑鲁c 一去妾+ 詈刀 wu ( nnmn ( 1 1 ，8 ) 这样由波场在，及t 一5 t 时刻的值，可近似推得波场在f + a t 时刻的佰。 二算法在计算机上的实现 力k啾啪 z x融峨 k | i 印一舐 查堡型三查兰堕主丝兰 堕堕塑堡坐型型鲨翌兰! ! 塑墨 用x 表示地表方向，z 表示地下深度。分别将x ，z 方向等间距地分成 n x ，n z 个点， x = m 越，小= 0 , 1 ，聊一l ， ( 1 19 a ) = n z ，h = 0 , 1 ，n z l ( 1 19 b ) ) ：，缸是网格间距。由抽样定理出，业应满足条件： m a x ( a x ，a z ) 害坠，( 1 1 1 0 ) o ，m “ 波数k x , k ：的具体采样值为： ( t ) ，= 要( ，) = 一 ；) ，s = 1 2 ，z x 2 1 ， ( 1 1 1 l a ) n x l x x ( 七：) ，= 兰 ( 尼：) 孵，= 一( t ：) ，= l 2 ，抛2 _ - ， ( 1 1 1 l b ) 为了提高计算效率，在做快速傅立叶变换( f f l 、) 时用到了实信号二维离散谱 的共轭性质【1 钔： p r ，j 2 p r ，一s 其中p 是压力在( x rz 。) 点的离散谱，n ，s 是p 的共轭。这样在进行傅 氏变换时只需计算全部数据的一半即可，而另一半数据的谱可取其共轭得 到。利用傅氏变换的线性性f l8 1 ，可以用一个复序列的傅氏变换计算两个实 序列的傅氏变换。另外，在做多个序列的傅氏变换时，许多运算是重复的， 如序列序号倒置的计算，运用类似向量化的运算，所有序列序号倒置的计 算只需计算一次即可，这样可以大大减少重复运算量，提高计算效率，节 省计算存储空间。 先用f f t 沿着x 方向对矩阵( p 。) ，( 。) ，( 。) 进行傅氏变换。当然， 由提到的傅氏变换的技巧只需计算一半的数据。变换的结果与水平方向的 波数相乘，再进行反傅氏变换，就得到了在所有x = x i 点的相对于x 的空 间导数。可用类似的过程计算z 方向的空间导数。再根据( 1 1 8 ) 式可以求得 波场在t + 出时刻的值p ( t + f ) 。 夫连理工大学顺 一论文曲网格伪谱注的原理发算法 三稳定性分析、误差估计 由( 1 ，1 5 ) 式，若x ，k ：为常数，则 妾= i k x z i , ( k 舻e 删 七， ( 1 11 2 a ) 笙o z 也i ( x , k z , t ) e x p l - - j ( 纠_ f 镢即， ( 1 1 1 2 b ) ，f 1i 1 2l k 先不考虑震源，将上式及c c 丑p 带入( 1 1 8 ) 中则有： p ( t + a t 卜2 烈。一烈卜) 十 _ 啦( - i 等) - i k = ( 一鲁) 印( f ) 出 = 2 p ( t ) 一p ( t 一r ) 一c 2 ( t ：+ k ；) p ( t ) a t 2 = 2 一c 2 ( 后；+ 霹) f2 】p ( f ) 一p ( t a t ) f 1 1 1 3 、 若以下标n 表示时刻，则n + l 表示时刻h 出，n 1 表示时刻卜a f 臼= c + 女弛， ( 1 1 “) 则上式可写成： p 月+ 1 = ( 2 0 2 ) p 月一p 月一1 ， 记为矩阵形式： ( 翟卅儿o 】) 。p - 1 j ， ；， 特征方程为 f 2 0 2 一r l | _ o l 1 一引 即：矛+ ( 口2 2 ) 五+ 1 = o 解得特征值五= - ( 2 - o z ) 【( 2 - 0 2 ) 2 _ 4 】1 ，2 ， ( 1 1 1 6 ) 当h l 时，算法就是稳定的。而曰2 时有( 2 - 口2 ) 2 4 ，故五叫以写 大连理工大学硕士论义 曲网格伪谱法的原理及算法 成 取模 = i ； r 2 一目2 ) 【4 一( 2 一0 2 ) 2 1 1 2 i 所以0 2 时h = 1 。 0 d ( 1 1 1 1 ) 式知波数 k x 墨兀x k z 曼冗a z 带入( 1 1 1 4 ) 式，并考虑到0 墨2 ，即得 即当 ( 1 缸) 2 + ( 1 a z ) 2 l 2 c a t 疡 出s 亟璺! 呈! 竺：竺! 时，算法稳定。 ( 1 1 ，1 7 ) 下面给出方程( 1 1 8 ) 的差分误差估计。令万表示p 的理论值，用台劳 级数给出： 酏地，2 器学 由( 1 1 6 ) ，( 1 ，1 8 ) 有： ( 1 1 1 8 ) 如+ a t ) 。2 p ( t ) 一p ( t + ，0 矿2 p 旷 ( 1 1 1 9 ) 于是可得 故 掣 矿 。 d + 出0 、j 一 帅等 十 + d 以 印 = i i “ 大连理 i 大学硕士论文 曲州格伪谱注的原理及算法 删一群等， s 对上式仅取一项时，可得较为实用的最低阶误差估计： 删z 一等譬， n 地。， 第二节曲网格伪谱法原理及数值网格的生成 尽管傅立叶伪谱法可以在大空间网格上得到精确的波场值，但如果界 面是弯曲的，离散化数值计算会将连续曲界面变成阶梯状界面，当在大间 距的直网格上进行傅氏变换时，阶梯状晃面会产生不必要的绕射。为了消 除这些绕射，产生了许多方法，将直网格改成曲网格是一种很自然的想法。 在曲网格上做数值计算的思想最早由f o r n b e r g 15 】提出。但是f o r n b e r g 的 映射技术只能使网格沿着主要的界面变化。n i e l s e np ta l l l 6 】采用了计算流 体力学中常用的分块映射技术使网格可以沿着所有界面变化，提高了算法 的精确性。本文采用了n i e l s e n 的方法生成网格。其基本思想是通过变换 与分块把地下所有曲界面映射成规则的界面，相应地产生一个曲坐标系， 然后在曲坐标系下进行伪谱法模拟计算。 曲的线性网格的生成与应用一度是计算流体力学的主要研究课题。在 过去十几年里出现了若干种网格生成方法，方法主要有两大类：基于插值 和优化的代数方法与基于某种偏微分方程解的微分方法。本文所采用的超 限插值法【1 6 】【2 0 ，是利用四条边界上的点通过插值产生区域内部的各个点， 从而完成直坐标系到曲坐标系的映射。 一基本变换方法 使数值网格沿着地下界面变化是n i e l s e n 的主要思想，为实现这一目 标，可将实际的几何构造映射成某一曲坐标系下的规则图形，比如矩形。 再利用插值法产生区域内部的各个点，从而生成所需的曲网格。具体实现 如f ： 6 查堡型王查堂堕主堡苎 些型堂堕堂鲨堕堕堡! 丝 如图l 一2 1 所示，设所选区域的四条边分别为a d ，b c ，a b ，c d 。将a d 边和b c 边分别对应于曲坐标系舌一叩域下的两条平行于t l 轴的直线a d 巾c ， 类似地，a b 边和c d 边对应于曲坐标系下的两条平行于的直线a b ，c d 。 这样，一个直角坐标系下的不规则区域就映射成了另一个坐标系下的矩形， 相应地建立了一个曲坐标系。然后在曲坐标系下建立一个规则网格，就可 以在规则网格上进行数值计算了。 图1 2 i 生成网格的过程( a ) 直角坐标系下的实际边界( b ) 曲坐标 系下的边界及网格( c ) 直角坐标系下与图b 所对应的网格 一般来说，常常将曲坐标系下矩形的网线的值取成整数，例如取 l = 1 ，2 ，；r l = l ，2 ，曲坐标系下网点的坐标为( 1 ，r 1 ) ，而其对应 的实际坐标系x z 下的值为x = x ( 亭，叩) ，z = z ( 毒，7 ) 。实际坐标值的计算可 以采用一种多维插值方法超限插值法，通过插值函数和区域的四条边 界产生区域内部的各个点的坐标值( 见图1 2 2 ) 。公式如下： 22 ， 22 蟛圳= 姒9 ( r 1 ) + 蹦枷( 翻一一傩砜( 咖( 厶川r ， ”= lm = 1月= lm = l ( 1 2 1 ) r = x e 、- t - z e 。，其中e 。，e z 是直角坐标系下的单位向量。丸，是插值函数，方 便起见，我们取线性插值函数： # 一i# 一# 5 嚣瞄) 2 蠢 ( 玎) = 竖旦，妒：( q ) ：旦丑， r h r , r l , 一玎l 人连理工大学硕士论文曲网格伪潜法的原理及算泣 ，( 岛，叩) 图1 - 2 2 超限插值法利用边界求 ，_ ) ，z ( f ，_ ) 的值 从图1 - 2 - 3 可以看到一个简单的单层模型的映射过程：先将直角坐标 系下的条曲线( 图卜2 1a ) 映射成曲坐标系下的一条直线( 图1 2 3 b ) ，把图 l 一2 - 3 b 用相应的直网格划分，再用超限插值公式( 1 i 2 ) 得到网格点( 等。，7 ，) 在直角坐标系下的坐标值( x ( 夤，玎) ，z ( 善f ，乃) ) ，从而得到了曲n n ( n 1 2 - 3 c ) 。 ( b ) ( c ) 剀卜2 - 3 一个简单的层模型( a ) 实际物理构造( b ) 曲坐标系f 的同一构造 ( c ) 映射后得到的网格 二分块方法 当地质结构比较复杂时，对整个区域直接做耿射有可能造成网格因过 于不规则或过于倾斜而产生畸变。为避免这种情况，并使网格沿所有的地 下界面变化，采用了在计算流体力学中应用广泛的分块的方法：按每一层 曲界面的单调性将曲线分成若干相邻的段，从而把整个区域分成若干块， 人连埋工大学硕士论文 曲网格伪倍法的原理发算法 将每一块单独映射生成网格后再把所有的块儿连起来，形成统一的整体。 连接时应注意网格的连续性。 图1 2 - 4 显示了一个分块过程。图1 2 - 4 a 是一个透镜体构造，将构造沿 地下界面及界面的单调性分块，将每块按前节介绍的方法映射成图1 2 4 b 中各对应子块，再将各子块相连得到图1 2 4 b ，把图l 一2 - 4 b 用直网格划分 后，通过超限插值得到网格1 2 4 c ，曲网格1 2 4 c 的网格线沿着所有界而 变化。 曲网格伪谱法就是将x z 坐标系下的声波方程转化到曲坐标系毒一叩 下，在曲坐标系下进行伪谱法计算，那么，在直网格( 1 2 - 4 b ) 上的计算等价 于在曲网格( 1 2 - 4 c ) 上的计算。 图l 一2 4 透镜模型的分块图 ( a ) 实际物理构造( b ) 曲坐标系下的同一构造 f c l 映射后的网格 大连理工碗士论文 曲网格伪谱法模型试算 第二章曲网格伪谱法模型试算 第一节曲网格伪谱法声波模型 在用伪谱法模拟地震模型时，离散化数值计算会将连续曲界面变成阶 梯状界面，大距离采样会在阶梯状界面产生不必要的人为绕射。为了消除 这些绕射，可以通过分块与变换把地下所有曲界面映射成规则的界面，相 应地产生个曲坐标系。在曲坐标系下地下界面是平的，离散网格也是规 则的，通过插值可以得到实际模型的离散网格，网格线随地下所有界面变 化。伪谱法模拟计算是在曲坐标系下的规则刚格( 直网格) 上进行的，因 为曲坐标系下的规则网格对应于直角坐标系下的曲网格，所以称这种方法 为曲网格伪谱法。 一声波方程 直角坐标系下声波动力学方程为 p 詈一署小刀， 。西叙， ，詈一票坞( 二一， a la z 。+， o p ：一五丝一五竺! a ta xa z ， ( 2 1 1 a ) ( 2 1 1 b ) f 2 1 ，2 ) 其中p = 一d 表示压力( o - 表示应力) ，“和v 表示位移速度分量，p 表示 密度，丑( = p c 2 ，c 表示波速) 表示【，a 吣系数。s x ，s z 代表震源。 s a c ( x ，z ) = 导厂( x ，z ) 矗( 出 s z ( x ，：) ：昙，( x ，= ) ( ，) o ( 2 13 a ) f 2 1 3 b ) 大连理工顺上论文 曲删格伪曹法模型试算 这里，h ( t ) 是雷克子波，f ( x ，z ) 是空间分布函数 f ( x , z ) = - e x p 一坚正等蔓】， (213。)2c o 一 式中。表示脉冲的半宽。( x o ，z o ) 表示震源位置。 因为计算是在曲坐标系掌一，7 下进行的，所以要将x z 坐标系下的 ( 2 11 ) 一( 2 1 3 ) 式转换到曲坐标系毒r l 下。有 ，塑一一卢塑 一o t b a f ，7 。罢+ s 。( 毒，7 7 ) d 7 7 p 祟= 一手：罢仉罢+ s z ( 毒，呷) ， p 百一乞壶仉翥+ ( 善，呷) ， _ o a p ，2 ( # a o 亨u + 玑丽o u + 关嚣+ 仉塑o r ) 系数最，乞，巩，砚由隐函数求导法则导出 x = z j 2 _ q f 一= 一鱼 。 j z r 仉。一_ 佗1 ，4 a ) ( 2 1 4 b ) ( 2 15 ) ( 2 ，1 6 a ) ( 2 16 b ) ( 2 17 a ) ( 2 17 b ) r 2 1 7 c ) v 历 缈纱 玑 仉 + 可鸳 望西 六 六 l l i l 野 玎 ，j # s f 人连理工硕士论文 曲唰格伪请法模型试算 x f 7 7 ：2 寸 这里，是j a c o b i a n 行列式 ( 2 ，i ，7 d ) 而x ，x 口，z # ，匆可在网格生成过程中通过对( 1 2 1 ) 式分别关于孝，誓求导数得： 毒吣莉2 喜毒以c 麴r c 己砌+ 耋y 翩，毒吣帆， 一景磊( 毒帆( ”) r ( 磊，繇) ， ( 2 i ，9 ) = lm = lo ， 南吣朋2 喜丸( # ) 南r ( 朋+ 嘉南( 咿( 锄j 喜耋珧) 南蹦护( 翻小 ( 2 1 1 0 ) 如图1 2 1 所示映射厝所在曲坐标系下的网格是一个嫂则网格，方程 ( 2 1 4 ) ( 2 1 5 ) 中的空间导数的值就可以在规则网格上的各个点上用标准 的傅立叶伪谱法计算了。方程的边界条件是傅立叶伪谱法的边界条件，即 周期边界条件，但周期性由x z 域转换到了毒一r 域。方程的时间导数项 可以用中心差分格式来计算，精度略差些。 可以看到曲坐标系下的方程( 2 1 4 ) 一( 2 ，1 5 ) 中的空问导数项比直角坐标 系下的方程( 2 1 1 ) 一( 2 1 2 ) 的空间导数项多一倍，相应的也会增加计算量， 但用传统的伪谱法在直网格上( 直接在直角坐标系下计算) 要想获得与曲 网格上( 映射后在蓝坐标系下计算) 相同的精度，就需要加密网格，也会 增加计算量。因此在处理大区域模型时，曲网格方法仍可以节省计算存储 量。 二稳定性分析 前一章推出了伪谱法的稳定性条件是 大连理工顺上论文 曲同格伪者法模型试算 蜓堕型垒：垒! z z m a x c ( 2 i 1 1 ) 其中a t ，a x ，a z 分别表示时闯与空间上的采样间隔。但因为曲网格伪谱 法的数值计算是在曲坐标系下进行的，波动方程由直角坐标系下的形式 ( 21 】) ( 2 i 2 ) 转变为曲坐标系下的形式( 2 ，1 4 ) ( 2 1 5 ) ，相应的差分格式不 周，故稳定性条件也不周。 由式( 1 ，1 ，8 ) 可推得曲坐标系下波动方程( 2 1 4 ) 一( 2 1 5 ) 的差分格式： 烈h 出) 2 烈f ) 呻o _ f ) _ ，叫鼻善卜古( 享。詈+ 矾器 + 仉南e 一击c 炙妻+ 玑考+ 彩，+ 参奏e 一万1e 吾：奏+ 毒+ 蹦， 叼= 南( - 专( f = 妻+ 仉赛 功 不考虑震源子波，由傅立时变换的微分性质( 1 1 1 2 ) 得： p ( f + f ) 2 p ( t ) 一p ( f 一f ) 一。 ( 鬈+ 亭z 2 0 2 + 2 ( 靠h + 乞砚) 吩吻+ ( 掘+ 谚) 孑 出2 p ( f ) = ( 2 0 2 ) p ( f ) 一p ( f 一f ) 其中c = 、仁表示波速。 vp 移= 罐暖十彰2 乍2 + 2 ( g r x + 乞琵) 嗨岛+ ( 砚2 + 磁2 j 勺2 2 盘， ( 2 1 1 2 ) 由第一章的讨论知：臼s 2 时算法稳定，将石考，sr e a t 代入 ( 2 1 1 2 ) 式，并令其满足条件0 2 。 若掌= 7 7 贝u 有： a t - 1 1 垒 万m a x 【( 炙+ 巩) 2 + ( + 魂) 2c 故曲网格伪谱法的稳定性条件是： 人连璀_ = 倾十论文 曲酬格伪谱法模型试箅 ， 2 m i n ( a 亭，a t ) a t 2 = = = = = = 兰= = 兰= = = = = 一 m a x ( 最+ 叩，) 2 + ( 亏：+ 叩二) 2c 显见，式( 2 1 1 1 ) 是式( 2 i 1 3 ) 的特例。 三数值算例 ( 2 1 1 3 ) 算例1 ：阶梯构造模型 这个例子通过在曲网格和直网格上对阶梯构造模型的模拟的比较可以 看出曲网格伪谱法很好地改善了直网格伪谱法中的绕射现象。图2 一l l 构 造了一个阶梯形的界面和一个平层，震源位置( x 0 ，z o ) = ( 6 4 0 m ，5 5 m ) ，检波点 在地表。运用第一章第二节介绍的方法将其转化成两个平层，从而把波动 方程( 2 1 1 ) 一( 2 1 2 ) 映射到曲坐标系下。数值计算用的曲坐标系下的网格 大小是( 叩f = 2 5 6 1 2 8 ) ，空间采样间隔孝= a t = 5 m ，时间采样间隔 a t = o 5 m s ，应当指出，在曲网格中空间采样间隔是随地层变化的，见图 2 一卜2 。震源脉冲半宽0 9 = 5 m ，取密度为单位密度。 图2 卜3 和图2 一卜5 是阶梯模型在直网格上的合成记录，可以清楚的 看到因为离散化数值计算造成的人为绕射。图2 - 1 - 4 和图2 一卜6 是阶梯模 型在曲网格上的合成记录，显然，人为绕射已经被很好地消除了。 算例2 透镜模型 图2 一l 一7 构造了一个透镜体模型。震源位置( z o ，：o ) = ( t 9 2 0 m ，1 9 5 m ) ，震 源脉冲半宽c o = 1 5 m ，检波点在地下3 1 5 米。雷克子波的频率是f = 2 0 h z ， 延时8 0 m s 。图2 1 8 是透镜体模型的曲网格。数值计算用的曲坐标系下的 网格大小是( 叩孝= 2 5 6 1 2 8 ) ，空间采样间隔蟛= a t = 1 5 m ，时间采样间隔 a t = l m s 。 图2 一卜9 和图2 一卜1 0 分别是模型在曲网格上和在直网格上的合成记 录，人为绕射的对比十分明显。图2 卜1 1 至图2 - 1 - 1 4 是曲网格上做的波 场快照，时间分别是l o o m s ，2 2 0 m s ，3 3 0 m s ，3 8 0 m s 。图2 卜1 5 给出了在 直网格上做的波场快照，时间是3 8 0 m s 。 大连理工顺士论文曲网格伪谱法模型试算 四讨论 因为曲网格伪谱法是在映射后的规则网格上进行计算，所以合成记录 也是在映射后的规则网格上得到的，而实际的检波点位置却不一定是在 眚一叩域网格点上。换句话说，计算得到的善一叩域下的合成记录，要返回 到z z 域下的整样点上，即要把数据从善一叩域反映射到x z 域，以得到 实际所需要的合成记录。这一问题，可以用插值的方法来解决。 最一般的情况是网格在x ，z 方向上都有变化( 如图2 一卜8 ) ，所以需要 在两个方向上插值，先在x 方向插值，获得x 方向上的整样点的波场值。 即：已知波场压力p ( x ，z ) 在( ，z o ) ，k = l ，2 ，n x 处的值，其中深度z o 是榆 波点的深度，要求构造一个次数不超过n x 一1 的代数多项式 芦月( x ) = “h x “+ d n l x “一1 + + “o 胛s ，? x l 使磊( x ) 在节点处满足p n ( 张) = p ( x k ) ，k = l ，2 ，。 i x 但在实际应用中，高次插值并不可取。尽管，节点的增多使插值多项 式在更多的点与| 口( x ) 一致，但两个插值节点之间，其逼近程度不一定很好， 通常是分段做低次插值。 选用二次代数插值的拉格朗日公式【2 1 】如下： 22 西( 加z p ( ，= 1 n 罴， 西( x ) 2 o ) o ( ，o ( 。) 2 l j 羔，r 2 1 1 4 1 j = o f = 0 。 。 i 由( 2 1 ，1 4 ) 可得x 方向上的整样点的波场值。再用类似的方法在z 方向上插