Joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking

Abstract

The present invention provides a joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking, which comprises the following steps: S1, taking the tracing result of the conventional backtracking method as the initial value of the segment-by-segment iteration method; S2, thinning the original grid to a preset sparsest level, and only retaining the path points on the boundary of the remaining original grid; S3, applying conventional segment-by-segment iteration to the thinned grid and changing the direction of the boundary of the thinned grid, and by changing the way of spatial discretization, making the refracted path become a path that satisfies Snell's law; S4, gradually densifying the thinned grid, repeating step S3 each time of densification until the current thinned level drops to 0, and then performing the last pass of conventional segment-by-segment iteration. The present invention solves the problems of slow convergence, inability to trace refracted ray paths and local convergence minimum existing in the conventional segment-by-segment iteration method. The present invention can not only adapt to complex velocity field changes, but also obtain high accuracy near the seismic source.

Description

Technical Field

[0001] The invention belongs to the technical field of seismic exploration and relates to a joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and reverse tracing. Background Art

[0002] The properties of seismic waves can be roughly divided into two aspects: kinematics and dynamics. Among the numerous seismic data processing and inversion methods, the requirements for numerical simulation of seismic waves are also different. Methods such as reverse time migration and full waveform inversion require numerical simulation of the full wave field, while methods such as conventional migration imaging, microseismic positioning, and tomography only need to simulate the kinematic laws of the wave field. The computational complexity of pure kinematic numerical simulation is much smaller than that of the full wave field.

[0003] Ray tracing is an important method for numerical simulation of seismic wave kinematics. At present, ray tracing mainly includes the following categories: test shooting method, bending method, pseudo bending method, segment-by-segment iteration method, reverse tracing method, shortest path method, interpolation method, etc. Among them, the theoretical basis of the test shooting method and segment-by-segment iteration method is Snell's law, the theoretical basis of the reverse tracing method, bending method and pseudo bending method is Fermat's principle, and the shortest path method and interpolation law are the combination of Fermat's principle and Huygens' principle.

[0004] The segment-by-segment iteration method first fixes the two end points of the ray, and then iteratively updates the ray path. The segment-by-segment iteration method does not cross multiple interfaces for overall update, so its stability is improved compared to the pseudo-bending method, but its effect is limited in media with drastic velocity changes. This type of method can be applied to grid models, but it cannot adapt to situations where there are drastic velocity changes in the grid model, that is, it has certain limitations in complex media.

[0005] The reverse tracking method is based on the travel time field and tracks from the detection point to the source. It has strong applicability to complex media and grid models, but the reverse tracking method tracks along the negative gradient direction of the travel time field, and the accuracy is closely related to the calculation of the gradient direction. This type of method does not have the step of iteratively updating the entire ray, so the calculation efficiency is usually higher than the segment-by-segment iteration method, but the accuracy of the corresponding calculation results cannot be effectively improved; on the other hand, since the tracking starts from the detection point, the error will accumulate near the source.

[0006] As seismic exploration develops towards high precision, ray tracing methods are increasingly required to adapt to the situation where high-density grids are used as medium models. Therefore, a ray tracing method is needed that can adapt to complex velocity field changes and obtain high accuracy near the earthquake source. Summary of the invention

[0007] To solve the above deficiencies, the present invention provides a combined ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking, which solves the problems existing in the conventional segment-by-segment iteration method: too slow convergence speed, possible entrapment in local minima, and inability to trace back-reflected waves. The present invention can not only adapt to complex velocity field changes but also obtain high accuracy near the seismic source.

[0008] The combined ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking includes the following steps:

[0009] S1. Use the tracing result of the conventional backtracking method as the initial value of the segment-by-segment iteration method;

[0010] S2. Thin the original grid to the preset most sparse level, and only retain the path points on the boundary of the remaining original grid;

[0011] S3. Apply the conventional segment-by-segment iteration to the thinned grid and change the direction of the boundary of the thinned grid. By changing the way of spatial discretization, the back-reflected path is changed into a path that satisfies Snell's law.

[0012] S4. Gradually densify the thinned grid. Repeat step S3 each time it is densified until the current thinned level drops to 0, and then perform the last pass of the conventional segment-by-segment iteration. In an embodiment of the present invention, the step S2 specifically includes the following processes:

[0013] 2-1. In the original grid, given the initial ray and the tracing point;

[0014] 2-2. Set the thinning level as n and thin the grid to the most sparse level, and only retain the path points on the boundary of the remaining grid.

[0015] In an embodiment of the present invention, in step 2-2, the thinning method is to retain one line every other line in the current grid model for each additional level of thinning. n-level thinning is equivalent to retaining one line every 2 n lines.

[0016] In an embodiment of the present invention, the step S3 specifically includes the following processes:

[0017] 3-1. Perform segment-by-segment iteration in the thinned grid and skip the tracing points that do not satisfy Snell's law;

[0018] 3-2. Rotate the thinned grid and the tracing points therein clockwise around the origin;

[0019] 3-3. Cover the rotated thinned grid with a horizontal-vertical grid having the same new grid spacing as the thinned grid;

[0020] 3-4. Use interpolation to calculate the intersection points of the rotated ray and the new grid boundary, that is, the tracing points in the new grid;

[0021] 3-5 Translate the new grid and the tracking points therein to the non-negative region;

[0022] 3-6 Iteratively update the rays in the new grid segment by segment, skipping the points that do not satisfy Snell's law;

[0023] 3-7 Transform the new grid and the corresponding tracking points back to the original grid and the corresponding tracking points.

[0024] In an embodiment of the present invention, in step 3-2, the rotation angle range is 0-90°.

[0025] In an embodiment of the present invention, in step 3-6, due to the change in the relative direction between the grid boundary and the ray, the points in the new grid that do not satisfy Snell's law are different from those in step 3-1.

[0026] In an embodiment of the present invention, step S4 gradually densifies the thinned grid into a new grid and facilitates interpolation to obtain new tracking points, while adding path points on the corresponding new grid boundary, and step S3 is performed at each thinning level.

[0027] In an embodiment of the present invention, step S4 specifically includes the following processes:

[0028] 4-1 Perform one densification on the current thinned grid, add grid lines in the middle of two adjacent grid lines, and then use interpolation to obtain new tracking points;

[0029] 4-2 Perform step S3 on the grid at the current thinning level;

[0030] 4-3 Repeat steps 4-1 and 4-2 until the current thinning level drops to 0, and then perform the last conventional segment-by-segment iteration.

[0031] In summary, the present invention provides a joint ray tracing method based on multi-scale, multi-direction segment-by-segment iteration and backtracking. The beneficial effects of the present invention are:

[0032] In view of the slow convergence of the conventional piecewise iteration method, the present invention proposes a multi-scale piecewise iteration method, that is, thinning the rectangular grid to reduce the number of interfaces in the model, thinning to the preset sparsest level and then gradually densifying the grid to the original state, effectively improving the convergence speed; secondly, in view of the problem that the conventional piecewise iteration method cannot track the existence of refracted ray paths, a multi-directional piecewise iteration method is proposed, that is, by changing the direction of the grid boundary, thereby changing the way of spatial discretization, so that the refracted path becomes a path that satisfies Snell's law, realizing the tracking of refracted waves; then, in view of the problem of local convergence minimum of the conventional piecewise iteration method, a joint ray tracing method based on multi-scale, multi-directional piecewise iteration and backtracking is proposed, that is, using the result of the conventional backtracking method as the initial value of the improved piecewise iteration method, effectively avoiding the iteration process from falling into local minimum. The present invention can not only adapt to complex velocity field changes, but also obtain high accuracy near the seismic source. Description of the Drawings

[0033] Figure 1 It is a schematic diagram of single-point tracking for backtracking.

[0034] Figure 2 It is a schematic diagram of the correction of the piecewise iteration method.

[0035] Figure 3 It is a numerical simulation result diagram of the backtracking method; among them, Figure 3 (a) is the forward model; Figure 3 (b) is the numerical simulation result diagram of the backtracking ray tracing method; Figure 3 (c) is the local magnification of the tracking result in the seismic source area in Figure 3 (b).

[0036] Figure 4 It is the error and statistics of the numerical simulation of the backtracking method example; among them, Figure 4 (a) is the forward model; Figure 4 (b) is a schematic diagram of the true ray (solid line) and the backtracking result (dashed line); Figure 4 (c) is the overall error statistics of a single ray; Figure 4 (d) is the distribution error statistics of a single ray.

[0037] Figure 5 It is a schematic diagram of the problems existing in the application of the piecewise iteration method in a rectangular grid; among them, Figure 5 (a) is a schematic diagram of the cross-grid correction process of the piecewise iteration method; Figure 5 (b) is a schematic diagram of local minimum during the correction process of the piecewise iteration method; Figure 5 (c) is a schematic diagram of ray paths such as refracted waves; Figure 5 (d) is a schematic diagram of the refracted wave ray path in the rectangular grid model.

[0038] Figure 6 It is the numerical simulation result diagram for problem (1); Figure 6 (a) is the forward model; Figure 6 (b) is the numerical simulation result diagram of the piecewise iteration method.

[0039] Figure 7 It is the numerical simulation result diagram for problem (2); Figure 7 (a) is the forward model; Figure 7 (b) is the numerical simulation result diagram of the piecewise iteration method.

[0040] Figure 8 It is the numerical simulation result diagram for problem (3); among them, Figure 8 (a) is the forward model; Figure 8 (b) is the numerical simulation result diagram of the piecewise iteration method (only the area where the ray path is located is shown).

[0041] Figure 9 It is the schematic diagram of the multi-scale piecewise iteration method; among them, Figure 9 (a) is the original grid, the initial ray path and the corresponding tracking points; Figure 9 (b) is for Figure 9 (a) is thinned at level 2; Figure 9 (c) is for Figure 9 (b) is subjected to conventional piecewise iteration; Figure 9 (d) is for Figure 9 (c) is encrypted; Figure 9 (e) is for Figure 9 (d) is subjected to conventional piecewise iteration.

[0042] Figure 10 It is the numerical simulation of the multi-scale piecewise iteration method; among them, Figure 10 (a) is the forward model; Figure 10 (b) is the numerical simulation result diagram of the multi-scale piecewise iteration method.

[0043] Figure 11 It is the schematic diagram of the multi-directional piecewise iteration method; among them, Figure 11 (a) is the schematic diagram of the back-refracted ray path; Figure 11 (b) is the schematic diagram of the path after the back-refracted ray is rotated 45 degrees clockwise.

[0044] Figure 12 It is the schematic diagram of the multi-directional piecewise iteration grid transformation; among them, Figure 12 (a) is the original grid; Figure 12 (b) is the rotated grid after being rotated 45 degrees clockwise; Figure 12 (c) is the new grid covered by the horizontal-vertical grid with the same grid spacing as the original grid; Figure 12 (d) is to translate the covered new grid to the non-negative region.

[0045] Figure 13 It is a graph of the numerical simulation results of multi-directional and multi-scale iterative segment by segment; among them, Figure 13 (a) is the forward model; Figure 13 (b) is the graph of the numerical simulation results of the multi-directional and multi-scale iterative method segment by segment.

[0046] Figure 14 It is a flowchart of the joint ray tracing method based on multi-scale, multi-directional iterative segment by segment and backtracking.

[0047] Figure 15 It is a graph of the numerical simulation results of the joint ray tracing method based on multi-scale, multi-directional iterative segment by segment and backtracking; among them, Figure 15 (a) is the forward model; Figure 15 (b) is the graph of the numerical simulation results of the joint ray tracing method based on multi-scale, multi-directional iterative segment by segment and backtracking; Figure 15 (c) is the locally enlarged graph of the horizontal segment of the ray.

[0048] Figure 16 It is the error and statistics of the numerical simulation of the example of the joint ray tracing method based on multi-scale, multi-directional iterative segment by segment and backtracking; among them, Figure 16 (a) is the graph of the numerical simulation results of the joint ray tracing method; Figure 16 (b) is for Figure 16 the local enlargement of the tracking results in the source region in (a); Figure 16 (c) is the overall error statistics of a single ray of the joint ray tracing method; Figure 16 (d) is the distribution error statistics of a single ray of the joint ray tracing method. Specific implementation manners

[0049] To make the purposes, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed present invention, but merely represents the selected embodiments of the present invention. Based on the embodiments in the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of protection of the present invention.

[0050] 1 Basic principle of the conventional ray tracing method

[0051] 1.1 Basic principle of the backtracking method

[0052] The basic theoretical basis of the backtracking method is Fermat's principle, that is: the ray direction is the normal direction of the wave front (isochronous surface), and the ray path is always perpendicular to the wave front (isochronous surface). Therefore, in the travel-time field, the gradient direction at a certain point is the advancing direction of the ray passing through that point.

[0053] Figure 1 It is a schematic diagram of single-point tracking for backtracking. Among them, the two squares represent two grids, point P is the current tracking point, PA and PB represent the partial derivatives in the vertical and horizontal directions of point P in the travel-time field, PC represents the gradient of point P, PD is the negative gradient direction, and D is the next tracking point.

[0054] Based on this, to obtain the ray path between two points, it can be implemented according to the following steps:

[0055] (1) Take the geophone point as the first tracking point on the path (as shown by point P in Figure 1 );

[0056] (2) Calculate the partial derivatives in the vertical and horizontal directions of this point in the travel-time field (as shown by PA and PB in Figure 1 ), and the gradient (as shown by PC in Figure 1 );

[0057] (3) Take the negative gradient direction (as shown by PD in Figure 1 ), and find the first intersection point with the grid along this direction (the intersection point with the side or exactly the grid point) (as shown by point D in Figure 1 ), as the next tracking point on the path;

[0058] (4) Repeat steps (2) and (3) until the source point is traced.

[0059] 1.2 Basic principle of the piecewise iteration method

[0060] The basic theoretical basis of the piecewise iteration method is Snell's law. Considering that any three consecutive points on the same ray path satisfy Snell's law, the source point and the geophone point can be connected by a straight line first. There are several intersection points between the connection line and the interfaces passed through by the path, which are used as the initial tracking results; then, according to Snell's law, the middle point of all three consecutive points on the ray is corrected, that is, the position relationship of all three consecutive points satisfies Snell's law; repeat this step until the sum of all correction amounts is less than a preset threshold, and it is considered that the ray tracing is completed.

[0061] As shown in Figure 2As shown in the figure, z = z2 represents an interface. The velocities above and below the interface are v1 and v2 respectively. P1 and P3 are tracking points on both sides of the interface and adjacent to the interface. Connect the two points with a straight line, which intersects the interface at point P2. P1, P2, and P3 are the initial tracking results. Let point P′2 be the true transmission point. Then the correction amount Δx from P2 to P′2 can be expressed as:

[0062]

[0063] In the formula, a = x2 - x1, b = x3 - x2.

[0064] 2. Error Analysis of Conventional Ray Tracing Method

[0065] 2.1 Error Analysis of Backward Tracing Method

[0066] The implementation of the backward tracing method includes two key points. One is the calculation of the travel time field, and the other is the calculation of the gradient of a certain point in the travel time field. For the first key point, the calculation of the travel time field itself is also a relatively complex problem. In all the examples of the present invention, the fast marching method is used for the numerical simulation of the travel time field. For the second key point, since both the velocity model and the travel time field are spatially discretized, the differential in the gradient calculation needs to be replaced by a difference. And replacing the differential with a difference actually implies the assumption that "the travel time field is linearly distributed locally".

[0067] According to the principle and process of the backward tracing method, the following error factors of this method are analyzed: One is that in addition to the positions of the source and the geophone, only the velocity model is input in the ray tracing process, but the final ray tracing is completed in the travel time field. Therefore, the error of the final ray path not only comes from the method itself, but also from the calculation result of the travel time field in the previous step, and the near-field error of the travel time field is greater than the far-field error. The other is that the assumption that "the travel time field is linearly distributed locally" does not always hold. First, the actual geological situation will cause deviations in the tracking path, and this kind of error will accumulate continuously during the tracking process, resulting in the greater the error accumulation closer to the source position. Second, the curvature of the wavefront is smaller in the far field and closer to a plane wave, so the assumption condition is easier to hold in the far field than in the near field. Based on the above error factors, the error of the ray tracing result obtained by this method near the source should be much larger than the area far from the source.

[0068] Use numerical simulation to test the backward tracing method. Adopt as Figure 3(a) The forward model shown: homogeneous model, with the velocity set to 5000 m / s, and only one source point is set, located at S(2500 m, 3000 m); 10 geophones (R1 - R10) are evenly arranged within the range of 50 m - 4950 m on the surface; a square grid is used for spatial discretization, with the grid size of 800*500 and the grid spacing of 10 m.

[0069] Figure 3 (b) shows the comparison of the tracing results between the backtracking (dashed line) and the analytical solution (solid line) on the entire model scale. There is no obvious error visible to the naked eye between the two in the figure. Zoom in on the area near the source for this result, as Figure 3 (c) shows that there is an obvious error between the backtracking result and the analytical solution, and it has two characteristics: one is that the ray shows a zigzag bend, and the other is that the paths of multiple rays merge.

[0070] The above example intuitively demonstrates the overall effect of the backtracking method and the error characteristics near the source. In order to further quantify the error distribution, on the basis of the forward model shown in Figure 3 (a), the geophones on the surface are increased to 500, and the other parameter settings remain unchanged, as Figure 4 (a) shows. The true ray (solid line) and the backtracking result (dashed line) from the source to a geophone in the homogeneous model are shown in Figure 4 (b). When quantifying the overall error of this tracing result, the average distance between these two lines is obtained, which is called the overall error of a single ray. When exploring the distribution of the error along the ray path, the distribution of the distance between these two lines is used as the quantification of this error, which is called the distribution error of a single ray. Based on the above two error quantification methods, error analysis is carried out on 500 rays, and the corresponding overall error of a single ray and the distribution error of a single ray statistics are obtained. As Figure 4 (c)(d) shows.

[0071] The overall error of a single ray of the backtracking method is shown in Figure 4 (c). The error distribution is roughly symmetric about the source position, and the error change is relatively smooth; in terms of relative magnitude, when the angle between the ray and the vertical direction is 0 degrees, ±30 degrees, there are obvious local minima in the error, and the error around ±30 degrees is the global minimum; in terms of absolute magnitude, the maximum error is close to 15 m, which is equivalent to the ray deviating from the true position by about 1.5 grid spacings as a whole. The average value of the errors of all geophones is about 9.2 m, which is already close to 1 grid spacing.

[0072] The distribution error of a single ray of the backtracking method is shown in Figure 4 (d). Due to the different total lengths of the rays, there is a blank area in the right part of the figure. Looking from the vertical axis direction, the pattern it shows is the same as that ofFigure 4 (c)Basically consistent; viewed from the horizontal axis direction, the error distribution characteristics conform to the feature that the near-field error is much larger than the far-field error; viewed from the absolute magnitude, the error in the local area of some rays has reached about 25 m, that is, 2.5 grid spacings.

[0073] 2.2 Error analysis of the step-by-step iteration method

[0074] According to the principle and process of the step-by-step iteration method, the following problems exist in this method:

[0075] (1) Since each correction only involves a single interface and the two velocity values on both sides of it, and different grids have different velocity values, this kind of correction cannot be carried out across grids. If cross-grid correction is to be realized, the tracking points close to the grid points need to be first reset to the grid points, and then the tracking points coinciding with the grid points are corrected. This process is as Figure 5 (a) shown, where the dashed line, solid line, and dotted line represent three ray paths. It is necessary to first artificially reset the dashed ray path to the solid ray path and then perform correction to obtain the dotted ray path. Even so, at most one grid can be crossed in one correction, which may lead to too slow convergence speed of this method.

[0076] (2) The grids involved in the correction process are always limited to the part within one to two grid spacings around the ray path and its vicinity. This results in that this method cannot obtain the velocity information of other regions. Therefore, the result of this method only represents a local minimum travel time rather than a global minimum travel time, which violates Fermat's principle. For example, in Figure 5 (b), assuming is much larger than, the ray path from point to point should be the refracted wave as shown by the dashed line. In the step-by-step iteration method, according to the process, it is necessary to first connect the two points with a straight line (the dotted line). For any three adjacent points on this straight line, Snell's law is already satisfied, so the correction amount is directly zero, the iteration cannot proceed, and the correct ray path cannot be obtained.

[0077] (3) The first arrival wave may have situations such as back-bending, as Figure 5 (c) shown, that is, the ray first propagates downward (or upward / left / right) and then upward (or downward / right / left). At this time, at the inflection point where the ray path changes from downward to upward, there must be a situation similar to Figure 5 (d) shown, that is, for three consecutive path points, the ordinate of the middle point is less than the minimum value of the ordinates of the two end points. Similar situations will also occur for other directions of back-bending. However, in Snell's law, when the interface is horizontal or vertical, the horizontal and vertical coordinates of the middle point must be between the horizontal and vertical coordinates of the two end points. That is to say, this situation does not satisfy Snell's law. Therefore, the step-by-step iteration method cannot track the rays of the back-bending wave or similar wave fields in the rectangular grid model.

[0078] Numerical simulations were carried out for the three problems existing in the above-mentioned successive-iteration method, and the backtracking results were used as approximate analytical solutions for comparative analysis. For the first problem, a forward model was constructed as shown in Figure 6 (a). This model Figure 3 divided the model shown in (a) into two layers. The depth of the bottom interface of the shallow layer was 2000 m, and the velocity was 2000 m / s. One geophone R(4950 m, 0 m) was set, and the other parameters remained unchanged. To adjust the computational workloads of the two methods to the same level to ensure the same experimental variables, the number of iterations of the successive-iteration method was set to 120 times. The numerical simulation results of the successive-iteration method are shown in Figure 6 (b). The results of 120 iterations of the successive-iteration method were almost no different from the initial value, so it was considered that the tracking failed.

[0079] For the second problem, a forward model was constructed as shown in Figure 7 (a). On the basis of the model shown in Figure 6 (a), the source and geophone positions were set at S(2000 m, 1500 m) and R(6000 m, 1500 m) respectively, and the number of iterations remained unchanged at 120 times. The numerical simulation results of the successive-iteration method are shown in Figure 7 (b). The first arrival wave between the source and the geophone should be a refracted wave, but the tracking result of the successive-iteration method remained the initial value unchanged, so it was considered that the tracking failed.

[0080] For the third problem, a forward model was constructed as shown in Figure 8 (a). This model was similar to the model shown in Figure 6 (a), changing the two-layer velocity model into a spatially uniformly varying model. The velocity increased linearly along the arrow direction. The velocity value at (0 m, 0 m) was 2000 m / s, and the velocity value at (5000 m, 5000 m) was 5000 m / s. The source and geophone positions were S(2000 m, 500 m) and R(6000 m, 500 m) respectively, and the number of iterations remained unchanged at 120 times. The numerical simulation results of the successive-iteration method are shown in Figure 8 (b). The first arrival wave between the source and the geophone should be a diffracted wave, but the result of the successive-iteration method still remained the initial value unchanged, so it was considered that the tracking failed.

[0081] Based on the three problems existing in the successive-iteration method and their corresponding numerical simulation analysis results, it can be seen that in the rectangular grid model, the conventional successive-iteration method cannot give the correct result within a reasonable time, and even in some cases, it cannot effectively update the iteration of the initial value.

[0082] 3 Joint ray tracing method based on multi-scale, multi-direction successive iteration and backtracking

[0083] In the rectangular grid model, the conventional step-by-step iteration method has three problems: slow convergence, local minimum, and inability to trace back the diffracted wave. The main reasons for these three problems are as follows: there are too many interfaces in the model (all grid boundaries are regarded as interfaces); the initial ray is quite different from the true ray; the interfaces are all horizontal and vertical, as well as the theoretical limitations of Snell's law. For the above reasons, the present invention will improve the conventional step-by-step iteration method.

[0084] 3.1 Multi-scale step-by-step iteration method

[0085] Since there are too many interfaces in the rectangular grid, resulting in slow convergence of the step-by-step iteration method. To solve this problem, based on the idea of thinning the rectangular grid to reduce the number of interfaces in the model and thus improve the calculation efficiency of the step-by-step iteration method, a multi-scale step-by-step iteration method is proposed. This method first thins the grid to the preset sparsest level and only retains the path points on the remaining grid boundaries. On this basis, the conventional step-by-step iteration is carried out; then the grid is gradually encrypted, and the path points on the corresponding grid boundaries are added, and the conventional step-by-step iteration is carried out at each sparse level; finally, the above steps are repeated until the grid is encrypted to the original state. The specific process is as follows:

[0086] (1) In the original grid, the initial ray and the tracking points are given. Taking Figure 9 (a) as an example, in a square grid model with 9*9 grid points, it is required to obtain the ray path between P1 and P 16 . First, connect P1 and P 16 with a straight line, and interpolate to find the tracking points P2 to P 15 on the path.

[0087] (2) Set the thinning level as n and thin the grid to the sparsest level, only retaining the path points on the remaining grid boundaries. The thinning method adopted in this method is to retain one line every other line in the current grid model for each additional level of thinning. n-level thinning is equivalent to retaining one line every 2 n lines. For example, for Figure 9 (a), perform 2-level thinning, and the thinning result is as shown in Figure 9 (b). Among them, the 9*9 grid lattice becomes 3*3, and only two of the 14 tracking points P2 to P 15 on the path, namely P8 and P9, remain. They are re-numbered as P2 and P3, and P4 is the previous endpoint P 16 .

[0088] (3) Apply the conventional step-by-step iteration method to the current grid. For example, perform the conventional step-by-step iteration on Figure 9 (b), and the result is as shown in Figure 9 (c).

[0089] (4)Perform one encryption on the current grid. The added grid lines are located in the middle of two adjacent grid lines, which is equivalent to changing the thinning level to level n - 1; then use interpolation to obtain new tracking points. For example, for Figure 9 (c)Perform one encryption, and the result is as shown in Figure 9 (d). The grid dot matrix becomes 5 * 5. The three points P2, P3, and P4 before encryption correspond to P4, P5, and P8 after encryption, while the four points P2, P3, P6, and P7 after encryption are new tracking points.

[0090] (5)Perform the conventional step-by-step iteration method on the grid at the current thinning level. For example, for Figure 9 (d)Perform the conventional step-by-step iteration, and the result is as shown in Figure 9 (e).

[0091] (6)Repeat steps (4) and (5) until the current thinning level drops to 0, that is, the same as the original grid model, and then perform the last pass of the conventional step-by-step iteration.

[0092] Use the multi-scale step-by-step iteration method to perform numerical simulation on the model shown in Figure 6 (also Figure 10 (a)). Considering that the initial ray is quite different from the true ray, a relatively large maximum thinning level is given, which is set to 6 here, and the number of iterations for each level is 20 times, and the total number remains 120 times unchanged. The tracking effect is as shown in Figure 10 (b). Compared with the conventional step-by-step iteration method, the multi-scale step-by-step iteration method has a significant improvement in the tracking effect, and there is a large overlap with the reverse tracking result. It can be considered that its accuracy has been improved to a level close to that of the reverse tracking method.

[0093] 3.2 Multi-direction step-by-step iteration method

[0094] Due to the existence mode of the interface in the rectangular grid model and the theoretical limitations of Snell's law, the step-by-step iteration method cannot track the ray path with a turning-back, as shown in Figure 11 (a). For this problem, the present invention considers the characteristics of the interface in the rectangular grid model, that is: the grid is only the result of spatial discretization of the medium model, and the grid boundary does not represent the real boundary in the actual medium. Its real boundary information is hidden in a large number of grids, and changing the way of spatial discretization will not destroy this information. According to the above characteristics, by changing the direction of the grid boundary, thus changing the way of spatial discretization, to make the turning-back path become a path that satisfies Snell's law, as shown in Figure 11 (b). For the convenience of observation, Figure 11 (b)rotates the ray, rather than the grid boundary, which is equivalent to rotating the grid. After the rotation process, the points that originally did not satisfy Snell's law can be updated using the step-by-step iteration method.

[0095] Based on the above idea, the multi-directional segment-by-segment iteration method is proposed, where the "direction" refers to the relative direction of the grid and the ray. The specific process of this method is as follows:

[0096] (1) Perform segment-by-segment iteration in the original grid (such as Figure 12 (a)), and skip the tracking points that do not satisfy Snell's law.

[0097] (2) Rotate the original grid and the tracking points in it clockwise around the origin by an arbitrary angle within the range of 0 - 90 degrees (in the text, taking a 45-degree rotation as an example), as shown in Figure 12 (b).

[0098] (3) Cover the rotated original grid with a horizontal-vertical grid having the same grid spacing as the original grid, as shown in Figure 12 (c).

[0099] (4) Use interpolation to calculate the intersection points of the rotated ray and the new grid boundary, that is, the tracking points in the new grid.

[0100] (5) Translate the new grid and the tracking points in it to the non-negative region, as shown in Figure 12 (d).

[0101] (6) Perform segment-by-segment iterative update on the rays in the new grid, and skip the points that do not satisfy Snell's law at the same time. Among them, due to the change in the relative direction of the grid boundary and the ray, the points that do not satisfy Snell's law in the new grid are different from those in step (1).

[0102] (7) Transform the new grid and the corresponding tracking points back to the original grid and the corresponding tracking points.

[0103] The above steps are a complete iteration of the multi-directional segment-by-segment iteration method, which actually includes two iterations corresponding to the original grid and the rotated grid.

[0104] Use the multi-directional segment-by-segment iteration method to perform numerical simulation on the Figure 8 model shown (also Figure 13 (a)). Select the original grid and the 45-degree grid as the parameters of the multi-directional method. To accelerate the convergence speed, the multi-scale method is added in this test, and the parameters are set as 6-level thinning and 20 iterations for each level (20 iterations for each of the original grid and the 45-degree grid), and the other parameters remain unchanged. The tracking results are shown in Figure 13 (b), where only the area where the ray path is located is shown. To ensure that the tracking results actually come from the multi-directional segment-by-segment iteration method, Figure 13(b) also gives the results obtained by only using the multi-scale method for comparison. The results of the multi-directional and multi-scale segment-by-segment iteration method have a high coincidence with the results of backtracking, which proves the effectiveness of this method; while simply using the multi-scale segment-by-segment iteration method has no effect at all because the initial ray is already perpendicular to all the grid boundaries it passes through. 3 Joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking

[0105] The multi-scale and multi-directional segment-by-segment iteration method can obtain an accuracy close to that of the backtracking method in the absence of local minima. In order to improve the accuracy of the initial ray of the improved segment-by-segment iteration method to solve the local minimum problem and further improve the convergence speed, the present invention uses the results of the backtracking method as the initial value of the multi-scale and multi-directional segment-by-segment iteration method, and proposes a joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking. This new joint ray tracing method can be considered as using the multi-scale and multi-directional segment-by-segment iteration method to perform a secondary correction on the tracking results of the backtracking method.

[0106] As Figure 14 shown, the specific process of the joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking includes the following steps:

[0107] (1) First, obtain the tracking results of the conventional backtracking method;

[0108] (2-1) Use the tracking results of the conventional backtracking method as the initial ray and tracking points of the segment-by-segment iteration method;

[0109] (2-2) Set the thinning level to n and thin the original grid to the sparsest level;

[0110] (2-3) Set the number of iterations to m. If m > 0 is satisfied, then proceed to step (2-4); otherwise, proceed to step (4-1)

[0111] (2-4) The number of iterations m becomes m - 1; proceed to (3-1)

[0112] (3-1) Apply the conventional segment-by-segment iteration method to the thinned grid and skip the points that do not satisfy Snell's law;

[0113] (3-2) Rotate the thinned grid and the tracking points in it clockwise by an arbitrary angle within the range of 0 - 90° around the origin;

[0114] (3-3) Cover the rotated thinned grid with a new grid having the same grid spacing as the thinned grid;

[0115] (3-4) Use interpolation to calculate the tracking points in the new grid (the intersection points of the rotated ray and the new grid boundary);

[0116] (3 - 5) Translate the new grid and the tracking points therein to the non - negative region for piece - by - piece iterative update and skip the points that do not satisfy Snell's law;

[0117] (3 - 6) Transform the new grid and the corresponding tracking points back to the original decimated grid and the corresponding tracking points;

[0118] (3 - 7) Set n = 0. If so, perform step 2 - 3; otherwise, perform (3 - 8);

[0119] (3 - 8) Set n = n - 1, encrypt the current decimated grid once (the added grid lines are in the middle of two adjacent grid lines) and obtain new tracking points by interpolation; continue with step (2 - 3);

[0120] (4 - 1) Set n = 0. If not, perform step (4 - 2); if so, perform step (4 - 5);

[0121] (4 - 2) Set the decimation level n to n - 1; continue with (4 - 3)

[0122] (4 - 3) Encrypt the current decimated grid once (the added grid lines are in the middle of two adjacent grid lines) and obtain new tracking points by interpolation;

[0123] (4 - 4) Perform the conventional piece - by - piece iterative method on the grid of the current decimation level; continue with step (4 - 1);

[0124] (4 - 5) Perform the conventional piece - by - piece iteration,

[0125] (4 - 6) Output the tracking result;

[0126] (4 - 7) End.

[0127] Use the joint ray - tracing method based on multi - scale, multi - direction piece - by - piece iteration and back - tracking to Figure 7 the model shown (also Figure 15 (a)) for numerical simulation. The results are as shown in Figure 15 (b, c). The results show that the joint ray - tracing method avoids local minima and tracks the refraction wave path. When the ray is in the horizontal segment, there are several small - amplitude bends in the back - tracking method, while the joint ray - tracing method has been significantly improved. To ensure that the reason for avoiding local minima comes from the joint ray - tracing method, Figure 15(b) also gives the results of using the multi-directional and multi-scale segment-by-segment iterative method. The tracking results still have no difference from the initial values, so it is considered that the method fails in tracking. In order to perform error analysis and statistics on the results of the combined ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking, the present invention initializes the rays and tracking points with the numerical simulation results of the backtracking method, sets the parameters as 45-degree rotation, 5-level decimation, and 10 iterations per level, and performs numerical simulation on the Figure 3 and Figure 4 shown models.

[0128] The tracking results of the combined ray tracing method are as shown in Figure 16 (a). Compared with the numerical simulation results of the backtracking method as shown in Figure 3 (c), even if the source area is locally enlarged in the results of the combined ray tracing method, as shown in Figure 16 (b), it still fits well with the analytical solution, effectively reducing the tracking error and significantly improving the accuracy. The overall error statistics of a single ray of the combined ray tracing method are as shown in Figure 16 (c). Compared with the backtracking method as shown in Figure 4 (c), the absolute magnitude of the overall error of a single ray of 500 rays is reduced to less than 5 m, the average error is reduced to 0.86 m, and the accuracy of the overall tracking results of a single ray is improved by 90.65%. The statistical results of the distribution error of a single ray of the combined ray tracing method are as shown in Figure 16 (d). In most areas, the error is close to 0, and in very few areas, the absolute magnitude of the error exceeds 5 m. There is basically no area where the absolute magnitude of the difference exceeds 10 m. Through the above error result analysis, it is verified that the combined ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking improves the accuracy of the conventional backtracking method by about one order of magnitude.

[0129] In summary, in order to solve the three problems existing in the application of the conventional segment-by-segment iterative method in rectangular grids: too slow convergence speed, possible entrapment in local minima, and inability to track the diffracted wave, the present invention proposes three targeted methods or strategies: (1) propose a multi-scale method, first use the decimated model to construct a rough ray path, and then encrypt the model to depict details, which effectively improves the convergence speed of the segment-by-segment iterative method; (2) propose a multi-directional method, by changing the direction of the grid boundary, adjust the rays such as diffracted waves that do not satisfy Snell's law in the original grid to a state where they can be updated using the segment-by-segment iterative method; (3) use the results of the backtracking method as the initial values of the segment-by-segment iterative method to avoid entrapment in local minima during the iterative process. From the perspective of the backtracking method, it can also be considered that the present invention uses the segment-by-segment iterative method as a post-processing process to further correct the results of the backtracking method. Therefore, the present invention effectively overcomes various disadvantages in the prior art and has high industrial utilization value.

[0130] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, the present invention may have various changes and modifications. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking, characterized in that, It includes the following steps: S1. Use the tracking result of the conventional backtracking method as the initial value of the piecewise iteration method; S2. Thin the original grid to the preset sparsest level, and only retain the path points on the boundary of the remaining original grid; S3. Apply conventional piecewise iteration to the thinned grid and change the direction of the boundary of the thinned grid. By changing the way of spatial discretization, make the folded path become a path that satisfies Snell's law; Specifically, it includes the following processes: 3-1. Conduct piecewise iteration in the thinned grid and skip the tracking points that do not satisfy Snell's law; 3-2. Rotate the thinned grid and the tracking points in it clockwise around the origin; 3-3. Cover the rotated thinned grid with a horizontal-vertical grid with the same new grid spacing as the thinned grid; 3-4. Use interpolation to calculate the intersection points of the rotated ray and the new grid boundary, that is, the tracking points in the new grid; 3-5. Translate the new grid and the tracking points in it to the non-negative region; 3-6. Update the piecewise iteration of the ray in the new grid, and skip the points that do not satisfy Snell's law at the same time; 3-7. Transform the new grid and the corresponding tracking points back to the original grid and the corresponding tracking points; S4. Gradually encrypt the thinned grid. Repeat step S3 each time it is encrypted until the current thinned level drops to 0, and then perform the last pass of conventional piecewise iteration.

2. The joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking according to claim 1, characterized in that The specific process of step S2 includes the following: 2-1. In the original grid, given the initial ray and tracking points; 2-2. Set the thinning level to n and thin the grid to the sparsest level, only retaining the path points on the boundary of the remaining grid.

3. The joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking according to claim 2, characterized in that In step 2-2, the thinning method is to retain one line out of every two lines in the current grid model for each additional level of thinning. n levels of thinning is equivalent to retaining one line out of every 2 n lines.

4. The joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking according to claim 1, characterized in that In step 3-2, the rotation angle range is 0-90°.

5. The joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking according to claim 4, characterized in that In step 3-6, due to the change in the relative direction of the grid boundary and the ray, the points that do not satisfy Snell's law in the new grid are different from those in step 3-1.

6. The joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking according to claim 1, characterized in that Step S4 gradually encrypts the thinned grid into a new grid and uses interpolation to obtain new tracking points, and at the same time adds the path points on the corresponding new grid boundary, and performs step S3 at each sparse level.

7. The joint ray tracing method based on multi-scale, multi-directional segment-by-segment iteration and backtracking according to claim 6, characterized in that The specific process of step S4 includes the following: 4-1. Encrypt the current thinned grid once. The added grid lines are located in the middle of two adjacent grid lines, and then use interpolation to obtain new tracking points; 4-2. Perform step S3 on the grid at the current thinned level; 4-3. Repeat steps 4-1 and 4-2 until the current thinned level drops to 0, and then perform the last pass of conventional piecewise iteration.

Images