An eikonal equation seismic wave travel time calculation method, system, computer storage medium and terminal

Abstract

The application belongs to the technical field of seismic wave travel time calculation, and relates to a method, a system, a computer storage medium and a terminal for calculating the travel time of a seismic wave based on an eikonal equation, the method comprising the following steps: obtaining a velocity model and relevant parameter information of a calculation area; initializing the calculation area, setting grid point attributes and initial source information; discretely solving the eikonal equation by using a high-order upwind difference scheme based on a weighted essential non-oscillatory format, and iteratively calculating the travel time value of a seismic wave at each grid point by using a narrow-band extension method; and outputting a calculation result when a termination condition is met. The application combines the weighted essential non-oscillatory format with the framework of the fast marching method, thereby keeping the unconditional stability of the algorithm, significantly improving the calculation accuracy of the travel time of a seismic wave under complex geological conditions, making the required calculation template more compact, reducing the calculation amount, improving the calculation efficiency, and providing reliable support for fine processing of tunnel seismic exploration data.

Description

Technical Field

[0001] This invention relates to the field of seismic wave travel time calculation technology, and in particular to a method, system, computer storage medium, and terminal for calculating seismic wave travel time using a process function equation. Background Technology

[0002] During tunnel construction, complex geological conditions, including well-developed faults and folds, frequently lead to sudden geological disasters such as collapses and water inrushes, severely impacting construction progress and even endangering the safety of workers. Tunnel Seismic Prediction (TSP) can effectively detect adverse geological conditions ahead of the tunnel face, guiding tunnel excavation and preventing construction accidents. Seismic wave travel time is a crucial parameter in seismic exploration data processing, and its calculation accuracy significantly affects the accuracy of subsequent data processing and seismic prediction.

[0003] In seismic wave travel time calculation methods, the equation-based method effectively avoids the problems of caustics and shadow areas inherent in traditional ray-based methods, making it the mainstream approach. The classic method within the equation-based seismic wave travel time algorithm is the fast-progression method, which uses an upwind difference scheme to solve the discretized equations and employs a narrowband extension to calculate seismic wave travel time. However, for complex tunnel conditions, the conventional fast-progression method has limitations in computational accuracy and cannot meet the requirements of practical, high-precision exploration. Higher-order schemes are a direct way to improve the computational accuracy of travel time algorithms. Compared to traditional higher-order finite difference schemes using Taylor expansion, the Weighted Essentially Non-Oscillatory (WENO) scheme offers higher accuracy and efficiency.

[0004] The journal *Advances in Exploration Geophysics*, Vol. 30, No. 5, 2007, published an article by Sun Zhangqing et al. entitled "Calculation of Seismic Wave Travel Time on Rippling Surfaces Using the Wavefront Rapid Advancement Method." The article details the specific implementation scheme of the conventional rapid advance method, explains the unconditional stability of the algorithm, and proposes a rapid advance method for rippled surfaces in actual data processing. The feasibility of the algorithm is verified through a theoretical model.

[0005] Chinese invention patent application 202010433784.2 discloses "A Method for Calculating Seismic Wave Travel Time." This invention first determines the basic form of the fundamental travel time calculation formula. Then, based on the Taylor series principle and the continuous fraction approximation formula for the square of travel time, it performs a higher-order approximation on the fundamental travel time calculation formula. Subsequently, it simultaneously reduces the powers of the fourth and sixth powers of the offset in the higher-order fundamental travel time calculation formula to output the seismic wave travel time. This algorithm improves the accuracy of seismic wave travel time calculation through higher-order calculations, but the multiple calculations increase the computation time.

[0006] Chinese invention patent application 202311556673.0 discloses "A method for calculating the travel time of seismic waves for tunnel exploration". This invention is based on the conventional rapid propagation method and improves the calculation efficiency of seismic wave travel time by optimizing the narrowband extension method and updating a set of points at a time. However, this invention uses a first-order difference operator when solving the equation, so there is still room for improvement in its calculation accuracy.

[0007] Therefore, how to provide a method that can effectively improve the accuracy of seismic wave travel time calculation for complex geological conditions (especially tunnel conditions) while ensuring computational efficiency is an urgent problem to be solved by those skilled in the art. Summary of the Invention

[0008] The purpose of this invention is to overcome the shortcomings of existing technologies in calculating seismic wave travel time with high accuracy and efficiency in complex tunnel conditions, and to provide a method, system, computer storage medium and terminal for calculating seismic wave travel time using the equation of process function.

[0009] In a first aspect, the present invention provides a method for calculating the travel time of seismic waves using a process function equation, comprising the following steps:

[0010] Step 1: Obtain the velocity model and related parameter information of the seismic wave travel time calculation area. The parameter information includes at least the grid size, grid spacing and source location. Step 2: Initialize the computational domain by setting grid point attributes and initial seismic source information; Step 3: Using a high-order upwind differential scheme based on a weighted essential non-oscillatory scheme, the equation is discretized and solved. The seismic wave travel time values ​​of each grid point are iteratively calculated using a narrowband extension method. Step 4: When the preset termination conditions are met, output the seismic wave travel time calculation results for all grid points within the calculation area.

[0011] In step 1 of the technical solution of this invention, the velocity model is a data volume that digitally describes the spatial distribution of seismic wave propagation velocity in the subsurface medium. It is typically represented using a gridded method, with each grid point assigned a corresponding velocity value or its reciprocal, i.e., a slowness value. The grid size refers to the number of grid points in the horizontal and vertical directions of the calculation area; the grid spacing refers to the actual distance between adjacent grid points; and the source location refers to the coordinates of the starting point that triggered the seismic wave. These parameters together constitute the basic data framework for calculating seismic wave travel time.

[0012] In step 2, the grid point attributes are categorized into three types: accepting points, narrowband points, and distant points. Accepting points are those whose travel times have already been calculated and will not participate in subsequent updates. Narrowband points are all grid points located on the approximate wavefront; these points have not yet been calculated and will still participate in subsequent updates. Distant points are all grid points except accepting points and narrowband points; these points have not yet participated in any updates. Source initialization specifically includes: setting all grid point attributes to distant points and travel times to infinity; setting the source point attribute to accepting points and travel times to zero; setting the attributes of adjacent grid points to narrowband points to form the initial narrowband, and calculating their travel times using analytical formulas. This initialization process lays the data foundation for subsequent wavefront extension.

[0013] In step 3, the weighted intrinsically oscillatory scheme is a high-precision numerical discretization method that approximates the derivative term by weighted combination of multiple candidate templates. The weights are adaptively determined based on the local smoothness of each template, achieving high-precision in smooth regions and automatically avoiding spurious oscillations in non-smooth regions. Specifically, based on the third-order weighted intrinsically oscillatory scheme, forward and backward difference schemes for the travel time gradient term in the equation are constructed. These schemes are then merged into the upwind difference scheme of the fast-margining method, resulting in the solution equation for the fast-margining method based on the third-order weighted intrinsically oscillatory scheme. The update calculation formula for the seismic wave travel time values ​​at the grid points is determined based on the solution equation.

[0014] Preferably, the equation is solved discretely, specifically including: Step 31: Construct the forward and backward difference schemes for the travel time gradient term in the equation based on the third-order weighted essential non-oscillatory scheme; Step 32: Combine the forward difference scheme and backward difference scheme of the travel time gradient term into the upwind difference scheme of the fast advance method to obtain the solution equation of the fast advance method based on the third-order weighted essential non-oscillating scheme; Step 33: Determine the update calculation formula for the seismic wave travel time values ​​at the grid points based on the solved equations.

[0015] The above process is the core technical aspect of this invention. By organically combining the weighted, inherently oscillatory scheme with the fast-progression method framework, it maintains the advantage of unconditional stability of the fast-progression method while significantly improving computational accuracy through higher-order schemes. Compared with traditional higher-order Taylor expansion schemes, the weighted, inherently oscillatory scheme requires fewer computation points under the same accuracy conditions, thereby effectively reducing the computational load and improving computational efficiency.

[0016] Preferably, step 31 specifically includes: At grid points in the computational domain, weighted interpolation is performed using the seismic wave travel time values ​​of three adjacent grid points to construct a difference expression for the travel time gradient term; wherein, the weights of the weighted interpolation are adaptively determined based on the local smoothness of each interpolation template.

[0017] This is the essential characteristic that distinguishes the weighted, oscillatory-free scheme from other higher-order schemes. The local interpolation template of the three grid points gives the scheme a compact template width, which is beneficial for boundary handling and parallel computation; the adaptive weighting mechanism ensures that the scheme will not produce non-physical numerical oscillations in discontinuous or large gradient regions, which is of great significance for dealing with complex geological structures such as faults and fracture zones.

[0018] Preferably, in step 2, the computational domain is initialized by setting grid point attributes and initial seismic source information. Specific steps include: Step 21: Initialize the properties of all grid points within the calculation area to be far away from the point, and set the seismic wave travel time value to infinity; Step 22: Set the source point's attributes to receive point and set the seismic wave travel time value to zero; Step 23: Set the grid points adjacent to the source point to narrowband points to form an initial narrowband, and calculate the seismic wave travel time value according to the analytical formula.

[0019] This initialization scheme strictly follows the basic framework of the fast-progression method, ensuring that the algorithm starts from the source point and gradually progresses outwards in a wavefront expansion manner. The accurate calculation of the initial narrowband provides a reliable starting value for subsequent iterative updates and is the foundation for the stable operation of the entire calculation process.

[0020] Preferably, step 3, which iteratively calculates the seismic wave travel time values ​​of each grid point using a narrowband extension method, specifically includes: Select the grid point with the smallest seismic wave travel time value from the current narrow band as the minimum travel time point, and change the attribute of this point to the receiving point; Perform attribute determination on the neighboring grid points of the minimum travel time point: If a neighboring point is a receiving point, no processing is performed; If the neighboring point is a narrow band point, the seismic wave travel time value is recalculated according to the updated calculation formula; If a neighboring point is a distant point, first modify its attributes to a narrowband point, move it into the narrowband, and then calculate the seismic wave travel time value according to the updated calculation formula.

[0021] This narrowband extension process is the core iterative mechanism of the rapid propagation method. By consistently selecting the point with the shortest travel time from the narrowband for extension, it ensures that the wavefront progresses according to causality; that is, points with smaller travel times are calculated first, and points with larger travel times are calculated later. This extension method conforms to the physical laws of seismic wave propagation, ensuring the rationality of the calculation results.

[0022] Preferably, the method further includes a boundary condition processing step: after each narrowband extension calculation is completed, a preset boundary condition is applied to the boundary grid points of the calculation region, and the grid point attribute for performing the boundary condition calculation is set as the accepting point, ensuring that all points in the calculation region have completed the update calculation.

[0023] Boundary conditions are a crucial component of boundary value problems in partial differential equations and must be considered in numerical calculations. For seismic wave travel time calculations, the treatment of boundary grid points directly affects the accuracy of the solution across the entire computational domain. This step ensures that seismic wave travel times can be accurately calculated even near the boundaries of the computational domain, preventing error accumulation or propagation due to improper boundary treatment.

[0024] Preferably, the preset termination condition is: all grid points within the calculation area have the attribute of receiving points; when this termination condition is met, the iterative calculation ends and the final seismic wave travel time calculation result is output.

[0025] This termination condition is based on the fundamental theory of the rapid advance method: when all grid points become receiving points, it means that the wavefront has extended to the entire computational domain, and the travel time value of each grid point has been calculated according to the equation and the given velocity model. At this point, the output travel time field is the final seismic wave travel time calculation result that satisfies the equation.

[0026] Secondly, the present invention provides a system for calculating the travel time of seismic waves using a functional equation, comprising: The model parameter input module is used to read in the velocity model and related parameter information, wherein the parameter information includes at least the grid size, grid spacing and earthquake source location; The initialization module is used to declare attributes and initialize seismic sources for grid points within the computational domain. The narrowband extended calculation module is used to select the minimum travel time point from the narrowband and perform attribute judgment and travel time calculation on the neighboring points of that point; The boundary condition processing module is used to execute boundary conditions after each narrowband extension calculation; The termination judgment module is used to check whether the attributes of all grid points in the calculation area are all receiving points. If so, the calculation ends and the final seismic wave travel time calculation results are output. The narrowband extended computing module uses a high-order upwind differential scheme based on a third-order weighted essentially oscillatory scheme for timekeeping calculation.

[0027] Thirdly, the present invention provides a computer storage medium storing a computer program, which, when executed by a processor, implements the steps of the above-described method.

[0028] Fourthly, the present invention provides a computer terminal, including a processor, a memory, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the above-described calculation method.

[0029] Compared with the prior art, the beneficial effects of the present invention are as follows: In the technical solution of this invention, a unified data benchmark is established for seismic wave travel time calculation by acquiring velocity models and related parameter information, ensuring that subsequent numerical calculations can accurately correspond to actual geographic coordinates. By initializing the calculation area, setting grid point attributes and initial source information, and strictly adhering to the basic framework of the fast-progression method, the algorithm is ensured to gradually advance outward from the source point in a manner consistent with physical laws. A high-order upwind difference scheme based on a weighted intrinsically non-oscillating scheme is used to discretize and solve the equations, organically combining the weighted intrinsically non-oscillating scheme with the fast-progression method framework. This maintains the advantage of unconditional stability of the fast-progression method while significantly improving the accuracy of seismic wave travel time calculation through the high-order scheme. Simultaneously, the weighted intrinsically non-oscillating scheme requires a more compact calculation template under the same accuracy conditions, effectively reducing the computational load for updating each grid point, improving overall computational efficiency, and solving the technical challenge of balancing high accuracy and high efficiency in complex tunnel conditions. Attached Figure Description

[0030] Figure 1 This is a flowchart of the method for calculating the travel time of seismic waves using the Cheng-Jun equation.

[0031] Figure 2 This is a schematic diagram of the local discrete point distribution of the third-order WENO scheme of the present invention, wherein (a) is a schematic diagram of the template of the forward difference scheme in the x direction, and (b) is a schematic diagram of the template of the backward difference scheme in the x direction.

[0032] Figure 3 This is a schematic diagram of the calculation process of the algorithm of the present invention.

[0033] Figure 4 This is a comparison chart of the absolute errors of the uniform model.

[0034] Figure 5This is a contour map of seismic wave travel time for a complex model. Detailed Implementation

[0035] The present invention will now be described in further detail with reference to specific embodiments. However, this should not be construed as limiting the scope of the present invention to the following embodiments; all technologies implemented based on the content of the present invention fall within the scope of the present invention.

[0036] Example 1 This embodiment provides a method for calculating the travel time of seismic waves using a process function equation, which specifically includes the following steps: Step 1: Obtain the velocity model and related parameter information of the seismic wave travel time calculation area.

[0037] A velocity model is a data volume that digitally describes the spatial distribution of seismic wave propagation velocities in the subsurface medium. It is typically represented using a gridded approach, with each grid point assigned a corresponding velocity value or its reciprocal, the slowness value. The grid size refers to the number of grid points in the calculation area both horizontally and vertically; the grid spacing refers to the actual distance between adjacent grid points; and the source location refers to the coordinates of the originating point that triggered the seismic wave. These parameters collectively constitute the basic data framework for calculating seismic wave travel times.

[0038] This embodiment uses a uniform model for illustration. The model size is 800×800, the grid spacing is 1m×1m, the source point is located at (400,400), and the velocity is 2000m / s.

[0039] Step 2: Initialize the computational domain by setting grid point attributes and initial seismic source information.

[0040] All grid points within the computational region are categorized into three types: accepting points, narrowband points, and distant points. Accepting points are those whose travel times have already been calculated and will not participate in subsequent updates. Narrowband points are all grid points located on the approximate wavefront; their calculations are not yet complete and they will still participate in subsequent updates. Distant points are all grid points except accepting points and narrowband points; these points have not yet participated in any updates. During the calculation, 't' represents the travel time value, and 'ch' represents the grid point attribute.

[0041] Initialization specifically includes the following sub-steps: Step 21: Initialize the properties of all grid points in the computational region to be far away points, and set their travel time value to infinity (1000 in this embodiment).

[0042] Step 22: Set the source point's properties to receive point and set its travel time value to zero.

[0043] Step 23: Set the grid points (up, down, left, right) adjacent to the source point as narrowband points to form the initial narrowband, and calculate its travel time value according to the analytical formula.

[0044] In this embodiment, the epicenter is located at (400, 400), and the attributes of its adjacent points (399, 400), (401, 400), (400, 399), and (400, 401) are modified to be narrowband points. The travel time values ​​are calculated using analytical formulas.

[0045] in, This represents the analytical time value at that point. Indicates grid spacing. This indicates the slowness at that point.

[0046] In this embodiment, the epicenter is located at (400, 400). First, the travel time of the epicenter is set to 0, and its grid point attribute is set to receiving point. Then, the travel times of the remaining points are set to the maximum value, which is 1000 in this embodiment, and their grid point attributes are set to far-away point. Finally, the attributes of the four grid points (399, 400), (401, 400), (400, 399), and (400, 401) near the epicenter are modified to narrowband points, and an initial narrowband is constructed based on these four points. The travel time value is calculated using an analytical formula.

[0047]

[0048] Step 3: Using a high-order upwind differential scheme based on a weighted essential oscillation-free scheme, the equation is discretized and solved. The seismic wave travel time values ​​of each grid point are iteratively calculated using a narrowband extension method.

[0049] This step is the core of the present invention, and specifically includes the following sub-steps: Step 31: Construct a time-lapse gradient term difference scheme based on the third-order WENO scheme. Equation of the process function in a two-dimensional isotropic medium: (1) in, Indicating departure, Let the slowness be represented by the reciprocal of the velocity. Using Godunov's upwind finite difference discretization, we can obtain the following: (2) in, , Indicates the grid spacing, and (3) The travel-time gradient term in the equation is represented using the third-order WENO scheme. The local interpolation template for the third-order WENO scheme is shown in Figure 2. Figure 2In the figure, (a) is a template diagram of the forward difference scheme in the x direction, and (b) is a template diagram of the backward difference scheme in the x direction. The grid points in the figure are used to illustrate the positional relationship of the adjacent grid points on which the calculation of the travel gradient of the current point depends.

[0050] Grid points within the calculation area The backward difference scheme for the travel time value in the x-direction is as follows: (4) in: (5) Grid points within the calculation area The forward difference scheme for the travel time value in the x-direction is as follows: (6) in: (7) in, This represents intermediate variables in the calculation process. It is a constant, and the calculation is set to the minimum value.

[0051] Similarly, grid points within the computational region Backward and forward difference schemes for travel time values ​​in the z-direction A similar method can also be used to obtain it.

[0052] Step 32: Incorporate the travel time gradient terms in the WENO format into the upwind difference scheme of the fast propulsion method to obtain the solution equation.

[0053] The higher-order travel gradient terms based on the WENO scheme in formulas (4)-(7) are merged into the scheme of the fast propulsion method upwind difference scheme: (8) (9) Based on the above, we can obtain the solution equations using the fast-progression method based on the third-order WENO solution scheme: (10) in: (11) in, Represents grid points Place The solution to be updated numerical values. Indicates the same grid point The current unupdated value.

[0054] Step 33: Determine the update calculation formula for seismic wave travel time values ​​at grid points.

[0055] For the point to be updated in the narrow band, its travel time T new It can be obtained by solving equation (10).

[0056] The solution to equation (10) is: (12) Where T old The values ​​are currently not updated, h is the grid spacing, and sij is the slowness at that point. In actual calculations, the appropriate difference value needs to be selected based on the windward direction.

[0057] Step 34: Iteratively calculate the travel time value of each grid point using the narrowband extension method.

[0058] The narrowband extension process is as follows: Select the grid point with the smallest travel time value from the current narrow band as the minimum travel time point P, and change the attribute of this point to accept point.

[0059] Perform attribute determination on the four neighboring points of the point with the shortest travel time P: If a neighboring point is a receiving point, no processing is performed; If the neighboring point is a narrowband point, then its travel time value is recalculated according to the update calculation formula in step 33; If a neighboring point is a distant point, first change its attribute to a narrowband point (move into the narrowband), and then calculate its travel time value according to the updated calculation formula.

[0060] Repeat the above process until all grid points have been calculated.

[0061] Figure 3 The black dots represent receiving points (travel time has been calculated), the gray dots represent narrowband points (points on the current wavefront, currently being calculated), and the white dots represent distant points (not yet being calculated). The distribution of dots of different colors illustrates the evolution of narrowband expansion.

[0062] Step 4: Boundary conditions and termination criteria.

[0063] Step 41: Boundary condition processing. After each narrowband extension calculation is completed, preset boundary conditions are applied to the boundary grid points of the computational domain, and the grid point attributes for boundary condition calculation are set to accept points to ensure that all points in the computational domain have completed the update calculation.

[0064] Specifically, narrowband extension calculation: Minimum travel time selection: Select the minimum travel time from the narrow band. And change the attribute of the minimum elapsed time point to the accept point (ch=2); Narrowband extension evolution: for minimum travel time The attributes of the four neighboring points are judged. If it is an accepting point, no processing is performed. If it is a narrow band point, its travel time is updated and calculated according to the solution formula (12). If it is a far away point, its attributes are first modified to narrow band point (moved into narrow band), and then its travel time value is calculated according to the solution formula (12).

[0065] In this embodiment, the boundary conditions can be expressed as: In the two-dimensional case, when using the fast-progression method of the third-order WENO solution scheme to calculate seismic wave travel time, boundary conditions must be considered. Assume the computational domain is: , For a point in the calculation region, The boundary conditions it satisfies are: (13) Step 42: Termination Check. Check if all grid points within the calculation area are receiver points. If so, end the iterative calculation and output the final seismic wave travel time calculation results; otherwise, return to Step 3 to continue the calculation.

[0066] Calculation result verification This embodiment uses a uniform model for verification, with model parameters as described above. The seismic wave travel time analysis values ​​are: Its seismic wave travel time resolution is:

[0067] in, This represents the analytical solution. Indicates slowness. Indicates the focal point, Indicates the grid spacing.

[0068] Comparison of seismic wave travel time and analytical values ​​calculated using the method of this invention, for example Figure 4 As shown, the results indicate that the calculation results of this invention are closer to the analytical values, and the travel time contour distribution conforms to the seismic wave propagation law. Compared with the conventional rapid advance method, the method of this invention significantly improves the calculation accuracy. (a) shows the absolute travel time error distribution obtained using the conventional rapid advance method, and (b) shows the absolute travel time error distribution obtained using the algorithm of this invention. The comparison demonstrates that the algorithm of this invention has smaller errors and higher accuracy.

[0069] For complex models such as Figure 5 As shown in the figure, the calculation results of the algorithm of the present invention under complex geological conditions are presented. The contour lines in the figure are smooth and continuous, and the distribution pattern conforms to the physical characteristics of seismic wave propagation, indicating that the algorithm has good adaptability to complex models.

[0070] The seismic wave travel time calculation method provided by this invention, for a uniform model, shows that the calculated results are closer to the analytical values, and the travel time contour distribution conforms to the seismic wave propagation law. Compared with the calculation results of conventional algorithms, the calculation results of this invention improve the accuracy of seismic wave travel time calculation. For complex models, the seismic wave travel time contour map in the calculation results of this invention conforms to the seismic wave propagation law, indicating that the algorithm has good adaptability to complex models.

[0071] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for calculating seismic wave travel time, characterized in that, Includes the following steps: Step 1: Obtain the velocity model and related parameter information of the seismic wave travel time calculation area. The parameter information includes at least the grid size, grid spacing and source location. Step 2: Initialize the computational domain by setting grid point attributes and initial seismic source information; Step 3: Using a high-order upwind differential scheme based on a weighted essential non-oscillatory scheme, the equation is discretized and solved. The seismic wave travel time values ​​of each grid point are iteratively calculated using a narrowband extension method. Step 4: When the preset termination conditions are met, output the seismic wave travel time calculation results for all grid points within the calculation area.

2. In the seismic wave travel time calculation method according to claim 1, step 3, discretizing the equation, specifically includes: Step 31: Construct the forward and backward difference schemes for the travel time gradient term in the equation based on the third-order weighted essential non-oscillatory scheme; Step 32: Combine the forward difference scheme and backward difference scheme of the travel time gradient term into the upwind difference scheme of the fast advance method to obtain the solution equation of the fast advance method based on the third-order weighted essential non-oscillating scheme; Step 33: Determine the update calculation formula for the seismic wave travel time values ​​at the grid points based on the solved equations.

3. The seismic wave travel time calculation method according to claim 2, step 31 specifically includes: At grid points in the computational domain, weighted interpolation is performed using the seismic wave travel time values ​​of three adjacent grid points to construct a difference expression for the travel time gradient term; wherein, the weights of the weighted interpolation are adaptively determined based on the local smoothness of each interpolation template.

4. The seismic wave travel time calculation method according to claim 1, characterized in that, Step 2 involves initializing the computational domain by setting grid point attributes and initial seismic source information. Specific steps include: Step 21: Initialize the properties of all grid points within the calculation area to be far away from the point, and set the seismic wave travel time value to infinity; Step 22: Set the source point's attributes to receive point and set the seismic wave travel time value to zero; Step 23: Set the grid points adjacent to the source point to narrowband points to form an initial narrowband, and calculate the seismic wave travel time value according to the analytical formula.

5. The seismic wave travel time calculation method according to claim 3, characterized in that, Step 3, which iteratively calculates the seismic wave travel time values ​​for each grid point using a narrowband extension method, specifically includes: Select the grid point with the smallest seismic wave travel time value from the current narrow band as the minimum travel time point, and change the attribute of this point to the receiving point; Perform attribute determination on the neighboring grid points of the minimum travel time point: If a neighboring point is a receiving point, no processing is performed; If the neighboring point is a narrow band point, the seismic wave travel time value is recalculated according to the updated calculation formula; If a neighboring point is a distant point, first modify its attributes to a narrowband point, move it into the narrowband, and then calculate the seismic wave travel time value according to the updated calculation formula.

6. The seismic wave travel time calculation method according to claim 1, characterized in that, The method also includes a boundary condition processing step: After each narrowband extension calculation is completed, preset boundary conditions are applied to the boundary grid points of the computational domain, and the grid point attributes for which boundary conditions are calculated are set to accept points, ensuring that all points in the computational domain are updated.

7. The seismic wave travel time calculation method according to claim 1, characterized in that, The preset termination condition is: all grid points within the calculation area are classified as receiving points; when this termination condition is met, the iterative calculation ends and the final seismic wave travel time calculation results are output.

8. A seismic wave travel time calculation system, characterized in that, include: The model parameter input module is used to read in the velocity model and related parameter information, wherein the parameter information includes at least the grid size, grid spacing and earthquake source location; The initialization module is used to declare attributes and initialize seismic sources for grid points within the computational domain. The narrowband extended calculation module is used to select the minimum travel time point from the narrowband and perform attribute judgment and travel time calculation on the neighboring points of that point; The boundary condition processing module is used to execute boundary conditions after each narrowband extension calculation; The termination judgment module is used to check whether the attributes of all grid points in the calculation area are all receiving points. If so, the calculation ends and the final seismic wave travel time calculation results are output. The narrowband extended computing module uses a high-order upwind differential scheme based on a third-order weighted essentially oscillatory scheme for timekeeping calculation.

9. A computer storage medium storing a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1-7.

10. A computer terminal, comprising a processor, a memory, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method as described in any one of claims 1-7.

Images