Joint absorption boundary method applied to low-frequency ship seismic wave numerical simulation

Abstract

The invention relates to a joint absorption boundary method applied to low-frequency ship seismic wave numerical simulation, and belongs to the field of ship seismic wave numerical simulation. In the numerical simulation process of the method, firstly, a simulation area is divided into a center wave field area and an artificial boundary area; in the central wave field area, taking the Ricker wavelet as a ship noise sound source, simulating sound source disturbance based on a Ricker wavelet expression, and performing ship seismic wave numerical simulation by adopting a staggered grid finite difference method; the ship seismic wave equation is described by a two-dimensional first-order speed-stress elastic wave equation; in a central wave field area, staggered grid finite difference discretization is carried out on an elastic wave equation, and an expansion form of time second-order and space 2N-order precision is adopted; in the artificial boundary area, L completely matched layers are added to serve as absorption boundaries, and a Gaussian attenuation function is adopted in a PML layer; and jointly using a second-order Higdon absorption boundary at the outer boundary. According to the method, the overall boundary absorption effect can be effectively improved under the condition that the storage consumption is not increased.

Description

Technical Field

[0001] The invention belongs to the field of numerical simulation of ship seismic waves, and in particular relates to a joint absorbing boundary method applied to numerical simulation of low-frequency ship seismic waves. Background Art

[0002] Numerical simulation of ship seismic waves is a method that simplifies the real ocean environment into a mathematical model and uses numerical calculations to simulate the propagation of ship seismic waves in the model. Through numerical simulations of various models, the propagation laws of ship seismic waves in different seabed media can be analyzed, thereby providing a theoretical basis for the actual underwater long-distance target detection based on ship seismic waves. In the computer-based finite difference numerical simulation of ship seismic waves, due to the limited memory of the computer, it is necessary to manually define the calculation area, and simply truncating the boundary will cause strong reflection wave interference at the boundary, significantly reducing the simulation accuracy. Therefore, it is usually necessary to properly process the artificial boundary to eliminate the false reflection at the calculation boundary, so that the numerical simulation results are more in line with the actual situation.

[0003] One-way wave absorption boundary and perfectly matched layer (PML) absorption boundary are the two main methods for artificial boundary processing at present. One-way wave absorption boundary is a boundary condition method that uses the one-way wave equation to absorb the energy of the incident wave in the boundary area, mainly including Clayton-Enquist (CE) absorption boundary and Higdon absorption boundary. Compared with CE absorption boundary, Higdon absorption boundary usually has better absorption effect, and its differential calculation is more convenient. PML absorption boundary achieves the absorption of false reflections by adding a perfectly matched layer in the boundary area and using an equation containing an attenuation function in the layer. This method can theoretically absorb incident waves at all angles, so it has received widespread attention in the industry in recent years. Meng Luwen, Zhu Xufang and others successively used PML absorption boundaries to realize the finite difference numerical simulation calculation of ship seismic waves, and obtained good simulation results on the basis of suppressing artificial boundary reflections.

[0004] However, the existing PML absorption boundary still has the following shortcomings in practical applications: First, the exponential attenuation function selected by the conventional PML absorption boundary increases too fast at the inner boundary, which leads to poor coupling between the central wave field and the inner boundary wave field, thereby generating false reflections at the inner boundary; second, the incident wave is usually not completely absorbed after passing through the PML attenuation layer, and some residual boundary reflections will still be generated at the outer boundary of the PML. In summary, the current conventional PML absorption boundary has not yet been able to solve the problem of false reflections at the inner and outer boundaries in the boundary layer, and its comprehensive absorption effect still needs to be further improved. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves. When the present invention performs numerical simulation of ship seismic waves based on the finite difference method, a Gaussian attenuation function is first used in the PML boundary area to optimize the coupling effect between the central wave field and the inner boundary, effectively reducing the inner boundary reflection; at the same time, a second-order Higdon absorbing boundary is used in combination at the outer boundary of the artificial boundary to further absorb the residual incident wave. Compared with the conventional PML boundary, the joint absorbing boundary used in the present invention can comprehensively suppress the false reflections of the inner and outer boundaries without increasing the storage capacity, thereby significantly improving the numerical simulation accuracy of ship seismic waves.

[0006] The present invention adopts the following technical solutions:

[0007] A joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves includes the following steps:

[0008] (1) In the numerical simulation process, the simulation area is first divided into a central wave field area and an artificial boundary area. In the central wave field area, the Ricker wavelet is used as the ship noise source, and the staggered grid finite difference method is used to perform numerical simulation of ship seismic waves; the ship seismic wave equation is described by a two-dimensional first-order velocity-stress elastic wave equation;

[0009] As a specific implementation, the two-dimensional first-order velocity-stress elastic wave equation is specifically expressed as follows:

[0010]

[0011] In formula (1), t is the time, ρ is the model medium density, w is the excitation sound source, x and z represent the horizontal and vertical directions in the space coordinates respectively, and v x 、v z Respectively represent the horizontal and vertical components of velocity, τ xx , τ zz , τ xz They represent the three components of stress respectively, λ and μ are the Lame coefficients, and the calculation formula is as follows:

[0012]

[0013] In formula (2), v p 、v s are the longitudinal wave velocity and the transverse wave velocity respectively.

[0014] (2) In the central wave field region, the elastic wave equation is discretized by staggered grid finite difference, using an expansion form with 2nd order accuracy in time and 2Nth order accuracy in space;

[0015] As a specific implementation method, the elastic wave equation is discretized by staggered grid finite difference as follows:

[0016]

[0017] In formula (3), Δt is the time sampling interval, b m is the spatial difference coefficient, k is the time discrete point number, i and j represent the sampling point numbers in the x and z directions in the spatial coordinates, respectively, and Δx and Δz represent the spatial grid step size, respectively.

[0018] (3) In the artificial boundary region, an L-layer fully matched layer is added as an absorption boundary. In order to optimize the coupling effect between the central wave field and the inner boundary and further suppress the inner boundary reflection, a Gaussian attenuation function is used in the PML layer. As a specific implementation, the Gaussian attenuation function is expressed as follows:

[0019]

[0020] In formula (4), d(h) is the value of the attenuation function that changes with distance, v p is the longitudinal wave velocity, h is the distance from the calculation point in the PML to the outermost boundary of the PML, and R is the theoretical reflection coefficient. The physical quantities in the two-dimensional first-order velocity-stress elastic wave equation are decomposed into two parts, the horizontal component and the vertical component, and the following form is obtained by finite difference discretization:

[0021]

[0022] In formula (5), i and j are the discrete point numbers in the horizontal and vertical directions of space, respectively, and k is the discrete point number in time. x 、v z are the horizontal and vertical components of velocity respectively. xx , σ zz are the normal stresses in the horizontal and vertical directions, σ xz is the shear stress and m represents the spatial number. Represents the horizontal component of velocity at the kth moment at the spatial position (i+1 / 2,j). Represents the horizontal component of velocity at the spatial position (i+1 / 2,j) at the k-1th moment. It represents the normal stress in the horizontal direction at the spatial position (i+m,j) at the k-1 / 2th moment. It represents the normal stress in the horizontal direction at the spatial position (i-m+1,j) at the k-1 / 2th moment. Represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j+m-1 / 2). It represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j-m+1 / 2). Represents the vertical component of the velocity at the kth moment at the spatial position (i+1 / 2,j). Represents the vertical component of the velocity at the k-1th moment at the spatial position (i+1 / 2,j). Represents the shear stress at the k-1 / 2th moment at the spatial position (i+m-1 / 2,j+1 / 2). Represents the shear stress at the k-1 / 2th moment at the spatial position (i-m+1 / 2,j+1 / 2). It represents the normal stress in the vertical direction at the spatial position (i,j+m) at the k-1 / 2th moment. It represents the normal stress in the vertical direction at the spatial position (i,j-m+1) at the k-1 / 2th moment. It represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1 / 2th moment. It represents the normal stress in the horizontal direction at the spatial position (i, j) at the k-1 / 2th moment. Represents the horizontal component of velocity at the kth moment at the spatial position (i+m-1 / 2,j). Represents the horizontal component of velocity at the kth moment at the spatial position (i-m+1 / 2,j). Represents the vertical component of the velocity at the kth moment at the spatial position (i,j+m-1 / 2). Represents the vertical component of the velocity at the kth moment at the spatial position (i,j-m+1 / 2). It represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1 / 2th moment. It represents the normal stress in the vertical direction at the spatial position (i, j) at the k-1 / 2th moment. Represents the shear stress at the k+1 / 2th moment at the spatial position (i+1 / 2,j+1 / 2). Represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j+1 / 2). Represents the horizontal velocity component at the kth moment at the spatial position (i+1 / 2,j+m). Represents the horizontal velocity component at the kth moment at the spatial position (i+1 / 2,j-m+1). Represents the vertical velocity component at the kth moment at the spatial position (i+m,j+1 / 2). represents the vertical velocity component at the kth moment at the spatial position (i-m+1,j+1 / 2). Similarly, Represents the horizontal velocity component at the spatial position (i, j) at the k+1th moment. Represents the vertical velocity component at the k+1th moment at the spatial position (i, j). Represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1th moment. Represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1th moment. represents the shear stress at the spatial position (i, j) at the k+1th moment. Represents the horizontal velocity component at the spatial position (i+1,j) at the k+1th moment. Represents the horizontal velocity component at the spatial position (i+2,j) at the k+1th moment.

[0023] (4) To absorb the residual incident wave at the outer boundary of the PML, a second-order Higdon absorbing boundary is used in combination at the outer boundary. As a specific implementation, the left boundary equation of the second-order Higdon absorbing boundary condition is expressed as:

[0024]

[0025] In formula (6), u x 、u z represent the wave field components in the x and z directions respectively, v p 、v s are the longitudinal wave velocity and the transverse wave velocity respectively; the Talyor formula is used to perform finite difference discretization on equation (6) and applied to the first-order velocity-stress elastic wave equation to obtain v x 、v z , τ xx , τ zz , τ xz The difference formats are:

[0026]

[0027] In formula (7), Represents the horizontal component of velocity at the spatial position (i, j) at the k+1th moment. Represents the horizontal component of velocity at the klth moment at the spatial position (i, j). Represents the vertical component of the velocity at the spatial position (i, j) at the k+1th moment. Represents the vertical component of the velocity at the klth moment at the spatial position (i, j). Represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1th moment. represents the normal stress in the horizontal direction at the spatial position (i, j) at the klth moment. Represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1th moment. Represents the normal stress in the vertical direction at the spatial position (i, j) at the klth moment. represents the shear stress at the spatial position (i, j) at the k+1th moment. represents the shear stress at the klth moment at the spatial position (i, j).

[0028] The coefficients in formula (7) can be obtained by the following identity:

[0029]

[0030] In formula (8), b is usually a constant between 0.3 and 0.5, and v p is the longitudinal wave velocity, Δt is the time sampling interval, and Δx is the spatial grid step. The right boundary, upper boundary, and lower boundary can be derived by analogy with the above derivation process.

[0031] The beneficial effects of the present invention compared with the prior art are as follows:

[0032] Compared with the exponential attenuation function used in the conventional PML absorption boundary, the Gaussian attenuation function changes more slowly at the inner boundary, which makes the coupling effect between the central wave field and the inner boundary better and significantly reduces the inner boundary reflection. In addition, the combined absorption boundary proposed in the present invention is based on the PML absorption boundary using the Gaussian attenuation function, and then combined with the second-order Higdon absorption boundary to perform secondary absorption on the incident wave that is not completely attenuated at the outer boundary of the PML, which can effectively improve the overall boundary absorption effect without increasing storage consumption. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 This is a schematic diagram of the numerical simulation area of ​​ship seismic waves;

[0034] Figure 2 The wave field snapshots of the numerical simulation of ship seismic waves using different artificial absorbing boundaries at the time of 500ms; (a) is the wave field snapshot of the artificial absorbing boundary without any treatment, (b) is the wave field snapshot of the conventional PML absorbing boundary with an exponential attenuation function, and (c) is the wave field snapshot using the joint absorbing boundary;

[0035] Figure 3 Comparison diagram of single-channel seismic records at (425m, 975m) under different artificial absorption boundary conditions. DETAILED DESCRIPTION

[0036] The technical solution of the present invention is further explained below by examples in combination with the accompanying drawings, but the protection scope of the present invention is not limited in any form.

[0037] Example 1

[0038] The present invention adopts a joint absorption boundary method for numerical simulation of low-frequency ship seismic waves. In the implementation process, the simulation area is divided into a central wave field area and an area within an artificial boundary, and a joint absorption boundary is constructed in the artificial boundary area, thereby effectively improving the boundary absorption effect of the numerical simulation. Figure 1 As shown, the specific implementation is as follows:

[0039] (1) The uniform medium velocity model is selected for numerical simulation. The number of central wave field grids is 391*391, and the spatial grid step is 5m. The longitudinal wave velocity is 2500m / s, the shear wave velocity is 1500m / s, and the density is 2000kg / m 3 , the time sampling interval is 0.0005s. A signal generator is used to obtain a Ricker wavelet with a main frequency of 30Hz, and it is excited at the center of the wave field. The numerical simulation of ship seismic waves is performed based on the staggered grid finite difference method. The ship seismic wave equation can usually be described by a two-dimensional first-order velocity-stress elastic wave equation, and its specific expression is as follows:

[0040]

[0041] In formula (1), t is the time, ρ is the model medium density, w is the excitation sound source, x and z represent the horizontal and vertical directions in the space coordinates respectively, and v x 、v z Respectively represent the horizontal and vertical components of velocity, τ xx , τ zz , τ xz They represent the three components of stress respectively, λ and μ are the Lame coefficients, and the calculation formula is as follows:

[0042]

[0043] In formula (2), v p 、v s are the longitudinal wave velocity and the transverse wave velocity respectively.

[0044] (2) In the central wave field region, the elastic wave equation is discretized by staggered grid finite difference, using an expansion form with 2nd order accuracy in time and 2Nth order accuracy in space;

[0045] As a specific implementation method, the elastic wave equation is discretized by staggered grid finite difference as follows:

[0046]

[0047] In formula (3), Δt is the time sampling interval, b mis the spatial difference coefficient, k is the time discrete point number, i and j represent the sampling point numbers in the x and z directions in the spatial coordinates, respectively, and Δx and Δz represent the spatial grid step size, respectively.

[0048] (3) In the artificial boundary region, L layers of perfectly matched layers are added as absorption boundaries. In this implementation, L=40. In order to optimize the coupling effect between the central wave field and the inner boundary and further suppress the inner boundary reflection, a Gaussian attenuation function is used in the PML layer. As a specific implementation method, the Gaussian attenuation function is expressed as follows:

[0049]

[0050] In formula (4), d(h) is the value of the attenuation function that changes with distance, v p is the longitudinal wave velocity, h is the distance from the calculation point in the PML to the outermost boundary of the PML, and R is the theoretical reflection coefficient. The physical quantities in the two-dimensional first-order velocity-stress elastic wave equation are decomposed into two parts, the horizontal component and the vertical component, and the following form is obtained by finite difference discretization:

[0051]

[0052] In formula (5), i and j are the discrete point numbers in the horizontal and vertical directions of space, respectively, and k is the discrete point number in time. x 、v z are the horizontal and vertical components of velocity respectively. xx , σ zz are the normal stresses in the horizontal and vertical directions, σ xz is the shear stress and m represents the spatial number. Represents the horizontal component of velocity at the kth moment at the spatial position (i+1 / 2,j). Represents the horizontal component of velocity at the spatial position (i+1 / 2,j) at the k-1th moment. It represents the normal stress in the horizontal direction at the spatial position (i+m,j) at the k-1 / 2th moment. It represents the normal stress in the horizontal direction at the spatial position (i-m+1,j) at the k-1 / 2th moment. Represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j+m-1 / 2). It represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j-m+1 / 2). Represents the vertical component of the velocity at the kth moment at the spatial position (i+1 / 2,j). Represents the vertical component of the velocity at the spatial position (i+1 / 2,j) at the k-1th moment. Represents the shear stress at the k-1 / 2th moment at the spatial position (i+m-1 / 2,j+1 / 2). Represents the shear stress at the k-1 / 2th moment at the spatial position (i-m+1 / 2,j+1 / 2). It represents the normal stress in the vertical direction at the spatial position (i,j+m) at the k-1 / 2th moment. It represents the normal stress in the vertical direction at the spatial position (i,j-m+1) at the k-1 / 2th moment. It represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1 / 2th moment. It represents the normal stress in the horizontal direction at the spatial position (i, j) at the k-1 / 2th moment. Represents the horizontal component of velocity at the kth moment at the spatial position (i+m-1 / 2,j). Represents the horizontal component of velocity at the kth moment at the spatial position (i-m+1 / 2,j). Represents the vertical component of the velocity at the kth moment at the spatial position (i,j+m-1 / 2). Represents the vertical component of velocity at the kth moment at the spatial position (i,j-m+1 / 2). It represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1 / 2th moment. It represents the normal stress in the vertical direction at the spatial position (i, j) at the k-1 / 2th moment. Represents the shear stress at the k+1 / 2th moment at the spatial position (i+1 / 2,j+1 / 2). Represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j+1 / 2). Represents the horizontal velocity component at the kth moment at the spatial position (i+1 / 2,j+m). Represents the horizontal velocity component at the kth moment at the spatial position (i+1 / 2,j-m+1). Represents the vertical velocity component at the kth moment at the spatial position (i+m,j+1 / 2). represents the vertical velocity component at the kth moment at the spatial position (i-m+1,j+1 / 2). Similarly, Represents the horizontal velocity component at the spatial position (i, j) at the k+1th moment. Represents the vertical velocity component at the k+1th moment at the spatial position (i, j). Represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1th moment. Represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1th moment. represents the shear stress at the spatial position (i, j) at the k+1th moment. Represents the horizontal velocity component at the spatial position (i+1,j) at the k+1th moment. Represents the horizontal velocity component at the spatial position (i+2,j) at the k+1th moment.

[0053] (4) To absorb the residual incident wave at the outer boundary of the PML, a second-order Higdon absorbing boundary is used in combination at the outer boundary. As a specific implementation, the left boundary equation of the second-order Higdon absorbing boundary condition is expressed as:

[0054]

[0055] In formula (6), u x 、u z represent the wave field components in the x and z directions respectively, v p 、v s are the longitudinal wave velocity and the transverse wave velocity respectively; the Talyor formula is used to perform finite difference discretization on equation (6) and applied to the first-order velocity-stress elastic wave equation to obtain v x 、v z , τ xx , τ zz , τ xz The difference formats are:

[0056]

[0057] In formula (7), Represents the horizontal component of velocity at the spatial position (i, j) at the k+1th moment. Represents the horizontal component of velocity at the klth moment at the spatial position (i, j). Represents the vertical component of the velocity at the spatial position (i, j) at the k+1th moment. Represents the vertical component of the velocity at the klth moment at the spatial position (i, j). Represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1th moment. represents the normal stress in the horizontal direction at the spatial position (i, j) at the klth moment. Represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1th moment. Represents the normal stress in the vertical direction at the spatial position (i, j) at the klth moment. represents the shear stress at the spatial position (i, j) at the k+1th moment. represents the shear stress at the klth moment at the spatial position (i, j).

[0058] The coefficients in formula (7) can be obtained by the following identity:

[0059]

[0060] In formula (8), b is usually a constant between 0.3 and 0.5, and v p is the longitudinal wave velocity, Δt is the time sampling interval, and Δx is the spatial grid step. The right boundary, upper boundary, and lower boundary can be derived by analogy with the above derivation process.

[0061] Figure 2 Comparison of wavefield snapshots at 500ms of numerical simulation of ship seismic waves using different artificial absorbing boundaries, where (a) is a wavefield snapshot of an artificial absorbing boundary without any processing, (b) is a wavefield snapshot of a conventional PML absorbing boundary with an exponential attenuation function, and (c) is a wavefield snapshot using a combined absorbing boundary. Figure 2 It can be seen that compared with the artificial absorption boundary without any treatment, the conventional PML absorption boundary has a better suppression effect on the boundary reflection, but the obvious residual boundary reflection can still be seen at the inner and outer boundaries of the PML; when the joint absorption boundary is used to realize the numerical simulation calculation of ship seismic waves, the reflected waves at the inner and outer boundaries of the artificial boundary area are effectively absorbed, and the overall simulation accuracy is significantly improved. Figure 3 It can also be seen from the comparison of the single-channel seismic record diagrams provided that when the joint absorption boundary is used, the amplitude of the internal and external reflections of the artificial absorption boundary is significantly smaller than the artificial absorption boundary without any processing and the conventional PML absorption boundary. In summary, the joint absorption boundary method for numerical simulation of low-frequency ship seismic waves proposed in the present invention effectively improves the boundary absorption effect, significantly improves the numerical simulation accuracy of ship seismic waves, and provides effective technical support for further improving the accuracy of underwater target detection.

Claims

1. A joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves, characterized in that: The method comprises the following steps: (1) In the numerical simulation process, the simulation area is first divided into a central wave field area and an artificial boundary area. In the central wave field area, the Ricker wavelet is used as the sound source of ship noise, and the sound source disturbance is simulated based on the Ricker wavelet expression. The staggered grid finite difference method is used to perform numerical simulation of ship seismic waves. The wave equation of ship seismic waves is described by a two-dimensional first-order velocity-stress elastic wave equation. (2) In the central wave field region, the elastic wave equation is discretized by staggered grid finite difference, using an expansion form with 2nd order accuracy in time and 2Nth order accuracy in space; (3) In the artificial boundary region, an L-layer perfectly matched layer is added as the absorption boundary, and a Gaussian attenuation function is used in the PML layer; (4) In order to absorb the residual incident wave at the outer boundary of the PML, a second-order Higdon absorbing boundary is used in combination at the outer boundary.

2. The joint absorbing boundary method according to claim 1, characterized in that: The specific expression of the two-dimensional first-order velocity-stress elastic wave equation in step (1) is as follows: In formula (1), t is the time, ρ is the model medium density, w is the excitation sound source, x and z represent the horizontal and vertical directions in the space coordinates respectively, and v x 、v z Respectively represent the horizontal and vertical components of velocity, τ xx , τ zz , τ xz They represent the three components of stress respectively, λ and μ are the Lame coefficients, and the calculation formula is as follows: In formula (2), v p 、v s are the longitudinal wave velocity and the transverse wave velocity respectively.

3. The joint absorbing boundary method according to claim 1, characterized in that: In the step (2), the elastic wave equation is discretized by staggered grid finite difference as follows: In formula (3), Δt is the time sampling interval, b m is the spatial difference coefficient, k is the time discrete point number, i and j represent the sampling point numbers in the x and z directions in the spatial coordinates, respectively, and Δx and Δz represent the spatial grid step size, respectively.

4. The joint absorbing boundary method according to claim 1, characterized in that: The Gaussian attenuation function in step (3) is expressed as follows: In formula (4), d(h) is the value of the attenuation function that changes with distance, v p is the longitudinal wave velocity, h is the distance from the calculation point in the PML to the outermost boundary of the PML, and R is the theoretical reflection coefficient. The physical quantities in the two-dimensional first-order velocity-stress elastic wave equation are decomposed into two parts, the horizontal component and the vertical component, and the following form is obtained by finite difference discretization: In formula (5), i and j are the discrete point numbers in the horizontal and vertical directions of space, respectively, and k is the discrete point number in time; v x 、v z are the horizontal and vertical components of velocity respectively; σ xx , σ zz are the normal stresses in the horizontal and vertical directions, σ xz is the shear stress, m represents the spatial number; represents the horizontal component of velocity at the kth moment at the spatial position (i+1 / 2, j); represents the horizontal component of velocity at the k-1th moment at the spatial position (i+1 / 2, j); represents the normal stress in the horizontal direction at the k-1 / 2th moment at the spatial position (i+m,j); represents the normal stress in the horizontal direction at the k-1 / 2th moment at the spatial position (i-m+1,j); represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j+m-1 / 2); represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2,j-m+1 / 2); represents the vertical component of velocity at the kth moment at the spatial position (i+1 / 2, j); represents the vertical component of velocity at the k-1th moment at the spatial position (i+1 / 2, j); represents the shear stress at the k-1 / 2th moment at the spatial position (i+m-1 / 2,j+1 / 2); represents the shear stress at the k-1 / 2th moment at the spatial position (i-m+1 / 2,j+1 / 2); represents the normal stress in the vertical direction at the k-1 / 2th moment at the spatial position (i, j+m); represents the normal stress in the vertical direction at the k-1 / 2th moment at the spatial position (i,j-m+1); represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1 / 2th moment; represents the normal stress in the horizontal direction at the k-1 / 2th moment at the spatial position (i, j); represents the horizontal component of velocity at the kth moment at the spatial position (i+m-1 / 2,j); represents the horizontal component of velocity at the kth moment at the spatial position (i-m+1 / 2,j); represents the vertical component of velocity at the kth moment at the spatial position (i,j+m-1 / 2); represents the vertical component of velocity at the kth moment at the spatial position (i,j-m+1 / 2); represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1 / 2th moment; represents the normal stress in the vertical direction at the k-1 / 2th moment at the spatial position (i, j); represents the shear stress at the k+1 / 2th moment at the spatial position (i+1 / 2, j+1 / 2); represents the shear stress at the k-1 / 2th moment at the spatial position (i+1 / 2, j+1 / 2); represents the horizontal velocity component at the kth moment at the spatial position (i+1 / 2,j+m); represents the horizontal velocity component at the kth moment at the spatial position (i+1 / 2,j-m+1); represents the vertical velocity component at the kth moment at the spatial position (i+m,j+1 / 2); represents the vertical velocity component at the kth moment at the spatial position (i-m+1,j+1 / 2); similarly, represents the horizontal velocity component at the spatial position (i, j) at the k+1th moment; represents the vertical velocity component at the k+1th moment at the spatial position (i, j); represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1th moment; represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1th moment; represents the shear stress at the spatial position (i, j) at the k+1th moment; represents the horizontal velocity component at the spatial position (i+1, j) at the k+1th moment; Represents the horizontal velocity component at the spatial position (i+2,j) at the k+1th moment.

5. The joint absorbing boundary method according to claim 1, characterized in that: The left boundary equation expression of the second-order Higdon absorbing boundary condition described in step (4) is: In formula (6), u x 、u z represent the wave field components in the x and z directions respectively, v p 、v s are the longitudinal wave velocity and the transverse wave velocity respectively; the Talyor formula is used to perform finite difference discretization on equation (6) and applied to the first-order velocity-stress elastic wave equation to obtain v x 、v z , τ xx , τ zz , τ xz The difference formats are: In formula (7), represents the horizontal component of velocity at the k+1th moment at the spatial position (i, j); represents the horizontal component of velocity at the klth moment at the spatial position (i, j); represents the vertical component of the velocity at the k+1th moment at the spatial position (i, j); represents the vertical component of velocity at the klth moment at the spatial position (i, j); represents the normal stress in the horizontal direction at the spatial position (i, j) at the k+1th moment; represents the normal stress in the horizontal direction at the klth moment at the spatial position (i, j); represents the normal stress in the vertical direction at the spatial position (i, j) at the k+1th moment; represents the normal stress in the vertical direction at the klth moment at the spatial position (i, j); represents the shear stress at the spatial position (i, j) at the k+1th moment; represents the shear stress at the klth moment at the spatial position (i, j); The coefficients in formula (7) can be obtained by the following identity: In formula (8), b is usually a constant between 0.3 and 0.5, and v p is the longitudinal wave velocity, Δt is the time sampling interval, and Δx is the spatial grid step size; the right boundary, upper boundary, and lower boundary can be derived by analogy with the above derivation process.

Images