地震层析成像概论.docx

1、地震层析成像概论地震层析成像概论大作业张义蜜，2012260301272016-01-041简述用于地震走时成像方法中的射线追踪算法及原理。1.1打靶法地震射线的打靶法原理简单，它是由射线的运动学方程(1式)定义的初值问题，可用来进行完全射线的追踪计算(存在速度界面时采用斯奈尔定律)。在寻找炮点和检波器点间的地震射线的两点问题是一个反演问题，其中未知参量是射线的初始方向，函数取最小值的原则是使射线终点到检波器间的距离为最小。这类算法的主要挑战是反演问题的非线性行为，且这种非线性行为随介质的复杂程度急剧增加(见图1所示)。（1）图2是打靶法原理图1.1.1近(旁)轴射线追踪射线理论的一个还未提及

2、的重要领域是近轴射线近似，它被广泛地用来进行数据预测。1.1.2完全非线性打靶算法在弱非均匀介质中，非线性迭代打靶算法解决边值问题是十分有效的，但当介质复杂程度增加时则变得不稳定。图3说明了这种情况，展示出具有不同复杂程度的两个速度模型中射线倾角与射线终点到检波器间距离的关系。给定非线性迭代打靶算法易出的错误，完全非线性算法至少是值得研究的。然而，相关文献较少，或许是近来太多的基于网格节点和波前算法出现就是用来克服上述局限。1.2弯曲(调整)法射线追踪调整算法的原理就是不断循环调整连接炮点与检波器间任意初始路径直至成为真实射线路径为止(即：满足费马原理所说的稳定时间路径，见图4所示)。通常的做

3、法是推导出可以循环求解的射线运动学方程的边值形式。1.2.1伪弯曲法伪弯曲法原理上与上述弯曲算法相似，但不对射线方程直接求解。该算法较早提出之一是Um和Thurber(1987)，主要基于射线可由一组线性插值点表示。给定某些任意路径，其目的是不断调整每个节点的位置，从而满足射线方程。利用确定射线路径的法方向，然后直接运用费马原理而有效地实现上述算法。 Zhao等(1992)修改了Um和Thurber(1987)的三点扰动算法，且容许速度界面的存在。这种情况下，定义射线的点序列包括射线穿过每个界面上的点。当更新至界面上这些点时，将沿界面扰动(此时，连接界面两侧的点不动)直到满足斯奈尔定律为止。K

4、oketsu和Sekine(1998)在球坐标系下提出了相似的算法。1.2.2其它弯曲算法除了弯曲和伪弯曲算法外，还有其他类似的弯曲算法也是值得一提的。Prothero等(1988)提出了基于函数极小值的三维弯曲算法。利用大量的搜寻找到炮点与检波器点间最小时间圆弧作为初始射线路径。利用三角谐函数的求和对射线进行扰动，去寻找产生走时最小的振幅系数。Debski和Ando(2004)开发了一种所谓“光谱射线追踪器”，其与Prothero等(1988)算法相近，不同的是射线由Chebyshev多项式定义。弯曲问题可形成函数求极小值问题，其中Chebyshev多项式分项系数变成可调整的变量直至两点间走

5、时为最小。其中不是利用线性化方法去寻求最小值，而是利用遗传基因算法去产生和选择多项式分项系数。射线追踪中大多数求解边值问题可归类为打靶法和弯曲法。然而，其它的算法也是存在的，较著名的是基于结构扰动理论(Cerveny,2001)。在这类算法中，假设在参考介质中已知两点间的射线路径，其目的是要通过对参考介质的扰动确定真实介质中相应两点的射线路径。1.3基于网格(节点)波前扩展的算法炮点与检波器间的射线追踪的另一类选择算法是计算由网格节点定义介质中所有节点上的走时和波前模拟。完整的走时场隐含了波前面是时间即：T(x)的等值线和所有可能路径(由T定义)的函数1.3.1快速行进法(FastMarchi

6、ngMethod)当波前自身相交(也就是说多值路径存在)时，初至波前自身包含绞缠或不连续，通常与后续波前相交。由于通常只计算初至波走时，这些信息被丢弃。为了保证计算的稳定性，计算新走时时，一般不用不连续介质两面不同的差分节点。克服这一问题方法之一是：求解粘性程函方程，其中对不连续进行了平滑。这种平滑解的局限是仅限于初至波波前的计算。业已证明，上述粘性局限解同样可利用熵满足算子求解方程(18)得到，其中估计T时考虑了方向性1.3.2最短路径算法最短路径算法是在网格化节点速度场中计算所有节点走时的另一种流行方法不需求解差分方程，而是由网格节点的连线作为具有走时的射线段。采用类似Dijkstra算法

7、来求取给定点到网格上所有节点的最短路径。根据费马原理两点间最短路径(最小走时路径)相应于实际射线路径。最短路径网格通常是由单元或中心节点来参数化。Nakanishi和Yamaguchi(1986)把速度场参数化为由常速度单元组成，而节点定义在单元边界上(参见图26a)。该算法的优势在于每对节点间的走时可以很容易地计算(t=ds)，这里d是两节点间的距离，s是含节点单元的波慢度。计算精度可通过缩小单元尺度或增加单元边界上节点的个数来提高。图26b刻画了在均匀介质中从某节点到其它几个节点间最小走时路径的选择。另一种形成网格节点的方法是利用规则的速度节点，并进行线性连接(Moser,1991)，如图

8、27a所示。两连接节点A和B间的走时可简单地由：t=d(sA+sB)/2得到，这里sA和sB分别为节点A和B处的波慢度值，而d则为两节点间的距离。图27a中自某一节点到相邻节点连线的角度为45是比较大的，但可以通过增加节点的连接减小角度(图27b)。这种连接方式有时称之为“forwardstar”(KlimesandKvasnicka,1994)。图27a中“forwardstar”有8种连接，而图27b中有16种连接。最短路径网格节点化方式的优势是可以准确地表示连续变化的速度场，同时界面的引入也更容易。一旦网格节点结构及节点间走时计算的方式确定之后，接下来就是计算整个走时场和相应射线路径。D

9、ijkstra(1959)最早提出原始的网格理论算法，其中计算时间正比于O(M2)，而M为所有节点总数。此算法的概念十分简单，即：有总数为Q的未知走时节点，起初Q含M个元素，而P是空集，将Q集节点的走时设置为任意大的数。算法将炮节点加入P集开始，然后由上述提及的“forwardstar”方式计算临近节点的走时(见图28)。这些组成了可能的走时，然后算法从中挑选最小走时，将其加入P集直至所有Q集内节点计算完毕为止。如果Q集内的节点在上次循环中已有计算走时，则选择具有最小走时使之更新。完整的走时场可通过M次迭代得到，射线路径可通过记录节点更新的顺序获得。图27由速度节点参数化的最短路径网格(Mos

10、er,1991)(a:由25个节点组成的网格,其中每个节点最多有8个连接;b:增加节点的连接可以更准确地表示小的路径偏离).尽管最短路径算法只能用来计算连续介质中的初至波路径，但同样也可作合适修改像解程函方程那样去追踪反射和折射波。Moser(1991)描述了在层状介质中利用带约束的最短路径算法计算这类震相地震射线，其中需要射线路径重访界面上特殊的节点。原则上来说，其算法与后面将要讨论的multistageFMM算法相似。最短路径算法在许多实际应用（要求在横向不均匀介质中计算大量的走时信息）中被证明是有效的。在最初的最短路径算法中，Nakanishi和Yamaguchi(1986)利用地方震走

11、时资料反演二维速度结构。Zhang和Toksoz(1998)利用更新的最短路径算法进行二维折射波走时成像，其中采用了波前面上均匀采样而除去网格中不必要的节点。在三维情形，Toomey等(1994)利用与Moser(1991)相似的最短路径算法来反演地壳结构。近来，Bai和Greenhalgh(2005)采用最短路径算法利用区域和地方震走时反演Rabaul(PapuaNewGuinea)火山三维速度结构。1.3.3改进型最短路径算法为了解决射线计算精度与CPU耗时之间的矛盾，Bai等(2004，2007)提出了MSPM算法，该方法与常规最短路径方法不同，通过在网格的主节点之间插入次级节点来增加射

12、线角度覆盖率，满足较大模型的计算(Bai,2004;Baietal,2007)。二维次级节点的情况在图29中给出，只在网格的边上进行插值。图30给出了三维情况下MSPM算法的网格单元模型图。从图中可以看出，一个LxLyLz的三维速度模型，可以划分为CxCyCz的小网格，在每个网格的六个面上插入次级节点(网格内部不进行插值)。为了保证走时计算的精度，可在不增加网格单元的情况下增加次级节点的密度，并能在射线走时计算的同时，给出射线的路径及走时场分布。射线追踪时首先找到炮点所在网格单元，计算此单元节点的走时，找到最小走时点作为次级震源向外进行扩展。计算最小走时的公式为:(72)1.4多次反射与透射波

13、射线追踪1.4.1分区多步快速行进法(MultistageFMM)Li和Ulrych(1993)提供了另一种选择算法来计算二维介质中的反射和折射波走时，其做法是采用局部节点重新划分，即对含有界面的单元重新用三角和矩形单元来较好地拟合界面。其中入射波走时场由Vidale(1988)方法计算，而反射波走时场则采用计算波前自界面走时最小节点重新开始上行波计算的方法得到。Rawlinson和Sambridge(2004b,a)提出的MultistageFMM算法与Li和Ulrych(1993)得算法相似，不同的是可用于更复杂的速度模型，同时可以计算任意多次反射、透射波路径。在计算反射和折射波时有下列两

14、个较难的问题需要解决：(1)准确地描述所含的界面；(2)折射、反射波前的正确传播。在不均匀层状介质中界面随深度变化，因此，界面时常与参数化节点不一致。1.4.2分区多步不规则最短路径算法(MultstageISPM)白超英等(2009)采用分区多步计算技术与上述改进型最短路径(MSPM)算法相结合，实现了模型规则网格划分下二维和三维复杂层状介质中多次波的追踪计算(唐小平、白超英，2009a,b,Baietal,2009)。在此基础上，吸取了MultstageISPM计算多次波时的优点，不同的是在模型参数化时采用一种不规则网格单元进行化分，在保证计算精度的前提下提高了运算速度和计算精度。1.5球

15、坐标系中MultistageISPM算法原理Snoke和Lahr(2001)的研究表明，当震中距大于200km时，则采用直角坐标系进行射线追踪将引入大约0.1s的计算误差，且误差随着震中距的增大几乎是指数型增加的(参见图46)。同样在区域成像研究中，Bijwaard和Spakman(1999)，ZhaoandLei(2004)的研究结果均表明由于采用一维和三维速度模型进行计算，其射线路径偏离近数公里，而走时差有数秒的差别。因此，研发一套可用于实际三维球坐标系下的多震相地震射线追踪的方法技术则显得尤为重要。另一方面上节中有关FMM算法和ISPM算法的对比中可知，无论是计算精度还是CPU时间，ISPM算法均好于FMM算法。鉴于此，有必要将计算高效的MultistageISPM算法推广至球坐标系下，使其能够满足上述要求。这即是本文中正演算法的主要目的之一。当然，这种推广并不是简单的，其中利用了费马的走时稳态路径的概念(见随后讨论)，而在直角坐标系下则仅仅是依据费马的最小走时原理。1.6多值波前(射线)追踪该算法的基本原理认为波前扩展可通过波前上一组点开始进行局部向外传播。对于某一时间步长T，新的波前由老的波前的终点所定义，如图66所示。为了避免波前扩展过程中的采样不足，由某种距离标准可通过插值函数得到相邻射线间新的插值点。因此为了实现该算法，所要求的是给出(38)式的某