CN105068133A - Ray tracking method used for complex geological structure model - Google Patents

Abstract

The invention provides a ray tracking method used for a complex geological structure model. The method comprises the steps of (a) constructing an angle space of initial incidence angles of trial rays into an angle grid topology comprising m*n grid points; (b) at excitation points, trying to emit a group of rays to a geological structure model by using corresponding angles ([Alpha]i, [Beta]j) of each grid point as the initial incidence angles; (c) according to the initial incidence angles of the trial rays, determining emergence points of trial rays in a receiving region after the trial rays are reflected by an object layer interface; (d) determining an emergence triangle of the receiving points; (e) calculating modified initial incidence angles in an iteration manner according to initial incidence angles of three trial rays corresponding to the emergence triangle of the receiving points; and (f) at the excitation points, trying to emit rays to the geological structure model by using the modified initial incidence angles as the initial incidence angles, determining the reflecting points of the trial rays reflected by the object layer interface, and connecting the excitation points, the reflecting points and the receiving points in sequence so as to three-dimensional ray tracking paths.

Description

For the ray-tracing scheme of the geological structure model of complexity

Technical field

All things considered of the present invention relates to technical field of geophysical exploration, more particularly, relates to a kind of ray-tracing scheme of the geological structure model for complexity.

Background technology

In method of seismic exploration, from seismic event geophone station receive the seismic event that excited by shot point in tectonic structure bed interface by reflection or refraction the refraction wave that produces and/or reflection wave, process by the raypath, travel-time etc. that the refraction wave detected and/or reflection wave are turned back to geophone station place and analyze, thus accurately can understand and determine architectonic form and position.Method of seismic exploration be widely used in petroleum gas seismic prospecting, the fields such as quality ore deposit.

Ray-tracing scheme in method of seismic exploration is as one wave field approximate calculation method fast and effectively, not only significant for theory of earthquake wave research, and directly applies to the process such as seismic data inversion and migration imaging.

In the tectonic structure (such as, the tectonic structure containing overthrust fault) of complexity, each big gun inspection between there are many rays.But, existing ray-tracing scheme cannot determine rapidly and accurately each big gun inspection between each bar ray.

In addition, traditional shot point and geophone station are all distributed in earth's surface, then there are following two kinds of situations: geophone station is distributed in earth's surface, and shot point is distributed in well in three-dimensional VSP recording geometry; Shot point is distributed in earth's surface, and geophone station is distributed in well.Because well is an infinitesimal line in three dimensions in logic, cause the ray-tracing scheme based on triangle area process of iteration to be difficult to be distributed in earth's surface at shot point, apply when geophone station is distributed in well.

Summary of the invention

Exemplary embodiment of the present invention is the ray-tracing scheme providing a kind of geological structure model for complexity, to solve existing ray-tracing scheme cannot carry out ray tracing rapidly and accurately problem to the tectonic structure of complexity.

According to exemplary embodiment of the present invention, a kind of ray-tracing scheme of the geological structure model for complexity is provided, comprise: the angular region of the initial incidence angle of shoot ray is configured to comprise the angle network topology of m × n net point by (a), wherein, each net point corresponding angle (α _i, β _j), α _ifor shoot ray in three-dimensional system of coordinate with angle formed by Z axis negative direction, β _jangle formed by the projection of shoot ray in XOY plane and X-axis positive dirction, i is greater than the integer that 0 is less than or equal to m, α _i> α _i-1, j is greater than the integer that 0 is less than or equal to n, β _j> β _j-1; (b) at shot point place, with the angle (α corresponding to each net point of described angle network topology _i, β _j) fire one group of ray for adjustment for geological structure model described in original incident angle; C according to the initial incidence angle of each shoot ray, () determines that each shoot ray is via the eye point in receiving area after destination layer boundary reflection; D () determines the outgoing triangle at acceptance point place, wherein, two the non-overlapping triangles of described outgoing triangle for being formed for summit with the eye point of four shoot ray corresponding to four net points of each angle grid of described angle network topology; E the initial incidence angle of () three shoot ray corresponding to the outgoing triangle at acceptance point place carrys out iteration and asks for correction initial incidence angle; F () is at shot point place, the reflection spot of described shoot ray via destination layer boundary reflection is determined to revise initial incidence angle for geological structure model shoot ray described in original incident angle, shot point, described reflection spot and acceptance point are connected successively, to obtain three-dimensional ray tracing path.

Alternatively, the focus being distributed in earth's surface is set to acceptance point, the geophone station be distributed in well is set to shot point to perform step (b) to (f).

Alternatively, described geological structure model is also included in the bed interface on destination layer interface, wherein, step (c) comprising: determine that each shoot ray is via the eye point in receiving area after the bed interface transmission on destination layer interface and destination layer boundary reflection according to the initial incidence angle of each shoot ray, step (f) comprising: at shot point place, determine that described shoot ray is via the transmission point of the bed interface transmission on destination layer interface and the reflection spot via destination layer boundary reflection to revise initial incidence angle for geological structure model shoot ray described in original incident angle, by shot point, described transmission point, described reflection spot and acceptance point connect successively, to obtain three-dimensional ray tracing path.

Alternatively, step (e) comprising: (e1) is according to the outgoing triangle T of acceptance point P at its place ₁t ₂t ₃in area coordinate and described outgoing triangle T ₁t ₂t ₃the initial incidence angle of three corresponding shoot ray, calculates new initial incidence angle; (e2) at shot point place, with new initial incidence angle for geological structure model shoot ray described in original incident angle, and determine that described shoot ray is via the eye point T in receiving area after destination layer boundary reflection; (e3) determine whether the distance between described eye point T and acceptance point P is less than predetermined threshold; (e4), when the distance between described eye point T and acceptance point P is less than described predetermined threshold, iteration is stopped, and using new initial incidence angle as correction initial incidence angle; (e5), when the distance between described eye point T and acceptance point P is more than or equal to described predetermined threshold, utilize described eye point T by described outgoing triangle T ₁t ₂t ₃be divided into three new triangles, and determine the new triangle at acceptance point P place; (e6) using the new triangle at acceptance point P place as outgoing triangle T ₁t ₂t ₃, return step (e1).

Alternatively, step (e1) comprising: (e11) calculates the direction vector (x of new initial incidence angle in three-dimensional system of coordinate according to following formula _new, y _new, z _new):

x _new＝u ₁x _a+u ₂x _b+u ₃x _c，

y _new＝u ₁y _a+u ₂y _b+u ₃y _c，

z _new＝u ₁z _a+u ₂z _b+u ₃z _c，

Wherein, (x _a, y _a, z _a), (x _b, y _b, z _b), (x _c, y _c, z _c) be respectively the outgoing triangle T at acceptance point P place ₁t ₂t ₃initial incidence angle (the α of three corresponding shoot ray _a, β _a), (α _b, β _b), (α _c, β _c) direction vector in three-dimensional system of coordinate, (u ₁, u ₂, u ₃) for acceptance point P is in the outgoing triangle T at its place ₁t ₂t ₃in area coordinate, wherein, s indicates described outgoing triangle T ₁t ₂t ₃area, S ₁, S ₂, S ₃indicate described outgoing triangle T respectively ₁t ₂t ₃three summits and three little leg-of-mutton areas forming of acceptance point P; (e12) new initial incidence angle is determined according to the direction vector of the new initial incidence angle calculated.

Alternatively, in step (d), according to the area coordinate (u of acceptance point in each outgoing triangle ₁, u ₂, u ₃) determine the outgoing triangle at acceptance point place, wherein, if u ₁>=0, u ₂>=0, u ₃>=0, and | u ₁|+| u ₂|+| u ₃|=1, then determine that acceptance point is in the outgoing triangle of correspondence, wherein, s indicates the leg-of-mutton area of outgoing, S ₁, S ₂, S ₃indicate three little leg-of-mutton areas that leg-of-mutton three summits of outgoing and acceptance point are formed respectively.

Alternatively, α _iscope be 0 ° to 90 °, β _jscope be 0 ° to 360 °.

In the ray-tracing scheme according to an exemplary embodiment of the present invention for the geological structure model of complexity, based on angle network topology carry out ray trial fire and gore area method iteration, can follow the trail of rapidly and accurately each big gun examine between many rays; In addition, can be distributed in earth's surface based on triangle area process of iteration to shot point, the VSP recording geometry that geophone station is distributed in well carries out ray tracing.

Part in ensuing description is set forth general plotting of the present invention other in and/or advantage, some will be clearly by describing, or can learn through the enforcement of general plotting of the present invention.

Accompanying drawing explanation

By below in conjunction with exemplarily illustrating the description that the accompanying drawing of embodiment carries out, the above and other object of exemplary embodiment of the present and feature will become apparent, wherein:

Fig. 1 illustrates according to an exemplary embodiment of the present invention for the process flow diagram of the ray-tracing scheme of the geological structure model of complexity;

Fig. 2 illustrates the schematic diagram of angle network topology according to an exemplary embodiment of the present invention;

Fig. 3 illustrates the schematic diagram of geological structure model according to an exemplary embodiment of the present invention;

Fig. 4 illustrates the schematic diagram of the spatial relationship that ray is crossing with gore according to an exemplary embodiment of the present invention;

Fig. 5 illustrates the schematic diagram of triangle area coordinate according to an exemplary embodiment of the present invention;

Fig. 6 illustrates the schematic diagram of six adjacent triangle of A point according to an exemplary embodiment of the present invention;

Fig. 7 illustrates the leg-of-mutton schematic diagram of outgoing according to an exemplary embodiment of the present invention;

Fig. 8 illustrates that iteration asks for the process flow diagram of the method revising initial incidence angle according to an exemplary embodiment of the present invention;

Fig. 9 illustrates the schematic diagram of the ray tracing of three-dimensional ray tracing method generation according to an exemplary embodiment of the present invention.

Embodiment

Now will in detail with reference to embodiments of the invention, the example of described embodiment is shown in the drawings, and wherein, identical label refers to identical parts all the time.Below by referring to accompanying drawing, described embodiment will be described, to explain the present invention.

Fig. 1 illustrates according to an exemplary embodiment of the present invention for the process flow diagram of the ray-tracing scheme of the geological structure model of complexity.Described method can have been come by the device for ray tracing, also realizes by computer program, thus when running this program, realizes said method.

As shown in Figure 1, in step S10, the angular region of the initial incidence angle of shoot ray is configured to the angle network topology comprising m × n net point, wherein, each net point corresponding angle (α _i, β _j), α _ifor shoot ray in three-dimensional system of coordinate with angle formed by Z axis negative direction, β _jangle formed by the projection of shoot ray in XOY plane and X-axis positive dirction, i is greater than the integer that 0 is less than or equal to m, α _i> α _i-1, j is greater than the integer that 0 is less than or equal to n, β _j> β _j-1.

Fig. 2 illustrates the schematic diagram of angle network topology according to an exemplary embodiment of the present invention.As shown in Figure 2, along the circumferential direction mesh lines instruction α _i, along the mesh lines instruction β of radial direction _j, both intersect the net point corresponding angle (α formed _i, β _j).

In addition, angle network topology also can be other forms, such as, and rectangular angular network topology.If rectangular angular network topology, then described rectangular angular network topology comprise that m is capable, n row, the net point corresponding angle (α on the i-th row jth row _i, β _j), and along with i larger, α _ilarger, along with j is larger, β _jlarger.

Exemplarily, α _iscope can be 0 ° to 90 °, β _jscope can be 0 ° to 360 °.Those skilled in the art also can arrange α according to the balance of efficiency and precision _iand β _jscope, and the density of grid.

In step S20, at shot point place, with the angle (α corresponding to each net point of described angle network topology _i, β _j) fire one group of ray for adjustment for geological structure model described in original incident angle.That is, respectively with the angle corresponding to each net point for initial incidence angle is from shot point to described geological structure model shoot ray.

In step S30, determine that each shoot ray is via the eye point in receiving area after destination layer boundary reflection according to the initial incidence angle of each shoot ray.

Particularly, when there are not other bed interfaces on destination layer interface in geological structure model, shoot ray is propagated to geological structure model from shot point along original incident angular direction, reflect according to Si Nieer (Snell) theorem when running into destination layer interface, then the shoot ray after reflection is propagated to direction, receiving area, arrive receiving area outgoing and namely complete a ray trial fire, and form corresponding eye point.

When also there are other bed interfaces on destination layer interface in geological structure model, shoot ray is propagated to geological structure model from shot point along original incident angular direction, transmission is carried out according to Si Nieer (Snell) theorem when running into other bed interfaces, reflect according to Si Nieer (Snell) theorem when running into destination layer interface, shoot ray after reflection is propagated to direction, receiving area and is carried out transmission when running into other bed interfaces, until arrive receiving area outgoing namely complete a ray trial fire, and form corresponding eye point.

Below, composition graphs 3 is illustrated how to determine the eye point of shoot ray in receiving area.Fig. 3 illustrates the schematic diagram of geological structure model according to an exemplary embodiment of the present invention.As shown in Figure 3, geological structure model can describe with face and block, block is formed by face closure, geological structure model shown in Fig. 3 has six closure bar, 27 bed interfaces (hereinafter referred to as sub-face), each closure bar is represented respectively with I, II, III, IV, V, VI, represent each sub-face of each closure bar with digital l to 27, wherein, sub-face 13 is destination layer interface.In each closure bar, medium is uniform, that is, the density in closure bar is identical, and ray has identical velocity of propagation in closure bar.Ray is rectilinear propagation in closure bar, and when running into the border of block and block, ray will occur to reflect or transmission, that is, reflect when running into destination layer interface 13, when running into other bed interfaces, transmission occurs.

Particularly, the tracking of the eye point T from shot point S to ray can be carried out as follows: the incidence from shot point S of (1) ray, determine that shot point S is in closure bar I, and according to the velocity of propagation of closure bar I correspondence and density parameter, and the sub-face 1 that the initial incidence angle determination ray of described ray intersects in closure bar I, and obtain intersecting point coordinate; (2) split closure bar I and closure bar II according to sub-face 1, determine that ray passes sub-face 1 and enters closure bar II from closure bar I transmission, and according to closure bar I and the velocity of propagation of closure bar II correspondence and the transmission direction of density parameter determination ray; (3) according to the transmission direction of the ray determined, determine the sub-face 2 that ray intersects in closure bar II, and obtain intersecting point coordinate; (4) split closure bar II and closure bar III according to sub-face 2, determine that ray passes sub-face 2 and enters closure bar III from closure bar II transmission, and according to closure bar II and the velocity of propagation of closure bar III correspondence and the transmission direction of density parameter determination ray; (5) according to the transmission direction of the ray determined, determine the sub-face 10 that ray intersects in closure bar III, and obtain intersecting point coordinate; (6) split closure bar III and closure bar IV according to sub-face 10, determine that ray passes sub-face 10 and enters closure bar IV from closure bar III transmission, and according to closure bar III and the velocity of propagation of closure bar IV correspondence and the transmission direction of density parameter determination ray; (7) according to the transmission direction of the ray determined, determine the sub-face 13 that ray intersects in closure bar IV, and obtain intersecting point coordinate; (8) because sub-face 13 is destination layer interface, calculate the reflection direction of ray, and determine with by the sub-face 5 of ray intersection reflected and intersecting point coordinate; (9) split closure bar IV and closure bar II according to sub-face 5, determine that ray passes sub-face 5 and enters closure bar II from closure bar IV transmission, and according to closure bar IV and the velocity of propagation of closure bar II correspondence and the transmission direction of density parameter determination ray; (10) according to the transmission direction of the ray determined, determine the sub-face 1 that ray intersects in closure bar II, and obtain intersecting point coordinate; (11) split closure bar II and closure bar I according to sub-face 1, determine that ray passes sub-face 1 and enters closure bar I from closure bar II transmission, and according to closure bar II and the velocity of propagation of closure bar I correspondence and the transmission direction of density parameter determination ray; (12) according to the transmission direction of the ray determined, determine the sub-face 18 that ray intersects in closure bar I, and obtain intersecting point coordinate, i.e. the coordinate of the eye point of ray.

Exemplarily, can according to transmission direction and the reflection direction calculating ray for the transmitting formula at three-dimensional fluctuating interface and Reflection formula.

The concrete mode of transmission direction calculating ray according to the transmitting formula for three-dimensional fluctuating interface is as follows: direction vector when setting ray incident is as d, d ₁, d ₂be respectively the durection component of d, the normal vector of point of intersection is n, the radioparent direction vector angle that to be r, θ be between r and n,

$r=cos\mathrm{\θ}\·\frac{d_1}{\left|d_1\right|}+sin\mathrm{\θ}\·\frac{d_2}{\left|d_2\right|}---\left(1\right),$

The concrete mode of reflection direction calculating ray according to the Reflection formula for three-dimensional fluctuating interface is as follows: direction vector when setting ray incident is as d, d ₁, d ₂be respectively the durection component of d, the normal vector of point of intersection is n, and the direction vector of reflected ray is f,

f＝d-2n·d(2)，

Calculate transmission direction and the reflection direction of ray by the way, can avoid solving multivariate linear equations, calculate accurately simple, be conducive to realizing rapid ray tracing.

About the coordinate of intersection point calculating ray and sub-face, because each sub-face can be made up of at least one gore, each gore is determined by three reference mark, and therefore, the intersection point calculating ray and sub-face is the intersection point of the gore calculated on ray and sub-face.In fact the intersection point calculating the gore on ray and sub-face comprises the problem of two aspects: one is the intersection point calculating ray and intersecting triangles face; Two is judge that whether institute's find intersection is the real intersection point (that is, whether institute's find intersection is positioned at this leg-of-mutton inside (comprising border)) of ray and gore.

About the intersection point calculation of ray and gore, Fig. 4 illustrates the schematic diagram of the spatial relationship that ray is crossing with gore according to an exemplary embodiment of the present invention, and as shown in Figure 4, ray is space line, known initial coordinate point P ₀(x ₀, y ₀, z ₀), the direction number of ray is respectively l, m, n, then ray equation can be expressed as:

$\left\{\begin{array}{c}x=x_0+lt\\ y=y_0+mt\\ z=z_0+nt\end{array}\right.---\left(3\right),$

Wherein, t is variable.For gore, three fixed point coordinate P of known gore ₁(x ₁, y ₁, z ₁), P ₂(x ₂, y ₂, z ₂), P ₃(x ₃, y ₃, z ₃), then consisting of the Representation Equation of gore be:

$\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(4\right),$

Now suppose that ray is crossing with gore, joining is Q (x _q, y _q, z _q), because intersection point is the point on ray, be again the point on gore, substitute into formula (3) respectively and formula (4) has:

$\left\{\begin{array}{c}x=x_0+lt\\ y_q=y_0+mt\\ z=z_0+nt\end{array}\right.---\left(3\right),$

$\left|\begin{array}{ccc}{x}_q-x_1& y_q-y_1& z_q-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(6\right),$

Formula (5) substitution formula (6) can be obtained the functional equation about variable t:

$\left|\begin{array}{ccc}{x}_0+lt-x_1& y_0+mt-y_1& z_0+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(7\right),$

For simplifying above formula equation, order:

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

Replace by formula (8), formula (7) becomes:

$\left|\begin{array}{ccc}lt-x_{10}& mt-y_{10}& 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(9\right),$

Launch above formula, the calculation expression that can obtain variable t is:

$t=\frac{a_1{x}_{10}+a_2{y}_{10}+a_3{y}_{10}}{{\mathrm{la}}_1+{\mathrm{ma}}_2+{\mathrm{na}}_3}---\left(10\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(11\right),$

Try to achieve t by formula (10), t is substituted into the intersection point Q (x that formula (5) can calculate ray and gore _q, y _q, z _q) coordinate figure.

The intersection point more than calculated only represents the intersection point of the plane at ray and triangle place, and whether intersection point also needs in triangle is determined further.Can judge intersection point whether in triangle according to triangle area coordinate, Fig. 5 illustrates the schematic diagram of triangle area coordinate according to an exemplary embodiment of the present invention, as shown in Figure 5, supposes that leg-of-mutton three coordinate points are P ₁(x ₁, y ₁, z ₁), P ₂(x ₂, y ₂, z ₂), P ₃(x ₃, y ₃, z ₃), intersection point is Q (x _q, y _q, z _q), form a large triangle Δ P like this ₁p ₂p ₃with three little triangles: Δ P ₁p ₂q, Δ P ₂p ₃q and Δ P ₃p ₁q.Large triangle Δ P ₁p ₂p ₃area is designated as S, and three little triangles are designated as S respectively ₁, S ₂and S ₃.Corresponding leg-of-mutton area value can be calculated respectively according to coordinate, be 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_q)+x_2(y_q-y_1)+x_q(y_1-y_2)\]\\ S_2=\frac{1}{2}\[x_2(y_3-y_q)+x_2(y_q-y_2)+x_q(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(12\right),$

Utilize above formula to calculate triangle area, if coordinate points by counterclockwise result be on the occasion of, if coordinate points in the direction of the clock result be negative value.For triangle area coordinate, at triangle Δ P ₁p ₂p ₃in every bit can represent with triangle area coordinate, for a Q, triangle area coordinate is expressed as (u ₁, u ₂, u ₃), each component is defined as follows:

$u_1=\frac{S_1}{S},u_2=\frac{S_2}{S},u_3=\frac{S_3}{S}---\left(13\right),$

Define from formula (13), three little triangle areas that three representation in components are split by Q point and the ratio of former triangle area, if fruit dot is at triangle interior (comprising border), then have: u ₁>=0, u ₂>=0, u ₃>=0 and | u ₁|+| u ₂|+| u ₃|=1, if Q point is at triangular exterior, then have: | u ₁|+| u ₂|+| u ₃| > 1.Thus, according to leg-of-mutton coordinate, ray and triangle intersection coordinate conversion are become corresponding triangle area coordinate, the character of recycling triangle area coordinate, judge intersection point whether in triangle.If intersection point is in triangle, then try to achieve the real intersection point in ray and sub-face, otherwise continued to find ray and the leg-of-mutton intersection point of the next one, until find out real intersection point.

The intersection point calculating ray and bed interface is very crucial to three-dimensional ray tracing, especially for the geological structure model of labyrinth, gore quantity included by bed interface is comparatively large, and calculates intersection point at every turn and all to need and all triangles compare, and causes calculated amount very large.Can counting yield be directly connected to the key problem that meet actual use, preferably, calculates the intersection point of the gore on ray and bed interface, to improve the speed of ray tracing by following layers interface subregion algorithm.Concrete steps are as follows:

First step: all triangles at constituting layer interface are divided into two sub regions (such as, on average dividing by triangle number) by triangle number.Calculate X, Y, Z coordinate extreme value of two sub regions under three-dimensional system of coordinate XYZ respectively (namely, maximal value in the absolute value of maximum coordinate value among and the absolute value of min coordinates value), X, Y, Z coordinate extreme value according to calculating constructs two space rectangular parallelepipeds respectively.Judge whether ray and bed interface intersect, namely judge whether the space rectangular parallelepiped of ray and structure intersects (that is, judging whether two space rectangular parallelepipeds of ray and structure intersect respectively).Particularly, ray and space rectangular parallelepiped are projected to XOY, XOZ and YOZ plane respectively, the asking of ray and space rectangular parallelepiped hands over judgement to be reduced to judge whether straight line and plane intersect, as long as meet ray and the projection of space rectangular parallelepiped on three plane X OY, XOZ and YOZ all intersects each other, just can determine that the subregion that ray and this space rectangular parallelepiped are corresponding exists intersection point, otherwise the subregion that ray is corresponding with this space rectangular parallelepiped is non-intersect.

Second step: the subregion intersected determined is divided into two sub regions again, the mode (coordinates computed extreme value, structure space rectangular parallelepiped etc.) identical with first step can be adopted, and judge whether ray and two sub regions again divided intersect respectively.

Third step: repeat second step, until when the triangle number in the subregion divided is less than specific threshold (such as, integer between specific threshold desirable 5 ~ 10, preferably desirable 8) time, no longer this subregion is divided, directly find intersection calculating is carried out to each triangle in ray and this subregion, obtain intersection point.Now only need carry out a small amount of find intersection to calculate, significantly can reduce calculated amount.

The essence of bed interface subregion algorithm is that bed interface is constantly divided into less bed interface, and the triangle of intersection point may be had progressively to be locked in a very little scope, and then carries out ray and leg-of-mutton find intersection calculates.Utilize bed interface subregion algorithm can get rid of rapidly in a large number without the triangle of intersection point, significantly can reduce calculated amount thus.

In triangle, the normal vector of arbitrfary point is very important for the reflection direction and projecting direction following the trail of ray, traditional way is the normal vector replacing triangle interior arbitrfary point with leg-of-mutton normal vector, if stratal surface fluctuation ratio is larger, then triangulation method vector will in nonlinearities change, causes normal vector that curved surface is put discontinuous.Therefore, the method effect in the ray tracing of relief surface calculates is not fine.Preferably, in bed interface, the normal vector of each point is obtained by this normal vector weighted interpolation of adjacent triangle in this bed interface, and the area of weight and adjacent triangle is directly proportional, and adjacent triangle center is inversely proportional to the distance of this point.Fig. 6 illustrates the schematic diagram of six adjacent triangle of A point according to an exemplary embodiment of the present invention, in bed interface as shown in Figure 6, adjacent 1 ~ No. 6 triangle of summit A, adjacent 3, No. 4 triangles of F point, adjacent 4, No. 5 triangles of E point etc. (more than two of the triangle that E point is adjacent in a model with F point).If six that summit A is adjacent leg-of-mutton planar process vectors are n _i(i=1,2 ... 6), area is designated as s _i, center is designated as d to the distance of A point _i, then the normal vector of summit A is:

$n_A=\underset{i=1}{\overset{6}{\Σ}}\left(\frac{s_j{n}_j}{d_i}\right)---\left(14\right),$

Other each summit, such as, the normal vector of summit B and C can be obtained by same method.

For the arbitrfary point N of triangle interior, this some area coordinate is in the triangles u _i, then the normal vector of N point is obtained by the normal vector linear interpolation on three summits,

$n_N=\underset{i=1}{\overset{3}{\Σ}}\left(u_i{n}_i\right)---\left(15\right),$

The normal vector of trying to achieve like this is single order continually varying in triangle, leg-of-mutton border is also continuous print.The normal vector such as calculating each point on AB limit respectively by triangle ABC and ABD is identical.In same bed interface, normal vector is single order continuous print, and on different borders, bed interface, normal vector is not necessarily continuous.This also tallies with the actual situation: vertical interface and the horizontal interface junction normal vector of model boundary should be not continuous.

Return Fig. 1, in step S40, determine the outgoing triangle at acceptance point place, wherein, two the non-overlapping triangles of described outgoing triangle for being formed for summit with the eye point of four shoot ray corresponding to four net points of each angle grid of described angle network topology.

In other words, four eye points of four shoot ray being initial incidence angle with the angle corresponding to four net points of same angle grid by formation two outgoing triangles, thus based on angle network topology carry out ray trial fire will form multiple outgoing triangle.Fig. 7 illustrates the leg-of-mutton schematic diagram of outgoing according to an exemplary embodiment of the present invention.As shown in Figure 7, diagonally the quadrilateral partition that four corresponding for same angle grid eye points are formed can be become two outgoing triangles.

Exemplarily, can according to the area coordinate (u of acceptance point in each outgoing triangle ₁, u ₂, u ₃) determine the outgoing triangle at acceptance point place, wherein, if u ₁>=0, u ₂>=0, u ₃>=0, and | u ₁|+| u ₂|+| u ₃|=1, then determine that acceptance point is in the outgoing triangle of correspondence, wherein, s indicates the leg-of-mutton area of outgoing, three little leg-of-mutton areas that S1, S2, S3 indicate leg-of-mutton three summits of outgoing and acceptance point to form respectively.Aforementionedly concrete account form to be described in detail, not to repeat them here.

In step S50, the initial incidence angle of three shoot ray corresponding to the outgoing triangle at acceptance point place carrys out iteration and asks for correction initial incidence angle.

Fig. 8 illustrates that iteration asks for the process flow diagram of the method revising initial incidence angle according to an exemplary embodiment of the present invention.Shown in Fig. 8, in step S501, according to the outgoing triangle T of acceptance point P at its place ₁t ₂t ₃in area coordinate and described outgoing triangle T ₁t ₂t ₃the initial incidence angle of three corresponding shoot ray, calculates new initial incidence angle.

Exemplarily, the direction vector (x of new initial incidence angle in three-dimensional system of coordinate can be calculated according to following formula _new, y _new, z _new):

x _new＝u ₁x _a+u ₂x _b+u ₃x _c(16)，

y _new＝u ₁y _a+u ₂y _b+u ₃y _c(17)，

z _new＝u ₁z _a+u ₂z _b+u ₃z _c(18)，

Wherein, (x _a, y _a, z _a), (x _b, y _b, z _b), (x _c, y _c, z _c) be respectively the outgoing triangle T at acceptance point P place ₁t ₂t ₃initial incidence angle (the α of three corresponding shoot ray _a, β _a), (α _b, β _b), (α _c, β _c) direction vector in three-dimensional system of coordinate, (u ₁, u ₂, u ₃) for acceptance point P is in the outgoing triangle T at its place ₁t ₂t ₃in area coordinate, wherein, s indicates described outgoing triangle T ₁t ₂t ₃area, S ₁, S ₂, S ₃indicate described outgoing triangle T respectively ₁t ₂t ₃three summits and three little leg-of-mutton areas forming of acceptance point P.

Then, new initial incidence angle is determined according to the direction vector of the new initial incidence angle calculated.Here, the direction vector of initial incidence angle (α, β) in three-dimensional system of coordinate (x, y, z) with the pass of (α, β) is: x=sin α cos β, y=sin α sin β, z=cos α.Said method obtains new initial incidence angle by carrying out correction to durection component, and in addition, also directly can carry out correction to obtain new initial incidence angle to angle, the present invention is not restricted this.

In step S502, at shot point place, with new initial incidence angle for geological structure model shoot ray described in original incident angle, and determine that described shoot ray is via the eye point T in receiving area after destination layer boundary reflection.

In step S503, determine whether the distance between described eye point T and acceptance point P is less than predetermined threshold.

Wherein, when determining that the distance between described eye point T and acceptance point P is less than described predetermined threshold, in step S504, stop iteration, and using new initial incidence angle as correction initial incidence angle.

When determining that the distance between described eye point T and acceptance point P is more than or equal to described predetermined threshold, in step S505, utilize described eye point T by described outgoing triangle T ₁t ₂t ₃be divided into three new triangles, and determine the new triangle at acceptance point P place.

In step S506, using the new triangle at acceptance point P place as outgoing triangle T ₁t ₂t ₃, return step S501.

Return Fig. 1, in step S60, when there are not other bed interfaces on destination layer interface in geological structure model, at shot point place, the reflection spot of described shoot ray via destination layer boundary reflection is determined to revise initial incidence angle for geological structure model shoot ray described in original incident angle, shot point, described reflection spot and acceptance point are connected successively, to obtain three-dimensional ray tracing path.

When also there are other bed interfaces on destination layer interface in geological structure model, step S60 can comprise: at shot point place, determine that described shoot ray is via the transmission point of the bed interface transmission on destination layer interface and the reflection spot via destination layer boundary reflection to revise initial incidence angle for geological structure model shoot ray described in original incident angle, shot point, described transmission point, described reflection spot and acceptance point are connected successively, to obtain three-dimensional ray tracing path.

Should be appreciated that, due to geological structure model, raypath non-linear, outgoing triangle may overlap, acceptance point may simultaneously in multiple outgoing triangle, therefore, in step s 40, can determine acceptance point simultaneously multiple outgoing triangles, in step S50, the initial incidence angle of three shoot ray respectively corresponding to each outgoing triangle at acceptance point place can carry out iteration and ask for and revise initial incidence angle accordingly (namely, utilize respectively each outgoing triangle come iteration ask for revise initial incidence angle accordingly), correspondingly, in step S60, corresponding three-dimensional ray tracing path can be obtained respectively according to the different correction initial incidence angle asked for, thus solve the multiresolution issue of raypath, obtain rapidly and accurately the inspection of big gun between many rays.Fig. 9 illustrates the schematic diagram of the ray tracing of three-dimensional ray tracing method generation according to an exemplary embodiment of the present invention.As shown in Figure 9, according to an exemplary embodiment of the present invention three-dimensional ray tracing method can obtain same big gun inspection between many rays.

Under normal circumstances, those skilled in the art using focus as shot point, using geophone station as acceptance point, but due to well be an infinitesimal line in three dimensions in logic, the ray-tracing scheme based on triangle area process of iteration is caused to be difficult to be distributed in earth's surface at shot point, apply when geophone station is distributed in well, the present invention considers the reversibility of ray, namely, even if the shot point of each ray and geophone station switch, the track of ray also can not change, according to exemplary embodiment of the present invention, the focus being distributed in earth's surface is set to acceptance point, the geophone station be distributed in well is set to shot point to perform step S20 to S60.Thus solve iteration in well and can not meet the problem of the pacing items of triangle area process of iteration, triangle area process of iteration can be applied to and be distributed in earth's surface to shot point, geophone station is distributed in the ray tracing of the VSP recording geometry in well.

According to an exemplary embodiment of the present invention for the ray-tracing scheme of the geological structure model of complexity, can follow the trail of rapidly and accurately each big gun inspection between many rays; In addition, can be distributed in earth's surface based on triangle area process of iteration to shot point, the VSP recording geometry that geophone station is distributed in well carries out ray tracing.

Although show and described exemplary embodiments more of the present invention, but those skilled in the art should understand that, when not departing from by the principle of the present invention of claim and its scope of equivalents thereof and spirit, can modify to these embodiments.

Claims (7)

1., for a ray-tracing scheme for the geological structure model of complexity, comprising:

A the angular region of the initial incidence angle of shoot ray is configured to the angle network topology comprising m × n net point by (), wherein, and each net point corresponding angle (α _i, β _j), α _ifor shoot ray in three-dimensional system of coordinate with angle formed by Z axis negative direction, β _jangle formed by the projection of shoot ray in XOY plane and X-axis positive dirction, i is greater than the integer that 0 is less than or equal to m, α _i> α _i-1, j is greater than the integer that 0 is less than or equal to n, β _j> β _j-1;

(b) at shot point place, with the angle (α corresponding to each net point of described angle network topology _i, β _j) fire one group of ray for adjustment for geological structure model described in original incident angle;

C according to the initial incidence angle of each shoot ray, () determines that each shoot ray is via the eye point in receiving area after destination layer boundary reflection;

D () determines the outgoing triangle at acceptance point place, wherein, two the non-overlapping triangles of described outgoing triangle for being formed for summit with the eye point of four shoot ray corresponding to four net points of each angle grid of described angle network topology;

E the initial incidence angle of () three shoot ray corresponding to the outgoing triangle at acceptance point place carrys out iteration and asks for correction initial incidence angle;

F () is at shot point place, the reflection spot of described shoot ray via destination layer boundary reflection is determined to revise initial incidence angle for geological structure model shoot ray described in original incident angle, shot point, described reflection spot and acceptance point are connected successively, to obtain three-dimensional ray tracing path.

2. method according to claim 1, wherein, is set to acceptance point by the focus being distributed in earth's surface, the geophone station be distributed in well is set to shot point to perform step (b) to (f).

3. method according to claim 1, wherein, described geological structure model is also included in the bed interface on destination layer interface, wherein, step (c) comprising: determine that each shoot ray is via the eye point in receiving area after the bed interface transmission on destination layer interface and destination layer boundary reflection according to the initial incidence angle of each shoot ray

Step (f) comprising: at shot point place, determine that described shoot ray is via the transmission point of the bed interface transmission on destination layer interface and the reflection spot via destination layer boundary reflection to revise initial incidence angle for geological structure model shoot ray described in original incident angle, shot point, described transmission point, described reflection spot and acceptance point are connected successively, to obtain three-dimensional ray tracing path.

4. method according to claim 1, wherein, step (e) comprising:

(e1) according to the outgoing triangle T of acceptance point P at its place ₁t ₂t ₃in area coordinate and described outgoing triangle T ₁t ₂t ₃the initial incidence angle of three corresponding shoot ray, calculates new initial incidence angle;

(e2) at shot point place, with new initial incidence angle for geological structure model shoot ray described in original incident angle, and determine that described shoot ray is via the eye point T in receiving area after destination layer boundary reflection;

(e3) determine whether the distance between described eye point T and acceptance point P is less than predetermined threshold;

(e4), when the distance between described eye point T and acceptance point P is less than described predetermined threshold, iteration is stopped, and using new initial incidence angle as correction initial incidence angle;

(e5), when the distance between described eye point T and acceptance point P is more than or equal to described predetermined threshold, utilize described eye point T by described outgoing triangle T ₁t ₂t ₃be divided into three new triangles, and determine the new triangle at acceptance point P place;

(e6) using the new triangle at acceptance point P place as outgoing triangle T ₁t ₂t ₃, return step (e1).

5. method according to claim 4, wherein, step (e1) comprising:

(e11) direction vector (x of new initial incidence angle in three-dimensional system of coordinate is calculated according to following formula _new, y _new, z _new):

x _new＝u ₁x _a+u ₂x _b+u ₃x _c，

y _new＝u ₁y _a+u ₂y _b+u ₃y _c，

z _new＝u ₁z _a+u ₂z _b+u ₃z _c，

Wherein, (x _a, y _a, z _a), (x _b, y _b, z _b), (x _c, y _c, z _c) be respectively the outgoing triangle T at acceptance point P place ₁t ₂t ₃initial incidence angle (the α of three corresponding shoot ray _a, β _a), (α _b, β _b), (α _c, β _c) direction vector in three-dimensional system of coordinate,

(u ₁, u ₂, u ₃) for acceptance point P is in the outgoing triangle T at its place ₁t ₂t ₃in area coordinate, wherein, s indicates described outgoing triangle T ₁t ₂t ₃area, S ₁, S ₂, S ₃indicate described outgoing triangle T respectively ₁t ₂t ₃three summits and three little leg-of-mutton areas forming of acceptance point P;

(e12) new initial incidence angle is determined according to the direction vector of the new initial incidence angle calculated.

6. method according to claim 1, wherein, in step (d),

According to the area coordinate (u of acceptance point in each outgoing triangle ₁, u ₂, u ₃) determine the outgoing triangle at acceptance point place, wherein, if u ₁>=0, u ₂>=0, u ₃>=0, and | u ₁|+| u ₂|+| u ₃|=1, then determine that acceptance point is in the outgoing triangle of correspondence,

Wherein, s indicates the leg-of-mutton area of outgoing, S ₁, S ₂, S ₃indicate three little leg-of-mutton areas that leg-of-mutton three summits of outgoing and acceptance point are formed respectively.

7. method according to claim 1, wherein, α _iscope be 0 ° to 90 °, β _jscope be 0 ° to 360 °.