CN102053258A - Self-adaptive three-dimensional ray tracing method based on complex geological structure - Google Patents

Abstract

A self-adaptive three-dimensional ray tracing method based on a complex geological structure model comprises the following steps: shooting a group of rays on the target layer in the full-angle range at the excitation point; binning the target layer, and counting the number of the incident rays on each bin of the target layer; encrypting the trial injection for the bin with the ray density lower than a preset threshold value; calculating the intersection point of the ray and the target layer according to the initial angle of the ray, and calculating the coordinates of the emergent point of the ray in the receiving area; searching for an emergent point of 3 initial rays closest to the receiving point in the receiving area, and correcting the initial angle of the rays at the exciting point; and calculating an intersection point of the ray and the target layer according to the corrected initial angle of the ray, and sequentially connecting the seismic source excitation point, the intersection point calculated according to the corrected initial angle and the receiving point to obtain the three-dimensional ray tracking path. The invention can rapidly eliminate a large number of triangles without intersection points, thereby reducing the calculation amount; through the area search, the ray search space is reduced, and the ray search speed is accelerated.

Description

Self-adaptive three-dimensional ray tracing method based on complex geological structure

Technical Field

The invention relates to seismic exploration data acquisition, in particular to a three-dimensional ray tracing method based on a complex geological structure model, which is used for the fields of petroleum and natural gas seismic exploration, geological and geological prospecting and the like.

Background

Seismic exploration is a method of receiving refracted and/or reflected waves of seismic waves excited at the earth's surface by shot points (i.e., excitation points) from seismic wave detection points, and processing and analyzing the detected ray paths, propagation times, etc. of the refracted and/or reflected waves generated by the seismic waves at the geologic formation layer interfaces back to the earth's surface in order to accurately understand and determine the morphology and location of the geologic formation. Seismic prospecting methods have been widely used in the fields of petroleum and gas seismic prospecting, geological and geological prospecting, and the like.

Seismology mainly comprises two major core theories, namely a wave theory for revealing the dynamic characteristics of seismic waves, and a ray theory for revealing the dynamic characteristics of the seismic waves. Therefore, two branches of forward modeling method based on fluctuation theory and forward modeling method based on ray theory are generated according to the two theories in the aspect of forward modeling.

The wave theory forward modeling method is to perform gridding on the established model, and iteratively calculate a wave field according to time step by using a difference method or a finite element method and the like according to wave equation approximation to obtain a forward wave field simulation result. The forward modeling method of the wave theory is characterized by being capable of accurately revealing the amplitude, frequency and phase changes of waves propagating in a medium, truly reflecting the dynamic characteristics of the waves and being an effective means for researching wave fields. Theoretically, the forward modeling method of the fluctuation theory can adapt to a very complex structure, but is limited by the size of a model grid in practical application, so that the very complex structure is limited to a certain extent.

The ray theory forward modeling method is that according to the model structure and the speed change, the movement track of seismic waves in an underground medium and the reflection and transmission (including refraction) changes of the waves through a speed interface are tracked, and wave field characteristics are approximately simulated by calculating the energy loss and phase change of the waves such as wave field diffusion, transmission coefficient, reflection coefficient and the like. The ray theory forward method is mainly used for researching wave kinematic characteristics such as wave propagation tracks, reflection point distribution, wave field distribution range and the like, and is also the basis for carrying out chromatographic inversion.

Snell's theorem shows that in the process of propagating seismic waves, when the seismic waves encounter a medium interface, a part of energy is reflected and returns to the original medium, called as reflected waves, and another part of energy penetrates through the interface and enters a new medium, called as transmitted waves, and the incident waves and the transmitted waves follow the following relationship:

wherein alpha is₁、α₂Respectively representing the angle of incidence and the angle of transmission, v₁、v₂Representing the incident medium velocity and the transmissive medium velocity, respectively.

The Fermat principle, the ray travel time minimum principle, means that the path of the seismic wave traveling in the medium is the minimum of travel time. Thus, the ray propagates as a straight line in a homogeneous medium, as a curved line in a heterogeneous medium, and perpendicular to the wavefront surface.

Huygens' principle indicates that a seismic source generates vibration, the vibration propagates around through mass points of a medium to form a wave front, each point on the wave front can be regarded as a new seismic source, the vibration continues to propagate through the mass points of the medium, and a new envelope surface is formed by superposition, namely a new wave front.

According to the Snell theorem, the Fermat principle and the Huygens principle, a plurality of ray tracing methods are formed, and the calculation methods of ray paths can be roughly divided into the following three types:

(1) shooting, or target shooting, is the most common ray tracing method that was first proposed and used. The ray tracing process is as follows: a series of ray parameter initial values are given at an excitation point, then tracking is carried out in sequence according to the Snell theorem, two rays which are closest to a receiving point are selected, the initial ray parameter values are adjusted through interpolation, and a satisfactory result can be obtained through multiple times of adjustment and modification. The method has the greatest advantages that the accurate tracking of rays is realized, and the tracking in a blind area can be avoided; iterative convergence is fast in a relatively simple model structure, but convergence is slow and time-consuming in a complex structure, which is a disadvantage of this method.

(2) And the bending method comprises the steps of gridding the model, starting from the position of a seismic source, calculating the travel time point by point according to the propagation direction of a wave field and the sequence of grid nodes, solving the travel time calculation of the grid by means of a difference equation method or a wave front method according to the Huygens principle and the like, comparing the nodes one by one from a receiving point according to the minimum travel time principle according to the Fermat principle, and finding back the seismic source point to obtain the minimum travel time ray. The bending method fully embodies the fluctuation characteristics of seismic waves and adapts to media with variable speed, but the method has the following defects: firstly, the model has certain limitation on the description of a complex structure through gridding approximation, namely, the requirement of a very complex geological structure cannot be met; secondly, the connecting lines of grid nodes are used for approximating ray paths (even if some methods adopt interpolation processing), and the approximation degree is determined by the size of the grid; searching ray paths according to the Fermat principle, possibly losing the shortest ray path, and judging difficultly when a plurality of ray paths exist, wherein ray blind areas cannot be eliminated; fourthly, a large amount of computer operation time and memory are consumed, and interactive calculation is not facilitated.

(3) And the iteration method is used for establishing a travel time equation according to the model structure, the known interface function and the medium speed, solving the partial derivatives of all unknown parameters and being equal to 0 according to the Fermat principle, approximately expanding the partial derivatives to form an equation set, and continuously iterating and converging the equation set to an accurate solution. The iteration method converts the ray tracing process into an iteration problem of a given initial value, so that the ray tracing process is simplified in form, and the convergence rate is high; the iteration speed is related to an initial ray path, the closer the initial ray path is to the iteration speed, the faster the iteration speed is, and otherwise, the slower the speed is; secondly, the method is only suitable for a layered structure model and a medium with uniform speed; solving multiple ray paths has certain difficulty.

As can be seen from ray theory, the propagation of waves in a medium satisfies the Snell theorem, the Fermat principle and the Huygens principle, and all current ray tracing methods basically expand around the three basic principles.

Disclosure of Invention

According to an aspect of the present invention, there is provided a three-dimensional ray tracing method based on a complex geological structure model, comprising the following steps: shooting a group of rays on a target layer of the geological structure model in a full-angle range at a seismic source excitation point; binning the target layer, and counting the number of the incident rays on each bin of the target layer; for a surface element with ray density lower than a preset threshold value, carrying out encryption trial injection on the surface element; calculating the intersection point of the pilot ray and the target layer according to the initial angle of the pilot ray, and then calculating the coordinates of the emergent point of the pilot ray in the receiving area; searching for the emergent points of 3 initial rays closest to the receiving point in the receiving area, and correcting the initial angle of the rays at the excitation point through triangular interpolation to ensure that the emergent points of the pilot rays in the receiving area gradually converge to the receiving point; and calculating an intersection point of the ray and the target layer according to the corrected initial angle of the ray, and sequentially connecting the seismic source excitation point, the intersection point calculated according to the corrected initial angle and the receiving point, thereby obtaining the three-dimensional ray tracking path.

Wherein, in the step of correcting the initial angles of the rays at the excitation point, the initial angles of the 3 initial rays are respectively (theta)₁，φ₁)、(θ₂，φ₂) And (theta)₃，φ₃) Wherein, theta₁、θ₂、θ₃Respectively is the included angle, phi, between the 3 initial rays and the Z coordinate axis in the three-dimensional space₁、φ₂、φ₃Respectively, the included angles between the projection on the horizontal plane of the 3 initial rays XOY and the X coordinate axis, and the corrected initial angle (theta) of the rays_p，φ_p) Comprises the following steps:

<math><mrow><msub><mi>&theta;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&theta;</mi><mn>1</mn></msub><msub><mi>S</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>2</mn></msub><msub><mi>S</mi><mn>3</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>3</mn></msub><msub><mi>S</mi><mn>1</mn></msub></mrow><mi>S</mi></mfrac><mo>=</mo><mfrac><msub><mi>S</mi><mn>2</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>1</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>3</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>2</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>1</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>3</mn></msub><mo>,</mo></mrow></math>

<math><mrow><msub><mi>&phi;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&phi;</mi><mn>1</mn></msub><msub><mi>S</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&phi;</mi><mn>2</mn></msub><msub><mi>S</mi><mn>3</mn></msub><mo>+</mo><msub><mi>&phi;</mi><mn>3</mn></msub><msub><mi>S</mi><mn>1</mn></msub></mrow><mi>S</mi></mfrac><mo>=</mo><mfrac><msub><mi>S</mi><mn>2</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>1</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>3</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>2</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>1</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>3</mn></msub><mo>,</mo></mrow></math>

wherein S is₁、S₂、S₃The areas of the 3 small triangles formed by the emitting points and the receiving points of the 3 initial rays are respectively, and S is the area of the triangle formed by the emitting points and the receiving points of the 3 initial rays.

The step of calculating the intersection of the shooting ray and the target slice may include: (1) dividing all triangles forming a target layer into two subregions according to the number of the triangles, respectively calculating X, Y, Z coordinate extreme values of the two subregions under a three-dimensional coordinate system XYZ, constructing a space cuboid according to the X, Y, Z coordinate extreme value obtained by calculation, respectively projecting a pilot ray and the cuboid to XOY, XOZ and YOZ planes, and determining that the pilot ray and the corresponding subregions have intersection points when the projection of the ray and the projection of the space cuboid on the three planes XOY, XOZ and YOZ are intersected with each other; (2) dividing the sub-region which is determined to be intersected with the shooting ray into two sub-regions again, and respectively judging whether the shooting ray is intersected with the sub-regions which are divided again or not in the same mode as the step (1); (3) and (3) repeatedly executing the step (2) until the number of the triangles in the divided sub-regions is less than a specific threshold, not dividing the sub-regions, and directly performing intersection calculation on the pilot ray and each triangle in the finally divided sub-regions to calculate the intersection.

The specific threshold value can be an integer between 5 and 10.

The specific threshold may be 8.

The step of calculating the intersection point may comprise: calculating the intersection point of the shooting ray and the plane of the triangle according to the initial angle of the shooting ray; converting the coordinates of the intersection points into triangle area coordinates according to the coordinates of three vertexes of the triangle; determining whether the intersection point is in the triangle according to the property of the area coordinate of the triangle; if the intersection point is in the triangle, the real intersection point of the shooting ray and the target layer is obtained; if the intersection point is not in a triangle, the intersection point of the ray with the next triangle continues to be calculated until the calculated intersection point is in a triangle.

The triangle area coordinate is defined as (u)₁，u₂，u₃) Wherein

s is the area of the triangle, S₁、S₂、S₃Respectively the three vertices of the triangle and the intersection point

u₁≥0，u₂≥0，u₃≥0；

The area of the formed 3 small triangles is if | u is satisfied₁|+|u₂|+|u₃If 1, the intersection point is determined to be in the triangle, and thus a true intersection point of the incident ray and the target slice is obtained.

Drawings

These and/or other aspects and advantages of the present invention will become apparent and more readily appreciated from the following description of the embodiments, taken in conjunction with the accompanying drawings of which:

FIG. 1 is a schematic diagram of determining whether a ray intersects a spatial cuboid;

FIG. 2 is a flow chart of target level sub-regionalization according to the present invention;

FIG. 3 is a flow chart of target level sub-regionalization using a non-recursive algorithm in accordance with the present invention;

FIG. 4 is a schematic diagram of the spatial relationship of ray intersection with a triangle;

FIG. 5 is a schematic illustration of a triangular area coordinate;

FIG. 6 is a flow chart illustrating a three-dimensional ray tracing method according to the present invention;

FIG. 7 is a schematic illustration of a trial ray initial angle correction calculation;

fig. 8 is a schematic diagram of a ray trace generated using a three-dimensional ray tracing method according to the present invention.

Detailed Description

Hereinafter, embodiments of the present invention are described in detail with reference to the accompanying drawings.

Before describing the three-dimensional ray tracing method based on the complex geological structure model, the hierarchical tree structure organization method of the sub-surface triangular gridding data and the method for rapidly calculating the intersection point of the ray and the target layer (the ground layer or the fault layer) of the geological structure model are described in detail.

Hierarchical tree structure organization method of sub-surface triangular gridding data

The calculation problem of the intersection point of the ray and the interface is very simple in mathematics, but is very critical to three-dimensional ray tracing, especially for the situations of complex structures and particularly large number of triangles of the interface, the calculation amount is very large, and the calculation efficiency is directly related to the core problem of whether the actual use can be met.

Although the basic algorithm for performing ray tracing on the triangulation network stratum model is simple, it is very inefficient and takes a long time to perform simple ray tracing on a huge (usually tens of thousands) number of complex triangulation network stratum models. Because each calculation of the intersection point needs to be compared with all triangles on the triangulation network, although the ray-triangle intersection algorithm is used, it is still very slow for large-scale ray tracing. In fact, a ray only intersects a few triangles in the triangular net, and the intersection algorithm only has the function of judging whether intersection points exist or not in many times. Therefore, if a large number of triangles without intersection points can be eliminated quickly, the ray tracing speed can be improved, and for this reason, the intersection judgment is carried out by introducing a layer sub-zoning algorithm.

Fig. 1 is a schematic diagram of determining whether a ray intersects a spatial cuboid, and fig. 2 is a flow chart of target slice sub-zoning according to the present invention.

The specific steps of the laminated region zoning algorithm are as follows:

the first step is as follows: and (3) uniformly dividing all triangles forming the target level into two sub-areas according to the number of the triangles (bisection method). X, Y, Z coordinate extreme values (namely, the maximum value of the absolute value of the maximum coordinate value and the absolute value of the minimum coordinate value) of the two subregions under the three-dimensional coordinate system XYZ are respectively calculated, and two space cuboids are respectively constructed according to the calculated X, Y, Z coordinate extreme values. It is determined whether the ray intersects the bedding plane, that is, whether the ray intersects the spatial cuboid of the structure (i.e., whether the ray intersects the two spatial cuboids of the structure, respectively). Fig. 1 shows only one spatial cuboid, and the other spatial cuboid is omitted. The ray and the space cuboid are respectively projected to XOY, XOZ and YOZ planes, intersection judgment of the ray and the space cuboid is simplified into judgment of whether a straight line and the plane are intersected or not, intersection points of the ray and the subregion can be determined only if the ray and the space cuboid are intersected on the three planes XOY, XOZ and YOZ, and otherwise, the ray and the subregion are not intersected.

The second step is as follows: and determining a sub-region intersected with the ray, averagely dividing the determined sub-region into two sub-regions again, and respectively judging whether the ray is intersected with the sub-region which is divided again by adopting the same method (calculating a coordinate extreme value, constructing a space cuboid and the like) as the first step.

The third step: and repeating the second step until the sub-region is not divided when the number of the triangles in the divided sub-region is less than a specific threshold (for example, the specific threshold may be an integer between 5 and 10, and preferably may be 8), and directly performing intersection calculation on the ray and each triangle in the sub-region to find the intersection. At the moment, only a small amount of intersection calculation is needed, and the calculation amount can be obviously reduced.

The essence of the bedding sub-zoning algorithm is to divide the bedding into smaller bedding, gradually lock the triangles with possible intersection points in a small range, and then perform intersection point calculation of the rays and the triangles. A large number of triangles without intersection points can be rapidly excluded by using the level sub-zoning algorithm, so that the calculation amount can be remarkably reduced.

Due to the low efficiency of the recursive algorithm, the recursive algorithm is implemented non-recursively in the programming implementation, and the specific algorithm flow is shown in fig. 3. The calculation is mainly carried out by using a stack data type in a circulating mode.

Method for rapidly calculating intersection point of ray and target layer (ground layer or fault layer) of geological structure model

The three-dimensional block structure model is described by points, triangles, faces and blocks, the blocks are formed by face closure, the faces are composed of triangular faces, each triangular face is determined by three control points, the medium in each closed block is uniform, namely, the density in the closed block is the same, and the rays have the same speed in the closed block. The ray travels straight within the closed block and is reflected or transmitted when encountering the boundary of the block, requiring the determination of the point at which the ray is reflected and transmitted on the boundary surface, i.e., the intersection of the ray with the triangular face on the face. Therefore, an important link in the ray tracing process of the three-dimensional block structure model is to calculate the intersection point of the ray and the triangular surface, and the speed and the efficiency of the ray tracing of the shooting are determined to a great extent. Calculating the intersection of a ray with a triangular face actually involves two problems: firstly, an intersection point of an intersected triangular surface is quickly calculated; secondly, whether the intersection point is the real intersection point of the ray and the triangular surface is judged rapidly (namely whether the intersection point is positioned inside the triangle (including the boundary)).

First, the problem of calculating the intersection of a ray with a triangle surface is discussed, and FIG. 4 is a schematic diagram of the spatial relationship of the ray intersection with the triangle. The ray is a spatial straight line, and an initial coordinate point P is known₀(x₀，y₀，z₀) The direction numbers of the rays are l, m and n respectively. The ray equation can be expressed as:

$\left\{\begin{array}{c}x=x_0+\mathrm{lt}\\ y={\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="HSA00000386477300051.png"\; state="\{\{state\}\}">y</figure-callout>}}_0+\mathrm{mt}\\ z={\mathrm{<figure-callout\; id="0"\; label="z"\; filenames="HSA00000386477300051.png"\; state="\{\{state\}\}">z</figure-callout>}}_0+\mathrm{nt}\end{array}\right.---\left(1\right)$

wherein t is a variable. For triangular faces, the three fixed-point coordinates P of the triangle are known₁(x₁，y₁，z₁)、P₂(x₂，y₂，z₂)、P₃(x₃，y₃，z₃) Then the triangular plane equation composed of it is expressed as:

$\left|\begin{array}{ccc}{x-x}_1& y-y_1& z-z_1\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=0---\left(2\right)$

now assume that the ray intersects the triangular face at a point P (x)_p，y_p，z_p) Since the intersection point is a point on the ray equation and a point on the triangular surface, the equations (1) and (2) are replaced by:

$\left\{\begin{array}{c}{x}_p=x_0+\mathrm{lt}\\ y_p={\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="HSA00000386477300051.png"\; state="\{\{state\}\}">y</figure-callout>}}_0+\mathrm{mt}\\ z_p={\mathrm{<figure-callout\; id="0"\; label="z"\; filenames="HSA00000386477300051.png"\; state="\{\{state\}\}">z</figure-callout>}}_0+\mathrm{nt}\end{array}\right.---\left(3\right)$

$\left|\begin{array}{ccc}{x_p-x}_1& y_p-y_1& z_p-z_1\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=0---\left(4\right)$

substituting equation (3) for equation (4) yields a functional equation for variable t:

$\left|\begin{array}{ccc}{x_0+\mathrm{lt}-x}_1& {\mathrm{<figure-callout\; id="0"\; label="y"\; filenames="HSA00000386477300051.png"\; state="\{\{state\}\}">y</figure-callout>}}_0+\mathrm{mt}-y_1& {\mathrm{<figure-callout\; id="0"\; label="z"\; filenames="HSA00000386477300051.png"\; state="\{\{state\}\}">z</figure-callout>}}_0+\mathrm{nt}-z_1\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=0---\left(5\right)$

to simplify the equation above, let:

$\left\{\begin{array}{c}{x}_{\mathrm{ij}}=x_i-x_j\\ y_{\mathrm{ij}}=y_i-y_j\\ z_{\mathrm{ij}}=z_i-z_j\end{array}\right.---\left(6\right)$

replacing with formula (5), formula (6) becomes:

$\left|\begin{array}{ccc}{\mathrm{lt}-x}_{10}& \mathrm{mt}-y_{10}& \mathrm{nt}-z_{10}\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=0---\left(7\right)$

by developing the above equation, the computational expression of the obtainable variable t is:

$t=\frac{a_1{x}_{10}+a_2{y}_{10}+a_3{y}_{10}}{{\mathrm{la}}_1+{\mathrm{<figure-callout\; id="2"\; label="ma"\; filenames="HSA00000386477300011.png,HSA00000386477300051.png"\; state="\{\{state\}\}">ma</figure-callout>}}_2+{\mathrm{na}}_3}---\left(8\right)$

wherein:

$\left\{\begin{array}{c}{a}_1=y_{21}{z}_{31}-y_{31}{z}_{21}\\ a_2=z_{21}{x}_{31}-z_{31}{x}_{21}\\ a_3=x_{21}{y}_{31}-x_{31}{y}_{21}\end{array}\right.---\left(9\right)$

t is obtained from the equation (7), and the intersection point P (x) of the ray and the triangular surface can be calculated by substituting t into the equation (6)_p，y_p，z_p) The coordinate values of (2).

In the above, the problem of calculating the intersection point of the ray and the triangle surface is discussed, and actually, the calculated intersection point only represents the intersection point of the ray and the plane where the triangle is located, and whether the intersection point is in the triangle needs to be further determined, and the problem of determining the intersection point will be discussed below.

By using threeThe intersection point is determined by the angular area coordinates, and as shown in FIG. 5, three coordinate points of the triangle are assumed to be P₁(x₁，y₁，z₁)、P₂(x₂，y₂，z₂)、P₃(x₃，y₃，z₃) The intersection point is P (x)_p，y_p，z_p). Thus forming a large triangle Δ P₁P₂P₃And three small triangles: delta P₁P₂P、ΔP₂P₃P and Δ P₃P₁And P. Large triangle delta P₁P₂P₃The area is marked as S, and the three small triangles are respectively marked as S₁、S₂And S₃. The area values of the corresponding triangles can be respectively calculated according to the coordinates, and are expressed as:

$\left\{\begin{array}{c}S=\frac{1}{2}[x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2)]\\ S_1=\frac{1}{2}[x_1(y_2-y_p)+x_2(y_p-y_1)+x_p(y_1-y_2)]\\ {\mathrm{<figure-callout\; id="2"\; label="S"\; filenames="HSA00000386477300011.png,HSA00000386477300051.png"\; state="\{\{state\}\}">S</figure-callout>}}_2=\frac{1}{2}[x_2(y_3-y_p)+x_2(y_p-y_2)+x_p(y_2-y_3)]\\ S_3=\frac{1}{2}[x_3(y_1-y_p)+x_1(y_p-y_3)+x_p(y_3-y_1)]\end{array}\right.---\left(10\right)$

the area of the triangle is calculated by the above formula, the result of the coordinate point in the counterclockwise direction is a positive value, and if the result of the coordinate point in the clockwise direction is a negative value, the position of the intersection point can be determined by using the characteristic. For triangular area coordinates, at triangle Δ P₁P₂P₃Can be represented by an area coordinate, which for point P is represented by (u)₁，u₂，u₃) Each component is defined as follows:

$u_1=\frac{S_1}{S},{\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000386477300011.png,HSA00000386477300051.png"\; state="\{\{state\}\}">u</figure-callout>}}_2=\frac{S_2}{S},u_3=\frac{S_3}{S}---\left(11\right)$

as can be seen from the definition of equation (11), the three components represent the ratio of the area of three small triangles divided by the point P to the area of the original triangle. The properties of the triangular area coordinates can be obtained very easily:

(1)u₁+u₂+u₃≡1

(2) when a point is inside a triangle (including the boundary vertices), there are:

u₁≥0，u₂≥0，u₃≥0；

|u₁|+|u₂|+|u₃|＝1

(3) when the P point is outside the triangle:

|u₁|+|u₂|+|u₃|＞1

therefore, the coordinates of the intersection points of the rays and the triangle are converted into the area coordinates of the triangle according to the coordinates of the triangle, and whether the intersection points are in the triangle or not can be quickly judged by utilizing the property of the area coordinates. If the intersection point is in the triangle, the real intersection point of the ray and the surface is obtained, otherwise, the intersection point of the ray and the next triangle is continuously searched until the real intersection point is found.

The three-dimensional ray tracing method based on the complex geological structure model of the invention is described in detail below. Fig. 6 shows a flow chart of a three-dimensional ray tracing method according to the invention.

Referring to FIG. 6, at step 601, a set of rays is shot at a full angular range at a source shot point onto a target horizon of a geological formation model. The bedding structure of the geological structure model is known.

When the next iteration is carried out, the requirement that the geological structures passed by the ray tracks participating in the iteration are basically the same or the regions of the ray reflection points participating in the iteration calculation are relatively flat is met, otherwise the iteration is not converged. Therefore, if the shooting rays are too dense, the calculation efficiency is reduced; if the shooting ray is too thin, the iteration may not converge. Therefore, a relatively sparse pilot can be performed first, and then an encrypted pilot can be performed on the basis of the sparse pilot.

In step 602, the target plane is binned, and the number of incident rays on each bin of the target plane is counted. In step 603, for a bin with a ray density lower than a predetermined threshold, the bin is subjected to encrypted trial injection. In order to facilitate ray density organization and rapid ray retrieval, the target layer is subjected to binning of a certain size, and a reflection point of each incident ray belongs to a certain bin. The organization method adopts area search, reduces the ray search space and accelerates the ray search speed.

In step 604, the intersection of the pilot ray and the target layer is calculated based on the initial angle of the pilot ray, and then the coordinates of the exit point of the pilot ray in the receiving area are calculated.

The calculation method of the intersection of a single ray and the target level has been explained in detail in the foregoing description (see "hierarchical tree structure organization method of triangular gridded data of a sub-surface" and "fast calculation method of the intersection of a ray and the target level (level or fault level) of the geological structure model"), and the description is not repeated here. By the calculation method, the calculation amount can be obviously reduced, and the real intersection point of the single ray and the target layer can be quickly obtained.

From the calculated position (coordinates) of the true intersection point, the coordinates of the exit point of the pilot ray in the receiving area can be obtained according to the reflection/transmission (refraction) theorem.

It should be appreciated that the geologic formation model may include one or more target layers, and thus, the incident radiation may be reflected/transmitted multiple times. Since the initial angle of the pilot ray and the layer structure of the geological structure model are known, the reflection points (intersection points) of the pilot ray and each target layer can be calculated in sequence.

Specifically, for a target slice, the intersection of the incident ray with the target slice can be directly calculated according to the initial angle of the incident ray.

For a plurality of target layers, calculating the intersection point of the ray and the uppermost target layer according to the initial angle of the pilot ray, reflecting and transmitting the ray at the intersection point, calculating the intersection point of the transmitted ray and the target layer which is immediately below the uppermost target layer according to the transmission angle, and so on until the intersection point of the transmitted ray and the lowermost target layer is calculated in the same way, and the transmitted ray is reflected at the intersection point; then, the intersection of the reflected ray with the target level immediately above the lowest target level is calculated from the reflection angle, the reflected ray is transmitted and reflected at this intersection, the intersection of the transmitted ray with the next target level is calculated from the transmission angle, and so on until the intersection of the transmitted ray with the highest target level is calculated in the same manner.

In step 605, the exit points of the 3 initial rays closest to the receiving point are searched in the receiving area, and the initial angle of the ray at the excitation point is corrected by triangle interpolation so that the exit point of the pilot ray in the receiving area gradually converges to the receiving point, and the ray path of the receiving point is obtained.

Specifically, a three-dimensional ray tracing iterative algorithm is adopted for ray path tracing.

Ray tracing is forward to obtain ray paths from a seismic source excitation point to a receiving point under the condition that the seismic source excitation point and the receiving point are known so as to obtain a wave field record, and various related applications and analysis can be carried out on the basis of the ray paths. Therefore, the most basic work of ray tracing is two-point ray tracing, and a three-dimensional shooting ray method has been discussed above, so that the tracing calculation of a single ray after a given initial ray angle is realized, and the three-dimensional ray two-point tracing is performed on the basis of a plurality of single ray tracing.

For the set of rays for the target slice in step 601, the initial angle of the set of rays at the excitation point is (θ)_i，φ_i)，θ_iIs the angle between the ray i and the Z coordinate axis in three-dimensional space, phi_iThe projection of the ray i on the XOY horizontal plane is at an angle to the X coordinate axis, i denoting the number from which the initial ray was emitted. Through angle theta_iAnd phi_iThe direction angle of the ray i in three-dimensional space can be determined. The method comprises the steps of forming an emergent point distribution in a receiving area on the ground through trial ray tracing, searching an emergent point of an initial ray closest to a receiving point according to the coordinates of the receiving point, correcting the emitting angle of the ray at an excitation point through interpolation, gradually converging the ray to the receiving point, and obtaining a ray path of the receiving point.

How to correct the ray initial angle is the most critical problem of three-dimensional incident ray convergence, and the problem directly influences the speed of the incident ray convergence. In the case of three-dimensional, the trial angle is corrected by triangle interpolation, and as shown in FIG. 7, the reception point is assumed to be P (x)_p，y_p，z_p) The three closest initial ray emitting points around the receiving point are respectively P₁(x₁，y₁，z₁)、P₂(x₂，y₂，z₂) And P₃(x₃，y₃，z₃) Corresponding to an initial angle of (θ)₁，φ₁)、(θ₂，φ₂) And (theta)₃，φ₃) The initial angle of P point correction to be solved is (theta)_p，φ_p). Using a one-dimensional linear interpolation equation, the azimuth angle theta_pAnd phi_pRespectively comprises the following steps:

<math><mrow><mfenced open='|' close='|'><mtable><mtr><mtd><msub><mi>&theta;</mi><mn>1</mn></msub><mo>-</mo><msub><mi>&theta;</mi><mi>p</mi></msub></mtd><mtd><msub><mi>&theta;</mi><mn>2</mn></msub><mo>-</mo><msub><mi>&theta;</mi><mi>p</mi></msub></mtd><mtd><msub><mi>&theta;</mi><mn>3</mn></msub><mo>-</mo><msub><mi>&theta;</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>x</mi><mn>1</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>2</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>3</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mn>1</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>2</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>3</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mn>0</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow></math>

<math><mrow><mfenced open='|' close='|'><mtable><mtr><mtd><msub><mi>&phi;</mi><mn>1</mn></msub><mo>-</mo><msub><mi>&phi;</mi><mi>p</mi></msub></mtd><mtd><msub><mi>&phi;</mi><mn>2</mn></msub><mo>-</mo><msub><mi>&phi;</mi><mi>p</mi></msub></mtd><mtd><msub><mi>&phi;</mi><mn>3</mn></msub><mo>-</mo><msub><mi>&phi;</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>x</mi><mn>1</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>2</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>3</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mn>1</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>2</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>3</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mn>0</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow></math>

the formula (12) and the formula (13) respectively represent the correction azimuth angle theta to be obtained_pAnd phi_pThe corresponding solution can be obtained by developing the two equations. Considering equation (12), we develop:

<math><mrow><mrow><mo>(</mo><msub><mi>&theta;</mi><mn>1</mn></msub><mo>-</mo><msub><mi>&theta;</mi><mi>p</mi></msub><mo>)</mo></mrow><mfenced open='|' close='|'><mtable><mtr><mtd><msub><mi>x</mi><mn>2</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>3</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mn>2</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>3</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd></mtr></mtable></mfenced><mo>+</mo><mrow><mo>(</mo><msub><mi>&theta;</mi><mn>2</mn></msub><mo>-</mo><msub><mi>&theta;</mi><mi>p</mi></msub><mo>)</mo></mrow><mfenced open='|' close='|'><mtable><mtr><mtd><msub><mi>x</mi><mn>3</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>1</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mn>3</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>1</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd></mtr></mtable></mfenced><mo>+</mo><mrow><mo>(</mo><msub><mi>&theta;</mi><mn>3</mn></msub><mo>-</mo><msub><mi>&theta;</mi><mi>p</mi></msub><mo>)</mo></mrow><mfenced open='|' close='|'><mtable><mtr><mtd><msub><mi>x</mi><mn>1</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd><mtd><msub><mi>x</mi><mn>2</mn></msub><mo>-</mo><msub><mi>x</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mn>1</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd><mtd><msub><mi>y</mi><mn>2</mn></msub><mo>-</mo><msub><mi>y</mi><mi>p</mi></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mn>0</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>14</mn><mo>)</mo></mrow></mrow></math>

to simplify equation (14), let:

$a_1=\left|\begin{array}{cc}{x}_2-x_p& x_3-x_p\\ y_2-y_p& y_3-y_p\end{array}\right|=x_2(y_3-y_p)+x_3(y_p-y_2)+x_p(y_2-y_3)$

$a_2=\left|\begin{array}{cc}{x}_3-x_p& x_1-x_p\\ y_3-y_p& y_1-y_p\end{array}\right|=x_3(y_1-y_p)+x_1(y_p-y_3)+x_p(y_3-y_1)$

$a_3=\left|\begin{array}{cc}{x}_1-x_p& x_2-x_p\\ y_1-y_p& y_2-y_p\end{array}\right|=x_1(y_2-y_p)+x_2(y_p-y_1)+x_p(y_1-y_2)---\left(15\right)$

the azimuth angle θ is obtained from the equations (14) and (15)_pExpression of (2):

<math><mrow><msub><mi>&theta;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&theta;</mi><mn>1</mn></msub><msub><mi>a</mi><mn>1</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>2</mn></msub><msub><mi>a</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>3</mn></msub><msub><mi>a</mi><mn>3</mn></msub></mrow><mrow><msub><mi>a</mi><mn>1</mn></msub><mo>+</mo><msub><mi>a</mi><mn>2</mn></msub><mo>+</mo><msub><mi>a</mi><mn>3</mn></msub></mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mrow></math>

referring to fig. 7, comparing equation (15) with equation (16), there are:

a₁＝2S₂，a₂＝2S₃，a₃＝2S₁ (17)

and:

S＝S₁+S₂+S₃ (18)

wherein S is₁、S₂、S₃The areas of the 3 triangles formed by the three initial ray emitting points and the receiving points are respectively, and S is the area of the triangle formed by the three initial ray emitting points. The area S can be calculated according to the coordinates of the three initial ray emitting points and the coordinates of the receiving points₁、S₂、S₃And S.

Equation (16) can be expressed as the area of the triangle:

<math><mrow><msub><mi>&theta;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&theta;</mi><mn>1</mn></msub><msub><mi>S</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>2</mn></msub><msub><mi>S</mi><mn>3</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>3</mn></msub><msub><mi>S</mi><mn>1</mn></msub></mrow><mi>S</mi></mfrac><mo>=</mo><mfrac><msub><mi>S</mi><mn>2</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>1</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>3</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>2</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>1</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>3</mn></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>19</mn><mo>)</mo></mrow></mrow></math>

equation (19) indicates that the azimuth angle (θ) of the point is known₁，θ₂，θ₃) And the azimuth angle theta to be found_pThe linear relation of (1) and the ratio of the area of three small triangles with coefficients determined by the position of the point P to be solved to the area of a triangle formed by known three points, the physical meaning is also very clear, for example, the closer the point P is to the point P₁Dot, S₂The larger the area of (A), P₁The specific gravity of the point is large, theta_pIs closer to theta₁. Similarly, the azimuth angle φ can be obtained similarly to the formula (19)_pThe calculation formula (c) is as follows:

<math><mrow><msub><mi>&phi;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&phi;</mi><mn>1</mn></msub><msub><mi>S</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&phi;</mi><mn>2</mn></msub><msub><mi>S</mi><mn>3</mn></msub><mo>+</mo><msub><mi>&phi;</mi><mn>3</mn></msub><msub><mi>S</mi><mn>1</mn></msub></mrow><mi>S</mi></mfrac><mo>=</mo><mfrac><msub><mi>S</mi><mn>2</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>1</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>3</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>2</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>1</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>3</mn></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>20</mn><mo>)</mo></mrow></mrow></math>

the initial ray correction angle (theta) at the receiving point is derived from the three exit points closest to the receiving point position_p，φ_p) The method solves the problem of two-point ray tracing convergence of the shooting ray.

In step 606, an intersection point of the ray and the target layer is calculated according to the corrected initial angle of the ray, and the seismic source excitation point, the intersection point calculated according to the corrected initial angle, and the receiving point are sequentially connected, thereby obtaining a three-dimensional ray tracing path.

Fig. 8 is a schematic diagram of a ray trace generated using a three-dimensional ray tracing method according to the present invention.

According to the three-dimensional ray tracing method based on the complex geological structure model, aiming at the relief surface and the three-dimensional block model containing the inverse fault complex geological structure, the ray tracing is completed by adopting the processes of three-dimensional shooting, encrypted shooting and iteration. According to the invention, a large number of triangles without intersection points can be rapidly eliminated by using a layer sub-zoning algorithm, so that the calculated amount is reduced; in the processes of shooting and encryption, through area search, the ray search space is reduced, and the ray search speed is accelerated.

While the present invention has been particularly shown and described with reference to exemplary embodiments thereof, it will be understood by those of ordinary skill in the art that various changes in form and details may be made therein without departing from the spirit and scope of the present invention as defined by the following claims.

Claims (7)

1. A three-dimensional ray tracing method based on a complex geological structure model comprises the following steps:

shooting a group of rays on a target layer of the geological structure model in a full-angle range at a seismic source excitation point;

binning the target layer, and counting the number of the incident rays on each bin of the target layer;

for a surface element with ray density lower than a preset threshold value, carrying out encryption trial injection on the surface element;

calculating the intersection point of the pilot ray and the target layer according to the initial angle of the pilot ray, and then calculating the coordinates of the emergent point of the pilot ray in the receiving area;

searching for the emergent points of 3 initial rays closest to the receiving point in the receiving area, and correcting the initial angle of the rays at the excitation point through triangular interpolation to ensure that the emergent points of the pilot rays in the receiving area gradually converge to the receiving point;

and calculating an intersection point of the ray and the target layer according to the corrected initial angle of the ray, and sequentially connecting the seismic source excitation point, the intersection point calculated according to the corrected initial angle and the receiving point, thereby obtaining the three-dimensional ray tracking path.

2. The three-dimensional ray tracing method of claim 1, wherein, in the step of correcting the initial angles of the rays at the excitation point, the initial angles of the 3 initial rays are respectively (θ [)₁，φ₁)、(θ₂，φ₂) And (theta)₃，φ₃) Wherein, theta₁、θ₂、θ₃Respectively is the included angle, phi, between the 3 initial rays and the Z coordinate axis in the three-dimensional space₁、φ₂、φ₃Respectively, the included angles between the projection on the horizontal plane of the 3 initial rays XOY and the X coordinate axis, and the corrected initial angle (theta) of the rays_p，φ_p) Comprises the following steps:

<math><mrow><msub><mi>&theta;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&theta;</mi><mn>1</mn></msub><msub><mi>S</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>2</mn></msub><msub><mi>S</mi><mn>3</mn></msub><mo>+</mo><msub><mi>&theta;</mi><mn>3</mn></msub><msub><mi>S</mi><mn>1</mn></msub></mrow><mi>S</mi></mfrac><mo>=</mo><mfrac><msub><mi>S</mi><mn>2</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>1</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>3</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>2</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>1</mn></msub><mi>S</mi></mfrac><msub><mi>&theta;</mi><mn>3</mn></msub><mo>,</mo></mrow></math>

<math><mrow><msub><mi>&phi;</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>&phi;</mi><mn>1</mn></msub><msub><mi>S</mi><mn>2</mn></msub><mo>+</mo><msub><mi>&phi;</mi><mn>2</mn></msub><msub><mi>S</mi><mn>3</mn></msub><mo>+</mo><msub><mi>&phi;</mi><mn>3</mn></msub><msub><mi>S</mi><mn>1</mn></msub></mrow><mi>S</mi></mfrac><mo>=</mo><mfrac><msub><mi>S</mi><mn>2</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>1</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>3</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>2</mn></msub><mo>+</mo><mfrac><msub><mi>S</mi><mn>1</mn></msub><mi>S</mi></mfrac><msub><mi>&phi;</mi><mn>3</mn></msub><mo>,</mo></mrow></math>

wherein S is₁、S₂、S₃The areas of the 3 small triangles formed by the emitting points and the receiving points of the 3 initial rays are respectively, and S is the area of the triangle formed by the emitting points and the receiving points of the 3 initial rays.

3. The three-dimensional ray tracing method of claim 1, wherein the step of calculating an intersection of the shooting ray with the target slice comprises:

(1) dividing all triangles forming the target level into two sub-areas according to the number of the triangles on average, respectively calculating X, Y, Z coordinate extreme values of the two sub-areas under an XYZ three-dimensional coordinate system, constructing a space cuboid according to the X, Y, Z coordinate extreme values obtained through calculation, respectively projecting the pilot rays and the cuboid to XOY, XOZ and YOZ planes, and determining that intersection points exist between the pilot rays and the corresponding sub-areas when the projections of the rays and the space cuboid on the XOY, XOZ and YOZ planes are intersected with each other;

(2) averagely dividing the sub-region which is determined to be intersected with the shooting ray into two sub-regions again, and respectively judging whether the shooting ray is intersected with the sub-regions which are divided again or not in the same mode as the step (1);

(3) and (3) repeatedly executing the step (2) until the number of the triangles in the divided sub-regions is less than a specific threshold, not dividing the sub-regions, and directly performing intersection calculation on the pilot ray and each triangle in the finally divided sub-regions to calculate the intersection.

4. The three-dimensional ray tracing method according to claim 3, wherein said specific threshold value is an integer between 5 and 10.

5. The three-dimensional ray tracing method of claim 4, wherein said specific threshold value is 8.

6. The three-dimensional ray tracing method of claim 3, wherein the step of calculating an intersection point comprises:

calculating the intersection point of the shooting ray and the plane of the triangle according to the initial angle of the shooting ray;

converting the coordinates of the intersection points into triangle area coordinates according to the coordinates of three vertexes of the triangle;

determining whether the intersection point is in the triangle according to the property of the area coordinate of the triangle;

if the intersection point is in the triangle, the real intersection point of the shooting ray and the target layer is obtained; if the intersection point is not in a triangle, the intersection point of the ray with the next triangle continues to be calculated until the calculated intersection point is in a triangle.

7. The three-dimensional ray tracing method of claim 6, wherein a triangle area coordinate is defined as (u)₁，u₂，u₃)，

Wherein,

s is the area of the triangle, S₁、S₂、S₃Three of the triangle respectivelyThe area of 3 small triangles formed by the vertex and the intersection point,

u₁≥0，u₂≥0，u₃≥0；

wherein if | u is satisfied₁|+|u₂|+|u₃If 1, the intersection point is determined to be in the triangle, and thus a true intersection point of the incident ray and the target slice is obtained.

Images