CN115358087A

CN115358087A - Combined ray tracing method based on multi-scale and multi-direction piecewise iteration and back tracing

Info

Publication number: CN115358087A
Application number: CN202211081546.5A
Authority: CN
Inventors: 乔汉青; 张耀阳; 吕琴音
Original assignee: Institute of Geophysical and Geochemical Exploration of CAGS
Current assignee: Institute of Geophysical and Geochemical Exploration of CAGS
Priority date: 2022-09-06
Filing date: 2022-09-06
Publication date: 2022-11-18

Abstract

The invention provides a combined ray tracing method based on multi-scale and multi-direction piecewise iteration and back tracing, which comprises the following steps: s1, taking a tracking result of a conventional reverse tracking method as an initial value of a piecewise iteration method; s2, thinning the original grid to a preset sparsest level, and only reserving path points on the boundary of the remaining original grid; s3, applying conventional piecewise iteration to the sparse grid, changing the direction of the boundary of the sparse grid, and changing the retraced path into a path meeting the snell' S law by changing a space discretization mode; and S4, gradually encrypting the thinning grid, repeating the step S3 for each encryption until the current thinning level is reduced to 0, and performing the last conventional piecewise iteration. The method solves the problems that the conventional piecewise iteration method has too slow convergence, cannot track the retraced ray path and has extremely small local convergence, can adapt to the complex velocity field change, and can obtain higher precision near the seismic source.

Description

Combined ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing

Technical Field

The invention belongs to the technical field of seismic exploration, and relates to a combined ray tracing method based on multi-scale and multi-direction piecewise iteration and back tracing.

Background

The properties of seismic waves can be broadly divided into two aspects, kinematics and dynamics. In many seismic data processing and inversion methods, the requirements for seismic wave numerical simulation are different. Methods such as reverse time migration, full waveform inversion and the like require numerical simulation of the full wavefield, while methods similar to conventional migration imaging, microseismic localization, tomography and the like only need to simulate the kinematics law of the wavefield, and the calculation amount of simple kinematics numerical simulation is far less than that of the full wavefield.

Ray tracing is an important method for seismic wave kinematics numerical simulation. Currently, ray tracing mainly includes the following categories: a fire-shooting method, a bending method, a pseudo-bending method, a segment-by-segment iteration method, a back-tracking method, a shortest path method, an interpolation method, and the like. The theoretical basis of the pilot method and the piecewise iteration method is the snell's law, the theoretical basis of the back tracking method, the bending method and the pseudo-bending method is the Fermat principle, and the shortest path method and the interpolation method are the ideas of the Wheatstone principle and the interpolation principle combined on the Fermat principle.

The piecewise iteration method firstly fixes two end points of the ray and then carries out iterative update on the ray path. The piecewise iteration method does not perform overall update across multiple interfaces, so that the stability of the method is improved compared with the pseudo-bending method, but the effect is limited in a medium with violent speed change. The method can be applied to the grid model, but cannot adapt to the situation that the speed in the grid model is changed violently, namely, a certain limitation exists in a complex medium.

The back tracking method is based on a travel time field and tracks from a wave detection point to a seismic source, and has strong applicability to complex media and grid models, but the back tracking method tracks along the negative gradient direction of the travel time field, and the accuracy has great relation with the calculation of the gradient direction. The method does not have the step of carrying out iterative update on the whole ray, so the calculation efficiency is generally higher than that of a piecewise iteration method, but the corresponding calculation result precision cannot be effectively improved; on the other hand, since the tracking starts from the detection point, errors accumulate near the source.

With the development of seismic exploration to high precision, the ray tracing method is also more and more required to adapt to the situation that a high-density grid is used as a medium model. Therefore, there is a need for a ray tracing method that can both accommodate complex velocity field variations and achieve higher accuracy near the source.

Disclosure of Invention

In order to solve the defects, the invention provides a combined ray tracing method based on multi-scale and multidirectional piecewise iteration and back tracing, which solves the problems of the conventional piecewise iteration method: the method has the advantages that the convergence speed is too low, the method can be trapped into local minimum and cannot track the inflection waves, the method can adapt to complex speed field change, and higher precision can be obtained near a seismic source.

The joint ray tracing method based on multi-scale and multi-direction piecewise iteration and back tracing comprises the following steps:

s1, taking a tracking result of a conventional reverse tracking method as an initial value of a piecewise iteration method;

s2, thinning the original grid to a preset sparsest level, and only reserving path points on the boundary of the remaining original grid;

and S3, applying conventional section-by-section iteration to the sparse grid, changing the direction of the boundary of the sparse grid, and changing the mode of space discretization to change the inflection path into a path meeting the snell' S law.

And S4, gradually encrypting the thinning grid, repeating the step S3 for each encryption until the current thinning level is reduced to 0, and performing the last conventional piecewise iteration. In an embodiment of the present invention, the step S2 specifically includes the following steps:

2-1 in the original grid, giving an initial ray and a tracking point;

2-2 set the thinning level to n and the grid to the most sparse level, only preserving the path points on the boundary of the remaining grid.

In one embodiment of the present invention, in step 2-2,the thinning mode is that one line is reserved in the current grid model every more stage of thinning, and n stages of thinning are equivalent to every 2 ⁿ One for each line.

In an embodiment of the present invention, the step S3 specifically includes the following steps:

3-1, performing piecewise iteration in the thinning grid, and skipping over tracking points which do not meet the snell's law;

3-2, rotating the thinning grid and the tracking points therein clockwise around the origin;

3-3, covering the rotated sparse grid by using a horizontal-vertical grid with the same new grid interval as the sparse grid;

3-4, calculating the intersection point of the rotated ray and the new grid boundary by utilizing interpolation, namely a tracking point in the new grid;

3-5 translating the new mesh and the tracking points therein to non-negative regions;

3-6, carrying out piecewise iteration updating on the rays in the new grid, and simultaneously skipping points which do not meet the snell's law;

3-7 transforming the new mesh and corresponding tracking points back to the original mesh and corresponding tracking points.

In one embodiment of the present invention, in step 3-2, the rotation angle ranges from 0 to 90 °.

In one embodiment of the present invention, in step 3-6, the points in the new grid that do not satisfy snell's law are different from step 3-1 due to changes in the relative orientation of the grid boundaries and the rays.

In an embodiment of the present invention, the step S4 gradually encrypts the sparse grid into a new grid and is favorable for interpolation to obtain a new tracking point, and meanwhile, path points corresponding to the new grid boundary are added, and the step S3 is performed at each sparse level.

In an embodiment of the present invention, the step S4 specifically includes the following steps:

4-1, encrypting the current thinning grid once, wherein the added grid line is positioned between two adjacent grid lines, and then solving a new tracking point by utilizing interpolation;

4-2, performing step S3 on the grid of the current rarefying level;

4-3 repeating the steps 4-1 and 4-2 until the current rarefaction level is reduced to 0, and then performing the last routine section-by-section iteration.

In summary, the present invention provides a joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing, and the beneficial effects of the present invention are as follows:

aiming at the problem that the conventional piecewise iteration method is too slow in convergence, the invention provides a multi-scale piecewise iteration method, namely, rectangular grids are thinned to reduce the number of interfaces in a model, the rectangular grids are thinned to a preset sparsest level and then are gradually encrypted to an original state, and the convergence speed is effectively increased; secondly, aiming at the problem that the conventional piecewise iteration method cannot track the path of the refracted ray, a multidirectional piecewise iteration method is provided, namely the direction of the grid boundary is changed, so that the space discretization mode is changed, the refracted path is changed into a path meeting the snell's law, and the refracted wave tracking is realized; then, aiming at the problem that the local convergence of the conventional piecewise iteration method is extremely small, a combined ray tracing method based on multi-scale and multidirectional piecewise iteration and back tracing is provided, namely the result of the conventional back tracing method is used as the initial value of the improved piecewise iteration method, and the iteration process is effectively prevented from falling into local minimum. The invention can adapt to the complex velocity field change and can obtain higher precision near the seismic source.

Drawings

FIG. 1 is a schematic diagram of backward tracking single-point tracking.

Fig. 2 is a schematic diagram of a segment-by-segment iterative correction method.

FIG. 3 is a diagram of the result of a numerical simulation of an inverse tracking method; wherein, fig. 3 (a) is a forward model; FIG. 3 (b) is a diagram showing the result of numerical simulation of the inverse ray tracing method; FIG. 3 (c) is a partial enlargement of the tracking result of the seismic source region in FIG. 3 (b).

FIG. 4 illustrates the error and statistics of an inverse tracking algorithm numerical simulation; wherein, fig. 4 (a) is a forward model; FIG. 4 (b) is a schematic of the true ray (solid line) and the back trace result (dashed line); FIG. 4 (c) is a single ray global error statistic; FIG. 4 (d) shows single ray distribution error statistics.

FIG. 5 is a schematic diagram of a problem with the application of a piecewise iterative method to a rectangular grid; wherein, fig. 5 (a) is a schematic diagram of a cross-grid correction flow by a piecewise iteration method; FIG. 5 (b) is a schematic diagram of local minima in a piecewise iterative correction process; FIG. 5 (c) is a schematic illustration of a ray path of the type of a switchback wave or the like; FIG. 5 (d) is a schematic diagram of the path of the refracted wave rays in the rectangular mesh model.

FIG. 6 is a graph of the results of numerical simulations for problem (1); FIG. 6 (a) is a forward model; fig. 6 (b) is a graph of the results of the numerical simulation by the stepwise iterative method.

FIG. 7 is a graph of the numerical simulation results for problem (2); FIG. 7 (a) is a forward model; fig. 7 (b) is a graph of the results of the numerical simulation by the stepwise iterative method.

FIG. 8 is a graph of the results of numerical simulations for problem (3); wherein, fig. 8 (a) is a forward model; fig. 8 (b) is a graph of the numerical simulation result of the piecewise iteration method (only the region where the ray path is located is shown).

FIG. 9 is a schematic diagram of a multi-scale piecewise iterative method; FIG. 9 (a) shows an original grid, an initial ray path and corresponding trace points; FIG. 9 (b) is a 2-stage pump-down of FIG. 9 (a); FIG. 9 (c) is a conventional piecewise iteration of FIG. 9 (b); FIG. 9 (d) is an encryption of FIG. 9 (c); FIG. 9 (e) is a conventional piece-wise iteration of FIG. 9 (d).

FIG. 10 is a multi-scale piecewise iterative numerical simulation; wherein, fig. 10 (a) is a forward model; fig. 10 (b) is a result diagram of numerical simulation of the multi-scale piecewise iteration method.

FIG. 11 is a schematic diagram of a multi-directional piecewise iterative method; FIG. 11 (a) is a schematic view of a folded ray path; FIG. 11 (b) is a schematic diagram of the path of the folded ray after rotating clockwise by 45 degrees.

FIG. 12 is a schematic diagram of a multi-directional piecewise iterative mesh transform; wherein, fig. 12 (a) is an original mesh; FIG. 12 (b) shows the rotating grid after rotating clockwise 45 degrees; FIG. 12 (c) is the new mesh after horizontal-vertical mesh overlay at the same mesh spacing as the original mesh; fig. 12 (d) shows the new mesh after coverage being translated to a non-negative area.

FIG. 13 is a multi-directional multi-scale piecewise iterative numerical simulation result diagram; wherein, fig. 13 (a) is a forward model; fig. 13 (b) is a multi-directional, multi-scale, piecewise iterative method numerical simulation result diagram.

FIG. 14 is a workflow diagram of a joint ray tracing method based on multi-scale, multi-directional piecewise iteration and back tracing.

FIG. 15 is a diagram of numerical simulation results of a joint ray tracing method based on multi-scale, multi-directional piecewise iteration and back tracing; wherein, fig. 15 (a) is a forward model; FIG. 15 (b) is a diagram of the numerical simulation results of the joint ray tracing method based on multi-scale, multi-directional piecewise iteration and back tracing; FIG. 15 (c) is a partial enlarged view of the ray level segment.

FIG. 16 is a block diagram of the error and statistics of an exemplary numerical simulation of a joint ray tracing method based on multi-scale, multi-directional piecewise iteration and back tracing; wherein, fig. 16 (a) is a numerical simulation result diagram of the joint ray tracing method; FIG. 16 (b) is a partial enlargement of the tracing result of the seismic source region in FIG. 16 (a); FIG. 16 (c) is a single ray ensemble error statistic for the combined ray tracing method; FIG. 16 (d) is the statistics of single ray distribution errors for the combined ray tracing method.

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention more apparent, the technical solutions of the embodiments of the present invention will be described clearly and completely with reference to the accompanying drawings of the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all embodiments of the present invention. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, are within the scope of the present invention. Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, are within the scope of the present invention.

1 basic principle of conventional ray tracing method

1.1 basic principle of the Back-tracking method

The basic theory behind the back-tracking method is the Fermat principle, namely: the ray direction is the normal direction of the wavefront surface (isochronal surface) and the ray path is always perpendicular to the wavefront surface (isochronal surface). Therefore, in the travel time field, the gradient direction at a certain point is the advancing direction of the ray passing through the point.

FIG. 1 is a schematic diagram of backward tracking single-point tracking. Wherein, two squares represent two grids, the point P is the current tracking point, PA and PB represent the vertical and horizontal direction partial derivatives of the point P in the travel time field, PC represents the gradient of the point P, PD is the negative gradient direction, and D is the next tracking point.

Based on this, to find the ray path between two points, the following steps can be performed:

(1) Using the demodulator probe as the first tracking point on the path (as indicated by point P in fig. 1);

(2) Calculating the vertical and horizontal derivatives (shown as PA and PB in fig. 1) of the point in the travel time field, and the gradient (shown as PC in fig. 1);

(3) Taking the direction of the negative gradient (as indicated by PD in fig. 1) and finding the first intersection point with the grid (the intersection point with the edge or exactly the grid point) along this direction (as indicated by point D in fig. 1) as the next tracking point on the path;

(4) Repeating (2) and (3) until the origin point is tracked.

1.2 basic principle of the piecewise iteration method

The basic theoretical basis for the piecewise iteration method is snell's law. Considering that any continuous three points on the same ray path meet the snell's law, a seismic source point and a detection point can be connected by a straight line, and a plurality of intersection points are formed between the connection line and an interface passing through the path and serve as an initial tracking result; then, correcting intermediate points of all continuous three points on the ray according to the snell's law, namely enabling the position relation of all continuous three points to meet the snell's law; the step is repeated until the sum of all the correction amounts is less than a preset threshold value, and then the ray tracing is considered to be completed.

As shown in fig. 2, z = z ₂ Indicates a certain oneThe velocity of the interface, the upper and lower of the interface are v ₁ And v ₂ ，P ₁ And P ₃ Two tracking points are connected by straight line and intersect with the interface at point P ₂ ，P ₁ 、P ₂ 、P ₃ I.e. the initial tracking result. Is P' ₂ The point is a real transmission point, then P ₂ To P' ₂ The correction amount Δ x of (a) can be expressed as:

wherein a = x ₂ -x ₁ ，b＝x ₃ -x ₂ ，

2. Error analysis of conventional ray tracing methods

2.1 error analysis of the Back-tracking method

The implementation of the back tracking method comprises two main points, namely the calculation of the travel time field and the calculation of the gradient of a certain point in the travel time field. For the first point, the calculation of the travel time field is a more complex problem, and all the calculation examples in the invention use the numerical simulation of the travel time field by using a fast propulsion method. For the second point, since the velocity model and the travel time field are both discretized spatially, the differentiation in the gradient calculation needs to be replaced by a difference, and the replacement of the differentiation by the difference actually implies the assumption that the travel time field is linearly distributed locally.

According to the principle and the flow of the back tracking method, the method is analyzed to have the following error factors: firstly, except the positions of a seismic source and a detector, only a velocity model is input in a ray tracing process, but final ray tracing is completed in a travel time field, so that the error of a final ray path not only comes from the method, but also comes from the calculation result of the travel time field in the previous step, and the near-field error of the travel time field is larger than the far-field error. Secondly, the assumption that the travel time field is locally linearly distributed is not always true, and one of the assumptions is that the actual geological condition causes deviation of a tracking path, and the error is continuously accumulated in the tracking process, so that the error accumulation is larger as the position of the seismic source is closer; the second is that the wave front has less curvature in the far field and is closer to a plane wave, so that the assumption is easier to be satisfied in the far field than in the near field. Based on the above error factors, the method should yield ray tracing results with much greater error near the source than in regions far from the source.

The back tracking method was tested using numerical simulation. Using the forward model as shown in fig. 3 (a): a uniform model, with a speed set to 5000m/S, and only one seismic source point set, located at S (2500m, 3000m); 10 detectors (R1-R10) are uniformly distributed in the range of 50m-4950m on the ground surface; spatial discretization was performed using a square grid with a grid size of 800 x 500 and a grid spacing of 10m.

Fig. 3 (b) shows a comparison of the tracking results for the inverse tracking (dashed line) and the analytical solution (solid line) at the entire model scale, with no apparent macroscopic error between the two. Amplifying the result in the area near the seismic source, as shown in fig. 3 (c), it can be observed that there is a significant error between the back-tracking result and the analytic solution, and there are two characteristics: firstly, the rays are zigzag bent, and secondly, paths of a plurality of rays are combined.

The above calculation example visually shows the overall effect of the back-tracking method and the characteristics of errors near the seismic source, and in order to further quantify the distribution of the errors, 500 detectors on the ground surface are added on the basis of the forward modeling shown in fig. 3 (a), and the settings of other parameters are unchanged, as shown in fig. 4 (a). In the homogeneous model, the real ray (solid line) from the source to a detection point and the back tracking result (dotted line) are shown in fig. 4 (b), when the overall error of the tracking result is quantified, the average distance between the two lines is obtained and is called the single-ray overall error, and when the distribution of the error on the ray path is explored, the distribution of the distance between the two lines is used as the quantification of the error and is called the single-ray distribution error. Based on the two error quantification methods, error analysis is carried out on 500 rays to obtain corresponding statistics of the single-ray integral error and the single-ray distribution error. As shown in fig. 4 (c) (d).

The single-ray overall error of the back tracking method is shown in fig. 4 (c), the error distribution is approximately symmetrical about the seismic source position, and the error change is smooth; from the relative size, when the included angle between the ray and the vertical direction is in the directions of 0 degree and +/-30 degrees, the error is obviously local minimum, wherein the error about +/-30 degrees is global minimum; from the absolute size, the error is up to 15m, which is equivalent to that the whole ray deviates from the real position by about 1.5 grid intervals, and the average value of the errors of all the detection points is about 9.2m, which is already close to 1 grid interval.

The error of the single ray distribution in the back tracking method is shown in fig. 4 (d), and the right part of the graph has a blank area due to the difference of the total lengths of the rays. Viewed from the longitudinal axis, the law is substantially the same as that shown in fig. 4 (c); when viewed from the direction of the transverse axis, the error distribution characteristic of the optical fiber is in accordance with the characteristic that the near-field error is far larger than the far-field error; the error of the local area of the partial ray is already about 25m in absolute size, i.e. 2.5 grid intervals.

2.2 error analysis by piecewise iteration

According to the principle and the flow of the piecewise iteration method, the method is analyzed to have the following problems:

(1) Since each correction only involves a single interface and two velocity values on both sides thereof, and different grids have different velocity values, the correction cannot be performed across the grids. If cross-grid correction is to be implemented, a tracking point close to a grid point needs to be returned to the grid point, and then a tracking point coincident with the grid point needs to be corrected. The process is shown in fig. 5 (a), wherein the virtual broken line, the real broken line and the dotted broken line represent three ray paths, and the virtual broken line path needs to be manually returned to the real broken line path and then corrected to obtain the dotted broken line path. Even so, a correction can only span one grid at most, which may result in too slow a convergence rate of the method.

(2) The mesh involved in the correction process is always limited to the ray path and the part of the periphery thereof within one to two mesh intervals, which results in that the method cannot obtain the speed information of other areas, and therefore the result of the method only represents one local minimum travel time and not a global minimum travel time, which violates the Fermat principle. For example, in fig. 5 (b), if the distance is far greater than the maximum distance, the ray path from the point to the point should be refracted as shown by a dashed line, and in the piecewise iteration method, according to the flow, two points need to be connected by a straight line (a dot-dash straight line), and any adjacent three point on the straight line satisfies the snell's law, so that the correction amount is directly zero, the iteration cannot be performed, and the correct ray path cannot be obtained.

(3) There may be a back-fold in the first-arrival wave, as shown in fig. 5 (c), i.e., the ray propagates downward (or upward/left/right) and then upward (or downward/right/left). At this time, at the inflection point of the ray path from the lower to the upper, there is necessarily a case similar to that shown in fig. 5 (d), that is, for the consecutive three path points, the ordinate of the middle point is smaller than the minimum value of the ordinates of the two end points. Similar situations can occur for other direction of inflection. However, in snell's law, the abscissa and ordinate of the intermediate point must be between the abscissas and ordinates of the two end points when the interface is horizontal or vertical. That is, this case is not satisfied by snell's law, so the piecewise iterative method cannot track the rays of the back-folded wave or similar wavefield in a rectangular mesh model.

And performing numerical simulation aiming at the three problems of the piecewise iteration method, and performing comparative analysis by taking the back tracking result as an approximate analytic solution. For the first problem, a forward model was constructed as shown in fig. 6 (a), which was divided into two layers, the depth of the bottom interface of the shallow layer was 2000m, the speed was 2000m/s, 1 demodulator probe R (4950m, 0m) was set, and the remaining parameters were unchanged. In order to adjust the calculated amount of the two methods to the same degree to ensure that the experimental variables are the same, the number of iterations of the piecewise iteration method is set to 120. The result of the numerical simulation by the piecewise iteration method is shown in fig. 6 (b), and the result of 120 iterations by the piecewise iteration method is almost not different from the initial value, so that the tracking is considered to be failed.

For the second problem, a forward model was constructed as shown in fig. 7 (a), which set the positions of the source and detector at S (2000m, 1500 m) and R (6000m, 1500 m), respectively, on the basis of the model shown in fig. 6 (a), and the number of iterations was kept constant for 120 times. The result of the numerical simulation of the piecewise iteration method is shown in fig. 7 (b), the previous first arrival wave of the source and the detector should be a refracted wave, and the tracking result of the piecewise iteration method remains unchanged from the initial value, so that the tracking is considered to be failed.

For the third problem, a forward model is constructed as shown in fig. 8 (a), which is similar to the model shown in fig. 6 (a), and the two-layer velocity model is changed into a spatial uniform velocity variation model, the velocity increases linearly in the direction of the arrow, the velocity values at (0 m,0 m) are 2000m/S, the velocity values at (5000 m ) are 5000m/S, the positions of the seismic source and the detector are S (2000 m, 500m) and R (6000 m, 500m), respectively, and the number of iterations is kept unchanged for 120 times. The result of the numerical simulation of the piecewise iteration method is shown in fig. 8 (b), where the previous first arrival wave of the source and detector should be a backward wave, and the result of the piecewise iteration method still remains the initial value, so that the tracking is considered to be failed.

By synthesizing three problems existing in the piecewise iteration method and the corresponding numerical simulation analysis results, it can be known that in the rectangular grid model, the conventional piecewise iteration method cannot give correct results within reasonable time consumption, and even cannot effectively iterate and update the initial values under certain conditions.

3 Joint ray tracing method based on multi-scale and multi-direction piecewise iteration and back tracing

In a rectangular grid model, a conventional piecewise iteration method has three problems of too slow convergence, local minimum and incapability of tracking a back folding wave, and the main reasons for the three problems are as follows: too many interfaces in the model (all grid boundaries are considered as interfaces); the difference between the initial ray and the real ray is large; the interfaces are both horizontal and vertical directions and the theoretical limits of snell's law, for which reason the present invention will improve upon conventional piecewise iterative methods.

3.1 Multi-Scale piecewise iterative method

In order to solve the problem, the rectangular grid is thinned to reduce the number of interfaces in a model and further improve the calculation efficiency of the piecewise iteration method, and a multi-scale piecewise iteration method is provided; then gradually encrypting the grids, adding path points on the corresponding grid boundary, and performing conventional piecewise iteration on each sparse level; and finally, repeating the previous step until the grid is encrypted to the original state. The specific process is as follows:

(1) In the original grid, an initial ray and a tracking point are given. Taking fig. 9 (a) as an example, in a square grid model with a grid point number of 9*9, P is required to be selected ₁ And P ₁₆ The ray path between, first, P ₁ And P ₁₆ The two points are connected by a straight line and interpolated to obtain a tracking point P on the path ₂ To P ₁₅ 。

(2) And setting the thinning level as n, thinning the grid to the most sparse level, and only reserving path points on the boundary of the residual grid. The method adopts a thinning mode that one partition line of the current grid model is reserved every more stage of thinning, and n stages of thinning are equivalent to every 2 ⁿ One for each line. For example, fig. 9 (a) is subjected to 2-stage thinning, and the result of thinning is shown in fig. 9 (b). The lattice of 9*9 is changed into 3*3, P on the path ₂ To P ₁₅ These 14 trace points leave only P ₈ And P ₉ Two, renumbering it as P ₂ And P ₃ ，P ₄ Is the previous endpoint P ₁₆ 。

(3) A conventional piece-wise iteration method is applied on the current grid. For example, a conventional piecewise iteration is performed on FIG. 9 (b), with the result shown in FIG. 9 (c).

(4) Encrypting the current grid for one time, wherein the added grid line is positioned between two adjacent grid lines, namely the thinning level is changed into n-1 level; then, a new tracking point is obtained by interpolation. For example, as a result of encrypting once in fig. 9 (c), the lattice becomes 5*5 as shown in fig. 9 (d). P before encryption ₂ 、P ₃ 、P ₄ P after three-point corresponding encryption ₄ 、P ₅ And P ₈ P after encryption ₂ 、P ₃ 、P ₆ 、P ₇ Four points are new tracking points.

(5) And performing a conventional section-by-section iteration method on the grid of the current thinning level. For example, a conventional piecewise iteration is performed on FIG. 9 (d), with the results shown in FIG. 9 (e).

(6) And (5) repeating the steps (4) and (5) until the current rarefaction level is reduced to 0, namely the current rarefaction level is the same as the original grid model, and performing the last conventional piecewise iteration.

The model shown in fig. 6 (also fig. 10 (a)) is numerically simulated by using a multi-scale piecewise iteration method, and a larger maximum thinning level is given considering that the initial ray is greatly different from the real ray, and is set to be 6 here, while the number of iterations of each level is 20, and the total number of iterations is kept unchanged at 120. The tracking effect is shown in fig. 10 (b), and compared with the conventional piecewise iteration method, the multi-scale piecewise iteration method has the advantages that the tracking effect is obviously improved, the tracking effect is greatly overlapped with the back tracking result, and the precision of the multi-scale piecewise iteration method is considered to be improved to a degree close to that of the back tracking method.

3.2 multidirectional piecewise iteration method

Due to the existence of the interface of the rectangular mesh model and the theoretical limitation of snell's law, the piecewise iteration method cannot track the ray path with the inflection, as shown in fig. 11 (a). For the problem, the invention considers the characteristics of the interface in the rectangular grid model, namely: the grid is only the result of the spatial discretization of the medium model, the grid boundary does not represent the real boundary in the actual medium, the real boundary information is hidden in a large number of grids, and the information is not damaged by changing the spatial discretization mode. According to the above feature, the direction of the mesh boundary is changed to change the spatial discretization, so that the inflection path is changed to a path satisfying snell's law, as shown in fig. 11 (b). For ease of viewing, fig. 11 (b) shows the ray rotated, rather than the grid boundary, which is equivalent to rotating the grid. After rotation processing, the points which originally do not satisfy snell's law can be updated by a piecewise iteration method.

Based on the above thought, a multi-direction piecewise iteration method is proposed, wherein the direction refers to the relative direction of the grid and the ray. The method comprises the following specific processes:

(1) The piecewise iteration is performed in the original grid (as in fig. 12 (a)) and tracking points that do not satisfy snell's law are skipped.

(2) The original grid and the tracking points therein are rotated clockwise by any angle in the range of 0-90 degrees around the origin (in this case, 45 degrees is taken as an example), as shown in fig. 12 (b).

(3) The rotated original mesh is overlaid with a horizontal-vertical mesh having the same mesh spacing as the original mesh, as shown in fig. 12 (c).

(4) And calculating the intersection point of the rotated ray and the boundary of the new grid, namely the tracking point in the new grid by utilizing interpolation.

(5) The new grid and the tracking points therein are translated to non-negative regions as shown in fig. 12 (d).

(6) And performing piecewise iteration updating on the rays in the new grid, and simultaneously skipping points which do not meet the snell's law. Wherein, because the relative direction of the grid boundary and the ray changes, the point in the new grid which does not satisfy snell's law is different from that in step (1).

(7) And transforming the new grid and the corresponding tracking points back to the original grid and the corresponding tracking points.

The above steps are a complete iteration of the multidirectional piecewise iteration method, wherein two iterations corresponding to the original grid and the rotating grid are actually included.

The model shown in fig. 8 (also fig. 13 (a)) was numerically simulated by a multidirectional piecewise iterative method. The original grid and the 45-degree grid are selected as parameters of a multi-direction method, in order to accelerate convergence speed, a multi-scale method is added in the test, the parameters are set to be 6-level rarefying and 20 times of each level of iteration (the original grid and the 45-degree grid are respectively 20 times), and the other parameters are unchanged. The tracking result is shown in fig. 13 (b), in which only the region where the ray path is located is shown. To ensure that the tracking results do come from a multi-directional piecewise iterative method, fig. 13 (b) also gives the results using only the multi-scale method for comparison. The result of the multidirectional and multi-scale piecewise iteration method has high coincidence with the result of the back tracking, and the effectiveness of the method is proved; while simply using the multi-scale piecewise iteration method is completely ineffective because the original ray is already perpendicular to all the grid boundaries that are traversed. 3, joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing

The multi-scale and multi-direction piecewise iteration method can obtain the precision close to that of a back tracking method under the condition that local minimum does not exist, in order to improve the initial ray precision of the improved piecewise iteration method, solve the problem of local minimum and further improve the convergence speed, the invention takes the result of the back tracking method as the initial value of the multi-scale and multi-direction piecewise iteration method, and provides a combined ray tracking method based on multi-scale and multi-direction piecewise iteration and back tracking.

As shown in fig. 14, a specific flow of the joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing includes the following steps:

(1) Firstly, obtaining a tracking result of a conventional back tracking method;

(2-1) taking the tracking result of the conventional back tracking method as an initial ray and a tracking point of a piecewise iteration method;

(2-2) setting the rarefying level to be n and rarefying the original grid to the sparsest level;

(2-3) setting the iteration number to be m, and if m is more than 0, performing the step (2-4); otherwise, go to step (4-1)

(2-4) the number of iterations m becomes m-1; carry out (3-1)

(3-1) applying a conventional piecewise iteration method to the thinning grid and skipping points which do not meet the snell's law;

(3-2) rotating the thinning grid and the tracking points therein clockwise around the origin by any angle within the range of 0-90 degrees;

(3-3) covering the rotated thinning grid with a new grid with the same grid interval as the thinning grid;

(3-4) calculating a tracking point (intersection point of the rotated ray and the new grid boundary) in the new grid by using interpolation;

(3-5) translating the new grid and the tracking points in the new grid to a non-negative area for carrying out section-by-section iterative updating and skipping points which do not satisfy the snell's law;

(3-6) converting the new grid and the corresponding tracking points back to the original thinning grid and the corresponding tracking points;

(3-7) setting n =0, if yes, performing step 2-3; otherwise, performing (3-8);

(3-8) setting n = n-1, carrying out primary encryption on the current thinning grid (the added grid line is positioned between two adjacent grid lines) and solving a new tracking point by utilizing interpolation; continuing to carry out the step (2-3);

(4-1) setting n =0, and if not, performing step (4-2); if yes, performing the step (4-5);

(4-2) setting the thinning grade n as n-1; continue (4-3)

(4-3) carrying out primary encryption on the current thinning grid (the added grid line is positioned between the two adjacent grid lines) and solving a new tracking point by utilizing interpolation;

(4-4) carrying out a conventional piecewise iteration method on the grid of the current thinning level; continuing to step (4-1);

(4-5) carrying out conventional piecewise iteration,

(4-6) outputting a tracking result;

and (4-7) finishing.

The model shown in fig. 7 (also fig. 15 (a)) is numerically simulated by using a joint ray tracing method based on multi-scale, multi-directional piecewise iteration and back tracing, and the results are shown in fig. 15 (b, c). The result shows that the combined ray tracing method avoids tracing the refracted wave path due to local minimum, and the backward tracing method has a plurality of small-amplitude bends in the ray horizontal segment, so that the combined ray tracing method is obviously improved. To ensure that the reason for avoiding local minima comes from the joint ray tracing method, the result of using the multi-directional, multi-scale piecewise iterative method, which is still not different from the initial value, is also given in fig. 15 (b), so the method is considered to fail in tracing. In order to carry out error analysis and statistics on the result of the combined ray tracing method based on multi-scale and multi-direction piecewise iteration and back tracing, the invention initializes rays and tracing points by the numerical simulation result of the back tracing method, sets parameters to be 45-degree rotation, 5-level thinning and 10 iterations of each level, and carries out numerical simulation on the models shown in the figures 3 and 4.

The tracking result of the combined ray tracing method is shown in fig. 16 (a), and compared with the numerical simulation result of the back tracing method, the result of the combined ray tracing method is good in fitting with an analytic solution even if the seismic source region is locally amplified as shown in fig. 16 (b), so that the tracking error is effectively reduced, and the precision is obviously improved. The statistics of the single-ray overall error of the combined ray tracing method are shown in fig. 16 (c), and compared with the back tracing method shown in fig. 4 (c), the absolute size of the single-ray overall error of 500 rays is reduced to less than 5m, the average value of the errors is reduced to 0.86m, and the accuracy of the single-ray overall tracing result is improved by 90.65%. The statistics of the single-ray distribution errors in the combined ray tracing method are shown in fig. 16 (d), in which most of the region errors are close to 0, few region errors have absolute values exceeding 5m, and there are substantially no regions having absolute differences of 10m or more. Through the analysis of the error results, the combined ray tracing method based on multi-scale and multidirectional piecewise iteration and back tracing is verified to improve the precision of the conventional back tracing method by about one order of magnitude.

In summary, the present invention is to solve three problems existing when the conventional piecewise iteration method is applied in the rectangular grid: the invention provides three methods or strategies with pertinence, wherein the convergence speed is too slow, the convergence speed possibly falls into local minimum and the inflection wave cannot be tracked: (1) A multi-scale method is provided, a rough ray path is constructed by using a sparse model, and then the model is encrypted to depict details, so that the convergence speed of a piecewise iteration method is effectively improved; (2) A multidirectional method is provided, and rays such as inflection waves which do not meet the snell's law in the original grid are adjusted to a state which can be updated by a piecewise iteration method by changing the direction of the grid boundary; (3) The result of the back tracking method is used as the initial value of the piecewise iteration method, and the iteration process is prevented from falling into local minimum. From the perspective of the back tracking method, it can also be considered that the present invention uses the piecewise iteration method as a post-processing procedure to further correct the result of the back tracking method. Therefore, the invention effectively overcomes various defects in the prior art and has high industrial utilization value.

The present invention is not limited to the above-described preferred embodiments, but various modifications and changes can be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. The joint ray tracing method based on multi-scale and multidirectional piecewise iteration and back tracing is characterized by comprising the following steps of:

And S4, gradually encrypting the thinning grid, repeating the step S3 for each encryption until the current thinning level is reduced to 0, and performing the last conventional piecewise iteration.

2. The joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing as claimed in claim 1, wherein said step S2 specifically comprises the following procedures:

2-1 in the original grid, giving an initial ray and a tracking point;

3. The multi-scale, multi-directional based on claim 2The joint ray tracing method of the piecewise iteration and the back tracing is characterized in that in the step 2-2, a thinning mode is that one line is reserved at the interval of the current grid model for every more stage of thinning, and n stages of thinning are equivalent to every 2 ⁿ One for each line.

4. The method for joint ray tracing based on multi-scale, multi-direction piecewise iteration and back tracing as claimed in claim 1, wherein said step S3 specifically comprises the following steps:

3-1, performing piecewise iteration in the thinning grid, and skipping over tracking points which do not satisfy the snell's law;

3-7 transform the new mesh and corresponding tracking points back to the original mesh and corresponding tracking points.

5. A joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing as claimed in claim 4, characterized in that in step 3-2, the angle of rotation is in the range of 0-90 °.

6. A joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing as claimed in claim 4, characterized in that in step 3-6, the points in the new grid that do not satisfy Snell's law are different from step 3-1 due to the change in the grid boundaries and the relative direction of the rays.

7. The joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing as claimed in claim 1, wherein said step S4 gradually encrypts the sparse grid into a new grid and facilitates interpolation to obtain new tracing points, and meanwhile, path points on the boundary of the corresponding new grid are added, and step S3 is performed at each sparse level.

8. The joint ray tracing method based on multi-scale, multi-direction piecewise iteration and back tracing as claimed in claim 7, wherein said step S4 specifically comprises the following procedures:

4-2, performing step S3 on the grid of the current rarefaction level;