Boundary reflected wave suppression method in infinite space dynamics numerical calculation and computer readable storage medium

Abstract

The invention relates to a boundary reflected wave suppression method in infinite space dynamics numerical calculation and a computer readable storage medium. The method comprises the following steps: acquiring a geometric dimension, a seismic source parameter and a basic geological parameter of a numerical model; obtaining a wave field parameter and a discrete grid size; based on the wave field parameters and the discrete grid size, the total analysis duration, the tail end damping area coverage range and the single-layer damping layer size are determined; based on the seismic source parameters, the wave field parameters, the coverage range of the tail end damping area and the size of a single damping layer, the number of damping layers and damping of each layer of the tail end damping area are obtained through calculation, and grids are divided; and setting an analysis step time increment, and finally endowing artificial damping in the numerical model. When the method is used, an elastic wave fluctuation equation does not need to be modified, a complex coordinate system does not need to be stretched and transformed, only artificial damping needs to be given to the tail end of a calculation area after grid subdivision, a seismic wave field at the boundary is gradually dissipated through damping of the tail end area, and the method can be applied to suppression of grazing incident waves without transformation of displacement-speed components.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation of seismic wave fields, and particularly to a method for suppressing boundary reflected waves in numerical calculation of infinite space dynamics and a computer-readable storage medium. Background Art

[0002] Due to the stress complexity and multi-field coupling of deep geological structures, numerical seismic wave field simulation is an important means for studying wave field characteristics in fields such as seismic exploration and structural damage detection. However, due to limitations in computing power and storage space, numerical models often need to be artificially truncated and the boundaries processed, and suppressing reflected waves at the truncated boundaries is the key to correctly simulating the seismic wave field.

[0003] Currently, there are mainly two types of methods for dealing with boundary reflected waves. One is a stress-type artificial boundary established using the approximate solution of the unilateral wave equation and the expression of the outgoing wave in the calculation region. However, the working efficiency of this scheme largely depends on the accurate expression of the outgoing wave, so the effect is poor when suppressing the reflection of spherical waves incident on the free boundary. The other type requires introducing stretched coordinates in the surrounding area of the model and defining an attenuation function in the stretched coordinates to achieve the propagation process of waves from the finite domain to the infinite space. However, there are limitations such as complex implementation, difficulty in selecting the parameters of the stretching function, or difficulty in directly applying it to commercial finite element platforms. Summary of the Invention

[0004] The present invention provides a method for suppressing boundary reflected waves in numerical calculation of infinite space dynamics and a computer-readable storage medium to solve the above technical problems.

[0005] The present invention is realized through the following technical solutions:

[0006] A method for suppressing boundary reflected waves in numerical calculation of infinite space dynamics, comprising the following steps:

[0007] Collect the geometric dimensions of the numerical model, source parameters, and basic geological parameters;

[0008] Calculate the wave field parameters based on the basic geological parameters;

[0009] Calculate the discrete grid size based on the wave field parameters;

[0010] Determine the total analysis duration, the coverage range of the end damping region, and the single-layer damping layer size based on the wave field parameters and the discrete grid size;

[0011] Calculate the number of damping layers and the damping of each layer in the end damping region based on the source parameters, wave field parameters, coverage range of the end damping region, and single-layer damping layer size;

[0012] Divide the numerical model into grids according to the discrete grid size based on the number of damping layers and the damping values of each layer;

[0013] The time increment of the analysis step is calculated based on the discrete grid size and wave field parameters;

[0014] In the numerical model, the grid is converted into an editable grid component, and artificial damping is assigned in the material properties.

[0015] Among them, the wave velocity is calculated by the following formula:

[0016]

[0017] In the above formula, C s represents the S-wave velocity, C p represents the P-wave velocity, E represents the elastic modulus, ρ represents the density of the geological body, and c represents the Poisson's ratio;

[0018] The wavelength is calculated by the following formula:

[0019]

[0020] In the above formula, λ s represents the S-wave wavelength, and f represents the source frequency.

[0021] Optionally, the discrete grid size is calculated by the following formula:

[0022]

[0023] In the above formula, Δx represents the discrete grid size, λ s represents the S-wave wavelength, C s represents the S-wave velocity, and f represents the source frequency.

[0024] Among them, the total analysis duration is calculated by the following formula:

[0025]

[0026] In the above formula, ΔT represents the total analysis duration, x represents the horizontal dimension of the numerical model, and y represents the vertical dimension of the numerical model.

[0027] Optionally, the size of the single-layer damping layer is consistent with the discrete grid size.

[0028] Optionally, the coverage range of the end damping region is 2 to 3 times the S-wave wavelength.

[0029] Optionally, the number of end damping layers is calculated by the following formula:

[0030]

[0031] In the above formula, m represents the number of end damping layers, H is the coverage range of the end damping region, and h is the size of a single damping layer.

[0032] Optionally, the damping of each layer in the end damping region is calculated using the following formula:

[0033]

[0034] ω = 2πf

[0035] In the above formula, D i is the damping of the end damping region, ω represents the circular frequency, f represents the source frequency; m represents the number of end damping layers; Δh i is the sum of the sizes of the first i damping layers; t represents the empirical coefficient of damping increment, and H represents the coverage range of the end damping region.

[0036] Optionally, the time increment is determined by the following formula:

[0037]

[0038] In the above formula, dt is the time increment, Δx represents the discrete grid size, and C p represents the P-wave velocity.

[0039] The computer-readable storage medium provided by the present application stores a computer program thereon, and when the program is executed by a processor, it implements the boundary reflected wave suppression method described above;

[0040] Optionally, the computer program includes a first subroutine, a second subroutine, and a third subroutine;

[0041] When the first subroutine is executed by a processor, it implements the acquisition of the geometric dimensions of the numerical model and the source parameters;

[0042] When the second subroutine is executed by a processor, it implements the acquisition of basic geological parameters, as well as the calculation of wave field parameters, discrete grid size, analysis step time increment, and total analysis duration;

[0043] When the third subroutine is executed by a processor, it implements the acquisition of wave field parameters and source parameters, as well as the calculation of the coverage range of the end damping region, the size of a single damping layer, the number of damping layers, and the damping of each layer in the end damping region;

[0044] Compared with the prior art, the present application has at least the following beneficial effects:

[0045] 1. The present application provides a new boundary treatment solution for the numerical simulation of seismic wave fields. Applying the present invention, there is no need to modify the elastic wave equation, nor to stretch and transform the complex coordinate system. Only after the mesh is divided, artificial damping is given at the end of the calculation region, and the seismic wave field at the boundary is gradually dissipated through the damping in the end region;

[0046] 2. This application is based on the finite element algorithm. When solving the wave equation using Newmark difference, displacement splitting is not required, so it is easy to implement, occupies less memory, and has high numerical efficiency.

[0047] 3. This application only defines the material damping value at the end of the calculation region and can be applied to suppress grazing incident waves without converting displacement-velocity components.

[0048] 4. This application can be extended to numerical simulations of seismic wave fields for internal explosion sources, multi-layer medium models, and three-dimensional models. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] To more clearly illustrate the technical solutions of the embodiments of the present invention, the following will briefly introduce the drawings required in the embodiments. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as limiting the scope. For those of ordinary skill in the art, without creative efforts, other related drawings can also be obtained based on these drawings.

[0050] Figure 1 It is a schematic program flow chart of the method for suppressing boundary reflection waves in the numerical calculation of infinite space dynamics in the present invention;

[0051] Figure 2 It is a diagram of a homogeneous, elastic, isotropic geological model under 2D grazing incidence in an embodiment of the present invention;

[0052] Figure 3 It is a snapshot diagram of the wave field of a homogeneous, elastic, isotropic geological model under 2D grazing incidence in an embodiment of the present invention;

[0053] Figure 4 It is a diagram of a 3-layer transversely isotropic medium geological model under 2D vertical incidence provided in an embodiment of the present invention;

[0054] Figure 5 It is a snapshot diagram of the wave field of a 3-layer transversely isotropic medium geological model under 2D vertical incidence with an artificial truncated boundary in an embodiment of the present invention;

[0055] Figure 6 It is a snapshot diagram of the wave field of a 3-layer transversely isotropic medium geological model under 2D vertical incidence after suppressing reflection waves using the method of the present invention in an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0056] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention.

[0057] It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other. It should also be noted that the various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments, and the same or similar parts among the various embodiments can be referred to each other.

[0058] As Figure 1 shown, the method for suppressing boundary reflection waves in numerical calculation of infinite space dynamics disclosed in this embodiment includes three sets of subroutines: namely, the first subroutine: source loading; the second subroutine: obtaining basic wave field parameters and grid discretization; the third subroutine: damping calculation.

[0059] In the first subroutine, it includes two parts: collecting the geometric dimensions of the numerical model and determining the source frequency, period, and acquisition period. The acquisition of these two parts of parameters is for modulating a suitable excitation source so that the elastic wave can carry sufficient formation information of the target detection area in the reflection characteristics.

[0060] 1.1 First, input the model horizontal coordinates (x0, x1) and vertical coordinates (y0, y1) into the first subroutine.

[0061] 1.2 Then, input the frequency, period, and acquisition period of the quasi-excited source.

[0062] In an exemplary embodiment, a sine function modulated by a Hanning window is used. Among them, for the field of seismic exploration, the number of cycles should be 1; for the field of non-destructive testing, especially for the detection of micro-defects, the number of source cycles should not exceed 3 to avoid being unable to cover formation details due to wavelength scale problems and diffraction phenomena.

[0063] The second step: In the second subroutine, it includes three parts: collecting geological parameters, collecting wave field parameters, and calculating the discrete grid size and total analysis duration. Among them, the geological parameters include elastic modulus, Poisson's ratio, and density, and the wave field parameters include wave velocity and wavelength.

[0064] The second subroutine mainly calculates the wave field parameters, discrete grid size, analysis step time increment, and total analysis duration of the geological body.

[0065] 2.1 In the second subroutine, it is necessary to input three parts of content in sequence: basic geological parameters, the number of grid points required to describe the S-wave wavelength, and the model horizontal coordinates (x0, x1) and vertical coordinates (y0, y1).

[0066] 2.2 Obtain the P- and S-wave velocities (C p , C s ) and S-wave wavelength (λ s ) of the current geological body (layer) from the basic geological parameters;

[0067] 2.3. By determining the ratio of the S-wave wavelength to the number of grid points, the discrete grid size Δx can be obtained; and through the horizontal and vertical dimensions of the numerical model, the required total analysis duration ΔT can be obtained.

[0068] Wave velocity and wavelength are the keys to determining the grid size and total analysis duration. The wave velocity is calculated using the following formula:

[0069]

[0070] In the above formula, C s represents the S-wave velocity, C p represents the P-wave velocity, E represents the elastic modulus, ρ represents the density of the geological body, and v represents the Poisson's ratio.

[0071] To restore the propagation process of seismic waves in the medium as much as possible and ensure the accuracy and integrity of numerical calculations, the grid size is calculated using the following formula:

[0072]

[0073] In the above formula, Δx represents the discrete grid size, λ s represents the S-wave wavelength, C s represents the S-wave velocity, and f represents the source frequency.

[0074] To fully simulate the seismic wave field in the calculation area and ensure the propagation range of the effective wave in the calculation area, the total analysis duration should be the ratio of twice the diagonal size of the model to the smaller wave velocity (S-wave velocity), as shown in the following formula:

[0075]

[0076] In the above formula, ΔT represents the total analysis duration, x represents the horizontal dimension, and y represents the vertical dimension.

[0077] Step 3: The third subroutine is mainly to solve the damping values of each layer in the end damping area, including: inputting wave field parameters, inputting source parameters, calculating the damping coverage, and calculating the number of damping layers and the required damping values of each layer.

[0078] 3.1. Subsequent damping calculations are carried out by obtaining wave field parameters and source parameters.

[0079] 3.2. Determine the coverage of the end damping area and divide the area in the model.

[0080] The end damping coverage area is at the truncated boundary on the periphery of the effective calculation area of the model, and a partial range is delimited as the damping area. The larger the coverage range of the damping area, the better the dissipation effect on seismic waves and the better the suppression effect on boundary reflected waves. On the basis of ensuring the suppression effect, considering the numerical convenience and efficiency, through a large number of numerical examples verification, the coverage range H of the end damping area can be 2 to 3 times of λ s times.

[0081] 3.3, Determine the size of the damping layer

[0082] When the end area damping coverage range H is certain, the thinner the damping layer, the more layers there are, the total damping value of the end area reaches the maximum, and the better the suppression effect on boundary reflected waves. In this embodiment, to ensure the suppression effect on boundary reflected waves and the convenience of later operations, through a large number of numerical examples verification, the size of the damping layer is preferably kept consistent with the discrete grid size, that is, h = Δx.

[0083] 3.4, Calculate the damping value of each layer in the end damping area.

[0084] The incremental damping scheme is adopted in the end damping area. The material damping value gradually increases to the t power from the inner layer to the outer layer of the end damping area. Through this scheme, the sudden reduction of wave velocity caused by excessive initial damping can be avoided, thus generating wave impedance differences and causing reflections.

[0085] The damping value is positively correlated with the dissipation effect on the outgoing wave, but it does not mean that the larger the damping value, the better the dissipation effect on the energy of the outgoing wave. When the damping value is small, the damping value in the end area will not be sufficient to dissipate the energy of the outgoing wave, resulting in boundary reflection; when the damping value is large, when the outgoing wave reaches the large damping material layer, the mechanical energy generated by the disturbance of the outgoing wave caused by the stress imbalance at the boundary node will generate reflections between the material layers. Therefore, a certain coordination must be maintained between the two. To characterize the relationship between the damping increasing characteristics in the damping area and the amplitude attenuation effect of the outgoing wave, a large number of numerical examples were carried out. The damping was fitted with the total amount of wavefront diffusion and absorption attenuation of the outgoing wave, and the relationship between the slope of the damping growth curve and the amplitude attenuation rate of the outgoing wave was obtained. The damping of each layer in the end damping area can be calculated by the following formula:

[0086]

[0087] ω = 2πf

[0088] In the above formula, D i is the damping of the end damping area, ω represents the circular frequency; f represents the source frequency; m represents the number of end damping layers; Δh iIt represents the sum of the damping sizes of the first i layers; t represents the empirical coefficient of the damping increment, which is obtained by fitting the damping growth curve and the peak curve of the node fluctuation energy. From the smallest to the largest seismic source frequency, its value range can be taken as 1.03 - 1.17; H represents the coverage range of the end damping area.

[0089] 3.5. After obtaining the number of damping layers and the damping values to be applied to each damping layer, the numerical model is meshed according to the size.

[0090] 3.6. Based on the discrete grid size and wave field parameters, the time increment of the analysis step is calculated.

[0091] To ensure calculation convergence and the integrity of seismic wave propagation in the medium, the time increment dt can be determined by the following formula:

[0092]

[0093] 3.7. Manually apply the damping assignment.

[0094] To facilitate the damping assignment to the end of the model, the discrete grid model is converted into a grid component, and then artificial material damping is applied layer by layer in the attribute definition. By introducing artificial damping at the unit nodes, the mechanical energy generated by the vibration (perturbation) transfer between nodes in the calculation area is gradually dissipated.

[0095] For the artificial material damping in the end damping area, in addition to the damping property, its mechanical parameters (such as elastic modulus, Poisson's ratio, and density) should be consistent with the effective calculation area. Avoid the change of wave impedance caused by the difference in mechanical properties or elastic parameters between the end damping area and the middle calculation area, so as to prevent the increase of the reflection coefficient.

[0096] In some embodiments, the implementation steps of the method for suppressing the boundary reflected wave in the infinite space dynamic numerical calculation are outlined as follows:

[0097] S1. Obtain a suitable seismic source and input the geometric size of the numerical model and the preset frequencies (f), periods (T), sampling frequencies (F), sampling periods (Ts), and seismic source durations (Tt) of the target seismic source in the first subroutine.

[0098] S2. Obtain the wave field parameters and the total analysis duration, and input the horizontal coordinates (x0, x1), vertical coordinates (y0, y1), target medium density (ρ), elastic modulus (E), Poisson's ratio (v) of the target area in the second subroutine, and calculate the P and S wave velocities (C p , C s ), and S wave wavelength (λ s ) of the target area medium.

[0099] S3. Calculate the coverage range of the end damping region and the size of a single damping layer. The coverage range of the end damping region can be 2 - 3 times of λ s , and the size of the damping layer should be consistent with the grid size (h = Δx).

[0100] S4. Calculate the damping to be applied to each damping layer. Substitute the source frequency, the coverage range of the end damping region, and the size of a single damping layer into the damping solution formula to calculate the damping to be applied to each damping layer. Calculate the total analysis duration ΔT according to the S-wave velocity and the target calculation region range.

[0101]

[0102] S5. Divide the grid, and divide the grid (Δx) according to the S-wave wavelength.

[0103] S6. Set the time increment of the analysis step.

[0104] S7. Assign damping. In the numerical model, convert the grid into an editable grid component and assign artificial damping in the material properties.

[0105] If seismic data without boundary reflections is to be obtained, the following conditions should be ensured during the simulation:

[0106] (1) Maintain the consistency of the elastic parameters of the medium in the calculation region and the end damping region, and avoid the increase in the absolute value of the reflection coefficient between the damping and the calculation region due to the wave impedance difference, resulting in multiple waves.

[0107] (2) The damping assignment scheme of the present invention adopts the incremental method. From the inner edge to the outer edge of the damping region, the damping increases in the t-th power to ensure that the peak value of the outgoing wave energy matches the damping value of the current damping layer.

[0108] (3) For the multi-layer medium model, the physical properties parameters of the outer damping layer of each layer of medium should be the same as those of this layer of medium.

[0109] Next, the suppression effect of the present invention on the boundary reflection wave field is verified through specific embodiments.

[0110] Embodiment 1

[0111] The designed embodiment is a homogeneous, elastic, isotropic geological model with a size of 100m × 100m, as Figure 2 shown. The source frequency is 100Hz, the source is located at the top free boundary, and it has a grazing incidence at an angle of 30° from northwest to southeast. The sampling interval is 2.5E-4s, and the calculation duration is 0.4s; C s is 1020m / s; λ s is 10.2m. For the end damping region, H is taken as 2λ s ; t is taken as 1.07; Δx is taken as The numerical results are as follows Figure 3 As shown, it can be seen that the present invention can better suppress the boundary reflection wave field and eliminate interference waves such as multiple waves.

[0112] Example 2

[0113] The design example is a three-layer transversely isotropic geological model with a size of 100m×100m. Figure 4 The source frequency is 100 Hz, the excitation point is located in the middle of the top free boundary, the sampling interval is 2.5E-4s, and the recording time is 0.35s; the formation medium is C from top to bottom s They are 1020m / s, 1455m / s, 1940m / s respectively; s They are 10.2m, 14.5m and 19.4m respectively. For the end damping area, H is 2λ s ; t is 1.07; Δx is The numerical results are as follows Figure 5 , Figure 6 As shown, Figure 5 6 is a snapshot of the wave field under the artificial truncation boundary, and 7 is a snapshot of the wave field after the reflection wave is suppressed by the present invention. It can be seen that the seismic wave propagates orderly in the geological model and diffuses into the infinite space without generating reflection on the truncation boundary. The present invention can still be well applied to the seismic wave field simulation of the multi-layer medium geological body.

[0114] The present application is a boundary reflection wave suppression method applicable to technical fields such as numerical calculation of infinite space dynamics problems, acoustic logging, and numerical simulation of non-destructive testing of underground structures. It is applicable to the simulation of seismic wave fields with vertical incidence, grazing incidence, and internal explosion sources at free boundaries under single-layer and multi-layer media. The use of the present invention does not require modification of the elastic wave equation, nor does it require stretching and transforming the complex coordinate system and performing displacement-wave velocity splitting. The reflected wave amplitude at the boundary is attenuated only by assigning incremental artificial material damping at the end of the calculation area. The implementation is simple, the memory usage is small, and the numerical efficiency is high. At the same time, the present invention can be extended to three-dimensional seismic wave field numerical simulation.

[0115] In several embodiments provided by the present application, it should be understood that the disclosed devices and methods can also be implemented in other ways. The device embodiments described above are merely illustrative. For example, the flowcharts and block diagrams in the accompanying drawings show the possible architectures, functions, and operations of devices, methods, and computer program products according to multiple embodiments of the present application. In this regard, each block in the flowchart or block diagram may represent a module, a program segment, or a part of code, and the module, program segment, or part of code contains one or more executable instructions for implementing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the blocks may occur in a different order from that marked in the accompanying drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in the reverse order, depending on the functions involved. It should also be noted that each block in the block diagram and / or flowchart, and the combination of blocks in the block diagram and / or flowchart, can be implemented by a dedicated hardware-based system that performs the specified functions or actions, or can be implemented by a combination of dedicated hardware and computer instructions.

[0116] In addition, the functional modules in each embodiment of the present application may be integrated together to form an independent part, or each module may exist separately, or two or more modules may be integrated to form an independent part.

[0117] When the above-mentioned functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in various embodiments of this application. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories (ROM, Read-Only Memory), random access memories (RAM, Random Access Memory), magnetic disks, or optical discs that can store program codes. It should be noted that in this embodiment, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variation thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or also includes elements inherent to such process, method, article or device. Without further limitation, an element defined by the statement "comprising an..." does not exclude the existence of additional identical elements in the process, method, article or device comprising the element.

[0118] The above are only the preferred embodiments of the present invention and are not used to limit the present invention. For those skilled in the art, the present invention can have various changes and modifications. Any modification, equivalent replacement, improvement, 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 method for suppressing boundary reflection waves in numerical calculation of infinite space dynamics, characterized in that: The following steps are involved: Collect numerical model geometry, source parameters and basic geological parameters; Wave field parameters are calculated based on basic geological parameters; Based on the wave field parameters, the discrete grid size is calculated; Determine the total analysis time, the coverage of the terminal damping area, and the size of the single damping layer based on the wave field parameters, model geometry, size, and discrete grid size; Based on the source parameters, wave field parameters, the coverage of the terminal damping area and the size of the single damping layer, the number of damping layers and the damping of each layer in the terminal damping area are calculated; Based on the number of damping layers and the damping value of each layer, the numerical model is meshed according to the discrete grid size; The time increment of the analysis step is calculated based on the discrete grid size and wave field parameters; In the numerical model, the meshed model is converted into an editable part and artificial damping is assigned in the material properties.

2. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: Source parameters include frequency, period and acquisition period.

3. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: The wave field parameters include S-wave velocity, P-wave velocity and S-wave wavelength; The wave speed is calculated using the following formula: In the above formula, C s represents the S wave velocity, C p represents the P-wave velocity, E represents the elastic modulus, ρ represents the density of the geological body, and v represents the Poisson's ratio; The wavelength is calculated using the following formula: In the above formula, λ s represents the wavelength of S wave, and f represents the source frequency.

4. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: The discrete grid size is calculated using the following formula: In the above formula, Δx represents the discrete grid size, λ s represents the wavelength of S wave, C s represents the S-wave velocity, and f represents the source frequency.

5. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: The total analysis time was calculated using the following formula: In the above formula, ΔT represents the total analysis time, x represents the horizontal dimension of the numerical model, and y represents the vertical dimension of the numerical model.

6. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: The size of the single damping layer is consistent with the discretization grid size.

7. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1 or 6, characterized in that: The end damping area covers 2 to 3 times the wavelength of the S wave; Optionally, the number of terminal damping layers can be calculated using the following formula: In the above formula, m represents the number of terminal damping layers, H is the coverage of the terminal damping area, and h is the size of a single damping layer.

8. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: The damping of each layer in the terminal damping area is calculated using the following formula: ω=2πf In the above formula, D i is the damping of the i-th layer in the terminal damping area, ω represents the circular frequency, f represents the source frequency; m represents the number of terminal damping layers; Δh i is the damping size and of the first i layers; t represents the empirical coefficient of damping increment; and H represents the coverage of the terminal damping area.

9. The method for suppressing boundary reflection waves in infinite space dynamics numerical calculation according to claim 1, characterized in that: The time increment is determined by the following formula: In the above formula, dt is the time increment, Δx represents the discrete grid size, and C p Indicates the P wave velocity.

10. A computer storable medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the boundary reflection wave suppression method according to any one of claims 1 to 9 is implemented; optionally, the computer program includes a first subroutine, a second subroutine and a third subroutine; When the first subroutine is executed by the processor, the acquisition of geometric dimensions of the numerical model and source parameters is realized; When the second subroutine is executed by the processor, the basic geological parameters are collected, and the wave field parameters, the discrete grid size, the analysis step time increment and the total analysis time are calculated; When the processor executes the third subroutine, the processor collects wave field parameters and source parameters, and calculates the coverage of the terminal damping area, the size of a single damping layer, the number of damping layers, and the damping of each layer in the terminal damping area.

Images