A method for quickly determining seismic wave travel time

Abstract

The application discloses a method for quickly determining seismic wave travel time, which comprises the following steps: step one, reading in relevant parameter files and a velocity model; and step two, using a computer to judge all grid points to obtain seismic wave travel time of all grid points. The method has simple steps and reasonable design, and the seismic wave travel time is obtained by multiplying the distance from a grid point to a seismic source point and a seismic wave travel time disturbance, so that the problem of large seismic source error in the calculation of seismic wave travel time by the original fast marching method is solved.

Description

A method for rapid determination of seismic wave travel time Technical Field

[0001] This invention belongs to the field of seismic wave travel time technology, specifically relating to a method for rapid determination of seismic wave travel time. Background Technology

[0002] Seismic wave travel time calculations are mainly divided into three categories: first, ray tracing methods based on kinematic equations, but these methods are prone to getting trapped in local solutions and have shadow zone problems when encountering complex media; second, the shortest path algorithm based on Fermat's principle, namely the SPM algorithm, which approximates seismic wave rays by using a large number of grid segments; to ensure accuracy, this method requires a much larger number of network nodes than other methods, which also leads to a larger computational space requirement for the SPM algorithm; and third, methods based on solving equations using finite difference calculus. These methods effectively avoid the shortcomings of kinematic equation methods and have significant advantages in terms of computational accuracy and efficiency.

[0003] The Fast Marching Method (FMM) is a typical method for solving equations using the finite difference approach. It calculates the global seismic wave travel time by employing a narrow-band approximate wavefront extension. However, due to the large wavefront curvature near the epicenter, significant errors occur in this area, and these errors propagate throughout the computational domain. Currently, source-point mesh refinement techniques are used. Refining the mesh near the epicenter in FMM can significantly reduce the computational error near the source, but the error near the epicenter remains large relative to the entire computational domain. Therefore, mesh refinement in the epicenter region cannot fundamentally solve the problem of large errors at the epicenter. Furthermore, mesh refinement in the epicenter region also increases the computational load of the algorithm. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a method for rapidly determining the travel time of seismic waves, which addresses the shortcomings of the prior art. The method is simple in steps and reasonable in design. It uses a fast travel method to obtain the travel time disturbance of seismic waves, and obtains the travel time of seismic waves by multiplying the distance from the grid point to the source point with the travel time disturbance. This solves the problem of large source error in the original fast travel method for calculating the travel time of seismic waves.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is: a method for rapid determination of seismic wave travel time, characterized in that the method includes the following steps:

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

[0007] Step 2: Use a computer to determine the seismic wave travel times of all grid points:

[0008] Step 201: Use a computer to divide all grid point attributes into far-away points, completed points, and narrow-band points to obtain the far-away point set, the completed point set, and the narrow-band point set; where, initially, the completed point set W is an empty set, the narrow-band point set Z is an empty set, and the far-away point set Y is all grid points;

[0009] Step 202: Use a computer to set the travel time of the seismic source point to 0 and add it to the narrowband point set Z. Set the travel time of all other grid points to infinity.

[0010] Step 203: Find the grid point with the smallest travel time in the narrowband point set Z and update it as the completion point. Then the source point is updated as the first completion point and removed from the narrowband point set Z and added to the completion point set W. Use a computer to remove the six grid points around the first completion point whose attributes are far away from the point from the far away point set Y, and then add them to the narrowband point set Z.

[0011] In the narrowband point set Z, the computer calculates the travel time of the six grid points around the first completed point, compares the calculation result with the travel time of the six grid points, and selects the smaller value for each grid point to assign to the travel time of that grid point, thus obtaining the updated travel time of the grid point;

[0012] Step 204: Based on the update time of the grid points, select the grid point with the smallest time in the narrowband point set Z and mark it as the second completion point, remove it from the narrowband point set Z, and add it to the completion point set W;

[0013] The computer determines whether grid points with the attribute of being farthest from the second completed point should be added to the narrowband point set Z. If they have already been added to the narrowband point set Z, no further processing is performed; if they have not been added, they are removed from the farthest point set Y and then added to the narrowband point set Z.

[0014] Step 205: In the narrowband point set Z, the computer calculates the travel time of I2 grid points around the second completion point, compares the calculation result with the travel time of I2 grid points, selects the smaller value for each grid point and assigns it to the travel time of that grid point to obtain the updated travel time of the grid point; where I2 is a positive integer not greater than 6.

[0015] Step 206: Repeat steps 204 and 205 multiple times. Based on the update time of the grid points, select the grid point with the smallest time in the narrowband point set Z and mark it as the (j-1)th completion point. Remove it from the narrowband point set Z and add it to the completion point set W.

[0016] The computer determines whether the grid points around the (j-1)th completed point whose attributes are farthest from the point are added to the narrowband point set Z. If they are already added to the narrowband point set Z, no further processing is performed; if they are not added, they are removed from the farthest point set Y and then added to the narrowband point set Z. Here, j is a positive integer and 3≤j.

[0017] In the narrowband point set Z, a computer is used to analyze the I around the (j-1)th completion point. j-1 The travel time of each grid point is calculated, and the calculation result is compared with that of the I j-1 The travel times of each grid point are compared, and the smaller value is selected for each grid point and assigned to its travel time to obtain the updated travel time of that grid point; where I j-1 It is a positive integer not greater than 6;

[0018] Step 207: Repeat step 206 multiple times until the narrowband point set is empty;

[0019] In step 206, a computer is used to analyze the I values ​​around the (j-1)th completion point. j-1 The process of calculating the travel time of each grid point is the same. Specifically, the computer is used to calculate the travel time of the i-th narrowband point surrounding the (j-1)-th completion point, where i and I... j-1 All are positive integers, and 1≤i≤I j-1 I j-1 The total number of grid points surrounding the (j-1)th completion point is given by the following process:

[0020] Step 2061: Set the travel time of the i-th narrowband point using a computer as follows:

[0021] T(x,y,z)=L0(x,y,z)×T1(x,y,z)(1); where T(x,y,z) represents the travel time of the i-th narrowband point (x,y,z), L0(x,y,z) represents the distance from the i-th narrowband point (x,y,z) to the source point; T1(x,y,z) represents the travel time perturbation of the i-th narrowband point (x,y,z);

[0022] Step 2062: Using a computer, obtain the equation of the i-th narrowband point (x,y,z) according to equation (1), as follows: Where S(x,y,z) represents the slowness of the seismic wave propagation at point (x,y,z). Denotes the gradient operator, || 2 This represents the square of the magnitude of the gradient;

[0023] Step 2063: Use a computer to solve the equation of the i-th narrowband point (x,y,z) using the fast travel method to obtain the seismic wave travel time disturbance T1(x,y,z) of the i-th narrowband point (x,y,z).

[0024] Step 2064: Using a computer, substitute the X-axis coordinates (x), Y-axis coordinates (y), and Z-axis coordinates (z) of the i-th narrowband point (x, y, z) in the spatial rectangular coordinate system into the formula. We obtain the distance L0(x,y,z) from the i-th narrowband point (x,y,z) to the source point; where x0 represents the X-axis coordinate of the source point in the spatial rectangular coordinate system, y0 represents the Y-axis coordinate of the source point in the spatial rectangular coordinate system, and z0 represents the Z-axis coordinate of the source point in the spatial rectangular coordinate system.

[0025] Step 2065: Use a computer to substitute T1(x,y,z) obtained in step 2063 and L0(x,y,z) obtained in step 2064 into equation (1) in step 2061 to obtain the travel time T(x,y,z) of the i-th narrowband point (x,y,z).

[0026] The above-mentioned method for rapidly determining seismic wave travel time is characterized in that: in step 2063, a computer is used to solve the equation of the i-th narrowband point (x,y,z) using the fast travel method to obtain the seismic wave travel time disturbance T1(x,y,z) from the source point to the i-th narrowband point (x,y,z). The specific process is as follows:

[0027] The equation of the function at the i-th narrowband point (x,y,z) is discretized as follows:

[0028] in, This represents the nth-order difference at point (x,y,z) in the negative x-axis direction. This represents the nth-order difference at point (x,y,z) in the positive x-axis direction. This represents the nth-order difference at point (x,y,z) in the negative y-axis direction. This represents the nth-order difference at point (x, y, z) in the positive y-axis direction. This represents the nth-order difference at point (x,y,z) in the negative z-axis direction. This represents the nth-order difference at point (x, y, z) in the positive z-axis direction.

[0029] The above-mentioned method for rapid determination of seismic wave travel time is characterized by:

[0030] When n takes the value of 1, the first-order difference is as follows:

[0031]

[0032]

[0033]

[0034]

[0035]

[0036] Where h represents the grid spacing, T1(xh,y,z) represents the travel time perturbation at grid point coordinates (xh,y,z), and T1(x+h,y,z) represents the travel time perturbation at grid point coordinates (x+h,y,z). T1(x,yh,z) represents the travel time perturbation at grid point coordinates (x,yh,z), and T1(x,y+h,z) represents the travel time perturbation at grid point coordinates (x,y+h,z). T1(x,y,zh) represents the travel time perturbation at grid point coordinates (x,y,zh), and T1(x,y,z+h) represents the travel time perturbation at grid point coordinates (x,y,z+h). This represents partial derivative operations;

[0037] When n takes the value of 2, the second-order difference is as follows:

[0038]

[0039]

[0040]

[0041]

[0042]

[0043] Where T1(x-2h,y,z) represents the travel time perturbation at grid point coordinates (x-2h,y,z), and T1(x+2h,y,z) represents the travel time perturbation at grid point coordinates (x+2h,y,z). T1(x,y-2h,z) represents the travel time perturbation at grid point coordinates (x,y-2h,z), and T1(x,y+2h,z) represents the travel time perturbation at grid point coordinates (x,y+2h,z). T1(x,y,z-2h) represents the travel time perturbation at grid point coordinates (x,y,z-2h), and T1(x,y,z+2h) represents the travel time perturbation at grid point coordinates (x,y,z+2h). This represents the partial derivative operation.

[0044] Compared with the prior art, the present invention has the following advantages:

[0045] 1. The method of the present invention has simple steps and reasonable design, and fundamentally solves the problem of large source error when calculating the travel time of seismic waves using the rapid travel method.

[0046] 2. This invention obtains the seismic wave travel time by multiplying the distance from the grid point to the source point by the seismic wave travel time perturbation. By employing a fast travel algorithm to obtain the seismic wave travel time perturbation of the i-th narrowband point, the seismic wave travel time of the i-th narrowband point is obtained, thereby improving the accuracy and efficiency of seismic wave travel time calculation.

[0047] In summary, the method of the present invention is simple in steps and reasonable in design. It obtains the seismic wave travel time by multiplying the distance from the grid point to the source point with the seismic wave travel time disturbance, thus solving the problem of large source error in the original fast travel method for calculating seismic wave travel time.

[0048] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0049] Figure 1 is a flowchart of the method of the present invention.

[0050] Figure 2 shows a 3D model of the spatial region of seismic wave propagation according to the present invention.

[0051] Figure 3 shows the error distribution of the first-order difference calculation method used in this invention.

[0052] Figure 4 shows the error distribution of the second-order difference calculation used in this invention. Detailed Implementation

[0053] Example 1

[0054] As shown in Figure 1, the present invention provides a method for rapid determination of seismic wave travel time, comprising the following steps:

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

[0056] Step 2: Use a computer to determine the seismic wave travel times of all grid points:

[0057] Step 201: Use a computer to divide all grid point attributes into far-away points, completed points, and narrow-band points to obtain the far-away point set, the completed point set, and the narrow-band point set; where, initially, the completed point set W is an empty set, the narrow-band point set Z is an empty set, and the far-away point set Y is all grid points;

[0058] Step 202: Use a computer to set the travel time of the seismic source point to 0 and add it to the narrowband point set Z. Set the travel time of all other grid points to infinity.

[0059] Step 203: Find the grid point with the smallest travel time in the narrowband point set Z and update it as the completion point. Then the source point is updated as the first completion point and removed from the narrowband point set Z and added to the completion point set W. Use a computer to remove the six grid points around the first completion point whose attributes are far away from the point from the far away point set Y, and then add them to the narrowband point set Z.

[0060] In the narrowband point set Z, the computer calculates the travel time of the six grid points around the first completed point, compares the calculation result with the travel time of the six grid points, and selects the smaller value for each grid point to assign to the travel time of that grid point, thus obtaining the updated travel time of the grid point;

[0061] Step 204: Based on the update time of the grid points, select the grid point with the smallest time in the narrowband point set Z and mark it as the second completion point, remove it from the narrowband point set Z, and add it to the completion point set W;

[0062] The computer determines whether grid points with the attribute of being farthest from the second completed point should be added to the narrowband point set Z. If they have already been added to the narrowband point set Z, no further processing is performed; if they have not been added, they are removed from the farthest point set Y and then added to the narrowband point set Z.

[0063] Step 205: In the narrowband point set Z, the computer calculates the travel time of I2 grid points around the second completion point, compares the calculation result with the travel time of I2 grid points, selects the smaller value for each grid point and assigns it to the travel time of that grid point to obtain the updated travel time of the grid point; where I2 is a positive integer not greater than 6.

[0064] Step 206: Repeat steps 204 and 205 multiple times. Based on the update time of the grid points, select the grid point with the smallest time in the narrowband point set Z and mark it as the (j-1)th completion point. Remove it from the narrowband point set Z and add it to the completion point set W.

[0065] The computer determines whether the grid points around the (j-1)th completed point whose attributes are farthest from the point are added to the narrowband point set Z. If they are already added to the narrowband point set Z, no further processing is performed; if they are not added, they are removed from the farthest point set Y and then added to the narrowband point set Z. Here, j is a positive integer and 3≤j.

[0066] In the narrowband point set Z, a computer is used to analyze the I around the (j-1)th completion point. j-1 The travel time of each grid point is calculated, and the calculation result is compared with that of the I j-1 The travel times of each grid point are compared, and the smaller value is selected for each grid point and assigned to its travel time to obtain the updated travel time of that grid point; where I j-1 It is a positive integer not greater than 6;

[0067] Step 207: Repeat step 206 multiple times until the narrowband point set is empty;

[0068] In step 206, a computer is used to analyze the I values ​​around the (j-1)th completion point. j-1 The process of calculating the travel time of each grid point is the same. Specifically, the computer is used to calculate the travel time of the i-th narrowband point surrounding the (j-1)-th completion point, where i and I... j-1 All are positive integers, and 1≤i≤I j-1 I j-1 The total number of grid points surrounding the (j-1)th completion point is given by the following process:

[0069] Step 2061: Set the travel time of the i-th narrowband point using a computer as follows:

[0070] T(x,y,z)=L0(x,y,z)×T1(x,y,z)(1); where T(x,y,z) represents the travel time of the i-th narrowband point (x,y,z), L0(x,y,z) represents the distance from the i-th narrowband point (x,y,z) to the source point; T1(x,y,z) represents the travel time perturbation of the i-th narrowband point (x,y,z);

[0071] Step 2062: Using a computer, obtain the equation of the i-th narrowband point (x,y,z) according to equation (1), as follows: (2); where S(x,y,z) represents the slowness of the seismic wave propagation at point (x,y,z). Denotes the gradient operator, || 2 This represents the square of the magnitude of the gradient;

[0072] Step 2063: Use a computer to solve the equation of the i-th narrowband point (x,y,z) using the fast travel method to obtain the seismic wave travel time disturbance T1(x,y,z) of the i-th narrowband point (x,y,z).

[0073] Step 2064: Using a computer, substitute the X-axis coordinates (x), Y-axis coordinates (y), and Z-axis coordinates (z) of the i-th narrowband point (x, y, z) in the spatial rectangular coordinate system into the formula. We obtain the distance L0(x,y,z) from the i-th narrowband point (x,y,z) to the source point; where x0 represents the X-axis coordinate of the source point in the spatial rectangular coordinate system, y0 represents the Y-axis coordinate of the source point in the spatial rectangular coordinate system, and z0 represents the Z-axis coordinate of the source point in the spatial rectangular coordinate system.

[0074] Step 2065: Use a computer to substitute T1(x,y,z) obtained in step 2063 and L0(x,y,z) obtained in step 2064 into equation (1) in step 2061 to obtain the travel time T(x,y,z) of the i-th narrowband point (x,y,z).

[0075] In this embodiment, step 2063 uses a computer to solve the equation of the i-th narrowband point (x,y,z) using the fast travel method, to obtain the seismic wave travel time disturbance T1(x,y,z) from the source point to the i-th narrowband point (x,y,z). The specific process is as follows:

[0076] The equation of the function at the i-th narrowband point (x,y,z) is discretized as follows:

[0077] in, This represents the nth-order difference at point (x,y,z) in the negative x-axis direction. This represents the nth-order difference at point (x,y,z) in the positive x-axis direction. This represents the nth-order difference at point (x,y,z) in the negative y-axis direction. This represents the nth-order difference at point (x, y, z) in the positive y-axis direction. This represents the nth-order difference at point (x,y,z) in the negative z-axis direction. This represents the nth-order difference at point (x, y, z) in the positive z-axis direction.

[0078] In this embodiment, when n is 1, the first-order difference is as follows:

[0079]

[0080]

[0081]

[0082]

[0083]

[0084] Where h represents the grid spacing, T1(xh,y,z) represents the travel time perturbation at grid point coordinates (xh,y,z), and T1(x+h,y,z) represents the travel time perturbation at grid point coordinates (x+h,y,z). T1(x,yh,z) represents the travel time perturbation at grid point coordinates (x,yh,z), and T1(x,y+h,z) represents the travel time perturbation at grid point coordinates (x,y+h,z). T1(x,y,zh) represents the travel time perturbation at grid point coordinates (x,y,zh), and T1(x,y,z+h) represents the travel time perturbation at grid point coordinates (x,y,z+h). This represents partial derivative operations;

[0085] When n takes the value of 2, the second-order difference is as follows:

[0086]

[0087]

[0088]

[0089]

[0090]

[0091] Where T1(x-2h,y,z) represents the travel time perturbation at grid point coordinates (x-2h,y,z), and T1(x+2h,y,z) represents the travel time perturbation at grid point coordinates (x+2h,y,z). T1(x,y-2h,z) represents the travel time perturbation at grid point coordinates (x,y-2h,z), and T1(x,y+2h,z) represents the travel time perturbation at grid point coordinates (x,y+2h,z). T1(x,y,z-2h) represents the travel time perturbation at grid point coordinates (x,y,z-2h), and T1(x,y,z+2h) represents the travel time perturbation at grid point coordinates (x,y,z+2h). This represents the partial derivative operation.

[0092] In this embodiment, first-order difference and second-order difference are used for comparative calculation. It should be noted that when using second-order difference, it is actually a mixture of first-order and second-order difference. When the narrowband point meets the calculation conditions of second-order difference, second-order difference is used for calculation. When the narrowband point does not meet the calculation conditions of second-order difference, first-order difference is used for calculation.

[0093] In this embodiment, in step one, relevant parameter files and velocity models are read in to obtain a 3D model of the spatial region for seismic wave propagation.

[0094] In this embodiment, the velocity model represents the propagation velocity of seismic waves at each grid point.

[0095] As shown in Figure 2, in this embodiment, a spatial rectangular coordinate system is established as follows: the upper left corner of the surface of the 3D model of the seismic wave propagation space region is taken as the origin, the horizontal line of the ground surface passing through the origin is the X-axis, the line pointing downwards through the origin is the Z-axis, and the line pointing outwards perpendicular to the X-axis and Z-axis through the origin is the Y-axis. Thus, the surface of the 3D model of the seismic wave propagation space region is the ground surface.

[0096] In this embodiment, the 3D model of the seismic wave propagation spatial region is a 3D uniform model, that is: the size of the 3D model of the seismic wave propagation spatial region is 10km×10km×10km, the seismic wave propagation speed is 2km / s, and the coordinates of the source point (x0,y0,z0) are (5km,5km,5km).

[0097] In this embodiment, the original fast travel algorithm uses a source point grid densification method to ensure accuracy, that is, the grid is densified in a spatial region of 1km near the source point, and the spacing of the densified grid is 0.005km.

[0098] In this embodiment, two indicators, mean absolute error and mean relative error, are used to evaluate the calculation accuracy. The formula for calculating the mean absolute error is: The formula for calculating the average relative error is: Where T exact (i) represents the actual travel time of the seismic wave at the i-th grid point, and T exact (i) is the ratio of the distance from the i-th grid point to the source point to the seismic wave propagation velocity at the i-th grid point; T cal (i) represents the seismic wave travel time of the i-th grid point obtained by the original fast travel algorithm or the method of this invention, and N is the total number of grid points, as shown in Table 1 below.

[0099] Table 1. Comparison of errors and CPU time in calculating seismic wave travel time between the original fast travel algorithm and the method of this invention in 3D models.

[0100]

[0101] In this embodiment, as shown in Table 1, for the uniform model, the calculation results of the method of the present invention are completely accurate (the error analysis results have eliminated the calculation errors caused by the computer's own decimal point retention). Furthermore, as the grid spacing decreases, the computational efficiency of both the original fast traversal algorithm and the method of the present invention decreases; however, at the same grid spacing, the computational efficiency of the method of the present invention is significantly higher than that of the original fast traversal algorithm.

[0102] Example 2

[0103] In this embodiment, the difference from Embodiment 1 is that the 3D model of the seismic wave propagation spatial region is a 3D velocity linearly increasing model, that is, the size of the 3D model of the seismic wave propagation spatial region is 10km×10km×10km, the seismic wave propagation velocity increases linearly from 2km / s at the surface (z=0) to 4km / s at the bottom interface (z=10km), and the coordinates of the source point (x0,y0,z0) are (5km,5km,0km).

[0104] In this embodiment, the formula for calculating the actual seismic wave travel time of the i-th grid point is as follows: Where α is the gradient of the seismic wave propagation velocity, v0 is the seismic wave propagation velocity at the source point, and v i Let L be the seismic wave propagation velocity at the i-th grid point. i,0 represents the distance of the i-th grid point from the source point, and arccosh(·) is the inverse hyperbolic cosine function.

[0105] In this embodiment, α is 0.2 and v0 is 2km / s.

[0106] In this embodiment, since the seismic source point of this model is located on the ground, the depth-direction densification area in the original fast travel algorithm is (z0, z0+1km).

[0107] In this embodiment, Table 2 is obtained as follows:

[0108] Table 2. Comparison of errors and CPU time in calculating seismic wave travel time between the original fast travel algorithm and the method of this invention in 3D models.

[0109]

[0110] In this embodiment, as shown in Table 2, for the 3D linearly increasing velocity model, as the grid spacing decreases, both the average absolute error and the average relative error of the original fast traversal algorithm and the method of this invention continuously decrease, resulting in improved computational accuracy. However, the decrease in grid spacing leads to an increase in the number of grid points, which in turn reduces computational efficiency. Nevertheless, when the grid spacing is the same, the method of this invention outperforms the original fast traversal algorithm in both computational accuracy and CPU computation time.

[0111] In this embodiment, the formula for calculating the absolute error is: |T exact (i)-T cal (i)|, the formula for calculating the relative error is: |T exact (i)-T cal (i)| / T exact (i).

[0112] In this embodiment, Figure 3 shows the error distribution diagram calculated using the first-order difference: Figure 3(a),(c),(e),(g) are the error distribution diagrams of each grid point on the 5km x-axis profile, and Figure 3(b),(d),(f),(h) are the results of the 5km y-axis profile. Among them, Figure 3(a),(b),(e),(f) are absolute error distribution diagrams; Figures (c),(d),(g),(h) are relative error distribution diagrams; Figures (a),(b),(c),(d) are obtained using the original fast travel algorithm, and Figures (e),(f),(g),(h) are obtained using the method of this invention.

[0113] Figure 4 shows the error distribution diagram of the second-order difference calculation: (a), (c), (e), and (g) in Figure 4 are the error distribution diagrams of each grid point of the 5km profile on the x-axis, and (b), (d), (f), and (h) are the results of the 5km profile on the y-axis. Among them, (a), (b), (e), and (f) in Figure 4 are the absolute error distribution diagrams; (c), (d), (g), and (h) are the relative error distribution diagrams. Figures (a), (b), (c), and (d) are obtained using the original fast travel algorithm, and Figures (e), (f), (g), and (h) are obtained using the method of this invention. It can be found that the seismic wave travel time calculated by the original fast travel algorithm has a large error in the region near the diagonal of the epicenter. When the first-order difference is used, the maximum absolute error can reach 20ms and the maximum relative error is 4%, while the maximum absolute error obtained by the method of this invention is only 1.5ms and the maximum relative error is only 0.04%. When using second-order finite difference, the computational accuracy of both algorithms is significantly improved. The maximum absolute error of the original fast marching algorithm is reduced to 1.5 ms, and the maximum percentage error is reduced to 2%. The method of this invention further reduces the maximum absolute error to 0.1 ms and the maximum percentage error to 0.02%. It is evident that the method of this invention, regardless of whether first-order or second-order finite difference is used, has higher computational accuracy than the original fast marching algorithm. Furthermore, the error distribution diagram shows that the method of this invention solves the problem of large errors at the seismic source point.

[0114] In summary, the method of the present invention has simple steps and reasonable design, and solves the problem of large source error when the original rapid travel method calculates the travel time of seismic waves.

[0115] The above description is merely a preferred embodiment of the present invention and does not constitute any limitation on the present invention. Any simple modifications, alterations, or equivalent structural changes made to the above embodiments based on the technical essence of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for rapid determination of seismic wave travel time, characterized in that, The method includes the following steps: Step 1: Read in the relevant parameter file and velocity model; wherein, the parameter file contains the number of grid points, grid spacing, and source location of the velocity model; Step 2: Use a computer to determine the seismic wave travel time of all grid points: Step 201: Use a computer to divide the attributes of all grid points into far-away points, completed points, and narrowband points, obtaining the far-away point set, completed point set, and narrowband point set; wherein, initially, the completed point set W is an empty set, the narrowband point set Z is an empty set, and the far-away point set Y is all grid points; Step 202: Use a computer to set the travel time of the source point to 0 and add it to the narrowband point set Z, and set the travel time of all other grid points to infinity; Step 2 03. Find the grid point with the smallest travel time in the narrowband point set Z and update it as the completion point. Then, the source point is updated as the first completion point, and it is removed from the narrowband point set Z and added to the completion point set W. Using a computer, remove the six grid points around the first completion point whose attributes are far away from the point from the far away point set Y, and then add them to the narrowband point set Z. In the narrowband point set Z, use a computer to calculate the travel time of the six grid points around the first completion point, compare the calculation result with the travel time of the six grid points, and select the smaller value for each grid point to assign to its travel time, thus obtaining the updated travel time of the grid point. Step 204. Select the narrowband point based on the updated travel time of the grid points. The grid point with the smallest travel time in set Z is designated as the second completion point and removed from the narrowband point set Z, then added to the completion point set W. A computer is used to determine whether grid points with the attribute of being farthest from the second completion point (up, down, left, right, front, back, center) should be added to the narrowband point set Z. If they are already added, no further processing is done; otherwise, they are removed from the farthest point set Y and then added to the narrowband point set Z. Step 205: In the narrowband point set Z, a computer calculates the travel time of I2 grid points surrounding the second completion point. The calculated result is compared with the travel time of these I2 grid points. For each grid point, the smaller value is selected and assigned to its travel time, thus obtaining the updated travel time of the grid point. Time; where I2 is a positive integer not greater than 6; Step 206: Repeat steps 204 and 205 multiple times. Based on the update time of the grid points, select the grid point with the smallest time in the narrowband point set Z and record it as the (j-1)th completion point. Remove it from the narrowband point set Z and add it to the completion point set W; Use a computer to determine whether the grid points around the (j-1)th completion point with the attribute of being far away from the point have been added to the narrowband point set Z. If they have been added to the narrowband point set Z, no processing is performed; if they have not been added, remove them from the far away point set Y and then add them to the narrowband point set Z; where j is a positive integer and 3≤j; In the narrowband point set Z, use a computer to analyze the I2 values ​​around the (j-1)th completion point. j-1 The travel time of each grid point is calculated, and the calculation result is compared with that of the I j-1 The travel times of each grid point are compared, and the smaller value is selected for each grid point and assigned to its travel time to obtain the updated travel time of that grid point; where I j-1 The integer is a positive integer not greater than 6; Step 207, repeat step 206 multiple times until the narrowband point set is empty; in step 206, the computer is used to analyze the I values ​​around the (j-1)th completed point. j-1 The process of calculating the travel time for each grid point is the same. Specifically, the computer is used to calculate the travel time of the i-th narrowband point surrounding the (j-1)-th completion point, where i and I... j-1 All are positive integers, and 1≤i≤I j-1 I j-1 The total number of grid points around the (j-1)th completion point is as follows: Step 2061: The computer sets the travel time of the i-th narrowband point as follows: T(x,y,z)=L0(x,y,z)×T1(x,y,z) (1); where T(x,y,z) represents the travel time of the i-th narrowband point (x,y,z), L0(x,y,z) represents the distance from the i-th narrowband point (x,y,z) to the source point; T1(x,y,z) represents the travel time disturbance of the i-th narrowband point (x,y,z); Step 2062: The computer obtains the equation of the i-th narrowband point (x,y,z) according to equation (1), as follows: (2); where S(x,y,z) represents the slowness of the seismic wave propagation at point (x,y,z). Denotes the gradient operator, | | 2 The square of the magnitude of the gradient is given; Step 2063: Using a computer, the equation of the i-th narrowband point (x,y,z) is solved using the fast travel method to obtain the seismic wave travel time disturbance T1(x,y,z) of the i-th narrowband point (x,y,z); Step 2064: Using a computer, the X-axis coordinate (x), Y-axis coordinate (y), and Z-axis coordinate (z) of the i-th narrowband point (x,y,z) in the spatial rectangular coordinate system are substituted into the formula. Obtain the distance L0(x,y,z) from the i-th narrowband point (x,y,z) to the source point; where x0 represents the X-axis coordinate of the source point in the spatial rectangular coordinate system, y0 represents the Y-axis coordinate of the source point in the spatial rectangular coordinate system, and z0 represents the Z-axis coordinate of the source point in the spatial rectangular coordinate system; Step 2065: Use a computer to substitute T1(x,y,z) obtained in step 2063 and L0(x,y,z) obtained in step 2064 into equation (1) in step 2061 to obtain the travel time T(x,y,z) of the i-th narrowband point (x,y,z).

2. The method for rapid determination of seismic wave travel time according to claim 1, characterized in that: In step 2063, a computer is used to solve the equation of the i-th narrowband point (x,y,z) using the fast travel method, to obtain the seismic wave travel time disturbance T1(x,y,z) from the source point to the i-th narrowband point (x,y,z). The specific process is as follows: The equation of the i-th narrowband point (x,y,z) is discretized, as follows: in, This represents the nth-order difference at point (x,y,z) in the negative x-axis direction. This represents the nth-order difference at point (x,y,z) in the positive x-axis direction. This represents the nth-order difference at point (x,y,z) in the negative y-axis direction. This represents the nth-order difference at point (x, y, z) in the positive y-axis direction. This represents the nth-order difference at point (x,y,z) in the negative z-axis direction. This represents the nth-order difference at point (x, y, z) in the positive z-axis direction.

3. The method for rapid determination of seismic wave travel time according to claim 2, characterized in that: When n takes the value of 1, the first-order difference is as follows: Where h represents the grid spacing, T1(xh,y,z) represents the travel time perturbation at grid point coordinates (xh,y,z), and T1(x+h,y,z) represents the travel time perturbation at grid point coordinates (x+h,y,z). T1(x,yh,z) represents the travel time perturbation at grid point coordinates (x,yh,z), and T1(x,y+h,z) represents the travel time perturbation at grid point coordinates (x,y+h,z). T1(x,y,zh) represents the travel time perturbation at grid point coordinates (x,y,zh), and T1(x,y,z+h) represents the travel time perturbation at grid point coordinates (x,y,z+h). This represents the partial derivative operation; when n is 2, the second-order difference is as follows: Where T1(x-2h,y,z) represents the travel time perturbation at grid point coordinates (x-2h,y,z), and T1(x+2h,y,z) represents the travel time perturbation at grid point coordinates (x+2h,y,z). T1(x,y-2h,z) represents the travel time perturbation at grid point coordinates (x,y-2h,z), and T1(x,y+2h,z) represents the travel time perturbation at grid point coordinates (x,y+2h,z). T1(x,y,z-2h) represents the travel time perturbation at grid point coordinates (x,y,z-2h), and T1(x,y,z+2h) represents the travel time perturbation at grid point coordinates (x,y,z+2h). This represents partial derivative operations.

Images