CN117892569B - 一种基于初至波走时约束的动态波场传播模拟方法 - Google Patents
一种基于初至波走时约束的动态波场传播模拟方法 Download PDFInfo
- Publication number
- CN117892569B CN117892569B CN202311644145.0A CN202311644145A CN117892569B CN 117892569 B CN117892569 B CN 117892569B CN 202311644145 A CN202311644145 A CN 202311644145A CN 117892569 B CN117892569 B CN 117892569B
- Authority
- CN
- China
- Prior art keywords
- simulation
- wave
- dynamic
- travel time
- time
- 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.)
- Active
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 82
- 238000000034 method Methods 0.000 title claims abstract description 53
- 238000005070 sampling Methods 0.000 claims description 12
- 238000005094 computer simulation Methods 0.000 claims description 11
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000005516 engineering process Methods 0.000 claims description 5
- 238000004806 packaging method and process Methods 0.000 claims description 4
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 5
- 238000007796 conventional method Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 239000011295 pitch Substances 0.000 description 1
- 230000008707 rearrangement Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明属于地震波传播模拟技术领域，具体涉及一种基于初至波走时约束的动态波场传播模拟方法，将求解程函方程所获得初至波走时利用time‑bin技术按时间方向重新排序，将模型网格点空间坐标封存到对应时刻的序列中，并筛选出每一时刻的波前范围，从而设计出每一时刻最合适的有限差分模拟的动态空间窗口范围，对每一时刻的波场模拟范围进行约束，通过有限差分法对地震波传播进行动态模拟，可有效降低波场模拟的计算量。在模拟任意时刻的地震波场时，高阶有限差分算子不需要遍历整个模型空间，仅需要在动态空间窗口内进行模拟。该方法不但兼顾波动方程法能有准确模拟复杂模型构造波场的优势，还可降低地震波传播模拟过程的计算量和内存消耗。
Description
技术领域
本发明属于地震勘探服务中的地震波传播模拟技术领域，具体涉及一种基于初至波走时约束的动态波场传播模拟方法。
背景技术
地震波传播模拟主要可以分为射线法和波动方程法两类。相比于射线类方法，波动方程法由于不需要高频近似假设条件，可以直接对波动方程进行求解，因此，在复杂构造中的地震波传播模拟结果更加准确，并且，可以得到全部的波场信息。有限差分法是众多波动方程模拟方法中最为常用的方法，具有占用内存小、计算效率高、易于编程实现及模拟精度高等优点。该方法首先，对需要求解的模型区域进行网格剖分；然后，根据Taylor展开式使用有限个网格点上的波场值对波场的空间导数进行离散表示；最后，离散接收时间，并且在时间方向上逐一时刻递推求解整个模型区域内的波场，从而实现利用计算机通过有限差分法对地震波场的传播过程进行数值模拟。
有限差分法进行波场模拟时，通常采用空间高阶有限差分算子对波场的空间导数进行表示，使用远小于地震子波周期的时间步长对接收时间进行离散(通常接收时间被离散成几千甚至上万个时刻)，可以有效地提高波场模拟的精度，并且，防止数值频散现象。然而，在地震波场模拟的过程中，计算机需要使用高阶有限差分算子计算每个时刻模型中每一个网格点上的波场值，这使得波场模拟过程计算量和存储量很高，大模型中的波场模拟甚至无法直接在计算机内存中进行存储，三维模型情况下的波场模拟更是难以实现。
地球物理学者为了提高波场模拟的计算效率，通过低秩近似、辛网格以及并行计算等方法进行优化，但是几乎鲜有学者关注到波场传播过程的稀疏性。实际上，当模型中的震源激发后，有效的模拟区域仅为波前通过的模型区域，而波前未到达的模型区域不存在波场值，同样也就不需要进行波场模拟。
而射线类方法虽然模拟复杂构造模型地精度远不及波动方程法，但是计算效率远高于波动方程法。射线类方法中通过求解程函方程可以快速且较高精度地获得地震波的初至波走时，有效模拟每一时刻波场中波前的空间位置。若能将射线法与波动方程法相结合，是一个新的尝试，也是实现高效模拟地震波传播的有效途径。
发明内容
本发明的目的就在于，提供一种基于初至波走时约束的动态波场传播模拟方法，在波场模拟过程中，设计一种动态的空间窗口，对每一时刻的波场模拟范围进行约束，从而通过有限差分法对地震波传播进行动态模拟，同时将射线法与波动方程法相结合，以解决有效降低波场模拟的计算量，高效模拟地震波传播的问题。该方法不但兼顾波动方程法能有准确模拟复杂模型构造波场的优势，而且可以降低地震波传播模拟过程的计算量和内存消耗。
本发明的目的是通过以下技术方案实现的：
一种基于初至波走时约束的动态波场传播模拟方法，包括以下步骤：
A、将给定的速度模型进行网格化，设置模型设置横向网格数nx，纵向网格数nz，以及相应的网格间距Δx和Δz，得到离散速度模型v(x,z)，x∈[0,nx],z∈[0,nz]；
B、设置差分阶数M和差分系数通过公式(1)中的有限差分稳定性计算公式计算采样率Δt，设置最大时间采样时刻nt；
C、将离散的速度模型v(x,z)代入公式(2)中的程函方程，计算初至波走时T(x,z)，初至波走时是一系列等时线；
D、将初至波走时T(x,z)中每一条等时间线上的点对应的空间坐标x和z通过time-bin技术按时间增加方向重新排序，封装到每一时刻的动态结构体数组D[c(t)]中，其中，c(t)代表对应t时刻的等时线上的点数；
E、定义基于初至波走时约束的动态模拟范围Ωt：
F、对波动方程进行离散，构建波动方程高阶有限差分格式，如公式(5)，其中，代表在t+1时刻网格点(x,z)处的波场；
G、利用离散波动方程的高阶有限差分格式，逐一时刻在动态的空间窗口范围Ωt内进行波场模拟，直至最大时间采样时刻nt；
H、随着传播时刻增加，当波场传播至模型边界时，加载PML边界条件进行处理，消除由于有限的模型边界产生的反射。
进一步地，步骤E，定义基于初至波走时约束的动态模拟范围，具体包括以下步骤：对每个时刻对应的D[c(t)]中的空间坐标x和z分别进行比较，确定每一个时刻数值模拟的动态窗口边界。
更进一步地，对于任意时刻t，最小的x坐标赋值给模拟窗口的左边界最大的x坐标赋值给模拟窗口的右边界/>最小的z坐标赋值给模拟窗口的上边界/>最大的z坐标赋值给模拟窗口的下边界/>
与现有技术相比，本发明的有益效果是：
1、本发明模拟过程利用波场快照的稀疏性，设计一种动态的空间窗口对每一时刻的波场模拟范围进行约束，从而有效降低波场模拟的计算量和内存消耗；
2、本发明的创新之处在于所使用的动态窗口与波前的拟合程度，使用的动态窗口是通过初至走时计算获得的，完全拟合波场中传播最快的波前，有效避免因为窗口过小而造成有效波场截断产生反射，或者窗口过大无法有效降低计算量和内存消耗，可以更加准确的界定计算区域，准确地模拟波场。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1结合走时约束的动态波场数值模拟方法流程图；
图2Marmousi速度模型；
图3Marmousi模型的初至走时图；
图4初至走时time bin封装示意图；
图5空间约束的波场动态模拟示意图；
图6震源波场动态模拟示意图；
图7a动态波场模拟的Marmousi模型中2s时震源波场快照；
图7b常规波场模拟的Marmousi模型中2s时震源波场快照；
图8a动态波场模拟的Marmousi模型中3s时震源波场快照；
图8b常规波场模拟的Marmousi模型中3s时震源波场快照；
图9吸收边界条件示意图。
具体实施方式
下面结合实施例对本发明作进一步说明：
下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部结构。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
本发明提供一种基于初至波走时约束的动态波场传播模拟方法，将求解程函方程所获得初至波走时利用time-bin技术按时间方向重新排序，将模型网格点空间坐标封存到对应时刻的序列中，并且筛选出每一时刻的波前范围，从而设计出每一时刻最合适的有限差分模拟的动态空间窗口范围，对每一时刻的波场模拟范围进行约束，从而通过有限差分法对地震波传播进行动态模拟，可以有效降低波场模拟的计算量。值得注意的是，这个动态的空间窗口设定尤为重要，偏小的窗口可能会导致截断有效波场，产生反射，偏大的窗口又会造成计算浪费，无法有效地减少波场模拟的计算量。在模拟任意时刻的地震波场时，高阶有限差分算子不需要遍历整个模型空间，仅需要在动态空间窗口内进行模拟。
如图1所示，本发明基于初至波走时约束的动态波场传播模拟方法，包括以下步骤：
A、将给定的速度模型进行网格化，设置模型设置横向网格数nx，纵向网格数nz，以及相应的网格间距Δx和Δz，得到离散速度模型v(x,z)，x∈[0,nx],z∈[0,nz]；
B、设置差分阶数M和差分系数通过公式(1)中的有限差分稳定性计算公式计算采样率Δt，设置最大时间采样时刻nt；
C、将离散的速度模型v(x,z)代入公式(2)中的程函方程，计算初至波走时T(x,z)，初至波走时是一系列等时线；
D、将初至波走时T(x,z)中每一条等时间线上的点对应的空间坐标x和z通过time-bin技术按时间增加方向重新排序，封装到每一时刻的动态结构体数组D[c(t)]中，其中c(t)代表对应t时刻的等时线上的点数；
E、定义基于初至波走时约束的动态模拟范围Ωt，见公式(3)。具体过程为：对每个时刻对应的D[c(t)]中的空间坐标x和z分别进行比较，确定每一个时刻数值模拟的动态窗口边界。例如，对于任意时刻t，最小的x坐标赋值给模拟窗口的左边界最大的x坐标赋值给模拟窗口的右边界/>最小的z坐标赋值给模拟窗口的上边界/>最大的z坐标赋值给模拟窗口的下边界/>
相比于常规静态模拟范围Ω，公式(4)，Ωt能够更加准确的表述波前区域，从而减少对波前以外的无效冗余区域的波场模拟。
Ω＝(0,nx,0,nz) (4)
F、对波动方程进行离散，构建波动方程高阶有限差分格式，如公式(5)。其中代表在t+1时刻网格点(x,z)处的波场；
G、利用离散波动方程的高阶有限差分格式，逐一时刻在动态的空间窗口范围Ωt内进行波场模拟，直至最大时间采样时刻nt；
H、随着传播时刻增加，当波场传播至模型边界时，加载PML边界条件进行处理，消除由于有限的模型边界产生的反射。
实施例1
A、将图2中的Marmousi速度模型(7km*12km)进行网格化，设置模型设置横向网格数nx＝1201，纵向网格数nz＝701，以及相应的网格间距Δx＝10m和Δz＝10m，得到离散速度模型v(x,z)；
B、设置差分阶数M和差分系数这里差分阶数M＝5，差分系数通过公式(1)中的有限差分稳定性计算公式计算采样率Δt＝1ms，设置最大时间采样时刻nt＝4000；
C、将离散的速度模型v(x,z)代入公式(2)中的程函方程，计算初至波走时场T(x,z)，如图3所示；
D、将初至波走时T(x,z)中每一条等时间线上的点对应的空间坐标x和z通过time-bin技术按时间增加方向重新排序，封装到每一时刻的动态结构体数组D[c(t)]中，如图4所示，其中c(t)代表对应t时刻的等时线上的点数；
E、定义则基于初至波走时约束的动态模拟范围Ωt，见公式(3)。对每个时刻对应的D[c(t)]中的空间坐标x和z分别进行比较，确定每一个时刻数值模拟的动态窗口边界。例如，对于任意时刻t，最小的x坐标赋值给模拟窗口的左边界最大的x坐标赋值给模拟窗口的右边界/>最小的z坐标赋值给模拟窗口的上边界/>最大的z坐标赋值给模拟窗口的下边界/>
图5中的虚线即代表每一时刻的动态模拟范围，相比于常规静态模拟范围Ω，公式(4)，Ωt能够更加准确的表述波前区域，从而减少对波前以外的无效冗余区域的波场模拟。
Ω＝(0,nx,0,nz) (4)
F、对波动方程进行离散，构建波动方程高阶有限差分格式，如公式(5)。其中代表在t+1时刻网格点(x,z)处的波场；
G、利用离散波动方程的高阶有限差分格式，逐一时刻在动态的空间窗口范围Ωt内进行波场模拟，直至最大时间采样时刻nt＝4000；
相比于常规的有限差分波场模拟过程，基于初至波走时约束的动态波场模拟每个时刻仅计算波前通过的有效区域Ωt，有效降低计算量和内存占用，对比示意图如图6所示。抽取两种模拟方法在2s即nt＝2000和3s即nt＝3000时的波场快照进行对比，如图7和图8所示。对比可以发现，本发明中涉及的动态波场模拟方法和常规波场模拟方法在复杂模型中均可以达到相同的模拟结果，证明使用初至波走时设计的动态窗口不影响模拟精度。但完成整个模拟过程所需要的计算时间，动态波场模拟方法仅需要2min,而常规方法需要3min，并且内存消耗也仅为常规方法的65％。
H、当波场传播至模型边界时，加载PML边界进行处理，加载方式如图9所示，消除人工边界产生的反射。
注意，上述仅为本发明的较佳实施例及所运用技术原理。本领域技术人员会理解，本发明不限于这里所述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行了较为详细的说明，但是本发明不仅仅限于以上实施例，在不脱离本发明构思的情况下，还可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。
Claims (3)
1.一种基于初至波走时约束的动态波场传播模拟方法，其特征在于，包括以下步骤：
A、将给定的速度模型进行网格化，设置模型设置横向网格数nx，纵向网格数nz，以及相应的网格间距Δx和Δz，得到离散速度模型v(x,z)，x∈[0,nx],z∈[0,nz]；
B、设置差分阶数M和差分系数通过公式(1)中的有限差分稳定性计算公式计算采样率Δt，设置最大时间采样时刻nt；
C、将离散的速度模型v(x,z)代入公式(2)中的程函方程，计算初至波走时T(x,z)，初至波走时是一系列等时线；
D、将初至波走时T(x,z)中每一条等时间线上的点对应的空间坐标x和z通过time-bin技术按时间增加方向重新排序，封装到每一时刻的动态结构体数组D[c(t)]中，其中，c(t)代表对应t时刻的等时线上的点数；
E、定义基于初至波走时约束的动态模拟范围Ωt：
F、对波动方程进行离散，构建波动方程高阶有限差分格式，如公式(5)，其中，代表在t+1时刻网格点(x,z)处的波场；
G、利用离散波动方程的高阶有限差分格式，逐一时刻在动态的空间窗口范围Ωt内进行波场模拟，直至最大时间采样时刻nt；
H、随着传播时刻增加，当波场传播至模型边界时，加载PML边界条件进行处理，消除由于有限的模型边界产生的反射。
2.根据权利要求1所述的一种基于初至波走时约束的动态波场传播模拟方法，其特征在于：步骤E，定义基于初至波走时约束的动态模拟范围，具体包括以下步骤：对每个时刻对应的D[c(t)]中的空间坐标x和z分别进行比较，确定每一个时刻数值模拟的动态窗口边界。
3.根据权利要求2所述的一种基于初至波走时约束的动态波场传播模拟方法，其特征在于：对于任意时刻t，最小的x坐标赋值给模拟窗口的左边界最大的x坐标赋值给模拟窗口的右边界/>最小的z坐标赋值给模拟窗口的上边界/>最大的z坐标赋值给模拟窗口的下边界/>
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311644145.0A CN117892569B (zh) | 2023-12-04 | 2023-12-04 | 一种基于初至波走时约束的动态波场传播模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311644145.0A CN117892569B (zh) | 2023-12-04 | 2023-12-04 | 一种基于初至波走时约束的动态波场传播模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117892569A CN117892569A (zh) | 2024-04-16 |
| CN117892569B true CN117892569B (zh) | 2024-06-25 |
Family
ID=90639956
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311644145.0A Active CN117892569B (zh) | 2023-12-04 | 2023-12-04 | 一种基于初至波走时约束的动态波场传播模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117892569B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109164488A (zh) * | 2018-10-10 | 2019-01-08 | 西安交通大学 | 一种梯形网格有限差分地震波场模拟方法 |
| CN109725345A (zh) * | 2018-11-15 | 2019-05-07 | 中国石油天然气集团有限公司 | 一种初至波正演模拟方法及装置 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111948708B (zh) * | 2019-10-18 | 2021-09-28 | 中国石油大学(北京) | 一种浸入边界起伏地表地震波场正演模拟方法 |
| US11281825B2 (en) * | 2020-06-30 | 2022-03-22 | China Petroleum & Chemical Corporation | Computer-implemented method for high speed multi-source loading and retrieval of wavefields employing finite difference models |
| CN114895351A (zh) * | 2022-04-29 | 2022-08-12 | 中国科学院地质与地球物理研究所 | 模拟地震波在任意不连续界面处传播的介质模型化的方法及装置 |
-
2023
- 2023-12-04 CN CN202311644145.0A patent/CN117892569B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109164488A (zh) * | 2018-10-10 | 2019-01-08 | 西安交通大学 | 一种梯形网格有限差分地震波场模拟方法 |
| CN109725345A (zh) * | 2018-11-15 | 2019-05-07 | 中国石油天然气集团有限公司 | 一种初至波正演模拟方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117892569A (zh) | 2024-04-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US7675818B2 (en) | Method for tomographic inversion by matrix transformation | |
| US20160202375A1 (en) | Eikonal Solver for Quasi P-Waves in Anisotropic Media | |
| US20220350042A1 (en) | Method and system for super resolution least-squares reverse time migration | |
| CN110221342A (zh) | 基于三维速度模型的震源定位方法、装置及存储介质 | |
| CN116774292B (zh) | 一种地震波走时确定方法、系统、电子设备及存储介质 | |
| Alkhalifah et al. | An eikonal-based formulation for traveltime perturbation with respect to the source location | |
| Hu et al. | Ray-illumination compensation for adjoint-state first-arrival traveltime tomography | |
| Weber | Seismic traveltime tomography: a simulated annealing approach | |
| CN108680968B (zh) | 复杂构造区地震勘探数据采集观测系统评价方法及装置 | |
| Guo et al. | Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique | |
| CN114460646A (zh) | 一种基于波场激发近似的反射波旅行时反演方法 | |
| CN117892569B (zh) | 一种基于初至波走时约束的动态波场传播模拟方法 | |
| CN111665556A (zh) | 地层声波传播速度模型构建方法 | |
| CN111007565B (zh) | 三维频率域全声波成像方法及装置 | |
| Zhou et al. | An iterative factored topography-dependent eikonal solver for anisotropic media | |
| Bregman | Tomographic inversion of crosshole seismic data | |
| Gillberg et al. | Accuracy and efficiency of stencils for the eikonal equation in earth modelling | |
| Nowack | Wavefronts and solutions of the eikonal equation | |
| Michelini | Fault zone structure determined through the analysis of earthquake arrival times | |
| CN102262789B (zh) | 一种利用拉普拉斯方程的用于地球物理领域的随机点数据网格化方法 | |
| Franklin et al. | A high-order fast marching scheme for the linearized eikonal equation | |
| Shin et al. | Laplace-domain full waveform inversion using irregular finite elements for complex foothill environments | |
| AlMuhaidib et al. | Finite difference elastic wave modeling including surface topography | |
| Sava et al. | Huygens wavefront tracing: A robust alternative to ray tracing | |
| CN116166988B (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 |