CN102495427B - Interface perception ray tracing method based on implicit model expression - Google Patents
Interface perception ray tracing method based on implicit model expression Download PDFInfo
- Publication number
- CN102495427B CN102495427B CN 201110414879 CN201110414879A CN102495427B CN 102495427 B CN102495427 B CN 102495427B CN 201110414879 CN201110414879 CN 201110414879 CN 201110414879 A CN201110414879 A CN 201110414879A CN 102495427 B CN102495427 B CN 102495427B
- Authority
- CN
- China
- Prior art keywords
- stratum
- ray
- node
- rule
- tetrahedron
- 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
Images
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses an interface perception ray tracing method based on implicit model expression and is applicable to a system for inverting, modeling and offsetting seismic data. The method comprises the following steps of: expressing a geological interface based on an implicit function; describing a formation depositional sequence by using a binary tree, and generating a three-dimensional geologic body construction model by using the geological interface which is implicitly expressed; on the basis of an implicitly-expressed construction model, assigning speed values to vertexes of a tetrahedron unit which is included in a formation body or intersected with a formation surface according to an appointed speed distribution rule of each formation body, and organizing and storing speed files by taking the formation body as a unit; and under the constraint of the implicitly-expressed construction model, performing ray tracing based on interface perception in the construction model in a segmented iteration mode according to a kinematical equation and a reflection law. The method has the advantages that: a complicated geological construction model can be subjected to automatic velocity modeling; high-accuracy and powerful ray tracing is realized; and important parameters are provided for seismic inversion, offset imaging and numerical simulation.
Description
Technical field
The present invention relates to automatically carry out the System and method for of velocity modeling and ray tracing, particularly a kind of interface perception ray tracing method based on implicit model expression under a kind of geological structure model constraint in implied expression.
Background technology
Because ray is carrying in a large number for seismic inversion, modeling and is being offset very Useful Information, so ray-tracing scheme is being played the part of very important role in Seismic Exploration Data Processing.If underground speed is known, the parameters such as raypath, whilst on tour, amplitude and emergence angle can obtain by the numerical evaluation of ray tracing.
Just proposed ray-tracing scheme as far back as the people such as Cerveny in 1977, Cerveny afterwards and Hron have carried out again replenishing and are perfect, and it is theoretical that these have consisted of classical ray tracing.Ray tracing is also very direct from the mathematical theory to the computer realization, but until in recent years, the application of three-dimensional ray tracing in actual seismic inverting and modeling is also very limited, that scientific research field carries out simulative example research mostly, for example in order to test new disposal route generating three-dimensional theogram, or check certain hypothesis and approximate rationality.So the major reason that this method can not be widely used at petroleum industrial circle, using for real complex geological structure area three-dimensional seismic prospecting is to expend very much computational resource, and the error of classic method is also very large.Existing many velocity modeling methods based on grid are used for ray tracing and tomographic inversion both at home and abroad.Chapman had also introduced a kind of velocity modeling algorithm based on ray in 2004, can perception aspect framework, and can be along the ray propagates physical attribute when ray runs into aspect.But these classic methods are usually based on tetrahedral grid, and the generation of grid is retrained by section and aspect, comprises unconformability and overlying strata etc., and in the area of complicated underground structure, grid generates often very difficult and generating mesh is of low quality.2006, the people such as Ruger propose a kind of method that automatically generates the formation velocity structured grid from the even speed grid of meticulous sampling, the method need not the constraint information (as tomography and stratum) that provides extra, the speed grid is embedded in the tetrahedral grid of Delaunay segmentation at discontinuity zone, but they fundamentally do not solve the difficult problem that the constraint grid generates.Traditional ray tracing is divided into workflow several different stages usually: pick up structural attitude, formation speed model and ray tracing, except the input of output meeting as next stage of previous stage, stages is often separate, there is no other positive connection.
The invention discloses a kind of interface perception ray tracing method based on implicit model expression, three Main Stage of ray tracing: structure modeling, velocity modeling and ray tracing have unified method basis---implied expression method, can not only the effective expression complex geological structure, and each stage organic connections, except being the last period the input of the latter half, also for the latter half provides support structure and computing method to support, consist of a unified integral body.
Summary of the invention
The technical problem to be solved in the present invention is: the method for velocity modeling and high precision ray tracing is carried out in research automatically in the complex geological structure model, the system that is used for seismic exploration data inverting, modeling and skew, solving the problem of existing ray-tracing scheme performance difficulty in complex geological structure, thereby overcome the difficulty that in prior art, precision is not high, efficient is lower and be difficult to effectively support seismic inversion, skew and simulation.
The present invention solves the technical scheme that its technical matters takes: a kind of interface perception ray tracing method based on implicit model expression is provided, it is under the tectonic model constraint of implicit model expression, carry out velocity modeling and interface perception ray tracing, comprise the following steps:
Steps A, geological interface and geologic body generation method construct three-dimensional geologic structure model based on implicit expression comprise symbolic distance field and stratum binary tree;
Step B, desin speed model automatically on the three-dimensional structure model basis of implied expression;
Step C carries out ray tracing according to tectonic model and rate pattern.
Symbolic distance field in steps A is the implied expression form of tectonic model, and the zero contour surface of symbolic distance field is in order to represent corresponding geological stratum; The stratum binary tree is the data structure of describing geological stratum and geologic body spatial topotaxy, and have following character: each non-leaf node represents an aspect, and each leaf node represents the area of space that aspect is cut apart, i.e. stratum body; Main stratum is father or the ancestor node on respective secondary stratum always; To arbitrary stratum binary tree node, its left sibling, oneself, its aspect corresponding to right node that always becomes sequence or reverse sequence arranges.
Step B further comprises:
Step B1, be that each stratum body is specified tetrahedral grid according to stratum binary tree and symbolic distance field, these tetrahedral grids comprise and are positioned at the tetrahedral grid that the stratum body is inner and intersect with stratum body outland aspect fully, crossing tetrahedral grid with its Different Strata body that intersects in respectively preserve a duplicate;
Step B2 is the tetrahedral grid summit tax velocity amplitude of each stratum body appointment according to the speed create-rule, and take the stratum body as unit organization speed file; Speed generates rule and is jointly determined by stratal surface base velocity function and velocity gradient, and wherein velocity gradient derives from sand smeller and geophysicist's priori or obtained by migration velocity analysis.
Step C further comprises:
Step C1 according to kinematical equation, adopts the method for cutting apart when waiting, and the raypath segmentation is calculated; Each calculates segmentation and is called a ray step, and to each ray step, if itself and the border triangle intersect of tetrahedron element, it will be blocked by this face so, and joining becomes this ray and goes on foot new terminal point.
Step C2, check whether the ray step is crossing with tomography: because fault surface adopts tetrahedral grid face Explicit Expression, and cancelled tomography both sides tetrahedron topological adjacency relation, the adjacent tetrahedron of sharing fault surface keeps a common sides duplicate separately, when ray step and tetrahedral border triangle intersect, check whether this face is fault surface, if, this ray step terminal point is labeled as and is positioned at tomography, namely this ray step intersects with tomography;
Step C3, check whether the ray step is crossing with stratal surface, because stratal surface has adopted the hidden function representation that shows, can not directly judge crossing, can use the rule of judgement spatial point layer position to determine in the binary tree of stratum: if starting point and the terminal point in ray step are positioned at the same leaf node of stratum binary tree, can determine that this ray step is positioned at certain stratum body fully and not crossing with stratal surface; Otherwise this ray step intersects with a stratal surface at least.Further determine the definite position of intersection point: check each node on the ancestral path of tracing back from this ray step stratum, starting point place body leaf node to root node, for nonleaf node, if the starting point in ray step and terminal point be opposite in sign in the symbolic distance field, can determine that ray step and this stratal surface intersect, the ratio of the starting point in ray step and the terminal point value in the symbolic distance field can be determined the definite position of intersection point.
Reflection and refraction when intersect on ray step and tomography or stratum, will occur in step C4, use the symbolic distance field, can calculate rapidly the normal vector at ray and aspect joining place, with conveniently calculating ray incident angle and emergence angle.When the ray step intersected at tomography, the reflected ray step still carried out in the current tetrahedron of current stratum body, but has new directions of rays; Refracted rays step will carry out in sharing the adjacent tetrahedron of fault surface with current tetrahedron, can determine stratum, new tetrahedron place body according to the rule of judgement spatial point layer position, if not at current stratum body, need to load corresponding stratum body speed file, and use the speed create-rule of new stratum body to go on foot point interpolation as refracted rays to go out velocity amplitude.When the ray step intersected at the stratum, the reflected ray step still carried out in the body of current stratum, but has new directions of rays; Refracted rays step will originate in new stratum body, and new stratum body uses the rule of the judgement spatial point layer position of revising to determine, and loads corresponding stratum body speed file, uses new tetrahedron duplicate to go on foot point interpolation as refracted rays and goes out velocity amplitude.
The rule of the judgement spatial point layer position in step C3 and C4 refers to: begin the selectivity traversal from stratum binary tree root node, gone out the value of given spatial point by interpolation in the symbolic distance field of aspect corresponding to current binary tree node, symbol by this value determines that the child node of next traversal is left child node or right child node, so continue traversal, namely provided layer position, spatial point place until arrive at leaf node; Referring to of the rule of the judgement spatial point layer position of revising: to the nonleaf node that intersects the stratum, if current stratum body is positioned at left subtree, continue traversal from right subtree when traversal, vice versa.
The advantage that the present invention compared with prior art has is:
1, the present invention because the implicit method that adopts is expressed tectonic model, can directly serve velocity modeling and ray tracing, is applicable to the complex geological structure environment;
2, in velocity modeling of the present invention, adopt multiple duplicate method, take the stratum body as unit formation speed file, can embed sand smeller and geophysicist to the understanding of geologic rule, can realize that the interruption of bed boundary place speed expresses, support for ray tracing has for strong;
3, the present invention utilizes symbolic distance field and stratum binary tree, has formed a whole set of rational ray and has intersected, reflects and the refraction rule, has effectively improved precision and the efficient of ray tracing.
Description of drawings
Fig. 1 is tectonic structure section and the corresponding stratum binary tree that generates;
Fig. 2 is the tetrahedron schematic diagram that intersects with stratal surface;
Fig. 3 is that the rate pattern with dual duplicate is expressed schematic diagram;
Fig. 4 is the complex geological structure model with a plurality of tomographies;
Fig. 5 is reflection, the refraction schematic diagram of injection line in the complex structure model;
Fig. 6 is the propagation schematic diagram that in the complex structure model, inner point source ray arbitrarily arrives certain section;
Fig. 7 is the complex geological structure model with the mushroom body;
Fig. 8 is inner point source ray arbitrarily at propagation (only the considering reflection) schematic diagram with the complex geological structure model of mushroom body;
Fig. 9 diagrammatic cross-section that to be ray propagate in mushroom body border has the tectonic model of high versus speed;
Figure 10 is rate pattern corresponding to actual seismic exploration geology tectonic model;
Figure 11 is propagation (only the consider refraction) schematic diagram of inner point source ray arbitrarily in actual seismic exploration geology tectonic model;
The structure wavefront surface schematic diagram that Figure 12 propagates in actual seismic exploration geology tectonic model for inner point source ray arbitrarily.
Embodiment
Below in conjunction with description of drawings the specific embodiment of the present invention.
A kind of interface perception ray tracing method based on implicit model expression comprises three steps, comprising:
Steps A, geological interface and geologic body generation method construct three-dimensional geologic structure model based on implicit expression comprise symbolic distance field and stratum binary tree;
wherein, in steps A based on geological interface and the geologic body generation method of implicit expression, can adopt " a kind of stratal surface and geologic body generation method based on level set " (to obtain national patent, license number is ZL200810112263.6), the mathematical method that also can take other generates corresponding symbol data field (as Frank for each aspect, T., Tertois, A.L., Mallet, J.L., 2007.3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data.Computer﹠amp, Geosciences 33:932-943, Calcagno, P., Chiles, J., Courrioux, G., Guillen, A., 2008, Geological modeling from field data and geological knowledge-part I:Modelling method coupling 3D potential-field interpolation and geological rules:Physics of the Earth and Planetary Interiors 171 (1-4), 147-157.) aspect is as the zero contour surface of symbol data field.
Implied expression method symbolization distance field is as representation, expresses implicit surface with the zero contour surface of symbolic distance field.In the present invention, be each stratal surface structure symbolic distance field, the zero contour surface of symbolic distance field is in order to express geological stratum.Simultaneously, defined a kind of data structure of describing geological stratum and geologic body spatial topotaxy---stratum binary tree, it has following character: each non-leaf node represents an aspect, and each leaf node represents the area of space that aspect is cut apart, i.e. stratum body; Main stratum is father or the ancestor node on respective secondary stratum always; To arbitrary stratum binary tree node, its left sibling, oneself, its aspect corresponding to right node that always becomes sequence or reverse sequence arranges.Fig. 1 is the schematic diagram of tectonic structure section and its corresponding stratum binary tree that generates, wherein, and Si (i=1,2,3,4) represent stratal surface, corresponding to the upper surface of stratum body, Ti (i=1,2 ..., 9) represent by the space subregion of stratigraphic division, in figure, stratum S2 is divided into T1 and T2 two subregions with area of space, and two sub regions are divided respectively again.
Step B, desin speed model automatically on the three-dimensional structure model basis of implied expression;
Step B1, be that each stratum body is specified tetrahedral grid according to stratum binary tree and symbolic distance field, these tetrahedral grids comprise and are positioned at the tetrahedral grid that the stratum body is inner and intersect with stratum body outland aspect fully, crossing tetrahedral grid with its Different Strata body that intersects in respectively preserve a duplicate;
Step B2 is the tetrahedral grid summit tax velocity amplitude of each stratum body appointment according to the speed create-rule, and take the stratum body as unit organization speed file; Speed generates rule and is jointly determined by stratal surface base velocity function and velocity gradient, and wherein velocity gradient derives from sand smeller and geophysicist's priori or obtained by migration velocity analysis.
Velocity modeling is a kind of attribute modeling in essence, and is very high to the rate pattern accuracy requirement when rate pattern is applied to the ray tracing of interface perception, and particularly at the interface speed is expressed on tomography, stratum etc.Existing many ray tracings that are used for the aspect perception based on the velocity modeling method of grid, but these main difficulties based on the interface perception velocities modeling algorithm of grid are the model constrained lower grid generation of complex structure.Be that based on a main challenge of the velocity modeling of implied expression the implicit expression stratum is not the tetrahedral grid interface, does not have grid vertex corresponding to stratal surface.Searching the affiliated stratum of each grid vertex and be their assignment, is a kind of consuming time and clumsy way by enumerating or comparing with structure at all levels.How the grid that those and aspect are intersected has different solutions for each summit tax velocity amplitude.Because not all four summits are all from same stratum body, the method by direct linear interpolation is that internal point is composed the result that velocity amplitude can lead to errors.The people such as Bargteil in 2007 have proposed a kind of method during a class Similar Problems in rebuilding solving geologic body, the tetrahedron that the method is positioned at above or below the stratum fully is set to constant, for and the tetrahedron that intersects of aspect adopt the approximate strategy of a kind of body according to the ratio between bi-material.For the solid material of plastic yield, this is a kind of simple and effective method, but and is not suitable for the velocity modeling of ray tracing.Because in the ray tracing process, the crossing calculating of a large amount of rays and aspect is arranged, the inaccurate calculating that can affect directions of rays of the speed that aspect is adjacent.
In isotropic medium, speed can be defined have the constant gradient, exponential function, conical function or analyzable feature.Because the deposition characteristics on stratum makes tectonic structure have transversely isotropic symmetrical feature, at vertical depth direction, speed may be no longer the constant gradient.Yet in dynamic construction zone, through long-term structural evolution, the stratum may run-off the straight, thereby causes the inclination of lateral symmetry axle.Most seismic prospecting research all is based on the layered geology model, and due to deposition characteristics, in layer, speed is linear change.Yet increase with the vertical degree of depth, always the variation of speed may be not continuous, when particularly running into the stratum.Therefore, need a kind of rate pattern of learning understanding of structure with meeting, make speed variation continuously in layer, and may be interrupted in the aspect junction.
Based on above geological knowledge, the rate pattern that is used for ray tracing and inverting must be consistent with the distribution of stratum media.Therefore, velocity modeling mainly comprises two parts: for the stratum body is specified tetrahedral grid, and be that grid vertex is composed velocity amplitude according to inclination transverse isotropy principle.The present invention has adopted a kind of velocity gradient method to calculate the speed of grid vertex in each stratum body.At first, for each stratum body round floor is specified a 2D base velocity function v (x, y), basic velocity amplitude on the bed interface, base area and velocity gradient (be generally sand smeller and geophysicist's priori or obtained by migration velocity analysis), use following formula to calculate the interior speed of stratum body:
v(P)＝v
0(P)+g*h(P)
In formula, the velocity amplitude of the tetrahedral grid summit P that v (P) indicates to ask for, v
0(P) for specifying the basic velocity amplitude of bed boundary, g is the velocity gradient of stratum body, and h (P) expression summit P is to the normal distance on stratum.Obviously, h (P) can be directly calculated by the symbolic distance field of stratigraphic structure model, tracks stratal surface according to symbolic distance field negative gradient direction from P equally, can obtain stratum base speed v
0(P).
Rate pattern is subdivided into tetrahedral grid, in the ray tracing process, need to use ray the velocity amplitude of each grid vertex of process.The velocity amplitude of tetrahedron element internal point is got by the velocity amplitude linear interpolation on four summits, but for the tetrahedral grid that intersects with aspect, this interpolation method is just no longer suitable.Because bed boundary does not have the constraint information that explicitly is defined as tetrahedral grid, and can not be expressed by tetrahedral triangular topological relations with the stratal surface of implicit method expression.Therefore, if only use the simple interpolations method, the velocity amplitude of grid and aspect intersection just may incorrectly maybe can't be explained.As shown in Fig. 2 a, intersect on tetrahedral grid unit and stratum, and summit A is positioned on stratal surface, and other three summit B, C and D are positioned under stratal surface, will calculate now the velocity amplitude of internal point E.If directly use the velocity amplitude on four summits to carry out linear interpolation, the E spot speed value that obtains so may be incorrect or be made us and can't explain, because the A point uses different velocity gradients calculate and get with B, C, D three spot speed values.The ray tracing of the mistake that the velocity amplitude of mistake will cause, and these bed boundaries places is the very responsive zone of ray tracing, very important to ray tracing, it will determine the direction of ray.
The present invention has adopted the method for dual duplicate to solve this problem, velocity amplitude for the point in the crossing grid cell on the calculating aspect, stratum body speed create-rule on aspect is extended to all summits (comprise on aspect and under) of this crossing grid cell, and the velocity amplitude on these summits is recorded in aspect overlying strata body speed file.For the summit under aspect of being positioned at of intersecting grid cell, make and use the same method, stratum body speed create-rule under aspect is extended to all summits of this crossing grid cell, and the velocity amplitude on these summits is recorded in the speed file of aspect sub-surface body.The situation that intersects for more complicated stage construction and grid cell, as Fig. 2 b, the method also is easy to expand to the method for multiple duplicate, for each stratum body generates the speed file that different duplicates are arranged, that is, if tetrahedron and a plurality of stratum body intersect, have a plurality of these tetrahedral friction speed value duplicates and exist, and be recorded in the speed file of corresponding stratum body, Fig. 3 describes the method for speed assignment.
Therefore, the velocity modeling process also needs respective extension.To each stratum body, need to search the tetrahedral grid that is positioned at the stratum body and intersects with its stratal surface, and the speed create-rule that uses this stratum body is tetrahedral grid summit assignment.Main difficulty in realization is how to search these tetrahedrons.Specific practice is: to each stratum body, travel through all tetrahedral grids, use symbolic distance field and stratum binary tree, determine whether tetrahedron is included in the stratum body interior or crossing with stratal surface, usually, the accessed path in the binary tree of stratum is for tracing back to the path of root node from the leaf node that represents this stratum body.For example, check whether tetrahedron C is included in stratum body S interior or crossing with it, from the leaf node of expression stratum body S, define the process that an iterator is indicated its traversal in the ergodic process in the binary tree of stratum.At first, iterator points to the father node (in the binary tree of stratum, nonleaf node represents geological stratum) of this leaf node, if stratum body S is positioned at a side on the indicated stratum of iterator, and all summits of tetrahedron C are positioned at opposite side, can determine that this tetrahedron is not included in the body S of stratum or intersect with it, traversal stops; Otherwise upgrade the position of iterator, along ancestors' path tracing, until root node, if trace back process is never violated Rule of judgment, can determine that tetrahedron S is included in the body S of stratum or with it intersects.When find all be included in stratum body S or with its tetrahedron that intersects after, be this stratum body tax velocity amplitude according to the speed create-rule, and generate its corresponding speed file.
Step C carries out ray tracing according to tectonic model and rate pattern.
Step C1 according to kinematical equation, adopts the method for cutting apart when waiting, and the raypath segmentation is calculated; Each calculates segmentation and is called a ray step, and to each ray step, if itself and the border triangle intersect of tetrahedron element, it will be blocked by this face so, and joining becomes this ray and goes on foot new terminal point.
The kinematical equation of ray tracing belongs to eikonal equation, and Cerveny use characteristic method has been described its concrete form, as shown in the formula:
du＝v
ndτ
Can find out from above-mentioned equation, along ray tracing x
iAnd distribution p
iIt is the function of monotone increasing independent variable u.x
i(u) be called as raypath, use this prescription journey can obtain the iterative computation formula of raypath.Ray tracing is from source point, and one group of analog wave front construction method emission is the ray of constant increment angle APPROXIMATE DISTRIBUTION.No matter when, when surpassing certain distance, a new ray can be inserted into to safeguard wavefront shape when the adjacent ray of topology.Cut apart segment iteration during raypath employing etc. and calculate, each calculates segmentation and is called a ray step.To each ray step, if its certain face with tetrahedron element intersects, it will be blocked by this face so, and joining becomes this ray and goes on foot new terminal point.Simultaneously, carry out mark to this intersection point, see whether its phase cross surface is positioned at tomography, because refraction or reflection will occur at the tomography place ray.
Step C2, check whether the ray step is crossing with tomography: because fault surface adopts tetrahedral grid face Explicit Expression, and cancelled tomography both sides tetrahedron topological adjacency relation, the adjacent tetrahedron of sharing fault surface keeps a common sides duplicate separately, when ray step and tetrahedral border triangle intersect, check whether this face is fault surface, if, this ray step terminal point is labeled as and is positioned at tomography, namely this ray step intersects with tomography;
Step C3, check whether the ray step is crossing with stratal surface, because stratal surface has adopted the hidden function representation that shows, can not directly judge crossing, can use the rule of judgement spatial point layer position to determine in the binary tree of stratum: if starting point and the terminal point in ray step are positioned at the same leaf node of stratum binary tree, can determine that this ray step is positioned at certain stratum body fully and not crossing with stratal surface; Otherwise this ray step intersects with a stratal surface at least.Further determine the definite position of intersection point: check each node on the ancestral path of tracing back from this ray step stratum, starting point place body leaf node to root node, for nonleaf node, if the starting point in ray step and terminal point be opposite in sign in the symbolic distance field, can determine that ray step and this stratal surface intersect, the ratio of the starting point in ray step and the terminal point value in the symbolic distance field can be determined the definite position of intersection point.
Therefore be not expressed as grid surface because stratal surface has explicitly, need further to check whether the ray step intersect with a certain stratal surface, the rule of judgement spatial point layer position can help to realize.If starting point and the terminal point in ray step are positioned at the same leaf node of stratum binary tree, can determine that this ray step is positioned at certain stratum body fully and does not intersect with stratal surface; Otherwise this ray step intersects with a stratal surface at least.In order to determine the definite position of intersection point, at first need to find out this ray step terminal point position.Therefore, in the binary tree of stratum, each node on the ancestral path of tracing back from this ray step stratum, starting point place body leaf node to root node may become checked object.For each nonleaf node (stratal surface), if the starting point in ray step and terminal point be opposite in sign in the symbolic distance field, can determine that ray step and this stratal surface intersect, simultaneously because the starting point in ray step and terminal point are put the symbol normal distance of aspect at the value representation in the symbolic distance field, therefore can be determined the definite position of intersection point by the ratio of distance value.In case after determining intersection point, the ray step will be truncated, intersection point will become new terminal point of ray step.
Reflection and refraction when intersect on ray step and tomography or stratum, will occur in step C4, use the symbolic distance field, can calculate rapidly the normal vector at ray and aspect joining place, with conveniently calculating ray incident angle and emergence angle.When the ray step intersected at tomography, the reflected ray step still carried out in the current tetrahedron of current stratum body, but has new directions of rays; Refracted rays step will carry out in sharing the adjacent tetrahedron of fault surface with current tetrahedron, can determine stratum, new tetrahedron place body according to the rule of judgement spatial point layer position, if not at current stratum body, need to load corresponding stratum body speed file, and use the speed create-rule of new stratum body to go on foot point interpolation as refracted rays to go out velocity amplitude.When the ray step intersected at the stratum, the reflected ray step still carried out in the body of current stratum, but has new directions of rays; Refracted rays step will originate in new stratum body, and new stratum body uses the rule of the judgement spatial point layer position of revising to determine, and loads corresponding stratum body speed file, uses new tetrahedron duplicate to go on foot point interpolation as refracted rays and goes out velocity amplitude.
The rule of the judgement spatial point layer position in step C3 and C4 refers to: begin the selectivity traversal from stratum binary tree root node, gone out the value of given spatial point by interpolation in the symbolic distance field of aspect corresponding to current binary tree node, symbol by this value determines that the child node of next traversal is left child node or right child node, so continue traversal, namely provided layer position, spatial point place until arrive at leaf node; Referring to of the rule of the judgement spatial point layer position of revising: to the nonleaf node that intersects the stratum, if current stratum body is positioned at left subtree, continue traversal from right subtree when traversal, vice versa.
If effective terminal point in ray step is positioned at stratal surface, reflection and refraction will occur in ray.The reflected ray step still originates in current stratum body but has had different directions, and the refracted rays step will originate in new stratum body, and starting point was the terminal point in a upper ray step, and the speed file that needs to load new stratum body calculates.So, how to find the new stratum body at refracted rays step place? directly the rule of use judgement spatial point layer position can run into some molehills, because starting point is located immediately on stratal surface.For this reason, the rule of judgement spatial point layer position is carried out simple modifications: to the nonleaf node that intersects the stratum, if current stratum body is positioned at left subtree, continue traversal from right subtree when traversal, vice versa.
Use the rule of revised judgement spatial point layer position, can determine rapidly the stratum body at refracted rays step place, and load corresponding speed file and calculate.Have dual duplicate owing to intersecting tetrahedron element with aspect, the velocity amplitude of the ground zero in refracted rays step will carry out interpolation according to the tetrahedron duplicate in the body speed file of new stratum.
If the ray step does not intersect with any stratum, need inspection ray step terminal point whether to be marked as and be positioned at tomography, if be positioned at tomography, reflection and refraction also will occur in the ray step so.as previously mentioned, owing to having cancelled tomography both sides tetrahedron topological adjacency relation, the tetrahedron that tomography is shared face keeps a duplicate separately, therefore the reflected ray step still carries out in the current tetrahedron of current stratum body, only had new directions of rays, and refracted rays step will carry out in sharing the tetrahedron of fault surface with current tetrahedron, can determine stratum, new tetrahedron place body according to the rule of judgement spatial point layer position, if not at current stratum body, need to load corresponding stratum body speed file, and use new tetrahedron duplicate to go on foot point interpolation as refracted rays to go out velocity amplitude.
Reflection and refraction can occur when ray runs into the interface, and the Snell law has been described this rule.According to the Snell law, the present invention uses following formula to calculate ray in the parameter at p point place:
V in formula (r) ray is in degree of depth r place velocity of propagation, and i is ray at degree of depth r place incident angle.
The tradition ray-tracing scheme is very difficult often when calculating the exit direction of ray at a certain intersection point place, particularly at the subterranean zone of complex structure.In the present invention, owing to adopting implicit method to express stratal surface, introduced symbolic distance field and stratum binary tree, can facilitate and easily calculate the normal vector at ray and aspect joining place, thereby calculate fast the ray exit direction in step, this is particularly favourable for extensive crossing calculating.And, when the crossing computing of ray and aspect each time, can obtain the parameters such as interfacial curvature matrix such as joining locus, interface normal vector and corresponding topical coordinate, these are all to calculate very important and sufficient information for ray reflection and refraction and kinematics and kinetic parameter.
In sum, provide following embodiment.Fig. 4～Figure 12 shows that ray tracing schematic diagram of the embodiment of the present invention.Fig. 4～Fig. 6 is that ray has the complex geological structure model tracking schematic diagram of multi-fault, and wherein, Fig. 4 is the 3 D complex geological structure model of input; Fig. 5 has shown the concrete situation that the wall scroll ray is propagated in this model, comprise reflection situation and refraction situation; Fig. 6 is for any point source of this model inside the propagation schematic diagram that the fan-shaped divergent-ray arrives a certain section.Fig. 7～Fig. 9 is that ray has the complex geological structure model tracking schematic diagram of mushroom body, and wherein, Fig. 7 is the 3 D complex geological structure model with the mushroom body; Fig. 8 is the propagation schematic diagram that any point source of model inside is the fan-shaped divergent-ray, the situation of the ray reflection that this schematic diagram is only considered; To be ray have on mushroom body border the diagrammatic cross-section of propagating under the situation of high versus speed to Fig. 9.Figure 10～Figure 12 implements the schematic diagram of example ray-tracing scheme of the present invention in the complex geological structure model in actual seismic exploration, wherein, Figure 10 uses the rate pattern corresponding to complex geological structure model of implied expression method construct in the actual seismic exploration; Figure 11 is the propagation schematic diagram that the inner point source arbitrarily of model is the fan-shaped divergent-ray, the situation of the broken line reflection that this schematic diagram is only considered; Figure 12 is the structure wavefront surface schematic diagram in this actual geologic model.
Certainly; the present invention can also have other various embodiments; in the situation that do not deviate from spirit of the present invention and essence thereof; those of ordinary skill in the art work as can make according to the present invention various corresponding changes and distortion, but these corresponding changes and distortion all should belong to the protection domain of the appended claim of the present invention.
Claims (4)
1. interface perception ray tracing method based on implicit model expression, under the three-dimensional geologic structure that implicit method is expressed was model constrained, the desin speed model, carried out ray tracing automatically, is used for the system of seismic exploration data inverting, modeling and skew; It is characterized in that, comprise the following steps:
Steps A, geological interface and geologic body generation method construct three-dimensional geologic structure model based on implicit expression comprise symbolic distance field and stratum binary tree;
Step B, desin speed model automatically on the three-dimensional structure model basis of implied expression;
Step C carries out ray tracing according to tectonic model and rate pattern;
Wherein, the symbolic distance field in steps A is the implied expression form of tectonic model, and the zero contour surface of symbolic distance field is in order to represent corresponding geological stratum; The stratum binary tree is the data structure of describing geological stratum and geologic body spatial topotaxy, and have following character: each non-leaf node represents an aspect, and each leaf node represents the area of space that aspect is cut apart, i.e. stratum body; Main stratum is father or the ancestor node on respective secondary stratum always; To arbitrary stratum binary tree node, its left sibling, oneself, its aspect corresponding to right node that always becomes sequence or reverse sequence arranges.
2. a kind of interface perception ray tracing method based on implicit model expression according to claim 1, is characterized in that, described step B further comprises:
Step B1, be that each stratum body is specified tetrahedral grid according to stratum binary tree and symbolic distance field, these tetrahedral grids comprise and are positioned at the tetrahedral grid that the stratum body is inner and intersect with stratum body outland aspect fully, crossing tetrahedral grid with its Different Strata body that intersects in respectively preserve a duplicate;
Step B2 is the tetrahedral grid summit tax velocity amplitude of each stratum body appointment according to the speed create-rule, and take the stratum body as unit organization speed file; Speed generates rule and is jointly determined by stratal surface base velocity function and velocity gradient, and wherein velocity gradient derives from sand smeller and geophysicist's priori or obtained by migration velocity analysis.
3. a kind of interface perception ray tracing method based on implicit model expression according to claim 1, is characterized in that, described step C further comprises:
Step C1 according to kinematical equation, adopts the method for cutting apart when waiting, and the raypath segmentation is calculated; Each calculates segmentation and is called a ray step, and to each ray step, if itself and the border triangle intersect of tetrahedron element, it will be blocked by this face so, and joining becomes this ray and goes on foot new terminal point;
Step C2, check whether the ray step is crossing with tomography: because fault surface adopts tetrahedral grid face Explicit Expression, and cancelled tomography both sides tetrahedron topological adjacency relation, the adjacent tetrahedron of sharing fault surface keeps a common sides duplicate separately, when ray step and tetrahedral border triangle intersect, check whether this face is fault surface, if, this ray step terminal point is labeled as and is positioned at tomography, namely this ray step intersects with tomography;
Step C3, check whether the ray step is crossing with stratal surface, because stratal surface has adopted the hidden function representation that shows, can not directly judge crossing, can use the rule of judgement spatial point layer position to determine in the binary tree of stratum: if starting point and the terminal point in ray step are positioned at the same leaf node of stratum binary tree, can determine that this ray step is positioned at certain stratum body fully and not crossing with stratal surface; Otherwise this ray step intersects with a stratal surface at least; Further determine the definite position of intersection point: check each node on the ancestral path of tracing back from this ray step stratum, starting point place body leaf node to root node, for nonleaf node, if the starting point in ray step and terminal point be opposite in sign in the symbolic distance field, can determine that ray step and this stratal surface intersect, the ratio of the starting point in ray step and the terminal point value in the symbolic distance field can be determined the definite position of intersection point;
Reflection and refraction when intersect on ray step and tomography or stratum, will occur in step C4, use the symbolic distance field, can calculate rapidly the normal vector at ray and aspect joining place, with conveniently calculating ray incident angle and emergence angle; When the ray step intersected at tomography, the reflected ray step still carried out in the current tetrahedron of current stratum body, but has new directions of rays; Refracted rays step will carry out in sharing the adjacent tetrahedron of fault surface with current tetrahedron, can determine stratum, new tetrahedron place body according to the rule of judgement spatial point layer position, if not at current stratum body, need to load corresponding stratum body speed file, and use the speed create-rule of new stratum body to go on foot point interpolation as refracted rays to go out velocity amplitude; When the ray step intersected at the stratum, the reflected ray step still carried out in the body of current stratum, but has new directions of rays; Refracted rays step will originate in new stratum body, and new stratum body uses the rule of the judgement spatial point layer position of revising to determine, and loads corresponding stratum body speed file, uses new tetrahedron duplicate to go on foot point interpolation as refracted rays and goes out velocity amplitude.
4. a kind of interface perception ray tracing method based on implicit model expression according to claim 3, is characterized in that, in the rule of the judgement spatial point layer position of the rule of the judgement spatial point layer position in described step C3 and C4 and correction:
The rule of judgement spatial point layer position refers to: begin the selectivity traversal from stratum binary tree root node, gone out the value of given spatial point by interpolation in the symbolic distance field of aspect corresponding to current binary tree node, symbol by this value determines that the child node of next traversal is left child node or right child node, so continue traversal, namely provided layer position, spatial point place until arrive at leaf node;
Referring to of the rule of the judgement spatial point layer position of revising: to the nonleaf node that intersects the stratum, if current stratum body is positioned at left subtree, continue traversal from right subtree when traversal, vice versa.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN 201110414879 CN102495427B (en) | 2011-12-10 | 2011-12-10 | Interface perception ray tracing method based on implicit model expression |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN 201110414879 CN102495427B (en) | 2011-12-10 | 2011-12-10 | Interface perception ray tracing method based on implicit model expression |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN102495427A CN102495427A (en) | 2012-06-13 |
| CN102495427B true CN102495427B (en) | 2013-06-19 |
Family
ID=46187267
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN 201110414879 Active CN102495427B (en) | 2011-12-10 | 2011-12-10 | Interface perception ray tracing method based on implicit model expression |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN102495427B (en) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR102042539B1 (en) | 2012-07-24 | 2019-11-08 | 삼성전자주식회사 | Method and apparatus for ray tracing |
| CN104077469B (en) * | 2014-05-28 | 2017-09-01 | 中国人民解放军海军航空工程学院 | Segment iteration remaining time method of estimation based on prediction of speed |
| US10267934B2 (en) * | 2015-01-13 | 2019-04-23 | Chevron U.S.A. Inc. | System and method for generating a depositional sequence volume from seismic data |
| CN104834005B (en) * | 2015-04-28 | 2017-07-07 | 中国石油天然气集团公司 | A kind of method for determining geological interface |
| CN105068133A (en) * | 2015-07-21 | 2015-11-18 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Ray tracking method used for complex geological structure model |
| EP3510564B1 (en) * | 2016-09-06 | 2022-11-02 | Virtamed AG | Ray-tracing methods for realistic interactive ultrasound simulation |
| US11471702B2 (en) * | 2016-12-23 | 2022-10-18 | Koninklijke Philips N.V. | Ray tracing for a detection and avoidance of collisions between radiotherapy devices and patient |
| CN106980139B (en) * | 2017-03-29 | 2018-10-30 | 安徽建筑大学 | A kind of seismic ray method for tracing based on direction vector |
| CN111210522B (en) * | 2020-01-14 | 2021-03-19 | 西南石油大学 | Method for tracking streamline distribution in three-dimensional unstructured grid flow field by using FEM (finite element modeling) |
| CN111968227B (en) * | 2020-06-02 | 2022-09-16 | 中南大学 | Three-dimensional geological fault network uncertainty analysis method, system and storage medium |
| CN113626894B (en) * | 2021-08-27 | 2024-01-19 | 北京航空航天大学 | Entity motion interface tracking method based on composite implicit boundary |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1997011390A2 (en) * | 1995-09-19 | 1997-03-27 | Exxon Production Research Company | Multiple suppression in geophysical data |
| CN1712991A (en) * | 2004-06-25 | 2005-12-28 | 中国石油化工股份有限公司 | Ray traction in earthquake prospection |
| CN101630014A (en) * | 2008-07-16 | 2010-01-20 | 中国石油天然气集团公司 | Method for imaging anisotropic medium through utilization of vertical seismic profile data |
-
2011
- 2011-12-10 CN CN 201110414879 patent/CN102495427B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO1997011390A2 (en) * | 1995-09-19 | 1997-03-27 | Exxon Production Research Company | Multiple suppression in geophysical data |
| CN1712991A (en) * | 2004-06-25 | 2005-12-28 | 中国石油化工股份有限公司 | Ray traction in earthquake prospection |
| CN101630014A (en) * | 2008-07-16 | 2010-01-20 | 中国石油天然气集团公司 | Method for imaging anisotropic medium through utilization of vertical seismic profile data |
Also Published As
| Publication number | Publication date |
|---|---|
| CN102495427A (en) | 2012-06-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102495427B (en) | Interface perception ray tracing method based on implicit model expression | |
| CN101303414B (en) | Method for generating ground layer surface and geologic body based on level set | |
| Jacquemyn et al. | Surface-based geological reservoir modelling using grid-free NURBS curves and surfaces | |
| EP3293552B1 (en) | System and method for editing geological models by switching between volume-based models and surface-based structural models augmented with stratigraphic fiber bundles | |
| CN101582173B (en) | Block model building method for complex geological structure | |
| CN102867330B (en) | Region-division-based spatial complex horizon reconstruction method | |
| CN104635262B (en) | A kind of positive reversed fault isopleth automatic generation method based on extended rectangular grid | |
| CN106934860A (en) | A kind of three-dimensional geological modeling method based on T battens | |
| CN106981093A (en) | A kind of three-dimensional formation parallel modeling method of subregion constraint coupling | |
| CN103413297A (en) | Cutting method based on integrated three-dimensional GIS model | |
| CN101533102B (en) | Triangular mesh ray tracing global method of two-dimensional complex construction | |
| CN103245972B (en) | A kind of method determining complex geological structure in two-dimensional space | |
| CN102194252A (en) | Geological-stratum-structure-based method for generating triangular lattice grids | |
| CN112394404B (en) | Progressive reservoir fine characterization method | |
| CN102867332B (en) | Based on the multistage subdivided meshes curved surface fitting method of complex boundary constraint | |
| CN101915088A (en) | Method and device for generating oil migration path | |
| Wu et al. | A 3D modeling approach to complex faults with multi-source data | |
| CN103698810A (en) | Hybrid network minimum travel time ray tracing tomography method | |
| US20210396897A1 (en) | Computer implemented method for correcting a reservoir model of a reservoir geological formation based on seismic images | |
| CN103543478A (en) | Geologic morphological interpolation KM (Kriging and Multiple-point geostatistics) method | |
| CN106814392B (en) | The velocity modeling method in three-dimensional secondary closely plast fusion face | |
| CN106846481B (en) | Geological profile generation method | |
| Gao et al. | Mathematical Interpolation and Correction of Three-Dimensional Modelling of High-Speed Railway. | |
| Shen et al. | Three-dimensional modeling of loose layers based on stratum development law | |
| CN1825139B (en) | Land bedding computer graph generating method based on deforming field |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20170721 Address after: 100191 Beijing City, Haidian District Zhichun Road Jinqiu International Building No. 6 A block 1501 Patentee after: Beijing grid world software technology Limited by Share Ltd Address before: 100191 Haidian District, Xueyuan Road, No. 37, Patentee before: Beihang University |