NL2020684B1 - Mixed 2D seismic wave travel time calculation method - Google Patents
Mixed 2D seismic wave travel time calculation method Download PDFInfo
- Publication number
- NL2020684B1 NL2020684B1 NL2020684A NL2020684A NL2020684B1 NL 2020684 B1 NL2020684 B1 NL 2020684B1 NL 2020684 A NL2020684 A NL 2020684A NL 2020684 A NL2020684 A NL 2020684A NL 2020684 B1 NL2020684 B1 NL 2020684B1
- Authority
- NL
- Netherlands
- Prior art keywords
- travel time
- node
- grid
- nodes
- seismic
- Prior art date
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 43
- 238000000034 method Methods 0.000 claims abstract description 37
- 238000010276 construction Methods 0.000 claims abstract description 13
- 230000014509 gene expression Effects 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 abstract description 3
- 241000826860 Trapezium Species 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
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/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
- G01V1/305—Travel times
-
- 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
- G01V2210/671—Raytracing
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)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a mixed 2D seismic travel time calculation method comprising inputting relevant parameters and a velocity model; emitting rays in different directions and calculating the information of central rays; calculating the seismic travel time of grid nodes within the ray range by using the wave front construction method; classifying the travel time attributes of grid nodes in the computational area and establishing an initial narrow band; and calculating the seismic travel time of the rest grid nodes by fast marching method. The fast marching method is connected with the wave front construction method through narrow band technology. Seismic travel time calculation precision of a small area near the source is improved by the wave-front construction method. Calculation precision of grid nodes in the rest area is increased by the fast marching method. A 2D seismic travel time calculation method with calculation efficiency and calculation precision is realized.
Description
Mixed 2D seismic travel time calculation method
Technical Field
The present invention relates to the field of seismic travel time calculation, in particular to a mixed 2D seismic travel time calculation method.
Background Art
No. 3, 2009 of Chinese Journal of Engineering Geophysics published “Analysis and Improvement on Calculation Accuracy of Fast Marching Method” by Zhang Shuangjie et al, which introduces two methods for improving the calculation efficiency of the fast marching method: firstly, performing calculating by adopting a high order difference scheme; secondly, performing local precise treatment on grid nodes near the source, so that the high order difference scheme is used for calculation near the source and a low order difference scheme is used for calculation in the rest area. The two methods are also used for calculating the seismic travel time of the homogeneous models, and beneficial experimental results are obtained.
No. 3, 2016 of Global Geology published “Calculation of Fast Marching Method on Seismic Travel-Time Based on Complete Ternary Tree” by Wang Qianlong et al, which introduces an improved fast marching method for calculating seismic travel time. By introducing a complete ternary tree heap sorting method into seismic travel time calculation, the improved method reduces the time to find the minimum travel time in narrow band extension and improves the calculation efficiency of the whole calculation method. Besides, the fast marching method for calculating seismic travel time based on a complete ternary tree is used to calculate the layer model, Marmousi model and Sigsbee 2a models, and beneficial experimental results are obtained.
As seen from the above examples, although through the fast marching method for calculating seismic travel time, the calculation accuracy can be improved to a certain extent, the realization process is complicated and improvement on the calculation accuracy is also limited.
Summary
The present invention aims to solve the technical problem of providing a mixed 2D seismic travel time calculation method by which a wave front construction method for calculating seismic travel time is seamlessly connected with a fast marching method for calculating seismic travel time together by flexibly using the narrow band technology, namely that the wave front construction method with high calculation accuracy is used in a small area near the source, and the fast marching method is used in the rest area to calculate the travel time. While the calculation accuracy of the new method is improved, the high efficiency is still remained.
In order to solve the technical problem, the present invention adopts the technical scheme: 1. a mixed 2D seismic travel time calculation method, comprising the following steps of: step 1: inputting a relevant parameter document and a velocity model, wherein the parameter document contains the grid node number, the grid space and the source position of the velocity model; step 2: emitting rays in different directions from the source and calculating the information of central rays; step 3: calculating the travel time of grid nodes within the ray range by a wave front construction method; step 4: classifying all the grid nodes and dividing the grid nodes into accepted nodes, narrow-band nodes or far-away nodes, namely that if the seismic travel time of one grid node is calculated and all the travel time of surrounding grid nodes is calculated already, then the node is a accepted node, if the seismic travel time of one grid node is calculated but not all the travel time of surrounding grid nodes is calculated, then the node is a narrow-band node, and if the travel time of one node is not calculated yet, the node is the far-away node; step 5: inserting all narrow-band nodes into the narrow band, and establishing an initial narrow band; and step 6: extending the narrow band until the narrow band is empty, wherein in the process, the travel time of the grid nodes is obtained by solving a 2D eikonal equation by the upwind difference method, and the 2D eikonal equation is as shown below:
wherein τ is the seismic travel time, 5 is slowness, V is gradient, and the upwind difference expression of the gradient term in the formula is as shown below:
wherein D~'t , D+*t , ZTjr and D^t are forward or backward difference expressions of travel time at the grid node (z, yj in an x or z direction respectively, and the concrete forms are
respectively.
Further, in step 2, the kinematical ray tracing equations are solved by a Runge-Kutta method to get the information of central rays, as shown below:
wherein x, represents spatial position, pt represents slowness, τ represents seismic travel time, and v represents the velocity value at a discrete node.
Compared with the prior art, the mixed 2D seismic travel time calculation method disclosed by the present invention has the beneficial effects that: because the seismic travel time calculation precision of grid nodes near the source is improved by the wave front construction method, the accuracy of the fast marching method for calculating grid nodes in the area is improved; and the calculation precision of the whole seismic travel time calculation method is greatly improved at the cost of reducing the calculation efficiency a little.
Brief Description of the Drawings
Fig. 1 shows the flow chart of the mixed 2D seismic travel time calculation method of the present invention.
Fig. 2 shows the schematic diagram of area division by the wave front construction method. Fig. 3 shows the relative error in a homogeneous medium by using the fast marching method. Fig. 4 shows the relative error in a homogeneous medium by using mixed seismic travel time calculation method.
Detailed Description
The invention will now be further described in detains in connection with the accompanying drawings and embodiments.
Fig. 1 shows the flow chart of the mixed 2D seismic travel time calculation method of the present invention. In the method, a reaction flow is shown, as shown below: step 1: inputting a relevant parameter document and a velocity model, wherein the parameter document contains the grid node number, the grid space and the source position; 2)emitting rays in different directions from a shot point and calculating the information of central rays, wherein the emission angle range of the rays is from -90 degrees to +90 degrees, the angle interval between every two adjacent rays is from 3 degrees to 10 degrees, and the kinematical ray tracing equations are solved by a Runge-Kutta method to get the information of central rays, as shown below:
wherein y represents spatial position, represents slowness, τ represents seismic travel time, and v represents the velocity value at a discrete node. step 3: calculating the travel time of grid nodes within the ray range by a wave front construction method, wherein the area is divided into multiple trapeziums by adjacent wave front points on adjacent central rays in the calculation process, and the grid nodes in each trapezium are acquired by interpolating vertex information of the corresponding trapezium; step 4: after the wave front construction method is completed, performing attributive classification on all of the grid nodes and specifying that if the seismic travel time of one grid node is calculated and all the travel time of surrounding grid nodes is calculated already, then the node is a accepted node, if the seismic travel time of one grid node is calculated but not all the travel time of surrounding grid nodes is calculated, then the node is a narrow-band node, and if the travel time of one node is not calculated yet, the node is the far-away node; step 5: moving all narrow-band nodes into a narrow band, and establishing an initial narrow band; and step 6: extending the narrow band until the narrow band is empty, wherein in the process, the
travel time of the grid nodes is obtained by solving a 2D eikonal equation by the upwind difference method, and the 2D eikonal equation is as shown below:
wherein τ is the seismic travel time, s is slowness, V is gradient, and the upwind difference expression of the gradient term in the formula is as shown below:
wherein
are forward or backward difference expressions of travel time at the grid node (z, j) in an x or z direction respectively, and the concrete forms are
respectively.
The calculation accuracy and stability of the method in disclosed by the present invention will be analyzed and verified below through a homogeneous medium model.
Fig. 3 and Fig. 4 show the relative errors in the homogeneous medium model respectively by the fast marching method and the mixed seismic travel time calculation method, the size of the homogeneous model is 601*601, the grid space is 10m*10m and the velocity is 1,000 m/s. The maximum length of a single tracked ray in the joint travel time method is 500m. As can be seen from the figure, the calculation precision of the mixed seismic travel time method is greatly improved than that of the fast marching method.
According to the mixed seismic travel time calculation method disclosed by the present invention, the fast marching method is connected with the wave front construction method through a narrow band technology, and the seismic travel time calculation accuracy in a small area near the source is improved by the wave front construction method, so that the accuracy of the fast marching method for calculating grid nodes in the rest area is improved, and a 2D seismic travel time calculation method considering the calculation efficiency and the calculation accuracy is realized.
Claims (2)
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810077621.8A CN108072897A (en) | 2018-01-23 | 2018-01-23 | It is a kind of to mix computational methods when two-dimensionally seismic wave is walked |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| NL2020684A NL2020684A (en) | 2018-04-20 |
| NL2020684B1 true NL2020684B1 (en) | 2018-11-14 |
Family
ID=61978378
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| NL2020684A NL2020684B1 (en) | 2018-01-23 | 2018-03-29 | Mixed 2D seismic wave travel time calculation method |
Country Status (4)
| Country | Link |
|---|---|
| CN (1) | CN108072897A (en) |
| BE (1) | BE1025488B1 (en) |
| LU (1) | LU100751B1 (en) |
| NL (1) | NL2020684B1 (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108957538A (en) * | 2018-06-21 | 2018-12-07 | 成都启泰智联信息科技有限公司 | A kind of virtual focus two dimension wavefront construction seimic travel time calculation method |
| CN115201901A (en) * | 2022-06-30 | 2022-10-18 | 中铁第四勘察设计院集团有限公司 | Method, device and equipment for determining tunnel wave front travel time and readable storage medium |
| CN117194855B (en) * | 2023-11-06 | 2024-03-19 | 南方科技大学 | Fitting analysis method and relevant equipment for weak anisotropy travel time |
| CN117607957B (en) * | 2024-01-24 | 2024-04-02 | 南方科技大学 | Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5724310A (en) * | 1996-10-15 | 1998-03-03 | Western Atlas International, Inc. | Traveltime generation in the presence of velocity discontinuities |
| US6018499A (en) * | 1997-11-04 | 2000-01-25 | 3Dgeo Development, Inc. | Three-dimensional seismic imaging of complex velocity structures |
| WO2002065159A1 (en) * | 2001-02-14 | 2002-08-22 | Hae-Ryong Lim | Method of seismic imaging using direct travel time computing |
| US8094513B2 (en) * | 2008-06-03 | 2012-01-10 | Westerngeco L.L.C. | Determining positioning of survey equipment using a model |
| CN106353793A (en) * | 2015-07-17 | 2017-01-25 | 中国石油化工股份有限公司 | Cross-well seismic tomography inversion method on basis of travel time incremental bilinear interpolation ray tracing |
| CN105425286A (en) * | 2015-10-30 | 2016-03-23 | 中国石油天然气集团公司 | Earthquake time-travelling acquisition method and crosshole earthquake time-travelling tomography method based on the earthquake time-travelling acquisition method |
-
2018
- 2018-01-23 CN CN201810077621.8A patent/CN108072897A/en active Pending
- 2018-03-29 NL NL2020684A patent/NL2020684B1/en not_active IP Right Cessation
- 2018-03-30 LU LU100751A patent/LU100751B1/en active IP Right Grant
- 2018-03-30 BE BE2018/0040A patent/BE1025488B1/en not_active IP Right Cessation
Also Published As
| Publication number | Publication date |
|---|---|
| BE1025488A1 (en) | 2019-03-18 |
| LU100751B1 (en) | 2018-10-01 |
| CN108072897A (en) | 2018-05-25 |
| NL2020684A (en) | 2018-04-20 |
| BE1025488B1 (en) | 2019-03-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| NL2020684B1 (en) | Mixed 2D seismic wave travel time calculation method | |
| US10062207B2 (en) | Method and system for reconstructing a three-dimensional model of point clouds | |
| CN104933064B (en) | The method and apparatus for predicting the kinematic parameter of destination object | |
| CN102692644B (en) | Depth domain common-image gather generation method | |
| NL2021354A (en) | 2-D Seismic Travel Time Calculation Method Based on Virtual Source Wavefront Construction | |
| Noichl et al. | " BIM-to-Scan" for Scan-to-BIM: Generating Realistic Synthetic Ground Truth Point Clouds based on Industrial 3D Models | |
| Gebreiter et al. | sbpRAY–A fast and versatile tool for the simulation of large scale CSP plants | |
| CN109255837A (en) | A kind of building method of the efficient B-spline surface for laser radar point cloud data processing | |
| EP2917768A1 (en) | System and method for analysis of designs of a seismic survey | |
| CN104502967B (en) | The method and device of quick obtaining seismic prospecting observation system bin information | |
| Quan et al. | Filtering LiDAR data based on adjacent triangle of triangulated irregular network | |
| CN108680968A (en) | Complex structural area seismic prospecting data collecting observation system evaluation method and device | |
| CN111076724B (en) | Three-dimensional laser positioning method and system | |
| AU739128B2 (en) | A method of seismic processing, and in particular a 3D seismic prospection method implementing seismic data migration | |
| CN110687598A (en) | Method and device for accelerating microseismic numerical simulation | |
| CN116381650A (en) | Laser radar point cloud position and intensity simulation and test method | |
| CN108051854B (en) | A kind of multiple dimensioned pseudo- bending ray tracing method | |
| US5265068A (en) | Efficient generation of traveltime tables for complex velocity models for depth migration | |
| Nunes et al. | Procedural generation of synthetic forest environments to train machine learning algorithms | |
| Schöler et al. | Towards an automated 3D reconstruction of plant architecture | |
| CN106556860B (en) | The method and apparatus for laying VSP observation systems | |
| CN112308134A (en) | Static fusion method based on Gaussian mixture probability hypothesis density filter | |
| Stefanelli et al. | Synthetic point cloud generation for class segmentation applications | |
| CN103020999A (en) | Method and system for identifying deformation area | |
| Tóth et al. | Tree growth simulation based on ray-traced lights modelling |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| MM | Lapsed because of non-payment of the annual fee |
Effective date: 20210401 |