CN117331119A - 一种面向隧道探测的快速地震波走时计算方法 - Google Patents

一种面向隧道探测的快速地震波走时计算方法 Download PDF

Info

Publication number
CN117331119A
CN117331119A CN202311556673.0A CN202311556673A CN117331119A CN 117331119 A CN117331119 A CN 117331119A CN 202311556673 A CN202311556673 A CN 202311556673A CN 117331119 A CN117331119 A CN 117331119A
Authority
CN
China
Prior art keywords
point
points
travel time
array
grid
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
CN202311556673.0A
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.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
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 Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN202311556673.0A priority Critical patent/CN117331119A/zh
Publication of CN117331119A publication Critical patent/CN117331119A/zh
Pending legal-status Critical Current

Links

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/303—Analysis for determining velocity profiles or travel times
    • 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/303—Analysis for determining velocity profiles or travel times
    • G01V1/305—Travel times
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

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

一种面向隧道探测的快速地震波走时计算方法
技术领域
本发明涉及地震波走时计算技术领域，具体涉及一种面向隧道探测的快速地震波走时计算方法。
背景技术
反射地震波成像方法是主要的隧道超前探测方法，地震波走时是其核心参数，常被广泛的应用于偏移与层析成像等技术中。现有的走时计算方法基本都是基于规则模型提出的，而现实中处于施工期的隧道常常拥有不规则的界面，此时要求地震波走时算法不仅要对复杂地表有良好的适应性而且在处理不规则边界时仍然有较高的计算效率，因此，研究面向隧道探测的快速地震波走时计算方法具有重要意义。
《勘探地球物进展》2007年第30卷第5期公开了孙章庆等“波前快速推进法起伏地表地震波走时计算”介绍了波前快速推进法(FMM)是一种准确而稳定的地震波走时计算方法，并且以FMM为基础，提出了一种复杂地表条件的计算地震波走时的计算方法，计算取得了良好的效果，其算法的计算复杂度为O(Nlog2N),N表示计算网格点总数。
中国发明专利申请201110397445.4中公开了《VTI介质中地震波的走时计算方法》，该算法通过输入地下地质速度模型以及观测系统参数，然后计算相应中间结果，主要采用多项式求解的方式，获取地震波走时，通过不同模型的验证计算结果取得了良好的效果。
中国发明专利申请201810077621.8中公开了《一种混合二维地震波走时计算方法》，该算法通过读入速度模型与相关计算参数，从炮点沿不同方向追踪一定距离的射线信息，利用波前构建法计算射线范围内地震波走时，然后利用快速推进法计算剩余网格节点的走时，通过数值模拟证明了该方法良好的计算效果。
通过以上的例子可以看出，现有的地震波走时计算方法针对常规水平模型具有良好的计算效果，但是这些方法在面对复杂起伏地表问题时为了获取较高的计算精度，常需要较长的计算时间，或者为了提高计算效率而相应降低了算法的计算精度。
发明内容
本发明要解决的技术问题是提供一种面向隧道探测的快速地震波走时计算方法，过采用迎风差分方案离散程函方程梯度项，求解程函方程从而获取地震波走时，基于快速推进法的窄带技术，采用一种更快速且更接近实际地震波前传播方式，即在地震波前扩展时每次在窄带中选取一组点作为扩展点，在改善原有走时计算方法的精度的同时，大大提高了计算效率。
为解决上述技术问题，本发明采用如下技术方案：
设计一种面向隧道探测的快速地震波走时计算方法，包括以下步骤：
S1、读入地质模型以及相关参数信息，所述参数包含速度模型的网格点数，网格间距，震源位置和起伏地表信息。
S2、对网格点进行属性划分：将所有网格点分为接收点、窄带点和远离点；接收点指所有网格点中被认为已经走时计算的点；窄带点是指所有窄带内的点；远离点是指除接收点和窄带点以外的所有网格节点，计算过程中用TT表示该计算点的走时值。
S3、震源初始化、窄带初始化和最小走时点选取；
S31、以震源点为起始点将震源点属性设置为idTT＝2，TT＝0，选取震源点的邻点，即上下左右四个点，判断其属性，idTT不为2，则求取其解析值，公式如下：
Ti±1,j＝Si±1,jh,Ti,j±1＝Si,j±1h
其中，h表示网格间距，S表示该点慢度；T表示该点的地震波走时值，i表示该点的横坐标，j表示该点的纵坐标；
S32、将震源临近网格点属性设置为idTT＝1，并将其纳入窄带构成初始窄带，同时将各邻近点的走时TT存入数组GT中；
S33、设置GT数组内最小走时值为TM。
S4、群组G选取：群组G表示在数组GT中走时属于(TM，TM+DT)范围内的点；群组G的所有点都作为窄带中的更新计算点，进行下一步更新计算。
所述步骤S4中群组G的选择标准为：
在二维情况下，设h＝Δx＝Δz,定义
其中，h表示网格间距，x、z分别表示横向与纵向的单位步长，W表示计算波前；表示波前W中的一点，/>表示慢度函数，/>表示走时函数；
则群组G的选择标准如下：
其中，DT表示时间，h表示网格间距，tW,min表示计算波前W中的最小走时点，sW,min表示计算波前W中最小走时点对应的慢度值。这样确定时间的数值，是为了保证，群组G中的元素个数大于1。
S5、GMM采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时。
所述步骤S5中，求取地震波走时的过程为：
网格节点走时是通过迎风差分法求解程函方程获得的，在二维情况下，地震波传播波前满足程函方程：
其中t(x,z)为走时函数，s(x,z)为模型介质慢度函数；
GMM采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时；
其中分别为x和z方向前向和后向差分算子，以一阶差分格式为例，有
S6、窄带扩展演化；
通过逆序遍历数组GT和顺序遍历数组GT，根据数组中的每个点(i,j)的走时，对其所有临近点(x,z)进行属性判断，更新该邻近点(x,z)的走时，将该临近点(x,z)纳入或移出GT数组，直至GT数组为空，输出地震波走时计算结果。
所述步骤S6中窄带扩展演化的具体过程为：
S61、设置
其中，DT表示时间，h表示网格间距，sW,min表示计算波前W中最小走时点对应的慢度值；
S62、设置TM＝TM+DT；
S63、以逆序遍历数组GT，对于数组中的每个点(i,j)如果其走时满足TT(i,j)＜＝TM，则对其所有临近点(x,z)进行属性判断，若idTT<＝1，则更新该邻近点(x,z)的走时；
S64、以顺序遍历数组GT，对于数组GT中的每个点(i,j)如果其走时满足TT(i,j)＜＝TM，则
(ⅰ)对其所有临近点(x,z)进行属性判断，若idTT<＝1，则更新该邻近点(x,z)的走时；
(ⅱ)若其临近点(x,z)属性为idTT＝0，则将其属性改为idTT＝1，并将其临近点(x,z)的走时纳入GT数组中；
(ⅲ)将该点(i,j)的属性设置为idTT＝2，即将该点(i,j)移出数组GT；
S7、判断数组GT是否为空，若为空则输出地震波走时计算结果，若不为空则返回S6继续计算。
本发明的有益效果在于：本发明提供的一种面向隧道探测的快速地震波走时计算方法，通过利用迎风有限差分方案求解程函方程来获取地震波走时，并且基于窄带技术，在地震波前扩展时选择一组点作为扩展点同时进行扩展计算，在保证了计算精度的同时，大大提高了地震波走时算法的计算效率。
隧道探测的环境通常比较复杂，对方法的稳定性要求较高，本发明具有无条件稳定性，有很强的适应能力。另外隧道探测对探测结果的时效性要求比较高，本方法具有很高的计算效率，可以提高隧道地震勘探处理的计算效率。
附图说明
图1为本发明面向隧道探测的快速地震波走时计算方法的流程图；
图2为本发明迎风差分方案示意图；
图3为本发明波前扩展示意图；
图4为本发明实施例中均匀模型的相对误差图；
图5、图6为本发明实施例中不同模型的地震波走时等值线图。
具体实施方式
下面结合实施例来说明本发明的具体实施方式，但以下实施例只是用来详细说明本发明，并不以任何方式限制本发明的范围。在以下实施例中所涉及的设备元件如无特别说明，均为常规设备元件。
实施例1：一种面向隧道探测的快速地震波走时计算方法，参见图1-图6，包括以下步骤：
S1、读入地质模型以及相关参数信息，其中，所述参数包含速度模型的网格点数，网格间距，震源位置和起伏地表信息；
本实施例中，计算相对误差使用均匀速度模型，起伏地表模型为自定义复杂模型,模型大小为1200×800，网格间距10m以下叙述中实施例数据为均匀速度模型，模型大小800×800，网格间距5m，速度为2000m/s震源点(2000m,2000m)。
S2、对网格点进行属性划分：对所有网格点进行分类，将它们分为接收点、窄带点和远离点。接收点指所有网格点中被认为已经走时计算的点(idTT＝2)；窄带点是指所有窄带内的点(idTT＝1)；远离点是指除接收点和窄带点以外的所有网格节点(idTT＝0)。计算过程中用TT表示该计算点的走时值。
S3、震源初始化、窄带初始化和最小走时点选取。
S31、以震源点为起始点将震源点属性设置为idTT＝2，TT＝0，选取震源点的邻点，即上下左右四个点，判断其属性，idTT不为2，则求取其解析值，公式如下：
Ti±1,j＝Si±1,jh,Ti,j±1＝Si,j±1h
其中，h表示网格间距，S表示该点慢度；T表示该点的地震波走时值，i表示该点的横坐标，j表示该点的纵坐标；
S32、将震源临近网格点属性设置为idTT＝1，并将其纳入窄带构成初始窄带，同时将各邻近点的走时TT存入数组GT中；
S33、设置GT数组内最小走时值TM。
S4、群组G选取。
群组G的选择标准：在二维情况下，设h＝Δx＝Δz,定义
其中，表示计算波前W中的一点，/>表示慢度函数，/>表示走时函数
则群组G的选择标准如下：
其中，DT表示时间，h表示网格间距，tW,min表示计算波前W中的最小走时点，sW,min表示计算波前W中最小走时点对应的慢度值。
在本实施例中：
S5、网格节点走时是通过迎风差分法求解程函方程获得的，在二维情况下，地震波传播波前满足程函方程：
其中t(x,z)为走时函数，s(x,z)为模型介质慢度函数。
GMM采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时。
其中分别为x和z方向前向和后向差分算子，以一阶差分格式为例，有
S6、窄带扩展演化；
S61、设置
其中，DT表示时间，h表示网格间距，sW,min表示计算波前W中最小走时点对应的慢度值。
S62、设置TM＝TM+DT；
S63、以逆序遍历数组GT，对于数组中的每个点(i,j)如果其走时满足TT(i,j)＜＝TM，则对其所有临近点(x,z)进行属性判断，若idTT<＝1，则更新该邻近点(x,z)的走时；
S64、以顺序遍历数组GT，对于数组GT中的每个点(i,j)如果其走时满足TT(i,j)＜＝TM，则
(ⅰ)对其所有临近点(x,z)进行属性判断，若idTT<＝1，则更新该邻近点(x,z)的走时；
(ⅱ)若其临近点(x,z)属性为idTT＝0，则将其属性改为idTT＝1，并将其临近点(x,z)的走时纳入GT数组中；
(ⅲ)将该点(i,j)的属性设置为idTT＝2，即将该点(i,j)移出数组GT。
S7、判断数组GT是否为空，若为空则输出地震波走时计算结果，若不为空则返回S6继续计算。图5、图6中的等时线图即为输出结果。
本实施例中针对隧道地质模型的地震波走时计算结果，从中可以看出：首先，本发明能够顺利完成全空间地震走时的计算，说明了其具有很强的适应能力。其次，本发明所用的计算时间是普通有限差分方法的三分之一左右，说明其具有更高的计算效率。分析计算结果可知：走时等值线的分布符合地震波在介质内部的传播规律。在面对速度模型中的断层、背斜时，走时等值线会出现弯曲以及分布逐渐变得稀疏等情况，同时计算结果也证明了本发明算法对复杂地表模型有良好的适应性。
上面结合实施例对本发明作了详细的说明，但是所属技术领域的技术人员能够理解，在不脱离本发明宗旨的前提下，还可以对上述实施例中的各个具体参数进行变更，形成多个具体的实施例，均为本发明的常见变化范围，在此不再一一详述。

Claims (4)

1.一种面向隧道探测的快速地震波走时计算方法，其特征在于，包括以下步骤：
S1、读入地质模型以及相关参数信息，所述参数包含速度模型的网格点数，网格间距，震源位置和起伏地表信息；
S2、对网格点进行属性划分：将所有网格点分为接收点、窄带点和远离点；接收点指所有网格点中被认为已经走时计算的点；窄带点是指所有窄带内的点；远离点是指除接收点和窄带点以外的所有网格节点，计算过程中用TT表示该计算点的走时值；
S3、震源初始化、窄带初始化和最小走时点选取；
S31、以震源点为起始点将震源点属性设置为idTT＝2，TT＝0，选取震源点的邻点，即上下左右四个点，判断其属性，idTT不为2，则求取其解析值，公式如下：
Ti±1,j＝Si±1,jh,Ti,j±1＝Si,j±1h
其中，h表示网格间距，S表示该点慢度；T表示该点的地震波走时值，i表示该点的横坐标，j表示该点的纵坐标；
S32、将震源临近网格点属性设置为idTT＝1，并将其纳入窄带构成初始窄带，同时将各邻近点的走时TT存入数组GT中；
S33、设置GT数组内最小走时值为TM；
S4、群组G选取：群组G表示在数组GT中走时属于(TM，TM+DT)范围内的点；群组G的所有点都作为窄带中的更新计算点，进行下一步更新计算；
S5、GMM采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时；
S6、窄带扩展演化；
通过逆序遍历数组GT和顺序遍历数组GT，根据数组中的每个点(i,j)的走时，对其所有临近点(x,z)进行属性判断，更新该邻近点(x,z)的走时，将该临近点(x,z)纳入或移出GT数组，直至GT数组为空，输出地震波走时计算结果；
S7、判断数组GT是否为空，若为空则输出地震波走时计算结果，若不为空则返回S6继续计算。
2.根据权利要求1所述的面向隧道探测的快速地震波走时计算方法，其特征在于，所述步骤S4中群组G的选择标准为：
在二维情况下，设h＝Δx＝Δz,定义
其中，h表示网格间距，x、z分别表示横向与纵向的单位步长，W表示计算波前；表示计算波前W中的一点，/>表示慢度函数，/>表示走时函数；
则群组G的选择标准如下：
其中，DT表示时间，h表示网格间距，tW,min表示计算波前W中的最小走时点，sW,min表示计算波前W中最小走时点对应的慢度值。
3.根据权利要求1所述的面向隧道探测的快速地震波走时计算方法，其特征在于，所述步骤S5中，求取地震波走时的过程为：
网格节点走时是通过迎风差分法求解程函方程获得的，在二维情况下，地震波传播波前满足程函方程：
|▽t(x,z)|＝s(x,z)
其中t(x,z)为走时函数，s(x,z)为模型介质慢度函数；
GMM采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时；
其中分别为x和z方向前向和后向差分算子，以一阶差分格式为例，有
4.根据权利要求1所述的面向隧道探测的快速地震波走时计算方法，其特征在于，所述步骤S6中窄带扩展演化的具体过程为：
S61、设置
其中，DT表示时间，h表示网格间距，sW,min表示计算波前W中最小走时点对应的慢度值；
S62、设置TM＝TM+DT；
S63、以逆序遍历数组GT，对于数组中的每个点(i,j)如果其走时满足TT(i,j)＜＝TM，则对其所有临近点(x,z)进行属性判断，若idTT<＝1，则更新该邻近点(x,z)的走时；
S64、以顺序遍历数组GT，对于数组GT中的每个点(i,j)如果其走时满足TT(i,j)＜＝TM，则
(ⅰ)对其所有临近点(x,z)进行属性判断，若idTT<＝1，则更新该邻近点(x,z)的走时；
(ⅱ)若其临近点(x,z)属性为idTT＝0，则将其属性改为idTT＝1，并将其临近点(x,z)的走时纳入GT数组中；
(ⅲ)将该点(i,j)的属性设置为idTT＝2，即将该点(i,j)移出数组GT。
CN202311556673.0A 2023-11-21 2023-11-21 一种面向隧道探测的快速地震波走时计算方法 Pending CN117331119A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311556673.0A CN117331119A (zh) 2023-11-21 2023-11-21 一种面向隧道探测的快速地震波走时计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311556673.0A CN117331119A (zh) 2023-11-21 2023-11-21 一种面向隧道探测的快速地震波走时计算方法

Publications (1)

Publication Number Publication Date
CN117331119A true CN117331119A (zh) 2024-01-02

Family

ID=89277649

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311556673.0A Pending CN117331119A (zh) 2023-11-21 2023-11-21 一种面向隧道探测的快速地震波走时计算方法

Country Status (1)

Country Link
CN (1) CN117331119A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117607957A (zh) * 2024-01-24 2024-02-27 南方科技大学 基于等效慢度的快速推进法的地震波走时求解方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117607957A (zh) * 2024-01-24 2024-02-27 南方科技大学 基于等效慢度的快速推进法的地震波走时求解方法及系统
CN117607957B (zh) * 2024-01-24 2024-04-02 南方科技大学 基于等效慢度的快速推进法的地震波走时求解方法及系统

Similar Documents

Publication Publication Date Title
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
Kim et al. 3-D traveltime computation using second-order ENO scheme
EP0750203B1 (en) Subsurface modeling from seismic data and secondary measurements
SG184803A1 (en) Artifact reduction in method of iterative inversion of geophysical data
CA2705505A1 (en) Forming a geological model
CN117331119A (zh) 一种面向隧道探测的快速地震波走时计算方法
CN110058307B (zh) 一种基于快速拟牛顿法的全波形反演方法
CN107817516B (zh) 基于初至波信息的近地表建模方法及系统
CN109459787B (zh) 基于地震槽波全波形反演的煤矿井下构造成像方法及系统
CN116774292B (zh) 一种地震波走时确定方法、系统、电子设备及存储介质
CN109444955A (zh) 三维地震射线追踪的双线性走时扰动插值方法
CN111580163B (zh) 一种基于非单调搜索技术的全波形反演方法及系统
CN103119472B (zh) 利用同时和顺序源方法进行全波形反演的混合方法
CN105573963B (zh) 一种电离层水平不均匀结构重构方法
Alkhalifah et al. An eikonal-based formulation for traveltime perturbation with respect to the source location
CN109946742A (zh) 一种TTI介质中纯qP波地震数据模拟方法
CN111665556A (zh) 地层声波传播速度模型构建方法
CN111948708B (zh) 一种浸入边界起伏地表地震波场正演模拟方法
CN103217715B (zh) 多尺度规则网格层析反演静校正方法
Chambers et al. A practical implementation of wave front construction for 3-D isotropic media
Sun et al. Joint 3D traveltime calculation based on fast marching method and wavefront construction
CN106443829A (zh) 一种近地表模型构建方法及装置
CN110568497B (zh) 一种复杂介质条件下地震初至波旅行时的精确求解方法
CN111914609B (zh) 井震联合叠前地质统计学弹性参数反演方法及装置
CN110376642B (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