CN114429047A - 一种基于三角网格的二次圆方程旅行时插值方法 - Google Patents
一种基于三角网格的二次圆方程旅行时插值方法 Download PDFInfo
- Publication number
- CN114429047A CN114429047A CN202210100554.3A CN202210100554A CN114429047A CN 114429047 A CN114429047 A CN 114429047A CN 202210100554 A CN202210100554 A CN 202210100554A CN 114429047 A CN114429047 A CN 114429047A
- Authority
- CN
- China
- Prior art keywords
- point
- vertex
- travel time
- equation
- triangular mesh
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000004364 calculation method Methods 0.000 claims abstract description 18
- 230000005284 excitation Effects 0.000 claims description 14
- 238000010276 construction Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Acoustics & Sound (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Geophysics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Geophysics And Detection Of Objects (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于三角网格的二次圆方程旅行时插值方法,通过设定当前所遍历的三角形(其中的两顶点A,B为已知点,第三个顶点C的坐标为待插值,点D为AB边上任意一点,点M为AB边的中点),并根据两顶点A、B计算其中点旅行时,然后结合两顶点以及中点坐标和旅行时,构造当前遍历三角形的圆方程函数,然后结合三角形内部旅行时表达式,并根据费马原理采用Ferrari法求解一元四次方程,即可获得当前遍历三角形顶点C的最短旅行时。本发明从波的传播特点出发,通过两个相邻节点之间的中心点构造辅助圆,并构造非线性旅行时插值,相比传统的LTI等方法可以获得更高的数值精度,具有更少的计算成本。
Description
技术领域
本发明涉及地球物理技术领域,具体涉及的是一种基于三角网格的二次圆方程旅行时插值方法。
背景技术
随着计算机技术的快速发展和进步,人们对地球内部的结构了解越发深入,层析成像技术也成为了进行研究的有力工具,而射线追踪方法则是地学层析成像中必不可少的方法之一。射线追踪方法是基于斯奈尔定律和惠更斯原理,来模拟地震波的传播运动学特征,其本质是在给定的震源激发点和接收点之间的两点射线追踪问题。
随着更加深入的研究,涌现了很多改进的新型算法。Sethian和Popovici(1999)提出了一种快速行进法(FMM),采用迎风差分格式求解局部eikonal方程。Asakawa和Kawanaka(1993)提出了线性旅行时插值(LTI)方法。其中,FMM方法因为其基于差分格式进行运算,所以适用于规则的矩形网格,易于计算高阶差分以提高其自身精度。但由于实际情况下地表起伏较大,矩形网剖分误差较大,边缘锯齿化,并不能很好的贴合外轮廓,因而不能灵活处理因复杂地形引起的不规则计算边界问题。
LTI是一种线性旅行时插值算法,其更容易适用于三角网格上进行计算,当实际情况下地表起伏较大时,三角网格可以更好的贴合地表外轮廓。由于射线追踪的线性旅行时插值对于快速模拟旅行时很重要,因此LTI在地球物理中也得到了广泛的应用。
目前,将LTI方法应用至震源激发点上的主要流程如下:
A1、首先计算包含震源激发点的三角网格的三个顶点的旅行;将三个顶点标记为固定点,并且将它们加入到可到达点(已经计算过旅行时的网格结点,但不确定求得了最小旅行时)表Q;
A2、判断表Q是否为空,若为空,算法结束;若不为空,则从表Q中找出旅行时最小的点Pi,判断Pi是否为固定点(已经确定求得了最小旅行时的网格结点);
A3、从表Q中找出旅行时最小的点Pi,若Pi不是固定点,将其标记为固定点;
A4、以Pi为子震源激发点(在每一轮计算中,从当前所有可到达的结点中找出的、未做过子震源且旅行时最小的结点),遍历所有与Pi点共边的邻接点Pj;
A5、若Pj是固定点,且ΔPiPjPk的另一点Pk不是固定点,则利用线性旅行时插值公式(LTI),以Pi—Pj为扩展边,插值计算Pk点的旅行时tk;
A6、如果tk小于Pk点已有旅行时Tk,则令Tk=tk;
A7、如果Pk原来不是可到达点,则将Pk加入到表Q中;
A8、将Pi点从表Q中去除,返回步骤A2。
LTI方法虽然可以用于三角网格插值,但因其是一种线性插值方法,因此在靠近震源激发点处,波前面呈球面型,即近场时误差较大,常常使得计算精度不高。
发明内容
本发明的目的在于提供一种基于三角网格的二次圆方程旅行时插值方法,用以替代现有的LTI线性旅行时插值方法,解决LTI算法中存在的近场时计算误差大、精度低且计算效率较低的问题。
为实现上述目的,本发明采用的技术方案如下:
一种基于三角网格的二次圆方程旅行时插值方法,包括以下步骤:
S1、设定一个包含有震源激发点的三角网格,其中,三角网格中任意两个顶点A、B的坐标为已知值,第三个顶点C的坐标为待插值,点D为AB边上任意一点,点M为AB边的中点;
S2、根据点A、B、M的坐标和旅行时构造当前三角网格的圆方程函数,用以表达点D的旅行时;
S3、根据点D的履行时,结合三角形函数,确定点D到顶点C的距离,最终根据费马原理并采用Ferrari法求解一元四次方程,即可插值计算顶点C的最短履行时。
具体地,所述步骤S2中,圆方程函数的构造过程如下:
(a)以三角网格中的顶点A为坐标原点建立坐标系,且该顶点A旅行时为纵轴,则该顶点A的坐标为A'(0,ta),三角网格中另一顶点B的坐标为B'(|AB|,tb);
(c)根据二维最大外接圆计算公式,确定三角形A'B'M'的最大外接圆O的圆心坐标(cx,cy),以及圆半径cr,从而构造出以下圆方程函数:
公式(1)中,x表示点D到A'的距离,即原坐标系中电D到顶点A的距离;td表示点D的旅行时。
具体地,所述步骤S3包括以下步骤:
(d)根据下列公式计算三角网格中第三个顶点C的旅行时:
(e)根据费马原理,将公式(2)转换为:
并进一步转换为关于x的四次方程:
ax4+bx3+cx2+dx+e=0 (4);
式(4)中系数a、b、c、d、e由下式表达;|AC|为AC边长;
(f)采用Ferrari法求解x,并将其代入公式(2)中的x,获得顶点C的最短旅行时tc。
与现有技术相比,本发明具有以下有益效果:
本发明采用圆弧的形式来拟合地震波前面,即通过三角网格中相邻两点及其中点构造出最大外接圆,并使用外接圆的下半圆弧来近似表示波前面。如此,本发明从波的传播特点出发,可以获得比LTI等方法更高的计算精度以及更少的计算成本,与以往的非线性插值方法相比,本发明具有更好的表现。模型数值实验也表明,本发明旅行时误差较小,不超过其他两种现有方法的十分之一(通常小于4%),这将非常有助于提高计算机层析成像和旅行时偏移的精度。
附图说明
图1为本发明-实施例的流程示意图。
图2为本发明-实施例中当前所遍历三角形的初始坐标系示意图。
图3为本发明-实施例中当前所遍历三角形的新坐标系示意图。
图4为本发明-实施例中采用的模型结构示意图。
图5为本发明-实施例中理论路径和理论接收点位置示意图。
图6为本发明-实施例中基于圆方程函数的计算结果示意图。
具体实施方式
本发明提供了一种新型的旅行时插值方法,应用于三角网格的计算上,可替代现有的LTI线性旅行时插值方法。本发明是在LTI等方法基础上改进的一种旅行时插值方法,可以实现高精度的旅行时计算。该方法的核心是通过相邻节点之间的辅助点对每个单元的旅行时间进行非线性插值。下面结合附图说明和实施例对本发明作进一步说明,本发明的实施包含但不限于以下实施例。
实施例
本实施例的设计思路在于,通过两个相邻节点之间的中心点构造辅助圆,构造非线性旅行时插值。本实施例的插值流程如图1所示,具体来说,其流程如下所述:
S1、设定当前所遍历的包含有震源激发点的三角网格,其中,三角网格中任意两个顶点A、B的坐标为已知值,第三个顶点C的坐标为待插值,点D为AB边上任意一点,点M为AB边的中点,如图2所示;
S2、以三角网格中的顶点A为坐标原点建立坐标系,且该顶点A旅行时为纵轴,则该顶点A的坐标为A'(0,ta),三角网格中另一顶点B的坐标为B'(AB,tb);
(c)根据二维最大外接圆计算公式,确定三角形A'B'M'的最大外接圆O的圆心坐标(cx,cy),以及圆半径cr,从而构造出以下圆方程函数:
公式(1)中,x表示点D到A'的距离,即原坐标系中电D到顶点A的距离;td表示点D的旅行时;
构造原理:由于射线追踪法中射线会从底边AB射入并到达顶点C,因此圆方程函数中的x范围为[0,|AB|];而根据波的传播特点,该圆方程只取下半圆,因此点D的旅行时可以表示为:
S5、根据费马原理,通过点D到点C的射线应满足如下要求:
因此,可将公式(2)转换为如下公式:
并进一步转换为关于x的四次方程:
ax4+bx3+cx2+dx+e=0 (5);
式(5)中系数a、b、c、d、e由下式表示;|AC|为AC边长;
S6、最后,采用Ferrari法求解x,并将其代入公式(2)中的x,即可获得顶点C的最短旅行时tc。
应用示例:
将上述插值方法应用到震源激发点中,其流程如下所述:
首先,建立一个长为2000,宽为400的双层地质模型,模型结构如图4所示,其中B1层电阻率为1500,B2层电阻率为2000。震源激发点坐标位于模型中(1000,0)处,在震源激发点所处水平面400处设置9个地震波接收点。地震波从震源激发点出发,以角度为10°的偏差设置9条关于震源激发点对称的追踪射线,并根据波传播理论计算这9条射线的理论传播路径以及9个理论接收点位置,如图5所示。
计算包含震源激发点的三角形三个顶点的旅行。将三个顶点标记为固定点,并且将它们加入到可到达点表Q。
判断表Q是否为空,若为空,则结束;若不为空,则从表Q中找出旅行时最小的点Pi,判断Pi是否为固定点。如果Pi不是固定点,则将其标记为固定点,并以Pi为子震源激发点,遍历所有与Pi点共边的邻接点Pj。
若Pj是固定点,且ΔPiPjPk的另一点Pk不是固定点,则利用本实施例设计的插值方法,以Pi—Pj为扩展边,以Pi、Pj及其中点M构建关于其扩展边上任意一点旅行时的最大外接圆圆方程函数,即上述公式(1)。然后结合三角形内部旅行时表达式求得Pk旅行时表达式,即上述公式(2),再转换为公式(4)。
最后,将公式(4)转换为一元四次方程,即上述公式(5),然后采用Ferrari'smethod(Fujii,2003)方法求解公式(5),即可获得当前遍历三角形顶点Pk的最短旅行时tk。
接着,判断tk与Pk点已有旅行时Tk的大小,如果tk小于Pk点已有旅行时Tk,则令Tk=tk;如果Pk原来不是可到达点,则将Pk加入到表Q中,并将Pi点从表Q中去除。
上述示例的插值效果如图6所示。根据模型数值计算结果可以看到,本发明方案相比传统的LTI等方法具有更高的数值精度和更少的计算成本。作为LTI方法的直接高阶扩展,本发明的方法在建立旅行时计算公式方面有完全不同的形式。因此,与以往的非线性插值方法相比,该方法具有更好的表现。
上述实施例仅为本发明的优选实施方式,不应当用于限制本发明的保护范围,凡在本发明的主体设计思想和精神上作出的毫无实质意义的改动或润色,其所解决的技术问题仍然与本发明一致的,均应当包含在本发明的保护范围之内。
Claims (3)
1.一种基于三角网格的二次圆方程旅行时插值方法,其特征在于,包括以下步骤:
S1、设定一个包含有震源激发点的三角网格,其中,三角网格中任意两个顶点A、B的坐标为已知值,第三个顶点C的坐标为待插值,点D为AB边上任意一点,点M为AB边的中点;
S2、根据点A、B、M的坐标和旅行时构造当前三角网格的圆方程函数,用以表达点D的旅行时;
S3、根据点D的履行时,结合三角形函数,确定点D到顶点C的距离,最终根据费马原理并采用Ferrari法求解一元四次方程,即可插值计算顶点C的最短履行时。
2.根据权利要求1所述的一种基于三角网格的二次圆方程旅行时插值方法,其特征在于,所述步骤S2中,圆方程函数的构造过程如下:
(a)以三角网格中的顶点A为坐标原点建立坐标系,且该顶点A旅行时为纵轴,则该顶点A的坐标为A'(0,ta),三角网格中另一顶点B的坐标为B'(|AB|,tb);
(c)根据二维最大外接圆计算公式,确定三角形A'B'M'的最大外接圆O的圆心坐标(cx,cy),以及圆半径cr,从而构造出以下圆方程函数:
公式(1)中,x表示点D到A'的距离,即原坐标系中电D到顶点A的距离;td表示点D的旅行时。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210100554.3A CN114429047B (zh) | 2022-01-27 | 2022-01-27 | 一种基于三角网格的二次圆方程旅行时插值方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210100554.3A CN114429047B (zh) | 2022-01-27 | 2022-01-27 | 一种基于三角网格的二次圆方程旅行时插值方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114429047A true CN114429047A (zh) | 2022-05-03 |
CN114429047B CN114429047B (zh) | 2023-08-22 |
Family
ID=81312961
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210100554.3A Active CN114429047B (zh) | 2022-01-27 | 2022-01-27 | 一种基于三角网格的二次圆方程旅行时插值方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114429047B (zh) |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20040122594A1 (en) * | 2002-12-23 | 2004-06-24 | Toshifumi Matsuoka | Methods for determining formation and borehole parameters using fresnel volume tomography |
CN1712991A (zh) * | 2004-06-25 | 2005-12-28 | 中国石油化工股份有限公司 | 一种用于地震勘探中射线追踪的方法 |
CN101533102A (zh) * | 2009-04-09 | 2009-09-16 | 长江工程地球物理勘测武汉有限公司 | 二维复杂结构三角网射线追踪全局方法 |
CN102830431A (zh) * | 2012-08-14 | 2012-12-19 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 真地表射线追踪自适应插值方法 |
CN108267781A (zh) * | 2017-12-15 | 2018-07-10 | 桂林理工大学 | 任意曲面非均匀介质快速行进程函方程求解射线追踪算法 |
CN108986218A (zh) * | 2018-06-06 | 2018-12-11 | 东南大学 | 一种基于pmvs的建筑物密集点云快速重建方法 |
CN112257241A (zh) * | 2020-10-15 | 2021-01-22 | 成都理工大学 | 一种三角网菲涅尔带时差层析反演方法 |
CN112596103A (zh) * | 2020-11-24 | 2021-04-02 | 中国地质科学院地球物理地球化学勘查研究所 | 射线追踪方法、装置和电子设备 |
-
2022
- 2022-01-27 CN CN202210100554.3A patent/CN114429047B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20040122594A1 (en) * | 2002-12-23 | 2004-06-24 | Toshifumi Matsuoka | Methods for determining formation and borehole parameters using fresnel volume tomography |
CN1712991A (zh) * | 2004-06-25 | 2005-12-28 | 中国石油化工股份有限公司 | 一种用于地震勘探中射线追踪的方法 |
CN101533102A (zh) * | 2009-04-09 | 2009-09-16 | 长江工程地球物理勘测武汉有限公司 | 二维复杂结构三角网射线追踪全局方法 |
CN102830431A (zh) * | 2012-08-14 | 2012-12-19 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 真地表射线追踪自适应插值方法 |
CN108267781A (zh) * | 2017-12-15 | 2018-07-10 | 桂林理工大学 | 任意曲面非均匀介质快速行进程函方程求解射线追踪算法 |
CN108986218A (zh) * | 2018-06-06 | 2018-12-11 | 东南大学 | 一种基于pmvs的建筑物密集点云快速重建方法 |
CN112257241A (zh) * | 2020-10-15 | 2021-01-22 | 成都理工大学 | 一种三角网菲涅尔带时差层析反演方法 |
CN112596103A (zh) * | 2020-11-24 | 2021-04-02 | 中国地质科学院地球物理地球化学勘查研究所 | 射线追踪方法、装置和电子设备 |
Non-Patent Citations (3)
Title |
---|
FRANCESC SORIGUERA等: "Requiem for Freeway Travel Time Estimation Methods Based on Blind Speed Interpolations Between Point Measurements", pages 1 - 7 * |
宋御杰: "地震波旅行时线性和非线性插值计算方法", pages 751 * |
王琦等: "VTI介质起伏界面混合网格旅行时线性插值计算方法", vol. 53, no. 6, pages 1175 - 1187 * |
Also Published As
Publication number | Publication date |
---|---|
CN114429047B (zh) | 2023-08-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Talmor | Well-spaced points for numerical methods | |
Lan et al. | A high‐order fast‐sweeping scheme for calculating first‐arrival travel times with an irregular surface | |
CN109636912B (zh) | 应用于三维声呐图像重构的四面体剖分有限元插值方法 | |
WO2019242045A9 (zh) | 一种虚拟震源二维波前构建地震波走时计算方法 | |
CN108108579B (zh) | 直流电阻率无单元法中耦合有限单元法的边界处理方法 | |
Le Bouteiller et al. | A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media | |
Ma et al. | Calculating Ray Paths for First‐Arrival Travel Times Using a Topography‐Dependent Eikonal Equation Solver | |
Zhang et al. | Eikonal solver in the celerity domain | |
CN116774292B (zh) | 一种地震波走时确定方法、系统、电子设备及存储介质 | |
CN115758938A (zh) | 面向粘性边界流场数值模拟的附面层网格生成方法 | |
JP2011039691A (ja) | メッシュモデル生成装置、プログラム及びメッシュモデル生成方法 | |
Cecil et al. | Simplex free adaptive tree fast sweeping and evolution methods for solving level set equations in arbitrary dimension | |
Acosta | High order surface radiation conditions for time-harmonic waves in exterior domains | |
Lavery | Shape-preserving, multiscale interpolation by bi-and multivariate cubic L1 splines | |
CN105931297A (zh) | 三维地质表面模型中的数据处理方法 | |
CN114429047A (zh) | 一种基于三角网格的二次圆方程旅行时插值方法 | |
Pandey et al. | Improved convergence of fast integral equation solvers for acoustic scattering by inhomogeneous penetrable media with discontinuous material interface | |
CN111158059A (zh) | 基于三次b样条函数的重力反演方法 | |
CN110568497B (zh) | 一种复杂介质条件下地震初至波旅行时的精确求解方法 | |
CN108802819B (zh) | 一种深度均匀采样梯形网格有限差分地震波场模拟方法 | |
CN111797551B (zh) | 一种面向试验应变数据三维云图显示的插值方法 | |
CN111898819B (zh) | 空间网格划分方法及装置 | |
Liu et al. | An improved Poisson surface reconstruction algorithm based on the boundary constraints | |
CN103914431A (zh) | 一种计算各向异性结构雷达横截面的无网格法 | |
CN105869209A (zh) | 三维地质表面模型中的畸形三角形数据处理方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |