最短路径法射线踪MATLAB实现.docVIP

最短路径法射线追踪的MATLAB实现 李志辉 刘争平 （西南交通大学土木工程学院 成都 610031） 摘 要：本文探讨了在MATLAB环境中实现最短路径射线追踪的方法和步骤，并通过数值模拟演示了所编程序在射线追踪正演计算中的应用。 关键词：最短路径法 射线追踪 MATLAB 数值模拟 利用地震初至波确定近地表介质结构，在矿产资源的勘探开发及工程建设中有重要作用。地震射线追踪方法是研究地震波传播的有效工具，目前常用的方法主要有有限差分解程函方程法和最小路径法。最短路径方法起源于网络理论，首次由Nakanishi和Yamaguchi应用域地震射线追踪中。Moser以及Klimes和Kvasnicha对最短路径方法进行了详细研究。通过科技人员的不断研究，最短路径方法目前已发展较为成熟，其基本算法的计算程序也较为固定。 被称作是第四代计算机语言的MATLAB语言，利用其丰富的函数资源把编程人员从繁琐的程序代码中解放出来。MATLAB用更直观的、符合人们思维习惯的代码，为用户提供了直观、简洁的程序开发环境。本文介绍运用Matlab实现最短路径法的方法和步骤，便于科研院校教学中讲授、演示和理解最短路径方法及其应用。 1 最短路径法射线追踪方法原理 最短路径法的基础是Fermat原理及图论中的最短路径理论。其基本思路是，对实际介质进行离散化，将这个介质剖分成一系列小单元，在单元边界上设置若干节点，并将彼此向量的节点相连构成一个网络。网络中，速度场分布在离散的节点上。相邻节点之间的旅行时为他们之间欧氏距离与其平均慢度之积。将波阵面看成式由有限个离散点次级源组成，对于某个次级源（即某个网格节点），选取与其所有相邻的点（邻域点）组成计算网格点；由一个源点出发，计算出从源点到计算网格点的透射走时、射线路径、和射线长度；然后把除震源之外的所有网格点相继当作次级源，选取该节点相应的计算网格点，计算出从次级源点到计算网格点的透射走时、射线路径、和射线长度；将每次计算出来的走时加上从震源到次级源的走时，作为震源点到该网格节点的走时，记录下相应的射线路径位置及射线长度。 图1 离散化模型（星点表示震源或次级震源，空心点为对应计算网格点） 根据Fermat原理逐步计算最小走时及射线方向。设Ω为已知走时点q的集合，p为与其相邻的未知走时点，tq分别和p点的最小走时，tqp为q至p最小走时。r为p的次级源位置，则 根据Huygens原理，q只需遍历Q的边界（即波前点），当所有波前邻点的最小走时都求出时，这些点又成为新的波前点。应用网络理论中的最短路径算法，可以同时求出从震源点传至所有节点之间的连线近似地震射线路径。 2 最短路径法射线追踪基本算法步骤 把网格上的所有节点分成集合p和q，p为已知最小旅行时的结点总数集合，q为未知最小旅行时的节点的集合。若节点总数为n，经过n次迭代后可为求出所有节点的最小旅行时。过程如下： 初始时 q集合包含所有节点，除震源s的旅行时已知为ts＝0外，其余所有节点的旅行时均为ti＝（i属于Q但不等于s）。P集合为空集。 在Q中找一个旅行时最小的节点i，它的旅行时为ti； 确定与节点i相连的所有节点的集合V； 求节点j（j属于V且j不属于P）与节点i连线的旅行时dtij； 求节点j（）的新旅行时tj（取原有旅行时tj与tj＋dtij的最小值）； 将i点从Q集合转到P集合； 若P集合中的节点个数小于总节点数N，转2，否则结束旅行时追踪； 从接收点开始倒推出各道从源点道接收点的射线路径，只要每个节点记下使它形成最小旅行时的前一个节点号，就很容易倒推出射线路径。 Matlab语言可以十分方便地构造用户自己的函数（可以同其它目录里的函数同名），供主函数调用，并且通过Matlab结构形式，根据属性名将节点上不同类型的数据组织起来，从而简洁地实现算法中有关节点的判断、调用和运算。 3 数值模拟 3.1 构造数值模型 图2 探测区域数值模型（空白区为低速介质区，黑色区为高速介质区） 如图2所示，数值模拟探测区域长Length=20m，宽Width=8m。假设已知背景介质速度为V低，存在三个黑色区高速介质区，其中区域1，2相对于水平中线对称，区域3相对水平居中，V高＝4V低。可离散化为间距为0.5m的17行，间距为1m的17行21列的网络。 3.2私有函数构造 （1）构造函数Vnew_jishuan（），用以确定子波源点所在的结点相连的所有结点集合V。 [V_i, V_j]=Vnew_jishuan(jiedian_i,jiedian_j ,m,n) 其中，输出V_I,V_j分别为子波源点所在的结点相连的所有结点的行向量与列向量；输入中jiedian_I，jiedian_j分别为子波源点所在的结点的网络行号与列号，m，n为离散网络行