CN108267781A - 任意曲面非均匀介质快速行进程函方程求解射线追踪算法 - Google Patents

任意曲面非均匀介质快速行进程函方程求解射线追踪算法 Download PDF

Info

Publication number
CN108267781A
CN108267781A CN201711354734.XA CN201711354734A CN108267781A CN 108267781 A CN108267781 A CN 108267781A CN 201711354734 A CN201711354734 A CN 201711354734A CN 108267781 A CN108267781 A CN 108267781A
Authority
CN
China
Prior art keywords
grid
mesh
arbitrary
focus
ray
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
Application number
CN201711354734.XA
Other languages
English (en)
Other versions
CN108267781B (zh
Inventor
薛霆虓
黄明曦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guilin University of Technology
Original Assignee
Guilin University of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Guilin University of Technology filed Critical Guilin University of Technology
Priority to CN201711354734.XA priority Critical patent/CN108267781B/zh
Publication of CN108267781A publication Critical patent/CN108267781A/zh
Application granted granted Critical
Publication of CN108267781B publication Critical patent/CN108267781B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/624Reservoir parameters

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Image Generation (AREA)

Abstract

本发明公开了在任意曲面非均匀介质快速行进法(FMM)程函方程求解射线追踪的算法。快速行进法在体波中的应用得到了很大的发展,在二维平面的面波射线追踪也得到了很好的使用,但是至今未能将FMM方法程函方程运用在带地形任意曲面中。传统的FMM方法在非正交坐标网格中的算法只有一阶精度,本发明的算法利用了上风网格中的方向导数,使得FMM在非正交网格的计算有更高的精度。此发明提出了在任意曲面、非均匀介质追踪的FMM程函方程算法,解决了在任意复杂地形、非均匀介质快速射线追踪的问题。此发明可用于任意起伏地形的面波勘探和面波层析成像。

Description

任意曲面非均匀介质快速行进程函方程求解射线追踪算法
技术领域:
本发明涉及任意曲面射线追踪,射线追踪是地震勘探领域和地震层析成像领域中的重要技术。
背景技术:
面波勘探和面波层析成像是地震勘探领域中要技术,不仅是研究大尺度地球科学中的重要手段,也在浅层地震勘探中得到了广泛地应用。地震勘探属于地球物理反演技术,其中效率最高的是走时反演方法。无论是体波勘探和面波勘探,都离不开正演过程,而正演过程中最关键的步骤是射线追踪。射线追踪技术非常丰富,包含无网格法和网格法。其中无网格法包含弯曲法、试射法、波前重建法等技术;而网格法包含了最短路径法、插值法和快速行进法等等。而快速行进法是最近发展起来的最为先进的方法,具有计算精确、算法无条件稳定等优点。快速行进法在体波层析成像中得到了很好的发展和应用(Rawlinson,etal.,2006;Rawlinson&Urvoy,2006等)。
面波层析成像技术来说,传统的方法是计算源和台站之间的大圆路径,既没有考虑面波速度造成的路径弯曲,更没有考虑研究区域地形的起伏变化。如果要考虑地形对面波传播的影响,根据费马原理,需要在曲面进行射线追踪。最易实现的方法是最短路径法,但是最短路径方法精度较差,计算效率也不够高。而FMM方法效率较高,而且精度更高。比如,文献(Saygin,E.2007.Seismic receiver and noise correlation based studies inAustralia,PhD thesis,Australian National University.)考虑了非均匀介质的影响,用FMM方法求解,得到了很好的结果。但是FMM方法是直接在计算程函方程,而程函方程的求解一般都是在正交坐标系下实现。而要实现曲面的射线追踪,不可能完全用正交坐标系。为此,需要发展出一种方法,能够在非正交坐标系中实现FMM程函方程算法,以便能较为精确地计算在复杂地形和非均匀介质中的射线追踪。
发明内容:
本发明完成了FMM算法求解程函方程不能在非正交坐标系下实现的困难,提出了解决这种困难的方法。实现了复杂地形和非均匀介质的射线追踪,为复杂地形和复杂介质的面波层析成像提供了正演计算。
为解决上述技术问题,本发明采用了以下技术方案:
步骤一:建立速度网格模型(Velocity grids),对曲面进行网格划分,给出每个网格点的速度值,划分得到的网格为三角形网格。一般来说,速度网格和传播网格(Propagation grids)相比较,划分较为粗略。
步骤二:首先建立起伏地形模型,进行传播网格化剖分,划分为三角形网格,划分三角形网格即为建立锐角三角形,具体分为在不考虑地表起伏的水平面建立等边三角形;而若考虑到地形的起伏,所建立的三角形内角在60度附近,为了保证算法的效率,对地形起伏进行一定程度的光滑处理。最后得到有一定光滑度的任意起伏地形模型(图2)。此外,对震源处附近的网格进行细化处理,保证计算的精度。然后根据传播网格周围的速度网格,通过插值方法得到传播网格节点的速度值。为了保证计算的正确,有必要对网格进行检查,并进行适当调整。
步骤三:对每一个三角形网格进行计算。
首先,从震源处所在的网格,开始计算。从震源开始,波前传至周围的网格点。刚开始,震源处的网格点是上风(Upwind)集合,Upwind集合是已经计算好的网格集合,这些正在计算得周围网格点的集合称为窄带(Narrow band),未计算的网格点集合是下风(Downwind)。如此,从震源处开始计算,一直到曲面上所有的网格点都遍历,都成为Upwind集合,如图2所示。
下面介绍每一个三角形网格点的计算算法。
传统的FMM算法都是在正交网格中进行计算程函方程。
式中,T是波前的走时,s(x)是每一个网格点的慢度值。在正交坐标系中可以求解。但是,对于任意曲面剖分成的任意三角型网格,并不是正交网格,因此,可以应用Sethian1999中的办法解决这一问题。
其中各参数见图3。不同于四边形网格,三角形是任意划分的,无法使用有限差分的高阶算法,因此这种方法只有一阶近似,精度不是很高。为了提高计算精度,我们的发明改进了这一算法。具体方法阐述如下:
如图4所示,已经三角形网格点i和j处的波前到时Ti和Tj,要求To的到时。传统方法中,假设了在三角形走时梯度相同,使用了O点的慢度。因此精度较差。为了增加计算精度,在计算Ti和Tj时,可以利用i和j节点处的走时梯度大小(慢度)和在方向导数 如果采取这种算法,需要在计算每个节点走时的同时,也要保存每个节点的走时梯度方向。下面的公式(2)和(3)表达了网格点i、j处a和b方向导数和三角形两条边a和b上的平均方向导数之间的关系。
其中分别表示网格点i、j和o处方向和方向的方向导数,是三角形a和b上的平均方向导数。利用这两个公式可以用已知网格点i和j上的方向导数和三角形边上的平均方向导数得到网格点O处方向的方向导数。
然后,我们就可以建立三角形中FMM算法。
图4显示了如何由网格点O处方向导数,求O处梯度的算法。建立如图4的直角坐标系,所求的O处走时梯度坐标为(x0,y0),那么建立下列方程组:
其中分别是网格点O处方向的方向导数,而且是网格点O到时To的函数(由公式(2)和(3)求得),SO是网格点O处的慢度。将方程组(4)求得的x0和y0代入公式(5)。最后得到:
由此公式求得网格点O处的到时TO。这样就完成了一个三角形网格中的计算。如此,循环下去,一直到所有网格点都变成上风(Upwind)集合。至此所有网格点的到时都已经计算完成。
步骤四:找出震源至台站的射线路径。从台站开始,在每个网格中求梯度,负梯度的方向就是射线方向,把射线沿途经过的网格点中负梯度都求出来。从接受台站开始,连接所有沿途网格中射线的线段,直至震源处,就得到了震源-台站的射线路径。此外,在计算每个网格点的到时的同时,同时计算每个速度网格节点速度的微小变化对走时的影响,也就是进行走时计算的同时,也计算出Frechet导数值,为面波反演作准备。保存走时信息和Frechet导数值。至此,FMM方法的射线追踪就完成了。
与现有技术相比,本发明具有以下优点:
(1)网格剖分灵活方便,无需特殊软件进行剖分。
(2)能够快速实现在非均匀速度介质中、波动在任意曲面的传播路径,同时计算Frechet导数,为面波反演提供基础。
(3)和最短路径算法相比,计算精度更高,计算效率也更高。
(4)本发明继承了FMM算法以下优点:a)算法无条件稳定性;b)满足波前传播规律;c)计算效率高的波前扩展方式。提高传统FMM算法在非正交网格的计算精度。
(5)本发明可以运用于二维曲面的波动研究,也可以用于面波相速度和群速度射线路径计算,继而为带地形面波反演提供高效的正演计算。
附图说明:
1、图1为任意曲面非均匀介质快速行进法射线追踪算法流程图;
2、图2为FMM算法上风、窄带和下风示意图;
3、图3为FMM传统算法在单个三角形网格中实现的示意图;
4、图4为本发明的FMM算法在单个三角形网格中实现的示意图;
具体实施例
步骤一:建立速度网格模型(Velocity grids),对曲面进行网格划分,给出每个网格点的速度值。一般来说,速度网格和传播网格(Propagation grids)相比较,划分较为粗略。
步骤二:首先建立起伏地形模型,进行传播网格化剖分,划分为三角形网格,划分三角形网格即为建立锐角三角形,具体分为在不考虑地表起伏的水平面建立等边三角形;而若考虑到地形的起伏,所建立的三角形内角在60度附近。为了保证算法的效率,对地形起伏进行一定程度的光滑处理。最后得到有一定光滑度的任意起伏地形模型(图2)。此外,对震源处附近的网格进行细化处理,保证计算的精度。然后根据传播网格周围的速度网格,通过插值方法得到传播网格节点的速度值。为了保证计算的正确,有必要对网格进行检查,并进行适当调整。
步骤三:对每一个三角形网格进行计算。
首先,从震源处所在的网格,开始计算。从震源开始,波前传至周围的网格点。刚开始,震源处的网格点是上风(Upwind)集合,Upwind集合是已经计算好的网格集合,这些正在计算得周围网格点的集合称为窄带(Narrow band),未计算的网格点集合是下风(Downwind)。如此,从震源处开始计算,一直到曲面上所有的网格点都遍历,都成为Upwind集合,如图2所示。
下面主要介绍每一个三角形网格点的计算算法。
传统的FMM算法都是在正交网格中进行计算程函方程。
式中,T是波前的走时,s(x)是每一个网格点的慢度值。在正交坐标系中可以求解。但是,对于任意曲面剖分成的任意三角型网格,并不是正交网格,因此,可以应用Sethian1999中的办法解决这一问题。
其中各参数见图3。不同于四边形网格,三角形是任意划分的,无法使用有限差分的高阶算法,因此这种方法只有一阶近似,精度不是很高。为了提高计算精度,我们的发明改进了这一算法。具体方法阐述如下:
如图4所示,已经三角形网格点i和j处的波前到时Ti和Tj,要求To的到时。传统方法中,假设了在三角形走时梯度相同,使用了O点的慢度。因此精度较差。为了增加计算精度,在计算Ti和Tj时,可以利用i和j节点处的走时梯度大小(慢度)和在方向导数 如果采取这种算法,需要在计算每个节点走时的同时,也要保存每个节点的走时梯度方向。下面的公式(2)和(3)表达了网格点i、j处a和b方向导数和三角形两条边a和b上的平均方向导数之间的关系。
其中分别表示网格点i、j和o处方向和方向的方向导数,是三角形a和b上的平均方向导数。利用这两个公式可以用已知网格点i和j上的方向导数和三角形边上的平均方向导数得到网格点O处方向的方向导数。
然后,我们就可以建立三角形中FMM算法。
图4显示了如何由网格点O处方向导数,求O处梯度的算法。建立如图4的直角坐标系,所求的O处走时梯度坐标为(x0,y0),那么建立下列方程组:
其中分别是网格点O处方向的方向导数,而且是网格点O到时To的函数(由公式(2)和(3)求得),SO是网格点O处的慢度。将方程组(4)求得的x0和y0代入公式(5)。最后得到:
由此公式求得网格点O处的到时TO。这样就完成了一个三角形网格中的计算。如此,循环下去,一直到所有网格点都变成上风(Upwind)集合。至此所有网格点的到时都已经计算完成。
步骤四:找出震源至台站的射线路径。从台站开始,在每个网格中求梯度,负梯度的方向就是射线方向,把射线沿途经过的网格点中负梯度都求出来。从接受台站开始,连接所有沿途网格中射线的线段,直至震源处,就得到了震源-台站的射线路径。此外,在计算每个网格点的到时的同时,同时计算每个速度网格节点速度的微小变化对走时的影响,也就是进行走时计算的同时,也计算出Frechet导数值,为面波反演作准备。保存走时信息和Frechet导数值。至此,FMM方法的射线追踪就完成了。
以上仅是本发明技术要点,凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理前提下的若干改进和润饰,应视为本发明的保护范围。

Claims (3)

1.一种在任意曲面实现射线追踪的任意曲面快速行进方法,其特征在于所述任意曲面快速行进方法的核心技术是在非正交网格实现走时的计算;
所述方法包括以下步骤:
(1)建立速度网格模型(Velocity grids),对曲面进行网格划分,给出每个网格点的速度值,划分得到的网格为三角形网格;
(2)建立起伏地形模型,进行传播网格化剖分,并对地形起伏进行光滑处理以保证算法的效率;然后得到有光滑度的任意起伏地形模型;对震源处附近的网格进行细化处理以保证计算的精度;然后根据传播网格周围的速度网格,通过插值方法得到传播网格节点的速度值;最后对网格进行检查,并进行适当调整以保证计算的正确;
(3)对每一个三角形网格进行计算;首先,从震源处所在的网格开始计算;即从震源开始,波前传至周围的网格点;从震源处开始计算,一直到曲面上所有的网格点都遍历,完成射线追踪;
(4)找出震源至台站的射线路径;从台站开始,在每个网格中求梯度,负梯度的方向就是射线方向,把射线沿途经过的网格点中负梯度都求出来;从接受台站开始,连接所有沿途网格中射线的线段,直至震源处,得到震源-台站的射线路径;此外,在计算每个网格点的到时的同时,计算每个速度网格节点速度的微小变化对走时的影响,也就是进行走时计算的同时,也计算出Frechet导数值为面波反演作准备;保存走时信息和Frechet导数值,完成FMM方法的射线追踪。
2.如权利要求1所述的任意曲面快速行进方法,其特征在于建立速度网格模型步骤中,将曲面划分三角形网格,并在三角形网格中进行射线追踪。
3.如权利要求2所述的任意曲面快速行进方法,其特征在于划分三角形网格即为建立锐角三角形,具体分为在不考虑地表起伏的水平面建立等边三角形;而若考虑到地形的起伏,所建立的三角形内角在60度附近。
CN201711354734.XA 2017-12-15 2017-12-15 任意曲面非均匀介质快速行进程函方程求解射线追踪算法 Active CN108267781B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711354734.XA CN108267781B (zh) 2017-12-15 2017-12-15 任意曲面非均匀介质快速行进程函方程求解射线追踪算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711354734.XA CN108267781B (zh) 2017-12-15 2017-12-15 任意曲面非均匀介质快速行进程函方程求解射线追踪算法

Publications (2)

Publication Number Publication Date
CN108267781A true CN108267781A (zh) 2018-07-10
CN108267781B CN108267781B (zh) 2020-06-05

Family

ID=62772217

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711354734.XA Active CN108267781B (zh) 2017-12-15 2017-12-15 任意曲面非均匀介质快速行进程函方程求解射线追踪算法

Country Status (1)

Country Link
CN (1) CN108267781B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110568496A (zh) * 2019-09-26 2019-12-13 核工业北京地质研究院 一种复杂介质条件下射线追踪方法
CN113221395A (zh) * 2021-03-16 2021-08-06 禁核试北京国家数据中心 基于分层介质的地震走时参数化网格模型构建方法及应用
CN114429047A (zh) * 2022-01-27 2022-05-03 成都理工大学 一种基于三角网格的二次圆方程旅行时插值方法
CN117607957A (zh) * 2024-01-24 2024-02-27 南方科技大学 基于等效慢度的快速推进法的地震波走时求解方法及系统
CN118011483A (zh) * 2024-02-04 2024-05-10 湖南工商大学 一种考虑地形效应的面波层析成像方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7315783B2 (en) * 2006-01-06 2008-01-01 Baker Hughes Incorporated Traveltime calculation in three dimensional transversely isotropic (3D TI) media by the fast marching method
CN103698810A (zh) * 2013-12-10 2014-04-02 山东科技大学 混合网最小走时射线追踪层析成像方法
CN105842735A (zh) * 2016-05-20 2016-08-10 四川大学 具有复杂速度分布的区域岩体微震震源定位方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7315783B2 (en) * 2006-01-06 2008-01-01 Baker Hughes Incorporated Traveltime calculation in three dimensional transversely isotropic (3D TI) media by the fast marching method
CN103698810A (zh) * 2013-12-10 2014-04-02 山东科技大学 混合网最小走时射线追踪层析成像方法
CN105842735A (zh) * 2016-05-20 2016-08-10 四川大学 具有复杂速度分布的区域岩体微震震源定位方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
CHENG CHENG ET AL.: "Three-dimensional pre-stack depth migration of receiver functions with the fast marching method: a Kirchhoff approach", 《GEOPHYS. J. INT.》 *
N. RAWLINSON AND M. SAMBRIDGE: "Multiple reflection and transmission phases in complex layered media using a multistage fast marching method", 《GEOPHYSICS》 *
PETER G. LELIEVRE ET AL.: "Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the FastMarchingMethod", 《GEOPHYS. J. INT.》 *
唐秋惠,等: "三维多步FMM射线追踪在正态分布速度地层模型中的应用", 《工程地球物理学报》 *
孟宪海,等: "基于三角域快速行进法的地震波走时计算", 《软件》 *
赵瑞,等: "复杂层状模型中多次波快速追踪-一种基于非规则网格的最短路径算法", 《地震学报》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110568496A (zh) * 2019-09-26 2019-12-13 核工业北京地质研究院 一种复杂介质条件下射线追踪方法
CN113221395A (zh) * 2021-03-16 2021-08-06 禁核试北京国家数据中心 基于分层介质的地震走时参数化网格模型构建方法及应用
CN113221395B (zh) * 2021-03-16 2023-10-20 禁核试北京国家数据中心 基于分层介质的地震走时参数化网格模型构建方法及应用
CN114429047A (zh) * 2022-01-27 2022-05-03 成都理工大学 一种基于三角网格的二次圆方程旅行时插值方法
CN114429047B (zh) * 2022-01-27 2023-08-22 成都理工大学 一种基于三角网格的二次圆方程旅行时插值方法
CN117607957A (zh) * 2024-01-24 2024-02-27 南方科技大学 基于等效慢度的快速推进法的地震波走时求解方法及系统
CN117607957B (zh) * 2024-01-24 2024-04-02 南方科技大学 基于等效慢度的快速推进法的地震波走时求解方法及系统
CN118011483A (zh) * 2024-02-04 2024-05-10 湖南工商大学 一种考虑地形效应的面波层析成像方法及系统

Also Published As

Publication number Publication date
CN108267781B (zh) 2020-06-05

Similar Documents

Publication Publication Date Title
CN108267781A (zh) 任意曲面非均匀介质快速行进程函方程求解射线追踪算法
CN108064348B (zh) 一种基于两点射线追踪的地震走时层析反演方法
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
CN102053258B (zh) 基于复杂地质构造的自适应三维射线追踪方法
CN110058307B (zh) 一种基于快速拟牛顿法的全波形反演方法
CN107843922A (zh) 一种基于地震初至波和反射波走时联合的层析成像方法
CN106814391A (zh) 基于菲涅尔体层析反演的地面微地震事件定位方法
WO2019242045A1 (zh) 一种虚拟震源二维波前构建地震波走时计算方法
US10162070B2 (en) Converting a first acquired data subset to a second acquired data subset
CN109444955A (zh) 三维地震射线追踪的双线性走时扰动插值方法
CN116774292B (zh) 一种地震波走时确定方法、系统、电子设备及存储介质
CN108303736B (zh) 各向异性ti介质最短路径射线追踪正演方法
CN105573963B (zh) 一种电离层水平不均匀结构重构方法
AU2015394432A1 (en) Model compression
Yu et al. A minimum traveltime ray tracing global algorithm on a triangular net for propagating plane waves
CN106353801B (zh) 三维Laplace域声波方程数值模拟方法及装置
CN109557588B (zh) 一种煤矿井下二维矿震波速反演降维方法
CN108983290B (zh) 一种三维横向各向同性介质中旅行时确定方法及系统
Yedlin et al. Uniform asymptotic conversion of Helmholtz data from 3D to 2D
Zhang et al. A new 3-D ray tracing method based on LTI using successive partitioning of cell interfaces and traveltime gradients
CN102998702A (zh) 保幅平面波叠前深度偏移方法
CN109655888B (zh) 地震数据处理中光滑浮动基准面的定量选择方法及系统
Li et al. Application of Snell's law in reflection raytracing using the multistage fast marching method
Martinelli An application of the level set method to underwater acoustic propagation
Padhi et al. Accurate quasi-p traveltimes in 3d transversely isotropic media using a high-order fast-sweeping-based eikonal solver

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