Rapid scanning seismic wave travel time calculation method based on process function equation

Abstract

The invention discloses a fast scanning seismic travel time calculation method based on a process function equation, and relates to the field of seismic travel time calculation, and the method comprises the steps: firstly, reading a speed model and related parameter information; then carrying out seismic source initialization, setting a slowness fixed item, calculating a travel time fixed item, and setting a travel time disturbance item according to the position of a seismic source point; next, Gauss-Seidel iterative calculation is carried out, in four alternating directions, a Gauss-Seidel type iterative method is used for solving an eikonal equation of factorization, a proper travel time disturbance term is selected through a causal condition, and finally, a travel time calculation result meeting a termination condition is output. According to the high-precision seismic travel time calculation method based on the factorization eikonal equation, the problem of seismic source singularity is well solved, high precision is guaranteed, meanwhile, the algorithm has high calculation efficiency, and the accuracy of subsequent seismic exploration data processing is improved.

Description

Technical Field

[0001] The present invention relates to the field of seismic wave travel time calculation, and particularly to a fast-scanning seismic wave travel time calculation method based on an improved eikonal equation. Background Art

[0002] The seismic wave travel time is an important parameter describing the kinematic characteristics of seismic waves, and is widely used in fields such as prestack migration, travel time tomography, and earthquake location. The calculation accuracy of the seismic wave travel time determines the accuracy of the processing results in the above fields. The seismic wave travel time solution method based on the eikonal equation is inevitably affected by the singularity of the seismic source. The source error will reach the entire calculation area along with the propagation of the wavefront, thus affecting the solution accuracy of the entire calculation area. In addition, for the processing of large model seismic data, the calculation efficiency of the travel time algorithm is also very important. Therefore, it is of great significance to study an efficient and high-precision seismic wave travel time calculation method for the seismic source singularity problem.

[0003] Chinese Patent Application No. 201810077621.8 discloses "A Hybrid Two-Dimensional Seismic Wave Travel Time Calculation Method". This algorithm reads in the velocity model and related calculation parameters, traces the ray information along different directions from the shot point for a certain distance, calculates the seismic wave travel time within the ray range using the wavefront construction method, and then calculates the travel time of the remaining grid nodes using the fast marching method. The good calculation effect of this method is proved by numerical simulation.

[0004] Chinese Patent Application No. 202010433784.2 discloses "A Seismic Wave Travel Time Calculation Method". This algorithm performs high-order approximation on the basic travel time calculation formula based on the Taylor series principle and the continuous fraction approximation time square formula, and then performs power reduction processing on the high-order travel time calculation formula to obtain the seismic wave travel time, and verifies the calculation effect of the algorithm through numerical simulation.

[0005] Chinese Patent Application No. 202410096376.0 discloses "A Seismic Wave Travel Time Solving Method and System Based on the Fast Marching Method of Equivalent Slowness". This invention provides a seismic wave travel time solving method and system based on the fast marching method of equivalent slowness. This method makes the space of the target grid formed by moving the discrete grid coincide with the coordinates of the seismic source point, establishes the relationship between the travel time of a certain point, the straight-line distance from the point to the seismic source point, and the equivalent slowness, and converts the equation based on the travel time into an equation based on the equivalent slowness according to this relationship. The equivalent slowness has a smaller radius of curvature, so that during the stage of calculating and solving the travel time of the equation, the calculated travel time has a higher accuracy. Finally, the travel time on the initial discrete grid points is calculated correspondingly, avoiding the systematic error caused by the incomplete coincidence of the seismic source point with the discrete grid, thereby improving the calculation accuracy of the seismic wave travel time.

[0006] As can be seen from the above examples, the conventional seismic wave travel-time algorithm has great limitations in terms of calculation accuracy and efficiency. The eikonal equation-based algorithm is a commonly used calculation method in travel-time calculation, but the source singularity will lead to large errors at the source point and will spread to the entire calculation area. Therefore, how to solve the source singularity and study an efficient and high-precision seismic wave travel-time calculation method is of great significance. Summary of the Invention

[0007] The technical problem to be solved by the present invention is to provide a fast-scan seismic wave travel-time calculation method based on an improved eikonal equation. The travel-time term and slowness term in the eikonal equation are decomposed into fixed terms and perturbation terms by using the multiplicative decomposition form. The fixed terms propagate in a spherical form and are known parameters, and the perturbation terms are correction factors and are parameters to be solved. The algorithm approximates the wavefront by spherical waves, well solves the source singularity problem in the travel-time calculation of complex models, uses the fast-scan method to solve the factored eikonal equation to obtain the seismic wave travel-time, and gives a new solution form of the local discrete equation of the seismic wave travel-time. Compared with the traditional fast-scan method, under the condition of similar calculation efficiency, the source error is effectively reduced and the accuracy of the seismic wave travel-time calculation is improved.

[0008] To solve the above technical problem, the technical solution adopted by the present invention is: a fast-scan seismic wave travel-time calculation method based on an improved eikonal equation, including the following steps:

[0009] Step 1: Determine the seismic wave travel-time calculation formula; the seismic wave travel-time calculation formula is an expression obtained by solving the factored eikonal equation using the fast-scan method.

[0010] Step 1.1 Eikonal equation in two-dimensional isotropic medium:

[0011]

[0012] Where: T represents the travel-time, and S represents the slowness, which is the reciprocal of the velocity.

[0013] To solve the source singularity problem, by performing multiplicative decomposition on Equation (1), the travel-time term and slowness term in the eikonal equation are decomposed into fixed terms and perturbation terms. Among them, the fixed terms propagate in a spherical form and are known parameters, and the perturbation terms are correction factors and are parameters to be solved. The algorithm approximates the wavefront by spherical waves, well solves the source singularity problem.

[0014] T(x,z) = T0(x,z)τ(x,z) (2)

[0015] S(x,z) = S0(x,z)α(x,z) (3)

[0016] Assume:

[0017]

[0018] Among them, \(T_0(x, z)\) represents the fixed travel time term, \(\tau(x, z)\) represents the travel time perturbation term, \(S_0(x, z)\) represents the fixed slowness term, and \(\alpha(x, z)\) represents the slowness perturbation term.

[0019] In formulas (2) and (3), both \(T_0(x, z)\) and \(S_0(x, z)\) are known. Using formulas (2), (3), and (4), the seismic wave travel time \(T(x, z)\) in equation (1) is converted into a solution problem in the eikonal equation \(\tau(x, z)\) to be decomposed.

[0020]

[0021] Step 1.2 Determine the solution formula; formula (5) is the seismic wave travel time calculation formula, and the following takes a two-dimensional rectangular grid as an example to illustrate the specific solution formula of the fast-scan seismic wave travel time calculation method based on the improved eikonal equation.

[0022] Figure 2 It shows an internal grid point C, and its adjacent points are W, E, N, and S. Discretize equation (5) on the four triangles: \(\Delta CEN\), \(\Delta CNW\), \(\Delta CWS\), and \(\Delta CSE\). For example, on the triangle \(\Delta CWS\), a discretized equation can be obtained:

[0023]

[0024] According to the quadratic formula for finding the roots of a quadratic equation, a new form of the seismic wave travel time solution for the discretized equation (6) can be obtained. Let \(\Delta = B\) 2 \(- 4AC\), and when \(\Delta\geq0\), solve for \(\tau\) C It is expressed as:

[0025]

[0026] Among them, \(\tau\) C \(= x\), \(\tau\) W \(= a\), \(\tau\) S \(= b\), \(h\) represents the grid spacing, \(T_0(x, z)\) represents the fixed travel time term, \(\tau(x, z)\) represents the travel time perturbation term, \(S_0(x, z)\) represents the fixed slowness term, \(\alpha(x, z)\) represents the slowness perturbation term, \(T\) 0x , \(T\) 0y respectively represent the gradients of \(T_0(x, z)\) in the \(x\) and \(z\) directions, and their mathematical expressions are:

[0027]

[0028] The solution \(\tau\) of formula (6) CThe causal condition needs to be satisfied, but Equation (6) may have no real roots or may have real roots that do not satisfy the causal condition. In this case, the characteristic method is used to transfer the information of τ along the edge from W to C and along the edge from S to C. By using the characteristic method on the edges respectively, two values of τ can be calculated: τ WC , τ SC . Assume that the edge (where δx = h, δy = 0, L = h)

[0029]

[0030] Step 2: Read in the velocity model and related parameter information, where the parameters include the grid size, grid spacing, and source position of the velocity model;

[0031] Step 3: Source initialization

[0032] Step 3.1 For the slowness fixed term S0(x, z), it is set to the velocity of a homogeneous medium, 1 m / s, then the travel time fixed term T0 is:

[0033]

[0034] For the source point (x0, y0), there is:

[0035] T0(x0, y0) = 0 (12)

[0036] Step 3.2 For the travel time perturbation term τ(x, z). Set τ(x0, y0) = 1 at the source point and assign a large positive value to all other grid points.

[0037] Step 4 Gauss - Seidel iterative calculation

[0038] In four alternating directions, use the Gauss - Seidel type iterative method to solve Equation (6), as Figure 3 shown. The scanning calculation directions are:

[0039] (1) i = 1:I; j = 1:J

[0040] (2) i = 1:I; j = J:1

[0041] (3) i = I:1; j = 1:J

[0042] (4) i = I:1; j = J:1

[0043] Step 4.1 As Figure 2As shown, at the grid point C, the discretized equations corresponding to the four triangles ΔCEN, ΔCNW, ΔCWS and ΔCSE composed of its neighboring points are solved respectively. For example, for the triangle ΔCWS, solve equation (6) to find two possible roots τ C,1 and τ C,2 .

[0044] If there are two real roots τ C,1 and τ C,2 , then first make a causal relationship judgment:

[0045] τ C,1 T0(C)≥τ W T0(W),τ C,1 T0(C)≥τ S T0(S)

[0046] τ C,2 T0(C)≥τ W T0(W),τ C,2 T0(C)≥τ S T0(S)

[0047] If τ C,1 and τ C,2 If both satisfy the causal relationship condition, then T WS =min{τ C,1 T0(C),τ C,2 T0(C)};

[0048] If τ C,1 If the causal relationship condition is met, then T WS =τ C,1 T0(C);

[0049] If τ C,2 If the causal relationship condition is met, then T WS =τ C,2 T0(C);

[0050] Otherwise, if both roots do not satisfy the causal condition, solve the characteristic equation according to formulas (9) and (10), using the edge The characteristic root τ is obtained by the characteristic method WC and τ SC .Execute causality τ WC T0(C)≥T W ,τ SC T0(C)≥T S .

[0051] If τ WC and τ SC If both satisfy the causal relationship, then T WS =min{τ WCT0(C), τ SC T0(C)};

[0052] If τ WC satisfies the causality, then T WS = τ WC T0(C);

[0053] If τ SC satisfies the causality, then T WS = τ SC T0(C);

[0054] If τ WC and τ SC both satisfy the causality, then T WS = ∞;

[0055] Step 4.2 In addition, when the equation (6) has no real roots, according to formulas (9) and (10), find the solutions of the characteristic equation, and use the method of characteristic sides to obtain τ WC and τ SC . Execute the causality τ WC T0(C) ≥ T W , τ SC T0(C) ≥ T S .

[0056] If τ WC and τ SC both satisfy the causality, then T WS = min{τ WC T0(C), τ SC T0(C)};

[0057] If τ WC satisfies the causality, then T WS = τ WC T0(C);

[0058] If τ SC satisfies the causality, then T WS = τ SC T0(C);

[0059] If τ WC and τ SC both satisfy the causality, then T WS = ∞;

[0060] Step 4.3 Select the minimum value from the results of the four triangle calculations, and calculate the value of τ(C) at the grid point C.

[0061]

[0062] Step 5: Termination condition

[0063] Set the threshold δ→0. If the travel time ∑|τ old -τ new | obtained in the previous iteration is less than δ, stop the iterative calculation and return the calculation result τ(x,z) of the perturbation factor. Otherwise, repeat step (4) until the termination condition is satisfied. Finally, the calculation result of the seismic wave travel time is T(x,z) = T0(x,z)τ(x,z).

[0064] The beneficial effects of the above technical solution are as follows: The present invention provides a fast-scanning seismic wave travel time calculation method based on an improved eikonal equation. By using the multiplicative decomposition form, the travel time term and the slowness term in the eikonal equation are decomposed into fixed terms and perturbation terms. The algorithm approximates the wavefront by spherical waves, which well solves the source singularity problem in the travel time calculation of complex models. The fast-scanning method is used to solve the decomposed eikonal equation to obtain the seismic wave travel time, and a new solution form of the local discrete equation of the seismic wave travel time is given. Compared with the traditional fast-scanning method, under the condition of similar calculation efficiency, the source singularity problem is solved, the source error is effectively reduced, and the accuracy of the seismic wave travel time calculation is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 A fast-scanning seismic wave travel time calculation method based on an improved eikonal equation according to the present invention

[0066] Figure 2 Schematic diagram of local solution of the algorithm of the present invention

[0067] Figure 3 Calculation scanning process of the algorithm of the present invention

[0068] Figure 4 Absolute error diagram of the gradient model in the embodiment of the present invention

[0069] Figure 5 Contour map of seismic wave travel time of complex model in the embodiment of the present invention DETAILED DESCRIPTION OF THE EMBODIMENTS

[0070] 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 the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0071] Please refer to Figures 1-5

[0072] The present invention provides a fast-scanning seismic wave travel time calculation method based on an improved eikonal equation, as Figure 1 shown, including the following steps:

[0073] Step 1: Determine the seismic wave travel time calculation formula; the seismic wave travel time calculation formula is an expression obtained by decomposing the eikonal equation and using the fast marching method to solve it.

[0074] Step 1.1 Eikonal equation in two-dimensional isotropic medium:

[0075] |▽T(x,z)| 2 = S 2 (x,z) (1)

[0076] Where: T represents the travel time, and S represents the slowness, which is the reciprocal of the velocity.

[0077] To solve the source singularity problem, by decomposing equation (1), the travel time term and the slowness term in the eikonal equation are decomposed into a fixed term and a perturbation term. Among them, the fixed term propagates in a spherical form and is a known parameter, and the perturbation term is a correction factor and is a parameter to be solved. The algorithm approximates the wavefront by a spherical wave and well solves the source singularity problem.

[0078] T(x,z) = T0(x,z)τ(x,z) (2)

[0079] S(x,z) = S0(x,z)α(x,z) (3)

[0080] Assume:

[0081]

[0082] Among them, T0(x,z) represents the travel time fixed term, τ(x,z) represents the travel time perturbation term, S0(x,z) represents the slowness fixed term, and α(x,z) represents the slowness perturbation term.

[0083] In formulas (2) and (3), both T0(x,z) and S0(x,z) are known. Using formulas (2), (3), and (4), the seismic wave travel time T(x,z) in equation (1) is converted into a solution problem in the factored eikonal equation τ(x,z).

[0084]

[0085] Step 1.2 Determine the solution formula; formula (5) is the seismic wave travel time calculation formula. The following takes a two-dimensional rectangular grid as an example to illustrate the specific solution formula of the fast marching seismic wave travel time calculation method based on the improved eikonal equation.

[0086] Figure 2Shows an internal grid point C, whose neighboring points are W, E, N, S. Discretize equation (5) on four triangles: ΔCEN, ΔCNW, ΔCWS, and ΔCSE. For example, on triangle ΔCWS, a discretized equation can be obtained:

[0087]

[0088] According to the quadratic formula for finding the roots of a quadratic equation, solve τ for equation (6) C Can be expressed as: Δ = B 2 - 4AC, and when Δ ≥ 0,

[0089]

[0090] where τ C = x, τ W = a, τ S = b, h represents the grid spacing, T0(x, z) represents the fixed travel time term, τ(x, z) represents the travel time perturbation term, S0(x, z) represents the fixed slowness term, α(x, z) represents the slowness perturbation term, T 0x , T 0y respectively represent the gradients of T0(x, z) in the x and z directions, and their mathematical expressions are:

[0091]

[0092] The solution τ of formula (6) C Needs to satisfy the causality condition, but formula (6) may have no real roots, or may have real roots that do not satisfy the causality condition. In this case, use the characteristic method to separately transmit the information of τ along the edge from W to C, and along the edge to transmit from S to C. By separately using the characteristic method on the edge , two values of τ can be calculated: τ WC , τ SC . Assume that the edge (where δx = h, δy = 0, L = h)

[0093]

[0094] Step 2: Read in the velocity model and related parameter information, where the parameters include the grid size, grid spacing, and source location of the velocity model;

[0095] In this embodiment, the gradient model is used to calculate the absolute error. The model size is 400 × 400, the grid spacing is 1 m × 1 m, the source point is located at (200, 0), and the velocity distribution of the gradient model and the corresponding seismic wave travel time analysis respectively satisfy (S0 = 0.001, G0 = 5):

[0096]

[0097] Among them, S(x) represents slowness, and T exact (x) represents the analytical solution, and

[0098] Step 3: Source initialization

[0099] Step 3.1 In this embodiment, for the slowness fixed term S0(x,z), it is set to the velocity of a homogeneous medium of 1 m / s, then the travel time fixed term T0 is:

[0100]

[0101] In this embodiment, for the source point (x0,y0), there is:

[0102] T0(x0,y0) = 0

[0103] Step 3.2 In this embodiment, for the travel time perturbation term τ(x,z). Set the source point τ(x0,y0) = 1, and assign a relatively large positive value (set to 10000) to all other grid points.

[0104] Step 4 Gauss - Seidel iterative calculation

[0105] In four alternating directions, use the Gauss - Seidel type iterative method to solve formula (6), as Figure 3 shown, the scanning calculation directions are:

[0106] (1) i = 1:I; j = 1:J

[0107] (2) i = 1:I; j = J:1

[0108] (3) i = I:1; j = 1:J

[0109] (4) i = I:1; j = J:1

[0110] Step 4.1 As Figure 2 shown, at grid point C, solve the discretized equations corresponding to the four triangles ΔCEN, ΔCNW, ΔCWS, and ΔCSE formed by its neighboring points respectively. For example, for triangle ΔCWS, solve equation (6) to find two possible roots τ C,1 and τ C,2 .

[0111] If there are two real roots τ C,1 and τ C,2 , then first perform a causality judgment:

[0112] τC,1 T0(C) ≥ τ W T0(W), τ C,1 T0(C) ≥ τ S T0(S)

[0113] τ C,2 T0(C) ≥ τ W T0(W), τ C,2 T0(C) ≥ τ S T0(S)

[0114] If τ C,1 and τ C,2 both satisfy the causality condition, then T WS = min{τ C,1 T0(C), τ C,2 T0(C)};

[0115] If τ C,1 satisfies the causality condition, then T WS = τ C,1 T0(C);

[0116] If τ C,2 satisfies the causality condition, then T WS = τ C,2 T0(C);

[0117] Otherwise, if neither root satisfies the causality condition, solve the characteristic equation according to formulas (9) and (10), and obtain the characteristic roots τ

[0118] and τ WC and τ SC . Execute the causality τ WC T0(C) ≥ T W , τ SC T0(C) ≥ T S .

[0119] If τ WC and τ SC both satisfy the causality, then T WS = min{τ WC T0(C), τ SC T0(C)};

[0120] If τ WC satisfies the causality, then T WS = τ WC T0(C);

[0121] If τ SC satisfies the causality, then T WS = τ SC T0(C);

[0122] If τ WC and τ SC both satisfy the causality, then T WS = ∞;

[0123] Step 4.2 In addition, when the equation (6) has no real roots, according to formulas (9) and (10), find the solutions of the characteristic equation, and use the method of characteristic sides to obtain τ WC and τ SC . Execute the causality τ WC T0(C) ≥ T W , τ SC T0(C) ≥ T S .

[0124] If τ WC and τ SC both satisfy the causality, then T WS = min{τ WC T0(C), τ SC T0(C)};

[0125] If τ WC satisfies the causality, then T WS = τ WC T0(C);

[0126] If τ SC satisfies the causality, then T WS = τ SC T0(C);

[0127] If τ WC and τ SC both satisfy the causality, then T WS = ∞;

[0128] Step 4.3 Select the minimum value from the results of the four triangle calculations, and calculate the value of τ(C) at the grid point C.

[0129]

[0130] Step 5: Termination condition

[0131] Set the threshold δ → 0. In this embodiment, set the termination threshold δ = 10 -9 . If the travel time ∑|τ old - τ new | < δ obtained in the previous cycle, stop the iterative calculation, and at the same time return the calculation result τ(x, z) of the perturbation factor. Otherwise, repeat Step 4 until the termination condition is satisfied. Finally, output the calculation result of the factorization eikonal equation fast sweeping method as T(x, z) = T0(x, z)τ(x, z).

[0132] In specific use, for a fast-scanning seismic wave travel time calculation method based on the improved eikonal equation of the present invention and for the gradient model, analyzing the calculation results shows that: the calculation results of the present invention are close to the analytical solution results, and the distribution of travel time isolines conforms to the seismic wave propagation law. Compared with the calculation results of the conventional fast-scanning method, the calculation results of the present invention well solve the source singularity problem, reduce the source error, and improve the calculation accuracy of seismic wave travel time. For the complex model, in the calculation results of the present invention, the seismic wave travel time isoline map conforms to the seismic wave propagation law, indicating that the algorithm has good adaptability and stability to the complex model.

[0133] Although the present invention has been described in detail with reference to the foregoing embodiments, for those skilled in the art, they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A fast-scanning seismic wave travel time calculation method based on an improved eikonal equation, comprising the following steps: S1: Determine the seismic wave travel time calculation formula, perform multiplicative decomposition, and determine the solution formula; S2: Read in the velocity model and related parameter information, where the parameters include the grid size, grid spacing, and source location of the velocity model; S3: Initialize the source; For the slowness fixed term S0(x,z), set it to the velocity of a homogeneous medium of 1 m / s, then the travel time fixed term T0 is: Formula 11: For the source point (x0,y0) there is: Equation 12: T0(x0,y0) = 0 For the travel time perturbation term τ(x,z); set τ(x0,y0) = 1 at the source point, and assign a relatively large positive value to all other grid points; S4: Perform Gauss-Seidel iterative calculation to solve the discretized equation, In the alternating 4 directions, use the Gauss-Seidel type iterative method to solve the discretized equation, and the scanning calculation directions are: (1) i = 1:I; j = 1:J (2) i = 1:I; j = J:1 (3) i = I:1; j = 1:J (4) i = I:1; j = J:1 S5: Termination condition. Set the threshold δ→0, and the travel time ∑|τ old -τ new | < δ, then stop the iterative calculation and return the calculation result τ(x,z) of the perturbation factor at the same time. Otherwise, repeat S4 until the termination condition is met; finally, the calculation result of the seismic wave travel time is T(x,z) = T0(x,z)τ(x,z).

2. The fast-scanning seismic wave travel time calculation method based on the improved eikonal equation according to claim 1, wherein: In the S1 step, the seismic wave travel time calculation formula is: Formula 1: |▽T(x,z)| 2 = S 2 (x,z) Where: T represents the travel time, S represents the slowness, that is, the reciprocal of the velocity.

3. A fast-scanning seismic wave travel time calculation method based on the improved eikonal equation according to claim 1, characterized in that: The decomposition method of the seismic wave travel time calculation formula is: By performing multiplicative decomposition on Equation 1, decompose the travel time term and the slowness term in the eikonal equation into fixed terms and perturbation terms. Among them, the fixed terms propagate in a spherical form and are known parameters, and the perturbation terms are correction factors and are parameters to be solved, used to solve the source singularity problem. Equation 2: T(x,z) = T0(x,z)τ(x,z) Equation 3: S(x,z) = S0(x,z)α(x,z) When: Formula 4: |▽T0(x,z)| 2 = S0 2 (x,z) Where, T0(x,z) represents the travel time fixed term, τ(x,z) represents the travel time perturbation term, S0(x,z) represents the slowness fixed term, and α(x,z) represents the slowness perturbation term. In Equations 2 and 3, both T0(x,z) and S0(x,z) are known. Using Equations 2, 3, and 4, convert the seismic wave travel time T(x,z) in Equation 1 into a solution problem in the decomposed eikonal equation τ(x,z).

4. A fast-scanning seismic wave travel time calculation method based on an improved eikonal equation according to claim 1, characterized in that: The method for determining the solution formula is: Taking a two-dimensional rectangular grid as an example, illustrate the specific solution formula of the fast-scanning seismic wave travel time calculation method based on the improved eikonal equation; Display an internal grid point C in the two-dimensional rectangular grid. The neighboring points of this point are W, E, N, S. Discretize the equation on the four triangles: ΔCEN, ΔCNW, ΔCWS, and ΔCSE; On the triangle ΔCWS, a discretized equation can be obtained: Formula 6: According to the quadratic formula for finding the roots of a quadratic equation, a new form of the solution for the travel time of seismic waves in the discretized equation formula 6 can be obtained. Let Δ = B 2 - 4AC, and when Δ ≥ 0, solve for τ C which is expressed as: Formula 7: where τ C = x, τ W = a, τ S = b, h represents the grid spacing, T0(x, z) represents the fixed travel time term, τ(x, z) represents the travel time perturbation term, S0(x, z) represents the fixed slowness term, α(x, z) represents the slowness perturbation term, T 0x , T 0y respectively represent the gradients of T0(x, z) in the x and z directions, and their mathematical expressions are: Formula 8: The solution τ of Equation 6 C To satisfy the causality condition, when Equation 6 has no real roots or has real roots that do not satisfy the causality condition, the characteristic method is used to transfer the information of τ along the edge from W to C, and along the edge from S to C; by using the characteristic method on the edges respectively, two values of τ can be calculated: τ WC , τ SC ; when, the edges where δx = h, δy = 0, L = h Formula 9: Formula 10:

5. A fast-scanning seismic wave travel time calculation method based on the improved eikonal equation according to claim 4, characterized in that: In the S5 step The actual method for solving the discretized equation is: At grid point C, the discretization equations corresponding to the four triangles ΔCEN, ΔCNW, ΔCWS, and ΔCSE formed by its neighboring points are solved respectively. For triangle ΔCWS, the solution formula 6 is used to find the two roots τ C,1 and τ C,2 . For two real roots τ C,1 and τ C,2 , first, a causality judgment is made: τ C,1 T0(C) ≥ τ W T0(W), τ C,1 T0(C) ≥ τ S T0(S) τ C,2 T0(C) ≥ τ W T0(W), τ C,2 T0(C) ≥ τ S T0(S) τ C,1 and τ C,2 both satisfy the causality condition, then T WS = min{τ C,1 T0(C), τ C,2 T0(C)}; τ C,1 If the causality condition is satisfied, then T WS = τ C,1 T0(C); τ C,2 If the causality condition is satisfied, then T WS = τ C,2 T0(C); Otherwise, neither root satisfies the causality condition. According to Formula 9 and Formula 10, solve the characteristic equation, and use the method of edge features to obtain the characteristic roots τ WC and τ SC ; Execute the causality relationship τ WC T0(C)≥T W , τ SC T0(C)≥T S ; τ WC and τ SC both satisfy the causal relationship, then T WS = min{τ WC T0(C), τ SC T0(C)}; τ WC If the causality is satisfied, then T WS = τ WC T0(C); τ SC If the causality is satisfied, then T WS = τ SC T0(C); τ WC and τ SC both satisfy the causal relationship, then T WS = ∞; Equation 6 has no real roots. According to Equations 9 and 10, the solutions of the characteristic equation are obtained, and τ is obtained by using the method of side characteristics WC and τ SC ; Execute the causal relationship τ WC T0(C) ≥ T W , τ SC T0(C) ≥ T S . τ WC and τ SC both satisfy the causal relationship, then T WS = min{τ WC T0(C), τ SC T0(C)}; τ WC Satisfy the causality, then T WS = τ WC T0(C); τ SC If the causality is satisfied, then T WS = τ SC T0(C); τ WC and τ SC both satisfy the causal relationship, then T WS = ∞; Finally, select the minimum value from the results calculated from the four triangles and calculate the τ(C) value of the grid point C.

Images