CN114879249B - 基于四面体单元走时扰动插值的地震波前走时计算方法 - Google Patents
基于四面体单元走时扰动插值的地震波前走时计算方法 Download PDFInfo
- Publication number
- CN114879249B CN114879249B CN202210384391.6A CN202210384391A CN114879249B CN 114879249 B CN114879249 B CN 114879249B CN 202210384391 A CN202210384391 A CN 202210384391A CN 114879249 B CN114879249 B CN 114879249B
- Authority
- CN
- China
- Prior art keywords
- travel time
- node
- nodes
- representing
- seismic
- 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
- 238000004364 calculation method Methods 0.000 title claims abstract description 33
- 230000001788 irregular Effects 0.000 claims abstract description 4
- 238000000034 method Methods 0.000 abstract description 16
- 230000007547 defect Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
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
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
技术领域
本发明涉及一种基于四面体单元走时扰动插值的地震波前走时计算方法，属于地震波前走时正演模拟技术领域。
背景技术
地震波前走时计算是地震射线追踪、层析成像、速度反演、震源定位和偏移成像等的重要环节。线性走时插值法是计算地震波前走时的常用方法之一。该方法假设波前走时在单元边界上呈线性变化，单元边界上任意点的波前走时可由已知节点走时的线性插值函数表示。然而，这种线性变化的假设与实际的非线性变化相差较大，会引起较大的走时计算误差。Zhang等（2015）提出了基于规则矩形单元的线性走时扰动插值法，该方法将实际走时分解为参考走时和走时扰动。参考走时为波在等效均匀介质中的走时，沿着单元边界非线性变化；走时扰动定义为实际走时和参考走时之差，并假设沿单元边界线性变化。由于走时扰动远小于参考走时，使得实际走时沿单元边界保持非线性变化，克服了线性走时插值法的缺陷，提高了波前走时的计算精度。但模型被离散为规则矩形单元，难以精确模拟起伏地形、起伏界面和倾斜断裂等复杂地质构造。
发明内容
针对现有技术存在的上述缺陷，本发明提出了一种基于四面体单元走时扰动插值的地震波前走时计算方法，其克服了基于规则矩形单元线性走时扰动插值法难以精确模拟起伏地形的问题。
本发明是采用以下的技术方案实现的：本发明所述一种基于四面体单元走时扰动插值的地震波前走时计算方法，包括如下步骤：
步骤一：采用不规则的四面体单元对速度模型进行离散化；
步骤二：给定初始条件；将距离震源大于100m的网格节点的走时设置为100000；计算距离震源小于和等于100m的网格节点的走时，计算公式如下：
（1）
其中，
t表示地震波从震源到达节点的走时，
s表示震源与节点的距离，
v表示震源点的速度；
步骤三：寻找一个待计算走时的节点
P；
步骤四：采用三个节点走时信息计算步骤三寻找的节点
P的走时；采用节点
P相邻节点中三个走时最小的节点信息计算节点
P的走时，计算公式如下：
（2）
其中，
t p 表示地震波从震源到达节点
P的走时，Δ
t K 表示震源到达节点
P的路径与四面体交点
K处的走时扰动量，
s eq 表示震源与交点
K间的等效慢度，
l SK 表示震源与交点
K间的距离，
s表示四面体单元的慢度，
l KP 表示交点
K与节点
P间的距离；
步骤五：采用两个节点走时信息计算步骤三寻找的节点
P的走时；采用的两个节点是节点
P相邻节点中三个走时最小的节点中的任意两个节点；
步骤六：计算地震波沿四面体的棱传播到达节点
P的走时；地震波经过节点
P1所在的棱到达节点
P的走时，采用下式计算：
（3）
其中，
t 1表示地震波先到达节点
P1再经过四面体的棱到达节点
P的走时，
t P1表示地震波从震源到达节点
P相邻三个走时最小的节点之一
P1的走时，
s表示四面体单元的慢度，
l P1
P 表示节点
P1和节点
P间的距离；地震波经过节点
P2所在的棱到达节点
P的走时，采用下式计算：
（4）
其中，
t 2表示地震波先到达节点
P2再经过四面体的棱到达节点
P的走时，
t P2表示地震波从震源到达节点
P相邻三个走时最小的节点之一
P2的走时，
s表示四面体单元的慢度，
l P2
P 表示节点
P2和节点
P间的距离；地震波经过节点
P3所在的棱到达节点
P的走时，采用下式计算：
（5）
其中，
t 3表示地震波先到达节点
P3再经过四面体的棱到达节点
P的走时，
t P3表示地震波从震源到达节点
P相邻三个走时最小的节点之一
P3的走时，
s表示四面体单元的慢度，
l P3
P 表示节点
P3和节点
P间的距离；
步骤七：分别判断步骤四、五和六计算的节点
P的走时是否满足因果性；如果满足，则保留，否则，舍弃；
步骤八：从步骤七保留的走时中选择最小值，作为节点
P最终的走时；
步骤九：判断是否完成所有节点的计算；如果完成，输出所有节点的走时，否则，返回步骤三。
进一步地，所述步骤四中，等效慢度
s eq 采用下式进行计算：
（6）
其中，
t P1、
t P2和
t P3分别表示地震波从震源到达节点
P相邻三个走时最小的节点
P1、
P2和
P3的走时，
l SP1、
l SP2和
l SP3分别表示震源与节点
P1、
P2和
P3间的距离。
进一步地，所述步骤四中，走时扰动量Δ
t K 采用下式进行计算：
（7）
其中，Δ
y和Δ
z分别表示交点
K与节点
P1、
P2和
P3的中点在
y和
z方向的距离，a、b和c是三个与走时扰动量相关的系数，通过节点
P1、
P2和
P3的走时扰动量确定。
本发明的有益效果是：采用本发明所述的基于四面体单元走时扰动插值的地震波前走时计算方法，通过对速度模型进行四面体单元离散，可以精确模拟起伏地形，其克服了矩形规则单元离散速度模型时难以准确模拟起伏地形的问题；本发明计算简单、易于实现，具有很强的适应性，且计算结果精度高。
附图说明
图1为本发明的流程图；
图2为采用规则网格离散化计算的Y方向1.5km处垂直切片的波前走时等值线图；
图3为本发明计算的Y方向1.5km处垂直切片的波前走时等值线图。
具体实施方式
为了使本发明目的、技术方案更加清楚明白，下面结合附图，对本发明作进一步详细说明。
为使本发明的上述和其他目的、特征和优点更明显易懂，下文结合实例，并配合所附图，作详细说明。本发明的流程图，如图1所示，包括如下步骤：
步骤一：采用不规则的四面体单元对速度模型进行离散化；
步骤二：给定初始条件；将距离震源大于100m的网格节点的走时设置为100000；计算距离震源小于和等于100m的网格节点的走时，计算公式如下：
（1）
其中，
t表示地震波从震源到达节点的走时，
s表示震源与节点的距离，
v表示震源点的速度；
步骤三：寻找一个待计算走时的节点
P；
步骤四：采用三个节点走时信息计算步骤三寻找的节点
P的走时；采用节点
P相邻节点中三个走时最小的节点信息计算节点
P的走时，计算公式如下：
（2）
其中，
t p 表示地震波从震源到达节点
P的走时，Δ
t K 表示震源到达节点
P的路径与四面体交点
K处的走时扰动量，
s eq 表示震源与交点
K间的等效慢度，
l SK 表示震源与交点
K间的距离，
s表示四面体单元的慢度，
l KP 表示交点
K与节点
P间的距离；
步骤五：采用两个节点走时信息计算步骤三寻找的节点
P的走时；采用的两个节点是节点
P相邻节点中三个走时最小的节点中的两个；
步骤六：计算地震波沿四面体的棱传播到达节点
P的走时；地震波经过节点
P1所在的棱到达节点
P的走时，采用下式计算：
（3）
其中，
t 1表示地震波先到达节点
P1再经过四面体的棱到达节点
P的走时，
t P1表示地震波从震源到达节点
P相邻三个走时最小的节点之一
P1的走时，
s表示四面体单元的慢度，
l P1
P 表示节点
P1和节点
P间的距离；地震波经过节点
P2所在的棱到达节点
P的走时，采用下式计算：
（4）
其中，
t 2表示地震波先到达节点
P2再经过四面体的棱到达节点
P的走时，
t P2表示地震波从震源到达节点
P相邻三个走时最小的节点之一
P2的走时，
s表示四面体单元的慢度，
l P2
P 表示节点
P2和节点
P间的距离；地震波经过节点
P3所在的棱到达节点
P的走时，采用下式计算：
（5）
其中，
t 3表示地震波先到达节点
P3再经过四面体的棱到达节点
P的走时，
t P3表示地震波从震源到达节点
P相邻三个走时最小的节点之一
P3的走时，
s表示四面体单元的慢度，
l P3
P 表示节点
P3和节点
P间的距离；
步骤七：分别判断步骤四、五和六计算的节点
P的走时是否满足因果性；如果满足，则保留，否则，舍弃；
步骤八：从步骤七保留的走时中选择最小值，作为节点
P最终的走时；
步骤九：判断是否完成所有节点的计算；如果完成，输出所有节点的走时，否则，返回步骤三。
其中，所述步骤四中，等效慢度
s eq 采用下式进行计算：
（6）
式中，
t P1、
t P2和
t P3分别表示地震波从震源到达节点
P相邻三个走时最小的节点
P1、
P2和
P3的走时，
l SP1、
l SP2和
l SP3分别表示震源与节点
P1、
P2和
P3间的距离。
其中，所述步骤四中，走时扰动量Δ
t K 采用下式进行计算：
（7）
式中，Δ
y和Δ
z分别表示交点
K与节点
P1、
P2和
P3的中点在
y和
z方向的距离，a、b和c是三个与走时扰动量相关的系数，通过节点
P1、
P2和
P3的走时扰动量确定。
实施例一：
下面结合具体实施方式，根据一个三层均匀层状速度模型进行解释和说明。
为了进一步说明本方法的实现思路及实现过程并证明方法的有效性，用一个三层均匀层状速度模型进行测试，并和常规方法计算结果进行比较。
S1:建立一个三层均匀层状速度模型。该模型地表起伏变化，大小为3000m×3000m×3000m，速度从上到下依次是2000m/s、3000m/s和4000m/s，震源点位于（1500m,1500m,84m）。
S2:分别采用边长约为15m和60m的四面体单元离散速度模型。
S3:采用公式（1）计算距离震源小于和等于100m的网格节点的走时。
S4:根据群扩展算法寻找一个待计算走时的节点
P。
S5:采用公式（2），根据节点
P相邻节点中三个走时最小的节点信息计算节点
P的走时。
S6:采用节点
P相邻节点中三个走时最小的节点中的两个节点信息计算节点
P的走时。
S7:采用公式（3）、（4）和（5），计算地震波沿四面体的棱传播到达节点
P的走时。
S8:分别判断步骤五、六和七计算的节点
P的走时是否满足因果性；如果满足，则保留，否则，舍弃。
S9:从步骤八保留的走时中选择最小值，作为节点
P最终的走时。
S10:判断是否完成所有节点的计算；如果完成，输出所有节点的走时，否则，返回步骤四。图3是根据上述步骤计算的Y方向1.5km处垂直切片的波前走时等值线图。
图2是采用规则网格离散化计算的Y方向1.5km处垂直切片的波前走时等值线图。从图2中可以明显看出，采用边长15m和60m的规则矩形单元离散速度模型计算的波前走时等值线不重合，表明规则网格离散速度模型会导致计算的波前走时存在明显的误差。然而，本发明采用不同大小四面体单元离散速度模型后计算的波前走时是重合的（详见图3），表明本发明计算的波前走时更加准确。
以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。
Claims (1)
1.一种基于四面体单元走时扰动插值的地震波前走时计算方法，其特征在于，包括如下步骤：
步骤一：采用不规则的四面体单元对速度模型进行离散化；
步骤二：给定初始条件；将距离震源大于100m的网格节点的走时设置为100000；计算距离震源小于和等于100m的网格节点的走时，计算公式如下：
其中，t表示地震波从震源到达节点的走时，l表示震源与节点间的距离，v表示震源点的速度；
步骤三：寻找一个待计算走时的节点P；
步骤四：采用三个节点走时信息计算步骤三寻找的节点P的走时；采用节点P相邻节点中三个走时最小的节点信息计算节点P的走时，计算公式如下：
tP＝ΔtK+seqlSK+slKP (2)
其中，tp表示地震波从震源到达节点P的走时，ΔtK表示震源到达节点P的路径与四面体交点K处的走时扰动量，seq表示震源与交点K间的等效慢度，lSK表示震源与交点K间的距离，s表示四面体单元的慢度，lKP表示交点K与节点P间的距离；走时扰动量ΔtK采用下式进行计算：
ΔtK＝aΔy+bΔz+c (3)
其中，Δy和Δz分别表示交点K与节点P相邻节点中三个走时最小的节点P1、P2和P3的中点在y和z方向的距离，a、b和c分别是与走时扰动量相关的三个系数，通过节点P1、P2和P3的走时扰动量确定；
步骤五：采用两个节点走时信息计算步骤三寻找的节点P的走时；采用的两个节点是节点P相邻节点中三个走时最小的节点中的两个；
步骤六：计算地震波沿四面体的棱传播到达节点P的走时；地震波经过节点P1所在的棱到达节点P的走时，采用下式计算：
t1＝tP1+slP1P (4)
其中，t1表示地震波先到达节点P1再经过四面体的棱到达节点P的走时，
表示地震波从震源到达节点P1的走时，s表示四面体单元的慢度，lP1P表示节点P1和节点P间的距离；地震波经过节点P2所在的棱到达节点P的走时，采用下式计算：
t2＝tP2+slP２P (5)
其中，t2表示地震波先到达节点P2再经过四面体的棱到达节点P的走时，
表示地震波从震源到达节点P2的走时，s表示四面体单元的慢度，lP2P表示节点P2和节点P间的距离；地震波经过节点P3所在的棱到达节点P的走时，采用下式计算：
t3＝tP3+slP3P (6)
其中，t3表示地震波先到达节点P3再经过四面体的棱到达节点P的走时，
表示地震波从震源到达节点P3的走时，s表示四面体单元的慢度，lP3P表示节点P3和节点P间的距离；
步骤七：分别判断步骤四、五和六计算的节点P的走时是否满足因果性；如果满足，则保留，否则，舍弃；
步骤八：从步骤七保留的走时中选择最小值，作为节点P最终的走时；
步骤九：判断是否完成所有节点的计算；如果完成，输出所有节点的走时，否则，返回步骤三。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210384391.6A CN114879249B (zh) | 2022-04-13 | 2022-04-13 | 基于四面体单元走时扰动插值的地震波前走时计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210384391.6A CN114879249B (zh) | 2022-04-13 | 2022-04-13 | 基于四面体单元走时扰动插值的地震波前走时计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114879249A CN114879249A (zh) | 2022-08-09 |
| CN114879249B true CN114879249B (zh) | 2023-04-28 |
Family
ID=82669498
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210384391.6A Active CN114879249B (zh) | 2022-04-13 | 2022-04-13 | 基于四面体单元走时扰动插值的地震波前走时计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114879249B (zh) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107843922A (zh) * | 2017-12-25 | 2018-03-27 | 中国海洋大学 | 一种基于地震初至波和反射波走时联合的层析成像方法 |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101533102B (zh) * | 2009-04-09 | 2011-07-20 | 长江工程地球物理勘测武汉有限公司 | 二维复杂结构三角网射线追踪全局方法 |
| AU2011320352B2 (en) * | 2010-10-29 | 2015-04-23 | Schlumberger Technology B.V. | Model based inversion of seismic response for determining formation properties |
| CN106353793A (zh) * | 2015-07-17 | 2017-01-25 | 中国石油化工股份有限公司 | 一种基于走时增量双线性插值射线追踪的井间地震层析反演方法 |
| CN106772577B (zh) * | 2016-06-29 | 2019-04-26 | 中国石油大学(华东) | 基于微地震数据和spsa优化算法的震源反演方法 |
| BR112019022031A2 (pt) * | 2017-04-21 | 2020-05-12 | Repsol Exploración, S.A. | Método de migração inversa de profundidade de kirchhoff pós-empilhamento para meios isotrópicos transversais inclinados (tti) e heterogêneos com base em rastreamento de raios em dados migrados |
| WO2019071504A1 (zh) * | 2017-10-12 | 2019-04-18 | 南方科技大学 | 一种基于两点射线追踪的地震走时层析反演方法 |
| CN109444955B (zh) * | 2019-01-09 | 2020-05-19 | 中国海洋大学 | 三维地震射线追踪的双线性走时扰动插值方法 |
| CN112099082B (zh) * | 2019-06-17 | 2021-08-20 | 中国海洋大学 | 一种共面元共方位角道集的地震回折波走时反演方法 |
-
2022
- 2022-04-13 CN CN202210384391.6A patent/CN114879249B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107843922A (zh) * | 2017-12-25 | 2018-03-27 | 中国海洋大学 | 一种基于地震初至波和反射波走时联合的层析成像方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114879249A (zh) | 2022-08-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US5812493A (en) | Method of representing the trajectory of at least one borehole in a space-time domain | |
| CN111948708B (zh) | 一种浸入边界起伏地表地震波场正演模拟方法 | |
| CN113671570B (zh) | 一种地震面波走时和重力异常联合反演方法与系统 | |
| CN112254798A (zh) | 一种预报海洋矢量声场的方法、系统及介质 | |
| CN114139335A (zh) | 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 | |
| Maleika | The influence of track configuration and multibeam echosounder parameters on the accuracy of seabed DTMs obtained in shallow water | |
| CN111158059A (zh) | 基于三次b样条函数的重力反演方法 | |
| CN111538079A (zh) | 基于全波形反演技术确定地质裂缝柔度参数的方法及装置 | |
| CN103698810A (zh) | 混合网最小走时射线追踪层析成像方法 | |
| CN108680968B (zh) | 复杂构造区地震勘探数据采集观测系统评价方法及装置 | |
| CN114879249B (zh) | 基于四面体单元走时扰动插值的地震波前走时计算方法 | |
| CN111665556A (zh) | 地层声波传播速度模型构建方法 | |
| Demin et al. | A numerical method of calculation of currents and sea surface topography in multiply connected domains of the ocean | |
| CN106353801A (zh) | 三维Laplace域声波方程数值模拟方法及装置 | |
| CN106950598B (zh) | 一种偏移速度场可靠性评价方法 | |
| CN107807392A (zh) | 一种自适应抗频散的分块时空双变逆时偏移方法 | |
| CN112257241B (zh) | 一种三角网菲涅尔带时差层析反演方法 | |
| CN110568497B (zh) | 一种复杂介质条件下地震初至波旅行时的精确求解方法 | |
| CN112649859B (zh) | 一种地震波速度自适应无网格场节点建立方法及系统 | |
| CN112596103A (zh) | 射线追踪方法、装置和电子设备 | |
| CN103425881A (zh) | 一种裂缝介质地震波响应的确定性数值模拟方法 | |
| Sava et al. | Huygens wavefront tracing: A robust alternative to ray tracing | |
| CN112346115A (zh) | 地下硐室群含空洞、复杂岩体波速环境下的微震震源定位方法 | |
| CN116660974B (zh) | 一种基于结构耦合约束的体波和面波三维联合反演方法 | |
| RU2611892C1 (ru) | Способ трехмерного моделирования заданного гидрогеологического объекта, реализуемый в вычислительной системе |
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 |