Finite difference seismic travel time calculation method based on Fibonacci heap sorting

Abstract

The invention discloses a finite difference seismic travel time calculation method based on Fibonacci heap sorting, and the method comprises the following steps: S1, reading parameters and a speed model, and reading related parameters and the speed model in a file; s2, a seismic source and a narrow band are initialized; s3, selecting and arranging a Fibonacci pile; s4, selecting a travel time minimum point in the narrow band; s5, judging attributes of adjacent points; s6, judging whether the narrow band is empty or not; and S7, outputting first-motion wave travel time, performing circulation until a narrow band is empty so as to complete travel time calculation of all grid nodes, and outputting a final travel time calculation result. According to the finite difference seismic wave travel time calculation method based on the Fibonacci heap sorting, the minimum travel time point in the narrow band is screened by adopting the Fibonacci heap sorting, so that the finite difference seismic wave travel time is calculated. The efficiency of grid point travel time calculation is improved, and the finite difference seismic travel time calculation method based on Fibonacci heap sorting is realized.

Description

Technical Field

[0001] The invention relates to the technical field of centering devices, and in particular to a finite difference seismic wave travel time calculation method based on Fibonacci heap sorting. Background Art

[0002] The 5th issue of "Progress in Geophysics" in 2022 published the article "High-precision travel time calculation of complex models using neural networks based on physical information constraints" by Shao Shuai et al., which combines neural networks based on physical information constraints and body-fitting grids to propose a method for effectively calculating the travel time of seismic waves in complex models; an arbitrary complex geological body is constructed using a body-fitting grid, and then a neural network based on physical information constraints is used to solve the travel time of seismic wave propagation on the complex geological body. The residual neural network is used to replace the original neural network architecture, which improves the learning ability of the neural network to a certain extent.

[0003] The 11th issue of "Chinese Journal of Geophysics" in 2009 published "Calculation of seismic wave travel time under complex surface conditions based on linear interpolation and narrowband technology", which proposed a new method for travel time calculation based on linear interpolation and narrowband technology. Among them, linear interpolation is used for local travel time calculation, and narrowband technology is used for local wavefront capture and tracking. In order to approximate the undulating surface, the velocity model is divided by combining triangulated and rectangular networks. The relevant programming practice and numerical experiments show that the new method can not only effectively and flexibly handle the drastic changes in surface elevation, but also has good adaptability and stability. The calculation results obtained meet the law of wavefront propagation.

[0004] The sixth issue of "Journal of Seismology" in 2004 published "Finite Difference Traveltime Calculation Method Based on Rectangular Grids" by Li Zhenchun et al., which derived a finite difference first-arrival wave traveltime calculation method based on arbitrary rectangular grids and local plane wavefront approximation. In addition, this method makes appropriate treatments for the first wave and scattered waves, and does not encounter problems such as shadow areas and caustics that exist in traditional ray methods. The results of simple model and complex Marmousi model trials show that this method has high accuracy and is suitable for complex media with strong longitudinal and lateral velocity changes. The Kirchhoff prestack depth migration based on this method basically achieves the position imaging effect of the wave equation method prestack depth migration in terms of the imaging effect of the main structure and target layer position.

[0005] The fifth issue of Progress in Exploration Geophysics in 2007 published the paper “Calculation of travel time of seismic waves on undulating surfaces by rapid wavefront advancement method” by Sun Zhangqing et al. Since conventional travel time calculation is basically based on horizontal surface, this limits the application of the method in undulating surface areas such as mountainous areas and hills. Based on FMM, the authors proposed a method for calculating the travel time of seismic waves on undulating surfaces and undulating interfaces under complex geological conditions. The calculation results of the theoretical model show that this method is feasible for travel time calculation under complex geological conditions.

[0006] From the above examples, it can be seen that the existing method can guarantee the calculation accuracy of seismic wave travel time to a certain extent, but the above method uses binary heap sorting to select the minimum travel time point in a narrow band. When targeting large seismic wave velocity models, it is difficult to guarantee the efficiency of the travel time calculation method. Summary of the invention

[0007] The purpose of the present invention is to provide a finite difference seismic wave travel time calculation method based on Fibonacci heap sorting to solve the problems raised in the above background technology.

[0008] To achieve the above object, the present invention provides the following technical solution: a finite difference seismic wave travel time calculation method based on Fibonacci heap sorting, comprising the following steps:

[0009] Step 1: Read in the relevant parameters and velocity model in the file, wherein the parameter file includes the number of grid points, grid spacing, and source location of the velocity model;

[0010] Step 2: Initialize the source and narrowband. The specific operation is as follows: Classify the grid nodes in the computational domain into: receiving points, narrowband points, and far-away points. The receiving point refers to the point in all grid points that has completed the travel time calculation. It is the source point during initialization, and its travel time value is: τ(i x ,j z )=0; narrow band points refer to all points within the narrow band. When initialized, they are 12 grid nodes near the earthquake source. The travel time values ​​of the four nodes above, below, left and right of the earthquake source are:

[0011]

[0012] The travel time values ​​of the four diagonal nodes are:

[0013]

[0014]

[0015] The travel time values ​​at the upper, lower, left and right grid spacings are:

[0016]

[0017] in:( For point (i x ,j z ), l is the grid spacing); the far-away points refer to all grid nodes except the receiving points and narrowband points. The value is set as large as possible during initialization (ensuring that the value is much larger than the travel time value that can be calculated for all points);

[0018] Step 3: Create an empty Fibonacci heap, and insert the travel time value of the initial narrow band midpoint as the initial element into the empty heap. The essence of the Fibonacci heap is a collection of minimum heaps, so when creating an empty heap, each tree in the Fibonacci heap must satisfy the property of the minimum heap: the keywords of all nodes in the heap are greater than the keywords of their parent nodes. The Fibonacci heap created in this way is very efficient for operations that require data merging, and the average time complexity of data operations is significantly better than that of a binary heap.

[0019] Step 4: Use Fibonacci heap sort to expand the narrow band, and extract the minimum node from the created Fibonacci heap, that is, the minimum travel time value in the narrow band. The operation is divided into three steps: 1. Put the child nodes of the minimum node into the root linked list of the heap, and delete the minimum node 2. Select the right sibling node of the minimum node as the minimum node 3. Maintain the properties of the heap; in this way, the point with the smallest travel time value in the narrow band is extracted and assigned as the acceptance point and removed from the narrow band.

[0020] Step 5: Use the Fibonacci heap sorting method to select the minimum travel time point τ(i min ,j min ), and set it as the acceptance point, and then min ,j min ) The four grid points above, below, left and right are used for judgment. The travel time value of the accepted point has been judged to be unchanged. The eikonal equation is established for the narrowband point to recalculate its travel time. The points that are still far away are included in the narrowband and the equation is established to calculate the travel time. The finite difference method is used when establishing the equation to calculate the travel time. The eikonal equation is as follows:

[0021]

[0022] Where τ is the travel time of the seismic wave and v is the propagation speed of the seismic wave.

[0023] The calculation of the first arrival travel time can be visualized as the expansion of a curve or surface (two-dimensional or three-dimensional) along the normal direction of a given velocity. We call such a curve or surface a wavefront. In the process of wavefront expansion, problems such as angle rotation and topological structure change will occur, and it is impossible to guarantee that the calculated travel time is the first arrival travel time. In order to solve this problem, we use the upwind difference method to solve the eikonal equation. The specific form is as follows:

[0024]

[0025] Where τ is the travel time of the seismic wave, v is the propagation speed of the seismic wave, Respectively represent the grid node (i, j), x The forward and backward difference operators of the travel time function τ in the direction and the forward and backward difference operators of τ in the z direction, among which are:

[0026]

[0027] The travel time of the unknown point can be solved by the travel time value of the known point through the above equation.

[0028] Step 6: Insert the recalculated value as a new node into the Fibonacci heap. The insertion method is to directly insert the new node before the minimum node of the original root linked list. If the inserted node is smaller than the minimum node of the original root linked list, the inserted node is updated to the minimum node; then repeat the operation in step 4 to extract the minimum node and assign it as the acceptance point until the narrow band is empty.

[0029] Step 7: Loop until the narrow band is empty, completing the travel time calculation of all grid nodes, and output the final travel time calculation result.

[0030] Compared with the prior art, the present invention has the following beneficial effects:

[0031] Taking into account the propagation law of seismic waves in the medium, the Fibonacci heap sorting method is used to select the minimum travel time point in a narrow band for the large seismic wave velocity model, which improves the efficiency of narrow band expansion by 26% on the basis of binary heap sorting, thereby improving the computational efficiency of the finite difference seismic wave travel time method. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 The present invention is a flow chart of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting.

[0033] Figure 2 Schematic diagram of narrowband extension.

[0034] Figure 3 Comparison of time complexity between Fibonacci heap and binary heap

[0035] Figure 4 Schematic diagram of the implementation process of extracting the minimum travel time point in a narrow band using the Fibonacci heap selection method.

[0036] Figure 5 The results of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting in homogeneous media are shown.

[0037] Figure 6 is the relative error of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting in homogeneous media.

[0038] Figure 7 The results of the finite difference seismic wave travel time calculation method based on Fibonacci heap ordering in layered media are given.

[0039] Figure 8The results of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting for the model containing abnormal geological body are shown in Figure 2.

[0040] Fig. 9 The results of finite difference seismic wave travel time calculation method based on Fibonacci heap sorting are calculated for complex stratigraphic models. DETAILED DESCRIPTION

[0041] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0042] See also Figure 1-9 The present invention provides a finite difference seismic wave travel time calculation method based on Fibonacci heap sorting, which is as follows:

[0043] S1, reading in the relevant parameters and velocity model in the file, wherein the parameter file includes the number of grid points, grid spacing, and source position of the velocity model;

[0044] S2. Initialize the source and narrowband. The specific operation is as follows: Classify the grid nodes in the computational domain into: receiving points, narrowband points, and far-away points. The receiving point refers to the point in all grid points that has completed the travel time calculation. It is the source point during initialization, and its travel time value is: τ(i x ,j z )=0; narrow band points refer to all points within the narrow band, which are the four grid nodes near the earthquake source when initialized, and their travel time values ​​are:

[0045]

[0046] in:( They are points (i x ±1,j z )、(i x ,j z ±1), l-bit grid spacing); the far-away points refer to all grid nodes except the receiving points and narrow-band points. When initializing, the value is set as large as possible (ensuring that the value is much larger than the travel time value that can be calculated for all points);

[0047] S3. Use the Fibonacci heap sorting method to select the point with the smallest travel time in the narrow band. The essence of the Fibonacci heap is a collection of minimum heaps, that is, each tree in the Fibonacci heap satisfies the properties of the minimum heap: the keywords of all nodes in the heap are greater than the keywords of their parent nodes. It has high efficiency for operations that require data merging, and the average time complexity of data operations is significantly better than that of a binary heap. The application of the Fibonacci heap in this problem is as follows:

[0048] First, an empty heap is created, and the known travel time values ​​in the narrow band are inserted into the empty heap as nodes to form a Fibonacci heap; the minimum node is extracted from the Fibonacci heap, that is, the minimum travel time value in the narrow band; the minimum node is deleted, that is, the minimum travel time value in the narrow band is set as the acceptance point.

[0049] Extracting the minimum node is mainly divided into three steps: 1. Put the child nodes of the minimum node into the root linked list of the heap and delete the minimum node 2. Select the right sibling node of the minimum node as the minimum node 3. Maintain the properties of the heap. In the process of maintaining the properties of the heap, the main thing is to maintain the degree of the tree, that is, the degree of each tree in the Fibonacci heap cannot be the same. If there are two trees with the same degree, the two trees will be merged. The principle of merging is to merge the root node with a larger weight into the tree with a smaller weight;

[0050] S4, use the Fibonacci heap sorting method to select the minimum travel time point τ(i min ,j min ), and set it as the acceptance point, and then min ,j min ) The four grid points above, below, left and right are used for judgment. The travel time value of the accepted point has been judged to be unchanged. The eikonal equation is established for the narrowband point to recalculate its travel time. The points that are still far away are included in the narrowband and the equation is established to calculate the travel time. The finite difference method is used when establishing the equation to calculate the travel time. The eikonal equation is as follows:

[0051]

[0052] Where τ is the travel time of the seismic wave and v is the propagation speed of the seismic wave.

[0053] The calculation of the first arrival travel time can be visualized as the expansion of a curve or surface (two-dimensional or three-dimensional) along the normal direction of a given velocity. We call such a curve or surface a wavefront. In the process of wavefront expansion, problems such as angle and topological structure change will occur, and it is impossible to guarantee that the calculated travel time is the first arrival travel time. In order to solve this problem, we use the upwind difference method to solve the eikonal equation. The specific form is as follows:

[0054]

[0055] Where τ is the travel time of the seismic wave, v is the propagation speed of the seismic wave, Respectively represent the grid node (i, j), x The forward and backward difference operators of the travel time function τ in the direction and the forward and backward difference operators of τ in the z direction, among which are:

[0056]

[0057] The travel time of the unknown point can be solved by the travel time value of the known point through the above equation.

[0058] S5. Insert the recalculated value as a new node into the Fibonacci heap. The insertion method is to directly insert the new node before the minimum node of the original root linked list. If the inserted node is smaller than the minimum node of the original root linked list, the inserted node is updated to the minimum node; then repeat the operation of step 4 to extract the minimum node and assign it as the acceptance point until the narrow band is empty.

[0059] S6. Loop until the narrow band is empty, that is, the travel time calculation of all grid nodes is completed, and the final travel time calculation result is output.

[0060] The calculation results and accuracy of the method of the present invention are analyzed and verified by four different models.

[0061] Figure 5 , Figure 6 The calculation results and relative errors of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting in homogeneous media are shown respectively. The model has 401 horizontal grid points and 401 vertical grid points. The horizontal and vertical grid spacings are both 1m. The velocity is 2000m / s, and the source is located at 201m horizontally. It can be seen from the figure that the calculation accuracy of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting is very reliable.

[0062] Figure 7 For layered media, Figure 8 For media containing abnormal bodies, Fig. 9 For complex stratigraphic media, the calculation results of the finite difference seismic wave travel time calculation method based on Fibonacci heap sorting are shown. The model has 401 horizontal grid points and 401 vertical grid points. The velocity is 0-4000m / s, and the horizontal and vertical grid spacing is 1m. The earthquake source is located at 201m horizontally. It can be seen from the figure that the calculation results of this method can well reflect the actual situation of the underground medium.

[0063] The present invention adopts Fibonacci heap sorting to screen the minimum travel time of grid points in a narrow band, improves the calculation efficiency of narrow band expansion for large earthquake wave velocity models, and realizes a Fibonacci heap sorting finite difference earthquake wave travel time calculation method.

[0064] Although the present invention has been described in detail with reference to the aforementioned embodiments, it is still possible for those skilled in the art to modify the technical solutions described in the aforementioned embodiments, or to make equivalent substitutions for some of the technical features therein. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A finite difference seismic wave travel time calculation method based on Fibonacci heap sorting, characterized in that: The following steps are involved: S1. Read in parameters and velocity model, and read in relevant parameters and velocity model in the file; S2. Initialization of the earthquake source and narrowband. To initialize the earthquake source and narrowband, first classify the grid nodes in the computational domain into: receiving points, narrowband points, and far-away points. The receiving points refer to the points in all grid points that have completed the travel time calculation. They are the earthquake source points during initialization, and their travel time values ​​are: τ(i x ,j z )=0; narrow band points refer to all points within the narrow band, which are 12 grid nodes near the earthquake source when initialized; S3. Fibonacci heap selection, create an empty Fibonacci heap, and insert the travel time value of the initial narrow band midpoint as the initial element into the empty heap. The essence of the Fibonacci heap is a collection of a series of minimum heaps, so when creating an empty heap, each tree in the Fibonacci heap must satisfy the properties of the minimum heap; S4, select the minimum travel time point in the narrow band, use Fibonacci heap sorting to expand the narrow band, and extract the minimum node from the created Fibonacci heap, that is, the minimum travel time value in the narrow band; S5. Determine the attributes of adjacent points and use the Fibonacci heap sorting method to select the minimum travel time point τ(i min ,j min ), and set it as the acceptance point, and then min ,j min ) The four grid points above, below, left and right are used for judgment. The travel time value of the accepted point has been judged to be unchanged. The eikonal equation is established for the narrow band point to recalculate its travel time. The points that are still far away are included in the narrow band and the equation is established to calculate the travel time. The finite difference method is used when establishing the equation to calculate the trend. The eikonal equation is as follows: Where τ is the travel time of the seismic wave and v is the propagation speed of the seismic wave. The calculation of the first arrival wave travel time can be visualized as the expansion of a curve or surface (two-dimensional or three-dimensional) along the normal direction of a given velocity. We call such a curve or surface a wavefront. During the expansion process, the wavefront will produce problems such as corners and topological structure changes, and it is impossible to guarantee that the calculated travel time is the first arrival travel time. In order to solve this problem, we use the upwind difference method to solve the eikonal equation. The specific form is as follows: Where τ is the travel time of the seismic wave, v is the propagation speed of the seismic wave, Respectively represent the grid node (i, j), x The forward and backward difference operators of the travel time function τ in the direction and the forward and backward difference operators of τ in the z direction, among which are: The travel time of the unknown point can be solved by the travel time value of the known point through the above equation; S6: Check if the narrow band is empty. Insert the recalculated value as a new node into the Fibonacci heap. The insertion method is to directly insert the new node before the minimum node of the original root list. If the inserted node is smaller than the minimum node of the original root list, update the inserted node to the minimum node. Then repeat the operation of S4 to extract the minimum node and assign it as the receiving point until the narrow band is empty. S7, output the travel time of the first arrival wave, and repeat until the narrow band is empty, that is, the travel time calculation of all grid nodes is completed, and the final travel time calculation result is output.

2. The finite difference seismic wave travel time calculation method based on Fibonacci heap sorting according to claim 1, characterized in that: The parameter file in step S1 includes the number of grid points, grid spacing, and source location of the velocity model.

3. The finite difference seismic wave travel time calculation method based on Fibonacci heap sorting according to claim 1, characterized in that: The travel time values ​​of the four nodes above, below, left and right of the earthquake source in step S2 are: The travel time values ​​of the four diagonal nodes are: The travel time values ​​at the upper, lower, left and right grid spacings are: in: For point (i x ,j z ), l is the grid spacing; far-away points refer to all grid nodes except receiving points and narrowband points, and the value is set as large as possible during initialization.

4. The finite difference seismic wave travel time calculation method based on Fibonacci heap sorting according to claim 1 is characterized in that: The property of the minimum heap in the S3 step is that the keywords of all nodes in the heap are greater than the keywords of their parent nodes. The Fibonacci heap created in this way has high efficiency for operations that require data merging, and the average time complexity of data operations is significantly better than that of a binary heap.

5. The finite difference seismic wave travel time calculation method based on Fibonacci heap sorting according to claim 1, characterized in that: The application of the Fibonacci heap in step S3 is to create an empty heap, insert the known travel time values ​​in the narrow band as nodes into the empty heap to form a Fibonacci heap; extract the minimum node from the Fibonacci heap, that is, the minimum travel time value in the narrow band; delete the minimum node, that is, set the minimum travel time value in the narrow band as the acceptance point.

6. The finite difference seismic wave travel time calculation method based on Fibonacci heap sorting according to claim 1, characterized in that: The method for extracting the minimum node in step S4 is: S4.

1. Put the child nodes of the smallest node into the root linked list of the heap and delete the smallest node; S4.2, select the right sibling node of the minimum node as the minimum node; S4.

3. Maintain the properties of the heap; thus, extract the point with the smallest travel time value in the narrow band, assign it as the acceptance point and move it out of the narrow band.

Images