A Joint Absorbing Boundary Method for Numerical Simulation of Low-Frequency Ship Seismic Waves

Abstract

This invention relates to a joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves, belonging to the field of ship seismic wave numerical simulation. The method first divides the simulation region into a central wavefield region and an artificial boundary region. Within the central wavefield region, the Ricker wavelet is used as the ship noise source, and the source disturbance is simulated based on the Ricker wavelet expression. The staggered-grid finite difference method is then used for ship seismic wave numerical simulation. The ship seismic wave wave equation is described by a two-dimensional first-order velocity-stress elastic wave equation. Within the central wavefield region, the elastic wave equation is discretized using staggered-grid finite difference, employing a second-order temporal and 2N-order spatial accuracy expansion. Within the artificial boundary region, an L-layer perfectly matched layer is added as an absorbing boundary, and a Gaussian attenuation function is used within the PML layer. A second-order Higdon absorbing boundary is used at the outer boundary. This invention can effectively improve the overall boundary absorption effect without increasing storage consumption.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation of ship seismic waves, and specifically relates to a joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves. Background Technology

[0002] Numerical simulation of ship seismic waves is a method that simplifies the real marine environment into a mathematical model and uses numerical calculations to simulate the propagation of ship seismic waves within the model. Through numerical simulations of various models, the propagation patterns of ship seismic waves in different seabed media can be analyzed, thus providing a theoretical basis for practical underwater long-range target detection based on ship seismic waves. In computer-based finite-difference numerical simulations of ship seismic waves, the computational domain needs to be manually defined due to the limited memory of the computer. However, simply truncating the boundaries will lead to strong reflected wave interference at the boundaries, significantly reducing simulation accuracy. Therefore, it is usually necessary to appropriately process the artificial boundaries to eliminate spurious reflections at the computational boundaries, thereby making the numerical simulation results more consistent with actual conditions.

[0003] One-way wave absorbing boundaries and perfectly matched layer (PML) absorbing boundaries are two main methods for handling artificial boundaries. One-way wave absorbing boundaries are boundary condition methods that use one-way wave equations to absorb incident wave energy within the boundary region. They mainly include Clayton-Enquist (CE) absorbing boundaries and Higdon absorbing boundaries. Compared to CE absorbing boundaries, Higdon absorbing boundaries generally have better absorption performance, and their differential calculations are more convenient. PML absorbing boundaries, on the other hand, absorb spurious reflections by adding a perfectly matched layer within the boundary region and using an equation containing attenuation functions within the layer. Theoretically, this method can absorb incident waves from various angles, and therefore has received widespread attention in recent years. Meng Luwen, Zhu Xufang, and others have used PML absorbing boundaries to achieve finite-difference numerical simulations of ship seismic waves, obtaining good simulation results while suppressing artificial boundary reflections.

[0004] However, existing PML absorbing boundaries still have the following shortcomings in practical applications: First, the exponential attenuation function selected by conventional PML absorbing boundaries increases too rapidly at the inner boundary, which leads to poor coupling between the central wavefield and the inner boundary wavefield, resulting in 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 still occur at the outer boundary of the PML. In summary, current conventional PML absorbing boundaries have not yet been able to solve the problem of false reflections at both the inner and outer boundaries of the boundary layer, and their overall absorption effect still needs to be further improved. Summary of the Invention

[0005] The technical problem this invention aims to solve is to provide a joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves. When performing numerical simulation of ship seismic waves based on the finite difference method, this invention first employs a Gaussian attenuation function within the PML boundary region to optimize the coupling effect between the central wavefield and the inner boundary, effectively reducing reflections from the inner boundary. Simultaneously, a second-order Higdon absorbing boundary is used at the outer boundary of the artificial boundary to further absorb residual incident waves. Compared to conventional PML boundaries, the joint absorbing boundary method employed in this invention can comprehensively suppress spurious reflections from both the inner and outer boundaries without increasing storage requirements, thereby significantly improving the accuracy of numerical simulation of ship seismic waves.

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

[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 wavefield region and an artificial boundary region. In the central wavefield region, 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 wave equation is described by the two-dimensional first-order velocity-stress elastic wave equation;

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

[0010]

[0011] In equation (1), t is time, ρ is the density of the model medium, w is the excitation sound source, x and z represent the horizontal and vertical directions in spatial coordinates, respectively, and v x v z τ represents the horizontal and vertical components of the velocity, respectively. xx τ zz τ xz These represent the three components of stress, with λ and μ being Lamé coefficients, calculated as follows:

[0012]

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

[0014] (2) Within the central wave field region, the elastic wave equation is discretized by staggered grid finite difference, and an expansion form with second-order time and 2Nth-order spatial accuracy is adopted.

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

[0016]

[0017] In equation (3), Δt is the time sampling interval, and b m denoted as spatial difference coefficient, k as the index of the time discrete point, i and j as the indexes of the sampling points in the x and z directions in the spatial coordinates, respectively, and Δx and Δz as the spatial grid step size.

[0018] (3) Within the artificial boundary region, an L-layer fully matched layer is added as an absorbing boundary. To optimize the coupling effect between the central wavefield and the inner boundary and further suppress inner boundary reflection, a Gaussian attenuation function is used within the PML layer. As a specific implementation method, the expression of the Gaussian attenuation function is as follows:

[0019]

[0020] In equation (4), d(h) is the numerical value of the decay function that varies with distance, and v p Let be the longitudinal wave velocity, h be the distance from the calculation point within the PML to the outermost boundary of the PML, and R be the theoretical reflection coefficient. The physical quantities in the two-dimensional first-order velocity-stress elastic wave equation are decomposed into horizontal and vertical components within the PML boundary layer, and then discretized using finite difference to obtain the following form:

[0021]

[0022] In equation (5), i and j are the discrete point indices in the horizontal and vertical directions of space, respectively, and k is the discrete point in time. x v z These represent the horizontal and vertical components of the velocity, respectively. σ xx σ zz These are the normal stresses in the horizontal and vertical directions, respectively, σ xz For shear stress, m represents the spatial index. This represents the horizontal component of the velocity at the k-th time point at spatial location (i+1 / 2,j). This represents the horizontal component of the velocity at the (k-1)th time at spatial location (i+1 / 2,j). This represents the normal stress in the horizontal direction at the (k-1 / 2)th time interval at spatial location (i+m,j). This represents the normal stress in the horizontal direction at the (k-1 / 2)th moment at the spatial location (i-m+1,j). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i+1 / 2, j+m-1 / 2). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i+1 / 2, j-m+1 / 2). This represents the vertical component of the velocity at the k-th moment at spatial location (i+1 / 2,j). This represents the vertical component of the velocity at the (k-1)th time at spatial location (i+1 / 2,j). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i+m-1 / 2,j+1 / 2). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i-m+1 / 2, j+1 / 2). This represents the normal stress in the vertical direction at the (k-1 / 2)th moment at spatial location (i,j+m). This represents the normal stress in the vertical direction at the (k-1 / 2)th moment at spatial location (i,j-m+1). This represents the normal stress in the horizontal direction at the (k+1 / 2)th time at spatial location (i,j). This represents the normal stress in the horizontal direction at the (k-1 / 2)th moment at spatial location (i,j). This represents the horizontal component of the velocity at the k-th time at spatial location (i+m-1 / 2,j). This represents the horizontal component of the velocity at the k-th time at spatial location (i-m+1 / 2,j). This represents the vertical component of the velocity at the k-th moment at spatial location (i,j+m-1 / 2). This represents the vertical component of the velocity at the k-th moment at spatial location (i, j-m+1 / 2). This represents the normal stress in the vertical direction at the (k+1 / 2)th moment at spatial location (i,j). This represents the normal stress in the vertical direction at the (k-1 / 2)th moment at spatial location (i,j). This represents the shear stress at the (k+1 / 2)th moment at spatial location (i+1 / 2, j+1 / 2). This represents the shear stress at the (k-1 / 2)th moment at spatial location (i+1 / 2, j+1 / 2). This represents the horizontal velocity component at the k-th time at spatial location (i+1 / 2, j+m). This represents the horizontal velocity component at the k-th time at spatial location (i+1 / 2, j-m+1). This represents the vertical velocity component at the k-th moment at spatial position (i+m,j+1 / 2). This represents the vertical velocity component at the k-th time point at spatial location (i-m+1, j+1 / 2). Similarly, This represents the horizontal velocity component at the (k+1)th time at spatial position (i,j). This represents the vertical velocity component at the (k+1)th time at spatial position (i,j). This represents the normal stress in the horizontal direction at the (k+1)th time at spatial location (i,j). This represents the normal stress in the vertical direction at the (k+1)th time at spatial location (i,j). This represents the shear stress at the (k+1)th time at spatial location (i,j). This represents the horizontal velocity component at the (k+1)th time at spatial location (i+1,j). This represents the horizontal velocity component at the (k+1)th time at spatial location (i+2,j).

[0023] (4) To absorb the residual incident wave at the outer boundary of the PML, a second-order Higdon absorbing boundary is used 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 equation (6), u x u z Let v represent the wave field components in the x and z directions, respectively. p v s Let v be the longitudinal wave velocity and the transverse wave velocity, respectively. Using the Talyor formula, equation (6) is discretized using finite difference and applied to the first-order velocity-stress elastic wave equation, yielding v. x v z τ xx τ zz τ xz The difference schemes are as follows:

[0026]

[0027] In equation (7), This represents the horizontal component of the velocity at the (k+1)th time at spatial location (i,j). This represents the horizontal component of the velocity at the kl-th time at spatial location (i,j). This represents the vertical component of the velocity at the (k+1)th time at spatial position (i,j). This represents the vertical component of the velocity at the kl-th time at spatial position (i,j). This represents the normal stress in the horizontal direction at the (k+1)th time at spatial location (i,j). This represents the normal stress in the horizontal direction at the kl-th moment at spatial location (i,j). This represents the normal stress in the vertical direction at the (k+1)th time at spatial location (i,j). This represents the normal stress in the vertical direction at the kl-th moment at spatial location (i,j). This represents the shear stress at the (k+1)th time at spatial location (i,j). This represents the shear stress at the kl-th moment at spatial location (i,j).

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

[0029]

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

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

[0032] Compared to the exponential decay function used in conventional PML absorbing boundaries, the Gaussian decay function exhibits a slower numerical change at the inner boundary. This results in better coupling between the central wavefield and the inner boundary, significantly reducing reflections at the inner boundary. Furthermore, the combined absorbing boundary proposed in this invention, based on the PML absorbing boundary employing a Gaussian decay function, further incorporates a second-order Higdon absorbing boundary to perform secondary absorption of the incompletely attenuated incident wave at the outer boundary of the PML. This effectively improves the overall boundary absorption effect without increasing storage consumption. Attached Figure Description

[0033] Figure 1 A schematic diagram of the area for numerical simulation of seismic waves from ships;

[0034] Figure 2 The following are wavefield snapshots at 500ms for numerical simulation of ship seismic waves using different artificial absorbing boundaries: (a) Wavefield snapshot of an artificial absorbing boundary without any processing, (b) Wavefield snapshot of a conventional PML absorbing boundary using an exponential decay function, and (c) Wavefield snapshot of a combined absorbing boundary.

[0035] Figure 3 Comparison of single-track seismic records at (425m, 975m) under different artificial absorption boundary conditions. Detailed Implementation

[0036] The technical solution of the present invention will be further explained below with reference to examples and accompanying drawings, but the scope of protection of the present invention is not limited in any way.

[0037] Example 1

[0038] This invention employs a joint absorbing boundary method for numerical simulation of low-frequency ship seismic waves. During implementation, the simulation region is divided into a central wavefield region and an artificial boundary region. A joint absorbing boundary is constructed within the artificial boundary region, thereby effectively improving the boundary absorption effect of the numerical simulation. A schematic diagram of the ship seismic wave numerical simulation region is shown below. Figure 1 As shown, its specific implementation method is as follows:

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

[0040]

[0041] In equation (1), t is time, ρ is the density of the model medium, w is the excitation sound source, x and z represent the horizontal and vertical directions in spatial coordinates, respectively, and v x v z τ represents the horizontal and vertical components of the velocity, respectively. xx τ zz τ xz These represent the three components of stress, with λ and μ being Lamé coefficients, calculated as follows:

[0042]

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

[0044] (2) Within the central wave field region, the elastic wave equation is discretized by staggered grid finite difference, and an expansion form with second-order time and 2Nth-order spatial accuracy is adopted.

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

[0046]

[0047] In equation (3), Δt is the time sampling interval, and b mdenoted as spatial difference coefficient, k as the index of the time discrete point, i and j as the indexes of the sampling points in the x and z directions in the spatial coordinates, respectively, and Δx and Δz as the spatial grid step size.

[0048] (3) Within the artificial boundary region, an L-layer fully matched layer is added as an absorbing boundary; in this implementation, L = 40. To optimize the coupling effect between the central wavefield and the inner boundary and further suppress inner boundary reflection, a Gaussian attenuation function is used within the PML layer. As a specific implementation method, the expression for the Gaussian attenuation function is as follows:

[0049]

[0050] In equation (4), d(h) is the numerical value of the decay function that varies with distance, and v p Let be the longitudinal wave velocity, h be the distance from the calculation point within the PML to the outermost boundary of the PML, and R be the theoretical reflection coefficient. The physical quantities in the two-dimensional first-order velocity-stress elastic wave equation are decomposed into horizontal and vertical components within the PML boundary layer, and then discretized using finite difference to obtain the following form:

[0051]

[0052] In equation (5), i and j are the discrete point indices in the horizontal and vertical directions of space, respectively, and k is the discrete point in time. x v z These represent the horizontal and vertical components of the velocity, respectively. σ xx σ zz These are the normal stresses in the horizontal and vertical directions, respectively, σ xz For shear stress, m represents the spatial index. This represents the horizontal component of the velocity at the k-th time point at spatial location (i+1 / 2,j). This represents the horizontal component of the velocity at the (k-1)th time at spatial location (i+1 / 2,j). This represents the normal stress in the horizontal direction at the (k-1 / 2)th time interval at spatial location (i+m,j). This represents the normal stress in the horizontal direction at the (k-1 / 2)th moment at the spatial location (i-m+1,j). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i+1 / 2, j+m-1 / 2). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i+1 / 2, j-m+1 / 2). This represents the vertical component of the velocity at the k-th moment at spatial location (i+1 / 2,j). This represents the vertical component of the velocity at the (k-1)th time at spatial location (i+1 / 2,j). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i+m-1 / 2,j+1 / 2). This represents the shear stress at the (k-1 / 2)th moment at the spatial location (i-m+1 / 2, j+1 / 2). This represents the normal stress in the vertical direction at the (k-1 / 2)th moment at spatial location (i,j+m). This represents the normal stress in the vertical direction at the (k-1 / 2)th moment at spatial location (i,j-m+1). This represents the normal stress in the horizontal direction at the (k+1 / 2)th time at spatial location (i,j). This represents the normal stress in the horizontal direction at the (k-1 / 2)th moment at spatial location (i,j). This represents the horizontal component of the velocity at the k-th time at spatial location (i+m-1 / 2,j). This represents the horizontal component of the velocity at the k-th time at spatial location (i-m+1 / 2,j). This represents the vertical component of the velocity at the k-th moment at spatial location (i,j+m-1 / 2). This represents the vertical component of the velocity at the k-th moment at spatial location (i, j-m+1 / 2). This represents the normal stress in the vertical direction at the (k+1 / 2)th moment at spatial location (i,j). This represents the normal stress in the vertical direction at the (k-1 / 2)th moment at spatial location (i,j). This represents the shear stress at the (k+1 / 2)th moment at spatial location (i+1 / 2, j+1 / 2). This represents the shear stress at the (k-1 / 2)th moment at spatial location (i+1 / 2, j+1 / 2). This represents the horizontal velocity component at the k-th time at spatial location (i+1 / 2, j+m). This represents the horizontal velocity component at the k-th time at spatial location (i+1 / 2, j-m+1). This represents the vertical velocity component at the k-th moment at spatial position (i+m,j+1 / 2). This represents the vertical velocity component at the k-th time point at spatial location (i-m+1, j+1 / 2). Similarly, This represents the horizontal velocity component at the (k+1)th time at spatial position (i,j). This represents the vertical velocity component at the (k+1)th time at spatial position (i,j). This represents the normal stress in the horizontal direction at the (k+1)th time at spatial location (i,j). This represents the normal stress in the vertical direction at the (k+1)th time at spatial location (i,j). This represents the shear stress at the (k+1)th time at spatial location (i,j). This represents the horizontal velocity component at the (k+1)th time at spatial location (i+1,j). This represents the horizontal velocity component at the (k+1)th time at spatial location (i+2,j).

[0053] (4) To absorb the residual incident wave at the outer boundary of the PML, a second-order Higdon absorbing boundary is used 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 equation (6), u x u z Let v represent the wave field components in the x and z directions, respectively. p v s Let v be the longitudinal wave velocity and the transverse wave velocity, respectively. Using the Talyor formula, equation (6) is discretized using finite difference and applied to the first-order velocity-stress elastic wave equation, yielding v. x v z τ xx τ zz τ xz The difference schemes are as follows:

[0056]

[0057] In equation (7), This represents the horizontal component of the velocity at the (k+1)th time at spatial location (i,j). This represents the horizontal component of the velocity at the kl-th time at spatial location (i,j). This represents the vertical component of the velocity at the (k+1)th time at spatial position (i,j). This represents the vertical component of the velocity at the kl-th time at spatial position (i,j). This represents the normal stress in the horizontal direction at the (k+1)th time at spatial location (i,j). This represents the normal stress in the horizontal direction at the kl-th moment at spatial location (i,j). This represents the normal stress in the vertical direction at the (k+1)th time at spatial location (i,j). This represents the normal stress in the vertical direction at the kl-th moment at spatial location (i,j). This represents the shear stress at the (k+1)th time at spatial location (i,j). This represents the shear stress at the kl-th moment at spatial location (i,j).

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

[0059]

[0060] In equation (8), b is usually a constant between 0.3 and 0.5, and v p Let Δt be the longitudinal wave velocity, Δt be the time sampling interval, and Δx be the spatial grid step size. 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 500 ms for numerical simulations of ship seismic waves using different artificial absorbing boundaries: (a) is the wavefield snapshot of the artificial absorbing boundary without any processing; (b) is the wavefield snapshot of the conventional PML absorbing boundary using an exponential decay function; and (c) is the wavefield snapshot of the combined absorbing boundary. Figure 2 It can be seen that, compared to untreated artificial absorbing boundaries, conventional PML absorbing boundaries have a better suppression effect on boundary reflections, but relatively obvious residual boundary reflections can still be seen at the inner and outer boundaries of the PML. When a combined absorbing boundary is used to realize the numerical simulation of ship seismic waves, the reflected waves at the inner and outer boundaries of the artificial boundary region are effectively absorbed, and the overall simulation accuracy is significantly improved. Simultaneously, combined with... Figure 3 The comparison of the provided single-channel seismic record images also shows that, when using the combined absorbing boundary, the amplitudes of the inner and outer reflections of the artificial absorbing boundary are significantly smaller than those of the unprocessed artificial absorbing boundary and the conventional PML absorbing boundary. In summary, the combined absorbing boundary method proposed in this invention for numerical simulation of low-frequency ship seismic waves effectively improves the boundary absorption effect, significantly enhances the accuracy of numerical simulation of ship seismic waves, and provides effective technical support for further improving the accuracy of underwater target detection.

Claims

1. A combined absorbing boundary method applied to numerical simulation of low-frequency seismic waves for ships, characterized in that, The method comprises the following steps: (1) in the numerical simulation process, first, the simulation area is divided into a central wave field area and an artificial boundary area; in the central wave field area, a Ricker wavelet is taken as a ship noise source, a source disturbance is simulated based on a Ricker wavelet expression, and an interlaced grid finite difference method is used for ship seismic wave numerical simulation; a ship seismic wave wave equation is described by a two-dimensional first-order velocity-stress elastic wave equation; (2) in the central wave field area, an interlaced grid finite difference is used for discrete dispersion of the elastic wave equation, and a time 2-order and space 2N-order accuracy expansion form is used; (3) in the artificial boundary area, L layers of completely matched layers are added as absorbing boundaries, and a Gaussian type attenuation function is used in the PML layer; the expression of the Gaussian type attenuation function is as follows: In formula (4), d(h) is a value of a distance-varying attenuation function, v p is a longitudinal wave velocity, h is a distance from a calculation point in the PML to an outermost boundary of the PML, and R is a theoretical reflection coefficient; (4) to absorb residual incident waves at the outer boundary of the PML, a second-order Higdon absorbing boundary is used at the outer boundary; 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 time, p is the density of the model medium, w is the excited sound source, x and z represent the horizontal and vertical directions in the spatial coordinates, respectively, v x , v z represent the horizontal and vertical components of the velocity, respectively, τ xx , τ zz , τ xz represent three components of stress, respectively, and λ and μ are Lame coefficients, the calculation formulas of which are as follows: In formula (2), v p , v s are respectively longitudinal wave velocity, transverse wave velocity; The step (3) decomposes each physical quantity in the two-dimensional first-order velocity-stress elastic wave equation into horizontal component and vertical component in the PML boundary layer, and obtains the following form by finite difference discretization: In formula (5), i, j are discrete point serial numbers in horizontal and vertical directions of space respectively, and k is a time discrete point serial number; v x , v z are horizontal and vertical components of velocity respectively; σ xx , σ zz are normal stresses in horizontal and vertical directions respectively, σ xz is a shear stress, and m represents a spatial serial number; represents a horizontal component of velocity at the kth moment at the spatial position (i+1 / 2, j); represents a horizontal component of velocity at the (k-1)th moment at the spatial position (i+1 / 2, j); represents a normal stress in horizontal direction at the (k-1 / 2)th moment at the spatial position (i+m, j); represents a normal stress in horizontal direction at the (k-1 / 2)th moment at the spatial position (i-m+1, j); represents a shear stress at the (k-1 / 2)th moment at the spatial position (i+1 / 2, j+m-1 / 2); represents a shear stress at the (k-1 / 2)th moment at the spatial position (i+1 / 2, j-m+1 / 2); represents a vertical component of velocity at the kth moment at the spatial position (i+1 / 2, j); represents a vertical component of velocity at the (k-1)th moment at the spatial position (i+1 / 2, j); represents a shear stress at the (k-1 / 2)th moment at the spatial position (i+m-1 / 2, j+1 / 2); represents a shear stress at the (k-1 / 2)th moment at the spatial position (i-m+1 / 2, j+1 / 2); represents a normal stress in vertical direction at the (k-1 / 2)th moment at the spatial position (i, j+m); represents a normal stress in vertical direction at the (k-1 / 2)th moment at the spatial position (i, j-m+1); represents a normal stress in horizontal direction at the (k+1 / 2)th moment at the spatial position (i, j); represents a normal stress in horizontal direction at the (k-1 / 2)th moment at the spatial position (i, j); represents a horizontal component of velocity at the kth moment at the spatial position (i+m-1 / 2, j); represents a horizontal component of velocity at the kth moment at the spatial position (i-m+1 / 2, j); represents a vertical component of velocity at the kth moment at the spatial position (i, j+m-1 / 2); Vz(k) = V(k, i, j - m + 1 / 2) represents the vertical component of velocity at spatial location (i, j - m + 1 / 2) at time k; σxx(k+1 / 2) = σ(k+1 / 2, i, j) represents the normal stress in the horizontal direction at spatial location (i, j) at time k + 1 / 2; σxx(k-1 / 2) = σ(k-1 / 2, i, j) represents the normal stress in the horizontal direction at spatial location (i, j) at time k - 1 / 2; τxy(k+1 / 2) = τ(k+1 / 2, i+1 / 2, j+1 / 2) represents the shear stress at spatial location (i+1 / 2, j+1 / 2) at time k + 1 / 2; τxy(k-1 / 2) = τ(k-1 / 2, i+1 / 2, j+1 / 2) represents the shear stress at spatial location (i+1 / 2, j+1 / 2) at time k - 1 / 2; Vx(k) = V(k, i+1 / 2, j+m) represents the horizontal component of velocity at spatial location (i+1 / 2, j+m) at time k; Vx(k) = V(k, i+1 / 2, j-m+1) represents the horizontal component of velocity at spatial location (i+1 / 2, j-m+1) at time k; Vz(k) = V(k, i+m, j+1 / 2) represents the vertical component of velocity at spatial location (i+m, j+1 / 2) at time k; Vz(k) = V(k, i-m+1, j+1 / 2) represents the vertical component of velocity at spatial location (i-m+1, j+1 / 2) at time k; Vx(k+1) = V(k+1, i, j) represents the horizontal component of velocity at spatial location (i, j) at time k + 1; Vz(k+1) = V(k+1, i, j) represents the vertical component of velocity at spatial location (i, j) at time k + 1; σxx(k+1) = σ(k+1, i, j) represents the normal stress in the horizontal direction at spatial location (i, j) at time k + 1; σzz(k+1) = σ(k+1, i, j) represents the normal stress in the vertical direction at spatial location (i, j) at time k + 1; τxy(k+1) = τ(k+1, i, j) represents the shear stress at spatial location (i, j) at time k + 1; Vx(k+1) = V(k+1, i+1, j) represents the horizontal component of velocity at spatial location (i+1, j) at time k + 1; Vx(k+1) = V(k+1, i+2, j) represents the horizontal component of velocity at spatial location (i+2, j) at time k + 1.

2. The federated absorbing boundary method of claim 1, wherein, The step (2) of discretely representing the elastic wave equation by staggered grid finite difference is expressed as: In formula (3), Δt is a time sampling interval, b m is a spatial difference coefficient, k is a time discrete point sequence number, i and j represent sampling point sequence numbers in x and z directions of the spatial coordinates respectively, and Δx and Δz represent spatial grid steps respectively.

3. The federated absorbing boundary method of claim 1, wherein, The left boundary equation expression of the second-order Higdon absorbing boundary condition in step (4) is as follows: In formula (6), u x , u z respectively represent wave field components in x and z directions, v p , v s respectively represent longitudinal wave velocity and transverse wave velocity; finite difference discretization is performed on formula (6) by using Taylor formula, and it is applied to a first-order velocity-stress elastic wave equation, so that difference formats of v x , v z , τ xx , τ zz , τ xz are respectively obtained: in formula (7), denotes the horizontal velocity component at spatial position (i,j) at time k+1; denotes the horizontal velocity component at spatial position (i,j) at time k-l; denotes the vertical velocity component at spatial position (i,j) at time k+1; denotes the vertical velocity component at spatial position (i,j) at time k-l; denotes the normal stress in horizontal direction at spatial position (i,j) at time k+1; denotes the normal stress in horizontal direction at spatial position (i,j) at time k-l; denotes the normal stress in vertical direction at spatial position (i,j) at time k+1; denotes the normal stress in vertical direction at spatial position (i,j) at time k-l; denotes the shear stress at spatial position (i,j) at time k+1; denotes the shear stress at spatial position (i,j) at time k-l; The coefficients in formula (7) can be obtained through the following identity: In formula (8), b is usually a constant between 0.3 and 0.5, v p is the longitudinal wave velocity, Δt is the time sampling interval, and Δx is the spatial grid step; the right boundary, the upper boundary, and the lower boundary can be derived by analogy with the above derivation process.

Images