CN112596103A - 射线追踪方法、装置和电子设备 - Google Patents

射线追踪方法、装置和电子设备 Download PDF

Info

Publication number
CN112596103A
CN112596103A CN202011334054.3A CN202011334054A CN112596103A CN 112596103 A CN112596103 A CN 112596103A CN 202011334054 A CN202011334054 A CN 202011334054A CN 112596103 A CN112596103 A CN 112596103A
Authority
CN
China
Prior art keywords
point
travel time
target
wave
points
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.)
Pending
Application number
CN202011334054.3A
Other languages
English (en)
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.)
Institute of Geophysical and Geochemical Exploration of CAGS
Original Assignee
Institute of Geophysical and Geochemical Exploration of CAGS
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 Institute of Geophysical and Geochemical Exploration of CAGS filed Critical Institute of Geophysical and Geochemical Exploration of CAGS
Priority to CN202011334054.3A priority Critical patent/CN112596103A/zh
Publication of CN112596103A publication Critical patent/CN112596103A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • 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
    • G01V1/282—Application of seismic models, synthetic seismograms
    • 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
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling

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)

Abstract

本发明提供了一种射线追踪方法、装置和电子设备，涉及地震勘探的技术领域，包括获取网格离散化的速度模型和射线追踪观测系统的系统参数；基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时；基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径，其中，射线路径中目标震源点的位置包括以下任一种：网格节点、网格边界、网格内部，目标震源点为射线路径中所有次级震源点中的任意一个。本发明方法在进行射线追踪时考虑了次级震源点出现在网格内部和网格边界的情况，有效的提高了射线路径的计算精度，从而缓解了现有技术中的射线追踪方法存在的计算精度低的技术问题。

Description

射线追踪方法、装置和电子设备
技术领域
本发明涉及地震勘探的技术领域，尤其是涉及一种射线追踪方法、装置和电子设备。
背景技术
射线追踪广泛应用于地震勘探的各个领域，为叠前偏移、速度分析和层析成像等提供了有效的手段，现有的射线追踪方法，在求解全局离散旅行时场的同时计算射线方向，这就需要计算网格节点的射线矢量或所有可能的射线路径，这种射线追踪方法一方面增加了存储空间，另一方面次级震源点也只能位于网格节点，导致射线路径的计算精度较差。
综上所述，现有技术中的射线追踪方法存在计算精度低的技术问题。
发明内容
本发明的目的在于提供一种射线追踪方法、装置和电子设备，以缓解了现有技术中的射线追踪方法存在的计算精度低的技术问题。
第一方面，本发明提供一种射线追踪方法，包括：获取网格离散化的速度模型和射线追踪观测系统的系统参数；其中，所述速度模型为射线追踪区域，所述系统参数包括震源点的空间位置和接收点的空间位置；基于所述网格离散化的速度模型和所述系统参数计算所述射线追踪区域内所有离散点的目标旅行时；其中，所述目标旅行时包括：初至波旅行时和后续波旅行时，所述后续波包括：一次反射波和多次反射波；基于所述所有离散点的目标旅行时和所述系统参数反向计算所述接收点到所述震源点的射线路径；其中，所述射线路径中目标震源点的位置包括以下任一种：网格节点、网格边界、网格内部，所述目标震源点为所述射线路径中所有次级震源点中的任意一个。
在可选的实施方式中，所述速度模型的模型参数包括：所属于透反射界面的若干离散点的空间位置；所述震源点的初至波旅行时为0；基于所述网格离散化的速度模型和所述系统参数计算所述射线追踪区域内所有离散点的目标旅行时，包括：利用窄带扩展方法依次计算所述网格离散化的速度模型内所有网格节点的初至波旅行时；基于所述所属于透反射界面的若干离散点的空间位置和所述所有网格节点的初至波旅行时确定界面初至旅行时；其中，所述界面初至旅行时包括所述透反射界面上所有离散点的初至波旅行时；基于所述界面初至旅行时和预设计算分区计算所述后续波旅行时；其中，所述预设计算分区为所述速度模型根据所述后续波的类型划分的计算区域；基于所述初至波旅行时和所述后续波旅行时确定所述目标旅行时。
在可选的实施方式中，利用窄带扩展方法依次计算所述网格离散化的速度模型内所有网格节点的初至波旅行时，包括：将所述网格离散化的速度模型中的所有网格节点进行分类，得到分类后的网格节点和窄带；其中，所述分类后的网格节点的节点类型包括：近点，窄带点和远点，所述震源点为近点，在未进行窄带扩展时所述窄带点仅包括震源点的上下左右四个邻点，所述远点为除所述近点和所述窄带点之外的网格节点；所述窄带为所有窄带点的集合；基于所述震源点的初至波旅行时计算所有窄带点的初至波旅行时；重复执行以下步骤，直至所有网格节点的节点类型均为近点：基于所有窄带点的初至波旅行时确定所述窄带中的最小旅行时点，并将所述最小旅行时点更新为近点；判断所述最小旅行时点的目标邻点的节点类型；其中，所述目标邻点为所述最小旅行时点的上下左右四个邻点中的任意一个；若所述目标邻点为近点，则所述目标邻点的初至波旅行时保持不变；若所述目标邻点为窄带点或远点，则基于所述最小旅行时点的初至波旅行时更新所述目标邻点的初至波旅行时，且将节点类型为远点的目标邻点移入所述窄带。
在可选的实施方式中，基于所述所属于透反射界面的若干离散点的空间位置和所述所有网格节点的初至波旅行时确定界面初至旅行时，包括：基于所述所属于透反射界面的若干离散点的空间位置和插值法离散化所述透反射界面，得到所述透反射界面上所有离散点的空间位置；利用算式tp＝tm+(tn-tm)L/|MN|计算目标离散点的初至波旅行时；其中，所述目标离散点为所述透反射界面上所有离散点中的任意一个，tp表示所述目标离散点的初至波旅行时，tm和tn分别表示所述目标离散点所在的网格边的两个端点m，n的初至波旅行时，且tn>tm，L表示所述目标离散点与端点m之间的距离，0≤L≤|MN|，|MN|表示端点m和端点n之间的距离；基于所述透反射界面上所有离散点的初至波旅行时确定所述界面初至旅行时。
在可选的实施方式中，基于所述界面初至旅行时和预设计算分区计算所述后续波旅行时，包括：判断后续波是否为初次经过目标透反射界面；其中，所述目标透反射界面为所有透反射界面中的任意一个界面；若是，则基于所述目标透反射界面的界面初至旅行时计算目标计算区域内所有网格节点的后续波旅行时；其中，所述目标计算区域为所述预设计算分区中所述后续波的下一传输方向对应的计算区域；若否，则基于所述目标透反射界面的目标界面旅行时计算所述目标计算区域内所有网格节点的后续波旅行时；其中，所述目标界面旅行时包括利用所述后续波上一计算区域内的所有网格节点的后续波旅行时插值得到的所述目标透反射界面上所有离散点的后续波旅行时。
在可选的实施方式中，基于所述所有离散点的目标旅行时和所述系统参数反向计算所述接收点到所述震源点的射线路径，包括：追踪目标波的射线路径时，针对所述射线路径中的目标射线，利用算式
Figure BDA0002796814550000041
计算所述目标射线的终点旅行时，以及，利用算式
Figure BDA0002796814550000042
计算所述目标射线的起点空间位置；其中，若已知所述接收点的目标波旅行时，则不使用所述终点旅行时的算式对所述接收点旅行时进行更新，所述目标波包括以下任一种：初至波，一次反射波，多次反射波，所述目标射线为所述射线路径中的任意一段射线，tc表示所述目标射线的终点C的目标波旅行时，ta表示所述目标射线所在的网格中已知离散点A的目标波旅行时，tb表示所述目标射线所在的网格中已知离散点B的目标波旅行时，
Figure BDA0002796814550000043
表示线段AB与线段AC之间的夹角，v表示波速，|AC|表示所述已知离散点A和所述目标射线的终点C之间的距离，|AB|表示所述已知离散点A和所述已知离散点B之间距离，(xa,za)表示所述已知离散点A的空间位置，(xb,zb)表示所述已知离散点B的空间位置，(xd,zd)表示所述目标射线的起点空间位置，r表示所述已知离散点A和所述目标射线的起点之间的距离，
Figure BDA0002796814550000044
按照接收点至震源点的方向将所有目标射线依次进行连接，得到所述目标波的射线路径。
在可选的实施方式中，所述网格离散化的速度模型中震源点附近的网格采用双重网格技术进行加密，以使所述速度模型划分为：加密网格区域和非加密网格区域。
第二方面，本发明提供一种射线追踪装置，包括：获取模块，用于获取网格离散化的速度模型和射线追踪观测系统的系统参数；其中，所述速度模型为射线追踪区域，所述系统参数包括震源点的空间位置和接收点的空间位置；第一计算模块，用于基于所述网格离散化的速度模型和所述系统参数计算所述射线追踪区域内所有离散点的目标旅行时；其中，所述目标旅行时包括：初至波旅行时和后续波旅行时，所述后续波包括：一次反射波和多次反射波；第二计算模块，用于基于所述所有离散点的目标旅行时和所述系统参数反向计算所述接收点到所述震源点的射线路径；其中，所述射线路径中目标震源点的位置包括以下任一种：网格边界、网格节点、网格内部，所述目标震源点为所述射线路径中所有次级震源点中的任意一个。
第三方面，本发明提供一种电子设备，包括存储器、处理器，所述存储器上存储有可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述前述实施方式中任一项所述的方法的步骤。
第四方面，本发明提供一种具有处理器可执行的非易失的程序代码的计算机可读介质，所述程序代码使所述处理器执行前述实施方式中任一项所述的方法。
本发明提供射线追踪方法，包括：获取网格离散化的速度模型和射线追踪观测系统的系统参数；其中，速度模型为射线追踪区域，系统参数包括震源点的空间位置和接收点的空间位置；基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时；其中，目标旅行时包括：初至波旅行时和后续波旅行时，后续波包括：一次反射波和多次反射波；基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径；其中，射线路径中目标震源点的位置包括以下任一种：网格节点、网格边界、网格内部，目标震源点为射线路径中所有次级震源点中的任意一个。
现有技术中的射线追踪方法在进行射线追踪时仅考虑所有次级震源点位于网格节点的情况，从而难以保证射线路径的计算精度。与现有技术相比，本发明提供的射线追踪方法，首先获取网格离散化的速度模型和射线追踪观测系统的系统参数；然后基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时；最后基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径。本发明方法在进行射线追踪时考虑了次级震源点出现在网格内部和网格边界的情况，有效的提高了射线路径的计算精度，从而缓解了现有技术中的射线追踪方法存在的计算精度低的技术问题。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种射线追踪方法的流程图；
图2为本发明实施例提供的一种标记有透反射界面和速度分区的速度模型示意图；
图3为本发明实施例提供的一种一次反射波和三次反射波的分区示意图；
图4为本发明实施例提供的一种波前窄带扩展示意图；
图5为本发明实施例提供的一种完全二叉树与数组的对应关系的示意图；
图6为本发明实施例提供的一种窄带扩展示意图；
图7为本发明实施例提供的一种传统Vidale差分和改进型Vidale差分网格节点示意图；
图8为本发明实施例提供的一种网格节点位置关系及模型区域内坐标方向示意图；
图9为本发明实施例提供的一种分区多步射线追踪的示意图；
图10为本发明实施例提供的一种由A和B两点的旅行时求C点旅行时及D点坐标的示意图；
图11为本发明实施例提供的一种由接收点向震源检索的射线路径追踪过程示意图；
图12为本发明实施例提供的一种双重网格技术的示意图；
图13为本发明实施例提供的一种网格间距为1m时不同差分格式下的射线路径与真实值对比；
图14为本发明实施例提供的一种不同网格间距不同差分格式计算误差的对比示意图；
图15为本发明实施例提供的一种高速夹层的模型参数、观测系统及射线追踪结果的示意图；
图16(a)为本发明实施例提供的一种使用本发明方法进行射线追踪的结果示意图；
图16(b)为本发明实施例提供的一种使用LTI技术进行射线追踪的结果示意图；
图16(c)为本发明实施例提供的一种本发明方法与LTI技术的旅行时对比示意图；
图17为本发明实施例提供的一种Marmousi模型参数及反射追踪模型分区示意图；
图18为本发明实施例提供的一种一次反射波追踪示意图；
图19为本发明实施例提供的一种多次反射波追踪示意图；
图20为本发明实施例提供的一种Marmousi模型单炮合成地震记录的示意图；
图21为本发明实施例提供的一种速度模型及反射界面的示意图；
图22(a)为本发明实施例提供的一种地表观测系统射线路径追踪结果示意图；
图22(b)为本发明实施例提供的一种地井观测系统射线路径追踪结果示意图；
图22(c)为本发明实施例提供的一种井地观测系统射线路径追踪结果示意图；
图22(d)为本发明实施例提供的一种井间观测系统射线路径追踪结果示意图；
图23为本发明实施例提供的一种不同观测系统下的三炮合成记录的示意图；
图24为本发明实施例提供的一种射线追踪装置的功能模块图；
图25为本发明实施例提供的一种电子设备的示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图，对本发明的一些实施方式作详细说明。在不冲突的情况下，下述的实施例及实施例中的特征可以相互组合。
实施例一
图1为本发明实施例提供的一种射线追踪方法的流程图，如图1所示，该方法具体包括如下步骤：
步骤S102，获取网格离散化的速度模型和射线追踪观测系统的系统参数。
具体的，在进行射线追踪时，首先要获取网格离散化的速度模型和射线追踪观测系统的系统参数，其中，速度模型为射线追踪区域，系统参数包括震源点的空间位置和接收点的空间位置。在本发明实施例中，射线追踪区域为网格离散化的速度模型，速度模型可以采用现有技术中已经离散化的速度模型(如公知的Marmousi模型、盐丘模型等)，或者在笛卡尔坐标系中，构建模型区域并采用矩形网格离散化模型区域，然后根据离散化后网格节点的空间位置给定速度来构建整个速度模型。震源点的空间位置和接收点的空间位置由用户给定，且接收点的数量不限，可以为一个或多个，用户可以根据实际需求进行设定，并利用震源点的空间位置和接收点的空间位置构建观测系统。
步骤S104，基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时。
在确定了射线追踪区域，震源点的空间位置以及接收点的空间位置后，基于以上数据即可计算出该射线追踪区域内所有离散点的目标旅行时，其中，目标旅行时包括：初至波旅行时和后续波旅行时，后续波包括：一次反射波和多次反射波。在本发明实施例中，离散点并不限定为矩形网格的节点(四个顶点)，如果透反射界面与网格的交点不是网格节点，则上述交点的目标旅行时也能计算出来。
步骤S106，基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径。
其中，射线路径中目标震源点的位置包括以下任一种：网格节点、网格边界、网格内部，目标震源点为射线路径中所有次级震源点中的任意一个。
在得到所有离散节点的目标旅行时之后，结合震源点的空间位置，接收点的空间位置，利用插值法对射线路径进行反向追踪，得到接收点到震源点的射线路径，且在进行反向追踪过程中，考虑次级震源点出现在网格边界和网格内部的情况，以使得到的射线路径更为精确。反向追踪的含义是，假设射线路径为A→B→C→D→E，那么在追踪上述射线路径时，首先要确定出D→E的这一段，然后依次是C→D，B→C，最后才是A→B。
现有技术中的射线追踪方法在进行射线追踪时仅考虑所有次级震源点位于网格节点的情况，从而难以保证射线路径的计算精度。与现有技术相比，本发明提供的射线追踪方法，首先获取网格离散化的速度模型和射线追踪观测系统的系统参数；然后基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时；最后基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径。本发明方法在进行射线追踪时考虑了次级震源点出现在网格内部和网格边界的情况，有效的提高了射线路径的计算精度，从而缓解了现有技术中的射线追踪方法存在的计算精度低的技术问题。
在一个可选的实施方式中，速度模型的模型参数包括：所属于透反射界面的若干离散点的空间位置；震源点的初至波旅行时为0；上述步骤S104，基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时，具体包括如下步骤：
步骤S1041，利用窄带扩展方法依次计算网格离散化的速度模型内所有网格节点的初至波旅行时。
具体的，在本发明实施例中，在震源点的初至波旅行时为0的情况下，为了求解网格离散化的速度模型内所有网格节点的初至波旅行时，具体使用窄带扩展法，这一过程可以想象为火焰在纸上的燃烧过程，一旦燃烧过的位置则认为是保持恒定的状态，不能再次被燃烧，转化到波前扩展过程则可将燃烧过的位置视为波前扩展过的区域，火焰前进的条带视为波前，未燃烧的位置视为待求旅行时的网格节点。根据波前传播过程中的单向性，地震波的波前由旅行时较小的网格节点扩展到旅行时较大的网格节点。采用窄带技术进行波前扩展，将波前面上的极小旅行时点作为波前扩展点，通过速度模型区域网格节点的属性来冻结“完成燃烧”的节点，并利用窄带的扩展更新“正在燃烧”的节点，由此计算所有网格节点的旅行时。
步骤S1042，基于所属于透反射界面的若干离散点的空间位置和所有网格节点的初至波旅行时确定界面初至旅行时。
在计算得到速度模型内所有网格节点的初至波旅行时之后，还需要进一步计算透反射界面上所有离散点的初至波旅行时，进而利用界面初至旅行时才能计算后续波旅行时，图2为本发明实施例提供的一种标记有透反射界面和速度分区的速度模型示意图，在获取到网格离散化的速度模型时，能够同时获得速度模型的模型参数，包括：所属于透反射界面的若干离散点的空间位置，一般地，透反射界面在给定时，往往是提供若干离散点(透反射界面上的部分离散点(控制点)，并非透反射界面上的全部离散点)的空间位置，为了得到透反射界面上所有离散点的初至波旅行时，首先要根据上述若干离散点的空间位置先确定透反射界面上所有离散点的空间位置，然后再进一步通过插值计算每个离散点的初至波旅行时，最后得到界面初至旅行时，其中，界面初至旅行时包括透反射界面上所有离散点的初至波旅行时。由于界面初至旅行时是由插值得到的，存在一定误差，因此，在实际操作中，可先由震源计算整个模型区域及各界面的初至旅行时，避免利用界面初至旅行时计算后续初至旅行时场而产生累计误差。
步骤S1043，基于界面初至旅行时和预设计算分区计算后续波旅行时。
通过上文中的描述可知，后续波包括一次反射波和多次反射波，后续波在传输过程中，其旅行时在速度模型中的不同速度分区需要分别计算，且在计算后续波旅行时的时候，需要根据界面初至旅行时和预设计算分区进行计算，上述预设计算分区为速度模型根据后续波的类型划分的计算区域，用户可以根据实际需求设置后续波的传播方向，也即，预设计算分区是用户预先给定的。图3为本发明实施例提供的一种一次反射波和三次反射波的分区示意图。
图3中(a)视图为一次反射模型分区示意图，图中箭头表示地震波在第二条透反射界面处发生一次反射，其下行波和上行波在穿过透反射界面时可根据需要选择是否发生透射和反射转换。按照分区多步思想，该反射波可分为四段进行计算，则相应的模型可分为四个计算区域，并在每个区域依次计算相应的下行波或上行波的旅行时。该反射波计算的具体实现过程分为四步：(1)计算震源到第一条透反射界面的下行波波前；(2)计算第一条透反射界面到第二条透反射界面的下行波波前；(3)计算第二条透反射界面到第一条透反射界面的上行波波前；(4)计算第一条透反射界面到地表的上行波波前。在各步计算过程中根据需要选用不同类型的速度，可实现第二条界面任意类型的P波(纵波)和S波(横波)组合得到的一次反射。若(1)、(2)中所用速度类型一致，则(1)、(2)可合并为计算震源到第二条透反射界面的下行波；若(3)、(4)中所用速度类型一致，则(3)、(4)可合并为计算第二条透反射界面到地表的上行波。
图3中(b)视图为多次反射波的模型分区示意图，图中箭头表示地震波在第二层内发生三次反射，其下行波和上行波在穿过透反射界面时可根据需要选择是否发生透射和反射转换。按照分区多步思想，该反射波可分为六步进行计算，则相应的模型可分为六个计算区域，并在每个区域依次计算相应的下行波或上行波的旅行时。该多次反射波计算的具体实现过程分为六步：(1)计算震源到第一条透反射界面的下行波波前；(2)计算第一条透反射界面到第二条透反射界面的下行波波前；(3)计算第二条透反射界面到第一条透反射界面的上行波波前；(4)计算第一条透反射界面到第二条透反射界面的下行波波前；(5)计算第二条透反射界面到第一条透反射界面的上行波波前；(6)计算第一条透反射界面到地表的上行波波前。在各步计算过程中根据需要选用不同类型的速度，可实现第二层内任意类型的P波和S波组合得到的多次反射。若(1)、(2)中所用速度类型一致，则(1)、(2)可合并为计算震源到第二条透反射界面的下行波；若(5)、(6)中所用速度类型一致，则(5)、(6)可合并为计算第二条透反射界面到地表的上行波。
一次反射波分区多步计算时，设模型含n条界面(不含模型上下边界)，即模型含n+1层，界面或层序号从上到下依次为1，2…n+1，则第i条界面的一次反射波的最多模型分区为2i，那么，按照分区多步可追踪的该界面一次反射波的类型为
Figure BDA0002796814550000131
(W表示P波或S波，下标表示界面或层号，上标1表示下行波，-1表示上行波)，若其中相邻两个分区内的波类型及传播方向一致，则这两个分区可合并为一个计算区。需要注意的是在发生透射或反射时P波能够转换为S波，但S波不可转换为P波。例如，第一条透反射界面可追踪的一次反射波的类型有3种：
Figure BDA0002796814550000132
而
Figure BDA0002796814550000133
是不存在的，其余透反射界面的一次反射波以此类推。
多次反射波分区多步计算时，图3中(b)视图的多次反射波其模型分区为6个，那么按照分区多步可追踪的该多次波的类型为
Figure BDA0002796814550000134
若其中相邻两个分区内的波类型及传播方向一致，则这两个分区可合并为一个计算区，其他类型多次波的追踪以此类推。需要注意的是在发生透射或反射时P波能够转换为S波，但S波不可转换为P波，
Figure BDA0002796814550000135
是不存在的。
步骤S1044，基于初至波旅行时和后续波旅行时确定目标旅行时。
上文中对计算目标旅行时的方法流程进行了概述，下面对其中涉及的各个步骤进行具体介绍。
在一个可选的实施方式中，上述步骤S1041，利用窄带扩展方法依次计算网格离散化的速度模型内所有网格节点的初至波旅行时，具体包括如下步骤：
步骤S10411，将网格离散化的速度模型中的所有网格节点进行分类，得到分类后的网格节点和窄带。
具体的，在计算网格离散化的速度模型内所有网格节点的初至波旅行时的时候，首先，要将所有网格节点进行分类，其中，分类后的网格节点的节点类型包括：近点，窄带点和远点，震源点为近点(已完成旅行时计算的点)，在未进行窄带扩展时窄带点仅包括震源点的上下左右四个邻点，远点为除近点和窄带点之外的网格节点；窄带为所有窄带点的集合。
步骤S10412，基于震源点的初至波旅行时计算所有窄带点的初至波旅行时。
在将网格节点分类完毕后，首先计算当前状态下窄带(初始窄带)中所有窄带点的初至波旅行时，如果震源点处于模型边界处，则初始窄带中不足4个窄带点。设震源点(is,js)的旅行时为t(is,js)＝0，那么震源点上下左右四个邻点的旅行时为，
Figure BDA0002796814550000141
其中，s表示地震波的慢度，即速度的倒数，h表示网格间距，初始状态下，远点的旅行时可以视为无穷大。
在进行窄带扩展过程中，近点逐渐增多，远点逐渐减少，图4为本发明实施例提供的一种波前窄带扩展示意图，当速度模型内所有的网格节点均为近点时，窄带扩展结束。所以为求得模型内所有网格节点的初至波旅行时，需要重复执行以下步骤，即重复执行步骤S10413至步骤S10416，直至所有：
步骤S10413，基于所有窄带点的初至波旅行时确定窄带中的最小旅行时点，并将最小旅行时点更新为近点。
窄带扩展的过程是近点逐渐增多的过程，新增近点为所有窄带点中旅行时最小的网格节点，也即，局部极小旅行时点，本发明实施例采用堆排序技术寻找窄带中最小旅行时t(imin,jmin)的网格节点(imin,jmin)(极小旅行时点)，用户也可以利用其他的技术手段来确定极小旅行时点。
步骤S10414，判断最小旅行时点的目标邻点的节点类型。
其中，目标邻点为最小旅行时点的上下左右四个邻点中的任意一个。
若目标邻点为近点，则执行步骤S10415；若目标邻点为窄带点或远点，则执行步骤S10416。
步骤S10415，目标邻点的初至波旅行时保持不变。
步骤S10416，基于最小旅行时点的初至波旅行时更新目标邻点的初至波旅行时，且将节点类型为远点的目标邻点移入窄带。
堆是一个完全二叉树的数据结构，即由n个元素构成的序列{a1，a2，…，an}的元素满足ai≤a2i，ai≤a2i+1或ai≥a2i，ai≥a2i+1，(i＝1，2，3…，n/2)，则该序列就称为堆，其最顶端的结点称为根结点，而每个结点作为父结点都有两个子结点。窄带扩展中使用的堆为小根堆，也称为最小堆，即在二叉树中每个父结点的关键值均小于其两个子结点的关键值。完全二叉树的数据结构和数组的对应关系如图5所示。
堆数据结构和堆排序技术的实现体现在窄带扩展过程中，其实现过程为：首先，初始化数据结构形成最小堆，即计算震源点周围网格节点的旅行时形成初始窄带，将窄带内网格节点信息(包括旅行时及节点坐标)按最小堆结构存贮，最小旅行时点位于堆顶；其次，提取根结点的信息，即将窄带内极小旅行时点调整为近点，也就是删除堆顶元素并重建堆数据结构。最后，加入新的元素或更新原有元素的关键值，即对上步中极小旅行时点的邻点属性进行判断，若其为近点则保持不变，若其为远点则计算旅行时并放入窄带，也就是在堆中加入新元素；若其为窄带点则更新其旅行时，也就是更新原有元素的关键值，并对更新后的堆重建。
图6为本发明实施例提供的一种窄带扩展示意图，通过上文中的描述可知，在得到窄带中的最小旅行时点后，将其节点类型更新为近点，接下来，要判断该新增近点上下左右四个临点的节点类型，针对近点类型的邻点，则保持其初至波旅行时不变，但是如果临点不是近点，则需要更新其旅行时，且针对远点类型的邻点，在更新完其旅行时后，还需要将其节点类型更新为窄带点，也即，将其移入窄带。
在更新窄带点和远点的初至波旅行时的时候，本发明实施例采用Vidale差分公式来进行旅行时的计算，Vidale差分方法含两种传统差分及一种改进型差分，图7中左侧和中间为两种传统Vidale差分的网格节点示意图，右侧为改进型Vidale差分的网格节点示意图，其中，tC为待求点旅行时，tA、tB、tD为已知网格旅行时，且tA为相对极小旅行时，网格间距h。
图7左侧的传统差分方法中，其差分格式表示为：
Figure BDA0002796814550000161
为求解该方程，将其带入二维情况下的程函方程
Figure BDA0002796814550000162
其中，t表示地震波的旅行时，s表示地震波的慢度，即速度的倒数。在笛卡尔坐标系中x表示水平方向，代表水平方向的位置；z表示垂直向下方向，代表深度方向的位置。旅行时函数t(x,z)是空间坐标(x,z)的函数，其在x、z方向上一阶偏微分算子分别为
Figure BDA0002796814550000163
在数值计算时将偏微分差分近似代替，用矩形网格将模型区域网格化，网格节点位置关系及模型区域内坐标方向示意图如图8所示。解方程，得到待求点C旅行时公式为：
Figure BDA0002796814550000171
其中，
Figure BDA0002796814550000172
sA、sB、sC和sD分别表示点A、B、C、和D处的慢度。
图7中间的传统差分方法中，其差分格式表示为：
Figure BDA0002796814550000173
为求解该方程，将其带入二维情况下的程函方程
Figure BDA0002796814550000174
解方程，得到待求点C旅行时公式为：
Figure BDA0002796814550000175
图7右侧的改进型差分方法中，其差分格式表示为：
Figure BDA0002796814550000176
为求解该方程，将其带入二维情况下的程函方程
Figure BDA0002796814550000177
解方程，得到待求点C旅行时公式为：
Figure BDA0002796814550000178
其中，
Figure BDA0002796814550000179
根据波前扩展规律，式(4)进一步推导为
Figure BDA0002796814550000181
Figure BDA0002796814550000182
利用上文中公式的公式(2)(3)(5)，求得待求点C的多个旅行时后，根据费马原理，选择极小旅行时作为待求点C的旅行时。按照以上方法更新邻点为窄带点类型和远点类型的网格节点的旅行时。
重复以上寻找窄带中极小旅行时点，判断其目标邻点的节点类型，更新邻点初至波旅行时的操作，直至窄带为空，也即，所有网格节点的节点类型均为近点，即完成了所有网格节点的初至波旅行时的计算。
在一个可选的实施方式中，上述步骤S1042，基于所属于透反射界面的若干离散点的空间位置和所有网格节点的初至波旅行时确定界面初至旅行时，具体包括如下步骤：
步骤S10421，基于所属于透反射界面的若干离散点的空间位置和插值法离散化透反射界面，得到透反射界面上所有离散点的空间位置。
具体的，通过上文中的描述可知，给定速度模型的透反射界面时仅给定透反射界面中的若干离散点，要确定透反射界面上所有离散点的空间位置，可以选择采用三次样条插值方法离散化透反射界面，本发明实施例不对插值方法进行具体的限定，用户也可以选择其他插值法求解透反射界面上所有离散点的空间位置。
步骤S10422，利用算式tp＝tm+(tn-tm)L/|MN|计算目标离散点的初至波旅行时；其中，目标离散点为透反射界面上所有离散点中的任意一个，tp表示目标离散点的初至波旅行时，tm和tn分别表示目标离散点所在的网格边的两个端点m，n的初至波旅行时，且tn>tm，L表示目标离散点与端点m之间的距离，0≤L≤|MN|，|MN|表示端点m和端点n之间的距离。
在求解透反射界面上所有离散点的初至波旅行时之前，已经确定了所有网格节点的初至波旅行时，所以采用上述算式表示的线性插值方法就可求得透反射界面上任一个离散点的初至波旅行时。
步骤S10423，基于透反射界面上所有离散点的初至波旅行时确定界面初至旅行时。
在一个可选的实施方式中，上述步骤S1043，基于界面初至旅行时和预设计算分区计算后续波旅行时，具体包括如下步骤：
步骤S10431，判断后续波是否为初次经过目标透反射界面；其中，目标透反射界面为所有透反射界面中的任意一个界面。
若是，则执行步骤S10432；若否，则执行步骤S10433。
步骤S10432，基于目标透反射界面的界面初至旅行时计算目标计算区域内所有网格节点的后续波旅行时；其中，目标计算区域为预设计算分区中后续波的下一传输方向对应的计算区域。
步骤S10433，基于目标透反射界面的目标界面旅行时计算目标计算区域内所有网格节点的后续波旅行时；其中，目标界面旅行时包括利用后续波上一计算区域内的所有网格节点的后续波旅行时插值得到的目标透反射界面上所有离散点的后续波旅行时。
为了计算后续波旅行时，在本发明实施例采用分区多步思想来计算，分区多步是指根据后续波的类型划分计算区域，分步计算各类波的旅行时场，以实现多波射线追踪。
计算后续波旅行时的时候，需要判断判断后续波是否为初次经过目标透反射界面，如果是，那么可以根据目标透反射界面的初至旅行时来计算后续波的下一传输方向对应的计算区域内网格节点的后续波旅行时；如果不是初次经过目标透反射界面，则首先需要根据后续波在上一计算区域内网格节点的旅行时通过插值法计算目标界面旅行时，然后再基于目标界面旅行时来计算后续波的下一传输方向对应的计算区域内所有网格节点的后续波旅行时。
图9提供了一种分区多步射线追踪的示意图，以图9中的两层速度模型为例，追踪透射波和反射波的实现过程可以描述为：
步骤(1)，离散化速度模型区域及透反射界面，按波类将划分计算区域Ⅰ和Ⅱ，由震源点开始计算区域Ⅰ中网格节点的初至波旅行时，如图9中的(a)视图；
步骤(2)，根据透反射界面所处的网格单元，利用插值法计算界面旅行时，如图9中的(b)视图；(计算方法可参考上述步骤S10422)。
步骤(3)，依次计算各个计算区域的上行波或下行波的波前，其处理方法为：①若波在透反射界面处发生透射，则由界面旅行时计算区域Ⅱ的下行波得到透射波波前，如图9中的(c)视图；②若波在透反射界面处发生反射，则由界面旅行时计算区域Ⅰ的上行波得到反射波或反射转换波的波前，如图9中的(d)视图。若①②中所追踪的波的类型发生转换，则在计算下一计算区域的上行波或下行波的波前时选用转换波的速度，则可得到透射波或反射波的波前。
步骤(4)，根据步骤(1)中的计算区域，重复步骤(2)、步骤(3)，即可实现任意次透射、反射及反射转换等各类型多次波的波前，多层的复杂速度模型亦然。
在一个可选的实施方式中，上述步骤S106，基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径，具体包括如下步骤：
步骤S1061，追踪目标波的射线路径时，针对射线路径中的目标射线，利用算式
Figure BDA0002796814550000211
计算目标射线的终点旅行时，以及，利用算式
Figure BDA0002796814550000212
计算目标射线的起点空间位置；其中，若已知接收点的目标波旅行时，则不使用终点旅行时的算式对接收点旅行时进行更新，目标波包括以下任一种：初至波，一次反射波，多次反射波，目标射线为射线路径中的任意一段射线，tc表示目标射线的终点C的目标波旅行时，ta表示目标射线所在的网格中已知离散点A的目标波旅行时，tb表示目标射线所在的网格中已知离散点B的目标波旅行时，
Figure BDA0002796814550000214
表示线段AB与线段AC之间的夹角，v表示波速，|AC|表示已知离散点A和目标射线的终点C之间的距离，|AB|表示已知离散点A和已知离散点B之间距离，(xa,za)表示已知离散点A的空间位置，(xb,zb)表示已知离散点B的空间位置，(xd,zd)表示目标射线的起点空间位置，r表示已知离散点A和目标射线的起点之间的距离，
Figure BDA0002796814550000213
步骤S1062，按照接收点至震源点的方向将所有目标射线依次进行连接，得到目标波的射线路径。
具体的，在本发明实施例中，反向计算接收点到震源点的射线路径时，采用步骤S1061中的算式来推算每一个次级震源点的坐标以及该次级震源点的下一次级震源点的旅行时，需要注意的是，由于是反向推算，所以如果接收点的目标波旅行时已经确定，则无需用步骤S1061中的算式进行计算，如果未知，则需要利用步骤S1061中计算目标射线的终点旅行时的算式求得接收点的目标波旅行时。本发明实施例为了使得追踪得到的射线路径更为精确，还考虑了次级震源点出现在网格内部和网格边界的情况，图10示出了由A和B两点的旅行时求C点旅行时及D点坐标的示意图，设点A、B处的坐标分别为(xa,za)、(xb,zb)及旅行时ta、tb，且点A为极小旅行时点，次级震源点C、D为待求点，射线方向由点D到点C，求解射线路径转化为求解次级震源点D的坐标(xd,zd)和点C的旅行时tc。
基于线性假设前提，网格边界A、B之间任意点D的旅行时td可由A、B两点的旅行时ta、tb通过线性插值得到，即td＝ta+(tb-ta)r/|AB|(6)，其中，点A为极小旅行时点，满足ta≤tb；r为点A、D之间的距离，满足0≤r≤|AB|，|AB|为点A、B之间的距离。由于点C处的射线来自点D，则点C的旅行时tc可表示为
Figure BDA0002796814550000221
v表示波速，
Figure BDA0002796814550000222
表示线段AB与线段AC之间的夹角，|CD|和|AC|分别表示对应点间的距离，将式(6)带入式(7)可得，
Figure BDA0002796814550000223
由式(8)可以看出，点C处的旅行时tc是关于r的函数，根据Fermat原理，由点D传到点C的射线旅行时最小，其旅行时梯度为零，即
Figure BDA0002796814550000224
将式(8)代入式(9)，可得：
Figure BDA0002796814550000225
求解式(10)可得到
Figure BDA0002796814550000226
进而可求得D点的坐标
Figure BDA0002796814550000227
同时。将式(11)代入式(8)可到点C的旅行时tc，即
Figure BDA0002796814550000228
也即，式(11)、(12)、(13)适合任意形状网格下的旅行时计算。
当点A、B位于网格节点时，如图10中左侧和中间的图所示，式(11)、(12)、(13)可进一步简化。
如图10左侧视图所示，当点A、B、D位于网格横边时，即za＝zb，则式(11)、(12)、(13)简化后对应的公式为：
Figure BDA0002796814550000231
Figure BDA0002796814550000232
Figure BDA0002796814550000233
如图10中间视图所示，当点A、B、D位于网格竖边时，即xa＝xb，则式(11)、(12)、(13)简化后对应的公式为：
Figure BDA0002796814550000234
Figure BDA0002796814550000235
Figure BDA0002796814550000236
由接收点向震源检索的射线路径追踪过程如图11所示，图中S为震源，R为接收点，虚线表示来自所有方向可能的射线路径，实线表示确定后的最小旅行时射线路径，其实现过程可分为以下三步：
①根据接收点的位置检索其所在网格单元所有节点的坐标及旅行时，若接收点不在网格节点上，则在确定其上级次级震源点的同时利用式(13)、(16)或(19)计算接收点旅行时。
②根据接收点的位置考虑来自所有方向的射线，利用式(16)或(19)计算各方向旅行时并选取最小旅行时射线段，利用式(15)或(18)计算最小旅行时方向次级震源点位置；若遇到透反射界面，则利用式(11)，(12)，(13)计算次级震源点位置。
③重复①②，直至次级震源点到达震源所在网格为止，依次连接震源→各次级震源点→接收点，得到震源点到接收点的整条射线路径，完成射线追踪。
在一个可选的实施方式中，网格离散化的速度模型中震源点附近的网格采用双重网格技术进行加密，以使速度模型划分为：加密网格区域和非加密网格区域。
具体的，现有技术中网格节点旅行时的计算误差只要集中在震源点附近，因此，本发明实施例中，网格离散化的速度模型中震源点附近的网格采用双重网格技术进行加密，加密后，速度模型划分为：加密网格区域和非加密网格区域，以此提高网格节点旅行时的计算精度。
图12为双重网格技术的示意图，双重网格技术的实现过程为：提取震源点附近某一范围内的网格，并进行加密，如图12所示，利用本发明实施例中的Vidale差分方法计算加密网格区域内所有网格节点处的旅行时，将完成计算的加密网格区域边界上主节点作为原区域的窄带点继续进行计算，直至完成所有网格节点旅行时的计算。
由双重网格技术的实现过程可知，该方法存在两个关键因子，即加密外扩层数和网格单元的加密点数。对于原网格进行加密可以提高计算精度，但势必会增加计算量，降低计算效率。如果加密外扩的层数过少，随着加密点数的增加虽然对计算效率影响不会太大，但对精度的提高不会有明显效果；如果加密外扩的层数过多，随着加密点数的增加，虽然计算精度会提高，但计算效率会越来越低，而且当达到一定的加密层数和点数后，计算精度的提高会越来越慢，不能兼顾效率。因此，要根据精度需要选择合适的加密层数和点数，原则上选择多的加密层数和少的加密点数，在实际操作中要根据设计的速度模型进行不同外扩层数和加密点数的试算，寻找最佳参数，本发明实施例不对双重网格技术中的参数进行具体限制。
射线追踪结束，可以输出各类波的旅行时场(旅行时和空间位置)，还可以利用所有接收点的旅行时场合成地震记录，以供统计或者其他研究使用。
本发明实施例提供的射线追踪方法，使用了波前窄带扩展技术、Vidale差分法和线性插值射线路径计算方法，且计算射线路径时考虑了次级震源点出现在网格内部和网格边界的情况，通过分区多步思想的运用，本发明实施例可以实现复杂模型地表、地井、井地、井间等任意观测系统下初至波及任意类型波任意次透射、反射、反射转换、透射转换的正演模拟，且本发明实施例对复杂模型具有较强的适应能力和稳定性，计算精度和效率较高，有效解决了现有的射线追踪方法存在的计算精度低的技术问题。
为了验证本发明实施例提供的射线追踪方法的计算精度、计算效率、稳定性以及对复杂模型的适应能力，发明人采用了多组典型模型进行测试。
为测试本发明方法的计算精度和效率，发明人设计了一均匀模型，模型区域为1000×1000m，波速为1000m/s。观测系统如图13中的(a)视图所示，▲表示震源，坐标(0，0)m，▼表示接收点，依次排列在模型右边界，道间距100m，接收排列共11道。在不同的网格间距在采用一阶、二阶迎风和Vidale差分，并使用双重网格技术，初至波射线路径追踪如图13所示，图14示出了不同网格间距不同差分格式计算误差的对比示意图，下表提供了不同网格间距下不同差分计算效率对比。
表1不同网格间距下不同差分计算效率对比
Figure BDA0002796814550000251
从图13可以看出，本发明实施例使用传统Vidale差分及改进型Vidale差分和插值法射线路径计算方法是可行的，使用小网格间距时，一阶、二阶迎风差分及Vidale差分格式下的射线路径与真解基本吻合，且射线路径垂直波前面。从图14与表1可以看出，网格间距越小，计算值与真解的之间的误差越小，而采用传统Vidale差分及改进型Vidale差分的计算精度优于一阶及二阶迎风差分，且计算效率与一阶迎风差分相当。从图14中的(b)视图可以看出，采用双重网格技术，在保证计算效率的前提下，计算精度得到明显改善。因此，本发明实施例采用Vidale差分及双重网格技术有效的提高了有限差分法射线追踪的计算精度，同时还保证了计算效率。
进一步的，为了验证本发明实施例中波前窄带扩展方式及射线路径计算的正确性和稳定性，发明人采用Asakawa和Kawanaka设计的高速夹层模型，高速夹层的模型参数、观测系统及射线追踪结果如图15所示，模型区域为1000×1000m，背景波速为1000m/s，高速夹层波速为1500m/s，模型a与模型b呈90°对称；▲表示震源，▼表示接收点，道间距100m，接收排列共11道，采用Vidale差分及与双重网格技术。
从图15可以看出，本发明得到的初至波旅行时等值线沿高速区弯曲，射线路径沿波前法向方向，并且图15中(a)、(b)视图中第3、4道接收的射线不是直接来自震源，而是来自高速层的折射，满足Fermat原理，与Asakawa和Kawanaka测试结果一致，表明本发明计算的旅行时和射线路径是正确的，波前窄带扩展方式稳定性较强，优于波前盒式扩展方式。针对两种高速夹层速度模型，利用波前窄带扩展方式和线性插值法准确地得到了正确结果，表明本发明采用的波前窄带扩展方式及射线路径计算方法是正确的、稳定的。
进一步的，为了验证本发明对复杂模型的适应能力和稳定性，发明人对Marmousi模型进行试算，对比分析本发明与LTI(Linear Traveltime Interpolation，旅行时线性插值法)射线追踪结果。模型区域3830m×1210m，网格间距为10m，网格节点数384×122(Nx×Nz)，速度参数参考图16(a)，图16(a)还示出了使用本发明方法进行射线追踪的结果示意图；震源▲位于模型底界面，坐标(2790，1210)m，接收排列位于模型顶部，道间距50m，共77道。本发明采用Vidale差分，并在震源周围进行100层×4点双重网格加密，LTI采用4×4的次生节点，图16(b)为本发明实施例提供的一种使用LTI技术进行射线追踪的结果示意图。
从图16(a)可以看出，本发明的射线追踪结果较为准确的反映了Marmousi模型的速度结构，高速区的波前较低速区超前，射线在高速区聚集而绕开低速区，且在速度突变点发生散射，射线总是沿最小旅行时方向传播，满足Fermat原理。对比图16(a)、图16(b)，本发明与LTI的旅行时场及射线路径基本吻合；同时，图16(c)为本发明实施例提供的一种本发明方法与LTI技术的旅行时对比示意图，从图16(c)中也可以看出本发明与LTI的旅行时高度吻合，两者残差几乎为零。由此可见，本发明对复杂模型的适应能力较强，且算法稳定性较好。
进一步的，本发明实施例采用了分区多步思想进行射线追踪，为了验证本发明方法对复杂模型多波射线追踪的能力，发明人采用Marmousi模型进行试算。模型参数与上一验证方法中的模型参数相同，震源▲与接收排列▼位于模型顶界面，震源坐标(1915，0)m，第一道接收点坐标(0,0)m，道间距383m，共11道接收。图17为本发明实施例提供的一种Marmousi模型参数及反射追踪模型分区示意图，如图17所示，在模型中给定两个透反射界面，并划分一次反射波及多次波的计算区域，其中一次反射波计算分区为Ⅰ-Ⅱ及Ⅰ-Ⅲ，多次波的计算分区为1–2–3–4。采用与上一验证方法中相同的差分格式及双重网格进行各类波射线追踪，其结果如图18至图20所示，图18为一次反射波追踪示意图，图19为多次反射波追踪示意图，图20为Marmousi模型单炮合成地震记录的示意图。
从图18至图19可以看出，本发明方法很好地解决了复杂模型反射波和多次波追踪问题，波前面及射线路径满足波的传播规律。结果表明，本发明方法成功实现了复杂模型各类波射线路径追踪，且对复杂模型的适应能力较强，稳定性较好。图20是Marmousi模型单炮合成地震记录的示意图，炮点位置不变，地表154道接收，道间距25m，合成地震记录中清晰显示了各接收道的初至波及来自两条界面的一次反射波及多次波。
进一步的，为了验证本发明方法对不同观测系统的适应能力，发明人设置了四种观测系统，即地表、地井、井地、井间观测系统进行正演模拟。
图21为本发明实施例提供的一种速度模型及反射界面的示意图，模型参数为：模型大小1000m×1500m，网格间距为2m，其速度和界面参数参考图21，模型的速度类型为P波，该模型共有三条透反射界面，第一条透反射界面为断层，第二条透反射界面为水平界面，第三条透反射界面为起伏界面。三条透反射界面将模型分为四层，第一层为均匀速度层，层速为2000m/s；第二层为含低速异常体的高速层，层速为2500m/s，透镜体的速度为2000m/s；第三层为含高速异常体的低速层，层速为1500m/s，正方块体的速度为2000m/s；第四层为均匀速度层，层速为2000m/s。
设计了地表、地井、井地、井间观测系统追踪各界面的一次反射波，各观测系统中的接收道位置固定。单炮模拟中道间距为50m，其中地表炮点坐标为(0,0)，接收排列为21道，起始坐标为(0,0)；井中炮点坐标为(0,150)，接收排列为31道，起始坐标为(1000,0)。多炮模拟中炮数为3，道间距为5m，其中地表炮间距为250m，炮点起始坐标为(0,0)，接收排列为201道，起始坐标为(0,0)；井中炮间距300m，炮点起始坐标为(0,150)，接收排列为301道，起始坐标为(1000,0)。
图22(a)为本发明实施例提供的一种地表观测系统射线路径追踪结果示意图；图22(b)为本发明实施例提供的一种地井观测系统射线路径追踪结果示意图；图22(c)为本发明实施例提供的一种井地观测系统射线路径追踪结果示意图；图22(d)为本发明实施例提供的一种井间观测系统射线路径追踪结果示意图。
图23提供了不同观测系统下的三炮合成记录示意图，图23中(a)视图为本发明实施例提供的一种地表观测系统合成记录示意图；图23中(b)视图为本发明实施例提供的一种地井观测系统合成记录示意图；图23中(c)视图为本发明实施例提供的一种井地观测系统合成记录示意图；图23中(d)视图为本发明实施例提供的一种井间观测系统合成记录示意图。
由不同观测系统下的模型测试可知，本发明提供的射线追踪方法适应地表、地井、井地、井间观测系统初至波与反射波的追踪，且对复杂模型具有较强的适应能力。
通过以上不同模型试算说明，本发明方法能够实现复杂模型地表、地井、井地、井间等任意观测系统下初至波及任意类型波任意次透射、反射、反射转换、透射转换的正演模拟，灵活性及稳定性较强，且本发明方法的计算精度较高，对复杂模型具有较强的适应能力。
实施例二
本发明实施例还提供了一种射线追踪装置，该射线追踪装置主要用于执行上述实施例一所提供的射线追踪方法，以下对本发明实施例提供的射线追踪装置做具体介绍。
图24是本发明实施例提供的一种射线追踪装置的功能模块图，如图24所示，该装置主要包括：获取模块10，第一计算模块20，第二计算模块30，其中：
获取模块10，用于获取网格离散化的速度模型和射线追踪观测系统的系统参数；其中，速度模型为射线追踪区域，系统参数包括震源点的空间位置和接收点的空间位置。
第一计算模块20，用于基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时；其中，目标旅行时包括：初至波旅行时和后续波旅行时，后续波包括：一次反射波和多次反射波。
第二计算模块30，用于基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径；其中，射线路径中目标震源点的位置包括以下任一种：网格边界、网格节点、网格内部，目标震源点为射线路径中所有次级震源点中的任意一个。
现有技术中的射线追踪方法在进行射线追踪时仅考虑所有次级震源点位于网格节点的情况，从而难以保证射线路径的计算精度。与现有技术相比，本发明提供的射线追踪装置，能够获取网格离散化的速度模型和射线追踪观测系统的系统参数；然后基于网格离散化的速度模型和系统参数计算射线追踪区域内所有离散点的目标旅行时；最后基于所有离散点的目标旅行时和系统参数反向计算接收点到震源点的射线路径。本发明装置在进行射线追踪时考虑了次级震源点出现在网格内部和网格边界的情况，有效的提高了射线路径的计算精度，从而缓解了现有技术中的射线追踪方法存在的计算精度低的技术问题。
可选的，速度模型的模型参数包括：所属于透反射界面的若干离散点的空间位置；震源点的初至波旅行时为0；第一计算模块20包括：
第一计算单元，用于利用窄带扩展方法依次计算网格离散化的速度模型内所有网格节点的初至波旅行时。
第一确定单元，用于基于所属于透反射界面的若干离散点的空间位置和所有网格节点的初至波旅行时确定界面初至旅行时；其中，界面初至旅行时包括透反射界面上所有离散点的初至波旅行时。
第二计算单元，用于基于界面初至旅行时和预设计算分区计算后续波旅行时；其中，预设计算分区为速度模型根据后续波的类型划分的计算区域。
第二确定单元，用于基于初至波旅行时和后续波旅行时确定目标旅行时。
可选的，第一计算单元具体用于：
将网格离散化的速度模型中的所有网格节点进行分类，得到分类后的网格节点和窄带；其中，分类后的网格节点的节点类型包括：近点，窄带点和远点，震源点为近点，在未进行窄带扩展时窄带点仅包括震源点的上下左右四个邻点，远点为除近点和窄带点之外的网格节点；窄带为所有窄带点的集合。
基于震源点的初至波旅行时计算所有窄带点的初至波旅行时。
重复执行以下步骤，直至所有网格节点的节点类型均为近点：基于所有窄带点的初至波旅行时确定窄带中的最小旅行时点，并将最小旅行时点更新为近点；判断最小旅行时点的目标邻点的节点类型；其中，目标邻点为最小旅行时点的上下左右四个邻点中的任意一个；若目标邻点为近点，则目标邻点的初至波旅行时保持不变；若目标邻点为窄带点或远点，则基于最小旅行时点的初至波旅行时更新目标邻点的初至波旅行时，且将节点类型为远点的目标邻点移入窄带。
可选的，第一确定单元具体用于：
基于所属于透反射界面的若干离散点的空间位置和插值法离散化透反射界面，得到透反射界面上所有离散点的空间位置。
利用算式tp＝tm+(tn-tm)L/|MN|计算目标离散点的初至波旅行时；其中，目标离散点为透反射界面上所有离散点中的任意一个，tp表示目标离散点的初至波旅行时，tm和tn分别表示目标离散点所在的网格边的两个端点m，n的初至波旅行时，且tn>tm，L表示目标离散点与端点m之间的距离，0≤L≤|MN|，|MN|表示端点m和端点n之间的距离。
基于透反射界面上所有离散点的初至波旅行时确定界面初至旅行时。
可选的，第二计算单元具体用于：
判断后续波是否为初次经过目标透反射界面；其中，目标透反射界面为所有透反射界面中的任意一个界面。
若是，则基于目标透反射界面的界面初至旅行时计算目标计算区域内所有网格节点的后续波旅行时；其中，目标计算区域为预设计算分区中后续波的下一传输方向对应的计算区域。
若否，则基于目标透反射界面的目标界面旅行时计算目标计算区域内所有网格节点的后续波旅行时；其中，目标界面旅行时包括利用后续波上一计算区域内的所有网格节点的后续波旅行时插值得到的目标透反射界面上所有离散点的后续波旅行时。
可选的，第二计算模块30具体用于：
追踪目标波的射线路径时，针对射线路径中的目标射线，利用算式
Figure BDA0002796814550000321
计算目标射线的终点旅行时，以及，利用算式
Figure BDA0002796814550000322
计算目标射线的起点空间位置；其中，若已知接收点的目标波旅行时，则不使用终点旅行时的算式对接收点旅行时进行更新，目标波包括以下任一种：初至波，一次反射波，多次反射波，目标射线为射线路径中的任意一段射线，tc表示目标射线的终点C的目标波旅行时，ta表示目标射线所在的网格中已知离散点A的目标波旅行时，tb表示目标射线所在的网格中已知离散点B的目标波旅行时，
Figure BDA0002796814550000324
表示线段AB与线段AC之间的夹角，v表示波速，|AC|表示已知离散点A和目标射线的终点C之间的距离，|AB|表示已知离散点A和已知离散点B之间距离，(xa,za)表示已知离散点A的空间位置，(xb,zb)表示已知离散点B的空间位置，(xd,zd)表示目标射线的起点空间位置，r表示已知离散点A和目标射线的起点之间的距离，
Figure BDA0002796814550000323
按照接收点至震源点的方向将所有目标射线依次进行连接，得到目标波的射线路径。
可选的，网格离散化的速度模型中震源点附近的网格采用双重网格技术进行加密，以使速度模型划分为：加密网格区域和非加密网格区域。
实施例三
参见图25，本发明实施例提供了一种电子设备，该电子设备包括：处理器60，存储器61，总线62和通信接口63，所述处理器60、通信接口63和存储器61通过总线62连接；处理器60用于执行存储器61中存储的可执行模块，例如计算机程序。
其中，存储器61可能包含高速随机存取存储器(RAM，RandomAccessMemory)，也可能还包括非不稳定的存储器(non-volatile memory)，例如至少一个磁盘存储器。通过至少一个通信接口63(可以是有线或者无线)实现该系统网元与至少一个其他网元之间的通信连接，可以使用互联网，广域网，本地网，城域网等。
总线62可以是ISA总线、PCI总线或EISA总线等。所述总线可以分为地址总线、数据总线、控制总线等。为便于表示，图25中仅用一个双向箭头表示，但并不表示仅有一根总线或一种类型的总线。
其中，存储器61用于存储程序，所述处理器60在接收到执行指令后，执行所述程序，前述本发明实施例任一实施例揭示的流过程定义的装置所执行的方法可以应用于处理器60中，或者由处理器60实现。
处理器60可能是一种集成电路芯片，具有信号的处理能力。在实现过程中，上述方法的各步骤可以通过处理器60中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器60可以是通用处理器，包括中央处理器(Central Processing Unit，简称CPU)、网络处理器(Network Processor，简称NP)等；还可以是数字信号处理器(Digital SignalProcessing，简称DSP)、专用集成电路(Application Specific Integrated Circuit，简称ASIC)、现成可编程门阵列(Field-Programmable Gate Array，简称FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。结合本发明实施例所公开的方法的步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器61，处理器60读取存储器61中的信息，结合其硬件完成上述方法的步骤。
本发明实施例所提供的一种射线追踪方法、装置和电子设备的计算机程序产品，包括存储了处理器可执行的非易失的程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
此外，术语“水平”、“竖直”、“悬垂”等术语并不表示要求部件绝对水平或悬垂，而是可以稍微倾斜。如“水平”仅仅是指其方向相对“竖直”而言更加水平，并不是表示该结构一定要完全水平，而是可以稍微倾斜。
在本发明的描述中，还需要说明的是，除非另有明确的规定和限定，术语“设置”、“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (10)

1.一种射线追踪方法，其特征在于，包括：
获取网格离散化的速度模型和射线追踪观测系统的系统参数；其中，所述速度模型为射线追踪区域，所述系统参数包括震源点的空间位置和接收点的空间位置；
基于所述网格离散化的速度模型和所述系统参数计算所述射线追踪区域内所有离散点的目标旅行时；其中，所述目标旅行时包括：初至波旅行时和后续波旅行时，所述后续波包括：一次反射波和多次反射波；
基于所述所有离散点的目标旅行时和所述系统参数反向计算所述接收点到所述震源点的射线路径；其中，所述射线路径中目标震源点的位置包括以下任一种：网格节点、网格边界、网格内部，所述目标震源点为所述射线路径中所有次级震源点中的任意一个。
2.根据权利要求1所述的方法，其特征在于，所述速度模型的模型参数包括：所属于透反射界面的若干离散点的空间位置；所述震源点的初至波旅行时为0；
基于所述网格离散化的速度模型和所述系统参数计算所述射线追踪区域内所有离散点的目标旅行时，包括：
利用窄带扩展方法依次计算所述网格离散化的速度模型内所有网格节点的初至波旅行时；
基于所述所属于透反射界面的若干离散点的空间位置和所述所有网格节点的初至波旅行时确定界面初至旅行时；其中，所述界面初至旅行时包括所述透反射界面上所有离散点的初至波旅行时；
基于所述界面初至旅行时和预设计算分区计算所述后续波旅行时；其中，所述预设计算分区为所述速度模型根据所述后续波的类型划分的计算区域；
基于所述初至波旅行时和所述后续波旅行时确定所述目标旅行时。
3.根据权利要求2所述的方法，其特征在于，利用窄带扩展方法依次计算所述网格离散化的速度模型内所有网格节点的初至波旅行时，包括：
将所述网格离散化的速度模型中的所有网格节点进行分类，得到分类后的网格节点和窄带；其中，所述分类后的网格节点的节点类型包括：近点，窄带点和远点，所述震源点为近点，在未进行窄带扩展时所述窄带点仅包括震源点的上下左右四个邻点，所述远点为除所述近点和所述窄带点之外的网格节点；所述窄带为所有窄带点的集合；
基于所述震源点的初至波旅行时计算所有窄带点的初至波旅行时；
重复执行以下步骤，直至所有网格节点的节点类型均为近点：
基于所有窄带点的初至波旅行时确定所述窄带中的最小旅行时点，并将所述最小旅行时点更新为近点；
判断所述最小旅行时点的目标邻点的节点类型；其中，所述目标邻点为所述最小旅行时点的上下左右四个邻点中的任意一个；
若所述目标邻点为近点，则所述目标邻点的初至波旅行时保持不变；
若所述目标邻点为窄带点或远点，则基于所述最小旅行时点的初至波旅行时更新所述目标邻点的初至波旅行时，且将节点类型为远点的目标邻点移入所述窄带。
4.根据权利要求2所述的方法，其特征在于，基于所述所属于透反射界面的若干离散点的空间位置和所述所有网格节点的初至波旅行时确定界面初至旅行时，包括：
基于所述所属于透反射界面的若干离散点的空间位置和插值法离散化所述透反射界面，得到所述透反射界面上所有离散点的空间位置；
利用算式tp＝tm+(tn-tm)L/|MN|计算目标离散点的初至波旅行时；其中，所述目标离散点为所述透反射界面上所有离散点中的任意一个，tp表示所述目标离散点的初至波旅行时，tm和tn分别表示所述目标离散点所在的网格边的两个端点m，n的初至波旅行时，且tn>tm，L表示所述目标离散点与端点m之间的距离，0≤L≤|MN|，|MN|表示端点m和端点n之间的距离；
基于所述透反射界面上所有离散点的初至波旅行时确定所述界面初至旅行时。
5.根据权利要求2所述的方法，其特征在于，基于所述界面初至旅行时和预设计算分区计算所述后续波旅行时，包括：
判断后续波是否为初次经过目标透反射界面；其中，所述目标透反射界面为所有透反射界面中的任意一个界面；
若是，则基于所述目标透反射界面的界面初至旅行时计算目标计算区域内所有网格节点的后续波旅行时；其中，所述目标计算区域为所述预设计算分区中所述后续波的下一传输方向对应的计算区域；
若否，则基于所述目标透反射界面的目标界面旅行时计算所述目标计算区域内所有网格节点的后续波旅行时；其中，所述目标界面旅行时包括利用所述后续波上一计算区域内的所有网格节点的后续波旅行时插值得到的所述目标透反射界面上所有离散点的后续波旅行时。
6.根据权利要求1所述的方法，其特征在于，基于所述所有离散点的目标旅行时和所述系统参数反向计算所述接收点到所述震源点的射线路径，包括：
追踪目标波的射线路径时，针对所述射线路径中的目标射线，利用算式
Figure FDA0002796814540000031
计算所述目标射线的终点旅行时，以及，利用算式
Figure FDA0002796814540000032
计算所述目标射线的起点空间位置；其中，若已知所述接收点的目标波旅行时，则不使用所述终点旅行时的算式对所述接收点旅行时进行更新，所述目标波包括以下任一种：初至波，一次反射波，多次反射波，所述目标射线为所述射线路径中的任意一段射线，tc表示所述目标射线的终点C的目标波旅行时，ta表示所述目标射线所在的网格中已知离散点A的目标波旅行时，tb表示所述目标射线所在的网格中已知离散点B的目标波旅行时，
Figure FDA0002796814540000041
表示线段AB与线段AC之间的夹角，v表示波速，|AC|表示所述已知离散点A和所述目标射线的终点C之间的距离，|AB|表示所述已知离散点A和所述已知离散点B之间距离，(xa,za)表示所述已知离散点A的空间位置，(xb,zb)表示所述已知离散点B的空间位置，(xd,zd)表示所述目标射线的起点空间位置，r表示所述已知离散点A和所述目标射线的起点之间的距离，
Figure FDA0002796814540000042
按照接收点至震源点的方向将所有目标射线依次进行连接，得到所述目标波的射线路径。
7.根据权利要求1所述的方法，其特征在于，所述网格离散化的速度模型中震源点附近的网格采用双重网格技术进行加密，以使所述速度模型划分为：加密网格区域和非加密网格区域。
8.一种射线追踪装置，其特征在于，包括：
获取模块，用于获取网格离散化的速度模型和射线追踪观测系统的系统参数；其中，所述速度模型为射线追踪区域，所述系统参数包括震源点的空间位置和接收点的空间位置；
第一计算模块，用于基于所述网格离散化的速度模型和所述系统参数计算所述射线追踪区域内所有离散点的目标旅行时；其中，所述目标旅行时包括：初至波旅行时和后续波旅行时，所述后续波包括：一次反射波和多次反射波；
第二计算模块，用于基于所述所有离散点的目标旅行时和所述系统参数反向计算所述接收点到所述震源点的射线路径；其中，所述射线路径中目标震源点的位置包括以下任一种：网格边界、网格节点、网格内部，所述目标震源点为所述射线路径中所有次级震源点中的任意一个。
9.一种电子设备，包括存储器、处理器，所述存储器上存储有可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述权利要求1至7中任一项所述的方法的步骤。
10.一种具有处理器可执行的非易失的程序代码的计算机可读介质，其特征在于，所述程序代码使所述处理器执行权利要求1至7中任一项所述的方法。
CN202011334054.3A 2020-11-24 2020-11-24 射线追踪方法、装置和电子设备 Pending CN112596103A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011334054.3A CN112596103A (zh) 2020-11-24 2020-11-24 射线追踪方法、装置和电子设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011334054.3A CN112596103A (zh) 2020-11-24 2020-11-24 射线追踪方法、装置和电子设备

Publications (1)

Publication Number Publication Date
CN112596103A true CN112596103A (zh) 2021-04-02

Family

ID=75183704

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011334054.3A Pending CN112596103A (zh) 2020-11-24 2020-11-24 射线追踪方法、装置和电子设备

Country Status (1)

Country Link
CN (1) CN112596103A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114429047A (zh) * 2022-01-27 2022-05-03 成都理工大学 一种基于三角网格的二次圆方程旅行时插值方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6035256A (en) * 1997-08-22 2000-03-07 Western Atlas International, Inc. Method for extrapolating traveltimes across shadow zones
CN1712991A (zh) * 2004-06-25 2005-12-28 中国石油化工股份有限公司 一种用于地震勘探中射线追踪的方法
CA2677135A1 (en) * 2007-02-06 2008-08-14 Alexander Kostyukevych Method of surface seismic imaging using both reflected and transmitted waves
CN106353793A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种基于走时增量双线性插值射线追踪的井间地震层析反演方法
CN106932821A (zh) * 2015-12-31 2017-07-07 上海青凤致远地球物理地质勘探科技有限公司 地震层析反演中的一种目标射线追踪技术
CN107843922A (zh) * 2017-12-25 2018-03-27 中国海洋大学 一种基于地震初至波和反射波走时联合的层析成像方法
CN109444955A (zh) * 2019-01-09 2019-03-08 中国海洋大学 三维地震射线追踪的双线性走时扰动插值方法
WO2019071504A1 (zh) * 2017-10-12 2019-04-18 南方科技大学 一种基于两点射线追踪的地震走时层析反演方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6035256A (en) * 1997-08-22 2000-03-07 Western Atlas International, Inc. Method for extrapolating traveltimes across shadow zones
CN1712991A (zh) * 2004-06-25 2005-12-28 中国石油化工股份有限公司 一种用于地震勘探中射线追踪的方法
CA2677135A1 (en) * 2007-02-06 2008-08-14 Alexander Kostyukevych Method of surface seismic imaging using both reflected and transmitted waves
CN106353793A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种基于走时增量双线性插值射线追踪的井间地震层析反演方法
CN106932821A (zh) * 2015-12-31 2017-07-07 上海青凤致远地球物理地质勘探科技有限公司 地震层析反演中的一种目标射线追踪技术
WO2019071504A1 (zh) * 2017-10-12 2019-04-18 南方科技大学 一种基于两点射线追踪的地震走时层析反演方法
CN107843922A (zh) * 2017-12-25 2018-03-27 中国海洋大学 一种基于地震初至波和反射波走时联合的层析成像方法
CN109444955A (zh) * 2019-01-09 2019-03-08 中国海洋大学 三维地震射线追踪的双线性走时扰动插值方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李永博: "VTI介质及复杂模型FMM射线追踪方法研究", 《中国优秀博硕士学位论文全文数据库（硕士） 基础科学辑，2013年第S2期》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114429047A (zh) * 2022-01-27 2022-05-03 成都理工大学 一种基于三角网格的二次圆方程旅行时插值方法
CN114429047B (zh) * 2022-01-27 2023-08-22 成都理工大学 一种基于三角网格的二次圆方程旅行时插值方法

Similar Documents

Publication Publication Date Title
US4736347A (en) Multiple stacking and spatial mapping of seismic data
US9291735B2 (en) Probablistic subsurface modeling for improved drill control and real-time correction
WO2019242045A9 (zh) 一种虚拟震源二维波前构建地震波走时计算方法
CN107843922A (zh) 一种基于地震初至波和反射波走时联合的层析成像方法
CN107817516B (zh) 基于初至波信息的近地表建模方法及系统
CN104570106A (zh) 一种近地表层析速度分析方法
Le Bouteiller et al. A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media
CN109444955A (zh) 三维地震射线追踪的双线性走时扰动插值方法
CN111948708B (zh) 一种浸入边界起伏地表地震波场正演模拟方法
US11397273B2 (en) Full waveform inversion in the midpoint-offset domain
CN106680870A (zh) 高精度地震波走时射线追踪方法
WO2014074173A1 (en) System and method for analysis of designs of a seismic survey
CN112596103A (zh) 射线追踪方法、装置和电子设备
CA2412995C (en) Seismic survey system
CN105204063B (zh) 地震数据速度模型建立方法和装置
CN105353406A (zh) 一种生成角道集的方法和装置
CN109581494B (zh) 叠前偏移方法及装置
CN117331119A (zh) 一种面向隧道探测的快速地震波走时计算方法
US5265068A (en) Efficient generation of traveltime tables for complex velocity models for depth migration
AU739128B2 (en) A method of seismic processing, and in particular a 3D seismic prospection method implementing seismic data migration
CN109188522B (zh) 速度场构建方法及装置
Li et al. A 3D reflection ray-tracing method based on linear traveltime perturbation interpolation
CN110824558A (zh) 一种地震波数值模拟方法
CN112257241B (zh) 一种三角网菲涅尔带时差层析反演方法
Pleshkevich et al. Kirchhoff-Type Implementation of Multi-Arrival 3-D Seismic Depth Migration with Amplitudes Preserved

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