CN101533102B - Triangular mesh ray tracing global method of two-dimensional complex construction - Google Patents

Abstract

The invention relates to a travel tomographic inversion of two-dimensional complex construction and a ray tracing method of wave field calculation based on Maslov ray theory. For two-dimensional complex external geometric boundary and complex space structure within domain, a medium is parameterized by using a triangle model with uniform slowness subblocks, and nodes are set on the boundary of a triangle unit. A traveling wave surface is used for replacing a wave front surface so that the expansion of the traveling wave surface and the transmission of the wave are transferred with the vibration of the triangle unit and the adjacent triangle unit, and can self-transfer backwards when encountering reverse branch without other shrinkage; the secondary source position of each node and the minimum travel time calculation are realized in the process of travel wave surface expansion; and the ray path from the receiving point to the source point is picked by using the travel time and direction information of each node secondary source and by means of the minimum travel time search. The calculation method is applied to any polygonal mesh.

Description

Triangular mesh ray tracing global method of two-dimensional complex construction

Technical field

The present invention relates to a kind of triangular mesh ray tracing global method of two-dimensional complex construction, tomographic inversion and based on the ray-tracing scheme that the wave field of Maslov ray theory calculates when specifically being the walking of a kind of two-dimensional complex construction belongs to the exploration geophysics technical field.

Background technology

Ray theory and ray method are the important channels of understanding wave field propagation law, are the important means of the wave field propagation problem in underground complex structure of research and the inhomogeneous medium, also are the bases of tomographic inversion when walking.

No matter be forward simulation or tomographic inversion, model parameterization is basic problem.Because engineering geology is surveyed the complicacy of field condition, be irregular shape under the most of situation of the external shape of institute's detection study object, as " worker " font structure; Inner geology anomalous body is complicated more on space structure, and what have is irregular column (as solution cavity, karst collapse col umn etc.), have be ribbon (as tomography, crack etc.), have be stratiform (as weak intercalated layer etc.).This proposes higher requirement with regard to model parameterization, surveys with the high resolving power that satisfies the labyrinth body.Model parameterization generally all adopts the rectangle net or the triangulation network at present, and wherein rectangle net subdivision is the most general form, and triangulation network subdivision also mainly is a simple trigonometric ratio on rectangle net basis.These simple rectangle nets, triangulation network subdivision do not take into full account regional complex boundary, also radiographic density, dielectric property are not combined with sizing grid, the ray tracing that is difficult to satisfy labyrinth is just being drilled the requirement with high resolving power chromatography imaging inverting.Must adopt irregular net, as the hybrid network of network of triangle, irregular quadrilateral net or network of triangle and irregular quadrilateral net, to model parameterization.Studies show that with respect to the square network parametrization, the triangle gridding parametrization has following advantage: 1. the dirigibility of model subdivision is strong; 2. to the description precision height of velocity discontinuity, rate pattern and INTERFACE MODEL have consistance; 3. the forward simulation memory space is little, computing time is few; 4. the forward simulation difficulty reduces, error reduces, efficient improves; Memory space is little when 5. being used for tomographic inversion, computing time is few, the equation condition is good, find the solution easy.

The theoretical foundation of ray tracing is that under the high-frequency approximation condition, the main energy of elasticity wave field is propagated along ray tracing.Traditional ray-tracing scheme comprises the bending method (Bending method) of the trial fire method (Shootingmethod) and the boundary value problem of initial-value problem on the ordinary meaning.Vidale (1988) points out that the subject matter of trial fire method and bending method is: 1. be difficult to velocity variations stronger in the treatment media when proposing the method for finite difference of eikonal equation; When the overall minimum in when 2. being difficult to obtain many-valued walk is walked; 3. counting yield is lower; 4. shadow region inner rays coverage density deficiency.

Ray tracing problem for triangulation network model, relevant document, done discussion, but its ray-tracing scheme that provides ranges the trial fire method and the bending method of the trial fire method and the boundary value problem of the initial-value problem on the ordinary meaning substantially, has its inherent defect inevitably.

For the ray tracing problem of rectangle pessimistic concurrency control, relevant document has been done more detailed summary to common ray-tracing scheme.In order to overcome this shortcoming of traditional trial fire method and bending method, in recent years, many geophysical research persons have carried out a large amount of work in this respect, have proposed the FRONT TRACKING method of the calculating first arrival whilst on tour that some precision are higher, efficient is higher and practical.Progress is mainly reflected in: the 1. improvement on the basis of traditional trial fire method and bending method, as all kinds of wavefront reconstruction method (Vinje, 1992; Sun, 1992; Lambare et al., 1996), during except that many-valued walk, also solved counting yield and shadow region preferably and covered not enough problem; The improvement of algorithm when 2. minimum being walked is calculated when making it to adapt to many-valued walk, and as slowness matching method (Symes, 1998), can think the popularization of shortest-path method; The combining of algorithm when 3. classic method is with minimum walk as HWT method (Sava, Fomel, 1998), then is to calculate raypath by front propagation.In these methods, when considering walking on all discrete points of space simultaneously and the global calculation method of raypath,, its higher precision and counting yield pay close attention to because having caused widely.Yet these are primarily aimed at the rule mesh ray-tracing scheme, are applied in the triangulation network and can complicate the issue, even be difficult to realize.

The advantage of the triangulation network makes that model triangulation network parametrization is the parameterized important trend of following geophysical model, but its wave field ray tracing problem lacks deep research, therefore, triangulation network model and wave field ray tracing problem thereof are that the important topic that must solve with inverting efficient, effect is just being drilled in present stage geophysical probing technique raising.

Summary of the invention

At above-mentioned existing problems, the object of the present invention is to provide a kind of Global Algorithm that can carry out ray tracing to the two-dimensional complex construction triangulation network.

Its technical solution is as follows:

Triangular mesh ray tracing global method of two-dimensional complex construction is to replace wavefront surface with the capable face of ripple, the expansion of the capable face of ripple, the transmission of ripple are undertaken by " vibration " transmission that a triangular unit is adjacent triangular unit, meet bow-tie, can transmit backward voluntarily, need not " contraction " separately, the secondary source of each node is determined and minimum is calculated when walking realize that its computing method are carried out according to the following steps in the capable face expansion process of ripple:

A) with the two-dimensional complex construction zone, with medium parameterization, on the border of triangular element, node is set with the uniform triangle model of slowness piecemeal,

Be zero when b) establishing the walking of source node, for infinitely great, include the triangle at source node place in ripple capable face triangular unit array or Priority Queues, and make the triangular unit at source node place be in state of activation during the walking of other nodes except that source node, make the sign that is present in the capable face array of ripple

C) in the capable face triangular unit of ripple array, search minimum in the triangular unit that is in state of activation node when walking; Determine this node place triangular unit; With this triangular unit as current triangular unit; If this node is the common node of a plurality of triangular units; Select one of them triangular unit to get final product; Other triangular unit can travel through voluntarily; Current triangular unit when this minimum is walked and adjacent triangular unit thereof are included the capable face triangular unit of ripple array in; And make the sign that is present in the capable face array of ripple; So that repeatedly do not deposit in the capable face triangular unit of the ripple array

Node when d) minimum of triangular unit is walked when minimum is walked calculates, hourage of each node of triangular unit when more current minimum is walked, and the secondary source of definite each node,

E) calculate again, during relatively more minimum walk the adjacent triangular unit node of triangular unit hourage and determine the secondary source of each node, when when adjacent triangular unit node is walked, changing, this adjacent triangular unit is made activation marker,

F) calculate again, triangular unit node hourage when more current minimum is walked, when each node of triangular unit is walked when current minimum is walked during no change, triangular unit was made " dormancy " sign when current minimum was walked, get back to c) step, when changing when each node of triangular unit is walked when current minimum is walked, triangular unit was made " activation " sign when current minimum was walked, get back to c) step, all triangular units all are in " dormancy " state when all no longer changing when the node of the capable face of all ripples is walked and in the capable face of ripple, and the capable face expansion of ripple finishes.

G), when walking and secondary source information, carry out picking up the raypath from the acceptance point to the source point to the source retrieval according to each node from acceptance point.

Described wavy surface is meant in the two dimensional surface, the plane that each node that is arrived by wave source point and a certain moment ripple constitutes.

Described source node is meant and excites source point.

Described adjacent triangular unit is meant and all relevant triangular units of a triangular unit summit or limit.

Owing to adopted above technical scheme, the present invention is directed to the ray tracing problem of two-dimensional complex construction, proposed to replace the triangular mesh ray tracing global method of wavefront surface with the capable face of ripple, numerical simulation shows that this algorithm has very high precision and reliability, overcome the weakness of the ray-tracing scheme bad adaptability of traditional rule-based net, this method makes rule mesh become a special case in the method, is applicable to the arbitrary polygon net, as network of quadrilaterals, quadrilateral and triangle hybrid network.For three-dimensional problem, transform suitable too a little.

That wavefront is portrayed as a curve is different with some, the characteristics of this method are: with the triangular unit that traveled through as a ripple battle array plane, triangular unit was adjacent triangular unit " vibration " transmission and carries out when on the one hand the transmission of the expansion of the capable face of ripple, ripple was walked by minimum, on the other hand when meeting bow-tie, transmit backward voluntarily, need not " contraction " separately, that is will before make the triangular unit " activation " of " dormancy " sign by the adjacent triangular unit of triangular unit, triangular unit when becoming minimum once more and walking.Like this with regard to the effective problem that has solved the expansion and the passback of ripple.

Two-dimensional complex regional model trigonometric ratio and triangular mesh ray tracing global method of two-dimensional complex construction will promote just drilling research work based on the wave field of ray theory, and make it just drill efficient greatly to improve.

Description of drawings

Accompanying drawing 1 rate pattern parametrization triangular unit

Accompanying drawing 2 gable-top neighborhood of a point triangle summits

The neighborhood point of accompanying drawing 3 nodes wherein, (a) for node is a triangle summit situation, (b) for node is the insertion point situation of three arms of angle, (c) is node situation between the insertion point of three arms of angle, (d) is that node is in the inner situation of triangular unit

The adjacent triangular unit of 4 one triangular units of accompanying drawing

When walking, calculates the wavefront neighborhood point of 5 one triangular units of accompanying drawing

The different approximate comparisons of accompanying drawing 6 wavefront wherein, (a) are the approximate synoptic diagram of point before the discrete wave, (b) are the approximate synoptic diagram of point before the continuous wave

The expansion process of the capable face of accompanying drawing 7 ripples from (a) to (h)

Accompanying drawing 8 each node secondary source calculated examples

The retrieval of accompanying drawing 9 node secondary sources wherein, (a) is node E situation in triangular unit, (b) is node E situation on triangle summit or insertion point, (c)～(d) is that node E is between the insertion point, limit

Accompanying drawing 10 2 flat seam dielectric model result of calculations wherein, (a) are divided result and primary wave raypath figure for triangular unit, are primary wave wave front figure (b), contrast when (c) walking for right margin, contrast when (d) walking for the coboundary

Accompanying drawing 11 horizontal layer model primary wave raypath and wave front figure wherein, (a) are the subdivision grid, (b) are primary wave raypath and wave front

Accompanying drawing 12 contains empty model primary wave raypath and wave front, wherein, (a) is the subdivision grid, (b) is primary wave raypath and wave front

Accompanying drawing 13 contains FAULT MODEL primary wave raypath and wave front, wherein, (a) is the subdivision grid, (b) is primary wave raypath and wave front

Embodiment

Below in conjunction with accompanying drawing and instantiation computing method of the present invention are described in further detail.

See accompanying drawing.

For ease of better understanding technical scheme of the present invention, at first the related notion that triangular mesh ray tracing global method of two-dimensional complex construction relates to is described.

As accompanying drawing 1, be the triangular unit of 2 dimensional region triangulation formation.Regulation: three summits of triangular unit are by counterclockwise arranging, and as A → B → C, three limits are directed line segment, as AB, BC, CA; The speed degree is even in each unit.Identical many nodes are set on 3 limits of each triangular unit, and nodal distance can uniform distribution, also can non-homogeneously distribute, and adopt uniform distribution here, and the node that the triangular unit common edge is provided with is identical.

Institute in the 2 dimensional region is called node a little, comprises source point, acceptance point, and the insertion point of the triangular unit summit and three arms of angle etc. accordingly, are referred to as source node, receiving node, triangular unit summit node and three arms of angle and insert node.

All triangular unit summit nodes adjacent with a summit of triangular unit are called the neighborhood triangle summit of this triangular unit summit node, and its direction is defined as from this point and points to neighborhood triangle summit, and agreement neighborhood triangle summit is around this some discharging counterclockwise.As accompanying drawing 2, P ₀Neighborhood triangle summit be respectively: P ₁, P ₂, P ₃, P ₄, P ₅, P ₆

Claim that any two the adjacent different nodes on triangular element inside or the limit are consecutive point, and arrange its directive property.All consecutive point of a node are called the neighborhood point of this node, and its direction is defined as from this node sensing neighborhood point, and arranges neighborhood point and discharge counterclockwise around this node, for reducing unnecessary neighborhood point, and on one side, a reservation and this node next-door neighbour's node.

For the regional triangulation network, the neighborhood point of node is not considering to mainly contain four kinds of situations under the border condition, as accompanying drawing 3, the neighborhood point rr of node ss (open circles) (all filled circles), that is: node on the triangular unit summit, node in the insertion point of three arms of angle, node between the insertion point of three arms of angle, node is in triangular unit inside.For a triangular unit, for (a) and (b), (c) three kinds of situations in the accompanying drawing 3, the neighborhood point of a node in this triangular unit only needs the neighborhood point of consideration in this triangular unit.Accordingly, the neighborhood of node point may comprise: the point between the insertion point of triangular unit summit, three arms of angle, the insertion point of three arms of angle, this four category node of triangular unit internal point.

Calculating when a known node is walked to the whilst on tour process of its neighborhood point, node is called the wavefront point of its neighborhood point during this known walks.As the ss in the accompanying drawing 3, its neighborhood point is called the wavefront neighborhood point of this node, as the rr in the accompanying drawing 3.The secondary source of node is ripple and propagates into this node by the source and experienced a last neighborhood node on the path during known walks.Therefore, put in the whilst on tour process of its wavefront neighborhood point at the calculating wavefront, each wavefront point can be regarded the secondary source of its wavefront neighborhood point in succession as, but has only when all wavefront points of traversal, and the secondary source of each node is only a last neighborhood node of minimum when this node is walked.

With all relevant triangular units of a triangular unit summit or limit, be called the adjacent triangular unit of this triangular unit.As accompanying drawing 4, triangular unit A ₀Adjacent triangular unit comprise A ₁～A ₁₂12 triangular units.The expansion of the capable face of ripple that will define below and the passback of bow-tie are expanded by adjacent triangular unit and are transmitted.

In two dimensional surface, the plane that each node that a certain moment is arrived by wave source point and ripple constitutes is called the capable face of ripple.In the triangular mesh ray tracing global algorithm, the form of expression of the capable face of ripple is the plane that the triangular unit of ripple traversal is formed.If regional subdivision is the polygon net of other form, as network of quadrilaterals, quadrilateral and triangle hybrid network, then the form of expression of the capable face of ripple is, the quadrilateral units or four deformation units of ripple traversal mix the plane of forming with triangular element.What deserves to be explained is, when calculating the walking of each node in the process, before calculating when the node of finishing whole zone is walked, not necessarily must be that minimum is when walking during the walking of capable each node of face of ripple.

By above-mentioned notion, just can obtain using in the triangular mesh ray tracing algorithm data structures such as more frequent directed edge, node, triangular unit.

The directed edge class

Class Edge // directed edge class

{

public：

Vertex*dest; // the second some position

Double cost; The cost on // limit can be a distance, also can be the time;

Triangle*lefttriangle; The left triangle on // limit

Triangle*righttriangle; The right triangle on // limit

}；

The node class

class?Vertex

{

public：

Double x; // node x coordinate

Double y; // node y coordinate

Dcllist＜Edge〉adj_tri_vex; // this puts adjacent triangle summit

Dcllist＜Edge〉adj_vex; // this vertex neighborhood point,

//dcllist is the double-linked circular list template

Vertex * prev_vex; The secondary source of // this point

Bool scratch; // search sign

Double density; // this triangular unit density control size

Double travel_time; // by the whilst on tour of source point to this node

}；

The triangular unit class

Class Triangle//triangle class

{

public：

Dcllist＜Vertex*〉tripoly; // three square rings that constitute by triangular unit summit and insertion point, limit

Vertex*a, * b, * c; Three summits of // triangular unit

Double v; // triangular unit velocity of wave

Bool scratch; // whether by traversal sign (dormancy, activation marker)

Bool insert; // whether enter ripple battle array plane landmark

Double min_time; When // minimum is walked

DCLLIST::iterator min_time_itr; Position when // minimum is walked,

//DCLLIST::iterator, triangle ring iterative position

Dcllist＜Triangle*〉adj_tri; // adjacent triangular unit

}；

Triangular mesh ray tracing global method of two-dimensional complex construction is to calculate when utilizing the expansion of the capable face of ripple to finish the node minimum to walk and the determining of secondary source, and the expansion of the capable face of ripple is adjacent unit " vibration " transmission by a unit and finishes.

For a complex region, adopt the triangulation method with regional medium parameterization, it is even to look each triangular unit slowness, and some topological relations of the node and the triangulation network are set.The basic step of triangular mesh ray tracing global computing method is as follows:

The first step is calculated when utilizing the expansion of the capable face of ripple to finish the node minimum to walk and secondary source definite.Replace wavefront surface with the capable face of ripple, the expansion of the capable face of ripple, the transmission of ripple are undertaken by " vibration " transmission that triangular unit is adjacent triangular unit, meet bow-tie, can transmit backward voluntarily, need not separately " contraction ".

Second step, the whilst on tour and the secondary source information pickup raypath of each node that utilization calculates.When promptly walking according to the minimum of each node, secondary source, from acceptance point, pick up raypath to the source.

The implementation procedure of triangular mesh ray tracing global method of two-dimensional complex construction is divided following several respects.

1 is provided with the topological relation of each node

With the two-dimensional complex construction zone, with medium parameterization, on the border of triangular element, node is set with the uniform triangle model of slowness piecemeal, the topological relation between each node is set simultaneously, mainly comprise the topological relation of following several category nodes:

(1) by counterclockwise, the neighborhood triangle summit of each node is set, the node here comprises the insertion point of source point, acceptance point, triangle summit and three arms of angle, the left and right triangle of the directed edge that is formed by node to neighborhood triangle summit is set simultaneously, and fundamental purpose is to be convenient to inquire about leg-of-mutton adjacent triangle.

(2) by counterclockwise, the neighborhood point of each node is set, the node here comprises the insertion point of source point, acceptance point, triangle summit and three arms of angle, the left and right triangle of the directed edge that is formed by node to its neighborhood point is set simultaneously, and fundamental purpose is to be convenient to calculate by the whilst on tour of wavefront point to wavefront neighborhood point.It is to make ray tracing to calculate separately source point and acceptance point that acceptance point and source point are equal to the advantage that the insertion point of the triangle summit and three arms of angle is provided with, and the path is picked up and only need be provided acceptance point coordinate or acceptance point numbering and get final product.

(3) by counterclockwise, the adjacent triangular unit of each triangular unit is set, is provided with according to gable-top neighborhood of a point triangle summit, fundamental purpose is to be convenient to the transmission of ripple between the triangular unit and the expansion of the capable face of ripple.

When 2 wavefront neighborhood point minimums are walked and the secondary source position calculation

The definite of calculating and secondary source only needed consideration to carry out in a triangular unit when wavefront neighborhood point minimum was walked, and this step Yun is contained in the expansion of the capable face of ripple, and the capable face of ripple is expanded a part of as follows description of details.

Accompanying drawing 5 is triangular units, and establishing the ss point is the wavefront node, and its wavefront neighborhood point represents that with filled circles wherein the rr node is a certain neighborhood point of ss, and their coordinate is respectively (x _Ss, y _Ss) and (x _Rr, y _Rr); (x _Rs, y _Rs) be the coordinate of rr node secondary source; t _SsDuring for the walking of ss node; t _RrDuring for the walking of rr node; V is the speed of triangular unit; t _SrFor ripple from the ss node with this triangular unit speed v along rectilinear propagation to the required time of rr node:

$t_{\mathrm{sr}}={[{(x_{\mathrm{ss}}-x_{\mathrm{rr}})}^2+{(y_{\mathrm{ss}}-y_{\mathrm{rr}})}^2]}^{\frac{1}{2}}/v---\left(1\right)$

So in the time of can obtaining the walking of rr node and corresponding secondary source node coordinate, promptly as (t _Ss+ t _Sr)＜t _RrThe time, have

$\left\{\begin{array}{c}{t}_{\mathrm{rr}}=t_{\mathrm{ss}}+t_{\mathrm{sr}}\\ x_{\mathrm{rs}}=x_{\mathrm{ss}}\\ y_{\mathrm{rs}}=y_{\mathrm{ss}}\end{array}\right.---\left(2\right)$

Computing node is walked constantly in a triangular unit, the node ss of minimum when this triangular unit is walked, allow ss travel through all nodes of this triangular unit then, calculate when the node of this triangular unit is walked and just finish, but not necessarily be exactly that minimum is when walking when each node is walked, has only when ss experiences all wavefront points t _RrWhen just becoming minimum and walking, (x _Rs, y _Rs) be only real secondary source coordinate.

The secondary source here is approximate with node, during the walking of wavefront neighborhood point be add when walking with the wavefront point wavefront put wavefront neighborhood point part when walking with represent, be equivalent to wavefront surface is regarded as series of discrete point secondary source, as accompanying drawing 6 (a).In fact the ripple that arrives each node might not just in time be that the wavefront point from the grid sends, and may be to send from the point between the net point, as accompanying drawing 6 (b).When the part is walked like this and secondary source position will be revised.

As accompanying drawing 6 (b), when the minimum of computing node E to walk and the position of secondary source, the secondary source of supposing E is between 3 of A, B, C, as the D point, establish certain a bit for local initial point (as the A point), at A, function representation is t (x) when walking between the B, 3 of C, and x is a local one-dimension coordinate on the AC direction.If the local coordinate of some E in the place ahead on the AC direction is x _E, the distance of vertical AC direction is d, when then any point x arrives walking of neighborhood point E on the AC line is

$t_E=t\left(x\right)+{[{(x-x_E)}^2+d^2]}^{\frac{1}{2}}/v---\left(3\right)$

V is the speed in the triangular unit, when making x between AC, change, and t _EGet minimumly, can obtain wavefront point is secondary source point D:

$\{x_D=x:t_E=\underset{x\∈(A,C)}{\min }[t\left(x\right)+{({(x-x_E)}^2+d^2)}^{\frac{1}{2}}]/v\}---\left(4\right)$

If the local coordinate that A, B, C are 3 is x _A, x _B, x _C, be t during 3 walk _A, t _B, t _C, when then upward any point is walked with hyperbolic curve approximate treatment AC

$t\left(x\right)=\left\{t_A^2\right[(x-x_B)(x-x_C)]/[(x_A-x_B)(x_A-x_C)]+$

$t_B^2\left[(x-x_A)(x-x_C)\right]/\left[(x_B-x_A)(x_B-x_C)\right]+---\left(5\right)$

$t_C^2{\left[(x-x_A)(x-x_B)\right]/\left[(x_C-x_A)(x_C-x_B)\right]\}}^{1/2}$

Easy proof, when homogeneous media, t (x) is walking constantly of causing of a certain point source, following formula is strict accurate.When medium was inhomogeneous, hyperbolic curve is approximate also to be that precision is the highest in 3 interpolation methods when walking.(4) formula is that a minimal value is found the solution problem, can adopt methods such as gradient method, golden search to try to achieve x _DAllow A, B, C go through all over the node on other both sides of removing limit, E point place when having finished that the minimum of E point in this triangle walked and the calculating of secondary source position by the limit.

The capable face expansion of 3 ripples

In global approach based on the ray tracing of rectangle net, when often utilizing the expansion of rectangle simulation wave front and shrinking the computing node minimum to walk and secondary source position, but in the triangulation network, the ripple front of drafting a certain moment will be very complicated, and the employing wave front, also be unfavorable for the transmission of bow-tie.

Way of the present invention is, utilizes the ripple parallel planes, and the expansion of ripple parallel planes is considered as a recursive procedure, and its arthmetic statement is as follows:

Be zero when (1) establishing the walking of source point, be infinity during the walking of other nodes except that source point.Include the triangle at source point place in the capable face triangular unit of ripple storage container, and make this triangular unit be in state of activation, make the sign that is present in the capable face of ripple.

(2) in the capable face triangular unit of ripple storage container, search minimum in the triangular unit that is in state of activation node when walking, determine this node place triangular unit, with this triangular unit as current triangular unit; If this node is the common point of a plurality of triangular units, select one of them triangular unit to get final product, other triangular unit can travel through voluntarily.Current triangular unit when this minimum is walked and adjacent triangular unit thereof are included the capable face triangular unit of ripple storage container in, and make the sign that is present in the capable face of ripple, are convenient to repeatedly not deposit in the capable face triangular unit of ripple storage container.

(3) utilize aforementioned " when wavefront neighborhood point minimum is walked and secondary source position calculation ", calculate, triangular unit node whilst on tour when more current minimum is walked (node when this triangular unit minimum is walked, same down); Calculate again, the whilst on tour of the adjacent triangular unit node of triangular unit during relatively more minimum walk, when when adjacent triangular unit node is walked, changing, this adjacent triangular unit is made activation marker; Calculate again, triangular unit node whilst on tour when more current minimum is walked, when each node of triangular unit is walked when current minimum is walked during no change, triangular unit was made " dormancy " sign when current minimum was walked, got back to for (2) step, when changing when each node of triangular unit is walked when current minimum is walked, triangular unit was made " activation " sign when current minimum was walked, and got back to for (2) step; Analogize in proper order, travel through all triangular units.When the node of the capable face of all ripples is walked, all no longer changing and the triangular unit in the capable face of all ripples all be in " dormancy " state, the expansion of the capable face of ripple finishes.

That wavefront is portrayed as a curve is different with some, and the characteristics of these computing method are: the triangular unit that traveled through as a ripple battle array plane, is replaced wavefront surface with the capable face of ripple; The capable face of expansion ripple is transmitted in " vibration " that be adjacent between triangular unit of triangular unit when walking by minimum on the one hand, on the other hand when meeting bow-tie, also incite somebody to action the triangular unit of " dormancy " " activation " before by the transmission between the adjacent triangular unit of triangular unit, and triangular unit when becoming minimum once more and walking need not separately " contraction ".Like this with regard to the effective problem that has solved the expansion and the passback of ripple.

Accompanying drawing 7 has been showed the expansion process of the capable face of ripple from (a) to (h), and source point is in the upper left corner, and the thick line triangular unit is represented current triangular unit.

Accompanying drawing 8 is a calculated examples, and source point demonstrates the secondary source of each node in the upper left corner, is B as the secondary source of A.

Picking up of 4 propagation path to the source

The arbitrfary point when the raypath of source point is walked by the node minimum and secondary source information to the source recourse finish.The present invention handles source point and acceptance point according to node is unified, thereby, in the time in the source path pick process, needn't walking, calculate separately acceptance point, directly recourse gets final product.The retrieval of the secondary source F of accompanying drawing 9 expression node E according to the position of E, has the secondary source retrieval of following several category nodes:

(1) when node E is positioned at triangular unit inside, as accompanying drawing 9 (a), search was found secondary source from triangular unit three limits when approximate and minimum was walked by hyperbolic curve, calculated this node when walking (this category node can be in the capable face expansion of ripple unified the processing, only need retrieval to get final product) here;

(2) when node E is the insertion point of the triangle summit or three arms of angle (9 (b)) can directly pick up secondary source;

(3) when node E is between three arm of angle insertion point A, C, the secondary source F of retrieval E will be according to the secondary source B of A and the secondary source D of C.At this moment mainly can run into as accompanying drawing 9 (c), (d) two kinds of situations: (c) situation be the secondary source B of A and C secondary source D in a triangular unit; (d) situation is with regard to more complicated, because the C point is the triangle summit, the secondary source D of C may be present in other triangular unit.To accompanying drawing 9 (c), (d) two kinds of situations, search was sought secondary source from the triangular unit corresponding edge when also approximate and minimum was walked by hyperbolic curve, and calculated this node when walking.

The check of 5 arithmetic accuracy

For precision and the reliability of checking the triangular mesh ray tracing algorithm, the same notional result of calculating according to the Snell theorem of the triangular mesh ray tracing Forward modelling result of accompanying drawing 10 (a) institute representation model is contrasted.The source is positioned at the left margin upper end, and acceptance point lays respectively at model coboundary and right margin, and acceptance point is apart from 2m.The 1st layer of medium velocity is 2000m/s, and the 2nd layer of medium velocity is 4000m/s, and model horizontal direction width is 50m, and the vertical direction degree of depth is respectively the 1st layer of 10m, the 2nd layer of 40m.When model being carried out the network division, employing density is that the function of 2m carries out triangulation.Accompanying drawing 10 (a) is triangulation network subdivision result and primary wave raypath figure (thick line); Accompanying drawing 10 (b) is primary wave wave front figure, and time line is 1ms at interval.Accompanying drawing 10 (c), (d) are respectively right margin and coboundary comparing result when walking, when solid line is walked for the theoretical primary wave that calculates with the analytical analysis formula among the figure, when round dot is walked for the primary wave of triangulation network tracing algorithm calculating.Relative error was 0.0009% when the right margin maximum was walked; Relative error only was 0.00028% when the coboundary maximum was walked.This shows, adopt the triangular mesh ray tracing algorithm to have quite high computational accuracy.

Specific embodiment

Embodiment one horizontal layer uniform dielectric model

Accompanying drawing 11 is the triangulation and the ray tracing result of horizontal layer model.Source point is positioned at the upper left corner, and acceptance point lays respectively at model coboundary, right margin and lower boundary.Model is shown in accompanying drawing 11 (a), the 1st layer of medium velocity is 2000m/s, the 2nd layer of medium velocity is 4000m/s, the 3rd layer of medium velocity is 3000m/s, the 4th layer of medium velocity is 5000m/s, model horizontal direction width is 50m, and 4 layer thicknesses of vertical direction are 10m, and it is that the function of 2m carries out triangulation that model is adopted density.Accompanying drawing 11 (a) is the triangulation result, different constantly wave front and the raypaths of accompanying drawing 11 (b) for adopting this algorithm to obtain.

Embodiment two contains empty model

Accompanying drawing 12 is for containing the triangulation and the ray tracing result of empty model.Source point is positioned at the upper left corner, and acceptance point lays respectively at model coboundary, right margin and lower boundary.Model is shown in accompanying drawing 12 (a), and background velocity is 4000m/s, is provided with: speed is the low speed rectangle region of 2000m/s, high speed rectangle region that speed is 5000m/s and the cavity of tunnel type shape, and mesh generation adopts the nonhomogeneous density control function.Accompanying drawing 12 (a) is the triangulation result, different constantly wave front and the raypaths of accompanying drawing 12 (b) for adopting this algorithm to obtain.

Embodiment three contains FAULT MODEL

Accompanying drawing 13 is for containing the triangulation and the ray tracing result of FAULT MODEL.Source point is positioned at the lower left corner, and acceptance point lays respectively at the model outer boundary.Model is shown in accompanying drawing 13 (a), is provided with tomography subdomain Ω in the zone ₂And low speed subdomain Ω ₅, be simulation Ω ₃To Ω ₅The graded that speed is bigger increases mesh refinement subdomain Ω ₄, mesh generation has adopted the nonhomogeneous density control function, the analog rate Ω of each subdomain ₁: v _p=3000m/s, Ω ₂: v _p=2000m/s, Ω ₃: v _p=4000m/s, Ω ₄: v _p=4000m/s, Ω ₅: v _p=3000m/s.Accompanying drawing 13 (a) is the triangulation result, different constantly wave front and the raypaths of accompanying drawing 13 (b) for adopting this algorithm to obtain.

Claims (1)

1. triangular mesh ray tracing global method of two-dimensional complex construction, it is characterized in that: triangular mesh ray tracing global method of two-dimensional complex construction is to replace wavefront surface with the capable face of ripple, the expansion of the capable face of ripple, the transmission of ripple are undertaken by " vibration " transmission that a triangular unit is adjacent triangular unit, meet bow-tie, can transmit backward voluntarily, need not " contraction " separately, the secondary source of each node is determined and minimum is calculated when walking realize that its computing method are carried out according to the following steps in the capable face expansion process of ripple:

A) with the two-dimensional complex construction zone, with medium parameterization, on the border of triangular element, node is set with the uniform triangle model of slowness piecemeal, the topological relation between each node is set simultaneously,

Be zero when b) establishing the walking of source node, for infinitely great, include the triangle at source node place in ripple capable face triangular unit array or Priority Queues, and make the triangular unit at source node place be in state of activation during the walking of other nodes except that source node, make the sign that is present in the capable face array of ripple

C) in the capable face triangular unit of ripple array, search minimum in the triangular unit that is in state of activation node when walking; Determine this node place triangular unit; With this triangular unit as current triangular unit; If this node is the common node of a plurality of triangular units; Select one of them triangular unit to get final product; Other triangular unit can travel through voluntarily; Current triangular unit when this minimum is walked and adjacent triangular unit thereof are included the capable face triangular unit of ripple array in; And make the sign that is present in the capable face array of ripple; So that repeatedly do not deposit in the capable face triangular unit of the ripple array

Node when d) minimum of triangular unit is walked when minimum is walked calculates, hourage of each node of triangular unit when more current minimum is walked, and the secondary source of definite each node, and computing formula is as follows:

In the triangular unit, establishing the ss point is the wavefront node when minimum is walked, and the rr node is a certain neighborhood point of ss, and their coordinate is respectively (x _Ss, y _Ss) and (x _Rr, y _Rr); (x _Rs, y _Rs) be the coordinate of rr node secondary source; t _SsDuring for the walking of ss node; t _RrDuring for the walking of rr node; V is the speed of triangular unit; t _SrFor ripple from the ss node with this triangular unit speed v along rectilinear propagation to the required time of rr node:

$t_{\mathrm{sr}}={[{(x_{\mathrm{ss}}-x_{\mathrm{rr}})}^2+{(y_{\mathrm{ss}}-y_{\mathrm{rr}})}^2]}^{\frac{1}{2}}/v---\left(1\right)$

So in the time of can obtaining the walking of rr node and corresponding secondary source node coordinate, promptly as (t _Ss+ t _Sr)＜t _RrThe time, have

$\left\{\begin{array}{c}{t}_{\mathrm{rr}}=t_{\mathrm{ss}}+t_{\mathrm{sr}}\\ x_{\mathrm{rs}}=x_{\mathrm{ss}}\\ y_{\mathrm{rs}}=y_{\mathrm{ss}}\end{array}\right.---\left(2\right)$

Computing node is walked constantly in a triangular unit, the node ss of minimum when this triangular unit is walked, allow ss travel through all nodes of this triangular unit then, calculate when the node of this triangular unit is walked and just finish, but not necessarily be exactly that minimum is when walking when each node is walked, has only when ss experiences all wavefront points t _RrWhen just becoming minimum and walking, (x _Rs, y _Rs) be only real secondary source coordinate; The node approximate representation of the secondary source here, the travel timetable of wavefront neighborhood point is shown and adds when the wavefront point is walked that wavefront puts wavefront neighborhood point part sum when walking, be equivalent to wavefront surface is regarded as series of discrete point secondary source, if the wavefront point that in fact arrives the ripple of each node and just in time be not from the grid sends, but send from the point between the net point, when the part is walked like this and secondary source position (3) and (4) is as follows revised;

When the minimum of calculating the node E of triangular element on is on one side walked and the position of secondary source, the secondary source of supposing E is at triangular element three adjacent node A, B, the D point between the C on certain limit on two limits in addition, if the A point is local initial point, D is ordered when walking, and function representation is t (x), x is a local one-dimension coordinate on the AC direction, and establishing the local coordinate of E on the AC direction is x _E, the distance of vertical AC direction is d, when then any point x arrives walking of node E on the AC line is

$t_E=t\left(x\right)+{[{(x-x_E)}^2+d^2]}^{\frac{1}{2}}/v---\left(3\right)$

V is the speed in the triangular unit, when making x between AC, change, and t _EGet minimumly, can obtain wavefront point is the local coordinate x of secondary source point D _D:

$\{x_D=x:t_E=\underset{x\∈(A,C)}{\min }[t\left(x\right)+{({(x-x_E)}^2+d^2)}^{\frac{1}{2}}]/v\}---\left(4\right)$

If the local coordinate that A, B, C3 are ordered is x _A, x _B, x _C, be t during 3 walk _A, t _B, t _C, when then upward any point is walked with hyperbolic curve approximate treatment AC

$t\left(x\right)=\left\{t_A^2\right[(x-x_B)(x-x_C)]/[(x_A-x_B)(x_A-x_C)]+$

$t_B^2\left[(x-x_A)(x-x_C)\right]/\left[(x_B-x_A)(x_B-x_C)\right]+---\left(5\right)$

$t_C^2\left[(x-x_A)(x-x_B)\right]/[(x_C-x_A)(x_C-x_B){\left]\right\}}^{1/2}$

E) calculate again, during relatively more minimum walk the adjacent triangular unit node of triangular unit hourage and determine the secondary source of each node, when when adjacent triangular unit node is walked, changing, this adjacent triangular unit is made activation marker,

F) calculate again, triangular unit node hourage when more current minimum is walked; When unchanged when each node of triangular unit is walked when current minimum is walked; Triangular unit was made " dormancy " sign when current minimum was walked; Get back to c) step; When changing when each node of triangular unit is walked when current minimum is walked; Triangular unit was made " activation " sign when current minimum was walked; Get back to c) step; All triangular units all are in " dormancy " state when all no longer changing when the node of the capable face of all ripples is walked and in the capable face of ripple; The capable face expansion of ripple is complete

G), when walking and secondary source information, carry out picking up the raypath from the acceptance point to the source point to the source retrieval according to each node from acceptance point.

2 triangular mesh ray tracing global method of two-dimensional complex construction as claimed in claim 1 is characterized in that: the capable face of described ripple is meant in the two dimensional surface, the plane that each node that is arrived by wave source point and a certain moment ripple constitutes.

3 triangular mesh ray tracing global method of two-dimensional complex construction as claimed in claim 1 is characterized in that: described source node is meant and excites source point.

4 triangular mesh ray tracing global method of two-dimensional complex construction as claimed in claim 1 is characterized in that: described adjacent triangular unit is meant and all relevant triangular units of a triangular unit summit or limit.

Images