基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法

1.本发明属于地震勘探技术领域，涉及一种基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法。

背景技术：

2.地震波的属性可以大致分为运动学和动力学两方面。在繁多的地震数据处理和反演方法中，其对于地震波数值模拟的要求也不尽相同。诸如逆时偏移、全波形反演等方法要求对全波场进行数值模拟，而类似常规偏移成像、微震定位、层析成像等方法则只需要模拟波场的运动学规律即可，单纯的运动学数值模拟的计算量要远小于全波场。
3.射线追踪是地震波运动学数值模拟的重要方法。目前，射线追踪主要包含以下几类：试射法，弯曲法，伪弯曲法，逐段迭代法，反向追踪法，最短路径法，插值法等。其中，试射法和逐段迭代法的理论基础是斯奈尔定律，反向追踪法、弯曲法和伪弯曲法的理论依据是费马原理，最短路径法和插值法则是在费马原理上又结合了惠更思原理的思想。
4.逐段迭代法先固定射线两端点，再对射线路径进行迭代更新。逐段迭代法不会跨越多个界面进行整体更新，故其稳定性相比伪弯曲法有所提高，但在速度变化剧烈的介质中效果有限。这一类方法能够被应用于网格模型中，但无法适应网格模型中存在速度剧烈变化的情况，即在复杂介质中存在一定的局限性。
5.反向追踪法是以走时场为基础，从检波点向震源进行追踪的，其对于复杂介质和网格模型具有较强的适用性，但反向追踪法是沿着走时场的负梯度方向进行追踪的，精度与梯度方向的计算有很大关系。这一类方法不存在对整条射线进行迭代更新的步骤，因此计算效率通常比逐段迭代法要高，但相应的计算结果精度就无法获得有效提升；另一方面，由于追踪是从检波点开始的，故误差会在震源附近累计。
6.随着地震勘探向高精度发展，射线追踪法也越来越需要适应以高密度网格作为介质模型的情况。因此，需要一种射线追踪方法既能够适应复杂的速度场变化，也能够在震源附近获得较高的精度。

技术实现要素：

7.为解决上述缺陷，本发明提供了一种基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，解决了常规逐段迭代法存在的问题：收敛速度过慢、可能陷入局部极小、无法追踪回折波，本发明既能够适应复杂的速度场变化，也能够在震源附近获得较高的精度。
8.基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其包括以下步骤：
9.s1、将常规反向追踪法的追踪结果作为逐段迭代法的初始值；
10.s2、将原始网格抽稀到预设的最稀疏级别，同时只保留剩余原始网格边界上的路径点；
11.s3、对抽稀网格应用常规逐段迭代，并改变抽稀网格边界的方向，通过改变空间离散化的方式，使得回折路径变为满足斯奈尔定律的路径。
12.s4、对抽稀网格逐渐加密，每次加密重复步骤s3，直到当前的抽稀级别降为0，再进行最后一遍常规逐段迭代。于本发明的一种实施方式中，所述步骤s2具体包括以下流程：
13.2-1在原始网格中，给定初始射线和追踪点；
14.2-2设定抽稀级别为n并将网格抽稀到最稀疏级别，只保留剩余网格边界上的路径点。
15.于本发明的一种实施方式中，在步骤2-2中，抽稀方式为每多一级抽稀就对当前网格模型中隔一条线保留一条，n级抽稀相当于每2n条线保留一条。
16.于本发明的一种实施方式中，所述步骤s3具体包括以下流程：
17.3-1在抽稀网格中进行逐段迭代，并跳过不满足斯奈尔定律的追踪点；
18.3-2将抽稀网格和其中的追踪点绕原点顺时针旋转；
19.3-3利用与抽稀网格相同新网格间距的水平-垂直网格覆盖旋转后的抽稀网格；
20.3-4利用插值，计算旋转后的射线与新网格边界的交点，即新网格中的追踪点；
21.3-5将新网格和其中的追踪点平移到非负区域；
22.3-6对新网格中的射线进行逐段迭代更新，同时跳过不满足斯奈尔定律的点；
23.3-7将新网格及对应的追踪点变换回原始网格及对应的追踪点。
24.于本发明的一种实施方式中，在步骤3-2中，旋转的角度范围为0-90
°
。
25.于本发明的一种实施方式中，在步骤3-6中，由于网格边界和射线的相对方向发生改变，新网格中不满足斯奈尔定律的点与步骤3-1不同。
26.于本发明的一种实施方式中，所述步骤s4对抽稀网格逐渐加密成新网格并利于插值求得新的追踪点，同时添加对应新网格边界上的路径点，并在每一个稀疏级别上都进行步骤s3。
27.于本发明的一种实施方式中，所述步骤s4具体包括以下流程：
28.4-1对当前抽稀网格进行一次加密，加入的网格线位于相邻两条网格线的中间，然后利用插值求得新的追踪点；
29.4-2在当前抽稀级别的网格上进行步骤s3；
30.4-3重复步骤4-1和4-2，直到当前的抽稀级别降为0，再进行最后一遍常规逐段迭代。
31.综上所述，本发明提供一种基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，本发明的有益效果是：
32.本发明针对常规逐段迭代法收敛过慢，提出了一种多尺度的逐段迭代法，即将矩形网格进行抽稀，以减少模型中界面的数量，抽稀到预设的最稀疏级别再逐渐加密网格到原始状态，有效提升收敛速度；其次针对常规逐段迭代法无法追踪存在回折射线路径的问题，提出了一种多方向的逐段迭代法，即通过改变网格边界的方向，从而改变空间离散化的方式，使得回折路径变满足斯奈尔定律的路径，实现了回折波追踪；然后针对常规逐段迭代法局部收敛极小的问题，提出了一种基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，即将常规反向追踪法的结果作为改进后的逐段迭代法的初始值，有效避免了迭代过程陷入局部极小。本发明既能够适应复杂的速度场变化，也能够在震源附近获得较高
的精度。
附图说明
33.图1为反向追踪单点追踪示意图。
34.图2为逐段迭代法修正示意图。
35.图3为反向追踪法数值模拟结果图；其中，图3(a)为正演模型；图3(b)为反向射线追踪法数值模拟结果图；图3(c)为对图3(b)中震源区域追踪结果局部放大。
36.图4为反向追踪法算例数值模拟的误差与统计；其中，图4(a)为正演模型；图4(b)为真实射线(实线)和反向追踪结果(虚线)示意图；图4(c)为单射线整体误差统计；图4(d)为单射线分布误差统计。
37.图5为逐段迭代法在矩形网格应用存在问题的示意图；其中，图5(a)为逐段迭代法跨网格修正流程示意图；图5(b)为逐段迭代法修正过程中局部极小示意图；图5(c)为回折波等类型的射线路径示意图；图5(d)为矩形网格模型中回折波射线路径示意图。
38.图6为针对问题(1)数值模拟结果图；图6(a)为正演模型；图6(b)为逐段迭代法数值模拟结果图。
39.图7为针对问题(2)数值模拟结果图；图7(a)为正演模型；图7(b)为逐段迭代法数值模拟结果图。
40.图8为针对问题(3)数值模拟结果图；其中，图8(a)为正演模型；图8(b)为逐段迭代法数值模拟结果图(只展示了射线路径所在的区域)。
41.图9为多尺度逐段迭代法示意图；其中，图9(a)为原始网格、初始射线路径及对应追踪点；图9(b)为对图9(a)进行2级抽稀；图9(c)为对图9(b)进行常规逐段迭代；图9(d)为对图9(c)进行加密；图9(e)为对图9(d)进行常规逐段迭代。
42.图10为多尺度逐段迭代法数值模拟；其中，图10(a)为正演模型；图10(b)为多尺度逐段迭代法数值模拟结果图。
43.图11为多方向逐段迭代法示意图；其中，图11(a)为回折射线路径示意图；图11(b)为回折射线顺时针45度旋转后路径示意图。
44.图12为多方向逐段迭代网格变换示意图；其中，图12(a)为原始网格；图12(b)为顺时针旋转45度后的旋转网格；图12(c)为与原始网格相同网格间距的水平-垂直网格覆盖后的新网格；图12(d)为将覆盖后的新网格平移到非负区域。
45.图13为多方向、多尺度逐段迭代数值模拟结果图；其中，图13(a)为正演模型；图13(b)为多方向、多尺度逐段迭代法数值模拟结果图。
46.图14为基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪法的工作流程图。
47.图15为基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法数数值模拟结果图；其中，图15(a)为正演模型；图15(b)为基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法数值模拟结果图；图15(c)为射线水平段局部放大结果图。
48.图16为基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法算例数值模拟的误差与统计；其中，图16(a)为联合射线追踪方法数值模拟结果图；图16(b)为对图16(a)中震源区域追踪结果局部放大；图16(c)为联合射线追踪方法单射线整体误差统计；图
16(d)为联合射线追踪方法单射线分布误差统计。
具体实施方式
49.为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
50.1常规射线追踪方法的基本原理
51.1.1反向追踪法的基本原理
52.反向追踪法的基本理论依据是费马原理，即：射线方向是波前面(等时面)的法线方向，射线路径总是垂直于波前面(等时面)。因此，在走时场中，某一点的梯度方向即是经过该点的射线的前进方向。
53.图1为反向追踪单点追踪示意图。其中，两个正方形表示两个网格，p点为当前追踪点，pa和pb表示p点在走时场中的垂直和水平方向的偏导，pc表示p点梯度，pd为负梯度方向，d为下一个追踪点。
54.基于此，要想求得两点之间的射线路径，可以按照以下步骤进行实现：
55.(1)将检波点作为路径上的第一个追踪点(如图1中的p点所示)；
56.(2)计算走时场中该点的垂直和水平方向的偏导(如图1中的pa和pb所示)，以及梯度(如图1中的pc所示)；
57.(3)取负梯度方向(如图1中的pd所示)，并沿该方向寻与网格的第一个交点(与边的交点或恰好是网格点)(如图1中的d点所示)，作为路径上的下一个追踪点；
58.(4)重复(2)和(3)，直到追踪到震源点为止。
59.1.2逐段迭代法的基本原理
60.逐段迭代法的基本理论依据是斯奈尔定律。考虑到同一射线路径上的任意连续三点都满足斯奈尔定律，可以先将震源点和检波点用直线相连，连线与路径上穿过的界面有若干交点，作为初始的追踪结果；之后，按照斯奈尔定律，对射线上的所有连续三点的中间点进行修正，即让所有连续三点的位置关系都满足斯奈尔定律；重复此步骤，直到所有修正量之和小于某个预设的阈值，即认为完成射线追踪。
61.如图2所示，z＝z2表示某一界面，界面上下方的速度分别为v1和v2，p1与p3为界面两侧且紧邻界面的追踪点，将两点用直线相连，与界面相交于点p2，p1、p2、p3即是初始追踪结果。设p
′2点为真实的透射点，则由p2到p
′2的修正量δx可以表示为：
[0062][0063]
式中，a＝x
2-x1，b＝x
3-x2，
[0064]
2.常规射线追踪方法的误差分析
[0065]
2.1反向追踪法的误差分析
[0066]
反向追踪法的实现包含两个要点，一是对走时场的计算，二是对走时场中某点梯度的计算。对于第一个要点，走时场的计算本身也是一个较为复杂的问题，本发明中所有算例均使用快速推进法进行走时场的数值模拟。对于第二个要点，由于速度模型与走时场都是经过空间离散化的，故梯度计算中的微分需要用差分来代替，而用差分来代替微分实际上隐含了“走时场在局部是线性分布的”这一假设。
[0067]
根据反向追踪法的原理和流程，分析该方法存在如下误差因素：一是除震源、检波器的位置外，射线追踪流程中只输入速度模型，但最终的射线追踪是在走时场中完成的，因此最终射线路径的误差不仅来源于方法本身，还来源于前一个步骤中走时场的计算结果，而走时场的近场误差要大于远场误差。二是“走时场在局部是线性分布的”这一假设并不总是成立，其一是实际地质情况会导致追踪路径的偏差，并且这种误差会在追踪的过程中不断地累积，致使越接近震源位置误差累积的越大；其二是波前面在远场时曲率更小，更接近平面波，致使假设条件在远场比在近场要更容易成立。基于上述误差因素，该方法得到的射线追踪结果在震源附近的误差应当远大于远离震源的区域。
[0068]
利用数值模拟对反向追踪法进行测试。采用如图3(a)所示的正演模型：均匀模型，速度设置为5000m/s，并仅设置了一个震源点，位于s(2500m,3000m)；在地表50m-4950m范围内均匀布设了10个检波器(r1-r10)；使用正方形网格进行空间离散，网格大小为800*500，网格间距为10m。
[0069]
图3(b)展示了反向追踪(虚线)和解析解(实线)在整个模型尺度下的追踪结果对比，图中二者之间并不存在明显的肉眼可见的误差。对该结果在震源附近的区域进行放大，如图3(c)所示，可以观察到反向追踪结果和解析解之间存在明显的误差，且具有两个特点：一是射线出现锯齿状弯折，二是多条射线的路径发生合并。
[0070]
上述算例直观地展示了反向追踪法的整体效果和震源附近的误差特点，为了更进一步地量化误差的分布，在图3(a)所示正演模型的基础上将地表的检波器增加为500个，其余参数设置不变，如图4(a)所示。在均匀模型中由震源到一个检波点的真实射线(实线)和反向追踪结果(虚线)如图4(b)，在量化该追踪结果的整体误差时，求取了这两条线之间的平均距离，称其为单射线整体误差，当探究误差在射线路径上的分布时，则以这两条线之间的距离的分布作为该误差的量化，称其为单射线分布误差。基于上述两种误差量化方法，对500条射线进行误差分析，得到对应的单射线整体误差和单射线分布误差统计。如图4(c)(d)所示。
[0071]
反向追踪法的单射线整体误差如图4(c)所示，误差分布大致关于震源位置对称，误差变化较为平滑；从相对大小来看，射线与垂直方向的夹角在0度、
±
30度方向时，误差存在明显的局部极小，其中
±
30度左右的误差为全局极小；从绝对大小来看，误差最高接近15m，相当于射线整体偏离了真实位置约1.5个网格间距，所有检波点的误差的平均值在9.2m左右，已经接近1个网格间距。
[0072]
反向追踪法的单射线分布误差如图4(d)所示，由于射线总长度不同，故图中右侧部分存在空白区域。从纵轴方向来看，其表现出的规律与图4(c)基本一致；从横轴方向来看，其误差分布特征符合近场误差远大于远场误差的特点；从绝对大小来看，部分射线的局部区域的误差已经高达25m左右，即2.5个网格间距。
[0073]
2.2逐段迭代法的误差分析
[0074]
根据逐段迭代法的原理和流程，分析该方法存在如下问题：
[0075]
(1)由于每次修正只涉及到单一界面及其两侧的两个速度值，而不同网格有不同的速度值，因此这种修正是无法跨网格进行。如果要实现跨网格修正，则需要先把靠近网格点的追踪点归位到该网格点上，再对与网格点重合的追踪点进行修正。此过程如图5(a)所示，其中虚折线、实折线、点划折线代表三条射线路径，需要先将虚折线路径人为归位到实折线路径，再进行修正，得到点划折线路径。即便如此，一次修正最多也只能跨过一个网格，这可能会导致该方法的收敛速度过慢。
[0076]
(2)修正过程中所涉及的网格始终局限在射线路径及其周围一到两个网格间距以内的部分，这导致该方法并不能获得其他区域的速度信息，因而该方法的结果只代表一个局部极小走时而非全局极小走时，其违背了费马原理。例如在图5(b)中，假设远大于，则点到点的射线路径应如虚折线所示折射波，而在逐段迭代法中，按照流程，需先将两点用直线相连(点划直线)，这条直线上的任意相邻三点都已经满足斯奈尔定律，故修正量直接为零，迭代无法进行，也无法得到正确的射线路径。
[0077]
(3)初至波可能存在回折等情况，如图5(c)所示，即射线先向下(或上/左/右)传播，再向上(或下/右/左)传播。此时，在射线路径由向下到向上的拐点处，必然有类似图5(d)所示的情况出现，即对于连续的三个路径点，中间点的纵坐标小于两端点纵坐标的最小值。对于其他方向的回折，也会出现类似的情况。然而，在斯奈尔定律中，当界面水平或垂直时，中间点的横纵坐标一定在两端点的横纵坐标之间。也就是说，这种情况是不满足斯奈尔定律的，因此逐段迭代法是无法在矩形网格模型中追踪回折波或类似波场的射线的。
[0078]
针对上述逐段迭代法存在的三个问题进行数值模拟并将反向追踪结果作为近似的解析解进行对比分析。针对第一个问题，构建正演模型如图6(a)所示，该模型是将图3(a)所示模型分为了两层，浅层的底界面深度为2000m，速度为2000m/s，设置1个检波点r(4950m,0m)，其余参数不变。为了将两种方法的计算量调整到相同程度，以保证实验变量相同，设定逐段迭代法的迭代次数为120次。逐段迭代法数值模拟结果如图6(b)所示，逐段迭代法120次迭代的结果几乎与初始值没有差别，故认为追踪失败。
[0079]
针对第二个问题，构建正演模型如图7(a)所示，该模型在图6(a)所示模型基础上将震源、检波器位置分别设定在s(2000m,1500m)和r(6000m,1500m)，迭代次数保持120次不变。逐段迭代法数值模拟结果如图7(b)所示，震源与检波器之前初至波应为折射波，而逐段迭代法的追踪结果却保持初始值未变，故认为追踪失败。
[0080]
针对第三个问题，构建正演模型如图8(a)所示，该模型在图6(a)所示模型类似，将二层速度模型变为空间匀速变化模型，速度沿箭头方向线性增加，(0m,0m)处的速度值为2000m/s，(5000m,5000m)处的速度值为5000m/s，震源、检波器位置分别为s(2000m,500m)和r(6000m,500m)，迭代次数保持120次不变。逐段迭代法数值模拟结果如图8(b)所示，震源与检波器之前初至波应为回折波，而逐段迭代法的结果仍然保持初始值未变，故认为追踪失败。
[0081]
综合逐段迭代法存在的三个问题与其对应的数值模拟分析结果，可知，在矩形网格模型中，常规的逐段迭代法无法在合理的耗时内给出正确结果，甚至在某些情况下无法对初始值做出有效的迭代更新。
[0082]
3基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法
[0083]
在矩形网格模型中，常规逐段迭代法存在收敛过慢、局部极小、不能追踪回折波这三个问题，而造成这三个问题的主要原因则分别为：模型中界面数量过多(所有网格边界都被当成界面)；初始射线与真实射线相差较大；界面均为水平和垂直方向以及斯奈尔定律的理论限制，针对上述原因，本发明将对常规逐段迭代法进行改进。
[0084]
3.1多尺度逐段迭代法
[0085]
由于矩形网格中界面过多，造成逐段迭代法收敛过慢，为了解决这一问题，基于对矩形网格进行抽稀，以减少模型中界面的数量，进而提高逐段迭代法计算效率的思路，提出了一种多尺度的逐段迭代法，该方法首先将网格抽稀到预设的最稀疏级别，同时只保留剩余网格边界上的路径点，在此基础上进行常规逐段迭代；然后逐渐加密网格，同时添加对应网格边界上的路径点，并在每一个稀疏级别上都进行常规逐段迭代；最后重复上一步骤直至网格加密到原始状态。具体流程如下：
[0086]
(1)在原始网格中，给定初始射线和追踪点。以图9(a)为例，在一个网格点数为9*9的正方形网格模型中，要求取p1和p
16
之间的射线路径，首先将p1和p
16
两点用直线相连，并插值求出路径上的追踪点p2到p
15
。
[0087]
(2)设定抽稀级别为n并将网格抽稀到最稀疏级别，只保留剩余网格边界上的路径点。本方法所采用的抽稀方式是每多一级抽稀就对当前网格模型中隔一条线保留一条，n级抽稀相当于每2n条线保留一条。例如，对图9(a)进行2级抽稀，抽稀结果如图9(b)所示。其中，9*9的网格点阵变成了3*3，路径上p2到p
15
这14个追踪点只剩下p8和p9两个，将其重新编号为p2和p3，p4则为之前的端点p
16
。
[0088]
(3)在当前的网格上应用常规逐段迭代法。例如，对图9(b)进行常规逐段迭代，结果如图9(c)所示。
[0089]
(4)对当前网格进行一次加密，加入的网格线位于相邻两条网格线的中间，即相当于抽稀级别变成n-1级；然后利用插值求得新的追踪点。例如，对图9(c)进行一次加密，结果如图9(d)所示，网格点阵变成了5*5。加密之前的p2、p3、p4三点对应加密后的p4、p5和p8，而加密后的p2、p3、p6、p7四点则为新的追踪点。
[0090]
(5)在当前抽稀级别的网格上进行常规逐段迭代法。例如，对图9(d)进行常规逐段迭代，结果图9(e)所示。
[0091]
(6)重复步骤(4)和(5)，直到当前的抽稀级别降为0，即与原始的网格模型相同，再进行最后一遍常规逐段迭代。
[0092]
利用多尺度逐段迭代法对图6所示模型(同为图10(a))进行数值模拟，考虑到初始射线与真实射线相差较大，故给定一个较大的最大抽稀级别，此处设定为6，而每个级别的迭代次数为20次，总次数保持120次不变。追踪效果如图10(b)所示，多尺度逐段迭代法与常规逐段迭代法相比，追踪效果有明显提升，与反向追踪结果有大幅重叠，可认为其精度提升到了与反向追踪法接近的程度。
[0093]
3.2多方向逐段迭代法
[0094]
由于矩形网格模型的界面存在方式以及斯奈尔定律的理论限制，逐段迭代法无法追踪存在回折的射线路径，如图11(a)所示。对于此问题，本发明考虑了矩形网格模型中界面的特点，即：网格只是介质模型进行空间离散化的结果，网格边界并不代表实际介质中的
真实边界，其真实边界信息是隐含在大量网格之中的，而改变空间离散化的方式并不会破坏这些信息。根据上述特点，通过改变网格边界的方向，从而改变空间离散化的方式，来使得回折路径变为满足斯奈尔定律的路径，如图11(b)所示。为了便于观察，图11(b)是对射线进行了旋转，而非网格边界，这与旋转网格是等价的。经过旋转处理后，原本不满足斯奈尔定律的点就可以利用逐段迭代法进行更新。
[0095]
基于上述思路提出了多方向逐段迭代法，其中“方向”是指网格和射线的相对方向。该方法的具体流程如下：
[0096]
(1)在原始网格中(如图12(a))进行逐段迭代，并跳过不满足斯奈尔定律的追踪点。
[0097]
(2)将原始网格和其中的追踪点绕原点顺时针旋转0-90度范围内任意角度(文中以旋转45度为例)，如图12(b)所示。
[0098]
(3)利用与原始网格相同网格间距的水平-垂直网格覆盖旋转后的原始网格，如图12(c)所示。
[0099]
(4)利用插值，计算旋转后的射线与新网格边界的交点，即新网格中的追踪点。
[0100]
(5)将新网格和其中的追踪点平移到非负区域，如图12(d)所示。
[0101]
(6)对新网格中的射线进行逐段迭代更新，同时跳过不满足斯奈尔定律的点。其中，由于网格边界和射线的相对方向发生改变，新网格中不满足斯奈尔定律的点与第(1)步中的不同。
[0102]
(7)将新网格及对应的追踪点变换回原始网格及对应的追踪点。
[0103]
以上步骤即是多方向逐段迭代法的一次完整迭代，其中实际上包含了原始网格和旋转网格对应的两次迭代。
[0104]
利用多方向逐段迭代法对图8所示模型(同为图13(a))进行数值模拟。选择原始网格和45度网格作为多方向方法的参数，为了加快收敛速度，在本次测试中加入了多尺度方法，设定参数为6级抽稀以及每级迭代20次(原始网格和45度网格各20次)，其余参数不变。追踪结果如图13(b)所示，其中只展示了射线路径所在的区域。为了保证追踪结果确实来自多方向逐段迭代法，图13(b)还给出了只使用了多尺度方法的结果作为对比。多方向、多尺度逐段迭代法的结果与反向追踪的结果有很高的重合度，证明了该方法的有效性；而单纯采用多尺度逐段迭代法则完全没有效果，因为初始射线已经与经过的所有网格边界垂直了。3基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法
[0105]
多尺度、多方向逐段迭代法在不存在局部极小的情况下能够获得与反向追踪法相接近的精度，为了提高改进后的逐段迭代法初始射线精度，以解决局部极小问题，同时进一步提高收敛速度，本发明将反向追踪法的结果作为多尺度、多方向逐段迭代法的初始值，提出了一种基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，这种新的联合射线追踪方法可被认为是利用多尺度、多方向逐段迭代法对反向追踪法追踪结果进行二次矫正。
[0106]
如图14所示，基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法的具体流程包括以下步骤：
[0107]
(1)首先，获得常规反向追踪法的追踪结果；
[0108]
(2-1)将常规反向追踪法的追踪结果作为逐段迭代法的初始射线及追踪点；
[0109]
(2-2)设定抽稀级别为n并将原始网格抽稀到最稀疏级别；
[0110]
(2-3)设定迭代次数为m，若满足m＞0，则进行步骤(2-4)；否则，则进行步骤(4-1)
[0111]
(2-4)迭代次数m成为m-1；进行(3-1)
[0112]
(3-1)对抽稀网格应用常规逐段迭代法并跳过不满足斯奈尔定律的点；
[0113]
(3-2)将抽稀网格和其中的追踪点绕原点顺时针旋转0-90
°
范围内任意角度；
[0114]
(3-3)用与抽稀网格相同网格间距的新网格覆盖旋转后的抽稀网格；
[0115]
(3-4)利用插值计算新网格中的追踪点(旋转后的射线与新网格边界的交点)；
[0116]
(3-5)将新网格和其中的追踪点平移到非负区域进行逐段迭代更新并跳过不满足斯奈尔定律的点；
[0117]
(3-6)将新网格及对应的追踪点变换回原始抽稀网格及对应的追踪点；
[0118]
(3-7)设定n＝0，若是，则进行步骤2-3；否则，则进行(3-8)；
[0119]
(3-8)设定n＝n-1，对当前抽稀网格进行一次加密(加入的网格线位于相邻两条网格线的中间)并利用插值求得新的追踪点；继续进行步骤(2-3)；
[0120]
(4-1)设定n＝0，若否，则进行步骤(4-2)；若是，则进行步骤(4-5)；
[0121]
(4-2)设定抽稀级别n为n-1；继续(4-3)
[0122]
(4-3)对当前抽稀网格进行一次加密(加入的网格线位于相邻两条网格线的中间)并利用插值求得新的追踪点；
[0123]
(4-4)对当前抽稀级别的网格上进行常规逐段迭代法；继续进行步骤(4-1)；
[0124]
(4-5)进行常规逐段迭代，
[0125]
(4-6)输出追踪结果；
[0126]
(4-7)结束。
[0127]
利用基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法对图7所示模型(同为图15(a))进行数值模拟，结果如图15(b、c)所示。其结果显示，联合射线追踪法避开了局部极小而追踪到了折射波路径，在射线水平段时反向追踪法出现若干小幅度的弯折，而联合射线追踪法得到了明显改善。为了确保避开局部极小的原因来自于联合射线追踪法，图15(b)中还给出使用多方向、多尺度逐段迭代法的结果，该追踪结果仍与初始值没有差异，故认为该方法追踪失败。为了对基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法的结果进行误差分析和统计，本发明以反向追踪法的数值模拟结果初始化射线和追踪点，设定参数为45度旋转，5级抽稀以及每级迭代10次，对图3、图4所示模型进行数值模拟。
[0128]
联合射线追踪法的追踪结果如图16(a)所示，相比于反向追踪法的数值模拟结果如图3(c)所示，联合射线追踪法的结果即使将震源区域局部放大，如图16(b)所示，也与解析解拟合良好，有效降低了追踪误差，精度得到了明显提升。联合射线追踪法的单射线整体误差统计如图16(c)所示，相比于反向追踪法如图4(c)所示，500条射线的单射线整体误差的绝对大小降低到5m以下，误差平均值降低到0.86m，单射线整体追踪结果精度提升了90.65％。联合射线追踪法的单射线分布误差统计如图16(d)所示，绝大多数区域误差接近0，有极少区域误差绝对大小超过5m，基本不存在差绝对大小在10m以上的区域。通过去上述误差结果分析，验证了基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法将常规反向追踪法的精度提高了一个数量级左右。
[0129]
综上所述，本发明为了解决常规逐段迭代法在矩形网格中应用时所存在的三个问题：收敛速度过慢、可能陷入局部极小、无法追踪回折波，本发明提出了具有针对性的三种方法或者策略：(1)提出多尺度方法，先利用抽稀的模型来构建大致的射线路径，再对模型进行加密，以刻画细节，这有效提高了逐段迭代法的收敛速度；(2)提出多方向方法，通过改变网格边界的方向，将回折波等在原始网格中不满足斯奈尔定律的射线调整到了可以利用逐段迭代法进行更新的状态；(3)将反向追踪法的结果作为逐段迭代法的初始值，避免了迭代过程陷入局部极小。从反向追踪法的角度出发，也可以认为本发明是将逐段迭代法作为一个后处理流程，对反向追踪法的结果进行进一步的矫正。因此，本发明有效克服了现有技术中的种种缺点而具高度产业利用价值。
[0130]
以上仅为本发明的优选实施方式而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：

1.基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，其包括以下步骤：s1、将常规反向追踪法的追踪结果作为逐段迭代法的初始值；s2、将原始网格抽稀到预设的最稀疏级别，同时只保留剩余原始网格边界上的路径点；s3、对抽稀网格应用常规逐段迭代，并改变抽稀网格边界的方向，通过改变空间离散化的方式，使得回折路径变为满足斯奈尔定律的路径。s4、对抽稀网格逐渐加密，每次加密重复步骤s3，直到当前的抽稀级别降为0，再进行最后一遍常规逐段迭代。2.根据权利要求1所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，所述步骤s2具体包括以下流程：2-1在原始网格中，给定初始射线和追踪点；2-2设定抽稀级别为n并将网格抽稀到最稀疏级别，只保留剩余网格边界上的路径点。3.根据权利要求2所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，在步骤2-2中，抽稀方式为每多一级抽稀就对当前网格模型中隔一条线保留一条，n级抽稀相当于每2
n
条线保留一条。4.根据权利要求1所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，所述步骤s3具体包括以下流程：3-1在抽稀网格中进行逐段迭代，并跳过不满足斯奈尔定律的追踪点；3-2将抽稀网格和其中的追踪点绕原点顺时针旋转；3-3利用与抽稀网格相同新网格间距的水平-垂直网格覆盖旋转后的抽稀网格；3-4利用插值，计算旋转后的射线与新网格边界的交点，即新网格中的追踪点；3-5将新网格和其中的追踪点平移到非负区域；3-6对新网格中的射线进行逐段迭代更新，同时跳过不满足斯奈尔定律的点；3-7将新网格及对应的追踪点变换回原始网格及对应的追踪点。5.根据权利要求4所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，在步骤3-2中，旋转的角度范围为0-90
°
。6.根据权利要求4所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，在步骤3-6中，由于网格边界和射线的相对方向发生改变，新网格中不满足斯奈尔定律的点与步骤3-1不同。7.根据权利要求1所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，所述步骤s4对抽稀网格逐渐加密成新网格并利于插值求得新的追踪点，同时添加对应新网格边界上的路径点，并在每一个稀疏级别上都进行步骤s3。8.根据权利要求7所述的基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其特征在于，所述步骤s4具体包括以下流程：4-1对当前抽稀网格进行一次加密，加入的网格线位于相邻两条网格线的中间，然后利用插值求得新的追踪点；4-2在当前抽稀级别的网格上进行步骤s3；4-3重复步骤4-1和4-2，直到当前的抽稀级别降为0，再进行最后一遍常规逐段迭代。

技术总结

本发明提供了一种基于多尺度、多方向逐段迭代和反向追踪的联合射线追踪方法，其包括以下步骤：S1、将常规反向追踪法的追踪结果作为逐段迭代法的初始值；S2、将原始网格抽稀到预设的最稀疏级别，同时只保留剩余原始网格边界上的路径点；S3、对抽稀网格应用常规逐段迭代，并改变抽稀网格边界的方向，通过改变空间离散化的方式，使得回折路径变为满足斯奈尔定律的路径；S4、对抽稀网格逐渐加密，每次加密重复步骤S3，直到当前的抽稀级别降为0，再进行最后一遍常规逐段迭代。本发明解决了常规逐段迭代法存在收敛过慢、无法追踪回折射线路径和局部收敛极小的问题，本发明既能够适应复杂的速度场变化，也能够在震源附近获得较高的精度。也能够在震源附近获得较高的精度。也能够在震源附近获得较高的精度。