(第九讲)地震层析成象

地震勘探新方法技术第九讲：地震层析成像技术一、地震层析成像研究发展概况二、地震层析成像方法面临的主要问题三、地震走时层析成像算法四、实例地震层析成像技术一、地震层析成像研究发展概况地震层析成像是地球物理学科的一个研究领域。在地球物理学研究中勘探地球物理是一个年青的学科，它起源于20世纪30年代。早期的地球物理勘探和地球物理方法从属于地质方法，即地质学家预测一个构造，地球物理学家用他们的原始的勘探技术去验证这一结构。60年代，地球物理学家获得了地下二维数据记录，使得地球物理学家从野外回到了室内，为从数据记录中获得地下结构的图像开始了精细的数据处理与解释研究。一、地震层析成像研究发展概况80年代，随着计算机工作站的发展，数据处理技术从二维向三维迅速发展，地球物理学家可以清晰的看清地下结构的图像。从而地球物理学完成了从受地质学驱使到驱动地质学发展的循环（Russell,1999），即地球物理学家通过使用地球物理的数据采集技术、数字处理技术和可视化技术看清三维地下结构。目前，勘探地球物理已经成为经济和生产领域的高技术之一。一、地震层析成像研究发展概况20世纪60年代初期，美国科学家Cormack从数学和实验结果证实了根据X射线的投影可以唯一地确定人体内部结构，从而奠定了医学诊断上图像重建的理论基础，即X射线CT(XRayComputerTomography).60年代中期和70年代中期,随着数学图像重建方法在射电天文学和电子显微学方面的应用和发展，在数学方法上出现了本质上与奥地利数学家1917年提出的Rndon逆变换方法相同的褶积投影方法，Chapman，1981）。此后，地学界借助医学CT思想，利用地震波的传播对地壳乃至上地幔结构开始进行半定量研究。从此，低着层析成像成为地球物理学研究的一个新领域。一、地震层析成像研究发展概况地震层析成像的研究在70年代首先以井间速度结构调查为研究对象（Boisetal.1972）。1979年，Dines和Lytle首先对地震层析成像坐了大量数值模拟，并公布了利用弯曲的地震射线进行地下地震波速度成像的结果，并首先将层析成像（ComputerizedGeophysicalTomography）这一名词用于论文的标题。1984年，美国的Anderson利用天然地震数据着手全球构造研究，并公布了全球三维速度结构。从而使人们对重力场变化、密度结构、地幔物质流动有了新的认识。一、地震层析成像研究发展概况80年代，地震层析成像发展到勘探地球物理学界，自从在亚特兰大（Atlanta）召开的第54届地球物理勘探学家协会（SEG(SocietyOfExplorationGeophysicists)）年会上设置了地震层析成像研究内容的主题之后，以Daily（1984），Somersten(1984),PrattandWorthington(1984),Bishop(1985)等人的研究为代表，利用人工地震发射与接收系统的地震层析成像理论、方法和技术以数值模拟的形式得到深入、广泛的研究。90年代，不论是利用天然地震数据还是人工地震数据的地震层析成像方法在认识地球的基础研究领域以及在资源勘探、工程勘探、环境保护、文物调查、防灾减灾等许多应用领域都得到实验性研究并取得有效的进展。二、地震层析成像方法面临的主要问题2.1地震波走时自动拾取问题在地震层析成像的研究中,可获得的观测数据是地震记录.从地震记录中可以获得地震波的走时、振幅和频率,其中最关键的是地震波走时.随着数字地震技术的发展,观测数据的数量迅速增加,准确地进行地震波走时的拾取越来越成为一项重要且繁重的工作.为此,走时的自动拾取成为人们研究与关注的对象.二、地震层析成像方法面临的主要问题2.1地震波走时自动拾取问题近年来,先后出现了相邻道互相关方法](Gelchinsky，1983)、能量比较方法(Coppens，1985)、改变褶积算子宽度的方法（Ramananant,1987）等等.90年代之后,地震初至波走时拾取的方法有了新的进展,如基于分形理论的Divider方法和Hurst方法(Boschetti,1996),特别是中国科学家提出的基于Hausdorff分维算法的地震波走时全自动拾取方法(Chang,1999)也取得了良好的效果.二、地震层析成像方法面临的主要问题2.1地震波走时自动拾取问题这些新方法以分形(Fractal)理论为依据,通过对地震记录时间序列分数维(Fractaldimension)的计算,实现了对地震波初至走时的自动拾取.分数维方法的最突出的优点是相邻道的无关,这一优点显然适用于地震层析成像方法以及天然地震的无规则观测方式,是一项非常有实用价值的方法.从文献中获悉,利用分形分维理论,根据地震波初至到达前后地震记录分维的差别,对反射地震记录、透射地震记录以及天然地震记录可以成功地进行初至走时的拾取.目前这一方法得到更多专业人员的关注和研究,可望获得更有效的研究成果.二、地震层析成像方法面临的主要问题2.2三维波动方程有限差分算法模拟地震波场的问题不论是天然地震还是人工地震(即使是二维观测方式)的观测数据都是在三维空间介质中形成.由于地下地质结构的千变万化,理论数据的正演计算只有在三维空间中实现才更具有实际意义.而目前大多采用二维计算,使得理论数据与观测数据之间的误差不仅由地质模型形成而且还由计算方法的数学模型形成.三维波动方程的有限差分解是获取地震波三维波场的有效方法.二、地震层析成像方法面临的主要问题2.2三维波动方程有限差分算法模拟地震波场的问题开展非弹性介质和完全弹性介质有限差分法三维地震波场的算法研究,采用多重网格算法求解差分方程,用粗网格的低频特性与细网格高频特性的互补提高差分方程的求解速度和精度.通过两种介质三维偏微分方程解的定量化对比,可得到由介质粘滞性引起的地震波振幅的衰减和频率的变化.从而提供地震波衰减特性层析成像方法中的基础数据——理论地震波场.这一问题的研究可为多分量地震层析成像方法奠定新的生长点.二、地震层析成像方法面临的主要问题2.3三维程函方程有限差分算法模拟地震时间场的问题三维程函方程的有限差分解是获取地震波三维空间走时的有效方法.但由于偏微分方程中走时对于空间位置的二阶导数在地质模型的突变点不连续,使得复杂地质模型三维时间场的计算出现畸变值.开展有限差分法三维程函方程突变型地质模型时间场的算法研究.对于突变形地质模型,用突变点线性震源方法以及多重网格算法进行计算.二、地震层析成像方法面临的主要问题2.3三维程函方程有限差分算法模拟地震时间场的问题这种方法在二阶导数的间断点考虑新的程函方程,可分别计算上行波和下行波(反射波和透射波),可望获得复杂地质模型的三维时间场函数(Mitcelletal,1980;Vidale,1988).这一方法的实现将比以往更为客观地提供地震层析成像方法中的重要基础数据——理论走时.此项研究可能取代传统的几何学射线追踪方法,成为地震层析成像研究中的一项至关重要的技术二、地震层析成像方法面临的主要问题2.4地震反演解的可靠性问题由于震源和检波器位置分布及连续问题的离散化,地震层析成像反演将遇到方程的不适定问题.若方程组是欠定的,解可能不存在,或者没有唯一解.当条件数很大时,反演问题将是不稳定的,所用算法也可能不稳定.若方程是超定的,说明方程组中的一个或几个方程是其它方程的线性组合,或者所有的方程中某些变量是其它变量的同一线性组合,这两种情况都得不出唯一解(杨文采,1997;刘福田,1989).二、地震层析成像方法面临的主要问题2.4地震反演解的可靠性问题尽管有些方程并非为彼此精确的线性组合,但它们可能很接近线性依赖,若在求解过程中,由于机器的舍入误差使它们成为线性依赖,计算将会失败.求解过程舍入误差的积累也会使结果和真实解相差甚远,尤其在解的数目很多时特别容易发生,而计算程序在方法上并没有错误.由于地形的限制,观测系统布置的不均匀引起射线分布的不均匀,反演的速度结果产生伪图像.因此应用一种适应于不适定问题的图像重建方法是确保得出地球物理真实解的重要前提.地球物理的反演问题总是不适定的,一般不存在经典意义下的解,只能给出某种意义下的广义解.对非唯一性条件下求出的解,给出解的评价显得尤其重要。二、地震层析成像方法面临的主要问题2.5地质解释问题地震层析成像方法的最终目的是要对成像结果做出切合实际的地质解释.成像结果中包含各类信息,有些与地质实际有关,有些则无关.须对地震波动在研究区域的覆盖程度与图像重建结果之间的内在联系进行分析,图像重建结果中的假象不容忽视,同时应该研究重建图像的地质解释方法.这将是地震层析成像方法真正地切实地用于解决实际问题时必须认真对待的一个极其重要的方面。三、地震走时层析成像算法地震层析成像技术大致可以分为2种类型：一种是基于射线理论的图像重建技术，包括地震走时层析成祥和地震衰减层析成像；另一种是基于波动方程反演的散射（或衍射）层析成像。从数学角度来看，后者是有一个函数的线性积分反求这个函数的问题，当射线是直线时，这种方法比较成熟。三、地震走时层析成像算法3．1离散图像重建设f(x,y)在区域外恒等于零，即f(x,y)0,(x,y),将区域分割成I个不重叠小区域（像元），I＝1，2，3，…,I。f(x,y)在上各点之值用它在上之平均值近似代替，即：其中表示小区域之面积，向量被称为图像向量。iSiSiSiSiSdxdyyxfif),(iSTIffff,...,,21三、地震走时层析成像算法3．1离散图像重建设射线与小区域相交部分之长度为，根据Radon变换，函数沿射线的投影函数为(1)其中j=1,2,…,J,为误差项。略去误差项，则可得离散图像重建的线性方程组为(2)其中J为射线总数。iLiSjiafiLIijijiLjefadsyxffj1),(jejeIiijijJjfaf1),...,2,1(三、地震走时层析成像算法在地震走时层析成像情况下，投影数据为地震波走时，图像向量f为像元内慢度的平均值，则式（2）可写成矩阵方程（3）式中，A为一个（IxJ）的值矩阵，其中J为穿过要讨论的区域的全部射线数，I要讨论的区域的全部单元数。A是一个相对松散的矩阵，因为任何一条射线通常只会穿过研究区中少部分单元。离散图像重建问题转化为给出一系列地震波走时τ，计算图像向量即地质体介质慢度向量f，只要矩阵A建立了，求它的逆A-1，则矢量f就可以很容易地求出。fjjτAfjia三、地震走时层析成像算法3．2重建算法方程组（3）中的系数矩阵A是极其稀疏的，因为它的每一行有J各元素，而每条地震波只通过所有I个像元中的一小部分，因此矩阵A中的大部分元素为零。根据系数矩阵稀疏的特点，对方程组（3）多采用迭代方法求解。三、地震走时层析成像算法3.2.1BPT算法（Backprojectiontechnique）求解方程组（3）的最简单和粗糙的离散图像重构方法之一是反投影法。将走时沿射线分配给每一个像元，分配时以第j条射线在像元内的长度与射线总长度之比为权，然后把通过i像元在加权后的走时对所有射线相加，并除以单元内总射线长度求得该单元的介质慢度。(4)式中：I为反演区域像元总数；J为地震波射线总数。Iijia1JjjiIijijIijiiaaaf111//~三、地震走时层析成像算法3.2.2ART算法（AlgebraicReconstructionTechnique）代数重建是按射线依次修改有关像元的图像向量的一类迭代算法。在方程（2）中令图像向量产生一增量，有(5)作为迭代算法要根据第j条射线的走时差求慢度的修改增量。由于方程（5）可能是欠定的或病态的，可以用它作为约束求的L2模的极小解。f),...,2,1(1JjafIijiijjifif三、地震走时层析成像算法由拉格朗日乘子法令目标函数为（其中λ为拉格朗日乘子），由得，代入式（5）有（6）（7）因此，再有式（6）便可写