Earthquake first arrival wave travel time calculation method for complex structure and undulating surface

Abstract

An earthquake first-motion wave travel time calculation method for a complex structure and an undulating surface belongs to the technical field of earthquake first-motion wave travel time calculation methods, and comprises the following steps: S1, reading a speed model with a terrain; s2, performing tetrahedron subdivision on the velocity model; s3, travel time calculation is carried out on the tetrahedron where the seismic source is located; s4, searching a vertex of travel time to be calculated; s5, carrying out attribute classification on all points; s6, calculating the travel time of the to-be-calculated point by adopting the travel time of the three known nodes; s7, carrying out travel time calculation on each tetrahedron by using a rapid propulsion method; and S8, outputting the file and displaying a result. According to the method, the wavefront expansion method on the tetrahedral domain is applied to travel time calculation, the tetrahedral mesh generation is flexible, the local adaptive capacity is high, the mesh size-variable generation can be performed according to the surface fluctuation and the construction complexity, the model boundary is consistent with the original interface, and the calculation precision is high.

Description

Technical Field

[0001] The invention belongs to the technical field of earthquake first-arrival wave travel time calculation methods, and in particular relates to an earthquake first-arrival wave travel time calculation method for complex structures and undulating surfaces. Background Art

[0002] Seismic exploration in areas with undulating surfaces and complex structures often faces the challenge of poor seismic data processing. The main reason is that the existing processing technology is based on the assumption of horizontal layered media. Therefore, it is necessary to study seismic data processing technology suitable for complex structural areas. In the process of seismic data processing, seismic travel time calculation is an indispensable part. Therefore, studying the seismic travel time calculation technology under complex structural conditions can lay the foundation for studying seismic processing data adapted to complex structural belts, thereby improving the success rate of seismic exploration in complex structural areas. The traditional basic assumptions of seismic data processing mainly include: uniform medium or horizontal layered medium, flat surface and isotropy, static correction is independent of time, etc. As seismic exploration advances into complex areas, the objects faced by seismic exploration are gradually becoming more complex. In areas with drastic terrain changes and complex structures, the traditional basic assumptions are no longer met. Therefore, the reliability of a series of methods based on traditional assumptions cannot be guaranteed, and we need to improve our methods to adapt to new problems. Seismic travel time field and seismic wave field are two important ways to describe the propagation of seismic waves. The research and development of seismic wave propagation time calculation methods under undulating surface conditions, the propagation time calculated from the surface is more consistent with the actual seismic wave propagation path and has high accuracy. Seismic travel time algorithms mainly include ray tracing method and finite difference method. For ray tracing method, there are mainly two-point ray tracing and shortest path ray tracing methods. Although this method has the advantages of intuitive principle, simple implementation, and numerical stability, it will produce shadows under complex velocity models and slow calculation speed, which is an inherent defect of this type of method. Finite difference method can be divided into two categories according to the different expansion methods of calculation methods: box expansion and wavefront expansion. There are mainly finite difference methods. The fundamental disadvantage of the box expansion method is that the expansion method of differential calculation does not conform to the propagation law of wave field, there is a phenomenon of causal inversion, and it cannot solve the problem of rotation wave. The wavefront expansion method follows the propagation law of wave field and has unconditional stability. At present, the wavefront expansion method for calculating travel time includes the rapid advancement method and the rapid scanning method. These two methods are completely different in the iterative method of solving the eikonal equation, but the difference format is the same. When calculating three-dimensional undulating surfaces and complex structural models, conventional methods use rectangles for calculation, which is difficult to guarantee accuracy. Summary of the invention

[0003] In order to solve the above problems, the present invention proposes: a method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces, comprising the following steps:

[0004] S1, read in the velocity model with terrain;

[0005] S2, tetrahedroning the velocity model;

[0006] S3, calculate the travel time of the tetrahedron where the earthquake source is located;

[0007] S4, finding a vertex for which travel time is to be calculated;

[0008] S5, classify attributes of all points;

[0009] S6. Calculate the travel time of the point to be calculated using the travel time of three known nodes;

[0010] S7, using the rapid marching method to calculate the travel time of each tetrahedron;

[0011] S8. Output file and display the results.

[0012] Furthermore, in step S2, the model is divided into tetrahedrons by tetrahedron partitioning. During the partitioning, it is necessary to ensure that the triangles on each face of the tetrahedron are all acute triangles and the speed in each tetrahedron is consistent.

[0013] Furthermore, in order to improve the accuracy, the vicinity of the earthquake source is densely divided, and finally a tetrahedral model of arbitrary undulating terrain is obtained.

[0014] Furthermore, in the steps S3-S4, the tetrahedron where the earthquake source is located is determined, and the distances between the earthquake source and the four vertices of the tetrahedron are determined. Based on this, the travel times of the four vertices of the tetrahedron where the earthquake source is located are calculated and marked as completion points. The points connected to the tetrahedron are marked as narrowband points, and the other points are marked as far-away points.

[0015] Furthermore, the calculation algorithm for the travel time of each tetrahedral grid point is as follows: the positions of the four vertices of tetrahedron ABCD are known, namely, the coordinates of point A are XA, YA, ZA, and the travel time is TA, the coordinates of point B are XB, YB, ZB, and the travel time is TB, the coordinates of point C are XC, YC, ZC, and the travel time is TC, the coordinates of point D are XD, YD, ZD, and the travel time is TD. Assume that TA, TB, TC have been calculated, TD is the point to be calculated, and the wave propagates from the surface where triangle ABC is located to point D.

[0016] Furthermore, assuming that the travel times of the four points have the following relationship, TA<=TB<=TC<=TD, then:

[0017] TD2((XB(YC-YA)+XC(YA-YB)+XA(YB-YC))2+(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))2+(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB)-ZC

[0018] +2TD(TC(XB(YC-YA)+XC(YA-YB)+XA(YB-YC))(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))+TA(XB(YC-YA)+XC(YA-YB)+XA(YB-Y) C))(XD(YC-YB)+XB(YD-YC)+XC(YB-YD))+TB(XB(YC-YA)+XC(YA-YB)+XA(YB-YC))(XC(YD-YA)+XD(YA-YC)+XA(YC-YD))+

[0019] TC(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XD(ZB-ZA)+XA(ZD-ZB)+

[0020] XB(ZA-ZD))+

[0021] TC(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))(YD(ZB-ZA)+YA(ZD-ZB)+

[0022] YB(FOR-ZD))+

[0023] TA(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XD(ZC-ZB)+XB(ZD-ZC)+

[0024] XC(ZB-ZD))+

[0025] TA(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))(YD(ZC-ZB)+YB(ZD-ZC)+

[0026] YC(ZB-ZD))+

[0027] <h2 style=";text-align:left;direction:ltr">TB(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XC(ZD-ZA)+XD(ZA-ZC)+XA(ZC-ZD))+TB(YB(ZC-ZA) )+YC(ZA-ZB)+YA(ZB-ZC))(YC(ZD-ZA)+YD(ZA-ZC)+YA(ZC-ZD)))+TC2(XD(YB-YA)+XA(YD-Y) B)+XB(YA-YD))2+TB2(XD(YC-YA)+XA(YD-YC)+XC(YA-YD))2+TA2(XD(YC-YB)+XB(YD-YC)+X C(YB-YD))2+TC2(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))2+TB2(XD(ZC-ZA)+XA(ZD-ZC)+XC(ZA) -ZD))2+TC2(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))2+TB2(YD(ZC-ZA)+YA(ZD-ZC)+YC(ZA-ZD) )2+TA2(XD(ZC-ZB)+XB(ZD-ZC)+XC(ZB-ZD))2+TA2(YD(ZC-ZB)+YB(ZD-ZC)+YC(ZB-ZD))2-F 2(-XBYCZA+XBYDZA+XAYCZB-XAYDZB+XD(YB(ZC_ZA)+YC(ZA-ZB)+YA(ZB-ZC))+XBYAZC-XAYB ZC+XAYDZC-XBYDZC+XC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))+(XB(YC-YA)+XA(YB-YC))ZD)2+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0028] <h2 style=";text-align:left;direction:ltr"> 2TATC(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))(XD(YC-YB)+XB(YD-YC)+XC(YB-YD))+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0029] <h2 style=";text-align:left;direction:ltr"> 2TBTC(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))(XC(YD-YA)+XD(YA-YC)+XC(YC-YD))+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0030] <h2 style=";text-align:left;direction:ltr"> 2TATB(XD(YC-YA)+XA(YD-YC)+XC(YA-YD))(XC(YD-YB)+XD(YB-YC)+XB(YC-YD))+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0031] <h2 style=";text-align:left;direction:ltr"> 2TATC(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))(XD(ZC-ZB)+XB(ZD-ZC)+XC(ZB-ZD))+<h2 style=";text-align:left;direction:ltr">

[0032] 2TATC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))(YD(ZC-ZB)+YB(ZD-ZC)+YC(ZB-ZD))+

[0033] 2TBTC(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))(XC(ZD-ZA)+XD(ZA-ZC)+XA(ZC-ZD))+

[0034] 2TATB(XD(ZC-ZA)+XA(ZD-ZC)+XC(ZA-ZD))(XC(ZD-ZB)+XD(ZB-ZC)+XB(ZC-ZD))+

[0035] 2TBTC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))(YC(ZD-ZA)+YD(ZA-ZC)+YA(ZC-ZD))+

[0036] 2TATB(YD(ZC-ZA)+YA(ZD-ZC)+YC(ZA-ZD))(YC(ZD-ZB)+YD(ZB-ZC)

[0037] +YB(ZC-ZD))=0 (2)

[0038] Formula (2) is a quadratic equation about TD. Formula (2) is used to calculate the travel time at the tetrahedron point D. In this way, the calculation in a tetrahedron is completed. In this way, the travel time in the narrowband point is calculated, and the point with the minimum travel time in the narrowband point is marked as the completed point. The tetrahedron where the minimum point is located is found from the far points, and the other vertices of the tetrahedron are marked as narrowband points;

[0039] Furthermore, the travel time within the narrowband point is calculated using formula (2). If the newly calculated travel time is smaller than or equal to the original travel time, the smaller one is retained. If the newly calculated travel time is larger than the original travel time, the travel time of the point does not need to be updated. After the calculation, the point with the smallest travel time within the narrowband point is marked as a completed point, and the tetrahedron where the smallest point is located is found from the far points, and the other vertices of the tetrahedron are marked as narrowband points.

[0040] Furthermore, determine whether a narrowband point exists. If so, return to the step of calculating the travel time within the narrowband point and repeat the cycle calculation. Otherwise, the calculation stops. This cycle continues until all grid points become completion points, and the travel time of all grid points has been calculated.

[0041] The beneficial effects of the present invention are as follows: compared with the prior art, the present invention has the following advantages: flexible and convenient grid generation, the present invention inherits the advantages of the FMM algorithm: the algorithm is unconditionally stable, satisfies the wavefront propagation law, and is a wavefront expansion method with high computational efficiency; the present invention has higher computational accuracy than the traditional orthogonal grid. In view of the defects of the existing calculations, the present invention applies the wavefront expansion method on the tetrahedron domain to the travel time calculation. The tetrahedron grid generation is flexible and has strong local adaptive capabilities. The grid size can be varied according to the surface undulation and the complexity of the structure, and the model boundary is consistent with the original interface. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a flowchart of the present invention;

[0043] Figure 2 It is a schematic diagram of the completed point, narrow band point and far away point after tetrahedron division;

[0044] Figure 3 It is a schematic diagram of the present invention implemented in a single tetrahedron. DETAILED DESCRIPTION

[0045] A method for calculating the travel time of the first arrival of earthquake waves for complex structures and undulating surfaces, relates to a method for calculating the travel time of earthquake waves for complex structures and undulating interfaces in three-dimensional space, and is a method for calculating the travel time of earthquake waves based on tetrahedron decomposition. The purpose of the present invention is to solve the problems of complex calculation and large interface error of the commonly used methods for calculating the travel time of the first arrival of earthquake waves in the case of undulating surfaces and complex structures, avoid the problem that the conventional travel time calculation method needs to smooth the complex interface when calculating the travel time of the first arrival of earthquake waves for the interface of complex structures on three-dimensional undulating surfaces, and improve the calculation accuracy.

[0046] Step 1: Use tetrahedron to divide the model into hexahedrons. When dividing, it is necessary to ensure that the triangles on each face of the tetrahedron are acute triangles and the velocity in each tetrahedron is consistent. In addition, in order to improve the accuracy, the vicinity of the earthquake source can be divided densely. Finally, the terrain model of arbitrary undulations of tetrahedron division is obtained ( Figure 2 );

[0047] Step 2, determine the tetrahedron where the earthquake source is located, determine the distances between the earthquake source and the four vertices of the tetrahedron, calculate the travel times of the four vertices of the tetrahedron where the earthquake source is located, and mark them as completion points, mark the points connected to the tetrahedron as narrowband points, and mark the other points as far-away points;

[0048] Step 3: The following introduces the calculation algorithm of the travel time of each tetrahedral grid point.

[0049] Traditional fast marching algorithms all calculate the eikonal equations within an orthogonal grid.

[0050] |▽T|=s (1)

[0051] In formula (1), T is the travel time of seismic waves, s is the slowness, and ▽ is the gradient operator. The upwind formula can be used to solve the equation in an orthogonal coordinate system, but for a tetrahedron, it cannot be directly applied because it is not an orthogonal grid. In view of this situation, the traditional method is to perform triangular travel time calculations on the four faces of the tetrahedron, and then take the minimum value as the travel time of the point to be calculated. For the specific calculation method, please refer to the solution proposed by Sethian et al. However, this method has only a first-order approximation and low accuracy. In order to improve the accuracy, the present invention improves this algorithm, which is specifically described as follows:

[0052] like Figure 3 As shown, the positions and travel times of the four vertices of tetrahedron ABCD are known, namely, the coordinates of point A are XA, YA, ZA, and the travel time is TA, the coordinates of point B are XB, YB, ZB, and the travel time is TB, the coordinates of point C are XC, YC, ZC, and the travel time is TC, the coordinates of point D are XD, YD, ZD, and the travel time is TD. Assume that TA, TB, and TC have been calculated, and TD is the point to be calculated. The wave propagates from the surface where triangle ABC is located to point D, then there is the following relationship, TA<=TB<=TC<=TD, then there is the following relationship:

[0053] TD2((XB(YC-YA)+XC(YA-YB)+XA(YB-YC))2+(XB(ZC-ZA)+XC(ZA-ZB)+

[0054] XA(ZB-ZC))2+(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))2)

[0055] +2TD(TC(XB(YC-YA)+XC(YA-YB)+XA(YB-YC))(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))+TA(XB(YC-YA)+XC(YA-YB)+XA(YB-Y C))(XD(YC-YB)+XB(YD-YC)+XC(YB-YD))+TB(XB(YC-YA)+XC(YA-YB)+XA(YB-YC))(XC(YD-YA)+XD(YA-YC)+XA(YC-YD))+

[0056] TC(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XD(ZB-ZA)+XA(ZD-ZB)+

[0057] XB(ZA-ZD))+

[0058] <h2 style=";text-align:left;direction:ltr">TC(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))(YD(ZB-ZA)+YA(ZD-ZB)+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0059] <h2 style=";text-align:left;direction:ltr"> YB(ZA-ZD))+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0060] <h2 style=";text-align:left;direction:ltr"> TA(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XD(ZC-ZB)+XB(ZD-ZC)+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0061] <h2 style=";text-align:left;direction:ltr"> XC(ZB-ZD))+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0062] <h2 style=";text-align:left;direction:ltr"> TA(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))(YD(ZC-ZB)+YB(ZD-ZC)+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0063] <h2 style=";text-align:left;direction:ltr"> YC(ZB-ZD))+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0064] <h2 style=";text-align:left;direction:ltr"> TB(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XC(ZD-ZA)+XD(ZA-ZC)+XA(ZC-ZD))+TB(YB(ZC-ZA) )+YC(ZA-ZB)+YA(ZB-ZC))(YC(ZD-ZA)+YD(ZA-ZC)+YA(ZC-ZD)))+TC2(XD(YB-YA)+XA(YD-Y) B)+XB(YA-YD))2+TB2(XD(YC-YA)+XA(YD-YC)+XC(YA-YD))2+TA2(XD(YC-YB)+XB(YD-YC)+X C(YB-YD))2+TC2(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))2+TB2(XD(ZC-ZA)+XA(ZD-ZC)+XC(ZA) -ZD))2+TC2(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))2+TB2(YD(ZC-ZA)+YA(ZD-ZC)+YC(ZA-ZD) )2+TA2(XD(ZC-ZB)+XB(ZD-ZC)+XC(ZB-ZD))2+TA2(YD(ZC-ZB)+YB(ZD-ZC)+YC(ZB-ZD))2-F 2(-XBYCZA+XBYDZA+XAYCZB-XAYDZB+XD(YB(ZC_ZA)+YC(ZA-ZB)+YA(ZB-ZC))+XBYAZC-XAYB ZC+XAYDZC-XBYDZC+XC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))+(XB(YC-YA)+XA(YB-YC))ZD)2+<h2 style=";text-align:left;direction:ltr"> <h2 style=";text-align:left;direction:ltr">

[0065] 2TATC(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))(XD(YC-YB)+XB(YD-YC)+XC(YB-YD))+

[0066] 2TBTC(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))(XC(YD-YA)+XD(YA-YC)+XC(YC-YD))+

[0067] 2TATB(XD(YC-YA)+XA(YD-YC)+XC(YA-YD))(XC(YD-YB)+XD(YB-YC)+XB(YC-YD))+

[0068] 2TATC(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))(XD(ZC-ZB)+XB(ZD-ZC)+XC(ZB-ZD))+

[0069] 2TATC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))(YD(ZC-ZB)+YB(ZD-ZC)+YC(ZB-ZD))+

[0070] 2TBTC(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))(XC(ZD-ZA)+XD(ZA-ZC)+XA(ZC-ZD))+

[0071] 2TATB(XD(ZC-ZA)+XA(ZD-ZC)+XC(ZA-ZD))(XC(ZD-ZB)+XD(ZB-ZC)+XB(ZC-ZD))+

[0072] 2TBTC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))(YC(ZD-ZA)+YD(ZA-ZC)+YA(ZC-ZD))+

[0073] 2TATB(YD(ZC-ZA)+YA(ZD-ZC)+YC(ZA-ZD))(YC(ZD-ZB)+YD(ZB-ZC)

[0074] +YB(ZC-ZD))=0 (2)

[0075] Use formula (2) to calculate the travel time at tetrahedron point D. This completes the calculation in a tetrahedron. In this way, the travel time in the narrowband point is calculated, and the point with the smallest travel time in the narrowband point is marked as the completed point. The tetrahedron where the minimum point is located is found from the far points, and the other vertices of the tetrahedron are marked as narrowband points;

[0076] Step 4: Use formula (2) to calculate the travel time within the narrowband point. If the newly calculated travel time is smaller than or equal to the original travel time, the smaller one is retained. If the newly calculated travel time is larger than the original travel time, the travel time of the point does not need to be updated. After the calculation, mark the point with the smallest travel time within the narrowband point as the completed point, find the tetrahedron where the smallest point is located from the far points, and mark the other vertices of the tetrahedron as narrowband points.

[0077] Step 5: Determine whether a narrowband point exists. If so, return to step 4 and repeat the calculation cycle. Otherwise, the calculation stops.

[0078] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical solutions and concepts of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A method for calculating the travel time of earthquake first arrival waves for complex structures and undulating surfaces, characterized in that: The steps include: S1, read in the velocity model with terrain; S2, tetrahedroning the velocity model; S3, calculate the travel time of the tetrahedron where the earthquake source is located; S4, finding a vertex for which travel time is to be calculated; S5, classify attributes of all points; S6. Calculate the travel time of the point to be calculated using the travel time of three known nodes; S7, using the rapid marching method to calculate the travel time of each tetrahedron; S8. Output file and display the results.

2. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 1, characterized in that: In step S2, the model is divided into tetrahedrons by tetrahedron partitioning. During the partitioning, it is necessary to ensure that the triangles on each face of the tetrahedron are all acute triangles and the speed in each tetrahedron is consistent.

3. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 2, characterized in that: In order to improve the accuracy, the area near the earthquake source is divided densely, and finally a tetrahedral terrain model with arbitrary undulations is obtained.

4. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 3, characterized in that: In the steps S3-S4, the tetrahedron where the earthquake source is located is determined, and the distances from the earthquake source to the four vertices of the tetrahedron are determined. Based on this, the travel times of the four vertices of the tetrahedron where the earthquake source is located are calculated and marked as completion points. The points connected to the tetrahedron are marked as narrowband points, and the other points are marked as far-away points.

5. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 4, characterized in that: The calculation algorithm for the travel time of each tetrahedral grid point is as follows: the positions of the four vertices of tetrahedron ABCD are known, namely, point A has coordinates XA, YA, ZA, and the travel time is TA, point B has coordinates XB, YB, ZB, and the travel time is TB, point C has coordinates XC, YC, ZC, and the travel time is TC, point D has coordinates XD, YD, ZD, and the travel time is TD. Assume that TA, TB, TC have been calculated, TD is the point to be calculated, and the wave propagates from the surface where triangle ABC is located to point D.

6. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 5, characterized in that: Assuming that the travel times of the four points have the following relationship, TA<=TB<=TC<=TD, then: TD2((XB(YC-YA)+XC(YA-YB)+XA(YB-YC))2+(XB(ZC-ZA)+XC(ZA-ZB) + +2TD(TC(XB(YC-YA)+XC(YA-YB)+XA(YB-YC))(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))+TA(XB(YC-YA)+XC(YA-YB)+XA(YB-Y C))(XD(YC-YB)+XB(YD-YC)+XC(YB-YD))+TB(XB(YC-YA)+XC(YA-YB)+XA(YB-YC))(XC(YD-YA)+XD(YA-YC)+XA(YC-YD))+ TC(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XD(ZB-ZA)+XA(ZD-ZB)+ XB(ZA-ZD))+ TC(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))(YD(ZB-ZA)+YA(ZD-ZB)+ YB(ZA-ZD))+ <h2 style=";text-align:left;direction:ltr">TA(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XD(ZC-ZB)+XB(ZD-ZC)+ <h2 style=";text-align:left;direction:ltr"> XC(ZB-ZD))+ <h2 style=";text-align:left;direction:ltr"> TA(YB(ZC-ZA)+YC(ZA-ZB)+YA(ZB-ZC))(YD(ZC-ZB)+YB(ZD-ZC)+ <h2 style=";text-align:left;direction:ltr"> YC(ZB-ZD))+ <h2 style=";text-align:left;direction:ltr"> TB(XB(ZC-ZA)+XC(ZA-ZB)+XA(ZB-ZC))(XC(ZD-ZA)+XD(ZA-ZC)+XA(ZC-ZD))+TB(YB(ZC-ZA) )+YC(ZA-ZB)+YA(ZB-ZC))(YC(ZD-ZA)+YD(ZA-ZC)+YA(ZC-ZD)))+TC2(XD(YB-YA)+XA(YD-Y) B)+XB(YA-YD))2+TB2(XD(YC-YA)+XA(YD-YC)+XC(YA-YD))2+TA2(XD(YC-YB)+XB(YD-YC)+X C(YB-YD))2+TC2(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))2+TB2(XD(ZC-ZA)+XA(ZD-ZC)+XC(ZA) -ZD))2+TC2(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))2+TB2(YD(ZC-ZA)+YA(ZD-ZC)+YC(ZA-ZD) )2+TA2(XD(ZC-ZB)+XB(ZD-ZC)+XC(ZB-ZD))2+TA2(YD(ZC-ZB)+YB(ZD-ZC)+YC(ZB-ZD))2-F 2(-XBYCZA+XBYDZA+XAYCZB-XAYDZB+XD(YB(ZC_ZA)+YC(ZA-ZB)+YA(ZB-ZC))+XBYAZC-XAYB ZC+XAYDZC-XBYDZC+XC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))+(XB(YC-YA)+XA(YB-YC))ZD)2+ <h2 style=";text-align:left;direction:ltr"> 2TATC(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))(XD(YC-YB)+XB(YD-YC)+XC(YB-YD))+ <h2 style=";text-align:left;direction:ltr"> 2TBTC(XD(YB-YA)+XA(YD-YB)+XB(YA-YD))(XC(YD-YA)+XD(YA-YC)+XC(YC-YD))+ <h2 style=";text-align:left;direction:ltr"> 2TATB(XD(YC-YA)+XA(YD-YC)+XC(YA-YD))(XC(YD-YB)+XD(YB-YC)+XB(YC-YD))+ 2TATC(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))(XD(ZC-ZB)+XB(ZD-ZC)+XC(ZB-ZD))+ 2TATC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))(YD(ZC-ZB)+YB(ZD-ZC)+YC(ZB-ZD))+ 2TBTC(XD(ZB-ZA)+XA(ZD-ZB)+XB(ZA-ZD))(XC(ZD-ZA)+XD(ZA-ZC)+XA(ZC-ZD))+ 2TATB(XD(ZC-ZA)+XA(ZD-ZC)+XC(ZA-ZD))(XC(ZD-ZB)+XD(ZB-ZC)+XB(ZC-ZD))+ 2TBTC(YD(ZB-ZA)+YA(ZD-ZB)+YB(ZA-ZD))(YC(ZD-ZA)+YD(ZA-ZC)+YA(ZC-ZD))+ 2TATB(YD(ZC-ZA)+YA(ZD-ZC)+YC(ZA-ZD))(YC(ZD-ZB)+YD(ZB-ZC)+YB(ZC-ZD))=0(2) Formula (2) is a quadratic equation about TD. Formula (2) is used to calculate the travel time at the tetrahedron point D. In this way, the calculation in a tetrahedron is completed. In this way, the travel time in the narrowband point is calculated, and the point with the minimum travel time in the narrowband point is marked as the completed point. The tetrahedron where the minimum point is located is found from the far points, and the other vertices of the tetrahedron are marked as narrowband points.

7. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 6, characterized in that: Formula (2) is used to calculate the travel time within the narrowband point. If the newly calculated travel time is smaller than or equal to the original travel time, the smaller one is retained. If the newly calculated travel time is larger than the original travel time, the travel time of the point does not need to be updated. After the calculation, the point with the smallest travel time within the narrowband point is marked as the completed point, and the tetrahedron where the minimum point is located is found from the far points, and the other vertices of the tetrahedron are marked as narrowband points.

8. The method for calculating the travel time of the first arrival wave of an earthquake for complex structures and undulating surfaces as claimed in claim 7, characterized in that: Determine whether a narrowband point exists. If so, return to the step of calculating the travel time within the narrowband point and repeat the cycle calculation. Otherwise, the calculation stops. This cycle continues until all grid points become completion points. In this case, the travel time of all grid points has been calculated.

Images