Method for rapidly tracking shortest path of seismic rays in anisotropic well

Abstract

The invention discloses an anisotropic well seismic ray shortest path fast tracking method comprising the following steps: S1, data preparation: S11, picking up a velocity gradient maximum value and a velocity gradient minimum value; s12, dividing the speed model into uniform grids at equal intervals; s13, initializing an observation system; s14, initializing a grid node; s2, processing a single grid obtained by subdivision, and constructing a distance table of a single grid node, an inclination angle table of ray paths of two grid nodes in the grid node and two trigonometric function tables of different ray inclination angles; determining an effective propagation range of the sub-wave source; and S3, carrying out anisotropic ray tracing. According to the method, the calculation speed of the in-well earthquake shortest path ray tracing algorithm of the anisotropic medium is greatly improved on the premise of ensuring that the precision is not influenced.

Description

Technical Field

[0001] The present invention belongs to the technical field of oil and gas exploration. Aiming at a complex structural model, a fast tracing method for the shortest path of anisotropic borehole seismic rays is proposed. Background Technique

[0002] Ray tracing technology is an important part of seismic research. Traditional methods include wavefront expansion methods based on the eikonal equation, Fermat's principle, and Huygens' principle, the symplectic geometric algorithm of seismic rays based on the Hamiltonian system, the shooting method, and the shortest path algorithm, etc. The shortest path algorithm is based on the method theory of graph theory. It uses the subdivided grid nodes to represent the depth model. By calculating the travel time of seismic waves between adjacent grid nodes, and then according to the calculated travel time information, it selects the grid nodes with the minimum travel time connecting the source point to the detector, and its connection line is used as the propagation path of the ray seismic wave. When the calculation is completed, a shortest path "tree" can be obtained, which records the shortest time received by all grid nodes and the seismic propagation path. The shortest path algorithm is stable and flexible, very suitable for various complex geological media, and can simultaneously trace the travel times of various types such as the first arrival wave, reflected wave, and refracted wave of multiple detectors. Therefore, this algorithm is applicable not only to isotropic media but also to ray tracing in anisotropic media.

[0003] The shortest path algorithm expresses the medium model with a series of regular spatial discrete grids and grid nodes. Then, according to Huygens' principle, the discrete grid nodes after spatial subdivision are regarded as sub-wave sources at one time, and the omnidirectional propagation of sub-wave rays is characterized by a finite number of discrete directions, so as to obtain the global minimum travel time of discrete points in the model space and the corresponding path of wave field propagation. Through a large number of research practices, it has been proved that the shortest path is unconditionally stable and can adapt to any dimension and complex medium model. However, the accuracy of the calculation results of this method is seriously affected by the grid accuracy. When the subdivided grid is sparse, the sparse space and direction discretization not only cause large deviations in the calculated travel time of seismic waves and the propagation path position, but also the incident angle and exit angle formed by the connection line of grid nodes are difficult to satisfy continuous changes. Especially when the subdivided grid is too large, the incident angle and exit angle calculated in the ray path are difficult to satisfy Snell's law. Reducing the size of the subdivided grid and increasing the number of grid nodes can significantly improve the accuracy of the results obtained by the shortest path method, but it will significantly increase the calculation time of seismic wave propagation ray tracing. More importantly, the energy of seismic waves propagates in the direction of group velocity in anisotropic media. Calculating the travel time between grid nodes based on the shortest path method mainly calculates the group velocity in the connection line direction of two grid nodes. Given the group angle, calculating the group velocity of a given group angle significantly increases the operation time. Therefore, for the travel time calculation required for migration imaging in anisotropic media, the existing ray tracing methods based on the shortest path are difficult to meet the requirements. Summary of the Invention

[0004] The object of the present invention is to overcome the deficiencies of the prior art and provide a fast anisotropic shortest path ray tracing method that greatly improves the calculation speed of the shortest path ray tracing algorithm for anisotropic media without affecting the accuracy.

[0005] The object of the present invention is achieved by the following technical solutions: A fast tracing method for the shortest path of anisotropic borehole seismic rays, comprising the following steps:

[0006] According to one aspect of an embodiment of the present invention, an anisotropic in-well seismic ray shortest path tracking method is provided, including the following steps: S1. Data preparation, including the following sub-steps: S11. Read the velocity model, calculate the velocity gradient in the velocity model, and pick up the maximum and minimum velocity gradients; S12. Based on the maximum velocity gradient that can be tolerated at the lowest within the same grid, divide the velocity model into equally spaced uniform grids; S13. Initialize the observation system, perform normalization processing on the coordinates of the observation system. If the grid length is m, then any point (x, y) on the grid is normalized to (x / m, y / m); S14. Initialize grid nodes: Establish uniformly spaced grid nodes on all sides of the grid, and number the grid nodes. The information stored in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives; S2. Process a single grid obtained by subdivision, including the following sub-steps: S21. Construct a distance table for single grid nodes: Calculate the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table; S22. Construct an inclination angle table for the ray paths between two grid nodes within the grid node: Calculate the inclination angles of the ray paths between any two grid nodes and store them in a two-dimensional data table; S23. Based on the inclination angle table, construct two trigonometric function tables for different ray inclination angles; S24. Determine the effective propagation range of the sub-wave source; S3. Perform anisotropic ray tracking to calculate the arrival time of the first arrival wave, including the following sub-steps: S31. Starting from the source point, query the velocity of the grid point where the source point is located, and look up the table to obtain the distances from the source point to the surrounding grid nodes and the inclination angles of the paths, and calculate the velocities of the surrounding grid nodes; S32. Perform anisotropic ray tracking on the constructed distance table for single grid nodes, the constructed inclination angle table for the ray paths between two grid nodes within the grid, the constructed trigonometric function tables for different ray inclination angles, and the rule for determining the effective propagation range of the sub-wave source, and save the first arrival information and path information in a new table; S33. Use a linked list to implement the storage structure, establish a queue to store grid nodes, and this queue is used to traverse all grid nodes by the breadth-first algorithm, and perform the enqueue operation on the grid nodes updated in the previous step; S34. According to the breadth-first search method, perform the dequeue operation on the grid nodes in the queue, and use each grid node as a new source point to calculate the travel time t to the adjacent grid nodes n , assuming that the arrival time of the first arrival wave at the current grid node is t, then the arrival time of the first arrival wave propagating through the current grid node to the surrounding grid nodes is t + t n , if this grid node already has an arrival time t n0 , and t n0 < t + t n, no operation is performed. Otherwise, update the first arrival time of the grid node in the table and enqueue the grid node; S35. Repeat the steps for the grid nodes in the queue until the queue is empty. At this time, the seismic wave propagation times and paths from the source point to all grid nodes are recorded in the table.

[0007] According to one aspect of the embodiments of the present invention, an anisotropic in-well seismic ray shortest path tracking device is provided, including: a data preparation module for data preparation, including the following sub-steps: S11. Read the velocity model, calculate the velocity gradient in the velocity model, and pick up the maximum and minimum velocity gradients; S12. Based on the maximum velocity gradient that can be tolerated at the lowest within the same grid, divide the velocity model into equally spaced uniform grids; S13. Initialize the observation system, perform normalization processing on the coordinates of the observation system. If the grid length is m, any point (x, y) on the grid is normalized to (x / m, y / m); S14. Initialize grid nodes: Establish uniformly spaced grid nodes on all sides of the grid and number the grid nodes. The information stored in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives; a single grid processing module for processing a single grid obtained by subdivision, including the following sub-steps: S21. Construct a distance table for the single grid node: Calculate the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table; S22. Construct an inclination angle table for the ray paths between two grid nodes in the grid node: Calculate the inclination angles of the ray paths between any two grid nodes and store them in a two-dimensional data table; S23. Based on the inclination angle table, construct two trigonometric function tables for different ray inclination angles; S24. Determine the effective propagation range of the sub-wave source; an anisotropic ray tracking module for performing anisotropic ray tracking and calculating the arrival time of the first arrival wave, including the following sub-steps: S31. Starting from the source point, query the velocity of the grid point where the source point is located, and look up the table to obtain the distances and path inclination angles from the source point to the surrounding grid nodes, and calculate the velocities of the surrounding grid nodes; S32. Perform anisotropic ray tracking on the constructed distance table for the single grid node, the constructed inclination angle table for the ray paths between two grid nodes in the grid, the constructed trigonometric function tables for different ray inclination angles, and the rules for determining the effective propagation range of the sub-wave source, and save the first arrival information and path information in a new table; S33. Use a linked list to implement the storage structure, establish a queue to store grid nodes. This queue is used for breadth-first algorithm to traverse all grid nodes, and perform the enqueue operation on the grid nodes updated in the previous step; S34. According to the breadth-first search method, perform the dequeue operation on the grid nodes in the queue, and use each grid node as a new source point to calculate the travel time t to the adjacent grid nodes n, assuming the arrival time of the current grid node is \(t\), then the arrival time of the wavefront propagating from the current grid node to the surrounding grid nodes is \(t + t\). n , if the grid node already has an arrival time \(t\). n0 , and \(t\). n0 < \(t + t\). n , then no operation is performed. Otherwise, update the arrival time of the grid node in the table and enqueue the grid node; S35. Repeat the above steps for the grid nodes in the queue until the queue is empty. At this time, the table records the seismic wave propagation time and path from the source point to all grid nodes.

[0008] According to one aspect of an embodiment of the present invention, an electronic device is provided, including: a processor; a memory for storing executable instructions of the processor; wherein, the processor is configured to execute the instructions to implement the anisotropic in - well seismic ray shortest path tracking method as described in any one of the above.

[0009] According to one aspect of an embodiment of the present invention, a computer - readable storage medium is provided. When the instructions in the computer - readable storage medium are executed by a processor of an electronic device, the electronic device can execute the anisotropic in - well seismic ray shortest path tracking method as described in any one of the above.

[0010] The beneficial effects of the present invention are:

[0011] (1) Through pre - calculations of forming a distance table for a single grid node, an inclination angle table for the ray paths between two grid nodes within a grid, and a trigonometric function table for different ray inclination angles, as well as improving the ray calculation method in grid nodes and the strategy of narrowing the effective propagation range of the wavelet source, a large amount of computing resources are saved by looking up tables during the subsequent ray - tracing process.

[0012] (2) Using a simplified anisotropic velocity equation and a linked - list - based storage structure to establish a queue to store grid nodes, and using this queue for breadth - first algorithm to traverse all grid nodes, which lays a foundation for further improving the efficiency of the shortest path ray - tracing algorithm in anisotropic media.

[0013] (3) The present invention greatly improves the calculation speed of the shortest path ray - tracing algorithm in anisotropic media without affecting the accuracy. The ray - tracing algorithm in anisotropic media is a basic tool for techniques such as anisotropic velocity inversion and anisotropic tomography. The improvement of this method will provide stronger technical support for such methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 is a flowchart of the anisotropic shortest path ray - tracing method of the present invention;

[0015] Figure 2 Schematic diagram of the information structure of the grid node of the present invention;

[0016] Figure 3 Schematic diagram of the grid node model of this embodiment;

[0017] Figure 4 Schematic diagram of the grid node update model of this embodiment;

[0018] Figure 5 Shortest path map of borehole seismic rays in a certain work area obtained by the method of the present invention. Detailed implementation manners

[0019] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0020] It should be noted that the terms "first", "second", etc. in the specification and claims of the present invention and the above drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged under appropriate circumstances so that the embodiments of the present invention described herein can be implemented in an order different from those illustrated or described herein. In addition, the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device comprising a series of steps or units is not necessarily limited to those steps or units clearly listed, but may include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0021] According to an embodiment of the present invention, an embodiment of a method for tracing the shortest path of anisotropic borehole seismic rays is provided. The method includes the following steps:

[0022] S1. Data preparation, including the following sub-steps:

[0023] S11. Read the velocity model, calculate the velocity gradient in the velocity model, and pick up the maximum value and the minimum value of the velocity gradient;

[0024] In this step, the velocity model is described. The velocity model refers to the distribution of the acoustic wave propagation velocity of the medium at different depths underground. The velocity gradient in the velocity model refers to the slope of the acoustic wave propagation velocity changing with depth in the underground medium, that is, the change rate of the medium density and elastic modulus with depth. The magnitude of the velocity gradient has an important influence on the propagation path and velocity of seismic waves underground, and has a great guiding role in the properties of underground media and the distribution of oil resources.

[0025] S12. Based on the maximum value of the minimum allowable velocity gradient within the same grid, divide the velocity model into equally spaced uniform grids;

[0026] In this step, the maximum value of the minimum allowable velocity gradient within the same grid is described. Since the maximum gradient value in the velocity model grid refers to the maximum value of the gradient (i.e., the velocity change rate) at the grid points in the velocity model. This usually indicates that there are drastic velocity changes in a certain area of the velocity model.

[0027] The velocity model is usually modeled in the form of grids, that is, the underground medium is divided into multiple grid cells, and each grid cell has a corresponding velocity value. Therefore, the maximum value of the minimum allowable velocity gradient within the same grid can be obtained.

[0028] Dividing the velocity model into equally spaced uniform grids means dividing the velocity model into grid cells of the same size and shape based on the maximum gradient value. This method can help better understand the underground geological structure and velocity distribution.

[0029] S13. Initialize the observation system and normalize the coordinates of the observation system. If the grid length is m, then any point (x, y) on the grid is normalized to (x / m, y / m);

[0030] In this step, the normalization of the coordinates of the observation system is described. By normalizing the coordinates, the value ranges of different data can be unified into the same range, making them comparable and facilitating comparison and analysis. Coordinate normalization can also reduce data overflow and gradient explosion phenomena, improve the stability and robustness of the system. Coordinate normalization can reduce the computational complexity and improve the efficiency and speed of the algorithm.

[0031] S14. Initialize grid nodes: Establish uniformly spaced grid nodes on all sides of the grid and number the grid nodes. The information stored in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives;

[0032] In this step, the grid nodes are numbered, which can make the subsequent processing more orderly. Moreover, the information stored in each grid node includes which edge of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives, which can enrich the stored information and enable further processing based on this sufficient information in the future.

[0033] S2. Process the individual grids obtained by meshing, including the following sub-steps:

[0034] S21. Construct a distance table for individual grid nodes: Calculate the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table;

[0035] In this step, a distance table for each grid is constructed to represent the distances from the current grid node to other grid nodes, clearly showing the distance relationships.

[0036] S22. Construct an inclination angle table for the ray paths between two grid nodes within a grid node: Calculate the inclination angles of the ray paths between any two grid nodes and store them in a two-dimensional data table;

[0037] In this step, an inclination angle table for every two grids is constructed to represent the inclination angles of the ray paths between two grid nodes, clearly showing the angle relationships.

[0038] S23. Based on the inclination angle table, construct two trigonometric function tables for different ray inclination angles;

[0039] In this step, the trigonometric function tables are determined. Among them, the trigonometric function tables can represent the values of trigonometric functions such as sine, cosine, and tangent at different angles. Such tables can be used for the calculation and analysis of trigonometric functions in mathematical and physical problems.

[0040] S24. Determine the effective propagation range of the wavelet source;

[0041] In this step, the effective propagation range of the wavelet source is determined. Determining this propagation range is very important for many applications. For example, in seismic exploration, determining the effective propagation range of the underground wavelet source can help determine the underground structure and resource distribution; in sonar systems, determining the effective propagation range of underwater acoustic sources can help with underwater detection and communication; in medical ultrasound imaging, determining the effective propagation range of ultrasonic sources can help with the diagnosis of tissue structures and lesions.

[0042] S3. Conduct anisotropic ray tracing and calculate the arrival time of the first arrival wave, including the following sub-steps:

[0043] S31. Starting from the seismic source point, query the velocity of the grid point where the seismic source point is located, look up the distances from the seismic source point to the surrounding grid nodes and the dip angles of the paths in a table, and calculate the velocities of the surrounding grid nodes.

[0044] In this step, the velocities of the surrounding grid nodes are determined. The velocities of the surrounding grid nodes can refer to the geological velocities at different positions and depths in the underground structure. These velocity values can be used to describe the propagation velocities of substances such as underground rocks and soils at different depths. Through these data, geological phenomena such as earthquakes, rock structures, and underground water flows can be predicted more accurately.

[0045] S32. Perform anisotropic ray tracing on the constructed single grid node distance table, the constructed table of the dip angles of the ray paths between two grid nodes within the grid, the constructed trigonometric function table of different ray dip angles, and the rules for determining the effective propagation range of the wavelet source, and save the first arrival information and path information in a new table.

[0046] In this step, anisotropic ray tracing is performed and the first arrival information and path information are saved in a new table, enabling a clear understanding of the relationships through the table.

[0047] S33. Use a linked list to implement the storage structure, establish a queue to store grid nodes. This queue is used for breadth - first algorithm to traverse all grid nodes, and perform the enqueue operation on the grid nodes updated in the previous step.

[0048] In this step, through the operations of dequeueing and enqueueing, it can be ensured that no grid nodes are missed. Moreover, the breadth - first algorithm is adopted. Since the breadth - first algorithm can search all nodes, it can find the shortest path or determine the relationships between nodes. The breadth - first algorithm can be used to find the shortest path between two nodes because it traverses nodes layer by layer to find the nearest node. And through the breadth - first algorithm, the relationships between nodes can be determined, such as whether a certain node can reach another node or the distance between two nodes. The breadth - first algorithm traverses nodes by level, so it can avoid getting into an infinite loop situation, thus better processing node information.

[0049] S34. According to the breadth - first search method, perform the dequeue operation on the grid nodes in the queue, and use each grid node as a new seismic source point to calculate the travel time t to the adjacent grid nodes. n Assume that the first arrival time of the current grid node is t, then the first arrival time of the wave propagating through the current grid node to the surrounding grid nodes is t + t. n If this grid node already has a first arrival time t n0 and t n0 <t + t n, no operation is performed; otherwise, update the first arrival time of the grid node in the table and enqueue the grid node;

[0050] In this step, the operations of enqueueing and dequeueing grid nodes are carried out in chronological order to ensure that the shortest path can be determined.

[0051] S35. Repeat the above steps for the grid nodes in the queue until the queue is empty. At this time, the seismic wave propagation time and path from the seismic source point to all grid nodes are recorded in the table.

[0052] In this step, by traversing all grid nodes, the determination of the seismic wave propagation time and path from the seismic source point to all grid nodes is realized.

[0053] As an optional embodiment, calculate the distance from all grid nodes in the grid to any other grid node in the same grid, including: respectively record two grid nodes in the grid as the incident point and the exit point of the seismic wave ray; adopt the following method to obtain the distance from the incident point of the seismic wave ray to the exit point of the seismic wave ray: where len in,out represents the distance from the incident point of the seismic wave ray to the exit point of the seismic wave ray, and x in , y in are the coordinate positions of the incident point of the seismic wave ray; x out , y out are the coordinate positions of the exit point of the seismic wave ray.

[0054] In this embodiment, the specific steps for calculating the distance from all grid nodes in the grid to any other grid node in the same grid are described. That is, the distance between two grid nodes can be determined by this method. By this method, the distance between two grid nodes can be accurately determined.

[0055] As an optional embodiment, determine the effective propagation range of the wavelet source, including: using Snell's law and the velocity of the adjacent grid of the grid node, calculate the direction θ2 of the seismic wave propagation path in the next grid with the ray incident angle θ1 of the previous grid path of the grid node; limit the inclination range of θ2±Δθ as the effective propagation area of the wavelet source, where Δθ is the set angular deviation.

[0056] In this embodiment, the specific steps for determining the effective propagation range of the wavelet source are described. Snell's law describes the refraction law of waves when propagating between media, that is, the relationship between the incident angle and the refraction angle. Determining the propagation range of the wavelet source based on the velocities of adjacent grids of grid nodes in the field of seismic exploration means that, according to Snell's law and the velocities of adjacent grid nodes, the propagation range of the wavelet in the subsurface medium is determined. This process can more accurately understand the structure and properties of the subsurface medium, thereby improving the efficiency and accuracy of exploration. Using Snell's law and the velocities of adjacent grids of grid nodes to determine the propagation range of the wavelet source can improve the accuracy of exploration. Moreover, by accurately determining the propagation range of the wavelet, the seismic exploration plan can be optimized, more appropriate exploration parameters and paths can be selected, and the exploration efficiency can be improved.

[0057] As an alternative embodiment, calculating the velocities of each surrounding grid node includes: querying the velocity of the grid point where the source point is located, and looking up a table to obtain the distance and the dip angle of the path from the source point to the surrounding grid nodes; using the following anisotropic velocity formula to calculate the velocities of each surrounding grid node and then calculating the arrival time of the first arrival wave: V ani = V0(1 + δsin 2 θcos 2 θ + εsin 4 θ); where: V0 is the seismic wave propagation velocity in the direction of the symmetry axis of the vertical transverse isotropic (VTI) medium; δ and ε are the weak anisotropic parameters of Thomas, and θ is the dip angle of the path from the source point to the surrounding grid nodes.

[0058] In this embodiment, the specific steps for calculating the velocities of each surrounding network node are described. Among them, the seismic wave propagation velocity in the direction of the symmetry axis of the VTI medium refers to the velocity of seismic waves propagating along the symmetry axis direction in a medium with vertical transverse isotropy (VTI). This velocity is usually affected by factors such as medium density, rock stiffness, and porosity. In seismic exploration, understanding the propagation velocity of seismic waves in a VTI medium can help to more accurately understand the subsurface structure and geological characteristics. Through the above method, the velocities of each grid node can be determined more accurately.

[0059] As an alternative embodiment, performing dequeue and enqueue operations on the grid nodes in the queue includes: performing a dequeue operation on the grid nodes in the queue, and using each grid node as a new source point to calculate the travel time t n to the adjacent grid nodes. Assuming that the arrival time of the first arrival wave at the current grid node is t, then the arrival time of the first arrival wave propagating through the current grid node to the surrounding grid nodes is t + t n . If this grid node already has an arrival time t n0 , and t n0 < t + t n, no operation is performed; otherwise, update the first arrival time of the grid node in the table and enqueue the grid node.

[0060] In this embodiment, the specific steps of dequeueing and enqueueing grid nodes in the queue are described. By determining the first arrival time to perform dequeueing and enqueueing, the recorded information can be valid and accurate. It can better determine the seismic wave propagation time and path from the seismic source point to all grid nodes.

[0061] As an alternative embodiment, the method is applied to oil and gas exploration.

[0062] In this embodiment, applying the method for determining seismic propagation velocity to oil and gas exploration, by determining the seismic propagation velocity, the location and depth of oil and gas reservoirs can be more accurately located, thereby improving exploration efficiency. Accurate seismic propagation velocity can help reduce the trial-and-error cost in the exploration process and reduce unnecessary exploration risks. Moreover, by determining the seismic propagation velocity, the formation structure and geological features can be more precisely interpreted, thereby increasing the success rate of oil and gas exploration. Accurate seismic propagation velocity can help determine the optimal oil production plan, improve oil and gas production efficiency and output. In addition, through precise determination of seismic propagation velocity, the impact on the environment caused by exploration activities can be reduced, and environmental risks can be reduced.

[0063] Based on the above content, the present application also proposes an alternative implementation manner, which is introduced in detail below:

[0064] Aiming at the problem of low computational efficiency of the shortest path ray tracing algorithm for anisotropic media, the present invention pre-computes the distance table of a single grid node, the dip angle table of the ray path between two grid nodes in the grid, and the trigonometric function table of different ray dip angles, and improves the ray calculation method in the grid node and reduces the effective propagation range of the sub-wave source, solves the difficult problem of low computational efficiency of the conventional ray tracing algorithm, and innovates the fast anisotropic shortest path ray tracing technology by means of a simplified anisotropic velocity calculation equation, realizing a great improvement in the computational speed of the shortest path ray tracing algorithm for anisotropic media on the premise of ensuring that the ray tracing accuracy of anisotropic media is not affected, and laying a foundation for high-precision velocity modeling and imaging of anisotropic media. The technical solution of the present invention is further described below with reference to the accompanying drawings.

[0065] Figure 1 is the flowchart of the anisotropic shortest path ray tracing method of the present invention. As Figure 1 shown, a fast anisotropic wellbore seismic ray shortest path tracing method of the present invention includes the following steps:

[0066] S1. Data preparation, including the following sub-steps:

[0067] S11. Read the velocity model, calculate the velocity gradient in the velocity model, and pick up the maximum and minimum values of the velocity gradient;

[0068] S12. Based on the maximum value of the velocity gradient that can be tolerated at least within the same grid, divide the velocity model into equally spaced uniform grids;

[0069] That is, set the maximum value of the velocity gradient that can be tolerated at least within the same grid in the velocity model, that is, the maximum rate of change of velocity. Based on this criterion, divide the grids in the velocity model into equally spaced uniform grids; In this way, both the accuracy of ray tracing can be guaranteed and the number of grids can be appropriately reduced.

[0070] S13. Initialize the observation system and normalize the coordinates of the observation system. If the grid length is m, then any point (x, y) on the grid is normalized to (x / m, y / m);

[0071] S14. Initialize the grid nodes: Establish uniformly spaced grid nodes on all sides of the grid and number the grid nodes. The information saved in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives. Figure 2 It is a schematic diagram of the information structure of the grid node of the present invention. As Figure 2 shown, node1, node2, and node3 are grid nodes;

[0072] S2. To improve the efficiency of the shortest path ray tracing, process the individual grids obtained by meshing. In this embodiment, taking 3 grid nodes on each side as an example, the grid node model is as Figure 3 shown. Actually, more grid nodes are divided in the calculation, and the calculation method is the same. It includes the following sub-steps:

[0073] S21. Construct a distance table for individual grid nodes: Calculate the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table; The two grid nodes are respectively recorded as the incident point and the exit point of the seismic wave ray, and len in,out is used to represent the distance from the incident point to the exit point:

[0074]

[0075] where, x in , y in are the coordinate positions of the incident point of the seismic wave ray; x out , y out are the coordinate positions of the exit point of the seismic wave ray; Figure 3 It is a schematic diagram of the grid node model of this embodiment. The distances from one grid node to the other grid nodes are as Figure 3As shown in (a), the distances from all grid nodes to the remaining grid nodes are as follows Figure 3 As shown in (b). The obtained distance table is shown in Table 1 (starting from Figure 3 (a), the grid node numbers are sequentially labeled in a clockwise or counterclockwise direction). Table 1 is a distance mapping table for calculating the distances from all grid nodes in the grid to any other grid node in the same grid. As shown in Table 1:

[0076] Table 1

[0077]

[0078] S22. Construct an inclination angle table for the ray paths between two grid nodes within a grid node: Calculate the inclination angles of the ray paths between any two grid nodes and save them in a two-dimensional data table; the inclination angle information can be used to calculate the anisotropic velocity, and the corresponding values can be directly located and queried through the grid node number information, and the query time complexity is O(1). Table 2 is an inclination angle mapping table for calculating the inclination angles of all grid nodes in the grid to any other grid node in the same grid. The obtained inclination angle table in this embodiment is shown in Table 2.

[0079] Table 2

[0080]

[0081]

[0082] S23. On the basis of the inclination angle table, construct two trigonometric function tables for different ray inclination angles;

[0083] S24. Determine the effective propagation range of the wavelet source: When calculating the propagation of seismic waves, the traditional method needs to calculate the travel times of the grid nodes connected by all discrete rays in all directions around the wavelet source. Based on Snell's law, when seismic waves propagate from one grid node to the next grid node, a large part of the directions within the grid node are unnecessary. The present invention uses Snell's law and the velocities of the adjacent grids of the grid node to calculate the path direction θ2 of the seismic wave propagation in the next grid based on the ray incident angle θ1 of the previous grid path of the grid node, and preliminarily determines the advancing direction of the seismic wave propagation ray in the next grid; the inclination angle range of θ2±Δθ is defined as the effective propagation area of the wavelet, and Δθ is a set angle threshold, so as to further improve the shortest path ray tracing algorithm.

[0084] S3. Perform anisotropic ray tracing and calculate the arrival time of the first arrival wave; including the following sub-steps:

[0085] S31. Starting from the seismic source point, query the velocity of the grid point where the seismic source point is located, and look up the table to obtain the distances from the seismic source point to the surrounding grid nodes and the inclination angles of the paths. Calculate the velocities of each surrounding grid node using the anisotropic velocity formula, and then calculate the arrival time of the first arrival wave. In the shortest ray tracing algorithm for anisotropic media, the group velocity is a function of the phase angle. There is often a problem that the group angle is known while the phase angle is unknown when calculating the travel time within the grid. Although the group angle can be expressed as an analytical function of the phase angle, it is difficult to obtain the expression of its inverse function. To simplify the calculation, the anisotropic velocity V ani The calculation formula is as follows:

[0086] V ani = V0(1 + δsin 2 θcos 2 θ + εsin 4 θ);

[0087] Where: V0 is the seismic wave propagation velocity in the direction of the symmetry axis of the transversely isotropic VTI medium with a vertical symmetry axis; δ and ε are the Thomas weak anisotropy parameters. δ indicates the difference between the longitudinal wave velocity in the vertical direction and the longitudinal wave velocity in the horizontal direction, expressing the anisotropy of the longitudinal wave. ε indicates the difference between the shear wave velocity in the vertical direction and the shear wave velocity in the horizontal direction, expressing the anisotropy of the shear wave. Due to a large number of trigonometric function calculations, the efficiency is reduced. Therefore, based on the ray inclination angle table, two trigonometric function tables with different ray inclination angles are constructed for subsequent calculation and query.

[0088] S32. Perform anisotropic ray tracing on the constructed single grid node distance table, the constructed inclination angle table of the ray paths between two grid nodes within the grid, the constructed trigonometric function tables with different ray inclination angles, and the rules for determining the effective propagation range of the wavelet source to obtain the first arrival information and path information, and save the first arrival information and path information in a new table;

[0089] S33. Use a linked list to implement the storage structure, establish a queue to store grid nodes. This queue is used for the breadth-first algorithm to traverse all grid nodes, and perform the enqueue operation on the grid nodes updated in the previous step;

[0090] S34. According to the breadth-first search method, perform the dequeue operation on the grid nodes in the queue, and use each grid node as a new seismic source point to calculate the travel time t n to the adjacent grid nodes. Assume that the first arrival time of the current grid node is t, then the first arrival time of the wave propagating through the current grid node to the surrounding grid nodes is t + t n . If this grid node already has a first arrival time t n0 , and t n0 < t + t n, no operation is performed; otherwise, the arrival time of the grid node is updated in the table, and the grid node is enqueued; Figure 4 is a schematic diagram of the grid node update model in this embodiment. The grid node update mode is as Figure 4 shown: The travel time of the node0-node1-node6-node7 path is shorter than that of the node0-node3-noe7 path. Therefore, according to Fermat's principle, the ray travels along the first path.

[0091] S35. Repeat step 2 until the queue is empty. At this time, the table records the seismic wave propagation time and path from the seismic source to all grid nodes.

[0092] Apply the method of the present invention to the data of a walkaway VSP work area in a certain complex work area to perform fast tracking of seismic rays. Figure 5 is the shortest path map of the borehole seismic rays in a certain work area obtained by the method of the present invention. The shortest travel time paths from each seismic source point to the receiving point reflected in the target layer are as Figure 5 shown.

[0093] According to an embodiment of the present invention, there is also provided an apparatus for implementing the above-mentioned anisotropic borehole seismic ray shortest path tracking method. The apparatus will be described in detail below.

[0094] The device includes: a data preparation module for data preparation, including the following sub-steps: S11, reading the velocity model, calculating the velocity gradient in the velocity model, and picking up the maximum and minimum values of the velocity gradient; S12, dividing the velocity model into equally spaced uniform grids based on the maximum value of the lowest allowable velocity gradient within the same grid; S13, initializing the observation system and normalizing the coordinates of the observation system. If the grid length is m, then any point (x, y) on the grid is normalized to (x / m, y / m); S14, initializing the grid nodes: establishing uniformly spaced grid nodes on all sides of the grid and numbering the grid nodes. The information stored in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives; a single grid processing module for processing a single grid obtained by dissection, including the following sub-steps: S21, constructing a distance table for single grid nodes: calculating the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table; S22, constructing an inclination angle table for the ray paths between two grid nodes within the grid node: calculating the inclination angles of the ray paths between any two grid nodes and storing them in a two-dimensional data table; S23, constructing two trigonometric function tables for different ray inclination angles based on the inclination angle table; S24, determining the effective propagation range of the wavelet source; an anisotropic ray tracing module for performing anisotropic ray tracing and calculating the arrival time of the first arrival wave, including the following sub-steps: S31, starting from the source point, querying the velocity of the grid point where the source point is located, and looking up the distances and path inclination angles from the source point to the surrounding grid nodes to calculate the velocities of the surrounding grid nodes; S32, performing anisotropic ray tracing on the constructed distance table for single grid nodes, the constructed inclination angle table for the ray paths between two grid nodes within the grid, the constructed trigonometric function tables for different ray inclination angles, and the rules for determining the effective propagation range of the wavelet source, and storing the first arrival information and path information in a new table; S33, using a linked list to implement the storage structure, establishing a queue to store grid nodes. This queue is used for breadth-first algorithm to traverse all grid nodes, and performing an enqueue operation on the grid nodes updated in the previous step; S34, according to the breadth-first search method, performing a dequeue operation on the grid nodes in the queue, and using each grid node as a new source point to calculate the travel time t to the adjacent grid nodes n , assuming the arrival time of the first arrival wave at the current grid node is t, then the arrival time of the first arrival wave propagating from the current grid node to the surrounding grid nodes is t + t n , if the grid node already has an arrival time t n0 , and t n0 < t + t n, no operation is performed; otherwise, update the first arrival time of the grid node in the table and enqueue the grid node; S35. Repeat the steps for the grid nodes in the queue until the queue is empty. At this time, the seismic wave propagation time and path from the seismic source point to all grid nodes are recorded in the table.

[0095] Optionally, the single-grid processing module is also used to calculate the distance from any grid node in the grid to any other grid node in the same grid by the following method: Denote two grid nodes in the grid as the seismic wave ray incident point and the seismic wave ray exit point respectively; Obtain the distance from the seismic wave ray incident point to the seismic wave ray exit point by the following method: where len in,out represents the distance from the seismic wave ray incident point to the seismic wave ray exit point, and x in , y in are the coordinate positions of the seismic wave ray incident point; x out , y out are the coordinate positions of the seismic wave ray exit point.

[0096] Optionally, the single-grid processing module is also used to determine the effective propagation range of the wavelet source by the following method: Use Snell's law and the velocities of the adjacent grids of the grid nodes to calculate the direction θ2 of the seismic wave propagation path in the next grid based on the ray incident angle θ1 of the previous grid path of the grid node; Limit the inclination range of θ2±Δθ as the effective propagation area of the wavelet source, where Δθ is the set angular deviation.

[0097] Optionally, the anisotropic ray tracing module is also used to calculate the velocities of the surrounding grid nodes by the following method: Query the velocity of the grid point where the seismic source point is located, and look up the table to obtain the distance from the seismic source point to the surrounding grid nodes and the inclination angle of the path; Calculate the velocities of the surrounding grid nodes and then calculate the first arrival time by the following anisotropic velocity formula: V ani =V0(1 + δsin 2 θcos 2 θ + εsin 4 θ); where: V0 is the seismic wave propagation velocity in the direction of the symmetry axis of the vertically transverse isotropic VTI medium; δ and ε are the Thomas weak anisotropy parameters, and θ is the inclination angle of the path from the seismic source point to the surrounding grid nodes.

[0098] It should be noted here that the above modules correspond to the steps of implementing the anisotropic borehole seismic ray shortest path tracing method. The instances and application scenarios implemented by multiple modules and the corresponding steps are the same, but are not limited to the content disclosed in the above embodiments.

[0099] According to another aspect of the embodiments of the present invention, an electronic device is further provided, including: a processor; and a memory for storing instructions executable by the processor, wherein the processor is configured to execute the instructions to implement the anisotropic in-well seismic ray shortest path tracing method according to any one of the above.

[0100] According to another aspect of the embodiments of the present invention, a computer-readable storage medium is further provided. When the instructions in the computer-readable storage medium are executed by the processor of the electronic device, the electronic device can execute the anisotropic in-well seismic ray shortest path tracing method according to any one of the above.

[0101] The serial numbers of the embodiments of the present invention above are only for description and do not represent the advantages or disadvantages of the embodiments.

[0102] In the above embodiments of the present invention, the descriptions of the respective embodiments have their own emphases. For the parts not detailed in a certain embodiment, reference may be made to the relevant descriptions of other embodiments.

[0103] In the several embodiments provided by the present application, it should be understood that the disclosed technical content can be implemented in other ways. Among them, the device embodiments described above are only illustrative. For example, the division of the units can be a logical function division. In actual implementation, there may be other division methods. For example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the displayed or discussed coupling or direct coupling or communication connection to each other can be through some interfaces. The indirect coupling or communication connection of the units or modules can be in an electrical or other form.

[0104] The units described as separate components may or may not be physically separated. The components displayed as units may or may not be physical units, that is, they can be located in one place or distributed to multiple units. Some or all of the units can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0105] In addition, the functional units in the various embodiments of the present invention can be integrated in a processing unit, or each unit can exist physically alone, or two or more units can be integrated in one unit. The above integrated units can be implemented in the form of hardware or in the form of software functional units.

[0106] When the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes: various media that can store program codes, such as USB flash drives, read-only memories (ROM, Read-Only Memory), random access memories (RAM, Random Access Memory), mobile hard disks, magnetic disks, or optical discs.

[0107] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. An anisotropic in-well seismic ray shortest path tracing method, characterized in that, It includes the following steps: S1. Data preparation, including the following sub-steps: S11. Read the velocity model, calculate the velocity gradient in the velocity model, and pick up the maximum and minimum values of the velocity gradient; S12. Divide the velocity model into equally spaced uniform grids based on the maximum value of the velocity gradient that can be tolerated at the lowest level within the same grid; S13. Initialize the observation system and normalize the coordinates of the observation system. If the grid length is m, then any point (x, y) on the grid is normalized to (x / m, y / m); S14. Initialize the grid nodes: Establish uniformly spaced grid nodes on all sides of the grid, number the grid nodes, and the information stored in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives; S2. Process a single grid obtained by subdivision, including the following sub-steps: S21. Construct a distance table for single grid nodes: Calculate the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table; S22. Construct an inclination angle table for the ray paths between two grid nodes within the grid node: Calculate the inclination angles of the ray paths between any two grid nodes and store them in a two-dimensional data table; S23. Based on the inclination angle table, construct two trigonometric function tables for different ray inclination angles; S24. Determine the effective propagation range of the wave source; S3. Perform anisotropic ray tracing to calculate the arrival time of the first arrival wave, including the following sub-steps: S31. Starting from the source point, query the velocity of the grid point where the source point is located, look up the table to obtain the distances and path inclination angles from the source point to the surrounding grid nodes, and calculate the velocities of the surrounding grid nodes; S32. Perform anisotropic ray tracing on the constructed distance table for single grid nodes, the constructed inclination angle table for the ray paths between two grid nodes within the grid, the constructed trigonometric function tables for different ray inclination angles, and the rules for determining the effective propagation range of the wave source, and save the first arrival information and path information in a new table; S33. Use a linked list to implement the storage structure, establish a queue to store grid nodes, and this queue is used for breadth-first algorithm to traverse all grid nodes, and perform the enqueue operation on the grid nodes updated in the previous step; S34. According to the breadth-first search method, perform dequeue and enqueue operations on the grid nodes in the queue, and use each grid node as a new seismic source point to calculate the travel time t to adjacent grid nodes n , assuming that the first arrival time of the current grid node is t, then the first arrival time of the wave propagating from the current grid node to the surrounding grid nodes is t + t n , if the grid node already has a first arrival time t n0 , and t n0 < t + t n , then do nothing, otherwise update the first arrival time of the grid node in the table and enqueue the grid node; S35. Repeat the steps of the grid nodes in the above queue until the queue is empty. At this time, the earthquake wave propagation times and paths from the source point to all grid nodes are recorded in the table.

2. The method according to claim 1, wherein The calculation of the distances from all grid nodes in the grid to any other grid node in the same grid includes: Denote two grid nodes in the grid as the incident point and the exit point of the seismic wave ray respectively; Obtain the distance from the incident point of the seismic wave ray to the exit point of the seismic wave ray in the following way: wherein, len in,out represents the distance from the incident point to the exit point of the seismic wave ray, and x in , y in are the coordinate positions of the incident point of the seismic wave ray; x out , y out are the coordinate positions of the exit point of the seismic wave ray.

3. The method according to claim 1, wherein The determination of the effective propagation range of the wave source includes: Using Snell's law and the velocities of the adjacent grids of the grid node, calculate the direction θ2 of the seismic wave propagation path in the next grid with the ray incident angle θ1 of the previous grid path on the grid node; The effective propagation region of the wavelet source is defined by the range of the inclination angle θ2±Δθ, where Δθ is the set angular deviation.

4. The method according to claim 1, wherein Calculating the velocities of each surrounding grid node includes: Querying the velocity of the grid point where the source point is located, and looking up the distance and the inclination angle of the path from the source point to the surrounding grid nodes in a table; Using the following anisotropic velocity formula to calculate the velocities of each surrounding grid node and then calculating the arrival time of the first arrival wave: V ani = V0(1 + δ sin 2 θ cos 2 θ + ε sin 4 θ); Where: V0 is the seismic wave propagation velocity in the direction of the symmetry axis of the vertical transverse isotropic VTI medium; δ and ε are the Thomas weak anisotropy parameters, and θ is the inclination angle of the path from the source point to the surrounding grid nodes.

5. The method according to claim 1, characterized in that Performing the enqueue and dequeue operations on the grid nodes in the queue includes: Dequeue the grid nodes in the queue, and use each grid node as a new source point to calculate the travel time t to adjacent grid nodes n , assuming the first arrival time of the current grid node is t, then the first arrival time propagated from the current grid node to the surrounding grid nodes is t + t n , if the grid node already has a first arrival time t n0 , and t n0 < t + t n , then do nothing, otherwise update the first arrival time of the grid node in the table and enqueue the grid node 6. The method according to any one of claims 1 to 5, characterized in that, The method is applied to oil and gas exploration.

7. An anisotropic in-well seismic ray shortest path tracking device, characterized in that Including: A data preparation module for data preparation, including the following sub-steps: S11. Reading the velocity model, calculating the velocity gradient in the velocity model, and picking up the maximum and minimum values of the velocity gradient; S12. Dividing the velocity model into equally spaced uniform grids based on the maximum value of the velocity gradient that can be tolerated at the lowest level within the same grid; S13. Initializing the observation system and normalizing the coordinates of the observation system. If the grid length is m, any point (x, y) on the grid is normalized to (x / m, y / m); S14. Initializing grid nodes: Establishing uniformly spaced grid nodes on all sides of the grid, numbering the grid nodes, and the information stored in each grid node includes which side of the current grid the corresponding grid node is on, the arrival time of the first arrival wave, and the previous grid node where the first arrival wave arrives; A single grid processing module for processing a single grid obtained by dissection, including the following sub-steps: S21. Constructing a distance table for single grid nodes: Calculating the distances from all grid nodes in the grid to any other grid node in the same grid to obtain a mapping table; S22. Constructing an inclination angle table for the ray paths between two grid nodes within the grid: Calculating the inclination angles of the ray paths between any two grid nodes and storing them in a two-dimensional data table; S23. Based on the inclination angle table, constructing two trigonometric function tables for different ray inclination angles; S24. Determining the effective propagation range of the wavelet source; An anisotropic ray tracing module for performing anisotropic ray tracing and calculating the arrival time of the first arrival wave, including the following sub-steps: S31. Starting from the source point, querying the velocity of the grid point where the source point is located, and looking up the distance and the inclination angle of the path from the source point to the surrounding grid nodes in a table, and calculating the velocities of each surrounding grid node; S32. Performing anisotropic ray tracing on the constructed distance table for single grid nodes, the constructed inclination angle table for the ray paths between two grid nodes within the grid, the constructed trigonometric function tables for different ray inclination angles, and the rule for determining the effective propagation range of the wavelet source, and storing the first arrival information and the path information in a new table; S33. Using a linked list to implement the storage structure, establishing a queue to store grid nodes, and using this queue for the breadth-first algorithm to traverse all grid nodes, and performing the enqueue operation on the grid nodes updated in the previous step; S34. According to the breadth-first search method, perform dequeue and enqueue operations on the grid nodes in the queue, and use each grid node as a new source point to calculate the travel time t to adjacent grid nodes n , assuming that the first arrival time of the current grid node is t, then the first arrival time propagated from the current grid node to the surrounding grid nodes is t + t n , if the grid node already has a first arrival time t n0 , and t n0 < t + t n , then do nothing, otherwise update the first arrival time of the grid node in the table and enqueue the grid node; S35. Repeat the steps for the grid nodes in the queue until the queue is empty. At this time, the table records the seismic wave propagation times and paths from the seismic source point to all grid nodes.

8. The device according to claim 7, wherein The single-grid processing module is further configured to calculate the distances from all grid nodes in the grid to any other grid node in the same grid by using the following method: Denote two grid nodes in the grid as the incident point and the exit point of the seismic wave ray respectively. Obtain the distance from the incident point of the seismic wave ray to the exit point of the seismic wave ray by using the following method: where len in,out represents the distance from the incident point to the exit point of the seismic wave ray, and x in , y in are the coordinate positions of the incident point of the seismic wave ray; x out , y out are the coordinate positions of the exit point of the seismic wave ray.

9. The device according to claim 7, characterized in that, The single-grid processing module is further configured to determine the effective propagation range of the sub-wave source by using the following method: Using Snell's law and the velocities of the adjacent grids of the grid nodes, calculate the direction θ2 of the seismic wave propagation path in the next grid with the ray incident angle θ1 of the path in the previous grid of the grid node. Limit the inclination angle range of θ2±Δθ as the effective propagation area of the sub-wave source, where Δθ is the set angular deviation.

10. The device according to claim 7, characterized in that, The anisotropic ray tracing module is further configured to calculate the velocities of the surrounding grid nodes by using the following method: Query the velocity of the grid point where the seismic source point is located, and look up the table to obtain the distances from the seismic source point to the surrounding grid nodes and the inclination angles of the paths. Calculate the velocities of the surrounding grid nodes by using the following anisotropic velocity formula, and then calculate the arrival time of the first arrival wave: V ani = V0(1 + δ sin 2 θ cos 2 θ + ε sin 4 θ); Where: V0 is the seismic wave propagation velocity in the direction of the symmetry axis of the vertical transverse isotropic (VTI) medium; δ and ε are the weak anisotropic parameters of Thomas, and θ is the inclination angle of the path from the seismic source point to the surrounding grid nodes.

11. An electronic device, characterized in that, Comprising: A processor; A memory for storing instructions executable by the processor; Wherein, the processor is configured to execute the instructions to implement the anisotropic borehole seismic ray shortest path tracing method according to any one of claims 1 to 6.

12. A computer-readable storage medium, characterized in that, When the instructions in the computer-readable storage medium are executed by the processor of the electronic device, the electronic device can execute the anisotropic borehole seismic ray shortest path tracing method according to any one of claims 1 to 6.

Images