CN102495427A - Interface perception ray tracing method based on implicit model expression - Google Patents

Abstract

The invention discloses an interface-aware ray tracing method based on implicit model expression, which is used in a seismic data inversion, modeling and migration system. The steps include: expressing the geological interface based on an implicit function; using a binary tree to describe stratum deposition sequence, combined with the implicitly expressed geological interface to realize the generation of the 3D geological body structure model; on the basis of the implicitly expressed structural model, according to the specified velocity distribution rules of each stratum body, Velocity values are assigned to the vertices of the tetrahedral unit intersecting layers, and the velocity files are organized and stored in units of formation volumes; under the constraints of the implicitly expressed structural model, according to kinematic equations and reflection laws, ray tracing based on interface perception is used in the construction The model is iteratively performed in sections. The invention can automatically carry out velocity modeling of complex geological structure models, complete high-precision and robust ray tracing, and obtain important parameters for seismic inversion, migration imaging and numerical simulation.

Description

An interface-aware ray tracing method based on implicit model representation

technical field

The invention relates to a system and method for automatically performing velocity modeling and ray tracing under the constraints of an implicitly expressed geological structure model, in particular to an interface-aware ray tracing method based on implicit model expression.

Background technique

Because rays carry a large amount of information that is useful for seismic inversion, modeling, and migration, ray tracing methods play a very important role in seismic exploration data processing. If the subsurface velocity is known, parameters such as ray path, travel time, amplitude, and exit angle can all be calculated numerically by ray tracing.

As early as 1977, Cerveny and others proposed the ray tracing method, and later Cerveny and Hron supplemented and perfected it, which constituted the classic ray tracing theory. Ray tracing is also very direct from mathematical theory to computer implementation. However, until recent years, the application of 3D ray tracing in actual seismic inversion and modeling is also very limited. Processing methods generate 3D synthetic seismic records, or test the plausibility of certain assumptions and approximations. An important reason why this method cannot be widely used in the petroleum industry is that the application of 3D seismic exploration in the real complex geological structure area is very computationally resource-intensive, and the error of the traditional method is also very large. Many grid-based velocity modeling methods have been used in ray tracing and tomographic inversion at home and abroad. In 2004, Chapman also introduced a ray-based velocity modeling algorithm that can perceive the layer frame and propagate physical properties along the ray when it encounters a layer. However, these traditional methods are usually based on tetrahedral grids, and the generation of grids is restricted by sections and layers, including unconformity and overlying strata. In areas with complex underground structures, grid generation is often very difficult and the quality of the generated grids is not good. high. In 2006, Ruger et al. proposed a method to automatically generate formation velocity structure grids from finely sampled uniform velocity grids. This method does not need to provide additional constraint information (such as faults and formations), and the velocity grids in discontinuous regions are embedded in Delaunay subdivided tetrahedral meshes, but they do not fundamentally solve the problem of constrained mesh generation. Traditional ray tracing usually divides the workflow into several different stages: picking up structural features, generating velocity models, and ray tracing. Except that the output of the previous stage will be used as the input of the next stage, each stage is often independent of each other. Necessary to contact.

The invention discloses an interface-aware ray tracing method based on implicit model expression. The three main stages of ray tracing: structure modeling, velocity modeling and ray tracing have a unified method basis—implicit expression method, which can not only The complex geological structure is effectively expressed, and each stage is organically connected. In addition to the input of the previous stage, it also provides structural support and calculation method support for the latter stage, forming a unified whole.

Contents of the invention

The technical problem to be solved by the present invention is to study the method of automatic velocity modeling and high-precision ray tracing in the complex geological structure model, and the system for seismic exploration data inversion, modeling and migration to solve the problem of existing ray The tracking method is difficult to implement in complex geological structures, so as to overcome the difficulties of low precision, low efficiency and difficulty in effectively supporting seismic inversion, migration and simulation in the existing technology.

The technical solution adopted by the present invention to solve the technical problem is to provide an interface-aware ray tracing method based on implicit model expression, which performs velocity modeling and interface-aware ray tracing under the constraints of the structural model expressed by the implicit model, Include the following steps:

Step A, constructing a three-dimensional geological structure model based on an implicit geological interface and geological body generation method, including a signed distance field and a stratum binary tree;

Step B, automatically constructing a velocity model based on the implicitly expressed 3D construction model;

Step C, performing ray tracing according to the construction model and the velocity model.

The symbolic distance field in step A is the implicit expression form of the structural model, and the zero isosurface of the symbolic distance field is used to represent the corresponding geological layer; the stratum binary tree is a data structure describing the spatial topological relationship between the geological layer and the geological body, which has the following Properties: each non-leaf node represents a layer, and each leaf node represents a spatial region divided by the layer, that is, the stratum body; the main stratum is always the parent or ancestor node of the corresponding auxiliary stratum; for any stratum binary tree node, its The level corresponding to the left node, itself, and its right node always becomes the sub-arrangement of the sequence or reverse sequence.

Step B further includes:

Step B1, specify tetrahedral grids for each stratigraphic volume according to the stratigraphic binary tree and the signed distance field, these tetrahedral grids include tetrahedral grids completely inside the stratum volume and intersecting with the boundary layer of the stratum volume, the intersecting tetrahedron A copy of the grid is kept in each of the different stratigraphic volumes it intersects;

Step B2, according to the velocity generation rules, assign velocity values to the vertices of tetrahedral grids specified by each stratum body, and organize velocity files in units of stratum bodies; the velocity generation rules are jointly determined by the stratum base velocity function and velocity gradient, where the velocity Gradients are derived from prior knowledge of geologists and geophysicists or obtained from migration velocity analysis.

Step C further includes:

Step C1, according to the kinematic equation, adopt the method of isochronous division to calculate the ray path in segments; each calculation segment is called a ray step, and for each ray step, if it intersects with the boundary triangle of the tetrahedron unit, It will then be truncated by that face, and the point of intersection becomes the new end point of the raystep.

Step C2, check whether the ray step intersects with the fault: Since the fault plane is expressed by a tetrahedral grid surface, and the topological adjacency relationship of tetrahedrons on both sides of the fault plane is canceled, the adjacent tetrahedrons sharing the fault plane each retain a common surface complex This, when the ray step intersects the boundary triangle of the tetrahedron, check whether the surface is a fault plane, if so, mark the end point of the ray step as being located on the fault, that is, the ray step intersects the fault;

Step C3, check whether the ray step intersects with the stratum layer. Since the stratum layer is expressed by an implicit function, the intersection cannot be directly determined. The rule for judging the spatial point horizon can be determined in the stratum binary tree: if the starting point and end point of the ray step are at If the same leaf node of the stratum binary tree is used, it can be determined that the ray step is completely located in a certain stratum body without intersecting the stratum; otherwise, the ray step intersects at least one stratum. Further determine the exact position of the intersection point: check each node on the ancestral path from the leaf node of the stratum where the starting point of the ray step is located to the root node, for non-leaf nodes, if the starting point and the ending point of the ray step have opposite signs in the signed distance , then it can be determined that the ray step intersects with the formation layer, and the exact position of the intersection point can be determined by the ratio of the value of the start point and end point of the ray step in the signed distance field.

Step C4, when the ray step intersects with the fault or formation, reflection and refraction will occur, and the normal vector at the intersection point of the ray and the layer can be quickly calculated by using the signed distance field, so as to facilitate the calculation of the incident angle and exit angle of the ray. When a raystep intersects a fault, a reflected raystep will still be performed in the current tetrahedron of the current formation volume, but with a new ray direction; a refracted raystep will be performed in an adjacent tetrahedron that shares the fault plane with the current tetrahedron, According to the rules for judging the horizon of the space point, the formation body where the new tetrahedron is located can be determined. If it is not in the current formation body, the corresponding formation body velocity file needs to be loaded, and the velocity value of the starting point of the refracted ray step is interpolated using the velocity generation rule of the new formation body. When a raystep intersects a stratum, a reflected raystep is still performed in the current stratum volume, but with a new ray direction; a refracted raystep will start in a new stratum volume, which uses the modified judgment space point horizon The rules are determined, and the velocity file of the corresponding layer is loaded, and the new tetrahedron copy is used to interpolate the velocity value for the starting point of the refraction ray step.

The rules for judging the spatial point horizon in steps C3 and C4 refer to: start selective traversal from the root node of the stratum binary tree, interpolate the value of a given spatial point from the signed distance field of the level corresponding to the current binary tree node, and use the value The sign of determines whether the next traversed child node is the left child node or the right child node, so continue traversing until it reaches the leaf node, which gives the level of the spatial point; the revised rule for judging the level of the spatial point refers to: when Traverse to non-leaf nodes of intersecting strata, if the current stratum is in the left subtree, continue traversing from the right subtree, and vice versa.

Compared with the prior art, the present invention has the following advantages:

1. The present invention can directly serve velocity modeling and ray tracing due to the implicit method used to express the structure model, and is suitable for complex geological structure environments;

2. In the velocity modeling of the present invention, the multi-repetition method is adopted to generate velocity files with formation bodies as units, which can embed the knowledge of geologists and geophysicists on geological laws, and can realize the discontinuous expression of the velocity at the formation interface. Strong support for ray tracing;

3. The present invention utilizes the signed distance field and the stratum binary tree to form a set of reasonable ray intersection, reflection and refraction rules, which effectively improves the accuracy and efficiency of ray tracing.

Description of drawings

Figure 1 is a geological structure section and the corresponding stratum binary tree generated;

Fig. 2 is a schematic diagram of a tetrahedron intersecting with the ground layer;

Fig. 3 is a schematic representation of a velocity model with double duplicates;

Figure 4 is a complex geological structure model with multiple faults;

Figure 5 is a schematic diagram of the reflection and refraction of a single ray in a complex structural model;

Fig. 6 is a schematic diagram of propagation of any internal point source ray to a section in the complex structural model;

Figure 7 is a complex geological structure model with a mushroom body;

Fig. 8 is a schematic diagram of the propagation (only reflection) of an internal arbitrary point source ray in a complex geological structure model with a mushroom body;

Fig. 9 is a schematic cross-sectional view of rays propagating in a structural model with high contrast velocity at the boundary of the mushroom body;

Fig. 10 is the velocity model corresponding to the actual seismic exploration geological structure model;

Fig. 11 is a schematic diagram of propagation (only considering refraction) of internal arbitrary point source rays in the actual seismic exploration geological structure model;

Fig. 12 is a schematic diagram of a structural wave front propagating in an actual seismic exploration geological structure model by an internal arbitrary point source ray.

Detailed ways

The specific implementation manner of the present invention will be described below in conjunction with the accompanying drawings.

An interface-aware ray tracing method based on implicit model representation consists of three steps, including:

Step A, constructing a three-dimensional geological structure model based on an implicit geological interface and geological body generation method, including a signed distance field and a stratum binary tree;

Among them, based on the implicit geological interface and geological body generation method in step A, "a level set-based formation layer and geological body generation method" (has obtained national patent, patent authorization number is ZL200810112263.6), also Other mathematical methods can be used to generate corresponding symbolic data fields for each level (such as Frank, T., Tertois, A.L., Mallet, J.L., 2007. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data. Computer & Geosciences 33: 932 -943; Calcagno, P., Chiles, J., Courrioux, G., Guillen, A., 2008, Geological modeling from field data and geological knowledge-part I: Modeling method coupling 3D potential-field interpolation and geological rules: Physics of the Earth and Planetary Interiors 171(1-4), 147-157.) plane as the zero isosurface of the symbolic data field.

The implicit expression method uses the signed distance field as the representation form, and uses the zero isosurface of the signed distance field to express the implicit surface. In the present invention, a signed distance field is constructed for each stratigraphic layer, and the zero isosurface of the signed distance field is used to express the geological layer. At the same time, a data structure describing the spatial topological relationship between the geological level and the geological body is defined—stratum binary tree, which has the following properties: each non-leaf node represents a level, and each leaf node represents a spatial region divided by the level. That is, the stratum body; the main stratum is always the parent or ancestor node of the corresponding auxiliary stratum; for any stratum binary tree node, its left node, itself, and the level corresponding to its right node always become the sub-arrangement of the sequence or reverse sequence . Fig. 1 is a schematic diagram of a geological structure profile and its corresponding stratum binary tree, wherein Si (i=1, 2, 3, 4) represents a stratum layer and corresponds to the upper surface of a stratum body, and Ti (i=1, 2, ..., 9) represent the spatial sub-regions divided by strata. In the figure, the stratum S2 divides the spatial region into two sub-regions T1 and T2, and the two sub-regions are divided respectively.

Step B, automatically constructing a velocity model based on the implicitly expressed 3D construction model;

Step B1, specify tetrahedral grids for each stratigraphic volume according to the stratigraphic binary tree and the signed distance field, these tetrahedral grids include tetrahedral grids completely inside the stratum volume and intersecting with the boundary layer of the stratum volume, the intersecting tetrahedron A copy of the grid is kept in each of the different stratigraphic volumes it intersects;

Step B2, according to the velocity generation rules, assign velocity values to the vertices of tetrahedral grids specified by each stratum body, and organize velocity files in units of stratum bodies; the velocity generation rules are jointly determined by the stratum base velocity function and velocity gradient, where the velocity Gradients are derived from prior knowledge of geologists and geophysicists or obtained from migration velocity analysis.

Velocity modeling is essentially a kind of attribute modeling. When the velocity model is applied to interface-aware ray tracing, the accuracy of the velocity model is very high, especially the velocity expression at the interface of faults and formations. Many mesh-based velocity modeling methods have been used for slice-aware ray tracing, but the main difficulty of these mesh-based interface-aware velocity modeling algorithms lies in the mesh generation under the constraints of complex structural models. A major challenge of velocity modeling based on implicit representations is that implicit formations are not tetrahedral mesh interfaces, and no mesh vertices correspond to formation layers. Finding which strata each mesh vertex belongs to and assigning them values, by enumerating or comparing with all strata, is time-consuming and clumsy. For meshes that intersect layers, there are different solutions for how to assign velocities to each vertex. Because not all four vertices are from the same formation volume, assigning velocities to interior points by direct linear interpolation can lead to erroneous results. In 2007, Bargteil et al. proposed a method to solve a similar problem in geological body reconstruction. This method sets the tetrahedron that is completely above or below the formation as a constant, and the tetrahedron that intersects with the layer is based on two materials The ratio between adopts a volume approximation strategy. This is a simple and effective approach for plastically deforming solid materials, but is not suitable for velocity modeling for ray tracing. Because in the process of ray tracing, there are a large number of intersecting calculations between rays and slices, the inaccuracy of the velocity of adjacent slices will affect the calculation of the ray direction.

In isotropic media, velocities can be defined with constant gradients, exponential functions, conic functions, or analyzable features. Due to the sedimentary characteristics of the formation, the geological structure has the symmetrical characteristics of lateral isotropy, and the velocity may no longer be a constant gradient in the vertical depth direction. However, in dynamic tectonic regions, strata may tilt after long-term tectonic evolution, resulting in a tilt of the transverse axis of symmetry. Most seismic exploration studies are based on layered geological models, where the velocity within a layer varies linearly due to sedimentary properties. However, the change in velocity may not always be continuous with increasing vertical depth, especially when encountering formations. Therefore, it is necessary to construct a velocity model that conforms to the understanding of geosciences, so that the velocity changes continuously within the layer, and discontinuity may occur at the junction of layers.

Based on the above geological knowledge, the velocity model used for ray tracing and inversion must be consistent with the distribution of formation medium. Therefore, velocity modeling mainly includes two parts: specifying tetrahedral grids for formation bodies, and assigning velocity values to grid vertices according to the principle of inclined transverse isotropy. The present invention adopts a velocity gradient method to calculate the velocity of the grid vertices in each formation volume. First, specify a 2D base velocity function v(x, y) for each subsurface of the formation body, according to the base velocity value and velocity gradient on the formation interface (usually the prior knowledge of geologists and geophysicists or by Migration velocity analysis), use the following formula to calculate the velocity in the formation volume:

v(P)=v ₀ (P)+g*h(P)

In the formula, v(P) represents the velocity value of the vertex P of the tetrahedral grid required to be taken, v ₀ (P) is the base velocity value of the specified formation interface, g is the velocity gradient of the formation body, and h(P) represents the vertex P The normal distance to the formation. Apparently, h(P) can be directly calculated from the signed distance field of the stratigraphic structure model, and also traced from P to the formation layer according to the negative gradient direction of the signed distance field, and the formation base velocity v ₀ (P) can be obtained.

The velocity model is subdivided into a tetrahedral mesh, and during ray tracing, the velocity value of each mesh vertex that the ray passes through is used. The velocities of points inside the tetrahedron are linearly interpolated from the velocities of the four vertices, but for tetrahedral meshes intersecting layers, this interpolation method is no longer appropriate. Because the formation interface is not explicitly defined as the constraint information of the tetrahedral mesh, and the formation layer expressed by the implicit method cannot be expressed by the triangular mesh surface of the tetrahedron. Therefore, velocity values where the grid intersects the slice may be incorrect or uninterpretable if only simple interpolation methods are used. As shown in Figure 2a, the tetrahedral grid unit intersects the formation, vertex A is above the formation, and the other three vertices B, C, and D are below the formation. Now the velocity value of the inner point E is calculated. If you directly use the velocity values of the four vertices for linear interpolation, the resulting velocity value at point E may be incorrect or unexplainable, because the velocity values at point A and points B, C, and D are calculated using different velocity gradients . Wrong velocity values will lead to wrong ray tracing, and these formation boundaries are extremely sensitive areas for ray tracing, which is very important for ray tracing, and it will determine the direction of the ray.

The present invention adopts the method of double replication to solve this problem. In order to calculate the velocity value of the point in the intersecting grid unit above the layer, the formation body velocity generation rule above the layer is extended to the intersecting grid unit. All vertices (including above and below the layer), and the velocity values of these vertices are recorded in the formation body velocity file on the layer. For the vertices of the intersecting grid unit located below the layer, use the same method to extend the formation body velocity generation rule below the layer to all vertices of the intersecting grid unit, and record the velocity values of these vertices under the layer in the velocity file for the formation body. For the more complex situation where multiple layers intersect with grid cells, as shown in Figure 2b, this method can also be easily extended to a multi-replica method, generating velocity files with different replicas for each stratigraphic body, that is, if the tetrahedron If it intersects with multiple formation bodies, there will be multiple copies of the tetrahedron with different velocity values, which will be recorded in the velocity file of the corresponding formation body. Figure 3 describes the velocity assignment method.

Therefore, the velocity modeling process also needs to be extended accordingly. For each formation volume, it is necessary to find the tetrahedron grids that are located in the formation volume and intersect with its formation plane, and use the velocity generation rules of the formation volume to assign values to the vertices of the tetrahedron grids. The main difficulty in implementation is how to find these tetrahedra. The specific method is: for each stratum body, traverse all tetrahedron grids, use the signed distance field and the stratum binary tree to determine whether the tetrahedron is contained in the stratum body or intersects with the stratum layer, usually, the search path in the stratum binary tree is from Indicates the path from the leaf nodes of the stratigraphic volume back to the root node. For example, to check whether a tetrahedron C is contained in or intersects with a stratum volume S, the traversal process in the stratum binary tree starts from the leaf node representing the stratum volume S, and an iterator is defined to indicate its traversal process. First, the iterator points to the parent node of the leaf node (the non-leaf node in the stratum binary tree represents the geological layer), if the stratum S is on one side of the stratum indicated by the iterator, and all vertices of the tetrahedron C are on the other side, then It can be determined that the tetrahedron is not included in or intersects with the formation body S, and the traversal is terminated; otherwise, the position of the iterator is updated, and traced along the ancestor path until the root node. If the tracing process has not violated the judgment condition, the tetrahedron S can be determined Contained in or intersected with the formation volume S. After finding all the tetrahedrons contained in or intersecting with the formation volume S, the velocity value is assigned to the formation volume according to the velocity generation rules, and its corresponding velocity file is generated.

Step C, performing ray tracing according to the construction model and the velocity model.

Step C1, according to the kinematic equation, adopt the method of isochronous division to calculate the ray path in segments; each calculation segment is called a ray step, and for each ray step, if it intersects with the boundary triangle of the tetrahedron unit, It will then be truncated by that face, and the point of intersection becomes the new end point of the raystep.

The kinematic equation of ray tracing belongs to the equation of the equation, and Cerveny uses the characteristic method to describe its specific form, as follows:

du=v ⁿ dτ

$\frac{{\mathrm{<span\; class="notranslate">\; dx</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; du</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; v</span>}}^{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}}}$

$\frac{{\mathrm{<span\; class="notranslate">\; dp</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; du</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{{\mathrm{<span\; class="notranslate">\; v</span>}}^{\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}}\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; v</span>}}{{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}$

From the above equations, it can be seen that along the ray trajectory x _i and its distribution p _i is a monotonically increasing function of the independent variable u. x _i (u) is called the ray path, and the iterative calculation formula of the ray path can be obtained by using this set of equations. Ray tracing starts at a source point, and the analog wavefront construction method fires a set of rays approximately distributed at constant incremental angles. Whenever topologically adjacent rays are separated by more than a certain distance, a new ray is inserted to maintain the wavefront shape. The ray path is iteratively calculated by isochronous segmentation, and each calculation segment is called a ray step. For each ray step, if it intersects with a face of the tetrahedron element, it will be truncated by the face, and the intersection point becomes the new end point of the ray step. At the same time, it is necessary to mark the intersection point to see if the intersecting surface is located on a fault, because the ray will be refracted or reflected at the fault.

Step C2, check whether the ray step intersects with the fault: Since the fault plane is expressed by a tetrahedral grid surface, and the topological adjacency relationship of tetrahedrons on both sides of the fault plane is canceled, the adjacent tetrahedrons sharing the fault plane each retain a common surface complex This, when the ray step intersects the boundary triangle of the tetrahedron, check whether the surface is a fault plane, if so, mark the end point of the ray step as being located on the fault, that is, the ray step intersects the fault;

Step C3, check whether the ray step intersects with the stratum layer. Since the stratum layer is expressed by an implicit function, the intersection cannot be directly determined. The rule for judging the spatial point horizon can be determined in the stratum binary tree: if the starting point and end point of the ray step are at If the same leaf node of the stratum binary tree is used, it can be determined that the ray step is completely located in a certain stratum body without intersecting the stratum; otherwise, the ray step intersects at least one stratum. Further determine the exact position of the intersection point: check each node on the ancestral path from the leaf node of the stratum where the starting point of the ray step is located to the root node, for non-leaf nodes, if the starting point and the ending point of the ray step have opposite signs in the signed distance , then it can be determined that the ray step intersects with the formation layer, and the exact position of the intersection point can be determined by the ratio of the value of the start point and end point of the ray step in the signed distance field.

Since the ground plane is not explicitly expressed as a grid surface, it is necessary to further check whether the ray step intersects with a certain ground plane, and the rules for judging the spatial point horizon can help realize this. If the starting point and end point of the ray step are located at the same leaf node of the stratum binary tree, it can be determined that the ray step is completely located in a stratum body without intersecting the stratum; otherwise, the ray step intersects at least one stratum. In order to determine the exact position of the intersection point, it is first necessary to find out the position of the end point of the ray step. Therefore, in the stratum binary tree, each node on the ancestral path from the leaf node of the stratum body where the starting point of the ray step is located to the root node may become the object to be checked. For each non-leaf node (ground layer), if the start and end points of the ray step are in the opposite sign in the signed distance field, it can be determined that the ray step intersects the layer, and because the start point and the end point of the ray step are in the signed distance field The value of represents the signed normal distance from the point to the slice, so the exact location of the intersection can be determined from the ratio of the distance values. Once the intersection point is determined, the ray step will be truncated, and the intersection point will become the new end point of the ray step.

Step C4, when the ray step intersects with the fault or formation, reflection and refraction will occur, and the normal vector at the intersection point of the ray and the layer can be quickly calculated by using the signed distance field, so as to facilitate the calculation of the incident angle and exit angle of the ray. When a raystep intersects a fault, a reflected raystep will still be performed in the current tetrahedron of the current formation volume, but with a new ray direction; a refracted raystep will be performed in an adjacent tetrahedron that shares the fault plane with the current tetrahedron, According to the rules for judging the horizon of the space point, the formation body where the new tetrahedron is located can be determined. If it is not in the current formation body, the corresponding formation body velocity file needs to be loaded, and the velocity value of the starting point of the refracted ray step is interpolated using the velocity generation rule of the new formation body. When a raystep intersects a stratum, a reflected raystep is still performed in the current stratum volume, but with a new ray direction; a refracted raystep will start in a new stratum volume, which uses the modified judgment space point horizon The rules are determined, and the velocity file of the corresponding layer is loaded, and the new tetrahedron copy is used to interpolate the velocity value for the starting point of the refraction ray step.

The rules for judging the spatial point horizon in steps C3 and C4 refer to: start selective traversal from the root node of the stratum binary tree, interpolate the value of a given spatial point from the signed distance field of the level corresponding to the current binary tree node, and use the value The sign of determines whether the next traversed child node is the left child node or the right child node, so continue traversing until it reaches the leaf node, which gives the level of the spatial point; the revised rule for judging the level of the spatial point refers to: when Traverse to non-leaf nodes of intersecting strata, if the current stratum is in the left subtree, continue traversing from the right subtree, and vice versa.

Rays will be reflected and refracted if the effective endpoint of the ray step is at the ground plane. The reflection ray step still starts from the current formation body but has a different direction, while the refraction ray step will start from the new formation body, the starting point is the end point of the previous ray step, and the velocity file of the new formation body needs to be loaded for calculation. So, how to discover the new stratigraphic body where the refracted raystep is located? There are some minor difficulties in directly using the rules for determining the horizon of a spatial point, since the starting point is located directly on the ground plane. Therefore, a simple modification is made to the rules for judging the horizon of spatial points: when traversing to non-leaf nodes of intersecting strata, if the current stratum is located in the left subtree, continue traversing from the right subtree, and vice versa.

Using the revised rules for judging the horizon of spatial points, the formation body where the refracted ray step is located can be quickly determined, and the corresponding velocity file can be loaded for calculation. Since the tetrahedral cells intersecting the bedding have double duplicates, the velocity values for the new start point of the refracted ray step will be interpolated according to the tetrahedral duplicates in the new formation volume velocity file.

If the raystep does not intersect any formations, it needs to check whether the raystep end point is marked as lying on a fault, if it is on a fault, then the raystep will also be reflected and refracted. As mentioned above, since the topological adjacency relationship of tetrahedrons on both sides of the fault is cancelled, each tetrahedron on the shared plane of the fault retains a copy, so the reflected ray step is still performed in the current tetrahedron of the current formation body, but with The new ray direction, and the refraction ray step will be performed in the tetrahedron that shares the fault plane with the current tetrahedron. According to the rules for judging the horizon of the space point, the layer body where the new tetrahedron is located can be determined. If it is not in the current layer body, it needs to be loaded Take the corresponding formation volume velocity file and use the new tetrahedron replica to interpolate the velocity values for the refraction ray step start point.

Snell's law describes the reflection and refraction of rays when they encounter an interface. According to Snell's law, the present invention uses the following formula to calculate the parameters of the ray at point p:

$\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; sin</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span>}$

where v(r) is the propagation velocity of the ray at the depth r, and i is the incident angle of the ray at the depth r.

Traditional ray tracing methods are often very difficult to calculate the outgoing direction of rays at a certain intersection point, especially in underground areas with complex structures. In the present invention, due to the use of an implicit method to express the stratum layer, the signed distance field and the stratum binary tree are introduced, and the normal vector at the intersection point of the ray and the layer can be calculated conveniently and easily, thereby quickly calculating the outgoing direction of the ray step, which is It is especially beneficial for large-scale intersection calculations. In addition, when each ray intersects with a layer, parameters such as the spatial position of the intersection point, the interface normal vector, and the interface curvature matrix of the corresponding local coordinates can be obtained, which are for ray reflection and refraction as well as kinematic and dynamic parameters Calculate very important and sufficient information.

In summary, the following examples are given. 4 to 12 are schematic diagrams of ray tracing according to the embodiment of the present invention. Figures 4 to 6 are schematic diagrams of ray tracing in a complex geological structure model with multiple faults, in which Figure 4 is the input 3D complex geological structure model; Figure 5 shows the specific situation of a single ray propagating in the model, including the reflection situation and refraction; Figure 6 is a schematic diagram of the propagation of rays emitted by any point source in the model in a fan shape to a certain section. Figures 7 to 9 are schematic diagrams of ray tracing in a complex geological structure model with a mushroom body, in which Figure 7 is a three-dimensional complex geological structure model with a mushroom body; Figure 8 is the propagation of rays emitted by any point source inside the model in a fan shape Schematic diagram, which only considers the situation of ray reflection; FIG. 9 is a schematic cross-sectional view of the ray propagating at the boundary of the mushroom body with high contrast velocity. Figures 10 to 12 are schematic diagrams of implementing the ray tracing method of the present invention in a complex geological structure model in actual seismic exploration, wherein Figure 10 is the velocity corresponding to the complex geological structure model constructed using the implicit expression method in actual seismic exploration Model; Fig. 11 is a schematic diagram of propagation of rays emitted by an arbitrary point source in the model in a fan shape, which only considers the situation of broken line reflection; Fig. 12 is a schematic diagram of the structural wave front in the actual geological model.

Certainly, the present invention can also have other multiple embodiments, without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and deformations according to the present invention, but these corresponding Changes and deformations should belong to the scope of protection of the appended claims of the present invention.

Claims (6)

1. An interface-aware ray tracing method based on implicit model expression. Under the constraint of the 3D geological structure model expressed by the implicit method, the velocity model is automatically constructed and ray tracing is performed, which is used for seismic exploration data inversion, modeling and biasing. The system of moving; It is characterized in that, comprises the following steps: Step A, constructing a three-dimensional geological structure model based on an implicit geological interface and geological body generation method, including a signed distance field and a stratum binary tree; Step B, automatically constructing a velocity model based on the implicitly expressed 3D construction model; Step C, performing ray tracing according to the construction model and the velocity model. 2. A kind of interface-aware ray tracing method based on implicit model expression according to claim 1, characterized in that, the signed distance field in the step A is an implicit expression form of the construction model, and the zero of the signed distance field Isosurfaces are used to represent corresponding geological layers. 3. A method of interface-aware ray tracing based on implicit model expression according to claim 1, characterized in that the stratum binary tree in the step A is a data structure describing the spatial topological relationship between geological layers and geological bodies, with The following properties: each non-leaf node represents a level, and each leaf node represents a space region divided by the level, that is, the stratum body; the main stratum is always the parent or ancestor node of the corresponding auxiliary stratum; for any stratum binary tree node, it The levels corresponding to the left node of , itself, and its right nodes always become sub-arrangements of the sequence or reverse sequence. 4. A method of interface-aware ray tracing based on implicit model expression according to claim 1, wherein the step B further comprises: Step B1, specify tetrahedral grids for each stratigraphic volume according to the stratigraphic binary tree and the signed distance field, these tetrahedral grids include tetrahedral grids completely inside the stratum volume and intersecting with the boundary layer of the stratum volume, the intersecting tetrahedron A copy of the grid is kept in each of the different stratigraphic volumes it intersects; Step B2, according to the velocity generation rules, assign velocity values to the vertices of tetrahedral grids specified by each stratum body, and organize velocity files in units of stratum bodies; the velocity generation rules are jointly determined by the stratum base velocity function and velocity gradient, where the velocity Gradients are derived from prior knowledge of geologists and geophysicists or obtained from migration velocity analysis. 5. A method of interface-aware ray tracing based on implicit model expression according to claim 1, wherein said step C further comprises: Step C1, according to the kinematic equation, adopt the method of isochronous division to calculate the ray path in segments; each calculation segment is called a ray step, and for each ray step, if it intersects with the boundary triangle of the tetrahedron unit, Then it will be truncated by the face, and the intersection point becomes the new end point of the raystep; Step C2, check whether the ray step intersects with the fault: Since the fault plane is expressed by a tetrahedral grid surface, and the topological adjacency relationship of tetrahedrons on both sides of the fault plane is canceled, the adjacent tetrahedrons sharing the fault plane each retain a common surface complex This, when the ray step intersects the boundary triangle of the tetrahedron, check whether the surface is a fault plane, if so, mark the end point of the ray step as being located on the fault, that is, the ray step intersects the fault; Step C3, check whether the ray step intersects with the stratum layer. Since the stratum layer is expressed by an implicit function, the intersection cannot be directly determined. The rule for judging the spatial point horizon can be determined in the stratum binary tree: if the starting point and end point of the ray step are at The same leaf node of the stratum binary tree, it can be determined that the ray step is completely located in a stratum body without intersecting the stratum layer; otherwise, the ray step intersects at least one stratum layer; further determine the exact position of the intersection point: check the stratum where the starting point of the ray step is located For each node on the ancestral path from the leaf node to the root node, for non-leaf nodes, if the start and end points of the ray step have opposite signs in the signed distance field, it can be determined that the ray step intersects the stratum layer, and the ray step’s The ratio of the start and end values in the signed distance field determines the exact location of the intersection point; Step C4, when the ray step intersects with the fault or formation, reflection and refraction will occur, and the normal vector at the intersection point of the ray and the layer can be quickly calculated by using the symbolic distance field, so as to facilitate the calculation of the incident angle and the outgoing angle of the ray; when the ray When the step intersects the fault, the reflected ray step will still be performed in the current tetrahedron of the current formation body, but with a new ray direction; the refracted ray step will be performed in the adjacent tetrahedron that shares the fault plane with the current tetrahedron, according to the judgment The rules of the spatial point horizon can determine the formation body where the new tetrahedron is located. If it is not in the current formation body, it is necessary to load the corresponding formation body velocity file, and use the velocity generation rule of the new formation body to interpolate the velocity value for the starting point of the refracted ray step; when the ray When the step intersects the formation, the reflection ray step is still carried out in the current formation body, but with a new ray direction; the refraction ray step will start from the new formation body, and the new formation body is determined using the modified rules for judging the spatial point horizon , and load the velocity file of the corresponding layer, and use the new tetrahedron replica to interpolate the velocity value for the starting point of the refracted ray step. 6. A method of interface-aware ray tracing based on implicit model expression according to claim 5, characterized in that, the rules for judging the spatial point horizon and the corrected judgment spatial point layer in the steps C3 and C4 Bit rules: The rule for judging the level of a spatial point refers to: start selective traversal from the root node of the stratum binary tree, interpolate the value of a given spatial point from the symbol distance field of the level corresponding to the current binary tree node, and determine the next traversal by the sign of the value Whether the child node is the left child node or the right child node, so continue to traverse until it reaches the leaf node, which gives the level of the spatial point; The revised rule for judging the horizon of a spatial point means: when traversing to a non-leaf node of an intersecting stratum, if the current stratum is located in the left subtree, continue traversing from the right subtree, and vice versa.

Images