（地球探测与信息技术专业论文）初至波层析成像应用研究.pdf

摘要 近地表层高程、岩性、速度、构造等对于地面地震勘探，会引起地震剖面相位的畸 变，极大地影响着叠加剖面的成像和构造形态，是制约地震勘探的重要因素。本文针对 近地表情况，利用初至波传播路径和走时信息，通过合理的走时反演得到近地表介质的 速度信息。 文章围绕着初至波层析成像，回顾了射线追踪和层析反演的发展历程；根据地震波 传播的基本原理，对初至波在表层模型中的传播规律及其特征进行了分析；编程实现了 基于f e r m a t 原理的最短路径( s p r ) 射线追踪方法，最短路径方法是用网络节点之间的 最小旅行时连线近似地震射线路径，这种方法可以同时计算出从震源到达空间所有点的 初至走时及相应的射线路径，并且不受射线理论约束，所得走时保证了全局最小的特点； 对比了两种使用的反演方法：代数重建技术( a r t ) 和联合迭代重建技术( s i r t ) 。并且实 现了编程。 文章最后建立了初至波层析反演的流程，通过理论模型试算，反演得到了可靠地表 介质速度，验证了最短路径射线追踪和s i r t 算法的正确性和可行性。对初至波与近地 表的对应关系有了明确的认识，能够更合理更准确的利用初至波信息来解决近地表问 题，为解决复杂地表地震勘探资料处理奠定了良好的基础。 关键字：初至波最短路径射线追踪层析成像 a b s t r a c t f o rt h es u r f a c es e i s m i c s u r v e y , e l e v a t i o n ，l i t h o l o g y , v e l o c i t y , s t r u c t u r e ，e t c o f n e a r - s u r f a c ew i l ll e a dt os e i s m i cp r o f i l e sp h a s ed i s t o r t i o n ，g r e a t l yi n f l u e n c e st h ei m a g eo f s t a c k e ds e c t i o na n dt h es t r u c t u r e ，i sa ni m p o r t a n tf a c t o rc o n s t r a i n i n gs e i s m i ce x p l o r a t i o n t h e s i su s e dt r a v e lt i m ea n dp a t ho fs e i s m i cf i r s tb r e a ka i m e da tc o m p l e xn e a r - s u r f a c e p r o b l e mi ns e i s m i cp r o s p e c t i n g ，p r o v i d i n gv e l o c i t yi n f o r m a t i o nt h r o u g hr e a s o n a b l et r a v e l t i m ei n v e r s i o n f i r s t l y , t h i sp a p e rb r i e f l yr e v i e w st h ed e v e l o p m e n tp r o c e s so fr a yt r a c i n ga n dt o m o g r a p h y i n v e r s i o na r o u n ds e i s m i ct o m o g r a p h yo ff i r s tb r e a k f u r t h e r m o r e ，t h ep a p e ra n a l y z e dt h e r e g u l a r i t yo fp r o p a g a t i o na n ds t y l ec h a r a c t e ro ff i r s tb r e a ki ns u r f a c em o d e lb a s e do nt h e r e g u l a r i t yo fs e i s m i cw a v ep r o p a g a t i o n p r o g r a m m i n gi nt h i st h e s i sr e a l i z e dt h em e t h o do f s h o r t e s tp a t hr a yt r a c i n gb a s e do nf e r m a t sp r i n c i p l e t h i sm e t h o di st ou s eb e e l i n e so ft h e m i n i m u mt r a v e lt i m eb e t w e e nt h en e t w o r kn o d e st oa p p r o x i m a t es e i s m i cr a yp a t h s ，a tt h e s a m et i m e ，i tc a l lc a l c u l a t et r a v e l t i m ea n dc o r r e s p o n d i n gr a yp a t h so fa l lp o i n t sf r o ms o u r c et o s p a c e t h ei m p o r t a n t i st h a t r a y - t h e o r y c a nn o tc o n s t r a i n ti ta n dt r a v e l t i m ea r e m i n i m u m f i n a l l y , w ed e s c r i b et w oi n v e r s i o nm e t h o d s ：a l g e b r a i cr e c o n s t r u c t i o nt e c h n i q u e ( a r t ) ，s i m u l t a n e o u si t e r a t i v er e c o n s t r u c t i o nt e c h n i q u e ( s i r t ) ，r e a l i z e dp r o g r a m m i n g a tt h ee n do fp a p e r ，w ei n t r o d u c et h ep r o c e s so ff i r s tb r e a kt o m o g r a p h yi n v e r s i o n ，b y u s i n gt h et h e o r e t i c a lm o d e lt ot e s ta n db ya p p l y i n gt h em e t h o do ft o m o g r a p h yi n v e r s i o no f f i r s tb r e a kt r a v e l - t i m ei ns u r f a c em o d e l ，w ea c q u i r e dr e l i a b l es u r f a c es p e e dm o d e la n dp r o v e d c o r r e c t n e s sa n df e a s i b i l i t yo fs h o r t e s tp a t hr a yt r a c i n ga n ds i r ta l g o r i t h m w j t l lac l e a r u n d e r s t a n d i n go fr e l a t i o n s h i pb e t w e e nf i r s tb r e a ka n dn e a r - s u r f a c e ，w ec o u l du s em o r e r e a s o n a b l ea n dm o r ea c c u r a t e l yt h ei n f o r m a t i o no ff i r s tb r e a kt os o l v en e a r - s u r f a c ep r o b l e m ， l a y i n gag o o df o u n d a t i o nt od e a lw i t hc o m p l e xs u r f a c es e i s m i cd a t ap r o c e s s i n g k e y w o r d s ：f i r s tb r e a k ，s h o r t e s tp a t h ，r a yt r a c i n g ，s e i s m i ct o m o g r a p h y i i 长安大学硕士学位论文 1 1 选题的依据 第一章绪论 我国西部的沙漠、戈壁、黄土塬、山地等复杂地表区，表层地质条件极为复杂，地 形起伏变化大。近地表介质的各向异性对波场造成不等量的延迟，导致反射波形态畸变， 影响叠加成像、信噪比、分辨率，以及成像位置的精度。 近地表从炮点出发，经过介质首先到达的初至波，可能包含直达波、回折波，透射 波、折射波或多个折射波的组合，这些波主要在近地表传播，往往能量很强，容易识别， 可追踪性好，其走时包含了近地表层速度信息，利用初至波层析成像方法构建近地表速 度模型成为研究热点。初至波走时反演层析成像是把介质离散成一系列小单元，利用射 线追踪方法追踪各个单元内的射线路径和走时信息，利用射线路径和走时反演每个单元 内介质的速度，最终得到近地表介质速度模型。由于该方法不需要对折射界面做任何地 质假设，利用全偏移距范围的初至时间，对起伏地表情况下的初至拟合好，能对复杂近 地表介质成像。 1 2 现有方法评述 1 2 1 射线追踪正演的方法现状 射线追踪方法是指给定激发点和接收点的坐标位置和地下介质的速度分布，计算激 发点到接收点的射线路径和走时。地震层析成像的兴起，引起了人们对射线追踪算法的 研究。射线追踪方法的计算速度和计算精度，直接影响层析成像的效果。下面对具有代 表性的射线追踪方法进行简述并对各种方法进行比较。 1 试射法和弯曲法是两种经典的射线追踪方法。试射法，即试着给出射线出射方 向，求射线路径，射线一般不能到达接收点的位置，需要根据实际到达点接收点之间的 偏差，修改射线出射方向，再次追踪射线路径。重复以上步骤，直到射线接近到达接收 点位置。弯曲法是运用波在介质中的曲线传播，经过整体计算和边值运用，解方程得出 路径。 2 v i d a l e ( 1 9 8 8 ) 提出了基于波前扩张的思想用有限差分求解程函数进行射线追踪。 首先用矩形网格将介质离散化，从激发点开始，地震波以方阵形式按网格节点次序逐点 向外扩展，用有限差分方法层层求解程函数方程，得到网格节点上的旅行时，然后根据 第一章绪论 地震波传播原理，射线路径与波前面垂直，从接收点沿旅行时最大梯度方向返回到激发 点，就得到激发点到接收点之间的最小旅行时射线路径，这种方法用已知求出的内圈网 格点上的旅行时来求相邻的外圈网格点上的旅行时，且按矩形波前面向外层层推进【3 引。 3 基于旅行时线性插值方法( l t i ) ，是把模型离散成均匀的正方形单元，旅行时 和射线路径的确定只与单元边界上的点有关。通过满足f e r m a t 的线性插值，使精确点 与边界点相交，单元内部的射线路径为直线，是一种精度较高的射线追踪方法。 4 最短路径方法起源于网络理论。最短路径射线追踪方法基于费马原理和网络中的 图论理论，它用一系列规则的空间离散网格和节点代替所要追踪的介质模型，介质的速 度用每个节点上的速度代替。并根据惠更斯原理依次将空间离散点作为子波震源，以有 限的离散方向来代替子波射线的全方位传播，逐步获得整个空间离散点的全局最小走时 和路径。通过不断研究，最短路径方法目前发展较为成熟，其基本算法的计算程序也较 为固定。 试射法适用于介质不复杂和满足一定误差条件的要求的情况，用最靠近接收点位置 的射线最为实际射线路径。当介质复杂时，很难确定震源点到接收点的射线路径，而且 计算量大。弯曲法相对于试射法效率高，但计算效率低，又是会陷入局部收敛。试射法 不能对折射波河射线理论不成立的射线路径进行追踪；弯曲法在震源和接收点距离较远 时效率会很低。 v i d a l e 法依据平面波理论，只适用于波前向外层层推进的情况，有时候初至波是 迂回传播的，这种情况下将得不到全局极小旅行时波前。因此，v a nt r i e r 和s y m e s ( 1 9 9 1 ) 修改了v i d a l e 的算法，实用逆向有限差分求解程函方程，避免了依赖模型的映射次序， 使其完全矢量化，这样提高了计算效率，但是，与v i d a l e 方法一样，仍不能保证得到最 小旅行时路径。 l t i 算法似v i d a l e 的一种高级形式，把模型离散成矩形网格单元，射线路径和旅行 时的确定只与单元边界上的节点有关系。该方法分为向前处理和向后处理，在向前处理 中，只计算各单元边界上的节点的旅行时，在向后处理中，根据f e r m a t 原理追踪射线 路径。这样确定的射线路径不是单元边界上的节点的连线，而是满足最小旅行时条件的 边界上的连线，因此，该算法会比其它算法精度更高。 最短路径射线追踪方法基于f e n n a t 原理和图论中的最短路径理论计算出最短射线 路径和最小接收点走时。该方法能够避免一些潜在的射线路径的丢失，当也有可能会有 几条射线路径同时到达接收点，很难选取，而且这种方法计算量大，需要很密的节点和 2 长安大学硕士学位论文 较大的搜索半径，才能得到精度高的射线路径。 1 2 2 近地表地层调查现状 近地表是指地表以下未成岩的低速介质，虽然厚度不是很大( 几十米，个别区域 达到几百米) ，但对地震波场改造很大，不利于地震勘探采集数据和处理。地表介质由 于沉积年代和沉积时间的差异、环境的温度变化、岩性不同、含水量的多少，压实程度 等，使表层介质地震波传播路径和传播速度变化很大。在大部分地区，由于地表受连续 或间断的沉积，低速带具有良好的成层性，有明显的速度界面。实际上低速带的划分没 有严格的定义，通常将表层按速度划分为低速层( 速度小于1 0 0 0 m s ) 、降速层( 速度在 1 0 0 0 2 0 0 0 m s 之间) 、高于2 0 0 0 m s 的介质归类为高速介质。近地表沉积的介质相对更 深层的地下介质，沉积年代相对较短，长年的分化作用使近地表的介质无胶结或半胶结， 地层中含水与不含水，含水量的多少都会引起地震波场的变化。另外，近地表厚度和速 度的各向异性、地表起伏都会对地震波场造成些影响，不利于地震数据的采集和处理。 产生对近地表的因素主要有： ( 1 ) 地表高程变化：由于地表的起伏引起的延迟时变化，地表高程变化是近地表异 常的主要因素。相对来说，地表的起伏造成的近地表延迟时往往大于近地表各向异性产 生的延迟时差，特别是在山地或山前地表高程变化大的地区。 ( 2 ) 近地表横向厚度的变化：由于近地表介质厚度变化，不同的地表条件低降速带 厚度差异很大。山前或山地表层变化很大，有些地方老地层出露地表，不存在低速介质， 有些地方又沉积很厚；戈壁砾石区通常低速层较薄，相对厚度变化不大；沙漠区则不同， 若存在潜水面，由于含水后沙层的厚度变高，可将潜水面作为低降速带界面，若没潜水 面或潜水面很低，这时低降速带在厚的区域可达几百米。 ( 3 ) 近地表横向速度的变化：在近地表由于沉积环境的变化，或沉积物的物理特征 发生变化，如是否含水，胶结不胶结，是砾石还是泥沙等，这些差异决定了地震波的传 播速度差别很大。 1 2 3 地震初至波层析成像方法现状 近地表速度信息的确定是解决复杂地区地震勘探问题的关键，利用初至波走时层析 成像重建近地表速度模型，是确定近地表速度信息的有效手段。近地表速度复杂模型接 收到的初至波可能是直达波、折射波、回折波等多种波的组合。传统的折射波成像，是 第一章绪论 基于层状介质模型来表示实际介质模型。当地表结构复杂时和初至波难以识别时，折射 波成像方法常常是不实际的和无效的。 随着基于f e r m a t 原理的多种射线追踪方法的使用，回折波层析成像成为层析成像 研究的热点。回折波层析成像能同时考虑直达波、回折波、透射波、折射波等初至波， 可以较好的模拟介质横向和纵向的速度变化。朱等1 提出了弯曲射线层析成像确定近 地表速度的层析成像方法假设介质速度虽深度线性变化，介质被离散成矩形单元，每个 单元速度为常速，用约束阻尼同时迭代重建技术反演。李录明等1 钉把近地表模型离散 成矩形单元，单元内的速度用单元节点速度的双线性函数表示，用最短路径射线追踪方 法计算射线路径和走时，用带阻尼的最小q r 分解算法反演。张建中等1 采用双线性 函数表示速度单元，用l s q r 算法反演，取得了一定的效果。 在反演问题中，在传统的线性迭代反演中，反演方法容易陷入局部极小解。非线性 反演局部优化反演方法逐渐得到了关注，但这种方法依赖于初始模型的选取，模拟退火 法、遗传算法等全局优化方法克服了上面方法的确定，不依赖于初始模型的选取，但这 类算法的效率太低，极大影响了全局优化的广泛应用。 1 3 本文的主要研究内容和研究思路 本文详细地分析了地震波传播规律、最短路径射线追踪原理、地震层析成像的基本 理论的基础、利用最短路径和s i r t 进行射线追踪和层析反演，对复杂的近地表模型进 行了正反演模拟。论文的主要内容包括： ( 1 ) 第一章绪论：对射线追踪、层析反演、初至波层析成像、近地表层的调查的发 展、方法、现状进行分析研究。 ( 2 ) 第二章：对地震波传播过程中的一些基本原理和初至波传播规律的理论基础进 行分析。 ( 3 ) 第三章：介绍最短路径射线追踪方法原理，包括速度建模和网格建立，节点旅 行时的计算和射线路径的确定，以及存在的问题及一些改进的方法。 ( 4 ) 第四章：对地震波旅行时层析成像的原理和数学公式进行推导。 ( 5 ) 第五章：设计包括水平界面、倾斜界面、弯曲界面、断层、地堑、起伏地表等 理论模型，进行正反演计算。 ( 6 ) 结论。 4 长安大学硕士学位论文 第二章地震波传播理论基础 地震勘探是以研究地震波在岩层中的传播规律为基础的，由于实际介质的非均匀性 和各项异性，波动方程的边界条件极为复杂，很难解出波动方程。当波速在一个范围内 变化很小时，射线理论是波动方程高频近似，这样就可以使数学问题简化，能方便的研 究地震波在介质中传播的路径、速度和旅行时等。地震波勘探理论一般把介质做一下几 种近似： 均匀介质模型：速度值不随空间坐标变化的介质。 层状介质模型：包括水平界面、倾斜界面、弯曲界面是地震学中最常用的地质模型， 每一层中地震波速度是相同的，不同层之间速度不同。 连续介质模型：地震波速度是空间连续变化的介质的函数，是层状介质的一种极限 情况，当层状介质无限增加，每一层厚度无限减少时，层状介质就变成连续介质。 层状连续介质模型：有些地区地下存在好几套岩性不同的地层，整体上是层状介质， 每一层中又是连续介质。 2 1 惠更斯( h u y g e n s ) 原理 地震波在介质中传播时，波的传播依赖于介质中各质点之间的相互作用。任一点的 震动，将引起相邻点介质的震动，介质中任一波面上的各点，都可看成是产生球面子波 的波源，在其后的任一时刻，这些子波的包络面构成新的波面。 2 2 费马( f e r m a t ) 原理 f e r m a t 原理是从射线角度描述波特点的。地震波在任意介质中沿射线的传播路径， 与其它路径相比，射线路径为所花时间最小的。根据这个原理我们可以得出，在均匀介 质中地震波沿直线传播，在非均匀介质中，地震波不再是直线，时间与波前面垂直的射 线路径。 2 3 斯奈尔( s n e l l ) 定理 根据s n e l l 定理，透射波在第二介质中的偏折，与两种介质的速度比有关。令k 代 表第一层介质，圪代表第二层介质，当k 屹时，透射波靠近法线偏折，当k 圪时， 第二章地震波传播理论基础 透射角大于入射角，透射角远离法线偏折。当入射角不断增大时，透射角也随之增大， 当入射角增大到一定角度时，透射角增大到9 0 0 ，上式可写为 粤：掣 ( 2 1 ) 一= 一 - z 1 - kk 、。 s i n f c = 甚 ( 2 2 ) 这时透射波就以圪速度沿界面滑行，并在第一介质中产生折射波。 2 4 互换原理 互换原理是指地震波在任意介质中传播时，激发点和接收点位置互换，同一种地震 波的传播路径相同。 2 5 直达波，折射波和回折波的传播规律 2 5 1 直达波的传播规律及特征 直达波：由震源直接传播到接收点的地震波，其传播规律如图2 1 。 激发地面 心 。 i 首迸 麓f浪? 矫懿 八 l i 分界面 滑行波 _ 口“：- 二、 2 5 2 折射波的传播规律及特征 图2 1 直达波传播示意图 地震波在传播中遇到下层介质速度大于上层介质速度的界面时，当入射角等于临界 角，透射角等于9 0 0 ，透过波将沿界面滑行，产生滑行波，同时产生滑行波时的入射角 称为临界角。根据斯奈尔( s n e l l ) 定理临界角f c 6 长安大学硕士学位论文 t ：s i 。n l 旦( 2 一3 ) 02 。 l j ) 吃 滑行波引起地层质点振动并传回地面，产生折射波，如图2 2 所示。折射波折向地 面的方向可由斯奈尔定理确定，有 s i n9 0 0 ：s i n r ( 2 4 ) ，2v 1 因此 s i n r = y a l = s i n c n r s l ns =2 哆 ( 2 5 ) 即折射角r 等于临界角i c 。折射波是一系列平行线，其波前是一系列直线。图2 3 表示出了多层介质中的折射波。 om s 图2 2 滑行波和折射波示意图 0m 1m 2m 3 图2 3 水平层状的折射波示意图 2 5 3 回折波的传播规律及特征 回折波时指在连续介质中产生的初至波，有时也叫做潜波，在地震勘探中较为常见， 这是因为在地球的浅部，速度虽深度多为连续增加的。在这种介质中，地震波不在是直 线，而是呈弧形逐步向上弯曲，传播到地下某一深度后开始返回，这种波就是地震勘探 7 第二章地震波传播理论基础 中的回折波。 地震回折波，也被称为连续介质中的“直达波”，当速度随深度线性增加时，地震 波射线是圆弧。地面观测时，会接收到一种和直达波相似的地震波，从震源出发，不经 过界面，直接到达接收点的。但是，它和直达波不同的地方是传播到某一深度以后，沿 圆弧状传播回地面观测系统的，所以，这种“直达波”又称为回折波，如图2 4 所示。 激发点地面 分界回 图2 4 回折波传播规律示意图 2 5 4 初至波与层状介质的对应关系 初至波是最先接收到的地震波，在均匀层状介质中，如果下层介质速度大于上层介 质的速度，那么在震源点附近，首先接收到直达波，如图2 5 中的o m l 段所示。 一段 距离后，将接收到在第一层面界的滑行波产生的折射波，如图2 5 中的m m ：段所示。 经过一段距离以后，会接收到来自第二层介质界面的折射波在第一层中产生的透射波如 图2 5 中的m ，s 段所示。层状介质中初至波的时距曲线如图2 6 所示 长安大学硕士学位论文 0 m im 2 滑行波 p 3 v 3 介质3 图2 5 层状介质初至波传播规律示意图 图2 6 层状介质初至波的时距曲线 第三章初至波射线追踪正演方法研究 第三章初至波射线追踪正演方法研究 很多领域中的问题都可以用网络中的最短路径理论来解决，城市中的交通问题就是 一个典型的例子，从一个地方到另一个地方，有很多条道路，由于车流量和道路的不同， 会产生不同的速度，我们可以用网络中的最短路径原理求取从一个地方到另外一个地方 所需时间最短的路径。在地震勘探中，根据f e r m a t 原理，地震波按最小旅行时路径传 播，我们可以把地下介质离散成一系列网格，每个网格上的节点用地下介质的速度表示， 节点与节点之间的连线用射线路径来表示，这样就可以借鉴网络理论中的最短路径来求 取震源点到接收点的旅行时最小的路径。解决最短路径问题适合用d i j k s t r a 算法，主要 特点是以起始点为中心向外层层扩展，直到扩展到终点为止。d i j k s t r a 算法能得出最短 路径的最优解，但由于其遍历的计算节点很多，所以效率很低。 最短路径射线追踪方法基于f e r m a t 原理和网络理论中的d ij k s t r a 算法，首先，把 地下介质离散成一系列网格，并在网格边界上设置一些节点，速度模型就可以用网格节 点上的速度值表示，网格中的每一个节点只能与相邻的节点连接。从震源点开始，波首 先到达震源点相邻的节点，从相邻节点中选出走时最小的节点，作为新的子波源，再向 相邻节点发出射线，进一步产生新的子波源，重复以上过程，直至获得所有节点的走时 和射线路径。 3 1 速度模型和网格的建立 本文采用的最短路径射线追踪中，把介质离散成一系列相同大小的矩形网格单元， 矩形网格大小可以任意设置，节点设置在矩形单元的四个角点上，用节点上的速度值代 表整个介质的速度，相邻节点之间都可以连接，相邻节点之间的旅行时等于它们之间的 距离与平均慢度之积。有效半径代表地震波的传播方向，有效半径越大，地震波传播方 向越多，计算相邻节点考虑的范围越大射线路径越接近实际地震波的传播特点。图 3 1 ( a ) 所示( 二维) ，一个节点与其相邻的8 个节点相连，这样就有8 个计算节点和8 个 与之对于那个的出射方向。3 1 ( b ) 所示( 二维) ，网格大小为l m * l m ，搜索半径为2 ，一 个节点与其相邻的1 6 个节点相连，这样就有1 6 个计算节点和1 6 个与之对于那个的出 射方向。节点设置在网格节点上，便于计算射线路径在每个网格单元内的大小，更重要 的是便于设置各种复杂模型，引入速度界面，进行地震波射线追踪。 1 0 长安大学硕士学位论文 ( a ) r = i 图3 1 模型网格结构示意图 3 2 节点最小旅行时的计算和射线路径的确定 l? 心钐 勿r。 ( b ) r - 2 从震源点开始，用d i j k s t r a 算法逐点找出从震源到各个网格节点的最小走时。在整 个搜索过程中，网络节点的集合n 分成四个子集p ：已经获得最小走时的节点，即已做 过子波源的节点的集合；q ：已计算出从p 中至少一个节点波传来的旅行时，但还没有 做子震源的节点集合；r ：在n 中除去p 、q 后剩余节点的集合；f ( i ) ：与震源或子震 源直接相连的节点集合，每一步仅计算这些节点的旅行时。具体算法为； ( 1 ) 初始化：p = m ，t ( s ) = 0 ，s 为震源，表示空集；= n ，t ( i ) = ，f n 。 ( 2 ) 选择：在q 中选择旅行时最小的节点i ，i q 。 ( 3 ) 更替：计算从i 点传到j 点的旅行时，若该值比原值小，则用该值取代原值，否 则保持原值不变。即，对f ( i ) n q ，( ) = m i n t ( j ) ，t ( i ) + d i j 】；对 j f ( i ) n r ，t ( j ) = f ( f ) + d i j 。将节点i 从q 转至p 。 ( 4 ) 迭代判断：如果p = n 或q = m 停止迭代，否则转向( 2 ) 确定震源点到接收点的波传播射线路径，是从接收点向震源点逆向搜索的过程。依 次找出前一级子波震源，直到搜索到震源点为止，把震源点和找出的子波震源依次连接 起来，得到震源点到接收点的射线路径。 在迭代过程中，最耗时的一步是每一步都要从q 中找出旅行时最小的点作为下一 个子波震源，包括子波震源最小走时的选取，子波震源的更新以及新子波震源的添加等。 本文采用快速排序法，对子波震源进行管理。 第三章初至波射线追踪正演方法研究 3 3 存在问题及一些改进方法 最短路径射线追踪算法，用网络节点之间的连线近似地震射线路径，存在如下一些 问题：当网络节点较少时，得到的地震射线路径呈之字形，计算走时比实际走时偏大； 需要大量的计算时间和计算机内存，特别是三维问题；在速度平缓区域，在两个节点之 间常常会有几条射线路径，不能从中选取一条接近实际射线的路径，在计算过程中，需 要增加一些约束条件，进行合理的选取。 针对上述存在的问题，一些学者提出了改进的方法：k l i m e s 和k v a s n i c h a 及z h a n g 和t o k s o z 等通过改变节点分布，m o s e r 等通过用其他方法进一步修正射线路径，以减 小计算误差；张建中等人引入了走时插值算法思想，提出了一种动态网络最短路径射线 追踪方法，消除了波前只能从一个节点到达另一个节点的局限，提高最小走时和射线路 径的计算精度。王辉等( 2 0 0 0 ) 对其在节点走时的计算中引入b r e s e n h a m 画线算法，在 走时节点查询中使用快速排序法和插入排序法，显著的提高了算法的精度和速度。 1 2 长安大学硕士学位论文 第四章地震初至波层析成像反演方法 照相及录像是利用光敏材料，把从物体表面反射的光波记录下来，这种成像方法属 于表面成像。另一种是在不切开物体或人体的情况下对其内部结构成像，这种成像的方 式又可分为全息成像和层析成像，全息成像技术主要用于对物体内部目标的轮廓成像， 而层析成像时对物体进行逐层剖析成像，相当于把物体切成片，然后照相。层析成像是 成像技术发展到目前的最高阶段，所谓c t 技术乃是计算机辅助层析成像的简称。 不“打开物体，而了解它的内部结构，常用的办法是用波穿透物体，让波带出关 于物体内部的信息。从物体内部图像重建的角度看，一张物体切片的图像两个空间变量 ( x ，y ) 的函数，称为图像函数，记为f ( x ，y ) ，用不同方向的入射波“照相”物体，观 测到的波场信息至少是入射波方向日和观测点位置p 两个变量的函数，称为投影函数 可以恢复唯一的图像函数f ( x ，y ) 一拉当变换，这个定理后来成为层析成像的理论基础。 总之，层析成像不仅利用波作为探测手段，而且要同时利用各种不同方向投射的波场信 息。 地震波层析成像可按输入资料的不同分为三大类：( 1 ) 输入地震波走时资料的成像 方法，对介质或波速成像；( 2 ) 输入某种波的最大振幅或某个波组的平均振幅，对介质 散射系数成像；( 3 ) 输入一定时窗内的接收地震波场和背景波速对介质波速扰动成像。 这三类方法可分别简称为走时反演成像，振幅反演成像及波场反演成像。按地震波理论 可分为射线层析成像和波动方程层析成像，波动方程层析成像能利用走时、振幅、相位 和频率等全波形信息、能提高分辨率，但是波动方程存在一些问题，没法消除各种干扰 因素的影响；射线层析成像虽然只利用了地震波旅行时信息，但干扰因素少，采用误差 较小的反演方法和先验信息，可以获得满意的反演效果。 4 1 模型的建立 针对本文采用最短路径射线追踪方法节点设置的特点，对于近地表模型，把介质离 散成一系列矩型小单元，给每个单元节点设置速度参数。如图4 1 所示，模型被离散成 ( m 一1 ) 木( 一1 ) 个矩形单元。单元节点排序按从左到右，从上到下顺序编号，射线路径 为矩形单元边或对角线，每个单元的射线长度等于矩形单元的边长或者对角线长度，这 样每个单元里边的射线长度很容易计算，当两个节点之间的连线速度突变时，射线的速 第四章地震初至波层析成像反演方法 度值为两个节点速度值的平均值。两个相邻节点之间的旅行时等于它们之间的距离与平 均慢度之积。 2 1 n z 上n 图4 1二维近地表模型离散化 4 2 射线路径和初至时间的计算 r 对于图4 1 的模型，用最短路径射线追踪方法进行射线追踪，给定激发点和接收 点坐标，以及地下介质模型大小和速度值，对模型介质离散化，确定网络大小和搜索半 径，然后进行正演计算，得到各个节点上的射线路径和走时，以及接收点的走时信息。 在反演过程中，利用射线追踪得到的节点上的射线路径长度和接收点走时，利用s i r t 层析成像的反演方法，得到反演结果。 4 3 反演方程的建立 根据地震射线理论，初至波旅行时与近地表速度有下列关系： f = f 知 ( 4 1 ) 月品) v ( x ，z ) 、。 其中，r ( v ) 为正演射线路径，速度v ( x ，z ) 为射线经过节点上的速度的平均。速 度用慢度可以表示为s ( x ，z ) = 1 v ( x ，z ) ，式4 1 可表示为： f = is ( x ，z ) d l ( 4 2 ) 尺i ) 假设第j 个节点的慢度是s ，第f 条射线的旅行时为，第i 条射线通过与其相连的 l 1 2 i j 长安大学硕士学位论文 两个节点的长度为射线长度乇，因此可以把式( 4 2 ) 可写成： 乃0 = t ， ( f = 1 ，2 ，聊) ( 4 3 ) j = t 式中，m 为正演的射线条数，n 为设计模型节点个数。用矩形( 或向量) 表示为： 其中，s = s l 是 a s = f ( 4 4 ) 彳= 厶。：，l 。 如，乞： ； 乞。 乙。乞：乙。 式( 4 4 ) 看似线性关系，由于慢度向量和射线路径都未知，射线路径又是慢度的非线 性函数，因此，该式为非线性关系。但慢度矢量s 和射线路径矩阵a 都不知道，而且射 线路径又是慢度的非线性函数，因此，该式实际上为非线性方程。 4 4 层析成像的算法 地震波走时层析成像方法非为线性方法和非线性方法，线性方法有反投影技术 ( b p t ) ，代数重建技术( a r t ) ，联合迭代重建技术( s i r t ) ，共轭梯度最小二乘法( c g l s ) ， 最小二乘正交分解法( l s q r ) 等，非线性反演的方法主要有速传算法、模拟退火法和 神经网络法等。 b p t 、a r t 和s i r t 是上世纪7 0 年代发展起来的三种线性反演算法，主要是由于当 时的计算机内存小设计的。b p t 方法简单，概念明确，但分辨率较低，只适用于粗略扫 描。a r t 方法计算速度快，但收敛性性不太好。s i r t 收敛性好，计算速度慢一些。 非线性反演方法是非线性问题线性化，主要问题是不稳定。非线性反演的方法为全 局搜索法，反演不依赖与初始模型的选取，当反演区域初始信息少时，可以为线性反演 提供初始信息，虽然计算速度慢，但对于复杂的非线性问题反演效果好。下面介绍两种 线性反演方法。 1 a r t 是a l g e b r a i cr e c o n s t r u c t i o nt e c h n i q u e 的简称，翻译为代数重建技术。 a r t 反演是个迭代过程。迭代过程为：首先，给定初始慢度向量s ? ( 产1 ，2 ，) ， 然后，计算第i 个方程( 即第i 条射线) 的估计值，用观测值和估计值之差计算慢度修正值， 1 5 第四章地震初至波层析成像反演方法 按方程组4 4 从第一个方程到最后一个方程依次对慢度向量进行修正，知道修正后的慢 度向量满足迭代次数要求或误差要求为止。 若用k 表示对慢度向量的修正次数，用f ( 江1 ，2 ，f ) 表示方程组( 4 4 ) 中的方程 的序号，则k 与i 的关系为 f - m o d ( k ，f ) ( 4 5 ) 若用s r ( ，= 1 ，2 ，) 表示第( k 一1 ) 次修正后的慢度向量，厶；( = 1 ，2 ，j ) 表示第 k 次慢度修正量，则第k 次慢度修正量应满足如下方程 ， 岛- z ( 砖卜1 + q ) = o ( 4 6 ) 其中，f 是由( 4 5 ) 式确定1 的。满足次方程的q ( 歹= l ，2 ，) 有无穷多，但须保 证s i f ，= 1 , 2 ，j 1 收敛，所以通常取 血，) = ， ， ( 6 ，- e s 扛1 ) j 口2 ， j = l ：# ( 1 ，2 ， = 1 ， 其中，= ( 包- z c l 9 8 ) k 一1 ) j = l ( 4 7 ) ( 4 8 ) 似是第( k 1 ) 次修正后的慢度向量对第k 次修正对应的第i 个方程的误差。此时， 第k 次修正后的慢度向量为 ：s ) k - + a s ) t ) ：s ( k - 1 ) 下丁c y r 一- ( k ) ( _ ，：1 ，2 ，力 2 j = l 综上所述，a r t 方法的具体实现步骤为： ( 1 ) 给定一组初始值拶( ，= 1 ，2 ，) ； ， ( 2 ) 计算第i 个方程( 即第i 条射线) 的估计值) 够- 1 j = l ( 3 ) 计算第i 个方程( 即第i 条射线) 的观测值与第i 个方程的估计值之差 1 6 ( 4 9 ) ( 4 1 0 ) ( 4 1 1 ) 长安大学硕士学位论文 ( 4 ) 计算第i 个方程的慢度修正量 q 。= 乃，：似 |j = l ( 4 1 2 ) ( 4 1 3 ) ( 5 ) 计算修正后的慢度值 够) - s 严d + 斛”= s y 。+ 乃似2 ( 4 1 4 ) ，1 = 1 并且根据其它资料对s s u = 1 ，2 ，j ) 加以约束： s k ) s 一，即 若。，令= s m i n 蔫s j 2 s 嘁，令s j 2 s 啦 ( 6 ) 将第一个方程的修正后的的代入第二个方程，依次进行下去，直至完成第一 轮。 ( 7 ) 笋l j 断迭代结果满足预定误差的要求时，则停止计算，不满足则重复以上步骤， 进入下一轮得循环，直到满足预定残差要求为止。 2 s i r t 是s i m u l t a n n e o u si t e r a t i v er e c o n s t r u c t i o n 的简称，翻译为联合迭代重建法， 最早由g i l b e r t ( 1 9 7 2 ) 提出。与a r t 的方法不同的地方是，s i r t 不是逐个方程( 或逐 条射线) 逐次地对慢度向量s i ( j = 1 ，2 ，) 进行修正，而是在每次迭代中所有方程同时 参入对慢度向量s j ( j = l ，2 ，j ) 的修正。在给定慢度向量的初值s ：( j 2 1 ，2 ，) 后，或 在第( k 1 ) 次迭代求出了修正的慢度向量s 他。1 ( j = 1 ，2 ，) 后，如( 4 7 ) 式所给出的， 为使单个的第i ( i = l ，2 ，i ) 个方程成立，第i ( i - 1 ，2 ，i ) 个方程对应的慢度修正量 l l 为s y = 巳，：2 ( ，= 1 ，2 ，) 。s i r t 方法是所有方程对同一个节点作平均后， 作为当前慢度的修正值，这样就可以压制一些干扰的因素，s i r t 的计算结果与观测的 数据使用顺序无关。 第k 次迭代的慢度修正量为 衅= 击喜c 叫) 姜叱 - 5 ， 、， 一 似， s 矿 口 j 一 一 6 1 l = r 1 第四章地震初至波层析成像反演方法 其中，m j 表示方程组( 4 4 ) 的系数矩阵a 的第j 列中非零元素的个数，第k 次代 修正后的慢度为 够) - 1 蟛”= 够1 ) + 古妻c 喜，( _ = 2 ，刀 ( 4 舶) 为加快收敛速度，可在( 4 1 6 ) 式中加入松弛因子，s i r t 方法的递推公式可写为 = 护+ 暑喜百a o r , ( k ) 俨墟，力o r 2 , 0 a 2 屯= 1 口；( 歹：1 ，2 ，j ) ，a ：ja 。2 m ( f _ l ，2 ，) 综上所述，s i r t 的具体实现步骤为： ( 1 ) 给出一组初始值巧= ( = 1 ，2 ，刀 j ( 2 ) 计算估计值 + s 5 。1 ( 江1 ，2 ，) j = l ( 3 ) 计算观测值匆“= 1 ，2 ，) 与估计值的差( ，= 1 ，2 ，) ( 4 17 ) ( 4 1 8 ) ( 4 1 9 ) “2 0 ) ( 4 ) 计算第j 个象素内的平均修正值厶于，假设在象素内有鸠条射线通过，则： s 尹= 击喜c ) 喜，c 歹= a 一，d c 4 m ， ( 5 ) 对第j 个象素的慢度值进行修正够= 砖扣1 + 衅，并且根据声波测井或者其 它地球物理资料对够( = 1 ，2 ，) 加以约束：够- - s m a x ，即 蒂s j s 幽，令s j 2 s 豳 若s j s 嗽，令s j = s 。 ( 6 ) 判断迭代结果是否满足预定误差的要求，满足停止，不满足则重复以上步骤， 进入下一轮迭代。 1 8 长安大学硕士学位论文 第五章模型计算 5 1 最短路径射线追踪计算值与真值对比 设计模型1 为均匀介质，参数设置：模型大小1 5 0 m * 1 0 0 m ，炮点( 0 ，o ) ，接收点位 置( o ，1 0 ) ( 0 , 2 0 ) ( 1 0 0 ，o ) ，网格大小为l m * l m ，搜索半径为1 0 和2 0 ，速度5 0 0 m s 。图 5 1 为搜索半径为1 0 的射线路径，图5 2 为搜索半径为2 0 的射线路径，表5 1 为最短路 径射线追踪计算值与真值的比较，从表中可以看出搜索半径增大，误差减小，精度提高， 最短路径射线追踪方法是一种有效、精确的射线追踪方法。 图5 1 射线路径，网格大小为1 1 ，搜索半径为1 0 1 9 第五章模型