（地球探测与信息技术专业论文）地震初至波走时层析成像与近地表速度建模.pdf

摘要 在整个地震勘探过程中，精确地求取地震波在地下介质中的传播速度，一直是地震 勘探的核心问题之一，尤其是近地表速度。 近地表速度建模的研究方法有很多，且各有其优劣，如折射法、面波法、层析法等。 初至波旅行时层析成像是指利用在地面某点放炮，在另一点接收的地震资料的初至波旅 行时信息反演炮点、接收点之间的近地表的地震波速度或地质构造。本文研究利用初至 波旅行时层析成像方法来建立近地表的速度模型，提供近地表地下介质的速度信息，进 一步为静校正或浅层工程勘探服务。 论文简要回顾了地震层析成像技术的发展历程，论述了地震初至波旅行时层析成像 的基本原理，介绍了初至波旅行时射线追踪正演算法和层析成像反演算法，进行了主要 算法推导，编程实现了初至波旅行时线性插值( l t i ) 射线追踪正演算法和联合重建迭 代( s 嬲) 层析反演算法。 在程序实现的基础上，着重试算分析了影响层析成像效果的3 个基本因素：初始速 度模型、速度约束限制和射线覆盖密度，并总结出各影响因素优选的原则和范围，以达 到提高初至层析反演精度及反演稳定性的目的；试算了包括水平界面、倾斜界面、弯曲 界面、断层等7 组层状介质理论模型、1 种加随机噪音的理论模型和两组野外实际数据， 较满意的试算结果表明：1 ) l t i 正演算法较之其他射线追踪算法编程简单、精度高、速 度快；2 ) s i r t 反演算法能得到稳定收敛的解；3 ) 通过对加随机噪音模型的试算可知， 文中的正、反演算法具有一定的抗噪能力；4 ) 比较满意的理论模型和实际资料处理结 果证明了程序的可靠性和实用性。 关键词：速度建模，层析成像，初至波，射线追踪 a b s t r a c t t h ea c c u r a t ev e l o c i t ym o d e lb u i l d i n gi st h ek e yf o rs e i s m i ce x p l o r a t i o n ，e s p e c i a l l y b u i l d i n gn e a r - s u r f a c ev e l o c i t ym o d e l t h e r ea r em a n ym e t h o d so fb u i l d i n gn e a r - s u r f a c ev e l o c i t ym o d e l ，s u c ha sr e f r a c t i o n m e t h o d ，s u r f a c ew a v em e t h o d ，t o m o g r a p h ym e t h o d ，a n ds oo n f i r s tb r e a kt o m o g r a p h yi ss h o t i no n ep l a c ea n dr e c e i v et h ef i r s tb r e a kt r a v e l - t i m ei na n o t h e rt oi n v e r s et h ev e l o c i t yo f s e i s m i cw a v eo rt h es t r u c t u r eb e t w e e ns h o ta n dr e c e i v e r t h r o u g ht h i sm e t h o d ，t h ep a p e r r e s e a r c h e sn e a r - s u r f a c ev e l o c i t ym o d e lb u i l d 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 nf o rf u r t h e r s t a t i cc o r r e c t i o no rs h a l l o w p r o j e c ts e i s m i ce x p l o r a t i o n t h ep 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 tc o u r s eo ft h es e i s m i ct o m o g r a p h yt e c h n o l o g y , a n dd i s c u s s e st h ep r i n c i p l eo ft h ef i r s tb r e a kt r a v e l - t i m et o m o g r a p h y w ei n t r o d u c et h er a y t r a c i n gm e t h o do ft h ef i r s tb r e a kt r a v e l t i m ea n di n v e r s i o na l g o r i t h mo ft o m o g r a p h y , d e d u c i n g p r i m a r ya l g o r i t h m ，p r o g r a m m i n gt h er a yt r a c i n gm e t h o do fl t ia n dt h es i r ti n v e r s i o n a l g o r i t h m o nt h i sb a s e ，w es t u d ya n ds u m m a r i z ef o r w a r d - i n v e r s i o ne f f e c ta n di t sp a r t i a li n f l u e n t i a l f a c t o r s ，s u c ha si n i t i a lv e l o c i t ym o d e l ，v e l o c i t yr e s t r a i n tc o n d i t i o n ，d e n s i t yo fr a yc o v e r t e s t s w i t hs e v e nt h e o r ym o d e l s ，o n et h e o r ym o d e lw i t hr a n d o mn o i s e ，a n dt w or e a ld a t as h o wt h a t ： i ) t h el t ir a yt r a c i n gm e t h o di se a s i e r , f a s t e ra n dm o r ep r e c i s et h a nt h eo t h e rm e t h o d s 2 ) s i r ta l g o r i t h mc a np r o v i d es t a b l ea n dc o n v e r g e ds o l u t i o n 3 ) t h r o u g ht h et e s t so ft h e o r y m o d e lw i t hr a n d o mn o i s e ，i tc a nb ek n o w nt h a tt h ef o r w a r d - i n v e r s i o na l g o r i t h mi nt h ep a p e r h a sa n t i - n o i s ea b i l i t y 4 ) t h ep r o g r a m sr e l i a b i l i t ya n du s a b i l i t yh a v eb e e np r o v e db yt e s t i n g w i t ht h e o r ym o d e l sa n dr e a ld a t a k e yw o r d s ：v e l o c i t ym o d e lb u i l d i n g ，t o m o g r a p h y , f i r s ta r r i v a l ，r a yt r a c i n g 论文独创性声明 本人声明：本人所呈交的学位论文是在导师的指导下，独立进行研究工 作所取得的成果。除论文中已经注明引用的内容外，对论文的研究做出重 要贡献的个人和集体，均已在文中以明确方式标明。本论文中不包含任何 未加明确注明的其他个人或集体已经公开发表的成果。 本声明的法律责任由本人承担。 论文作者签名：景司红 砷年，夕月彳日 论文知识产权权属声明 本人在导师指导卞所完成的论文及相关的职务作品，知识产权归属学 校。学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权 利。本人离校后发表或使用学位论文或与该论文直接相关的学术论文或成 果时，署名单位仍然为长安大学。 ( 保密的论文在解密后应遵守此规定) 论文作者签名： 导师签名： 景坷红 象谣 似7 年f 月7 日 加1 年j 其缉b 长安大学硕士学位论文 1 1 研究目的和意义 第一章绪论 地震勘探是利用人工在地表激发和接收地震波，再对地震波作分析处理以及解释而 得到地下构造信息和岩性信息的一种方法。在整个地震勘探过程中，精确的求取地震波 在地下介质中的传播速度，一直是地震勘探的核心问题之一【1 1 。尤其在地表条件较复杂 的区域，地表速度的横向剧烈变化会严重影响中深层目的层的成像效果。近地表速度不 准确，将会直接影响到速度分析、偏移成像的质量以及静校正的精度等地震勘探的各个 环节和最终的勘探成果。因此，建立准确的近地表深度速度模型已经成为地下准确成 像的迫切需要【2 】。由于直接测量速度场只能在很小的局部进行，无法得到整个地下的速 度场结构，因此地球物理学家们一直致力于利用地震波本身的信息，在各种假设前提下 进行地下速度结构的反演。 通常使用的速度求取方法有很多：如叠加速度分析和偏移速度分析，这两种速度分 析方法由于覆盖次数不足、缺少近地表反射，一般很难取得近地表层的有效速度信息； 又如折射波法，人们以往对于浅层一般都是应用该方法进行速度分析，但该方法由于受 近地表速度横向变化、实际地震记录中折射波识别困难、折射界面不稳定等因素的影响， 其结果往往难以令人满意；而微测井等直接测量近地表速度的方法，虽然能够较准确地 确定测点处表层速度随深度变化的规律，但在低速带厚度和速度变化较大的地区，必须 加密测点以保证测量速度的准确性，从而使得生产效率降低，成本提高。 几年来，在近地表速度建模的问题中，利用初至波的层析成像方法越来越受到人们 重视。地震勘探中的初至波是指从炮点出发，经过介质首先到达接收点的地震波，可能 是直达波、回转波( 连续折射波) 、折射波或多个层折射波的组合：这些波主要在近地表 层传播，因此，其走时包含了近地表层的速度信息。而且，直达波主要体现均匀介质模 型，回折波主要体现连续介质模型，折射波主要体现层状介质模型，因此，利用初至波 走时信息对近地表速度进行精细的分析最为有利；又因初至波容易获得，往往能量强， 可追踪性好，给波的识别和到达时拾取带来方便性和可靠性；因此，人们一直对用初至 波旅行时层析反演方法进行近地表速度建模情有独钟。本文主要研究利用初至走时的层 析反演方法建立近地表速度模型，提供近地表地下介质的速度信息，进一步为静校正或 浅层工程勘探服务。 第一章绪论 1 2 层析成像发展史 层析成像的英文t o m o g r a p h y 源于希腊语t o m o s ，本意是断面或切片的意思。在实际 应用中，层析成像就是一种利用在探测目标表面观测到的信号来求取目标内部信息的反 演方法，最早应用于医学领域，随后扩展到其它领域。地震层析成像的主要目标是确定 地球内部的精细结构和局部不均匀性。这不仅可以促进地球物理科学的发展，而且还可 以解决许多地质勘探中和矿产资源开发中的难题【3 一钉。 地震层析成像的研究自7 0 年代以井间速度结构调查为研究对象在地学界出露头角 ( b o i s e ta l ，1 9 7 2 ) p j 。d i n e s 和l y t l e ( 1 9 7 9 ) 首先对地震层析成像做了大量数值模拟，并公布 了利用弯曲的地震波射线进行地下地震波速度成像的结果，并首先将层析成像 ( c o m p u t e r i z e dg e o p h y s i c a lt o m o g r a p h y ) 这- - 名词用于论文的标题【6 】。8 0 年代，地震层析 成像的研究在地震勘探领域得以发展，自从在亚特兰大( a t l a n t a ) 召开的第5 4 届地球物理 勘探学家协会( s e gs o c i e t yo fe x p l o r a t i o ng e o p h y s i c i s t s ) 年会上设置了地震层析成像研究 内容的专题之后，以d a i l y t 7 1 ( 1 9 8 4 ) ，s o m e r s t e i n 8 1 ( 1 9 8 4 ) ，b i s h o p 9 1 ( 1 9 8 5 ) ，d y e r 和 w o r t h i n g t o n 1 0 1 ( 1 9 8 8 ) 等人的研究为代表，利用人工地震发射与接收系统的地震层析成像 的理论、方法和技术以数值模拟的形式得到广泛研究。9 0 年代，不论是利用人工地震数 据还是天然地震数据的地震层析成像方法在工程勘察、资源勘探、环境保护、文物调查、 岩体结构研究等许多应用领域得到实验性研究【1 1 1 2 1 。 我国的地震层析成像工作始于8 0 年代初期。其中，金安蜀【l3 】等( 1 9 8 0 ) 首先对北 京地区速度结构作反演并进行了地震层析成像方法研究，使地震层析成像方法得以开始 在中国应用。应用之初，多为对区域结构的层析成像：刘福田【1 5 】等( 1 9 8 6 ，1 9 8 9 ) 利 用近震和远震的资料对中国多个地区和整个中国大陆及邻近地区进行层析成像，取得了 一系列的成果；朱露培【l6 】( 1 9 9 0 ) 提出一种在频率域展开的模型参数的描述方法，并应 用于北京遥测台网的远震层析成像中，发现京津唐地区的上地幔有一低速区；等等。9 0 年代以来，随着层析理论的日趋成熟，地震层析成像在工程勘察、资源勘探、岩体结构 研究等许多领域得到实验性研究并取得有效的进展。 1 3 射线初至走时层析 在地震勘探问题中，层析成像就是利用在地表或井中接收记录到的地震波的旅行 时、振幅和波形等信息，来重建地下介质的速度等参数的分布技术。 2 长安大学硕士学位论文 从不同的角度，地震层析成像可以进行不同的分类【1 1 7 】1 ) 根据地震勘探中的布局 方式不同，可以分为地表反射( 或折射) 层析成像、井间层析成像、v s p 层析成像和反 向v s p 层析成像；2 ) 根据所利用的数据类型不同，可分为反射层析、透射层析和折射 层析；3 ) 根据所利用的地震剖面或炮集上的地震属性的不同，可以分为走时层析、振 幅层析和波形层析；4 ) 根据层析成像方法所依据的理论基础不同，可以分为以射线理 论为基础的射线层析成像和以波动理论为基础的绕射层析成像。在实际中，层析方法所 依据的理论基础和所使用的地震信息常常是紧密结合在一起的，比如通常所说的旅行时 层析成像就是在射线理论基础上利用地震走时信息的层析成像方法；而以波动理论为基 础绕射层析成像利用的则是地震波的波形信息。 所谓基于射线理论的初至走时层析成像是指利用地震波初至走时和其传播的射线 路径来反演地下介质速度的技术。近些年来，利用基于射线理论的地震初至波层析成像 方法重构近地表速度模型成为人们的研究热点。国内许多学者作过此方面研究：刘振宽 1 8 1 等( 1 9 9 4 ) 将近地表模型离散成三角形单元，每个单元内的速度线性变化，用试射法 的两点射线追踪技术完成炮点至接收点的射线追踪，并用最小二乘共扼梯度法进行反演 求解，理论模型试算证明方法具有可行性；李录明【1 9 】等( 2 0 0 0 ) 把近地表模型离散成矩 形单元，每个单元内的速度用由单元网格点速度得到的双线性函数表示，用网络法进行 射线追踪正演，用带阻尼的最小o r 分解迭代算法求解大型稀疏矩阵方程，理论模型和 实际资料都得到很好的效果；李满树【2 0 】等( 2 0 0 4 ) 将初至走时层析反演用于复杂地区近 地表速度建模，正演采用变步长扰动法，反演采用求解最小二乘问题的豪斯荷尔德变化 法，理论模型和实际资料试算均取得理想的效果；陈爱萍【2 l 】等( 2 0 0 6 ) 采用改进的最短 路径射线追踪算法进行正演，利用改进的约束带阻尼联合迭代重建技术求解大型、稀疏 方程组来反演近地表速度场，四川复杂地区某条测线的处理结果表明该方法能稳定、快 速、准确地重建复杂地区的近地表速度场；刘玉柱【2 2 】等( 2 0 0 7 ) 通过正则化手段将各种 先验信息纳入到初至波层析反演方程中，替代传统的外部约束模式，理论模型和实际资 料试算表明，正则化方法比传统外部约束方法具有更高的分辨率和反演精度；韩晓丽【2 习 ( 2 0 0 8 ) 等采用非显示射线追踪方式在全偏移距内进行反演，降低了速度深度模型的 不确定性，提高了表层模型的精度，适合任意介质，很好的解决了复杂地表引起的的静 校正问题；等等。 本论文主要着眼于基于射线方法的初至波走时层析反演。 1 1 射线追踪 第一章绪论 射线旅行时层析成像是利用地震波走时和其传播的射线路径来反演地下介质速度 的技术，它是目前在地震勘探领域应用最广泛的层析成像方法。射线追踪正演算法的精 度和速度直接影响着层析反演的精度和效率【2 4 1 。 地震射线追踪的方法很多，不同的方法适用的条件不同，计算效率也不同，而研究 适应任意构造和速度变化较大的情况，且计算效率高的射线追踪方法对于地震层析成像 技术的发展有着极其重要的意义。地震层析成像研究中，通常用由许多常速矩形单元组 成的二维单元模型来描述速度结构。对于这种模型，常规试射法不能很好地解决问题， 主要问题一是它无法处理折射波( 首波) ，二是速度结构复杂时可能造成射线的漏失， 三是计算速度慢。因此，前人提出了适应复杂条件的多种射线追踪或地震旅行时计算方 法，如用有限差分求解程函方程的方法【2 孓2 7 1 、最短路径法 2 8 - 3 2 】、波前法【3 3 1 等。1 9 9 3 年 a s a k a w a 提出了一种称为旅行时线性插值( l t i ：l i n e a rt i m ei n t e r p o l a t i o n ) 的新射线追 踪算法【3 4 】，通过试验表明用该算法计算旅行时和追踪射线路径比其它常规方法( 如有限 差分法) 更为快速、精确。因此，本文采用旅行时线性插值( l t i ) 射线追踪正演算法 2 ) 地震走时 地震走时是地震图中最有用的信息之一。与地震图的振幅和波形信息相比，地震走 时具有以下特点：1 ) 它主要取决于地下介质的波速结构，而与激发接收条件和仪 器性能关系不密切，因此通常具有较高的信噪比；2 ) 对球面波和柱面波的传播，其振 幅衰减的规律是不同的，但走时规律相同，因此利用走时反演地球波速结构的方法可同 时用于二维、三维及点源加上二维介质的所谓二维半的情况。地震走时的这些特点使利 用走时的层析成像方法具有简单、直观和通用的优点。 3 ) 初至走时 理论上，不同类型的波均可用于层析成像速度函数的反演。 由于浅层反射的地震信息分布于近道，且被面波严重干扰，难以用其求准速度，因 此，我们以往对于浅层一般都是应用折射波进行速度分析。但是，由于用折射波对表 层速度进行分析比较适合于介质为层状分布的地区；在实际地震记录中往往很难将折 射波信息识别出来；对于新测区，浅层折射界面的位置也是很难确定的，所以我们认 为用初至波走时信息来反演近地表速度更为合适【2 0 】。 我们定义初至波为从炮点出发，经过介质最先到达接收点的地震波，有可能是直达 波、回折波或者折射波；在均匀介质中主要出现的可能是直达波；在连续介质中主要出 现的可能是回折波；而在层状介质中主要表现出的就可能是折射波了。通过三者的组合 4 长安大学硕士学位论文 以及层析方法对介质横向变化的适应性，从而使得初至波层析成像方法适应任意表层模 型的反问趔1 9 】；另外，由于在浅层工程地震勘探中，初至容易获得，且往往能量强，可 追踪性好，给波的识别和到达时拾取带来方便性和可靠性；因此，利用初至波走时信息 对近地表速度进行精细的分析最为有利。 1 4 本文主要研究内容 地震层析成像方法主要分为两大部分，即正演方法和反演方法。正演方法一般从两 方面入手，一是射线追踪方法，二是对波场进行数值模拟；反演方法可以分为两类：第 一类是基于算子的线性化或拟线性反演方法，又称为“确定性反演方法”，如a r t 法、 s i r t 法、l s q r 法、最速下降法、拟牛顿法、共轭梯度法等；另一类是基于模型的完全 非线性反演方法，又称为“随机反演方法 ，如m o n t ec a r l o 法、模拟退火( s a ) 法、神经 网络法、遗传算法( g a ) 、单纯形法以及混沌算法等【3 5 1 。 结合层析成像正演的基本理论和层析成像反演的基本算法，本文主要对初至波旅行 时线性插值( l t i ) 正演射线追踪方法和联合迭代重建( s i r t ：s i m u l t a n n e o u si t e r a t i v e r e c o n s t r u c t i o n ) 层析反演算法展开研究工作。 研究的主要内容包括： 1 ) 对地震初至波旅行时层析成像的理论基础r a d o n 变换做了简单介绍，在此基础 上分析了地震初至波旅行时层析成像实现的基本原理，并绘制实现流程图；详细推导了参 初至波旅行时线性插值射线追踪算法，介绍了层析成像三种常用的反演算法，包括代数 重建技术( a r t ) 、联合迭代重建法( s i r t ) 、最小平方q r 因子分解法( l s q r ) 。 2 ) 利用f o r t r a n 语言自行编写并实现了初至波l t i 正演射线追踪程序；该程序可以 适用于多种观测系统，如地面放炮地面接收、地面放炮井中接收、井中放炮井中接收、 井中放炮地面接收、左边放炮右边接收、右边放炮左边接收，可以单边放炮，也可以双 边放炮；正演程序中的速度建模模块可快速精确地构建常用的地质模型，如断层、透镜 体等；l t i 初至波旅行时追踪程序在网格剖分足够小的条件下可以适用于包括水平层状、 弯曲层状等任意构造和速度变化较大的情况的初至波旅行时射线追踪，程序运行时间 短，并能得到精确的结果； 3 ) s i r t 属于迭代算法，该迭代算法能稳定收敛，收敛速度快，本文利用f o r t r a n 语言自行编写并实现了s i r t 层析成像反演程序； 4 ) 影响地震层析反演效果的因素有很多，本文以理论模型为基础详细分析了初始 第一章绪论 速度模型、速度约束限制、射线覆盖密度三个影响地震层析成像效果的因素，并总结出 各影响因素优选的原则和范围，达到提高初至层析反演精度及反演稳定性的目的； 5 ) 理论模型初至波层析成像反演结果分析：结合实际地质情况，设计了包括水平 界面、倾斜界面、弯曲界面、断层等7 组层状介质理论模型，进行正、反演方法试算及 效果分析，验证方法的正确性和有效性； 6 ) 噪音影响层析成像结果的试算：利用f o r t r a n 语言对初至波加随即噪音，之后对 加噪音后的初至波提取初至时间进行层析反演试算和结果分析，检验方法的抗噪能力； 7 ) 野外实际数据处理和反演结果分析：人工拾取实际初至时间，根据工区的叠加 速度资料和已有的其他地球物理资料( 如面波) 确定初始速度模型和速度约束范围，在 此基础上进行层析反演，获得近地表速度模型信息，并与其他反演结果( 如面波) 对比， 验证方法的实用性和可靠性性。 地震初至波射线层析成像方法实现的流程图如图1 1 所示： 图1 1 初至波旅行时地震层析成像流程 6 长安大学硕士学位论文 第二章地震初至波旅行时层析成像的基本原理 本章主要介绍地震初至波旅行时层析成像的理论基础和层析成像实现的基本原理。 2 1 地震初至波旅行时层析成像的理论基础【4 ，1 1 ，1 2 ，3 6 l _ 求a d o n 变换 各种层析成像的数学理论基础都是经典r a d o n 变换或广义r a d o n 变换。经典r a d o n 正变换是沿直线簇对某已知函数的积分，广义r a d o n 正变换是经典r a d o n 正变换的推广， 它是沿某种曲线簇对某已知函数的积分。r a d o n 反变换是重建r a d o n 正变换中所涉及的 函数，或者恢复该函数的某种特性。在地震勘探中，当地下介质的层与层之间的速度变 化不大时，即波阻抗差异较小时，初至波射线路径可以用直线近似，此时，可将初至波 旅行时看作是地下介质的地震波慢度( 速度的倒数) 函数的经典r a d o n 正变换，通过经 典r a d o n 反变换可由初至波旅行时重建地下介质的慢度函数。当地下介质的波阻抗差异 较大时，初至波射线路径是曲线，不再能用直线近似，这时，应将初至波旅行时看作是 地下介质的地震波慢度函数的广义r a d o n 正变换，由初至波旅行时重建地下介质的慢度 函数的问题属于广义r a d o n 反变换问题。地震层析成像是以经典r a d o n 变换为基础的， 只有在少数地震条件简单的地区，可以直接应用经典r a d o n 变换，而在大多数情况下， 岩性分布很不均匀，不能用直射线进行层析成像，只能利用象素划分和射线追踪技术实 现曲射线层析成像。 下面给出二维经典r a d o n 变换。 设u ( x ，z ) 是在x o z 平面上定义的充分光滑的二元函数，在充分大的区域之外为零，三 为x o z 平面内的任意一条直线，则称u ( x ，z ) 沿直线三的积分u ( x ，z ) 的r a d o n 正变换，记 为 r u = p ( 础) d s 工 ( 2 1 ) 其中，凼表示线长元素。通常“为图像，称r “为图像的投影。 x o z 平面内任意直线( 如图2 1 所示) 可表示为 三：p = x c o s o + z s i n o( 2 2 ) 其中，( c o s o ，s i n 0 ) 为三的单位法向量，口为三的法线与x 轴正向的夹角，p 为坐标原点 到三的距离。由此可见，任意直线三对应确定的数对( p ，秒) 对应x o z 平面内的确定的一条 7 第二章地震初至波旅行时层析成像的基本原理 直线( 当0 0 , 2 7 r ) 时，这种对应是一一对应) 。所以u ( x ，z ) 的r a d o n 正变换也可以表示 为 r u ( p ，口) = j u ( x , z ) d s ( 2 3 ) p = x s o + z s i n o 二维r a d o n 反变换，是指根据函数u ( x ，z ) r a d o n 正变换的值( 即许多的积分值) 来求 被积函数u ( x ，z ) 。二维r a d o n 反变换，也称为二维r a d o n 反演公式，公式如下： r u ( p ，9 ) 嘶) - - 专p 岛去力 ( 2 4 ) 图2 1直射线投影表示示意图图2 2 象素划分与射线穿过象素的长度示意图 2 2 地震初至波旅行时层析成像实现的基本原理 初至波旅行时层析成像反演就是给定初始模型并网格化，计算初至走时并将其与实 际观测的初至走时比较，根据比较结果修改模型，再计算走时，然后再比较，如此反复 比较修改直至迭代达到满意的精度，最终得到的模型就是反演结果【3 7 】。 具体实现过程如下【3 8 】： 在地震初至波旅行时层析成像中，首先将要重建地震波慢度图像的区域进行网格划 分，每个网格称为一个象素，每个象素可以是正方形、长方形，或者三角形的。如图2 2 所示。假设用长方形的象素覆盖要重建图像的区域，按一定的顺序排列后，共有j 个象 素，且每个象素内地震波的慢度是常数，用s ，( = 1 , 2 ，) 表示。对于每个 j ( j = 1 , 2 ，j ) ，定义一个如下的图像基函数 州础) = 熙篙多g 外j - ? - 羹蓁 ( 2 5 ) g ，( 础) 2 饭掷) 叠 象素 【2 5 ) 长安大学硕士学位论文 于是，地下介质的地震波慢度图像s ( x ，z ) 可用s ，和g ，( 戈，z ) ( = 1 , 2 ，力的线性组合 来近似，即 - ， s ( x ，z ) = s ，g 小，z ) j = l ( 2 6 ) 假设某次地震实际观测共记录了，道，每道对应一条从炮点到接收点的初至波射线 路径，用厶g = 1 , 2 ，d 表示，每道也对应一个初至波旅行时，用饥o = 1 , 2 ，) 表示。 而仇等于地震波慢度j ( x ，z ) 沿第f 条射线的曲线积分，因此以o = 1 , 2 ，) 可看作s ( 石，z ) 的广义r a d o n 正变换，即 包= 扣，z ) 凼= i 喜俐础肛蔷jj ，胁础油( 2 ，d ( 2 7 ) 根据g - ，( 工，z ) 的定义可知，g ，( z ，z ) 凼为第f 条射线在第，个象素内的射线长度，用 a u ( f = 1 , 2 ，歹= 1 , 2 ，刀表示这个积分。于是，通过( 2 7 ) 式可得如下的方程组 a l l口1 2 a 2 1a 2 2 口na i 2 a l la t 2 b l 6 2 i b j ； b i ( 2 8 ) 或写为a s = 6 ，其中，彳= ( ) 是由各条射线在各个象素内的长度构成的矩阵；s 是由 各象素内的地震波慢度构成的向量，s = ( 黾，j ：，曲) z ；b 是由各地震接收道对应的初 至波旅行时构成的向量，即b = ( b l ，6 ：，6 ，) r 。 当第七条射线不穿过第，个象素时，a 盯= 0 ，因此，一般地，系数矩阵彳是大型的 稀疏矩阵。 地震初至波旅行时层析成像实现的基本原理就是根据已知的地震初至波旅行时 岛( f = 1 , 2 ，d 设法从方程组( 2 8 ) 中解出慢度向量s ，( 歹= 1 , 2 ，j ) 。 在直射线层析成像中，方程组( 2 8 ) 的系数矩阵很容易求出来，且是不变的，因此该 方程组成为一个线性方程组，可用后面章节将要介绍的层析成像反演算法直接求解此方 程组，得到慢度向量s ，( = 1 , 2 ，) 。 9 屯；0；曲 矿iiiiiiijiii皿 u ； u； 口 口 口 口 口 卯鲥；砌；砌 第二章地震初至波旅行时层析成像的基本原理 在曲射线层析成像中，慢度向量j ，( 歹= 1 , 2 ，) 为未知，射线路径也是未知的，即 方程组( 2 8 ) 系数矩阵中的每个元素口甜都是未知的，因此，从方程组( 2 8 ) 中解出慢度向量 j ，( ，= 1 , 2 ，) 的问题是一个非线性问题。解决这个非线性问题的一般做法是：先给定 t ! 曼度向l i s ，( ，= 1 , 2 ，) 的初始值s ；( j = 1 , 2 ，刀，通过正演射线追踪求出方程组( 2 8 ) 的系数矩阵么的近似矩阵，及初始慢度下的初至波旅行时酽( f = 1 , 2 ，i ) 与已知的初至 波旅行时岛( f = 1 , 2 ，i ) 的差向量0 ”= b j 一印”( i = 1 ，2 ，i ) ，然后，仅将慢度向量 s ，( 歹= 1 , 2 ，) 看成未知量，方程组( 2 8 ) 就成为了线性方程组，用层析成像反演算法解 此线性方程组得到慢度修正量血；，从而得到新的慢度向量s ；+ 1 = s ；+ 厶；( 歹= 1 , 2 ，) ， 将修正后的慢度向量作为新的初始值，再重复以上的过程，直到初至波旅行时误差减 到足够小为止。具体实现流程图见图2 3 ： 图2 3 地震初至走时层析反演实现流程图 1 0 长安大学硕士学位论文 第三章地震初至波旅行时线性插值射线追踪方法 射线追踪在地震勘探中被广泛应用于地震正演模拟、a v o 分析、基于模型的地震速 度分析、叠前深度偏移和地震层析成像之中。地震射线追踪的方法有很多，不同的方法 适用的条件不同，计算效率也不同，而研究适应任意构造和速度变化较大的情况，且计 算效率高的射线追踪方法对于叠前深度偏移及地震层析成像技术的发展有着极其重要 的意义。 在地震勘探研究中，特别是在叠前深度偏移和地震层析成像研究中，通常用由许多 常速矩形( 或正方形、三角形) 单元组成的二维单元模型来描述地下介质的速度结构。 对于这种模型，常规试射法不能很好地解决问题，主要问题是：1 ) 它无法处理折射波 ( 首波) ，2 ) 速度结构复杂时可能造成射线的漏失，3 ) 计算速度慢。因此，前人提出 了适应复杂条件的多种射线追踪或地震旅行时计算方法，如用有限差分求解程函方程的 方法【2 5 之7 1 、最短路径法【2 8 。2 1 、波前法【3 3 】等。 1 9 9 3 年a s a k a w a 提出了一种称为旅行时线性插值( l t i ) 的新射线追踪算法【3 4 1 ，它 是以正反向方法为基础，并假定旅行时变化呈线性。在正向处理( 亦称为向前处理) 过 程中，计算网格单元边界上各个离散点处的初至波旅行时；在反向处理( 亦称为向后处 理) 过程中，根据费马原理追踪所要寻找的射线路径，同时求出射线路径与其穿过的单 元边界的交点。旅行时线性插值( l t i ) 射线追踪方法主要适用于初至波旅行时的计算 和射线路径的追踪。 旅行时线性插值( l t i ) 射线追踪方法通过对模型网格进行剖分，能有效计算包括 直达波、折射波、透射波在内的多种初至波的旅行时间和射线路径。该算法可以适应任 意变速介质，且不受“屏蔽层 的限制；另外，l t i 算法还有一个很大的优势，即它不 是简单地将边界节点的连线作为射线路径，而是通过线性插值，在满足f e r m a t 原理的 条件下的精确点处与边界相交，因此常速单元内的射线路径总是直线；同时边界处的折 射角度也是随着入射角而连续变化。 本论文在运用旅行时线性插值( l t i ) 正演算法计算地质模型初至波走时的同时， 将其与常规v i d a l e 法以及传统的打靶方法做了比较。通过试算表明，用旅行时线性插值 ( l t i ) 算法计算初至波旅行时和追踪射线路径比其它常规方法( 如有限差分法等) 更 为快速、精确。 第三章地震初至波旅行时线性插值射线追踪方法 3 1 模型描述 正演模拟技术在地震勘探领域发挥着非常重要的作用，而地质模型的建立又是正演 模拟技术的基础和前提。地质模型的建立归结为对地质模型结构的数学描述，本文通过 界面参数和地层速度参数进行地质模型描述。此模型描述方法更注重对地下结构的描 述，它适用于地下波场特征正演分析与研究。 3 1 1 界面方程的确定 图3 1 二维层状结构模型 地下沉积地层是通过漫长的地质过程形成的，普遍具有成层的特点。我们把地下介 质当作一个二维非均匀性、具有弯曲界面的层状介质。图3 1 所展示的是一个具有三个 分界面的二维层状结构模型。 二维层状结构模型可以描述为由贯穿模型左右边界的地层分界线分隔成许多单连 通的条带状区块，并且同一区块内具有统一属性的地层。其界面可用一次、二次、三次 或多次函数表示如下： z ( x ) = 口f 石 f = o 在上式中，z ) 是表示界面深度的函数，x 表示水平方向的空间距离，a i 是地层分 界线的拟合系数，以是拟合次数( 一般的，三次拟合就可以很好地描述界面的形态，地 质情况也比较吻合) ，具体意义如下： n = 0 ：表示界面为水平界面； r t = 1 ：表示界面为线性拟合( 线性非水平界面) ； 1 2 长安大学硕士学位论文 刀2 ：表示界面为弯曲界面； 本文在编写程序过程中，通过给定界面的控制点坐标来求取拟合系数a ；，对于弯曲 光滑界面利用界面控制点通过三次样条函数来建立界面方程。图3 1 中三个界面的控制 点坐标如下( 可参照图3 2 ( a ) ) ： 第一个界面：2 个控制点( o ，6 ) 、( 2 0 0 ，6 ) 第二个界面：3 个控制点( o ，1 2 ) 、( 9 0 ，1 6 ) 、( 2 0 0 ，1 0 ) 第三个界面：4 个控制点( 0 ，2 6 ) 、( 5 0 ，2 2 ) 、( 1 3 0 ，2 8 ) 、( 2 0 0 ，2 0 ) 3 1 2 速度界面的近似 本论文是在每个网格内的速度均为常速的基础上进行程序编写的，因此建立地质模 型时，必须保证此条件的成立。为此，采取以下方法对界面进行近似。 假设地质模型由n 个界面组成，用网格( 本文中采用正方形或者矩形) 进行剖分； 将第m 个界面的界面方程与每一条纵向的网格线求交点，如( x 1 ，y 1 ) 、( x 2 ，y 2 ) ，计算相 邻两交点连线的中点坐标( x 3 , y 3 ) ，根据此相邻两交点判断所属网格的位置，并求出该网 格中心坐标( x 4 ，y 4 ) ，判断：若y 3 y 4 ，即网格中心刚好位于真实界面上或位于真实界 面下方，则该网格属于下一层，若y 3 y 4 ，即网格中心位于真实界面上方，则该网格 属于上一层，依据所有这样的交点进行判断，则得到了满足网格内速度为常速的近似界 面；如图3 2 ( a ) 所示，其中虚线表示真实界面，实线表示近似界面。为了说明方便，图 3 2 ( a ) 中网格剖分所选择的矩形网格较大，导致近似界面与真实界面相差很大；当网格 剖分所选择的矩形网格足够小时，近似界面与真实界面相互吻合，见图3 2 ( b ) 。为方便 比较，不同网格剖分大小得到的近似界面在同一网格剖分下绘图。 ，弟一层v l ，v v 、 ”1 叩q 7 ( x y )k ￥z ) 一f 一 二r 。一一一j第二层： ，x4 - 一- 一一l 一一一 一 1 j 一一。一一 一5 ：第三层 v 3山衅二# 。一王7 。 一，一t 一： 一1 、一、。i j _ z r r i ，1 。7i 、- 一- 一一：二，卜一 o2 040(101 1 01 0 01 2 01 4 0100伯o2 0 0 3 2 ( a ) 网格剖分及界面近似示意图( 大网格剖分：1 0 2 ) 第三章地震初至波旅行时线性插值射线追踪方法 弟一一层 v l r 一，一 一一 - 、 第二二屡 v 9 _ l 一一 - _ 一 - _ 一 - - _ _ 昂工层 3 ，r、i 一r ， 1 、，_ 一一 h 一一_ 一一一 、 o2 04 01 1 08 0 1 0 01 2 01 4 01 1 0t 1 1 02 0 0 3 2 ( b ) 网格剖分及界面近似示意图( 小网格剖分：0 5 x 0 1 ) 3 2 旅行时线性插值射线追踪的基本公式 3 4 , 3 8 - 4 1 如图3 3 所示，n ：a ( x 。、z 。) ，b ( x b 、乙) ，c ( x 。、z e ) 为已知，彳，曰两点的旅行时f 。，乙也 是已知的，射线从么，曰之间的d ( 为、乃) 点通过到达c 点。求d 点的位置和c a n n , j , 旅行时。 以( 矗，z 曩)4 ( ，z 正) 占( 而，乃) a ( ，z 且) d ( 占( 而 乃) d ( x a ，乃) 乙)b ( 而 a j( b ) ( c ) 图3 3 由a 和b 两点的旅行时插值求c 点旅行时的几何关系 设a b 与a c 线段之间的夹角为伊，彳到d 距离为，么，b 之间的d 点的旅行时可根 据彳，b 两点的旅行时t 。，乙通过线性插值可得到，即 气= f 。+ ( 乙- t 。) r i a b i ( 3 1 ) 这时c 点的旅行时f 。可表示为f 。= 屯+ y 。1 ，2 + l a c2 - 2 ，i a c l c o s o ( 3 2 ) 其中v 为介质的波传播速度，i 彳c | ，阻b in i b c l 可由两点距离公式求得，c o s 伊n 在- - n 形a b c 中用余弦定理计算得出。将( 3 1 ) 式代a ( 3 2 ) 式得 1 4 长安大学硕士学位论文 乞：乙+ o 。一乞) ，i a s i + y 一1 ；_ ：下i 苫f i = j 云i 丽 ( 3 3 ) ( 3 3 ) 式将c 点的旅行时表达成了，的函数，根据f e r m a t 原理，初至波的射线从d 点 传到c 点应满足 o t 。o r = 0 即( t b - t , ) l a s + y 一- ( ，一i 彳c l c 。s 驴) ；_ ：f i 苫广二j 云雨= 0 ( 3 4 ) 从方程。4 ，可解出，= a c , c o s o o - 了主崭， 。5 ， 这时，d 点的坐标为 ( 3 6 ) c 点的最小旅行时为 乙= 乙+ i i 留 c o s c , o ( t b - t , ) + s i l l 伊f 五i f 二7 j 而y 】 ( 3 7 ) 需要注意的是，确保按( 3 5 ) 式算出的，是方程( 3 4 ) 的解的条件是t b 乙，所以在根 据已知两点的旅行时按( 3 6 ) 式、( 3 7 ) 计算c 点的最小旅行时* r i d 点的坐标时，应将旅行 时小的点设为a 点，这样才能保证，是彳到插值点d 的距离。还需要注意的是，由于射 线从4 ，b 之间通过，所以必须满足约束条件o ，l a b l 。 以上的旅行时线性插值射线追踪的基本公式既适应于三角形单元或多边形单元划 分的情况，也适应于矩形单元划分的情况。在用矩形单元划分模型时，由于其单元边要 z , 是水平的( 如图3 3 ( b ) 所示) ，要么是垂直的( 如图3 3 ( c ) 所示) ，且c ，d 点都在单元 边上，所以相应的旅行时线性插值射线追踪的基本公式要简单些。 对于如图3 3 ( b ) 所示，a b 在单元水平边界上，即z 。= 乙时，相应的基本公式为 r = c 工c 一工a ，一：