CN1712991A - Ray traction in earthquake prospection - Google Patents
Ray traction in earthquake prospection Download PDFInfo
- Publication number
- CN1712991A CN1712991A CN 200410049719 CN200410049719A CN1712991A CN 1712991 A CN1712991 A CN 1712991A CN 200410049719 CN200410049719 CN 200410049719 CN 200410049719 A CN200410049719 A CN 200410049719A CN 1712991 A CN1712991 A CN 1712991A
- Authority
- CN
- China
- Prior art keywords
- ray
- point
- grid
- path
- tour
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 claims abstract description 80
- 238000009499 grossing Methods 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 238000004458 analytical method Methods 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 4
- 238000013213 extrapolation Methods 0.000 claims description 3
- 238000013459 approach Methods 0.000 claims description 2
- 238000001914 filtration Methods 0.000 claims description 2
- 238000003384 imaging method Methods 0.000 abstract 1
- 238000005266 casting Methods 0.000 description 13
- 238000010586 diagram Methods 0.000 description 11
- 238000004364 calculation method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000007664 blowing Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Images
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
A method for tracking ray in seismic prospecting picks up three points from rectangular net and solves out position of middle point based on principle of Fermat minimum hourage to raise two time calculating speed of disturbance method. The method can be widely used in seismic prospecting activities such as shift, inversion and imaging.
Description
Technical field:
The invention belongs to oil seismic exploration data processing method technical field, relate in particular to the method that a kind of minimum whilst on tour seismic ray is followed the trail of.
Background technology:
Method of seismic exploration, promptly each seismic event check point of arranging with mode according to certain rules receive by according to certain rules or mode arrange and move respectively blow out refraction wave and/or the reflection wave that a little excites the seismic event that is produced on the face of land, handle and analyze by time and waveform that the detected refraction wave that at the interface this seismic event is produced by each geology synthem and/or reflection wave are returned the face of land, so that can accurately understand and determine architectonic form and position, be a kind of the using always and effective method of in explorations such as oil and natural gas, instructing position, prospect drilling hole to dispose at present.Because underground each geology synthem is in the difference of the aspects such as degree of uniformity of composition and density, distribution, the seismic event velocity of propagation therein that a little excites generation of blowing out is corresponding also to have nothing in common with each other, and each check point is detected to be produced at the interface by each geology synthem and return also corresponding difference to some extent such as the time of the corresponding refraction wave on the face of land and/or reflection wave and waveform.
Seismic event is in the process of underground propagation, and the path that ripple passes and the speed of medium have determined the hourage of ripple.According to wave theory, ripple sets out from shot point (focal point, the earthquake centre during as earthquake), at underground propagation, through after a while, arrives acceptance point (record website) with medium velocity, the path of being walked always hourage the shortest.This shortest path (path that the ray whilst on tour is the shortest) is exactly the raypath of ripple.How determining path and whilst on tour that ripple is propagated from the shot point to the acceptance point, is exactly the main contents of ray tracing research.In seismic migration method, tomography and inversion method, extensively adopt ray-tracing scheme seismic wave whilst on tour definitely, so ray-tracing scheme is an important research field.The method that has several ray tracings in the prior art is as two point ray tracing methods (Feregra, 1980), finite difference method is separated eikonal equation method (Vidale, 1988), wavefront reconstruction estimation whilst on tour and amplitude (Vinje etc., 1993), wavefront back tracking method (Bube etc., 2001).
But two problems that all ray tracing techniques mainly face are:
1, the blind area that in tracing process, produces of ray (be raypath walk less than the zone) and be not the time shortest path.The blind area of ray mainly occurs in the low velocity band on stratum, when ray runs into low velocity layer (LVL), walks around low velocity layer (LVL), passes from high-speed layer, thereby has produced the blind area.
2, some ray-tracing scheme might can not find shortest path (minimum hourage), as cubic grid shooting method.Arithmetic speed also is the problem that must consider in the practical application.Too much algorithm consuming time although can overcome the blind zone problem of ray, also can't be used widely in actual applications.
In order to overcome the problem of above-mentioned existence, 3 disturbances of a kind of rectangular node (intermediate point is done the slip of horizontal level) ray casting has been proposed in the prior art, this method is to propose at the discrete situation of rate pattern, no matter which type of velocity field distributes or the existence of weathering zone, this method can find the path of shortest time, and having avoided the ray blind area and having followed the trail of the path is not the time shortest route problem.Its concrete operation method is:
1, trace regions is made rectangular nodeization, distribution speed value on each grid.Connect with straight line to acceptance point at focal point, obtain the initial ray path, the intersection point of initial path and grid is respectively p
1, p
2... p
NCalculating path p
1p
2P
NTotal whilst on tour.
2, get three intersection point p of initial path and rectangular node
I-1, p
iAnd p
I+1, fix two end points p
I-1And p
I+1, to intermediate point p
iRevise, provide a disturbance quantity (increment of disturbance), calculate corresponding hourage its disturbance.Might be that the position of minimum whilst on tour is all carried out such position and slided, the whilst on tour after calculating location slides, whilst on tour is relatively determined minimum whilst on tour, thereby is determined to have the position q of minimum whilst on tour then
i, use q
iReplace p
iThe position.
3, get q
i, p
I+1And p
I+2Three points are fixed two end points q
iAnd p
I+2, repeat 2 disturbance and whilst on tour calculating.Determine q
I+1Replace p
I+1The position.Get 3 corrections so successively up to being adapted to p
NTill, and calculate new route p
1q
1q
2P
NWhilst on tour T.
4, p
1q
1q
2P
NLine be used as new initial ray path, repeated for 2,3 steps.
5, repeat the 1-4 step repeatedly, obtain a series of curved rays path and corresponding whilst on tour sequence T
0＞=T
1＞=... if, certain whilst on tour T
MSatisfy T
M-1-T
M＜=ε, then tracking stops.
This technical characterstic is, can avoid ray blind area and ray tracing path is not the time shortest route problem, but there are 2 technical matterss in the method:
1, ray tracing speed is slow: in concrete operations, in making the process of disturbance, must carry out disturbance successively to all possible point, so just increase operand, greatly reduce the speed of ray tracing.
2, the selection problem of disturbance quantity: in the operation and application of reality, too small if disturbance quantity is selected, precision generally can reach requirement, but obvious increase consuming time.Otherwise excessive if disturbance quantity is chosen, therefore the requirement that does not then reach minimum whilst on tour need carry out repeated tests to disturbance quantity, to determine effective and rational span.To reach the requirement that promptly can satisfy computational accuracy on the one hand, do not increase meaningless operand on the other hand again.In the application of reality, rationally choosing of disturbance quantity is very difficult.
Summary of the invention:
3 disturbance ray castings of rectangular node of the prior art are to propose at the discrete situation of rate pattern, no matter the influence of which type of velocity field, weathering zone etc., can find the shortest time path, having avoided the ray blind area and having followed the trail of the path is not the time shortest route problem, is a comparatively desirable ray-tracing scheme.But there are two hang-ups in this method, and first speed is too slow.In making the process of disturbance, must carry out disturbance successively to all possible point.Calculate whilst on tour then, carry out the comparison of minimum whilst on tour, determine the position in final path.It two is choice problems of disturbance quantity size.Disturbance quantity is excessive, and speed can be accelerated, but can not find the shortest time path; Disturbance quantity is too small, and after computational accuracy reached, minimum whilst on tour did not have the change of any practical significance, but speed obviously slows down.Have only a suitable disturbance quantity could promptly satisfy computational accuracy, do not increase meaningless operand again.This disturbance quantity in addition of making us perplexing most need be adjusted with the difference of compute depth.When the degree of depth increased, disturbance quantity need suitably reduce.In the actual computation process, realize quite difficulty of this requirement.
The present invention is based on 3 disturbance ray castings of rectangular node, utilize the minimum whilst on tour principle of Fermat, directly obtain the position of disturbance intermediate point, do not need intermediate point is carried out disturbance successively, thereby improved 3 disturbance ray castings of rectangular node computing velocity, simultaneously because therefore the position of calculating minimum whilst on tour under the minimum whilst on tour principle of Fermat has avoided disturbance quantity to choose problem.
Concrete method is as follows:
A kind of method that is used for the seismic prospecting ray tracing, it is characterized in that: after this method carry out matrix gridization with the ray tracing field, utilize 3 of rectangular nodes and FERMAT principle, one by one the raypath point is followed the trail of, ray tracing path when forming the shortest walk;
Its concrete grammar is:
A, the ray trace regions is carried out gridding operation;
B, set up the discrete velocity field of the two-dimension speed model of Depth Domain;
C, discrete velocity field is done smoothing processing: the grid at each degree of depth interface, get the odd number point along horizontal axis, as 5,7,9 grades are carried out the smoothing processing of mathematic(al) mean or medium filtering;
The intersection point p of D, setting initial ray path and this initial path and described grid
1, p
2... p
N::
T hourage in E, calculating initial ray path
0
Raypath is followed the trail of in F, pointwise:
(1) gets a little: get three order intersection point (p on the initial ray path
1, p
2, p
3);
(2) calculate: fixing two-end-point p
1And p
3, to intermediate point p
2, calculate with the Fermat principle
Obtain the intermediate point q of minimum whilst on tour correspondence
1, and with new intermediate point q
1Replace original
Intermediate point p
2, make p
1q
1p
3Be p
1To p
3The shortest curved rays path of time;
(3) get q again
1, p
3, p
43 points are revised it, repeat the operation of (2), like this
Get 3 corrections successively up to the terminal point p that is adapted to raypath
NTill;
(4) through behind step (1)-(3), form a new raypath p
1q
1q
2P
N,
(5) calculate hourage: according to p
1q
1q
2P
NT hourage of path computing new route
1
(6) threshold ratio: with front and back (T twice hourage
0And T
1) difference and threshold value carry out
Relatively;
(7) repeating step: if threshold ratio result does not meet accuracy requirement, then with p
1q
1q
2P
N
As initial path, repeat the tracing process of (1)-(6) among the above-mentioned steps F,
T hourage makes new advances
2, and and T
1Carry out threshold ratio, do not satisfy accuracy requirement then repeatedly
The tracing process of (1) among the repeating step F-(6) constantly forms new raypath,
Raypath when pointwise approaches the shortest walk;
G, output raypath step:
(1) if (6) threshold ratio result in the step F meets accuracy requirement, i.e. travelling
The difference of time meets accuracy requirement: | T
M-1-T
M|＜=ε, tracking stops;
(2) output: corresponding raypath is the seismic ray of minimum hourage and follows the trail of
This raypath is exported in the path.
In concrete application, transverse grid at interval greater than longitudinal grid at interval in the gridding of the described steps A operation.
In concrete application, the grid interval is that transverse grid is 1: 1 with longitudinal grid ratio at interval at interval in the gridding of the described steps A operation.
In concrete operation, the discrete velocity field of the two-dimension speed model of setting up Depth Domain among the described step B is: from seismic data, obtain the stack velocity of the shear wave of the compressional wave and the corresponding degree of depth with velocity analysis method.With Dick's Si formula stack velocity is scaled interval velocity, and with time-the speed territory is transformed into the degree of depth-speed territory, forms the rate pattern of Depth Domain.Between focal point-geophone station, velocity field is carried out the discrete velocity field that interpolation and extrapolation just generate the two-dimension speed model that the Depth Domain ray tracing needs.
In concrete operation, the discrete velocity field of the two-dimension speed model of setting up Depth Domain among the described step C: make smoothing processing in the operating process, get the level and smooth of 5-9 road.
In concrete application, the FERMAT principle in the described step F, its core formula:
X wherein: grid intersection point P
I+2-P
iLevel interval;
R: grid intersection point P
I+2-P
I+1Level interval;
h
1: grid intersection point P
I+1-P
iVertical interval;
h
2: grid intersection point P
I+2-P
I+1Vertical interval;
v
1, v
2Be respectively P
I+1Point is grid speed up and down;
T: from P
i～P
I+2Whilst on tour;
Dt/dr: the time is to the differentiate of r.
In the application of reality, recursion formula is followed the trail of in the pointwise in the described step F:
x＝p
3-p
1、q
1＝p
1+(x-r)。
x＝p
i+2-q
i、q
i+1＝q
i+(x-r)，(i＝1，2，……，N-1)
The r position with minimum whilst on tour of finding the solution wherein for the Fermat principle.
P in the formula
1, p
2... p
NSeries is the horizontal coordinate of initial mesh intersection point; q
1Be to p
2The amended horizontal coordinate of point; Then p series is made amendment successively, and utilize the amended q value in front, form new q series, also just obtained p
1To p
NNew raypath.
According to above-mentioned technical scheme, 3 Fermat ray castings of the application's rectangular node are based on a kind of method that improves computing velocity of 3 methods of perturbation of rectangular node.Get three points of rectangular node, under the minimum whilst on tour principle of Fermat, ask for the position of intermediate point, and do not resemble disturbance successively the method for perturbation.Therefore, computing velocity improves more than 2 times than method of perturbation, the puzzlement that not selected by the disturbance quantity size in the perturbation process.This method has been inherited the advantage of 3 methods of perturbation of rectangular node, and the velocity field to discrete arbitrarily can find the shortest time path, and having avoided the ray blind area and having followed the trail of the path is not the time shortest route problem.In seismic migration, inverting with become to analyse formation method and be with a wide range of applications in using.
Description of drawings:
Fig. 1 for the present invention from the S point to the R point, pass the ray travelling synoptic diagram at interface;
Fig. 2 is 3 Fermat ray tracings of rectangular node of the present invention specific algorithm synoptic diagram;
Fig. 3 is the process flow diagram of 3 Fermat ray-tracing schemes of rectangular node of the present invention;
Fig. 4 is the interval velocity model of P ripple Depth Domain among the present invention;
Fig. 5 is 3 Fermat ray tracings of grid synoptic diagram of the present invention 50 * 50;
Fig. 6 is the minimum whilst on tour raypath of the PP ripple figure of 3 Fermat ray castings of rectangular node of the present invention;
Fig. 7 is the present invention and 3 on-the-spot operation interface synoptic diagram and comparison diagrams consuming time of method of perturbation ray tracing operation of rectangular node;
Fig. 8 is the superimposed displayed map in minimum whilst on tour path of the present invention and 3 method of perturbation ray-tracing schemes of rectangular node.
Concrete embodiment:
Fig. 1 for the present invention from the S point to the R point, pass the ray travelling synoptic diagram at interface.On the basis of 3 disturbance back tracking methods of rectangular node, the present invention proposes 3 Fermat ray castings of rectangular node, gets three intersection point x of rectangular node and raypath equally
I-1, x
iAnd x
I+1, to intermediate point x
i, be not to provide a disturbance quantity, but, determine to have the position (r) of minimum whilst on tour according to the minimum whilst on tour principle of Fermat to its disturbance.So once just can determine the position of minimum whilst on tour.Not the intermediate point of minimum travelling, disturbance successively.Substitute x with new intermediate point position
i, order is got next intersection point, said process, to the last an intersection point repeatedly, obtain new raypath, calculate the whilst on tour of new route, with new raypath as new initial ray path, repeat said process, less than given precision, stop following the trail of and calculate up to the change amount of whilst on tour.Faster on its computing velocity than method of perturbation, and inherited the advantage of method of perturbation, adapt to various complicated velocity fields.
Wherein, x: grid intersection point P
I+2-P
iLevel interval; R: grid intersection point P
I+2-P
I+1Level interval;
h
1: grid intersection point P
I+1-P
iVertical interval; h
2: grid intersection point P
I+2-P
I+1Vertical interval;
v
1, v
2Be respectively P
I+1Point is grid speed up and down; T: from P
i～P
I+2Whilst on tour;
Dt/dr: the time is to the differentiate of r.
As shown in Figure 1, x=x in the following formula
I+1-x
I-1Be fixing two-end-point S (x
I-1) and R (x
I+1) horizontal range; R=x
I+1-x
iBe intermediate point (x
i) to R (x
I+1) horizontal range; V1, v2 are respectively the disturbance point speed of grid up and down; H1, h2 are the disturbance point height of grid up and down.T (r) is through intermediate point x from S
i, the ray hourage of arrival R,, redefine intermediate point x if wish this section whilst on tour minimum from S to R
iThe position, 3 methods of perturbation are constantly to change x with a disturbance quantity on the interface
iCorresponding whilst on tour t (r) is calculated in the position, relatively each x
iWhilst on tour, obtain minimum whilst on tour, thereby determine x
iThe position as the raypath of minimum whilst on tour.3 Fermat methods of invention are not to x
iCarry out disturbance successively, but whilst on tour t (r) is asked derivative about r, when derivative is zero, solves and separate corresponding x about r
iIt is exactly the raypath of minimum whilst on tour.Fermat principle that Here it is.
Fig. 2 is 3 Fermat ray tracings of rectangular node of the present invention specific algorithm synoptic diagram.
Concrete grammar: among the figure, a is an initial ray, and b is hour ray.1, connect S, the straight line of R is obtained the intersection point of this straight line and mesh lines, P as the initial ray path
1, P
2..., P
N, and calculate whilst on tour T
0
2, use P
1, P
2, P
3Revise initial ray path, wherein P at 3
1, P
3Keep motionless,, directly calculate Q with the minimum whilst on tour formula of Fermat
1The position of point makes P
1Q
1P
3Be P
1To P
3The shortest curved rays path of time.
3, get Q
1, P
3, P
4Revise for 3,, get 3 corrections so successively up to being adapted to P with the 2nd step
NTill, and calculate new route P
1Q
1Q
2P
NWhilst on tour T
1
4, P
iQ
1Q
2P
NLine be used as new initial ray path, repeated for 2,3 steps.
5, repeat the 1-4 step repeatedly, obtain a series of curved rays path and corresponding whilst on tour sequence T
0＞=T
1＞=... if, certain whilst on tour T
MSatisfy T
M-1T
M＜=ε, then tracking stops.
Fig. 3 is the process flow diagram of 3 Fermat ray-tracing schemes of rectangular node of the present invention, and its concrete steps are the same.
Fig. 4 is the interval velocity model figure of P ripple Depth Domain among the present invention.
The zone of getting ray tracing is a degree of depth 0-5000 rice, and maximum focal point-geophone interval is from 4000 meters.According to the calculation process of Fig. 3,
1, with 40 * 50 meters grid gridding is carried out in the zone of ray tracing;
2, from seismic data, can obtain the stack velocity of compressional wave (P ripple) and shear wave (S ripple) with velocity analysis method.With Dick's Si formula stack velocity is scaled interval velocity, and with time-the speed territory is transformed into the degree of depth-speed territory, forms the rate pattern of Depth Domain.Between focal point-geophone station (acceptance point), velocity field is carried out the discrete velocity field (see figure 4) that interpolation and extrapolation just generate the two-dimension speed model that the Depth Domain ray tracing needs.The speed of discrete velocity field changes obviously in vertical, horizontal, accompanies high-velocity bed in the velocity layering, the low velocity layer is arranged, whole velocity field more complicated under the high-velocity bed;
3, provide the initial ray path, calculate initial ray path and grid intersection point, and calculate the whilst on tour in initial ray path;
4, circulation intersection point is got the intersection point of three orders on the ray, and fixing two-end-point calculates the intermediate point of minimum whilst on tour correspondence to intermediate point with the Fermat principle, replaces old intermediate point with new intermediate point;
5, do you judge whether intersection point arrives a last point? not, get the ray intersection downwards, form three new points, the repeating step step 3 with new point and order.Be to calculate the whilst on tour of new route;
6, do you judge that whilst on tour satisfies computational accuracy?, be not used as new initial ray path with new path, repeating step 3,4 and 5.Be to stop ray tracing.
Fig. 5 is 3 Fermat ray tracings of grid synoptic diagram of the present invention 50 * 50;
Grid with 50 * 50, mesh spacing (50 meters), velocity field is made horizontal and vertical linear change at grid.Fix 2 points (focal point-geophone station), carry out initial ray tracking and minimum whilst on tour and follow the trail of.Fig. 5 is the ray diagram of two point ray tracing of focal point-geophone station.Straight line is an initial ray, and curve is minimum whilst on tour ray.The whilst on tour minimum of curved rays, and direct rays are not to be the path of minimum whilst on tour.Because the velocity field linear change is so curved rays is very slick and sly.Wherein, a is an initial ray, and b is hour ray.
Fig. 6 is the minimum whilst on tour raypath of the PP ripple figure of 3 Fermat ray castings of rectangular node of the present invention.From the depth point 2000 meters to 4800 meters, fixing focal point-reflection spot-geophone station, 3 Fermat ray tracings of rectangle PP wave ray path profile.From Fig. 6 as seen, ray the blind area do not occur except obvious bending, and all rays can both pass through complicated rate pattern, arrives the ground acceptance point with minimum whilst on tour.
Fig. 7 is the present invention and 3 on-the-spot operation interfaces of method of perturbation ray tracing operation of rectangular node.3 method of perturbation ray tracing operations of rectangular node scene (on), the input disturbance amount is 0.05, the operation start time is 8: 58: 41, the depth point is 2500 meters, the minimum whilst on tour of PP ripple (the PP ripple is a compressional wave) is 1814 milliseconds, the minimum whilst on tour of PS ripple (the PS ripple is a transformed wave) is 2522 milliseconds, and the termination time is 8: 59: 24, calculates 43 seconds consuming time; 3 Fermat ray casting operations of rectangular node of the present invention on-the-spot (descending), the operation start time is 9: 4: 9, the depth point is 2500 meters, the minimum whilst on tour of PP ripple is 1812 milliseconds, the minimum whilst on tour of PS ripple is 2520 milliseconds, and the termination time is 9: 4: 29, calculates 20 seconds consuming time.
On the SGI workstation, with rectangular node (grid interval 40 * 50) 3 methods of perturbation, disturbance quantity selects 0.05, and (disturbance quantity was less than 0.05 o'clock, and the change amount of minimum whilst on tour is less than 4 milliseconds, but obvious increase consuming time.Disturbance quantity is excessive, as gets 0.1, does not then reach the requirement of minimum whilst on tour), therefore, disturbance quantity 0.05 is the amount of determining behind overtesting, can promptly satisfy computational accuracy, does not increase meaningless operand again.At focal point, the P ripple incides subsurface interface from three angles, and 2500 meters boundary reflection P ripple and conversion S ripple are to the ground geophone station in the depth point, the minimum whilst on tour of PP ripple is 1814 milliseconds, the minimum whilst on tour of PS ripple is 2522 milliseconds, and calculating always consuming time is 43 seconds, and computing is on-the-spot as on Fig. 7.Same computational accuracy (even than method of perturbation accurate 2 milliseconds, under Fig. 7), same depth point and focal point-geophone station, the calculating of 3 Fermat ray castings of rectangular node total consuming time be 20 seconds.About 2.15 times of raising speed.The raypath of the minimum whilst on tour of the PP ripple of two kinds of methods and PS ripple is seen Fig. 8, and solid line is the result of the compressional wave ray tracing of 3 methods of perturbation of rectangular node, and open circles is the result of 3 Fermat methods of rectangular node compressional wave ray casting; Little square frame is the result of the transformed wave ray tracing of 3 methods of perturbation of rectangular node, and ten sons are the result of 3 Fermat methods of rectangular node transformed wave ray casting.The descending incident P wave path of PP and PS ripple is identical, observes the Snall law on the interface, and reflected P ripple and conversion S ripple go upward to the geophone station position.The raypath of two kinds of methods is overlapping substantially.From the result of minimum whilst on tour as can be seen, how many difference are 2 milliseconds difference do not reflect substantially on raypath yet.But the arithmetic speed difference of two kinds of methods is very big, and at 2500 meters of the degree of depth, both computing velocitys differ from 2.15 times, and at shallow-layer, as 1000 meters, result of calculation shows that both computing velocitys differ from 2.67 times.Promptly under identical disturbance quantity condition, the interface is more shallow, and the computing velocity difference of two kinds of methods is bigger.
For 3 disturbance algorithms of rectangular node,, fix the requirement that a disturbance quantity can't satisfy minimum whilst on tour along with the increase of the degree of depth.Need suitably reduce disturbance quantity according to the increase of the degree of depth, make the consuming time further increase of this method like this.Therefore 3 Fermat ray castings of rectangular node improve more than 2 times than 3 method of perturbation arithmetic speeds, are widely used in the application of reality.
Fig. 8 is the superimposed displayed map in minimum whilst on tour path of two kinds of ray-tracing schemes.
Pp-r among the figure: the compressional wave ray of 3 methods of perturbation; Pp-f: the compressional wave ray of 3 Fermat methods; Ps-r: the transformed wave ray of 3 methods of perturbation; Ps-f: the transformed wave ray of 3 Fermat methods.
Claims (7)
1, a kind of method that is used for the seismic prospecting ray tracing is characterized in that:
Said method comprising the steps of:
A, the ray trace regions is carried out gridding operation;
B, set up the discrete velocity field of the two-dimension speed model of Depth Domain;
C, discrete velocity field is done smoothing processing:, get the smoothing processing that the odd number point carries out mathematic(al) mean or medium filtering along horizontal axis for the grid at each degree of depth interface;
The intersection point p of D, setting initial ray path and definite this initial path and described grid
1, p
2... p
N
T hourage in E, calculating initial ray path
0
Raypath is followed the trail of in F, pointwise:
(1) gets a little: get three order intersection point (p on the initial ray path
1, p
2, p
3);
(2) calculate: fixing two-end-point p
1And p
3, to intermediate point p
2, calculate the intermediate point q of minimum whilst on tour correspondence with the Fermat principle
1, and with new intermediate point q
1Replace original intermediate point p
2, make p
1q
1p
3Become p
1To p
3The curved rays path that time is the shortest;
(3) get q again
1, p
3, p
43 points are to p
3Revise, repeat the operation of (2), get 3 corrections so successively up to the terminal point p that is adapted to raypath
NTill;
(4) through behind step (1)-(3), form a new raypath p
1q
1q
2P
N
(5) calculate hourage: according to p
1q
1q
2P
NT hourage of path computing new route
1
(6) threshold ratio: with front and back (T twice hourage
0And T
1) difference and threshold value compare;
(7) repeating step: if threshold ratio result does not meet accuracy requirement, then with p
1q
1q
2P
NAs initial path, repeat the tracing process of (1)-(6) among the above-mentioned steps F, T hourage must make new advances
2, and and T
1Carry out threshold ratio, do not satisfy the then tracing process of (1)-(6) among the repeating step F repeatedly of accuracy requirement, constantly form new raypath, the raypath when pointwise approaches the shortest walk;
G, output raypath step:
(1) if (6) threshold ratio result in the step F meets accuracy requirement, promptly the difference of hourage meets accuracy requirement: | T
M-1-T
M|＜=ε, tracking stops;
(2) output: corresponding raypath is the seismic ray of minimum hourage and follows the trail of the path, exports this raypath.
2, method according to claim 1 is characterized in that: transverse grid at interval greater than longitudinal grid at interval in the gridding operation of described steps A.
3, method according to claim 1 and 2 is characterized in that: transverse grid is 1: 1 with longitudinal grid ratio at interval at interval in the gridding operation of described steps A.
4, method according to claim 1 is characterized in that: the discrete velocity field of the two-dimension speed model of setting up Depth Domain among the described step B: from seismic data, obtain the stack velocity of the shear wave of the compressional wave and the corresponding degree of depth with velocity analysis method; With Dick's Si formula stack velocity is scaled interval velocity again, and with time-the speed territory is transformed into the degree of depth-speed territory, forms the rate pattern of Depth Domain; Between focal point-geophone station, velocity field is carried out the discrete velocity field that interpolation and extrapolation just generate the two-dimension speed model that the Depth Domain ray tracing needs then.
5, method according to claim 1 is characterized in that: described step C does to get in the smoothing processing the level and smooth of 5-9 road to discrete velocity field.
6, method according to claim 1 is characterized in that: the Fermat principle in the described step F, and its application of formula is:
Wherein, x: grid intersection point P
I+2-P
iLevel interval;
R: grid intersection point P
I+2-P
I+1Level interval;
h
1: grid intersection point P
I+1-P
iVertical interval;
h
2: grid intersection point P
I+2-P
I+1Vertical interval;
v
1, v
2Be respectively P
I+1Point is grid speed up and down;
T: from P
i～P
I+2Whilst on tour;
Dt/dr: the time is to the differentiate of r.
7, method according to claim 1 is characterized in that: recursion formula is followed the trail of in the pointwise in the described step F:
x＝p
3-p
1、q
1＝p
1+(x-r)；
x＝p
i+2-q
i、q
i+1＝q
i+(x-r)，(i＝1，2，……，N-1)；
The r position with minimum whilst on tour of finding the solution wherein for the Fermat principle;
p
1, p
2... p
NSeries is the horizontal coordinate of initial mesh intersection point;
q
1Be to p
2The amended horizontal coordinate of point; Then p series is made amendment successively, and utilize the amended q value in front, form new q series, also just obtained from p
1To p
NNew raypath.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN 200410049719 CN1292263C (en) | 2004-06-25 | 2004-06-25 | Ray traction in earthquake prospection |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN 200410049719 CN1292263C (en) | 2004-06-25 | 2004-06-25 | Ray traction in earthquake prospection |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN1712991A true CN1712991A (en) | 2005-12-28 |
| CN1292263C CN1292263C (en) | 2006-12-27 |
Family
ID=35718681
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN 200410049719 Expired - Lifetime CN1292263C (en) | 2004-06-25 | 2004-06-25 | Ray traction in earthquake prospection |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN1292263C (en) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102495427A (en) * | 2011-12-10 | 2012-06-13 | 北京航空航天大学 | Interface perception ray tracing method based on implicit model expression |
| CN103105622A (en) * | 2011-11-11 | 2013-05-15 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Homomorphous wave time difference positioning method based on data base technology |
| CN102216810B (en) * | 2008-10-06 | 2013-11-06 | 雪佛龙美国公司 | System and method for deriving seismic wave fields using both ray-based and finite-element principles |
| CN106569279A (en) * | 2015-10-12 | 2017-04-19 | 中国石油化工股份有限公司 | Method for correcting three-dimensional shortest ray tracking path and device thereof |
| CN108051854A (en) * | 2018-02-01 | 2018-05-18 | 西安科技大学 | A kind of multiple dimensioned pseudo- bending ray tracing method |
| CN111324968A (en) * | 2020-03-06 | 2020-06-23 | 西南大学 | Laying method of microseismic monitoring sensors for inclined stratum tunnel engineering |
| CN112596103A (en) * | 2020-11-24 | 2021-04-02 | 中国地质科学院地球物理地球化学勘查研究所 | Ray tracing method and device and electronic equipment |
| CN114429047A (en) * | 2022-01-27 | 2022-05-03 | 成都理工大学 | Quadratic equation travel time interpolation method based on triangular mesh |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101533102B (en) * | 2009-04-09 | 2011-07-20 | 长江工程地球物理勘测武汉有限公司 | Triangular mesh ray tracing global method of two-dimensional complex construction |
| CN103076627B (en) * | 2011-10-26 | 2015-10-07 | 中国石油化工股份有限公司 | A kind of rate pattern smooth optimization method |
-
2004
- 2004-06-25 CN CN 200410049719 patent/CN1292263C/en not_active Expired - Lifetime
Cited By (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102216810B (en) * | 2008-10-06 | 2013-11-06 | 雪佛龙美国公司 | System and method for deriving seismic wave fields using both ray-based and finite-element principles |
| CN103105622A (en) * | 2011-11-11 | 2013-05-15 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Homomorphous wave time difference positioning method based on data base technology |
| CN103105622B (en) * | 2011-11-11 | 2015-08-12 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Based on the homotype ripple time difference positioning method of database technology |
| CN102495427B (en) * | 2011-12-10 | 2013-06-19 | 北京航空航天大学 | Interface perception ray tracing method based on implicit model expression |
| CN102495427A (en) * | 2011-12-10 | 2012-06-13 | 北京航空航天大学 | Interface perception ray tracing method based on implicit model expression |
| CN106569279B (en) * | 2015-10-12 | 2018-08-07 | 中国石油化工股份有限公司 | Method and apparatus for correcting three-dimensional most short ray tracing path |
| CN106569279A (en) * | 2015-10-12 | 2017-04-19 | 中国石油化工股份有限公司 | Method for correcting three-dimensional shortest ray tracking path and device thereof |
| CN108051854A (en) * | 2018-02-01 | 2018-05-18 | 西安科技大学 | A kind of multiple dimensioned pseudo- bending ray tracing method |
| CN108051854B (en) * | 2018-02-01 | 2019-09-03 | 西安科技大学 | A kind of multiple dimensioned pseudo- bending ray tracing method |
| CN111324968A (en) * | 2020-03-06 | 2020-06-23 | 西南大学 | Laying method of microseismic monitoring sensors for inclined stratum tunnel engineering |
| CN111324968B (en) * | 2020-03-06 | 2023-03-28 | 西南大学 | Laying method of microseismic monitoring sensors for inclined stratum tunnel engineering |
| CN112596103A (en) * | 2020-11-24 | 2021-04-02 | 中国地质科学院地球物理地球化学勘查研究所 | Ray tracing method and device and electronic equipment |
| CN114429047A (en) * | 2022-01-27 | 2022-05-03 | 成都理工大学 | Quadratic equation travel time interpolation method based on triangular mesh |
| CN114429047B (en) * | 2022-01-27 | 2023-08-22 | 成都理工大学 | Triangular mesh-based interpolation method for travel time of quadratic circle equation |
Also Published As
| Publication number | Publication date |
|---|---|
| CN1292263C (en) | 2006-12-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN100351650C (en) | Method for inversion constituting virtual well data using before-folded seismic wave form | |
| CN107843922A (en) | One kind is based on seismic first break and the united chromatography imaging method of Travel time | |
| CN1797037A (en) | Method for carrying out inversion for wave impedance of earthquake wave | |
| CN106353793A (en) | Cross-well seismic tomography inversion method on basis of travel time incremental bilinear interpolation ray tracing | |
| CN1748157A (en) | Methods for determining formation and borehole parameters using fresnel volume tomography | |
| CN106842295A (en) | The waveform inversion method of logging information constrained | |
| CN105093277A (en) | Shallow-medium-deep strata velocity fusion method in seismic modeling | |
| CN105093301B (en) | The generation method and device of common imaging point angle of reflection angle gathers | |
| CN1292263C (en) | Ray traction in earthquake prospection | |
| CN104749617A (en) | Multi-scale fractured reservoir forward model establishing method | |
| CN101487898A (en) | Method for oil gas water recognition by employing longitudinal wave seismic exploration post-stack data | |
| CN101561512A (en) | Multi-scale crosshole SIRT tomography method | |
| CN105301639A (en) | Speed field inversion method and device based on VSP double-weight travel time tomography | |
| CN102495427A (en) | Interface perception ray tracing method based on implicit model expression | |
| CN106199704B (en) | A kind of Three-dimendimal fusion submarine cable seismic data velocity modeling method | |
| CN107817520A (en) | The pressure coefficient Forecasting Methodology and system of marine facies mud shale stratum | |
| CN109425900A (en) | A kind of Seismic Reservoir Prediction method | |
| CN102901984B (en) | Method for constructing true earth surface dip angle trace gathers of seismic data | |
| CN103576198A (en) | Method for rapidly predicting two-dimensional offshore earthquake data free surface multiple | |
| CN102338887B (en) | Irregular-size space-variant grid tomography imaging statics correction method | |
| CN107817516A (en) | Near surface modeling method and system based on preliminary wave information | |
| CN111722283B (en) | Stratum velocity model building method | |
| CN104991272A (en) | Earthquake speed disturbance modeling method for well-free earthquake reversion | |
| CN100429530C (en) | Observation method of between well earthquake excitation and reception interchange reflection wave | |
| CN103852789B (en) | Nonlinear chromatography method and device for seismic data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| CX01 | Expiry of patent term |
Granted publication date: 20061227 |