（地球探测与信息技术专业论文）折射波层析成像在西安地裂缝勘探中的应用研究.pdf

摘要 地震层析成像，是在计算机层析成像( c o m p u t e r i z e dt o m o g r a p h y ，简称c t ) 基础 上发展起来的一种新方法，是调查和认知地下地质结构的有效手段，目前已广泛应用于 工程与地质灾害等方面的地球物理勘探领域。地震层析成像是利用地震数据来反演地下 结构的物质属性，可确定地球内部的精细结构和局部的不均匀性，成像过程包括模型参 数化、正演计算、反演计算，其中反演是地震层析成像的关键部分。本文在广泛查阅国 内外文献的基础上，介绍了地震层析成像技术的正反演理论基础，并对地震层析成像三 种反演算法进行了对比。文中采用l s q r 算法，对西安地裂缝理论地质模型进行反演模 拟，确定了该方法的正确可行性。文章最后，将折射地震层析成像技术应用于西安地裂 缝探测中，取得了成功。 论文绪论部分叙述了c t 技术的发展历史，简单介绍了地震层析成像的基本概念、 分类及目前的研究现况；第二章，介绍了地震波传播的基本原理，对折射波走时曲线方 程进行了论述；第三章，叙述了地震层析成像反演的理论基础，通过对理论模型正反演 的结果，对三种反演方法进行了对比分析；第四章介绍了射线追踪的方法及理论基础， 针对西安地裂缝的实际地质状况设计了理论模型，并进行了正反演计算；文章最后，对 实际采集的西安地裂缝勘探地震资料进行了折射地震层析成像处理，得到了近地表地层 速度场，地裂缝位置处速度异常明显，处理结果与实际吻合。 本文首次将折射地震层析成像应用于西安地裂缝勘探中，通过上述研究工作，得出 以下研究成果： 1 通过对a r t 、s i r t 、l s q r 算法反演比较，证明l s q r 算法具有反演速度快、 稳定性好、反演精度高等优点； 2 折射地震层析成像技术可成功探测西安地裂缝，在其速度剖面上能判断地裂缝的 位置、倾向、倾角等参数； 3 在不增加传统的反射波地震勘探外业工作量基础上，利用折射波信息，可给出近 地表几十米内的地层剖面速度场，刚好弥补了西安地裂缝中反射波勘探在近地表的空 白。 关键词：西安地裂缝，折射波，层析成像，反演 a b s t r a c t t h es e i s m i ct o m o g r a p h yi san e wm e t h o df o u n d e do nt h ec o m p u t e r i z e dt o m o g r a p h y ( c t )a n di se f f e c t i v ef o ri n v e s t i n ga n du n d e r s t a n d i n gt h eg e o l o g i cf r a m e w o r k i ti s e x t e n s i v e l ya p p l i e dt ot h ef i e l do fg e o l o g yd i s a s t e ra n dg e o - e x p l o r a t i o na n ds oo n i ti sa t e c h n i q u ew h i c hu s e ss e i s m i cd a t et oi n v e r tt h em a t e r i a la t t r i b u t eo fu n d e r g r o u n ds t r u c t u r ea n d d e t e r m i n et h ef i n es t r u c t u r eo ft h ee a r t hi n t e r i o ra n dr e g i o n a li n - h o m o g e n e i t y i ti n c l u d e s m o d e lp a r a m e t e r i z a t i o n , f o r w a r dm o d e l i n gc a l c u l a t i o n ，a n di n v e r tc a l c u l a t i o n ，i t sk e yi s i n v e r t i o nt e c h n o l o g y o nt h eb a s i so fe x t e s i v e l yr e f e r r i n gt od o m e s t i ca n df o r e i g nd o c u m e n t s ， t h i sa r t i c l ei n t r o d u c et h ef o r w a r da n di n v e r t i o nc a l c u l a t i o nt h e o r yo fs e i s m i ct o m o g r a p h ya n d c o n t r a s tt h et h r e ei n v e r t i o na l g o r i t h m t h ea r t i c l ei n t r o d u c e dl s q rw h i c hs i m u l a t et o i n v e r t i o nx i a ng r o u n df i s s u r et h e o r e t i cg e o l o g ym o d e la n dm a k es u r et h i sm e t h o d s p r e c i s e n e s sa n df e a s i b i l i t y i nt h ee n d o ft h ep a p e r , i ti ss u c c e s s f u lt oa p p l yt h et e c h n o l o g yo f t h er e f r a c t i o ns e i s i m i ct o m o g r a p h yt ot h ee x p l o r a t i o no fx i a ng r o u n df i s s u r e i nt h i sp a p e r ，t h em a i nc o n t e n t so fc h a p t e r sa sf o l l o w s ： t h ep a r to fp r o l e g o m e n o nd i s c u s s e dt h ed e v e l o p m e n to fc tt e c h n o l o g y , a n di n t r o d u c e d t h eb a s i cc o n c e p t i o no ft h ec t , i t ss o r ta n dt h ep r e s e n tr e s e a r c h c h a p t e r2i n t r o d u c e dt h e b a s i cp r i n c i p l eo ft h es e i s m i cw a v ea n dd i s c u s s e dt h ee q u a t i o nt ot h et r e n dc u r v eo fr e f r a c t e d w a v e c h a p t e r3i n t r o d u c e dt h eb a s i ci n v e r t i o nm o d e lt h e o r yo fc t ，b yt h ef o r w a r da n d i n v e r t i o nt h e o r ym o d l eo fc t , a n dc o m p a r e dt h et h r e ei n v e r t i o nc a l c u l a t i o n c h a p t e r4 i n t r o d u c e dt h em e t h o do fm a k i n gt r a c k sf o rr a d i a la n di t sb a s i ct h e o r y , d e s i g n e dt h et h e o r y m o d e l ，a c c o r d i n gt ot h ep r a c t i c a lg e o l o g i c a ls i t u a t i o no fx i a ng r o u n df i s s u r et h e nm a d ea f o r w a r da n di n v e r t i o nc a l c u l a t i o n a tl a s ti td e a lt h ep r a c t i c a le x p l o r i n gs e i s m i cd a t eo ft h e x i a ng r o u n df i s s u r ew i t hc t ，w h i c hg e n e r a l i z ev e l o c i t yf i e l da r o u n dr a p i d i t ys e c t i o np l a n e t h ev e l o c i t yw i t h i ng r o u n df i s s u r ei s a p p a r e n t l yd i f f e r e n t t h er e s u l tc o i n c i d e dw i t hr e a l s i t u a t i o n c tf i r s t l ya p p l yt ot h ep h y s i c a lp r e s p e c t i o no fx i a ng r o u n df i s s u r e b yt h ea b o v e r e s e a r c h ，t h ea c h i e v e m e n ta sf o l l o w ： 1a f t e rc o m p a r a t i o no ft h ea r t , s i r t , l s q r , i ti sp r o v e dt h a tl s q rm e t h o di sg o o da t i n v e r s i o ns p e e d ，s t a b i l i t ya n dp r e c i s i o n 2i t i ss u c c e s s f u lt op r o bt h ex i a ng r o u n df i s s u r e ，a c c u r a t l yd e t e r m i n et h ep o s i t i o n ， t e n d e n c y , a z i m u t ho ft h eg r o u n df i s s u r ei nr e f r a c t i o ns e i s m i ct o m o g r a p h yp r o f i l e 3o nt h eb a i s eo fn o ti n c r e a s i n gt h ec o n v e n t i o n a lr e f l e c t e dw a v ew o r k l o a df r o ms i s m i c e x p l o r i n gt ou s er e f r a c t i v ew a v ei n f o r m a t i o n ，w h i c hc a nd e t e r m i n ev e l o c i t yf i e l do fs t r a t u m s e c i o n i tj u s to f f s e tv a c a n c yo fx i a ng r o u n df i s s u r er e f r a c t i v ew a v es e i s m i cs u r v e y k e yw o r d ：x i a ng r o u n df i s s u r e ；r e f r a c t e dw a v e ；s e i s m i ct o m o g r a p h y ；i n v e r t i o n i i 论文独创性声明 本人声明：本人所呈交的学位论文是在导师的指导下，独立进行研究 工作所取得的成果。除论文中已经注明引用的内容外，对论文的研究做出 重要贡献的个人和集体，均已在文中以明确方式标明。本论文中不包含任 何未加明确注明的其他个人或集体已经公开发表的成果。 本声明的法律责任由本人承担。 彦怍 训7 年6 月l 日 论文知识产权权属声明 本人在导师指导下所完成的论文及相关的职务作品，知识产权归属学 校。学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权 利。本人离校后发表或使用学位论文或与该论文直接相关的学术论文或成 果时，署名单位仍然为长安大学。 ( 保密的论文在解密后应遵守此规定) 论文作者签名： 1 1 垮 q 导师签名： 叼引胪日 r 肌日 长安大学硕士学位论文 1 1ct 成像技术的背景 第一章绪论 c o m p u t e r i z e dt o m o g r a p h y 简称c t ，即计算机层析成像，是由投影重建图像的应用 技术之一。层析成像的英文( t o m o g r a p h y ) 源于希腊语t o m o s ，本意是断面或切片的意思。 c t 技术是根据物体外部所测量的数据，按照一定的物理关系和数学关系反演物体内部 物理量的分布，最后得到清晰分布图像的技术 6 1 。 c t 技术最早应用于医学领域之中，之后迅速被引入到各个科学领域。如在天文学、 核磁共振、显微技术、无损检测、磁流体力学、物质结构、化学和地球物理等方面。其 数学基础是1 9 1 7 年奥地利数学家j r a d o n 发表的“关于函数沿某些回路积分确定该函 数”的论文，它的物理思想是1 9 2 1 年法国医生a b mb o c a g e 提出的焦平面层析技术 ( f o c a lp l a n et o m o g r a p h y ) 。1 9 7 2 年第一台x 射线c t 装置诞生了，使得c t 技术在医学上 的应用逐渐成熟。为此，1 9 7 9 年的诺贝尔生理和医学奖授予了c t 的发明者物理学家 a m c o r m a c k 和工程师gn h o u n s f i l e d 。层析成像技术的意义在于通过在物体外部的 非破坏性探测，能获得物体内部物性分布的图像，就好像在某种领域“光线”照射下，物 体成为透明体，使物体内部信息得到充分显示。 而将层析成像技术应用到地学领域，对探查地下异常体的空间展布起到很大的作 用。通过探查工程岩体或矿层中的异常体，可有效地指导工程建设和采矿生产的快速、 高效运行，避免重大人员和财产损失。 1 2 地震层析成像技术的基本概念和分类 上世纪8 0 年代，层析成像技术进入到地学研究领域，被称为地学层析成像 ( g e o t o m o g r a p h y ) 。它是利用医学c t 的理论对地球内部进行探测方法的总称，是一种详 细调查地下物性参数分布状况的地球物理探测技术。根据传播物质的不同，地学层析成 像的研究内容又可以分为地震层析成像、电磁波层析成像和电阻率层析成像。其中，地 震层析成像又是地学层析成像中发展最成熟、应用最广和研究最深入的技术。 地震层析成像技术( 下称地震c t ) 是由布置在孔内的激发点和接受点所组成的观 测系统，由激发点产生地震波，经地下地质体的折射、反射等物理过程到达接收点，并 用接收点所记录的直达波、反射波以及面波等震相的走时或振幅资料，在计算机上通过 第一章绪论 不同的数学处理方法得出地震波速度场分布图像，并通过速度的分布特征以及速度大小 重建各种地质异常体的位置和分布情况。 地震c t 如果从数据的观测形式来分，可分为地面c t 、井间c t 、v s pc t 和联合 c t ；若从成像所用的震相分类，可分为反射c t 、透射c t 和绕射c t ；根据地震波传播 理论，地震层析又可分为两种：一种是射线c t ，它忽略地震波动力学特征，是在射线 路径上将旅行时反投影；另一种是波动方程c t ，它是在波程上把微分波场反投影。 1 3 地震层析成像方法在勘探中的进展 最早进行地震层析成像是b o i s 等人( 1 9 7 1 、1 9 7 2 ) ，利用地震波动旅行时推测两井 问的主要构造边界。m a s o n ( 1 9 8 1 年) 在煤层上利用首波做了代数重建实验。在地震勘探 中，真正开展层析成像技术的研究工作是从上世纪8 0 年代开始，1 9 8 3 至1 9 8 5 年，美国 地球物理联合会( a g u ) 和美国勘探地球物理协会( s e g ) 多次组织c t 专题研究，这是地学 界c t 研究广泛展开的一个标志。 随后，地震层析成像正、反演理论及应用研究取得了丰硕成果。在正演算法方面， v i d a l e ，周兵等，黄联捷、李幼铭等，l a n g g a n 等提出了快速、精确的计算走时和射线 路径的射线路径的射线追踪法，p a r t 提出了一种适合于地震层析成像的频率域中的弹性 波正演模拟的有限差分法，c a v i c c h i 提出了一种时间域内波动方程正演算法。在反演算 法上，s c a l e 提出了共扼梯度法和快速l p 方法，d v a n e y 及w u 等( 1 9 8 7 ) 提出了用于波动 层析成像的滤波反传算法，d e v a n e y 等( 1 9 9 1 ) 研究了衍射层析成像广义a r t 法，b e y l k i n 和栾文贵研究了广义r a d o n 变换进行地震散射波场的层析成像的方法技术，王宏伟提出 菲涅耳特地震层析成像技术，l i u 研究了将脉冲谱技术( p s t ) 应用于二维波动方程反演 中，李世雄等研究了广义r a d o n 变换在散射层析成像法中应用，吴律等对有限观测角层 析进行了研究，张霖斌等研究了层状介质中衍射层析成像算法，p h i l i p s 等比较了旅行时 层析成像中的常用算法，l u o 等提出了波动方程旅行时的层析成像方法，袁晓晖研究了 密度与压缩系数的散射层析成像。在利用反射波资料方面，黄鑫、吴建成、e l a n d a ，c h j r s w a r d 、左海道等在反演方面也有所不同。近年来，又出现了仿真淬火法、遗传算法 和神经网络法，在重建图像时收到很好的效果。这些理论研究在数值模拟中取得了令人 满意的结果，为应用打下了良好的基础。 地震c t 成像在地球物理领域的应用非常广泛。a n d e r s o n 等公布了他们对地球三 维速度结构成像的结果，用地震c t 技术揭示地球内部精细结构的层析图，日本物理学 2 长安大学硕士学位论文 家柳谷俊等利用c t 技术研究了岩石破裂局部化，岩石断层形成过程以及浸水岩石的波 速与吸收的变化情况，n e l s a n 等用旅行时层析成像确定地震位置，h a r r i s 等利用井间层 析成像检测井间的断裂构造和孔岩层，宋正宗利用井间层析成像在二滩电站对岩体质量 进行分级，黄金莉等利用地震透射波走时层析成像技术对地下巷道揭露的断层点的连通 情况进行了探测，获得较好的效果。 随着高度集成化多道数字地震仪广泛投入使用，地震c t 技术进入了地球物理勘探 的实用阶段，在许多大型构造探查中，如洋中脊张裂断层、板块隐伏挤压带一级油气藏 构造等，地震c t 技术都发挥出它的优越性。 1 4 本文主要研究目的及意义 西安地裂缝是西安特有的一种地质灾害，目前在西安市区内已发现具有一定长度规 模的地裂缝( 带) 达1 4 条之多，地表出露总长度超过1 2 0 k m 5 3 】。在西安地裂缝的地 震勘探中，主要采用反射地震方法，由于反射波地震方法很难做到勘探目的层深浅兼顾， 其勘探目标通常为几十米到几百米，其勘探目的是对该深度范围内的地裂缝性质和特征 进行确定。这样一来，经资料处理后得到的反射时间剖面上，近地表几十米范围内一般 为空白。而折射波，作为反射波资料处理中通常被切除掉的干扰波，却携带了丰富的地 下波场信息，特别是近地表的信息。在反射波勘探中，偏移距多在3 0 8 0 m ，这时得到 的现场记录中折射波非常突出清晰，其反映的地层信息恰好在近地表数十米范围内。 在西安地裂缝的变形带内，由于地裂缝附近的土体结构遭到破坏，造成结构强度下 降，已有的研究表明，地裂缝变形带内地层剪切波速及等效剪切波速明显低于缝两侧正 常值，这种岩土体强度降低的影响范围约在1 0 m 左右【l 引。 因此本文拟在不增加地裂缝勘探外业工作量的基础上，利用折射波震相信息，反演 近地表地层，根据反演得到的速度场信息，结合西安地裂缝变形带内岩土体波速降低这 一性质，得到地裂缝及其变形带内的低速异常，结合反射勘探结果，给出地裂缝从地表 到深部数百米范围内的位置、倾向、倾角等参数。 如果利用折射层析成像方法成功实现上述目的，则将对西安地裂缝研究和工程勘察 具有直接的现实意义，主要体现在如下几方面。 1 地震反射波勘探在西安地裂缝勘察中，主要从地层深部定性判断地裂缝，遗漏近 地表详细资料，而折射地震层析成像可反映出近地表地层速度差异，弥补了反射波勘探 第一章绪论 在近地表的空白区，折射波层析成像与反射波勘探的有效结合，可为地裂缝研究和工程 勘探提供更多的对比判断信息。 2 折射地震层析成像结果可以辅助或部分取代近地表的的钻探、槽探，这将大大节 省研究和勘察经费，降低建设费用。同时，由于折射层析成像方法无需额外的外业工作 量，可缩短勘察施工周期，为工程项目争取时间。 3 折射层析成像方法如能在西安地裂缝勘探领域中得到成功应用，将为指导实际生 产提供相应的理论依据。 1 5 论文的研究内容 本论文主要进行如下几方面内容研究： 1 对几种地震层析成像方法进行分析，编制相应计算程序，用这些方法对有代表性 的理论模型进行反演，对反演结果进行对比，得到各反演方法的优缺点。 2 根据折射波射线追踪理论，编写计算程序，设计符合西安地裂缝的理论模型，通 过对模型正反演计算分析，确定利用折射波震相进行层析成像方法在西安地裂缝地震勘 探中的可行性。 3 利用上述反演方法和射线追踪理论，对西安地裂缝的实际地震采集资料进行处理 分析，结合已掌握的地质资料，如已知地裂缝的位置、倾向、倾角、走向等，判断该反 演方法的实际效果，验证方法的正确性。 4 长安大学硕士学位论文 第二章地震波传播理论基础 地震波是在地下岩层中传播的弹性波，其传播规律满足波动方程，但由于实际岩层 的非均匀性，使得波动方程的边界条件非常复杂，以至于无法解出该波动方程。在一定 条件下，用射线理论来近似波动理论，可使数学问题大大简化，从而能方便的研究地震 波在地下介质中的传播路径、传播速度、旅行时间等运动学特征【2 0 】。地震波在岩层中的 传播速度与岩层的性质，诸如弹性参数、岩石成分、密度、埋藏深度、地质年代、孔隙 度等因素有关，这导致了地震速度分布规律及波传播规律的复杂性，给研究问题造成了 困难。因此，地震勘探理论一般把地下介质做以下几种近似： 均匀介质：这种介质假设地震波在岩层各处的传播速度都相同。 层状介质：假设岩层是成层分布的，地震波在每一层中的传播速度都是相同的，而 各层之间的地震速度不同。 连续介质：假设地震波在岩层中的传播速度是随空间位置连续变化的。连续介质的 物性随空间坐标连续变化，显然是非均匀介质。层状介质实际上也是一种非均匀介质， 只不过物性的不均匀只出现在分界面处，是一种比较特殊的非均匀介质。 岩层介质结构不同，地震波的传播规律也会不同。对地震波在不同介质中传播特征 的研究，是地震勘探的重要内容，也是进行定性和定量解释的依据。下面首先介绍地震 波传播过程中的一些基本原理和定律。 2 1 惠更斯( h u y g e n s ) 原理2 7 】 h u y g e n s 原理是利用波前的概念来处理问题。地震波向外传播时，是由震源产生的 振动，通过地层介质的质点依次传播出去的。h u y g e n s 原理指出，波前面上的每一点都 可以看作是新的震源，而这些新震源发出的子波波前的包络面，就是新的波前面。根据 这个原理，只要我们知道某一时刻的波前位置，就能够确定出地震波在各种不同时间的 波前位置。当波遇到另一种介质时，这个原理可以形象地说明波的反射、透射、折射和 绕射等现象。 2 2 费马( f e r m a t ) 原理【2 7 】 f e r m a t 原理是从射线角度描述波传播特点的。它指出，波在各种介质中地震波沿射 线传播的时间和沿其它任何路径传播的时间相比为最小，或者说波沿所花时间为最小的 5 第二章地震波传播的理论基础 路径传播。根据这个原理可以确定地震波在己知传播速度的介质中的射线形状，如在均 匀介质中射线路径为直线，而在非均匀介质中，射线路径不再是直线，是与波前面垂直 的曲线。 2 3 斯耐尔( s n e l l ) 定理【2 7 】 根据s n e l l 定理，透射波在第二介质中的偏折，与两种介质的速度比有关。令v 1 代表第一层介质，v 2 代表第二层介质，当v i v 2 时，透射波靠近法线偏折，当v i j、lrj 矿一锄可一钞 rl，。1l 足 足 第三章地震c t 的理论基础及反演算法比较 即得 厂( q ，哆) = f ( c o c o s o ，c o s i n o ) = 妙( 细s 日+ s s i no 舢i n - sc o so ) e - e n d s d t r 2 = ，r f ( t ，p ) p f 耐d t = r f ( c o ，a ) r l 利用r a d o n 变换与f o u r i e r 变换之间的关系式( 3 5 ) 及f o u r i e r 变换的反演公式 厂( y ) 2 西孬1 旷( q ，哆) p f ( q ，+ 喇) d q d 哆 即可求得r a d o n 变换的反演公式 ( 3 5 ) m 川= 南了硼h ! 矽( t , o ) e - i n d t 卜刚州叫砒 6 ， 上面的反演公式还可以进一步化为 弛：一南m s 毒篙 ( 3 7 ) 由此可见，r a d o n 变换的反演公式包括：求导，h i l b e r t 变换及对0 求平均三种运算。 因此，当f ( x ，y ) 满足适当条件时( 这在实际问题中一般总是满足的) ，由r f ( t ，日) 可以唯 一确定厂( x ，y ) ，也就是说r a d o n 变换的反演问题是存在且唯一的，但是要成为数学上 适定的问题，还必须满足稳定性的要求。另外，反演公式( 3 6 ) 或( 3 7 ) 在理论上给出了重 建连续函数厂( x ，y ) 的图像的方法，它要求射线入射角0 从0 到丌连续变化，而且投影函 数r f ( t ，0 ) 的观测完整和准确，而在实际测量中这些假设通常是不能满足的，因此以上 反演公式只有理论上的意义，并不适宜具体数值计算。 在地震c t 反演时，是通过炮点和接收点来实现对反演区域内的信息的获取。设s ( r ) 为射线穿过介质的慢度( 速度的倒数) ，名分别表示炮点、检波器点坐标，f ( ，名) 表 示旅行时，则射线积分可表示为： t ( 5 ，名) = j - j ( r ) d l 川，s l 1 0 ( 3 8 ) 长安大学硕士学位论文 其中，尺 ，乞，s 】表示积分曲线，称为射线。根据f e r m a t 原理，射线不但与炮点、检波 器点的位置有关，而且与介质慢度有关，即积分函数s ( ，) 与积分路径研，名，j ( ，) 】有关。 可见，地震波层析成像是个非线性的反问题。对于这样的反问题，理论上可以用广 义r a d o n 逆变换来直接求解，但在实际地震勘探中，由于观测角度有限，无法提供上述 反算子需要的全部数据。因此，只能采用正问题来解反问题的办法，一般用迭代法求解， 其基本过程为：给定正演时波的旅行时或从实际资料中拾取旅行时，即观测走时，然后 给定一个初始速度模型，用射线追踪来计算合成旅行时，把该旅行时进行比较，得到一 个差值，用这个差值通过适当的算法修改模型。反复迭代上述过程，直到计算出的旅行 时与真实旅行时相差在允许误差范围内为止，处理流程如图3 2 所示。 3 2 离散图像重建技术 图3 2 地震层析成像处理流程图 在实际c t 反演中，首先确定成像区域，把成像区域划分成有限不重叠的象元，对 应每一个象元可以定义它的基，则图象的基表示为多维向量。根据有限维向量空间元的 表现定理可知，图象的数字化版为基的线性组合在每一象元内的平均值。如果象元足够 第三章地震c t 的理论基础及反演算法比较 小，则任意一条射线图象的投影函数可根据r a d o n 变换可知它在任意一个象元内的长 度。于是我们便得到了投影函数的级数展开式，它也给出了重建图象数字化版的方程组。 设f ( x ，y ) 在区域q ( 图3 3 ) 外恒等于零，即f ( x ，y ) 三0 ，( x ，y ) 萑q ，将区域分 割成1 个不重叠小区域s ，( 像元) ，i = 1 ，2 ，3 ，i 。f ( x ，y ) 在s ，上各点之值用它 ks ，上之平均值近似代替h 0 1 ，即 z ：f , f 矿( x , y ) d x d y ( 3 9 ) 其中m 表示小区域s ，之面积。l f i - j 量f = ( f l ，厶，乃) 7 被称为图像向量。 y 口 图3 3 反演区域离散示意图 x 设射线三与小区域s ，相交部分之长度为，根据r a d o n 变换，函数沿射线的 投影函数为 吼，厂= f ( 墨y ) 凼= 善i 彳， ( 3 1 0 ) 其中为射线总数。 在地震走时层析成像情况下，投影数据孵，f 为地震波走时f ，图像向量厂为像元内 慢度的平均值，则式( 3 1 0 ) 可写成矩阵方程 a j = 气 ( 3 1 1 ) 离散图像重建问题转化为给出一系列地震波走时r ，计算图像向量即地质体介质慢 1 2 长安大学硕士学位论文 度向量厂。 方程组( 3 1 1 ) 中的系数矩阵a 是极其稀疏的，因为它的每一行有f 个元素，而每条 地震波只通过所有1 个像元中的一小部分，因此矩阵a 中的大部分元素为零。根据系数 矩阵稀疏的特点，对方程组( 3 1 1 ) 多采用迭代方法求解。 3 3 反演算法 根据投影数据重建图像，它是一个产生具有某种物理性质的二维分布图像的过程， 其数学基础就是前面介绍的r a d o n 变换及逆变换。但是，r a d o n 逆变换公式应用到实际 中存在不少困难，例如： ( 1 ) 在c t 中，我们仅能得到有限的测量值，从有限的线积分去决定此图像，也就是说， 不完全数据投影是我们首先遇到的困难。 ( 2 ) c t 中得到的测量值，只能用以估算线积分值，这些估算中的不准确，是射线束的 宽度、散射及检波器的不准确性等等造成的，而r a d o n 逆变换对这些不准确性是敏感的。 ( 3 ) 对于r a d o n 逆变换公式不便作数值计算，我们需要一个有效的算法去代替它，并 希望它不是繁琐的，有利于提高数据计算的速度。 来解决上述困难，多年来人们作了大量的工作去研究和寻找各种算法。这些算法不 但应在计算机上完成得快，而且在所获得的数据有限或不准确时，能产生可接受的重建 图像，目前常用的用正问题解决反问题的迭代法，可统一用下式表达 b = a s ( 3 1 2 ) 式中，b = ( ，t 2 ，) r ，代表i 条射线走时残差；s = ( 而，s 2 ，s j ) r ，代表j 个慢度修 正量；a 即为式( 3 1 1 ) 中的系数矩阵。 通过非线性化处理，地震波旅行时层析成像实现的关键变为求解线性方程组b = a s 的问题。通过前面的分析，方程b = a s 具有以下特点： ( 1 ) 由于网格的离散细化，每条射线只通过少量网格，即只有相对少的矩阵元素是 非零的，因此矩阵a 是一个大型稀疏矩阵。 ( 2 ) 不完全投影决定了问题的欠定性，即模型参数不能唯一确定。 ( 3 ) 大多数情况下方程个数i 多于模型参数j ，即i j ，使得问题超定。 ( 4 ) 测量数据的误差决定了问题的不相容性，即很难找到能精确满足方程b = a s 的 s 矢量。 第三章地震c t 的理论基础及反演算法比较 由于b = a s 通常没有精确解，人们常常借助最小二乘解算法求取s ，它使得i ia s b 取极小值，其中i 表示欧几里得范数。该s 称为( 3 1 2 ) 式的最小范数最小二乘解，这也 就是人们通常要寻找的解。但用最d , - - 乘法直接求取具有大型病态稀疏系数矩阵的方程 组( 3 1 2 ) 时，要占去大量内存，耗时过多，而且矢量b 的误差严重影响最小二乘解。 基于上述原因，诸多专家经过研究，给出了一些有效的求解大型稀疏线性方程组 ( 4 1 3 ) 的算法，主要有代数重建技术( a r t ) 、联合迭代重建技术( s i r t ) 、正交分解最 d , - - 乘法( l s q r ) 等等，以下逐一介绍这些算法。 3 3 1a r t 算法 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 s 的简称，中文译为代数重建技术。用 a r t 求解方程组( 3 1 2 ) 是个迭代过程。 1 a r t 的递推公式 求解线性方程组( 4 1 2 ) 的解，从某个初值s o 开始，在迭代q 次之后近似解为s 训， 记残差为，- ( 引，于是 ，9 = 6 一a s g ( 3 1 3 ) 上式是一个递推方程，在递推过程中选择修正量厶( 鲋，使得残差为零，即 ，；叮+ 1 = 包一i s 9 + f 9 】= 0 ( 3 1 4 ) i 是q 除以i 的余数，本质上是逐条射线即逐个方程修正。 在特定的递推公式下，j 叮收敛于它的极限。这里对厶- 给定p 一范数。 渺= 睁州p r 列 因此得到修正量 掣= 警舻士p - l ( 3 1 6 ) i 鸣r 因为在地震层析成像中射线通过某个象素的长度4 ，0 。于是上式可写为 1 4 长安大学硕士学位论文 厶卜苎 鸣叶1_ 一v ( 3 1 7 ) 当p = 2 ，= 1 时有 一s j ( q ，+ 掣 ( 3 1 8 ) 鸣2 ( 3 1 8 ) 式就是a r t 算法的表达式。 式中i = l 2 i 为射线号，j = l ，2 j 为单元号，为系数矩阵a 中的元素，q = 0 ，1 为迭代次数。 2 a r t 实现步骤 ( 1 ) 选定一组初值芦，j = l ，2 ，j 。 ( 2 ) 计算第i 个方程( 即第i 条射线) 的估计值 b = 以_ ( 3 1 9 ) j = l ( 3 ) 计算观测值与估计值之差_ r j = 屯一b ( 3 2 0 ) ( 4 ) 计算第i 个方程的修正值 钙：乒 ( 3 2 1 ) 厶2 ( 5 ) 修正射线所通过的各个象素的慢度值 s j 2 sj 七沁j 并且根据约束条件 s m i 。 即 s j s 凇，令s j = s 呲 1 5 ( 3 2 2 ) ( 3 2 3 ) 第三章地震c t 的理论基础及反演算法比较 ( 6 ) 将( 3 2 1 ) 中修正后的的值代入( 3 1 9 ) 中，循环下去，直到完成i 条射线。 ( 7 ) 迭代结果满足预定残差的要求后，则停止，否则进入下一轮迭代。直到满足预订 残差要求为止，即完成第一轮迭代。 1r l 停止迭代 l 是 图3 4 a r t 算法流程图 1 6 长安大学硕士学位论文 3 3 2s i r t 算法 s i 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 的简称，中文译为联合迭代 重建技术。 1 s i r t 的递推公式 s i r t 的计算结果与使用数据的次序无关。这是计算方式上与a r t 法的不同之处。 s i r t 利用重建图像中的某一个像素所有通过它的射线的修正值来确定这个像素的平均 修正值，这样的平均值可以消除某些干扰因素。 从数学上考虑，我们不一定要求某个方程的残差消失，而是要使残差量呈递减趋势。 从物理上说，是对穿过图像重建中某一个像素的射线取平均修正量，即 掣2 古军嚣 ( 3 2 4 ) 式中，q 为循环次数，m j 表示a 的第j 列中非零元素的个数。从效果上，n n ( 3 2 4 ) 式可以消除某些干扰和随机测量误差。在计算上，该方法比较稳定。 2 实现步骤 ( 1 )n g - - 缈j j n0 ，j = l ，2 ，j 。s i r t 初值选择与a r t 的选择相同。 ( 2 )计算估计值 ( 3 )计算观测值与估算值的差 = 包一6 f ( 3 2 5 ) ( 3 2 6 ) ( 4 ) 求出第j 个像素内平均修正值衅，假设在像素内共有m 条射线通过，则有： 衅2 古；，；，弘 2 7 ) ( 5 ) 用平均修正值缸夕对第j 个像素的慢度值- 进行修正： 同时要加约束条件 哆“) _ 妒+ 衅 妒“ 1 7 ( 3 2 8 ) ( 3 2 9 ) r l 乙l = z g乃 j 一 = 岛 当 第三章地震c t 的理论基础及反演算法比较 s j s 啦，令s j = s 蕊 ( 6 )设e 为某个残差值，当 s j ( q 一妒“i 口，j = l 2 ，j 时停止迭代，否则进入下一轮迭代，重复( 2 ) ( 6 ) 。 长安大学硕士学位论文 1r l 停止迭代 l 是 图3 5s i r t 算法流程图 3 3 3l s q r 算法1 4 0 i l s q r 法最早是p a i g e 和s a u n d e r s 于19 8 2 年提出的，它是利用l a n c z o s 方法求解最 小二乘乘问题的一种投影法，由于在求解过程用到q r 因子分解法，故这种方法叫做 一1 9 - 第三章地震c t 的理论基础及反演算法比较 l s q r 方法( l e a s ts q u a r e sq r - f a c t o r i z a t i o n ) 。 l s q r 法是利用迭代方法求解舡_ b ，m i ne i a s - b l l ：的一种算法。这种方法基于g a l u b 和k a h a n 提出的双对角处理方法，通过分析此方法与标准共轭梯度法相等价，但是具有 数值计算上的优点。其基本原理如下：将i x j 矩阵a 用l a n c z o s 方法双对角化矩阵，即 用u 和v 正交阵使得 a = u b v l ( 3 3 0 ) 其中： b = a 1 0 0 卢2a 2 00 00 9 j 伐j ( 3 3 1 ) 再经过q r 分解求出( 3 3 0 ) 式的迭代向量序列，将其应用到( 3 1 2 ) 式就可得出l s q r 方法重建的慢度图像。下面分别介绍l a n c z o s 方法、q r 分解和l s q r 算法。 1 l a n c z o s 方法 l a n c z o s 方法是一种求解具有结称系数矩阵的大型线性方程组的一种投影方法。该 方法的基本过程是从包含n 个未知量的具有对称系数矩阵的线性方程组本身出发，递推 构造一列n 维空间尺”中的标准正交向量，同时形成r ”的一列单调子空间，且每个都能 方便地求出方程组在相应子空间内的投影近似解，当子空间的维数( 标准正交向量的个 数) 等于方程组的秩时，相应近似解就是所要求的一个解。 对方程 b x = b ( 3 3 2 ) 其中，b 为n * n 矩阵。设矿们，y n ，v ”1 为n 维空间中m 个无关向量。令 圪= i v 们，y ，v ”一1 ) 为n * m 矩阵，k m = 印口门 y ，y n ，v ( m - 1 ) ) 为y 们，y n ，v 川张 成的子空间。投影法的基本思想是寻找( 3 3 2 ) 式的近似解x ”使得： (。)1xy(m(，)，e，k：=m(bx b 。，1 ，2 ，m 1 ( 3 3 3 ) i 一) 上y ，= o ，2 ，一 r 7 令x ”= 圪y ，其中y ”为m 宰1 的实向量，则( 3 3 3 ) 式等价于： 长安大学硕士学位论文 ( v f b 圪) y ”= 吃6 ( 3 3 4 ) 假设( 3 3 4 ) 式有唯一的解，则通过求解( 3 3 4 ) 得到y ，从而得到x 。如何选择 y 们，v n ，v ”1 使得( 3 3 4 ) 式中的矩阵曙b 具有较简单的形式，使求解( 3 3 4 ) 式在计算 机上容易实现。下面给出如何构造l a n c z o s 向量y 们，y n ，v ”n ，使得矩阵曙b 圪具有 三对解形式。 v o = b l f l o ，风= 1 1bi i , v 一1 = 0 研= 0 ，1 ，2 ， ”= b v 一卢。v ”一1 ( 3 3 5 ) a 。= ( y ，”) v ”+ 1 = ( c o 一a 。v ”) 卢。+ l 上式a 。，成+ 。为标量，选取尾+ ，使得v 埘+ n1 1 = 1 ，如不再存在这样的凡+ 。则迭代停止。 按