CN108072897A - It is a kind of to mix computational methods when two-dimensionally seismic wave is walked - Google Patents

It is a kind of to mix computational methods when two-dimensionally seismic wave is walked Download PDF

Info

Publication number
CN108072897A
CN108072897A CN201810077621.8A CN201810077621A CN108072897A CN 108072897 A CN108072897 A CN 108072897A CN 201810077621 A CN201810077621 A CN 201810077621A CN 108072897 A CN108072897 A CN 108072897A
Authority
CN
China
Prior art keywords
mrow
point
msup
walked
mfrac
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
CN201810077621.8A
Other languages
Chinese (zh)
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 CN201810077621.8A priority Critical patent/CN108072897A/en
Priority to NL2020684A priority patent/NL2020684B1/en
Priority to BE2018/0040A priority patent/BE1025488B1/en
Priority to LU100751A priority patent/LU100751B1/en
Publication of CN108072897A publication Critical patent/CN108072897A/en
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30Analysis
    • G01V1/303Analysis for determining velocity profiles or travel times
    • G01V1/305Travel times
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/67Wave propagation modeling
    • G01V2210/671Raytracing

Abstract

Computational methods when two-dimensionally seismic wave is walked the invention discloses a kind of mixing, comprise the following steps:Read in relevant parameter, rate pattern;Ray is sent to different directions and calculate central ray information along focus;The seimic travel time of radiation range internal net point is calculated using wavefront construction method;Attribute is classified and establishes initial narrow band on this basis when being walked to total space mesh point;Remaining mesh point seimic travel time is calculated using fast-marching algorithm.By the present invention in that it is connected to fast-marching algorithm and wavefront construction method with narrow-band technologies, the seimic travel time computational accuracy of smaller area near focus is improved by wavefront construction method, so as to improving computational accuracy of the fast-marching algorithm in remaining area grid node, realize when a kind of two-dimensionally seismic wave for taking into account computational efficiency and computational accuracy is walked and calculate.

Description

It is a kind of to mix computational methods when two-dimensionally seismic wave is walked
Technical field
It is particularly a kind of to mix computational methods when two-dimensionally seismic wave is walked the present invention relates to seimic travel time calculating field.
Background technology
《Journal of Engineering Geophysics》3rd phase in 2009 disclose Zhang Shuanjie etc. " fast-marching algorithm computational accuracy analyze and Improve ", the method for describing two kinds of raising fast-marching algorithm computational accuracies, first, being calculated using higher-order wave equation;Two It is that local fine processing is carried out to grid node near focal point, is calculated near focus using higher-order wave equation, Remaining area is calculated using low order differential form.And it carries out seismic wave to uniform dielectric model by both approaches to walk When calculate, experimental result has obtained relatively good effect.
《International Geology》3rd phase in 2016 discloses " the fast-marching algorithm seismic wave based on complete ternary tree such as Wang Qianlong Calculated when walking ", describe it is a kind of it is improved quickly propel seimic travel time computational methods, complete ternary tree sort method is introduced In being calculated to seimic travel time, reduce and the time that minimum walks duration is found in narrowband extension, improve the calculating of entire algorithm Efficiency.And by quickly propelling seimic travel time computational methods to stratified model, Marmousi moulds based on complete ternary tree Type, Sigsbee 2a models are calculated, and experimental result has obtained relatively good effect.
By example above as can be seen that the existing seimic travel time computational methods that quickly propel can carry to a certain extent Computational accuracy is risen, but realizes that process is complicated, the computational accuracy of promotion is also limited.
The content of the invention
The technical problems to be solved by the invention are to provide a kind of mixing computational methods when two-dimensionally seismic wave is walked, by flexible By wavefront construction seimic travel time computational methods and the seamless connection of seimic travel time computational methods is quickly propelled with narrow-band technologies Get up, i.e., by near focus a small range using the higher wavefront construction method of computational accuracy, used in remaining area quick When propelling method calculates, while original fast-marching algorithm computational accuracy is improved, the characteristics of its is efficient is still remained.
In order to solve the above technical problems, the technical solution adopted by the present invention is:
It is a kind of to mix computational methods when two-dimensionally seismic wave is walked, comprise the following steps:
Step 1:Relevant parameter file, rate pattern are read in, wherein, the Parameter File includes the mesh point of rate pattern Number, grid spacing and hypocentral location;
Step 2:Along different directions divergent-ray and central ray information is calculated from shot point;
Step 3:When being walked using grid node in wavefront construction method calculating ray coverage;
Step 4:Classify to all mesh points, divide them into receiving station, narrowband point or point of distance, i.e., if one The seimic travel time of mesh point had been computed, and that is put around it all calculated when walking, then this point is receiving station;Such as The seimic travel time of one mesh point of fruit calculated, and put around it when walking and it is not all calculated, then this point is narrowband Point;If point is not calculated when walking, this point is point of distance;
Step 5:All narrowband points are moved into narrowband, build initial narrow band;
Step 6:Narrowband is extended, until narrowband is sky;In this course, when mesh point is walked asked by upwind difference method The two-dimentional eikonal equation of solution obtains, and the two dimension eikonal equation is:
| ▽ τ |=s
Wherein, τ is seimic travel time, and s is slowness, and ▽ is gradient signs, the upwind difference table of gradient terms in above-mentioned formula It is up to formula:
Wherein,Be respectively when walking at the grid node (i, j) on x or z directions Difference expression forward or backward, concrete form are respectively:
Further, in step 2, solve kinematics ray tracing equation group using Long Gekutafa and obtain central ray Information is shown below:
Wherein, xiRepresentation space position, piRepresent slowness component, τ represents seimic travel time, and v represents the speed at discrete point Value.
Compared with prior art, the beneficial effects of the invention are as follows:It is improved as a result of wavefront construction method near focus Area grid node seimic travel time computational accuracy, so as to improve the computational accuracy of fast-marching algorithm zoning mesh point; The computational accuracy of entire seimic travel time algorithm greatly improved to reduce less computational efficiency as cost.
Description of the drawings
Fig. 1 is that the present invention mixes computational methods flow chart when two-dimensionally seismic wave is walked.
Fig. 2 wavefront construction method region segmentation schematic diagrames.
Fig. 3 is fast-marching algorithm relative error in uniform dielectric.
Fig. 4 is that seimic travel time computational methods relative error is mixed in uniform dielectric.
Specific embodiment
The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.
Fig. 1 is to mix computational methods flow chart when two-dimensionally seismic wave is walked, and the realization flow of the method for the present invention is shown in figure, It is specific as follows:
1) relevant parameter file, rate pattern are read in, wherein, Grid dimension of the Parameter File comprising rate pattern, Grid spacing and hypocentral location.
2) from shot point along different directions divergent-ray, the launch angle scope of ray is:- 90 ° to+90 °, between ray Angle at intervals of 3 ° to 10 °, using Long Gekutafa solve kinematics ray tracing equation group obtain central ray information, such as Shown in following formula:
Wherein, xiRepresentation space position, piRepresent slowness component, τ represents seimic travel time, and v represents the speed at discrete point Value.
3) when being walked using grid node in wavefront construction method calculating ray coverage.This region quilt in calculating process It is put before adjacent wave on adjacent center ray and is divided into multiple irregular quadrilaterals (as shown in Figure 2), in each quadrangle Grid node is obtained by quadrangle vertex information interpolation.
4) complete after wavefront construction method calculates, attributive classification is carried out to total space grid node, it is specified that:An if grid The seimic travel time of point had been computed, and that is put around it all calculated when walking, then this point is receiving station;If one The seimic travel time of a mesh point calculated, and put around it when walking and it is not all calculated, then this point is narrowband point;Such as One point of fruit is not calculated when walking, then this point is point of distance.
5) all narrowband points are moved into narrowband, builds initial narrow band.
6) narrowband is extended, until narrowband is sky.In this course, it is to solve two by upwind difference method when mesh point is walked Tie up what eikonal equation obtained, the two dimension eikonal equation is:
| ▽ τ |=s
Wherein, τ is seimic travel time, and s is slowness, and ▽ is gradient signs, the upwind difference table of gradient terms in above-mentioned formula It is up to formula:
Wherein,Be respectively when walking at the grid node (i, j) on x or z directions Difference expression forward or backward, concrete form are respectively:
Analysis verification is carried out to the computational accuracy and stability of the method for the present invention below by uniform dielectric model.
Fig. 3, Fig. 4 are respectively that fast-marching algorithm is opposite in uniform dielectric model with mixing seimic travel time computational methods Error, the size of homogeneous model is 601 × 601, and grid spacing is 10m × 10m, speed 1000m/s.When combining to walk in method The maximum length for finding out single ray is 500m.As can be seen from the figure seimic travel time computational methods are mixed compared with quickly pushing away It is greatly improved into method computational accuracy.
By the present invention in that being connected to fast-marching algorithm and wavefront construction method with narrow-band technologies, improved by wavefront construction method The seimic travel time computational accuracy of the neighbouring smaller area of focus, so as to improving fast-marching algorithm in remaining area grid node Computational accuracy, realize when a kind of two-dimensionally seismic wave for taking into account computational efficiency and computational accuracy is walked and calculate.

Claims (2)

1. a kind of mix computational methods when two-dimensionally seismic wave is walked, which is characterized in that comprises the following steps:
Step 1:Relevant parameter file, rate pattern are read in, wherein, Grid dimension of the Parameter File comprising rate pattern, Grid spacing and hypocentral location;
Step 2:Along different directions divergent-ray and central ray information is calculated from shot point;
Step 3:When being walked using grid node in wavefront construction method calculating ray coverage;
Step 4:Classify to all mesh points, divide them into receiving station, narrowband point or point of distance, i.e., if a grid The seimic travel time of point had been computed, and that is put around it all calculated when walking, then this point is receiving station;If one The seimic travel time of a mesh point calculated, and put around it when walking and it is not all calculated, then this point is narrowband point;Such as One point of fruit is not calculated when walking, then this point is point of distance;
Step 5:All narrowband points are moved into narrowband, build initial narrow band;
Step 6:Narrowband is extended, until narrowband is sky;In this course, it is to solve two by upwind difference method when mesh point is walked Tie up what eikonal equation obtained, the two dimension eikonal equation is:
| ▽ τ |=s
Wherein, τ is seimic travel time, and s is slowness, and ▽ is gradient signs, the upwind difference expression formula of gradient terms in above-mentioned formula For:
<mrow> <mo>|</mo> <mo>&amp;dtri;</mo> <mi>&amp;tau;</mi> <mo>|</mo> <mo>=</mo> <msup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>max</mi> <msup> <mrow> <mo>(</mo> <msubsup> <mi>D</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> <mrow> <mo>-</mo> <mi>x</mi> </mrow> </msubsup> <mi>&amp;tau;</mi> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <mi>min</mi> <msup> <mrow> <mo>(</mo> <msubsup> <mi>D</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> <mrow> <mo>+</mo> <mi>x</mi> </mrow> </msubsup> <mi>&amp;tau;</mi> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mi>max</mi> <msup> <mrow> <mo>(</mo> <msubsup> <mi>D</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> <mrow> <mo>-</mo> <mi>z</mi> </mrow> </msubsup> <mi>&amp;tau;</mi> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <mi>min</mi> <msup> <mrow> <mo>(</mo> <msubsup> <mi>D</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> <mrow> <mo>+</mo> <mi>z</mi> </mrow> </msubsup> <mi>&amp;tau;</mi> <mo>,</mo> <mn>0</mn> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mrow> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </msup> </mrow>
Wherein,Be respectively when walking at the grid node (i, j) on x or z directions forward Or backward difference expression formula, concrete form are respectively:
2. a kind of as described in claim 1 mix computational methods when two-dimensionally seismic wave is walked, which is characterized in that in step 2, makes Kinematics ray tracing equation group is solved with Long Gekutafa and obtains central ray information, is shown below:
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <msub> <mi>dx</mi> <mi>i</mi> </msub> </mrow> <mrow> <mi>d</mi> <mi>&amp;tau;</mi> </mrow> </mfrac> <mo>=</mo> <msup> <mi>v</mi> <mn>2</mn> </msup> <mo>&amp;CenterDot;</mo> <msub> <mi>p</mi> <mi>i</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <msub> <mi>dp</mi> <mi>i</mi> </msub> </mrow> <mrow> <mi>d</mi> <mi>&amp;tau;</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> </mrow> </mfrac> <mrow> <mo>(</mo> <mfrac> <mn>1</mn> <mi>v</mi> </mfrac> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <msup> <mi>v</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>v</mi> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
Wherein, xiRepresentation space position, piRepresent slowness component, τ represents seimic travel time, and v represents the velocity amplitude at discrete point.
CN201810077621.8A 2018-01-23 2018-01-23 It is a kind of to mix computational methods when two-dimensionally seismic wave is walked Pending CN108072897A (en)

Priority Applications (4)

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
NL2020684A NL2020684B1 (en) 2018-01-23 2018-03-29 Mixed 2D seismic wave travel time calculation method
BE2018/0040A BE1025488B1 (en) 2018-01-23 2018-03-30 Method for calculating 2D mixed seismic propagation time
LU100751A LU100751B1 (en) 2018-01-23 2018-03-30 Method for calculating 2D mixed seismic propagation time

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 (1)

Publication Number Publication Date
CN108072897A true CN108072897A (en) 2018-05-25

Family

ID=61978378

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810077621.8A Pending CN108072897A (en) 2018-01-23 2018-01-23 It is a kind of to mix computational methods when two-dimensionally seismic wave is walked

Country Status (4)

Country Link
CN (1) CN108072897A (en)
BE (1) BE1025488B1 (en)
LU (1) LU100751B1 (en)
NL (1) NL2020684B1 (en)

Cited By (3)

* Cited by examiner, † Cited by third party
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
CN117194855A (en) * 2023-11-06 2023-12-08 南方科技大学 Fitting analysis method and relevant equipment for weak anisotropy travel time
CN117607957A (en) * 2024-01-24 2024-02-27 南方科技大学 Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002065159A1 (en) * 2001-02-14 2002-08-22 Hae-Ryong Lim Method of seismic imaging using direct travel time computing
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
CN106353793A (en) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 Cross-well seismic tomography inversion method on basis of travel time incremental bilinear interpolation ray tracing

Family Cites Families (3)

* Cited by examiner, † Cited by third party
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
US8094513B2 (en) * 2008-06-03 2012-01-10 Westerngeco L.L.C. Determining positioning of survey equipment using a model

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002065159A1 (en) * 2001-02-14 2002-08-22 Hae-Ryong Lim Method of seismic imaging using direct travel time computing
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

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
丁鹏程等: "《基于三维多模板快速推进算法的复杂近地表射线追踪》", 《石油物探》 *
孙章庆等: "《三维复杂山地多级次群推进迎风混合法多波型走时计算》", 《地球物理学报》 *
李永博: "《快速行进法射线追踪提高旅行时计算精度和效率的改进措施》", 《石油地球物理勘探》 *

Cited By (5)

* Cited by examiner, † Cited by third party
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
CN117194855A (en) * 2023-11-06 2023-12-08 南方科技大学 Fitting analysis method and relevant equipment for weak anisotropy travel time
CN117194855B (en) * 2023-11-06 2024-03-19 南方科技大学 Fitting analysis method and relevant equipment for weak anisotropy travel time
CN117607957A (en) * 2024-01-24 2024-02-27 南方科技大学 Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method
CN117607957B (en) * 2024-01-24 2024-04-02 南方科技大学 Seismic wave travel time solving method and system based on equivalent slowness rapid propulsion method

Also Published As

Publication number Publication date
NL2020684B1 (en) 2018-11-14
BE1025488A1 (en) 2019-03-18
BE1025488B1 (en) 2019-03-25
NL2020684A (en) 2018-04-20
LU100751B1 (en) 2018-10-01

Similar Documents

Publication Publication Date Title
CN107392875A (en) A kind of cloud data denoising method based on the division of k neighbours domain
CN104318622B (en) Triangular mesh modeling method of indoor scene inhomogeneous three dimension point cloud data
Ben‐Chen et al. Conformal flattening by curvature prescription and metric scaling
CN106528740B (en) Road axis extracting method based on Delaunay triangulation network
CN107610223A (en) Power tower three-dimensional rebuilding method based on LiDAR point cloud
CN104200212B (en) A kind of building external boundary line drawing method based on airborne LiDAR data
CN105005995B (en) A kind of method for calculating three-dimensional point cloud model bone
CN108072897A (en) It is a kind of to mix computational methods when two-dimensionally seismic wave is walked
CN108734728A (en) A kind of extraterrestrial target three-dimensional reconstruction method based on high-resolution sequence image
Alharthy et al. Detailed building reconstruction from airborne laser data using a moving surface method
CN107886528A (en) Distribution line working scene three-dimensional rebuilding method based on a cloud
CN104700451A (en) Point cloud registering method based on iterative closest point algorithm
CN103729872B (en) A kind of some cloud Enhancement Method based on segmentation resampling and surface triangulation
CN107843922A (en) One kind is based on seismic first break and the united chromatography imaging method of Travel time
CN108711194A (en) A kind of three-dimensional grid model joining method based on cubic Bézier curves
CN105975974A (en) ROI image extraction method in finger vein identification
CN105303612B (en) A kind of extract digital network method based on Triangulated irregular network model
CN105354881B (en) Distortion of the mesh optimized algorithm based on Category Attributes data
CN107424166A (en) Point cloud segmentation method and device
CN114332291A (en) Oblique photography model building outer contour rule extraction method
CN110095775A (en) The platform SAR fast time-domain imaging method that jolts based on mixed proportion
CN108957538A (en) A kind of virtual focus two dimension wavefront construction seimic travel time calculation method
CN110147575A (en) A kind of calculation method that the two-phase stream interface based on single layer particle levels collection captures
CN110176061A (en) Human body surface reconstructing method in a kind of three-dimensional reconstruction
CN115294287A (en) Laser SLAM mapping method for greenhouse inspection robot

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20180525

RJ01 Rejection of invention patent application after publication