1、地震勘探新方法技术地震勘探新方法技术第九讲：地震层析成像技术第九讲：地震层析成像技术一、地震层析成像研究发展概况一、地震层析成像研究发展概况 二、地震层析成像方法面临的主要问题二、地震层析成像方法面临的主要问题 三、地震走时层析成像算法三、地震走时层析成像算法四、实例四、实例地震层析成像技术地震层析成像技术 一、一、地震层析成像研究发展概况地震层析成像研究发展概况 地震层析成像是地球物理学科的一个研究领域。在地震层析成像是地球物理学科的一个研究领域。在地球物理学研究中勘探地球物理是一个年青的学科，地球物理学研究中勘探地球物理是一个年青的学科，它起源于它起源于2020世纪世纪3030年代。早期的
2、地球物理勘探和地年代。早期的地球物理勘探和地球物理方法从属于地质方法，即地质学家预测一个球物理方法从属于地质方法，即地质学家预测一个构造，地球物理学家用他们的原始的勘探技术去验构造，地球物理学家用他们的原始的勘探技术去验证这一结构。证这一结构。6060年代，地球物理学家获得了地下二年代，地球物理学家获得了地下二维数据记录，使得地球物理学家从野外回到了室内，维数据记录，使得地球物理学家从野外回到了室内，为从数据记录中获得地下结构的图像开始了精细的为从数据记录中获得地下结构的图像开始了精细的数据处理与解释研究。数据处理与解释研究。一、一、地震层析成像研究发展概况地震层析成像研究发展概况 8080年
3、代，随着计算机工作站的发展，数据处理技术年代，随着计算机工作站的发展，数据处理技术从二维向三维迅速发展，地球物理学家可以清晰的从二维向三维迅速发展，地球物理学家可以清晰的看清地下结构的图像。从而地球物理学完成了从受看清地下结构的图像。从而地球物理学完成了从受地 质 学 驱 使 到 驱 动 地 质 学 发 展 的 循 环地 质 学 驱 使 到 驱 动 地 质 学 发 展 的 循 环（Russell,2019Russell,2019），即地球物理学家通过使用地球），即地球物理学家通过使用地球物理的数据采集技术、数字处理技术和可视化技术物理的数据采集技术、数字处理技术和可视化技术看清三维地下结构。目
4、前，勘探地球物理已经成为看清三维地下结构。目前，勘探地球物理已经成为经济和生产领域的高技术之一。经济和生产领域的高技术之一。一、一、地震层析成像研究发展概况地震层析成像研究发展概况 2020世纪世纪6060年代初期，美国科学家年代初期，美国科学家CormackCormack从数学和实从数学和实验结果证实了根据验结果证实了根据X X射线的投影可以唯一地确定人体射线的投影可以唯一地确定人体内部结构，从而奠定了医学诊断上图像重建的理论内部结构，从而奠定了医学诊断上图像重建的理论基础，即基础，即X X射线射线CT(X Ray Computer Tomography).CT(X Ray Computer
5、 Tomography).6060年代中期和年代中期和7070年代中期年代中期,随着数学图像重建方法在随着数学图像重建方法在射电天文学和电子显微学方面的应用和发展，在数射电天文学和电子显微学方面的应用和发展，在数学方法上出现了本质上与奥地利数学家学方法上出现了本质上与奥地利数学家19171917年提出年提出的的RndonRndon逆变换方法相同的褶积投影方法，逆变换方法相同的褶积投影方法，ChapmanChapman，19811981）。此后，地学界借助医学）。此后，地学界借助医学CTCT思想，利用地震思想，利用地震波的传播对地壳乃至上地幔结构开始进行半定量研波的传播对地壳乃至上地幔结构开始进
6、行半定量研究。从此，低着层析成像成为地球物理学研究的一究。从此，低着层析成像成为地球物理学研究的一个新领域。个新领域。一、一、地震层析成像研究发展概况地震层析成像研究发展概况 地震层析成像的研究在地震层析成像的研究在7070年代首先以井间速度结构年代首先以井间速度结构调查为研究对象（调查为研究对象（Bois et al.1972Bois et al.1972）。）。19791979年，年，DinesDines和和LytleLytle首先对地震层析成像坐了大量数值模首先对地震层析成像坐了大量数值模拟，并公布了利用弯曲的地震射线进行地下地震波拟，并公布了利用弯曲的地震射线进行地下地震波速 度 成
7、像 的 结 果，并 首 先 将 层 析 成 像速 度 成 像 的 结 果，并 首 先 将 层 析 成 像（Computerized Geophysical TomographyComputerized Geophysical Tomography）这一名）这一名词用于论文的标题。词用于论文的标题。19841984年，美国的年，美国的AndersonAnderson利用利用天然地震数据着手全球构造研究，并公布了全球三天然地震数据着手全球构造研究，并公布了全球三维速度结构。从而使人们对重力场变化、密度结构、维速度结构。从而使人们对重力场变化、密度结构、地幔物质流动有了新的认识。地幔物质流动有了新的
8、认识。一、一、地震层析成像研究发展概况地震层析成像研究发展概况 8080年代，地震层析成像发展到勘探地球物理学界，年代，地震层析成像发展到勘探地球物理学界，自从在亚特兰大（自从在亚特兰大（AtlantaAtlanta）召开的第）召开的第5454届地球物理届地球物理勘探学家协会（勘探学家协会（SEG(Society Of Exploration SEG(Society Of Exploration Geophysicists)Geophysicists)）年会上设置了地震层析成像研究）年会上设置了地震层析成像研究内 容 的 主 题 之 后，以内 容 的 主 题 之 后，以 D a i l yD
9、a i l y（1 9 8 41 9 8 4），），Somersten(1984),Pratt and Worthington Somersten(1984),Pratt and Worthington(1984),Bishop(1985)(1984),Bishop(1985)等人的研究为代表，利用人工等人的研究为代表，利用人工地震发射与接收系统的地震层析成像理论、方法和地震发射与接收系统的地震层析成像理论、方法和技术以数值模拟的形式得到深入、广泛的研究。技术以数值模拟的形式得到深入、广泛的研究。9090年代，不论是利用天然地震数据还是人工地震数据年代，不论是利用天然地震数据还是人工地震数据的
10、地震层析成像方法在认识地球的基础研究领域以的地震层析成像方法在认识地球的基础研究领域以及在资源勘探、工程勘探、环境保护、文物调查、及在资源勘探、工程勘探、环境保护、文物调查、防灾减灾等许多应用领域都得到实验性研究并取得防灾减灾等许多应用领域都得到实验性研究并取得有效的进展。有效的进展。二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.12.1地震波走时自动拾取问题地震波走时自动拾取问题在地震层析成像的研究中在地震层析成像的研究中,可获得的观测数据是地可获得的观测数据是地震记录震记录.从地震记录中可以获得地震波的走时、振从地震记录中可以获得地震波的走时、振幅和频率幅和频率
11、,其中最关键的是地震波走时其中最关键的是地震波走时.随着数字随着数字地震技术的发展地震技术的发展,观测数据的数量迅速增加观测数据的数量迅速增加,准确准确地进行地震波走时的拾取越来越成为一项重要且繁地进行地震波走时的拾取越来越成为一项重要且繁重的工作重的工作.为此为此,走时的自动拾取成为人们研究与走时的自动拾取成为人们研究与关注的对象关注的对象.二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.12.1地震波走时自动拾取问题地震波走时自动拾取问题 近 年 来近 年 来 ,先 后 出 现 了 相 邻 道 互 相 关 方先 后 出 现 了 相 邻 道 互 相 关 方法法(Ge
12、lchinsky(Gelchinsky，1983)1983)、能量比较方法、能量比较方法(Coppens(Coppens，1985)1985)、改变褶积算子宽度的方法（、改变褶积算子宽度的方法（Ramananant,1 Ramananant,1 987987）等等）等等.90.90年代之后年代之后,地震初至波走时拾取的方地震初至波走时拾取的方法有了新的进展法有了新的进展,如基于分形理论的如基于分形理论的DividerDivider方法和方法和HurstHurst方法方法(Boschetti,1 996),(Boschetti,1 996),特别是中国科学家特别是中国科学家提出的基于提出的基于H
13、ausdorffHausdorff分维算法的地震波走时全自动分维算法的地震波走时全自动拾取方法拾取方法(Chang,1 999)(Chang,1 999)也取得了良好的效果也取得了良好的效果.二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.12.1地震波走时自动拾取问题地震波走时自动拾取问题这些新方法以分形这些新方法以分形(Fractal)(Fractal)理论为依据理论为依据,通过对通过对地震记录时间序列分数维地震记录时间序列分数维(Fractaldimension)(Fractaldimension)的计的计算算,实现了对地震波初至走时的自动拾取实现了对地震波初至
14、走时的自动拾取.分数维分数维方法的最突出的优点是相邻道的无关方法的最突出的优点是相邻道的无关,这一优点显这一优点显然适用于地震层析成像方法以及天然地震的无规则然适用于地震层析成像方法以及天然地震的无规则观测方式观测方式,是一项非常有实用价值的方法是一项非常有实用价值的方法.从文献从文献中获悉中获悉,利用分形分维理论利用分形分维理论,根据地震波初至到达根据地震波初至到达前后地震记录分维的差别前后地震记录分维的差别,对反射地震记录、透射对反射地震记录、透射地震记录以及天然地震记录可以成功地进行初至走地震记录以及天然地震记录可以成功地进行初至走时的拾取时的拾取.目前这一方法得到更多专业人员的关注目前
15、这一方法得到更多专业人员的关注和研究和研究,可望获得更有效的研究成果可望获得更有效的研究成果.二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.2 2.2 三维波动方程有限差分算法模拟地震波场的问三维波动方程有限差分算法模拟地震波场的问题题不论是天然地震还是人工地震不论是天然地震还是人工地震(即使是二维观测方即使是二维观测方式式)的观测数据都是在三维空间介质中形成的观测数据都是在三维空间介质中形成.由于由于地下地质结构的千变万化地下地质结构的千变万化,理论数据的正演计算只理论数据的正演计算只有在三维空间中实现才更具有实际意义有在三维空间中实现才更具有实际意义.而目前大而
16、目前大多采用二维计算多采用二维计算,使得理论数据与观测数据之间的使得理论数据与观测数据之间的误差不仅由地质模型形成而且还由计算方法的数学误差不仅由地质模型形成而且还由计算方法的数学模型形成模型形成.三维波动方程的有限差分解是获取地震三维波动方程的有限差分解是获取地震波三维波场的有效方法波三维波场的有效方法.二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.2 2.2 三维波动方程有限差分算法模拟地震波场的问三维波动方程有限差分算法模拟地震波场的问题题 开展非弹性介质和完全弹性介质有限差分法三开展非弹性介质和完全弹性介质有限差分法三维地震波场的算法研究维地震波场的算法研究
17、,采用多重网格算法求解差采用多重网格算法求解差分方程分方程,用粗网格的低频特性与细网格高频特性的用粗网格的低频特性与细网格高频特性的互补提高差分方程的求解速度和精度互补提高差分方程的求解速度和精度.通过两种介通过两种介质三维偏微分方程解的定量化对比质三维偏微分方程解的定量化对比,可得到由介质可得到由介质粘滞性引起的地震波振幅的衰减和频率的变化粘滞性引起的地震波振幅的衰减和频率的变化.从从而提供地震波衰减特性层析成像方法中的基础数而提供地震波衰减特性层析成像方法中的基础数据据理论地震波场理论地震波场.这一问题的研究可为多分量这一问题的研究可为多分量地震层析成像方法奠定新的生长点地震层析成像方法奠
18、定新的生长点.二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.3 2.3 三维程函方程有限差分算法三维程函方程有限差分算法 模拟地震时间场的问题模拟地震时间场的问题三维程函方程的有限差分解是获取地震波三维空间走三维程函方程的有限差分解是获取地震波三维空间走时的有效方法时的有效方法.但由于偏微分方程中走时对于空间位但由于偏微分方程中走时对于空间位置的二阶导数在地质模型的突变点不连续置的二阶导数在地质模型的突变点不连续,使得复杂使得复杂地质模型三维时间场的计算出现畸变值地质模型三维时间场的计算出现畸变值.开展有限差开展有限差分法三维程函方程突变型地质模型时间场的算法研分法
19、三维程函方程突变型地质模型时间场的算法研究究.对于突变形地质模型对于突变形地质模型,用突变点线性震源方法以用突变点线性震源方法以及多重网格算法进行计算及多重网格算法进行计算.二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.3 2.3 三维程函方程有限差分算法三维程函方程有限差分算法 模拟地震时间场的问题模拟地震时间场的问题 这种方法在二阶导数的间断点考虑新的程函方这种方法在二阶导数的间断点考虑新的程函方程程,可分别计算上行波和下行波可分别计算上行波和下行波(反射波和透射反射波和透射波波),),可望获得复杂地质模型的三维时间场函数可望获得复杂地质模型的三维时间场函数 (
20、Mitcelletal,1980;Vidale,1988).(Mitcelletal,1980;Vidale,1988).这一方法的实现这一方法的实现将比以往更为客观地提供地震层析成像方法中的重要将比以往更为客观地提供地震层析成像方法中的重要基础数据基础数据理论走时理论走时.此项研究可能取代传统的几此项研究可能取代传统的几何学射线追踪方法何学射线追踪方法,成为地震层析成像研究中的一项成为地震层析成像研究中的一项至关重要的技术至关重要的技术二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.42.4地震反演解的可靠性问题地震反演解的可靠性问题由于震源和检波器位置分布及连续问
21、题的离散化由于震源和检波器位置分布及连续问题的离散化,地震层析成像反演将遇到方程的不适定问题地震层析成像反演将遇到方程的不适定问题.若方若方程组是欠定的程组是欠定的,解可能不存在解可能不存在,或者没有唯一解或者没有唯一解.当条件数很大时当条件数很大时,反演问题将是不稳定的反演问题将是不稳定的,所用算所用算法也可能不稳定法也可能不稳定.若方程是超定的若方程是超定的,说明方程组中说明方程组中的一个或几个方程是其它方程的线性组合的一个或几个方程是其它方程的线性组合,或者所或者所有的方程中某些变量是其它变量的同一线性组合有的方程中某些变量是其它变量的同一线性组合,这两种情况都得不出唯一解这两种情况都得
22、不出唯一解(杨文采杨文采,1 997;,1 997;刘福刘福田田,1989).,1989).二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.42.4地震反演解的可靠性问题地震反演解的可靠性问题尽管有些方程并非为彼此精确的线性组合尽管有些方程并非为彼此精确的线性组合,但它们但它们可能很接近线性依赖可能很接近线性依赖,若在求解过程中若在求解过程中,由于机器由于机器的舍入误差使它们成为线性依赖的舍入误差使它们成为线性依赖,计算将会失败计算将会失败.求解过程舍入误差的积累也会使结果和真实解相差求解过程舍入误差的积累也会使结果和真实解相差甚远甚远,尤其在解的数目很多时特别容易发
23、生尤其在解的数目很多时特别容易发生,而计而计算程序在方法上并没有错误算程序在方法上并没有错误.由于地形的限制由于地形的限制,观观测系统布置的不均匀引起射线分布的不均匀测系统布置的不均匀引起射线分布的不均匀,反演反演的速度结果产生伪图像的速度结果产生伪图像.因此应用一种适应于不适因此应用一种适应于不适定问题的图像重建方法是确保得出地球物理真实解定问题的图像重建方法是确保得出地球物理真实解的重要前提的重要前提.地球物理的反演问题总是不适定的地球物理的反演问题总是不适定的,一般不存在经典意义下的解一般不存在经典意义下的解,只能给出某种意义下只能给出某种意义下的广义解的广义解.对非唯一性条件下求出的解
24、对非唯一性条件下求出的解,给出解的给出解的评价显得尤其重要评价显得尤其重要 。二、二、地震层析成像方法面临的主要问题地震层析成像方法面临的主要问题 2.5 2.5 地质解释问题地质解释问题地震层析成像方法的最终目的是要对成像结果做出地震层析成像方法的最终目的是要对成像结果做出切合实际的地质解释切合实际的地质解释.成像结果中包含各类信息成像结果中包含各类信息,有些与地质实际有关有些与地质实际有关,有些则无关有些则无关.须对地震波动须对地震波动在研究区域的覆盖程度与图像重建结果之间的内在在研究区域的覆盖程度与图像重建结果之间的内在联系进行分析联系进行分析,图像重建结果中的假象不容忽视图像重建结果中
25、的假象不容忽视,同时应该研究重建图像的地质解释方法同时应该研究重建图像的地质解释方法.这将是地这将是地震层析成像方法真正地切实地用于解决实际问题时震层析成像方法真正地切实地用于解决实际问题时必须认真对待的一个极其重要的方面必须认真对待的一个极其重要的方面 。三、地震走时层析成像算法三、地震走时层析成像算法 地震层析成像技术大致可以分为地震层析成像技术大致可以分为2 2种类型：一种是种类型：一种是基于射线理论的图像重建技术，包括地震走时层析基于射线理论的图像重建技术，包括地震走时层析成祥和地震衰减层析成像；另一种是基于波动方程成祥和地震衰减层析成像；另一种是基于波动方程反演的散射（或衍射）层析成
26、像。从数学角度来看，反演的散射（或衍射）层析成像。从数学角度来看，后者是有一个函数的线性积分反求这个函数的问题，后者是有一个函数的线性积分反求这个函数的问题，当射线是直线时，这种方法比较成熟。当射线是直线时，这种方法比较成熟。三、地震走时层析成像算法三、地震走时层析成像算法3 31 1离散图像重建离散图像重建设设 f(x,y)f(x,y)在 区 域在 区 域 外 恒 等 于 零，即外 恒 等 于 零，即 f(x,y)f(x,y)0,(x,y),0,(x,y),将区域分割成将区域分割成I I个不重叠小区域个不重叠小区域 （像（像元），元），I I1 1，2 2，3 3，,I,I。f(x,y)f(
27、x,y)在在 上各点之值上各点之值用它在上之平均值近似代替，即：用它在上之平均值近似代替，即：其中其中 表示小区域表示小区域 之面积，之面积，向量向量 被称为图像向量。被称为图像向量。iSiSiSiSiSdxdyyxfif),(iSTIffff,.,21三、地震走时层析成像算法三、地震走时层析成像算法3 31 1离散图像重建离散图像重建设射线设射线 与小区域与小区域 相交部分之长度为相交部分之长度为 ，根据，根据RadonRadon变换，函数变换，函数 沿射线沿射线 的投影函数为的投影函数为 (1)(1)其中其中j=1,2,J,j=1,2,J,为误差项。略去误差项为误差项。略去误差项 ，则可得
28、离散图像重建的线性方程组为则可得离散图像重建的线性方程组为 (2)(2)其中其中J J为射线总数。为射线总数。iLiSjiafiLIijijiLjefadsyxffj1),(jejeIiijijJjfaf1),.,2,1(三、地震走时层析成像算法三、地震走时层析成像算法在地震走时层析成像情况下，投影数据在地震走时层析成像情况下，投影数据 为地震为地震波走时波走时 ，图像向量，图像向量f f为像元内慢度的平均值，则为像元内慢度的平均值，则式（式（2 2）可写成矩阵方程）可写成矩阵方程 （3 3）式中，式中，A A为一个（为一个（I IxJxJ）的）的 值矩阵，其中值矩阵，其中J J为穿过为穿过要
29、讨论的区域的全部射线数，要讨论的区域的全部射线数，I I要讨论的区域的全部要讨论的区域的全部单元数。单元数。A A是一个相对松散的矩阵，因为任何一条射是一个相对松散的矩阵，因为任何一条射线通常只会穿过研究区中少部分单元。线通常只会穿过研究区中少部分单元。离散图像重建问题转化为给出一系列地震波走离散图像重建问题转化为给出一系列地震波走时时，计算图像向量即地质体介质慢度向量，计算图像向量即地质体介质慢度向量f，只要矩，只要矩阵阵A建立了，求它的逆建立了，求它的逆A-1，则矢量，则矢量f就可以很容易地就可以很容易地求出。求出。fjjA fjia三、地震走时层析成像算法三、地震走时层析成像算法3 32
30、 2重建算法重建算法 方程组（方程组（3 3）中的系数矩阵）中的系数矩阵A A是极其稀疏的，因为它是极其稀疏的，因为它的每一行有的每一行有J J各元素，而每条地震波只通过所有各元素，而每条地震波只通过所有I I个个像元中的一小部分，因此矩阵像元中的一小部分，因此矩阵A A中的大部分元素为零。中的大部分元素为零。根据系数矩阵稀疏的特点，对方程组（根据系数矩阵稀疏的特点，对方程组（3 3）多采用迭）多采用迭代方法求解。代方法求解。三、地震走时层析成像算法三、地震走时层析成像算法 3.2.1 BPT3.2.1 BPT算法（算法（Back projection techniqueBack projec
31、tion technique）求解方程组（求解方程组（3 3）的最简单和粗糙的离散图像重构方）的最简单和粗糙的离散图像重构方法之一是反投影法。将走时沿射线分配给每一个像法之一是反投影法。将走时沿射线分配给每一个像元，分配时以第元，分配时以第j j条射线在像元内的长度条射线在像元内的长度 与射线与射线总长度之比为权，然后把通过总长度之比为权，然后把通过i i像元在加权后的走时像元在加权后的走时对所有射线相加，并除以单元内总射线长度求得该对所有射线相加，并除以单元内总射线长度求得该单元的介质慢度。单元的介质慢度。(4)(4)式中：式中：I I为反演区域像元总数；为反演区域像元总数；J J为地震波射
32、线总数。为地震波射线总数。Iijia1JjjiIijijIijiiaaaf111/三、地震走时层析成像算法三、地震走时层析成像算法 3.2.2 ART3.2.2 ART算法（算法（Algebraic Reconstruction Algebraic Reconstruction TechniqueTechnique）代数重建是按射线依次修改有关像元的图像向量的代数重建是按射线依次修改有关像元的图像向量的一类迭代算法。在方程（一类迭代算法。在方程（2 2）中令图像向量产生一增）中令图像向量产生一增量量 ，有，有 (5)(5)作为迭代算法要根据第作为迭代算法要根据第j j条射线的走时差条射线的走时
33、差 求慢度求慢度的修改增量的修改增量 。由于方程（。由于方程（5 5）可能是欠定的或病）可能是欠定的或病态的，可以用它作为约束求态的，可以用它作为约束求 的的L L2 2模的极小解。模的极小解。f),.,2,1(1JjafIijiijjifif三、地震走时层析成像算法三、地震走时层析成像算法由拉格朗日乘子法令目标函数为由拉格朗日乘子法令目标函数为（其中（其中为拉格朗日乘子），由为拉格朗日乘子），由 得得 ，代入式（，代入式（5 5）有）有 （6 6）（7 7）因此，再有式（因此，再有式（6 6）便可写出对第）便可写出对第j j条射线及第条射线及第i i个像个像元求波慢修改增量得公式。元求波慢修
34、改增量得公式。IijjiiiaffQ12min)(0)(/ifQjiiaf2Iijija122Iijija12/2三、地震走时层析成像算法三、地震走时层析成像算法当然，不一定非取当然，不一定非取 的的L L2 2模极小不可，也可以取任意模极小不可，也可以取任意阶的模，如阶的模，如L L2p2p模极小来求慢度的修改增量。其中模极小来求慢度的修改增量。其中P P1 1，2 2，。但是当。但是当P1P1时涉及开方运算，速度太慢，一时涉及开方运算，速度太慢，一般很少采。只有在令般很少采。只有在令 时，可导出时，可导出ARTART迭代的最迭代的最简单修正公式简单修正公式 （8 8）ifpIijijiaf
35、1/三、地震走时层析成像算法三、地震走时层析成像算法 上式说明走时差平均地分配给每一条射线上式说明走时差平均地分配给每一条射线j j通过地单通过地单元，而不考虑像元内射线地长短。由式（元，而不考虑像元内射线地长短。由式（6 6）和式（）和式（8 8）可得可得ARTART方法迭代公式方法迭代公式 jiIijikjjkikiiaafafff12)()()1()0(),(初始任取Iijikjjkikiiafafff1)()()1()0(),(初始任取式中：式中：为系数矩阵为系数矩阵的分量；的分量；k k=0,1,=0,1,为迭代次数。当为迭代次数。当时，取时，取 ，其，其中中为给定允许误为给定允许误
36、差差。jaTjaaaA),.,(21)()1(kkffffk)1(三、地震走时层析成像算法三、地震走时层析成像算法3 3.2.3 SIRT2.3 SIRT算法（算法（simultaneous iterative simultaneous iterative reconstruction techniquereconstruction technique）SIRTSIRT算法就是联合迭代重建技术，用一个或多个优化算法就是联合迭代重建技术，用一个或多个优化准则使得的解估计唯一。按照二次最优化准则，得到准则使得的解估计唯一。按照二次最优化准则，得到SIRTSIRT典型的迭代修正公式典型的迭代修正公式
37、 （1111）式中：式中：，为非负实数；为非负实数；u u为松弛系数；为松弛系数；W W为对为对称矩阵；称矩阵；A A为系数矩阵；为系数矩阵；B B为非负定矩阵；为非负定矩阵；C C为对称正为对称正定矩阵；定矩阵；B B和和C C又称为平滑矩阵；又称为平滑矩阵；I I为单位矩阵；为单位矩阵；为为与与 同维的向量，它与图像数值的某种先验知识同维的向量，它与图像数值的某种先验知识有关。有关。)()(1fffuffkokTkkkICBAWCAf0fk三、地震走时层析成像算法三、地震走时层析成像算法矩阵矩阵W,B,CW,B,C和向量和向量 要根据实际问题来选择，如对最要根据实际问题来选择，如对最小二乘
38、法和最小范数解估计，选小二乘法和最小范数解估计，选 ，B=I,B=I,及及 ,可作为初始猜测（可作为初始猜测（为数据协方差矩阵为数据协方差矩阵的逆、的逆、为为图像向量的协方差矩阵）。图像向量的协方差矩阵）。f0CW 1CCff0C1Cf三、地震走时层析成像算法三、地震走时层析成像算法可以看出可以看出，SIRTSIRT重建算法的特点是在某一轮迭代中，重建算法的特点是在某一轮迭代中，所有像元的值（即图像向量所有像元的值（即图像向量 ）都用前一轮的迭代）都用前一轮的迭代结果结果 来修正。为了得到最简单情况下的来修正。为了得到最简单情况下的SIRTSIRT迭代迭代修正公式，可令式（修正公式，可令式（1111）中）中 0 0，适当选取矩，适当选取矩阵阵W W，C C和重新定义不随迭代变化的松弛系数和重新定义不随迭代变化的松弛系数u u，则式，则式（1111）可写为）可写为 fk)1(fk)(JjjiJjJjjikjjjikikiaafaafff111)()()1(0/),(初始任取四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例四、实例
