一种基于三角网格的二次圆方程旅行时插值方法

1.本发明涉及地球物理技术领域，具体涉及的是一种基于三角网格的二次圆方程旅行时插值方法。


背景技术：

2.随着计算机技术的快速发展和进步，人们对地球内部的结构了解越发深入，层析成像技术也成为了进行研究的有力工具，而射线追踪方法则是地学层析成像中必不可少的方法之一。射线追踪方法是基于斯奈尔定律和惠更斯原理，来模拟地震波的传播运动学特征，其本质是在给定的震源激发点和接收点之间的两点射线追踪问题。
3.随着更加深入的研究，涌现了很多改进的新型算法。sethian和popovici(1999)提出了一种快速行进法(fmm)，采用迎风差分格式求解局部eikonal方程。asakawa和kawanaka(1993)提出了线性旅行时插值(lti)方法。其中，fmm方法因为其基于差分格式进行运算，所以适用于规则的矩形网格，易于计算高阶差分以提高其自身精度。但由于实际情况下地表起伏较大，矩形网剖分误差较大，边缘锯齿化，并不能很好的贴合外轮廓，因而不能灵活处理因复杂地形引起的不规则计算边界问题。
4.lti是一种线性旅行时插值算法，其更容易适用于三角网格上进行计算，当实际情况下地表起伏较大时，三角网格可以更好的贴合地表外轮廓。由于射线追踪的线性旅行时插值对于快速模拟旅行时很重要，因此lti在地球物理中也得到了广泛的应用。
5.目前，将lti方法应用至震源激发点上的主要流程如下：
6.a1、首先计算包含震源激发点的三角网格的三个顶点的旅行；将三个顶点标记为固定点，并且将它们加入到可到达点(已经计算过旅行时的网格结点，但不确定求得了最小旅行时)表q；
7.a2、判断表q是否为空，若为空，算法结束；若不为空，则从表q中找出旅行时最小的点pi，判断pi是否为固定点(已经确定求得了最小旅行时的网格结点)；
8.a3、从表q中找出旅行时最小的点pi，若pi不是固定点，将其标记为固定点；
9.a4、以pi为子震源激发点(在每一轮计算中，从当前所有可到达的结点中找出的、未做过子震源且旅行时最小的结点)，遍历所有与pi点共边的邻接点pj；
10.a5、若pj是固定点，且δpipjpk的另一点pk不是固定点，则利用线性旅行时插值公式(lti)，以pi—pj为扩展边，插值计算pk点的旅行时tk；
11.a6、如果tk小于pk点已有旅行时tk，则令tk＝tk；
12.a7、如果pk原来不是可到达点，则将pk加入到表q中；
13.a8、将pi点从表q中去除，返回步骤a2。
14.lti方法虽然可以用于三角网格插值，但因其是一种线性插值方法，因此在靠近震源激发点处，波前面呈球面型，即近场时误差较大，常常使得计算精度不高。


技术实现要素：

15.本发明的目的在于提供一种基于三角网格的二次圆方程旅行时插值方法，用以替代现有的lti线性旅行时插值方法，解决lti算法中存在的近场时计算误差大、精度低且计算效率较低的问题。
16.为实现上述目的，本发明采用的技术方案如下：
17.一种基于三角网格的二次圆方程旅行时插值方法，包括以下步骤：
18.s1、设定一个包含有震源激发点的三角网格，其中，三角网格中任意两个顶点a、b的坐标为已知值，第三个顶点c的坐标为待插值，点d为ab边上任意一点，点m为ab边的中点；
19.s2、根据点a、b、m的坐标和旅行时构造当前三角网格的圆方程函数，用以表达点d的旅行时；
20.s3、根据点d的履行时，结合三角形函数，确定点d到顶点c的距离，最终根据费马原理并采用ferrari法求解一元四次方程，即可插值计算顶点c的最短履行时。
21.具体地，所述步骤s2中，圆方程函数的构造过程如下：
22.(a)以三角网格中的顶点a为坐标原点建立坐标系，且该顶点a旅行时为纵轴，则该顶点a的坐标为a'(0,ta)，三角网格中另一顶点b的坐标为b'(|ab|,tb)；
23.(b)根据顶点a、b的坐标，确定三角网格中ab边的中点m的坐标为且点m的旅行时为：
24.(c)根据二维最大外接圆计算公式，确定三角形a'b'm'的最大外接圆o的圆心坐标(c
x
,cy)，以及圆半径cr，从而构造出以下圆方程函数：
[0025][0026]
公式(1)中，x表示点d到a'的距离，即原坐标系中电d到顶点a的距离；td表示点d的旅行时。
[0027]
具体地，所述步骤s3包括以下步骤：
[0028]
(d)根据下列公式计算三角网格中第三个顶点c的旅行时：
[0029][0030]
公式(2)中，表示点d到顶点c的距离，根据公式(1)、x及三角函数得出；α为边ab与ac的夹角；s为慢度；
[0031]
(e)根据费马原理，将公式(2)转换为：
[0032][0033]
并进一步转换为关于x的四次方程：
[0034]
ax4 bx3 cx2 dx e＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)；
[0035]
式(4)中系数a、b、c、d、e由下式表达；|ac|为ac边长；
[0036][0037]
(f)采用ferrari法求解x，并将其代入公式(2)中的x，获得顶点c的最短旅行时tc。
[0038]
与现有技术相比，本发明具有以下有益效果：
[0039]
本发明采用圆弧的形式来拟合地震波前面，即通过三角网格中相邻两点及其中点构造出最大外接圆，并使用外接圆的下半圆弧来近似表示波前面。如此，本发明从波的传播特点出发，可以获得比lti等方法更高的计算精度以及更少的计算成本，与以往的非线性插值方法相比，本发明具有更好的表现。模型数值实验也表明，本发明旅行时误差较小，不超过其他两种现有方法的十分之一(通常小于4％)，这将非常有助于提高计算机层析成像和旅行时偏移的精度。
附图说明
[0040]
图1为本发明-实施例的流程示意图。
[0041]
图2为本发明-实施例中当前所遍历三角形的初始坐标系示意图。
[0042]
图3为本发明-实施例中当前所遍历三角形的新坐标系示意图。
[0043]
图4为本发明-实施例中采用的模型结构示意图。
[0044]
图5为本发明-实施例中理论路径和理论接收点位置示意图。
[0045]
图6为本发明-实施例中基于圆方程函数的计算结果示意图。
具体实施方式
[0046]
本发明提供了一种新型的旅行时插值方法，应用于三角网格的计算上，可替代现有的lti线性旅行时插值方法。本发明是在lti等方法基础上改进的一种旅行时插值方法，可以实现高精度的旅行时计算。该方法的核心是通过相邻节点之间的辅助点对每个单元的旅行时间进行非线性插值。下面结合附图说明和实施例对本发明作进一步说明，本发明的实施包含但不限于以下实施例。
[0047]
实施例
[0048]
本实施例的设计思路在于，通过两个相邻节点之间的中心点构造辅助圆，构造非线性旅行时插值。本实施例的插值流程如图1所示，具体来说，其流程如下所述：
[0049]
s1、设定当前所遍历的包含有震源激发点的三角网格，其中，三角网格中任意两个顶点a、b的坐标为已知值，第三个顶点c的坐标为待插值，点d为ab边上任意一点，点m为ab边的中点，如图2所示；
[0050]
s2、以三角网格中的顶点a为坐标原点建立坐标系，且该顶点a旅行时为纵轴，则该顶点a的坐标为a'(0,ta)，三角网格中另一顶点b的坐标为b'(ab,tb)；
[0051]
s3、根据顶点a、b的坐标，确定三角网格中ab边的中点m的坐标为且点m的旅行时为：如图3所示；
[0052]
(c)根据二维最大外接圆计算公式，确定三角形a'b'm'的最大外接圆o的圆心坐标(c
x
,cy)，以及圆半径cr，从而构造出以下圆方程函数：
[0053][0054]
公式(1)中，x表示点d到a'的距离，即原坐标系中电d到顶点a的距离；td表示点d的旅行时；
[0055]
构造原理：由于射线追踪法中射线会从底边ab射入并到达顶点c，因此圆方程函数中的x范围为[0,|ab|]；而根据波的传播特点，该圆方程只取下半圆，因此点d的旅行时可以表示为：
[0056][0057]
s4、基于点d旅行时的表达式，以及点a和点d的距离x，使用三角函数可得点d到顶点c的距离为
[0058]
假设慢度(与速度相对，是速度的倒数)为s，α为边ab与ac的夹角，则点d到三角形顶点c的旅行时为结合点d旅行时，顶点c的旅行时可以表示为：
[0059][0060]
s5、根据费马原理，通过点d到点c的射线应满足如下要求：
[0061][0062]
因此，可将公式(2)转换为如下公式：
[0063][0064]
并进一步转换为关于x的四次方程：
[0065]
ax4 bx3 cx2 dx e＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)；
[0066]
式(5)中系数a、b、c、d、e由下式表示；|ac|为ac边长；
[0067][0068]
s6、最后，采用ferrari法求解x，并将其代入公式(2)中的x，即可获得顶点c的最短旅行时tc。
[0069]
应用示例：
[0070]
将上述插值方法应用到震源激发点中，其流程如下所述：
[0071]
首先，建立一个长为2000，宽为400的双层地质模型，模型结构如图4所示，其中b1层电阻率为1500，b2层电阻率为2000。震源激发点坐标位于模型中(1000，0)处，在震源激发点所处水平面400处设置9个地震波接收点。地震波从震源激发点出发，以角度为10
°
的偏差设置9条关于震源激发点对称的追踪射线，并根据波传播理论计算这9条射线的理论传播路
径以及9个理论接收点位置，如图5所示。
[0072]
计算包含震源激发点的三角形三个顶点的旅行。将三个顶点标记为固定点，并且将它们加入到可到达点表q。
[0073]
判断表q是否为空，若为空，则结束；若不为空，则从表q中找出旅行时最小的点pi，判断pi是否为固定点。如果pi不是固定点，则将其标记为固定点，并以pi为子震源激发点，遍历所有与pi点共边的邻接点pj。
[0074]
若pj是固定点，且δpipjpk的另一点pk不是固定点，则利用本实施例设计的插值方法，以pi—pj为扩展边，以pi、pj及其中点m构建关于其扩展边上任意一点旅行时的最大外接圆圆方程函数，即上述公式(1)。然后结合三角形内部旅行时表达式求得pk旅行时表达式，即上述公式(2)，再转换为公式(4)。
[0075]
最后，将公式(4)转换为一元四次方程，即上述公式(5)，然后采用ferrari's method(fujii,2003)方法求解公式(5)，即可获得当前遍历三角形顶点pk的最短旅行时tk。
[0076]
接着，判断tk与pk点已有旅行时tk的大小，如果tk小于pk点已有旅行时tk，则令tk＝tk；如果pk原来不是可到达点，则将pk加入到表q中，并将pi点从表q中去除。
[0077]
上述示例的插值效果如图6所示。根据模型数值计算结果可以看到，本发明方案相比传统的lti等方法具有更高的数值精度和更少的计算成本。作为lti方法的直接高阶扩展，本发明的方法在建立旅行时计算公式方面有完全不同的形式。因此，与以往的非线性插值方法相比，该方法具有更好的表现。
[0078]
上述实施例仅为本发明的优选实施方式，不应当用于限制本发明的保护范围，凡在本发明的主体设计思想和精神上作出的毫无实质意义的改动或润色，其所解决的技术问题仍然与本发明一致的，均应当包含在本发明的保护范围之内。