二维复杂结构三角网射线追踪全局方法.pdf

本发明涉及二维复杂结构的走时层析成像反演和以Maslov射线理论为基础的波场计算的射线追踪方法。对于二维复杂外部几何边界及区域内部复杂空间结构，用慢度分块均匀的三角形模型将介质参数化，在三角形单元的边界上设置节点。以波行面代替波前面，波行面的扩展、波的传递通过三角单元与其相邻三角单元的“振动”传递进行，遇回转波，会自行向后传递，无需另行“收缩”，各节点的次级源位置和最小走时计算在波行面扩展过程中实现；利用各节点次级源走时与方向信息，通过最小走时搜索，拾取从接收点到源点的射线路径。本计算方法适用于任意多边形网。

权利要求书

二维复杂结构三角网射线追踪全局方法
技术领域
本发明涉及一种二维复杂结构三角网射线追踪全局方法，具体是一种二维复杂结构的走时层析成像反演和以Maslov射线理论为基础的波场计算的射线追踪方法，属于勘探地球物理技术领域。
背景技术
射线理论和射线方法是认识波场传播规律的重要途径，是研究地下复杂构造和不均匀介质中的波场传播问题的重要手段，也是走时层析成像反演的基础。
无论是正演模拟还是层析反演，模型参数化是基本的问题。由于工程地质探测现场条件的复杂性，所探测研究对象的外部几何形状大部分情况下是不规则的形状，如“工”字型结构；内部地质异常体在空间结构上更加复杂，有的呈不规则柱状体(如溶洞、陷落柱等)，有的呈条带状(如断层、裂缝等)、有的呈层状(如软弱夹层等)。这就模型参数化提出更高的要求，以满足复杂结构体的高分辨率探测。目前模型参数化一般都采用矩形网或三角网，其中矩形网剖分是最普遍的形式，而三角网剖分也主要是在矩形网基础上简单的三角化。这些简单的矩形网、三角网剖分没有充分考虑区域复杂边界，也未将射线密度、介质特性同网格大小有机结合起来，难以满足复杂结构的射线追踪正演和高分辨率层析成像反演的要求。必须采用不规则网，如三角形网、不规则四边形网、或三角形网和不规则四边形网的混合网，对模型参数化。研究表明，相对于矩形网络参数化，三角网格参数化具有如下优点：①模型剖分的灵活性强；②对速度间断面的描述精度高，速度模型和界面模型具有一致性；③正演模拟存储量小、计算时间少；④正演模拟难度降低、误差减少、效率提高；⑤用于层析反演时存储量小、计算时间少、方程性态好、求解容易。
射线追踪的理论基础是，在高频近似条件下，弹性波场的主能量沿射线轨迹传播。传统的射线追踪方法，通常意义上包括初值问题的试射法(Shootingmethod)和边值问题的弯曲法(Bending method)。Vidale(1988)在提出程函方程的有限差分法时，指出试射法和弯曲法的主要问题在于：①难于处理介质中较强的速度变化；②难于求出多值走时中的全局最小走时；③计算效率较低；④阴影区内射线覆盖密度不足。
对于三角网模型的射线追踪问题，有关文献，作了探讨，但其给出的射线追踪方法基本归类于通常意义上的初值问题的试射法和边值问题的试射法和弯曲法，不可避免地存在其固有缺点。
对于矩形网模型的射线追踪问题，有关文献对常见的射线追踪方法做了比较详细的总结。为了克服传统试射法和弯曲法的这缺点，近年来，许多地球物理研究者在这方面进行了大量的工作，提出了一些精度较高、效率较高而且实用的计算初至旅行时的波前追踪方法。研究进展主要体现在：①在传统的试射法及弯曲法的基础上的改进，如各类波前重建方法(Vinje，1992；Sun，1992；Lambare et al.，1996)，除多值走时外，还较好地解决了计算效率及阴影区覆盖不足的问题；②对最小走时算法的改进，使之可适应多值走时计算，如慢度匹配法(Symes，1998)，可认为是最短路径方法的推广；③传统方法与最小走时算法的结合，如HWT方法(Sava，Fomel，1998)，则是通过波前传播计算射线路径。这些方法中，同时考虑空间所有离散点上的走时和射线路径的全局计算方法，由于其较高的精度和计算效率引起了广泛的关注。然而，这些主要针对规则网射线追踪方法，应用到三角网中会使问题复杂化，甚至难以实现。
三角网的优势使得模型三角网参数化是未来地球物理模型参数化的重要趋势，但其波场射线追踪问题却缺乏深入的研究，因此，三角网模型及其波场射线追踪问题是现阶段地球物理探测技术提高正演与反演效率、效果必须要解决的重要课题。
发明内容
针对上述存在问题，本发明的目的在于提供一种可对二维复杂结构三角网进行射线追踪的全局算法。
其技术解决方案如下：
二维复杂结构三角网射线追踪全局方法是以波行面代替波前面，波行面的扩展、波的传递通过一个三角单元与其相邻三角单元的“振动”传递进行，遇回转波，会自行向后传递，无需另行“收缩”，各节点的次级源确定和最小走时计算在波行面扩展过程中实现，其计算方法按以下步骤进行：
a)将二维复杂结构区域，用慢度分块均匀的三角形模型将介质参数化，在三角形单元的边界上设置节点，
b)设源节点的走时为零，除源节点外的其他节点的走时为无穷大，将源节点所在的三角形纳入波行面三角单元数组或优先队列，并使源节点所在的三角单元处于激活状态，作出存在于波行面数组的标志，
c)在波行面三角单元数组中查找处于激活状态的三角单元中的最小走时节点，确定这一节点所在三角单元，将此三角单元作为当前三角单元，如果这一节点是多个三角单元的公共节点，选其中一个三角单元即可，其它三角单元会自行遍历，将这个最小走时的当前三角单元及其相邻三角单元纳入波行面三角单元数组，并做出存在于波行面数组的标志，以便于不多次存入波行面三角单元数组中，
d)从最小走时三角单元的最小走时节点开始，计算、比较当前最小走时三角单元每个节点的旅行时间，并确定各节点的次级源，
e)再计算、比较最小走时三角单元的相邻三角单元节点的旅行时间及确定各节点的次级源，当相邻三角单元节点走时有变化时，将这一相邻三角单元作出激活标志，
f)再一次计算、比较当前最小走时三角单元节点旅行时间，当当前最小走时三角单元各节点走时无变化时，对当前最小走时三角单元作出“休眠”标志，回到第c)步，当当前最小走时三角单元各节点走时有变化时，对当前最小走时三角单元作出“激活”标志，回到第c)步，当所有波行面的节点走时都不再变化时且波行面中所有的三角单元均处于“休眠”状态，波行面扩展完毕。
g)从接收点开始，根据各节点走时及次级源信息，进行向源检索，拾取从接收点到源点的射线路径。
所述的波形面是指二维平面内，由波源点及某一时刻波到达的每个节点构成的平面。
所述的源节点是指激发源点。
所述的相邻三角单元是指与一个三角单元顶点或边相关的所有三角单元。
由于采用了以上技术方案，本发明针对二维复杂结构的射线追踪问题，提出了以波行面代替波前面的三角网射线追踪全局方法，数值模拟表明该算法具有很高的精度和可靠性，克服了传统的基于规则网的射线追踪方法适应性差的弱点，该方法使规则网成为此方法中的一个特例，适用于任意多边形网，如四边形网、四边形与三角形混合网。对于三维问题，稍加改造也同样适用。
与一些将波前刻画成一条曲线不同，本方法的特点是：将所遍历的三角单元作为一个波阵平面，一方面波行面的扩展、波的传递通过最小走时三角单元与其相邻三角单元“振动”传递进行，另一方面当遇回转波时，自行向后传递，无需另行“收缩”，亦即通过三角单元的相邻三角单元将先前已经作出“休眠”标志的三角单元“激活”，再次成为最小走时三角单元。这样就有效的解决了波的扩展和回传的问题。
二维复杂区域模型三角化及二维复杂结构三角网射线追踪全局方法，将能促进以射线理论为基础的波场正演研究工作，并使其正演效率极大提高。
附图说明
附图1速度模型参数化三角单元
附图2三角顶点的邻域三角顶点
附图3节点的邻域点，其中，(a)为节点是三角顶点情形，(b)为节点是三角边的插入点情形，(c)为节点在三角边的插入点之间情形，(d)为节点在三角单元内部情形
附图4一个三角单元的相邻三角单元
附图5一个三角单元的波前邻域点走时计算
附图6波前不同近似的比较，其中，(a)为离散波前点近似示意图，(b)为连续波前点近似示意图
附图7波行面从(a)到(h)的扩展过程
附图8各节点次级源计算实例
附图9节点次级源检索，其中，(a)为节点E在三角单元内情形，(b)为节点E在三角顶点或插入点上情形，(c)～(d)为节点E在边插入点之间
附图10二水平层介质模型计算结果，其中，(a)为三角单元划分结果和初至波射线路径图，(b)为初至波波阵面图，(c)为右边界走时对比，(d)为上边界走时对比
附图11水平层状模型初至波射线路径及波阵面图，其中，(a)为剖分网格，(b)为初至波射线路径及波阵面
附图12含空洞模型初至波射线路径及波阵面，其中，(a)为剖分网格，(b)为初至波射线路径及波阵面
附图13含断层模型初至波射线路径及波阵面，其中，(a)为剖分网格，(b)为初至波射线路径及波阵面
具体实施方式
下面结合附图及具体实例对本发明的计算方法作进一步详细说明。
见附图。
为便于更好理解本发明的技术方案，首先对二维复杂结构三角网射线追踪全局方法涉及到的相关概念进行描述。
如附图1，为二维区域三角剖分形成的三角单元。规定：三角单元的三个顶点按逆时针方向排列，如A→B→C，三边均为有向线段，如AB、BC、CA；每个单元内速度度均匀。在每个三角单元的3个边上设置相同多的节点，节点距离可以均匀分配，也可以非均匀分配，这里采用均匀分配，且三角单元公共边设置的节点相同。
二维区域内的所有点称为节点，包括源点、接收点，三角单元顶点及三角边的插入点等，对应地，称之为源节点、接收节点、三角单元顶点节点及三角边插入节点。
与三角单元的一个顶点相邻的所有三角单元顶点节点称为该三角单元顶点节点的邻域三角顶点，其方向定义为从该点指向邻域三角顶点，并约定邻域三角顶点绕该点逆时针排放。如附图2，P₀的邻域三角顶点分别为：P₁、P₂、P₃、P₄、P₅、P₆。
称一个三角形单元内部或边上的任意两个相邻不同节点为相邻点，并约定其指向性。一个节点的所有相邻点称为这个节点的邻域点，其方向定义为从该节点指向邻域点，并约定邻域点绕该节点逆时针排放，为减少不必要的邻域点，同一边上，只保留与该节点紧邻的节点。
对区域三角网而言，节点的邻域点在不考虑边界情况下主要有四种情形，如附图3，节点ss(空心圆)的邻域点rr(所有实心圆)，即：节点在三角单元顶点、节点在三角边的插入点、节点在三角边的插入点之间、节点在三角单元内部。对一个三角单元而言，对于附图3中(a)、(b)、(c)三种情形，一个节点在该三角单元中的邻域点只需考虑在该三角单元中的邻域点。相应的，一个节点的邻域点可能包括：三角单元顶点、三角边的插入点、三角边的插入点之间的点、三角单元内部点这四类节点。
在计算从一个已知节点走时到其邻域点的旅行时过程中，该已知走时节点称为其邻域点的波前点。如附图3中的ss，其邻域点称为该节点的波前邻域点，如附图3中的rr。一个已知走时节点的次级源即为波由源传播到该节点所经历路径上的上一个邻域节点。因此，在计算波前点到其波前邻域点的旅行时过程中，各波前点可以相继看做是其波前邻域点的次级源，但只有当遍历所有波前点时，各节点的次级源才是使该节点走时最小的上一个邻域节点。
与一个三角单元顶点或边相关的所有三角单元，称为该三角单元的相邻三角单元。如附图4，三角单元A₀的相邻三角单元包括A₁～A₁₂12个三角单元。在下面将要定义的波行面的扩展与回转波的回传就是通过相邻三角单元来扩展和传递的。
在二维平面内，某一时刻由波源点及波到达的每个节点构成的平面，称为波行面。在三角网射线追踪全局算法中，波行面的表现形式为，波遍历的三角单元组成的平面。若区域剖分为其它形式的多边形网，如四边形网、四边形与三角形混合网，则波行面的表现形式为，波遍历的四边形单元、或四变形单元与三角形单元混合组成的平面。值得说明的是，计算各节点的走时过程中，在完成整个区域的节点走时计算之前，波行面各节点的走时不一定必须是最小走时。
由上述概念，就可得到三角网射线追踪算法中运用较频繁的有向边、节点、三角单元等数据结构。
有向边类
class Edge              //有向边类
{
public：
    Vertex*dest；           //第二个点位置
    doublecost；           //边的代价，可以是距离，也可以是时间；
    Triangle*lefttriangle；//边的左三角形
    Triangle*righttriangle；//边的右三角形
}；
     节点类
class Vertex
{
public：
   double x；                     //节点x坐标
   double y；                     //节点y坐标
   dcllist<Edge>adj_tri_vex；//该点相邻三角顶点
   dcllist<Edge>adj_vex；      //该点邻域点，
                                 //dcllist为双向循环链表模板
   Vertex *prev_vex；            //该点的次级源
   bool scratch；                //搜索标志
   double density；              //该点三角单元密度控制大小
   double travel_time；          //由源点到该节点的旅行时
}；
   三角单元类
class Triangle//三角形类
{
public：
    dcllist<Vertex*>tripoly；//由三角单元顶点及边插入点构成的三角环
    Vertex*a，*b，*c；          //三角单元三个顶点
    double v；                      //三角单元波速
    bool scratch；              //是否被遍历标志(休眠、激活标志)
    bool insert；                   //是否进入波阵平面标志
    double min_time；               //最小走时
    DCLLIST::iteratormin_time_itr；//最小走时位置，
                                   //DCLLIST::iterator，三角环迭代位置
   dcllist<Triangle*>adj_tri；//相邻三角单元
}；
二维复杂结构三角网射线追踪全局方法是利用波行面的扩展完成节点最小走时计算和次级源的确定，而波行面的扩展是通过一个单元与其相邻单元的“振动”传递完成。
对于一个复杂区域，采用三角剖分方法将区域介质参数化，视每个三角单元慢度均匀，并设置节点及三角网的一些拓扑关系。三角网射线追踪全局计算方法的基本步骤如下：
第一步，利用波行面的扩展完成节点最小走时计算和次级源的确定。以波行面代替波前面，波行面的扩展、波的传递通过三角单元与其相邻三角单元的“振动”传递进行，遇回转波，会自行向后传递，无需另行“收缩”。
第二步，利用计算出的各节点的旅行时和次级源信息拾取射线路径。即根据每个节点的最小走时、次级源，从接收点出发，向源拾取射线路径。
二维复杂结构三角网射线追踪全局方法的实现过程分以下几方面。
1 设置各节点的拓扑关系
将二维复杂结构区域，用慢度分块均匀的三角形模型将介质参数化，在三角形单元的边界上设置节点，同时设置各节点之间的拓扑关系，主要包括以下几类节点的拓扑关系：
(1)按逆时针方向，设置各节点的邻域三角顶点，这里的节点包括源点、接收点、三角顶点及三角边的插入点，同时设置由节点至邻域三角顶点形成的有向边的左、右三角形，主要目的是便于查询三角形的相邻三角形。
(2)按逆时针方向，设置各节点的邻域点，这里的节点包括源点、接收点、三角顶点及三角边的插入点，同时设置由节点至其邻域点形成的有向边的左、右三角形，主要目的是便于计算由波前点至波前邻域点的旅行时。将接收点和源点等同于三角顶点及三角边的插入点进行设置的一个优点是使得射线追踪不必对源点和接收点另行计算，路径拾取只需给出接收点坐标或接收点编号即可。
(3)按逆时针方向，设置各三角单元的相邻三角单元，根据三角顶点的邻域三角顶点设置，主要目的是便于三角单元之间波的传递和波行面的扩展。
2 波前邻域点最小走时及次级源位置计算
波前邻域点最小走时计算及次级源的确定只需考虑在一个三角单元内进行，本步骤藴含于波行面的扩展中，波行面扩展细节见下一部分的描述。
附图5是一个三角单元，设ss点是波前节点，其波前邻域点用实心圆表示，其中rr节点是ss的某一邻域点，它们的坐标分别为(x_ss，y_ss)和(x_rr，y_rr)；(x_rs，y_rs)为rr节点次级源的坐标；t_ss为ss节点的走时；t_rr为rr节点的走时；v为三角单元的速度；t_sr为波从ss节点以本三角单元速度v沿直线传播到rr节点所需的时间：
t sr = [ ( x ss - x rr ) 2 + ( y ss - y rr ) 2 ] 1 2 / v - - - ( 1 ) ]]>
于是可以求出rr节点的走时及相应的次级源节点坐标，即当(t_ss+t_sr)<t_rr时，有
t rr = t ss + t sr x rs = x ss y rs = y ss - - - ( 2 ) ]]>
在一个三角单元内计算节点走时时，从这个三角单元走时最小的节点ss开始，然后让ss遍历本三角单元的所有节点，本三角单元的节点走时计算才算完毕，但各节点走时不一定就是最小走时，只有当ss经历所有波前点时，t_rr才变成最小走时，(x_rs，y_rs)才是真正次级源坐标。
这里的次级源是用节点近似，波前邻域点的走时是用波前点走时加上波前点到波前邻域点局部走时的和表示的，相当于把波前面看成是一系列离散点次级源，如附图6(a)。实际上到达每一节点的波并不一定正好是从网格上的波前点发出，而可能是从网格点之间的点发出，如附图6(b)。这样局部走时和次级源位置就要进行修正。
如附图6(b)，要计算节点E的最小走时和次级源的位置，假定E的次级源在A、B、C三点之间，如D点，设某一点为局部原点(如A点)，在A，B，C三点之间走时函数表示为t(x)，x为AC方向上局部一维坐标。设前方点E在AC方向上的局部坐标为x_E，垂直AC方向的距离为d，则经AC线上任一点x到达邻域点E的走时为
t E = t ( x ) + [ ( x - x E ) 2 + d 2 ] 1 2 / v - - - ( 3 ) ]]>
v为三角单元内的速度，令x在AC之间变化时，t_E取极小，可求出波前点即次级源点D：
{ x D = x : t E = min x&Element; ( A , C ) [ t ( x ) + ( ( x - x E ) 2 + d 2 ) 1 2 ] / v } - - - ( 4 ) ]]>
若A、B、C3点的局部坐标为x_A、x_B、x_C，3点的走时为t_A，t_B，t_C，则用双曲线近似计算AC上任一点走时
t ( x ) = { t A 2 [ ( x - x B ) ( x - x C ) ] / [ ( x A - x B ) ( x A - x C ) ] + ]]>
    t B 2 [ ( x - x A ) ( x - x C ) ] / [ ( x B - x A ) ( x B - x C ) ] + - - - ( 5 ) ]]>
    t C 2 [ ( x - x A ) ( x - x B ) ] / [ ( x C - x A ) ( x C - x B ) ] } 1 / 2 ]]>
容易证明，当介质均匀时，t(x)是某一点源引起的走时时，上式是严格精确的。当介质不均匀时，双曲线近似也是3点走时插值方法中精度最高的。(4)式是一个极小值求解问题，可采用梯度法、黄金搜索等方法求得x_D。按边让A、B、C历遍除E点所在边的其它两边的节点就完成了E点在这个三角形中的最小走时和次级源位置的计算。
3 波行面扩展
在基于矩形网的射线追踪的全局方法中，常利用一矩形模拟波阵面的扩张与收缩来计算节点最小走时和次级源位置，但在三角网中，要拟定某一时刻的波阵线将会非常复杂，而且采用波阵面，也不利于回转波的传递。
本发明的做法是，利用波行平面，将波行平面的扩展视为一个递归过程，其算法描述如下：
(1)设源点的走时为零，除源点外的其他节点的走时为无穷大。将源点所在的三角形纳入波行面三角单元存储容器，并使这个三角单元处于激活状态，作出存在于波行面的标志。
(2)在波行面三角单元存储容器中查找处于激活状态的三角单元中的最小走时节点，确定这一节点所在三角单元，将此三角单元作为当前三角单元；如果这一节点是多个三角单元的公共点，选其中一个三角单元即可，其它三角单元会自行遍历。将这个最小走时的当前三角单元及其相邻三角单元纳入波行面三角单元存储容器，并做出存在于波行面的标志，便于不多次存入波行面三角单元存储容器。
(3)利用前述“波前邻域点最小走时及次级源位置计算”，计算、比较当前最小走时三角单元节点旅行时(从本三角单元最小走时节点开始，下同)；再计算、比较最小走时三角单元的相邻三角单元节点的旅行时，当相邻三角单元节点走时有变化时，将这一相邻三角单元作出激活标志；再一次计算、比较当前最小走时三角单元节点旅行时，当当前最小走时三角单元各节点走时无变化时，对当前最小走时三角单元作出“休眠”标志，回到第(2)步，当当前最小走时三角单元各节点走时有变化时，对当前最小走时三角单元作出“激活”标志，回到第(2)步；以次类推，遍历所有三角单元。当所有波行面的节点走时都不再变化时且所有波行面中的三角单元均处于“休眠”状态，波行面扩展完毕。
与一些将波前刻画成一条曲线不同，本计算方法的特点是：将所遍历的三角单元作为一个波阵平面，以波行面代替波前面；一方面通过最小走时三角单元的与其相邻三角单元间的“振动”传递扩展波行面，另一方面当遇回转波时，也通过三角单元的相邻三角单元之间的传递将先前已经“休眠”的三角单元“激活”，而再次成为最小走时三角单元，无需另行“收缩”。这样就有效的解决了波的扩展和回传的问题。
附图7展示了波行面从(a)到(h)的扩展过程，源点在左上角，粗线三角单元表示当前三角单元。
附图8为一计算实例，源点在左上角，显示出各节点的次级源，如A的次级源为B。
4 波传播路径的向源拾取
任意点到源点的射线路径通过节点最小走时和次级源信息向源追索完成。本发明将源点和接收点按照节点统一处理，因而，在向源路径拾取过程中不必对接收点走时另行计算，直接追索即可。附图9表示节点E的次级源F的检索，根据E的位置，有以下几类节点的次级源检索：
(1)当节点E位于三角单元内部时，如附图9(a)，通过双曲线近似和最小走时搜索从三角单元三边上找到次级源，算出该节点走时(这一类节点可在波行面扩展中统一处理，这里只需检索即可)；
(2)当节点E是三角顶点或三角边的插入点时(9(b))，可直接拾取次级源；
(3)当节点E位于三角边插入点A、C之间时，检索E的次级源F要依据A的次级源B与C的次级源D。这时主要会遇到如附图9(c)、(d)两种情况：(c)情形是A的次级源B与C的的次级源D在一个三角单元内；(d)情形就比较复杂，由于C点是三角顶点，C的次级源D可能存在于其它三角单元内。对附图9(c)、(d)两种情况，也通过双曲线近似和最小走时搜索从三角单元相应边上寻找次级源，并算出该节点走时。
5 算法精度检验
为了检验三角网射线追踪算法的精度和可靠性，对附图10(a)所示模型的三角网射线追踪正演模拟结果同根据Snell定理计算的理论结果进行了对比。源位于左边界上端，接收点分别位于模型上边界和右边界，接收点距2m。第1层介质速度为2000m/s，第2层介质速度为4000m/s，模型水平方向宽度为50m，垂直方向深度分别是第1层10m，第2层40m。在对模型进行网络划分时，采用密度为2m的函数进行三角剖分。附图10(a)为三角网剖分结果和初至波射线路径图(粗线)；附图10(b)为初至波波阵面图，等时线间隔1ms。附图10(c)、(d)分别为右边界和上边界走时对比结果，图中实线为用解析分析公式计算的理论初至波走时，圆点为三角网追踪算法计算的初至波走时。右边界最大走时相对误差为0.0009％；上边界最大走时相对误差仅为0.00028％。由此可以看出，采用三角网射线追踪算法具有相当高的计算精度。
具体实施例
实施例一  水平层状均匀介质模型
附图11为水平层状模型的三角剖分及射线追踪结果。源点位于左上角，接收点分别位于模型上边界、右边界及下边界。模型示于附图11(a)，第1层介质速度为2000m/s，第2层介质速度为4000m/s，第3层介质速度为3000m/s，第4层介质速度为5000m/s，模型水平方向宽度为50m，垂直方向4个层厚度均为10m，对模型采用密度为2m的函数进行三角剖分。附图11(a)为三角剖分结果，附图11(b)为采用本算法得到的不同时刻的波阵面和射线路径。
实施例二  含空洞模型
附图12为含空洞模型的三角剖分及射线追踪结果。源点位于左上角，接收点分别位于模型上边界、右边界及下边界。模型示于附图12(a)，背景速度为4000m/s，设置：速度为2000m/s的低速矩形区、速度为5000m/s的高速矩形区及隧道型状的空洞，网格剖分采用不均匀密度控制函数。附图12(a)为三角剖分结果，附图12(b)为采用本算法得到的不同时刻的波阵面和射线路径。
实施例三  含断层模型
附图13为含断层模型的三角剖分及射线追踪结果。源点位于左下角，接收点分别位于模型外边界。模型示于附图13(a)，在区域内设置了断层子域Ω₂及低速子域Ω₅，为模拟Ω₃到Ω₅速度较大的梯度变化，增加网格加密子域Ω₄，网格剖分采用了不均匀密度控制函数，各子域的模拟速度Ω₁∶v_p＝3000m/s，Ω₂：v_p＝2000m/s，Ω₃∶v_p＝4000m/s，Ω₄∶v_p＝4000m/s，Ω₅∶v_p＝3000m/s。附图13(a)为三角剖分结果，附图13(b)为采用本算法得到的不同时刻的波阵面和射线路径。