Fast calculation method for seafloor seismic wave field induced by axisymmetric volume sound source in water

Abstract

The present invention relates to a method for rapidly calculating the submarine seismic wave field induced by an axisymmetric volume sound source in water, comprising the steps of: constructing an equivalent sound source intensity - sound pressure relationship equation; constructing a series form of the Green's function; calculating the reflection coefficient under an elastic seabed; constructing a vibration velocity - equivalent sound source relationship equation; calculating the equivalent sound source intensity; constructing an axisymmetric finite element coupling calculation model for the equivalent source intensity; constructing an axisymmetric finite element coordinate system in cylindrical coordinates; establishing an axisymmetric finite element corresponding to the structure; establishing an axisymmetric finite element corresponding to the seawater; establishing an axisymmetric finite element corresponding to the acoustic - solid coupling boundary surface; establishing a matching infinite element for truncating the infinite outer domain; establishing an acoustic - vibration coupling dynamic equation under axisymmetric conditions; constructing an axisymmetric finite element coupling calculation model for the equivalent source intensity; and calculating the submarine seismic wave induced by an axisymmetric body target in water. The present invention realizes rapid and high - precision calculation of the submarine seismic wave field induced by an axisymmetric volume sound source in water.

Description

Technical Field

[0001] The present invention relates to the technical field of underwater seismic wave field calculation, and particularly to a fast calculation method for underwater seismic wave field induced by an axisymmetric volume sound source in water. Background Art

[0002] In addition to propagating in water, the radiated noise of an underwater target can also propagate through the water body to the interface between the sea water and the seabed, that is, the acoustic-solid coupling surface, which is the so-called underwater seismic wave, especially obvious in shallow waters. Different from the acoustic wave in sea water which only contains longitudinal wave components, the underwater seismic wave contains different forms of energy components such as longitudinal waves, transverse waves, lateral waves, and Scholte waves. They not only carry the effective information of the underwater target, but are also little affected by the underwater acoustic channel and have the ability of remote detection, and have received extensive attention from various countries in recent years.

[0003] Currently, the prediction of underwater seismic waves induced by underwater target sound sources at home and abroad mainly focuses on point sound sources and pulsed sound sources, and there is relatively little attention paid to volume sound sources, and there is no attention to the fast prediction of underwater seismic waves induced by axisymmetric volume sound sources in water, which is the premise for quickly and accurately locating axisymmetric volume sound sources in water;

[0004] However, there is no any applied research on the above technical problems in the existing technology, which belongs to a technical vacuum. Summary of the Invention

[0005] In view of the above problems, the present invention provides a fast calculation method for underwater seismic wave field induced by an axisymmetric volume sound source in water, and its purpose is to realize fast and high-precision calculation of the underwater seismic wave field induced by an axisymmetric volume sound source in water.

[0006] To solve the above problems, the technical solution provided by the present invention is as follows:

[0007] A fast calculation method for underwater seismic wave field induced by an axisymmetric volume sound source in water, comprising the following steps:

[0008] S100. Calculate the equivalent sound source intensity; specifically, it comprises the following steps:

[0009] S110. Construct an equation for the relationship between equivalent sound source intensity and sound pressure;

[0010] S120. Construct a series form of the Green's function; then calculate the reflection coefficient under the elastic seabed;

[0011] S130. Construct an equation for the relationship between vibration velocity and equivalent sound source; then use the equation for the relationship between vibration velocity and equivalent sound source to calculate the equivalent sound source intensity;

[0012] S200. Construct an axisymmetric finite element coupling calculation model for the equivalent source intensity; specifically, it comprises the following steps:

[0013] S210. Construct an axisymmetric finite element coordinate system in cylindrical coordinates;

[0014] S220. Establish an axisymmetric finite element corresponding to the structure;

[0015] S230. Establish an axisymmetric finite element corresponding to the seawater;

[0016] S240. Establish an axisymmetric finite element corresponding to the acoustic-solid coupling boundary surface;

[0017] S250. Establish a matching infinite element for truncating the infinite outer domain;

[0018] S260. Establish the acoustic-vibration coupling dynamic equation under axisymmetric conditions; and establish the axisymmetric finite element coupling calculation model for the equivalent source strength according to the acoustic-vibration coupling dynamic equation under the axisymmetric conditions;

[0019] S300. Calculate the submarine seismic waves induced by the axisymmetric body target according to the equivalent sound source strength and the axisymmetric finite element coupling calculation model for the equivalent source strength.

[0020] Preferably, constructing the equivalent sound source strength-pressure relationship equation in S110 includes the following steps:

[0021] S111. Establish a sound radiation calculation model of the structure under the shallow sea channel; the sound radiation calculation model of the structure under the shallow sea channel includes the following parts:

[0022] The area where the structure is located under the shallow sea channel, the area where the shallow sea fluid medium is located, the area where the solid seabed is located, the coupling boundary between the elastic structure and the fluid, the coupling boundary between the seabed and the seawater, the outer normal direction of the structure, the outer normal direction of the seabed, the virtual surface where the equivalent sound source is located, the coordinate origin, the position vector coordinates of the equivalent sound source, the relative position vector coordinates between the equivalent sound source and the point where the sound field is to be obtained, and the seawater depth;

[0023] S112. Use the sound radiation calculation model of the structure under the shallow sea channel to construct the equivalent sound source strength-pressure relationship equation;

[0024] The equivalent sound source strength-pressure relationship equation is expressed as follows:

[0025]

[0026] Where: p(r S ) is the sound pressure on the surface of the structure under the shallow sea channel; r S is the surface of the structure under the shallow sea channel; r 0 = is the position vector coordinates of the equivalent sound source, denoted as r 0 =(x0 , y 0 , z 0 ), located above the virtual plane where the equivalent sound source is located; k is the acoustic wave number, expressed by the following formula:

[0027] k = ω / c f

[0028] where: ω is the angular frequency; c f is the acoustic wave propagation speed;

[0029] ρ f is the density of the fluid medium; i is the unit imaginary number; c is the acoustic wave speed; s(r 0 ) is the equivalent sound source intensity; G k (r S , r 0 ) is the Green's function of the waveguide space, expressed by the following formula:

[0030]

[0031] where: δ(r - r 0 ) is the Delta function, is the Laplace operator;

[0032] M is the number of equivalent sources; m is the number marking symbol of the equivalent source, and m ∈ [0, M];

[0033] Preferably, the series form of the Green's function is expressed by the following formula:

[0034]

[0035] where: l is the marking symbol, and l ∈ [0, ∞); R l1 , R l2 , R l3 and R l4 The calculation formulas are as follows; e is the natural logarithm;

[0036] In the formula:

[0037]

[0038] where: x, y, z are the relative position vector coordinates between the equivalent sound source and the point to be solved in the sound field, and are uniformly expressed as r = (x, y, z); h is the depth of the seawater; h 1 is the distance from the equivalent sound source to the seawater surface; h 2 is the distance from the equivalent sound source to the seabed; γ 1 = -1; γ is the reflection coefficient under the elastic seabed, expressed by the following formula:

[0039]

[0040] Wherein: Z is an intermediate variable; ρ s2 is the seabed density; θ is both the incident wave angle and the reflected wave angle; β is the refracted longitudinal wave angle; α is the refracted shear wave angle; θ, β, and α satisfy the following relationships:

[0041]

[0042] Wherein: c p2 is the shear wave velocity of the seabed medium; c s2 is the longitudinal wave velocity of the seabed medium.

[0043] Preferably, the vibration velocity - equivalent sound source relationship equation in S130 is expressed as follows:

[0044]

[0045] Wherein: u(r S ) represents the normal displacement of the surface of the structure; is the Laplace operator; n S is the outer normal direction of the structure, is the normal gradient operator.

[0046] Preferably, in the axisymmetric finite element coordinate system in cylindrical coordinates in S210, the cylindrical coordinate system is expressed as (r, θ, z); the overall seawater and seabed coupled solution domain takes r = 0 as the axis of symmetry; the equivalent sound source is located on r = 0.

[0047] Preferably, in the axisymmetric finite element corresponding to the structure in S220, the partial differential equation of the structure in the cylindrical coordinate system on the axis of symmetry is expressed as follows:

[0048]

[0049] Wherein: σ rr is the normal stress of the axisymmetric element in the r direction; σ zz is the normal stress of the axisymmetric element in the z direction σ θθ is the normal stress of the axisymmetric element in the θ direction; τ rz is the shear stress; u is the displacement along the r direction; w is the displacement along the z direction; the axisymmetric element includes an acoustic axisymmetric element and a structural axisymmetric element;

[0050] The strain matrix in the cylindrical coordinate system is expressed as follows:

[0051]

[0052] Wherein: ε rr is the normal strain of the axisymmetric element in the r direction; ε θθis the normal strain of the axisymmetric element in the z direction; ε zz is the normal strain of the axisymmetric element in the θ direction; γ rz is the shear strain;

[0053] The corresponding linear elastic stress matrix is expressed as follows:

[0054] σ = [σ rr σ θθ σ zz τ rz T = Dε

[0055] where: D is the matrix coefficient, expressed as follows:

[0056]

[0057] where: E is the elastic modulus; v is the Poisson's ratio;

[0058] The displacement on each axisymmetric element is expressed as follows:

[0059]

[0060] where: q is the nodal displacement vector of the axisymmetric element, and q = [u 1 w 1 … u m w m T ; N s is the matrix composed of the shape function N(ξ, η), expressed as follows:

[0061]

[0062] where: N(ξ, η) is the shape function; ξ is the interpolation function in the ξ direction of the local coordinate system; η is the interpolation function in the η direction of the local coordinate system;

[0063] The strain matrix of the axisymmetric element is expressed as follows:

[0064] ε = Bq

[0065] where: B is the stress-displacement matrix, and B = [B 1 , B 2 , …, B m ; B m is expressed as follows:

[0066]

[0067] The stress matrix corresponding to the axisymmetric element is expressed as follows:

[0068] σ = DBq ​​

[0069] The stiffness matrix corresponding to the axisymmetric element is expressed by the following formula:

[0070]

[0071] The mass matrix corresponding to the axisymmetric element is expressed by the following formula:

[0072]

[0073] Where: J is the Jacobian matrix for converting from the global coordinates (r, z) to the local coordinates (ξ, η), and is expressed by the following formula:

[0074]

[0075] det(J) is the determinant corresponding to the matrix J;

[0076] The elastic seabed response matrix is obtained by calculating the matrix and the matrix respectively using Gauss integration, and then combining each axisymmetric element according to each discrete node.

[0077] Preferably, each axisymmetric element corresponds to 9 calculation nodes; in the matrix N s formed by the shape functions N(ξ, η), m = 9, and each shape function corresponds to one of the said calculation nodes; each shape function N(ξ, η) is expressed as follows:

[0078]

[0079]

[0080]

[0081] Preferably, in the axisymmetric finite element corresponding to the seawater in S230:

[0082] When the equivalent sound source is located in the seawater, the control equation of the acoustic axisymmetric element is expressed by the following formula:

[0083]

[0084] Where: p is the sound pressure on each acoustic axisymmetric element, and is expressed by the following formula:

[0085] p = N f p

[0086] Where: N f is the matrix formed by the shape functions N(ξ, η), and N f = [N 1 (ξ, η) N 2(ξ,η) … N m (ξ,η)]; p is the acoustic pressure vector of the nodes of the axisymmetric acoustic element, and p = [p 1 p 2 … p m T ;

[0087] The stiffness matrix corresponding to the axisymmetric acoustic element is expressed by the following formula:

[0088]

[0089] The mass matrix corresponding to the axisymmetric acoustic element is expressed by the following formula:

[0090]

[0091] Where:

[0092] The seawater corresponding matrix is obtained by calculating the matrix and using Guass integration, and then combining each structural axisymmetric element according to each discrete node;

[0093] In the axisymmetric finite element corresponding to the acoustic-solid coupling boundary surface in S240:

[0094] On the acoustic-solid interface, the stress continuity condition is expressed by the following formula:

[0095] p| Γ =(σ·n)| Γ

[0096] On the acoustic-solid interface, the displacement continuity condition is expressed by the following formula:

[0097]

[0098] Where: Γ is the acoustic-solid interface; p| Γ is the acoustic pressure at the interface of the axisymmetric acoustic element; (σ·n)| Γ is the normal stress at the interface of the structural axisymmetric element; is the projection of the acoustic pressure gradient of the axisymmetric acoustic element in the normal direction of the interface; (u T ·n)| Γ is the normal displacement of the structural axisymmetric element at the interface;

[0099] The coupling matrix is expressed by the following formula:

[0100]

[0101] Where: n is the outer normal direction of the acoustic-solid coupling interface structure; ​

[0102] In the matching infinite elements of the truncated infinite domain in S250:

[0103] The mapping relationship of the coordinate positions of the matching infinite elements from the global coordinate system to the local coordinate system is expressed by the following formula:

[0104] r = M 1 (ξ)r 1 + M 2 (ξ)r 2

[0105] Where: M 1 (ξ) and M 2 (ξ) are both shape functions and are expressed by the following formula:

[0106]

[0107] The shape function M 1 (ξ) and the shape function M 2 (ξ) match the matching infinite elements in the global coordinate system and are expressed by the following formula:

[0108]

[0109] Where: the position of the pole is the position of the axis of symmetry; r 1 is the connection point between the matching infinite element and the finite element, and r 3 is the place where the matching infinite element extends to infinity.

[0110] Preferably, S260 specifically includes the following steps:

[0111] S261. Assemble all the structural axisymmetric elements into a structural axisymmetric element stiffness matrix and a structural axisymmetric element mass matrix; assemble all the acoustic axisymmetric elements into an acoustic axisymmetric element stiffness matrix and an acoustic axisymmetric element mass matrix; assemble all the fluid-structure interaction elements into a fluid-structure interaction element stiffness matrix and a fluid-structure interaction element mass matrix; assemble all the matching infinite elements into a matching infinite element stiffness matrix and a matching infinite element mass matrix;

[0112] S262. According to the corresponding relationship between the nodes, assemble the structural axisymmetric element stiffness matrix, the acoustic axisymmetric element stiffness matrix, the fluid-structure interaction element stiffness matrix, and the matching infinite element stiffness matrix into a coupled stiffness matrix;

[0113] S263. According to the corresponding relationship between the nodes, assemble the structural axisymmetric element mass matrix, the acoustic axisymmetric element mass matrix, the fluid-structure interaction element mass matrix, and the matching infinite element mass matrix into a coupled mass matrix;

[0114] S264. According to the pressure release boundary condition on the sea surface and the displacement boundary condition related to the axis of symmetry, the coupled stiffness matrix and the coupled mass matrix are respectively trimmed to obtain the acoustic-vibration coupling dynamic equation under the axisymmetric condition;

[0115] S265. According to the acoustic-vibration coupling dynamic equation under the axisymmetric condition, the sea water and the seabed are respectively truncated by using the infinite element for the infinite space, so as to establish the axisymmetric finite element coupling calculation model for the equivalent source strength.

[0116] Preferably, S300 specifically includes the following steps:

[0117] S310. Calculate the intensity corresponding to the equivalent sound source by using the mirror image method and the equivalent source method simultaneously;

[0118] S320. Substitute the position and intensity of the equivalent sound source into the axisymmetric finite element coupling calculation model for the equivalent source strength, calculate the sound pressure corresponding to the seabed surface in the near field and the sound pressure corresponding to the seabed surface in the far field respectively, and then superimpose the sound pressure corresponding to the seabed surface in the near field and the sound pressure corresponding to the seabed surface in the far field, and finally obtain the submarine seismic wave induced by the axisymmetric body target.

[0119] Compared with the prior art, the present invention has the following advantages:

[0120] 1. Since the present invention adopts the method combining the equivalent source and the finite element, based on the wave superposition principle, the radiation sound field of the volume sound source is equivalent to the superposition of the radiation sound fields of a number of equivalent point sources, and the finite element equation considering the acoustic-vibration coupling effect is further considered, so that the submarine seismic wave induced by the axisymmetric volume in water can be calculated quickly;

[0121] 2. Since the present invention calculates the sound pressure / stress at the sea water-seabed interface generated by each equivalent source under different distribution positions and intensity conditions, and linearly superimposes them according to the linear superposition principle, the submarine seismic wave induced by the axisymmetric volume in water can be calculated with high precision;

[0122] 3. Since compared with the volume sound source, the point sound source can be regarded as a force and directly applied to the corresponding water body unit, when applying the method of the present invention to calculate the submarine seismic wave, it is not necessary to consider the interface coupling effect between the volume sound source and the sea water unit, so that the modeling method is simple and the calculation amount is small;

[0123] 4. Since the present invention adopts the method combining the equivalent source and the finite element, as long as the axisymmetric volume target is in water, the seismic wave induced by it can be calculated, whether the target is partially or fully located in water, and whether it is in shallow sea or deep sea, so that the application range is very wide. Description of the Drawings

[0124] Figure 1 Schematic diagram of the method flow for a specific embodiment of the present invention;

[0125] Figure 2 Schematic diagram of the local calculation model of structural acoustic radiation in a shallow - sea channel for a specific embodiment of the present invention;

[0126] Figure 3 Schematic diagram of the reflection and refraction of plane acoustic waves on the elastic seabed for a specific embodiment of the present invention;

[0127] Figure 4 Schematic diagram of the axisymmetric element and the corresponding stress tensor in the cylindrical coordinate system for a specific embodiment of the present invention;

[0128] Figure 5 Schematic diagram of the unit corresponding to the second - order shape function for a specific embodiment of the present invention;

[0129] Figure 6 Schematic diagram of the node - mapped infinite element for a specific embodiment of the present invention;

[0130] Figure 7 Schematic diagram of the 2D infinite element in the local coordinate system for a specific embodiment of the present invention;

[0131] Figure 8 Schematic diagram of the axisymmetric finite - element calculation model for a specific embodiment of the present invention;

[0132] Figure 9 Schematic diagram of the spherical sound source and its coordinate system under the shallow sea for a specific embodiment of the present invention;

[0133] Figure 10 Schematic diagram of the distribution of field points and equivalent point sources for a specific embodiment of the present invention;

[0134] Figure 11-a Schematic diagram of the sound - pressure propagation loss at the sea - seabed interface excited by a pulsating sphere with a hard seabed at an analysis frequency of 30 Hz for a specific embodiment of the present invention;

[0135] Figure 11-b Schematic diagram of the sound - pressure propagation loss at the sea - seabed interface excited by a pulsating sphere with a hard seabed at an analysis frequency of 50 Hz for a specific embodiment of the present invention;

[0136] Figure 11-c Schematic diagram of the sound - pressure propagation loss at the sea - seabed interface excited by a pulsating sphere with a soft seabed at an analysis frequency of 30 Hz for a specific embodiment of the present invention;

[0137] Figure 11-d Schematic diagram of the sound - pressure propagation loss at the sea - seabed interface excited by a pulsating sphere with a soft seabed at an analysis frequency of 50 Hz for a specific embodiment of the present invention;

[0138] Figure 12-a Schematic diagram of the sound pressure propagation loss at the seawater-seabed interface excited by a laterally vibrating sphere on a hard seabed at an analysis frequency of 30 Hz in a specific embodiment of the present invention;

[0139] Figure 12-b Schematic diagram of the sound pressure propagation loss at the seawater-seabed interface excited by a laterally vibrating sphere on a hard seabed at an analysis frequency of 50 Hz in a specific embodiment of the present invention;

[0140] Figure 12-c Schematic diagram of the sound pressure propagation loss at the seawater-seabed interface excited by a laterally vibrating sphere on a soft seabed at an analysis frequency of 30 Hz in a specific embodiment of the present invention;

[0141] Figure 12-d Schematic diagram of the sound pressure propagation loss at the seawater-seabed interface excited by a laterally vibrating sphere on a soft seabed at an analysis frequency of 50 Hz in a specific embodiment of the present invention. Detailed implementation manners

[0142] The following further clarifies the present invention in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and not to limit the scope of the present invention. After reading the present invention, various equivalent forms of modification by those skilled in the art fall within the scope defined by the appended claims of this application.

[0143] As Figure 1 shown, a fast calculation method for the seabed seismic wave field induced by an axisymmetric volume sound source in water includes the following steps:

[0144] S100. Calculate the equivalent sound source intensity; specifically, it includes the following steps:

[0145] S110. According to the wave superposition principle, construct an equation for the relationship between the equivalent sound source intensity and the sound pressure.

[0146] In this specific embodiment, constructing the equation for the relationship between the equivalent sound source intensity and the sound pressure includes the following steps:

[0147] S111. Establish a sound radiation calculation model for the structure in a shallow sea channel; the sound radiation calculation model for the structure in a shallow sea channel includes the following parts:

[0148] The area where the structure is located in the shallow sea channel, the area where the shallow sea fluid medium is located, the area where the solid seabed is located, the coupling boundary between the elastic structure and the fluid, the coupling boundary between the seabed and the seawater, the outer normal direction of the structure, the outer normal direction of the seabed, the virtual surface where the equivalent sound source is located, the coordinate origin, the coordinate of the position vector of the equivalent sound source, the relative position vector coordinate between the equivalent sound source and the point where the sound field is to be obtained, and the seawater depth.

[0149] As Figure 2As shown, in this specific embodiment, the area where the structure is located in the shallow - sea channel is denoted by the symbol Ω S The area where the shallow - sea fluid medium is located is denoted by the symbol Ω W The area where the solid seabed is located is denoted by the symbol Ω B The coupling boundary between the elastic structure and the fluid is denoted by the symbol Γ SW The coupling boundary between the seabed and the sea water is denoted by the symbol Γ BW The outer normal direction of the structure is denoted by the symbol n S The outer normal direction of the seabed is denoted by the symbol n B The virtual surface where the equivalent sound source is located is denoted by the symbol Γ′, the coordinate origin is denoted by the symbol O, and the position - vector coordinate where the equivalent sound source is located is denoted by r 0 =(x 0 ,y 0 ,z 0 ) The relative - position - vector coordinate between the equivalent sound source and the point where the sound field is to be solved is denoted by r=(x,y,z), and the sea - water depth is denoted by the symbol d

[0150] It should be noted that the coordinate origin O is taken near the centroid of the structure and is preset artificially

[0151] It should be noted that according to the wave - superposition principle, the sound pressure p(r S ) on the surface r S of the structure in the shallow - sea channel can be regarded as the linear superposition of the radiated sound fields of all equivalent sources located at the position r 0 on the virtual surface Γ′ where the equivalent sound source is located. Thus, it can be expressed by Equation (1):

[0152] p(r S ) = ikρ f c∫ Γ′ q(r 0 )G k (r S ,r 0 )dS(r 0 ) (1)

[0153] Where: q(r 0 ) is the source strength of the virtual sound source distributed at r 0 ; ρ f is the density of the fluid medium; i is the unit imaginary number; c is the sound - wave speed; s(r 0 ) is the equivalent - sound - source strength; G k (r S ,r 0 ) is the Green's function of the waveguide space and is expressed by Equation (2):

[0154]

[0155] where: δ(r - r 0 ) is the Delta function, and is the Laplace operator.

[0156] S112. Construct an equivalent sound source intensity - sound pressure relationship equation using the sound radiation calculation model of the structure under the shallow - sea channel; the equivalent sound source intensity - sound pressure relationship equation is expressed by Equation (3):

[0157]

[0158] It should be noted that Equation (3) is transformed from Equation (2), and its principle lies in that it can be theoretically proved that Equation (2) is equivalent to the Helmholtz boundary integral equation. Thus, according to Equation (2), the sound source intensity can be regarded as the superposition form of a finite number of point sound sources.

[0159] where: p(r S ) is the sound pressure on the surface of the structure under the shallow - sea channel; r S is the surface of the structure under the shallow - sea channel; r 0 = is the vector coordinate of the position where the equivalent sound source is located, denoted as r 0 = (x 0 , y 0 , z 0 ), which is above the virtual surface where the equivalent sound source is located; k is the wave number of the sound wave, expressed by Equation (4):

[0160] k = ω / c f (4)

[0161] where: ω is the angular frequency; c f is the sound wave propagation speed.

[0162] M is the number of equivalent sources; m is the number - marking symbol of the equivalent source, and m ∈ [0, M].

[0163] S120. Use the mirror image method to construct the series form of the Green's function; then calculate the reflection coefficient under the elastic seabed.

[0164] In this specific embodiment, the series form of the Green's function is expressed by Equation (5):

[0165]

[0166] where: l is the marking symbol, and l ∈ [0, ∞); R l1 , R l2 , R l3 and R l4 are calculated as follows; e is the natural logarithm.

[0167] In the formula:

[0168]

[0169] Wherein: x, y, and z are the relative position vector coordinates between the equivalent sound source and the point to be determined in the sound field, and are uniformly expressed as r = (x, y, z); h is the depth of the seawater; h 1 is the distance from the equivalent sound source to the sea surface; h 2 is the distance from the equivalent sound source to the seabed; γ 1 = -1; γ is the reflection coefficient under the elastic seabed, and is expressed by Equation (6):

[0170]

[0171] Wherein: Z is an intermediate variable; ρ s2 is the seabed density; θ is both the incident wave angle and the reflected wave angle; β is the refracted longitudinal wave angle; α is the refracted transverse wave angle; θ, β, and α satisfy the relationship of Equation (7):

[0172]

[0173] Wherein: c p2 is the transverse wave velocity of the seabed medium; c s2 is the longitudinal wave velocity of the seabed medium.

[0174] It should be noted that, as Figure 3 shown, when the solid seabed parameters are determined, the calculation of the seabed reflection coefficient changes according to different sound wave incident angles; therefore, here, the theory of Brekhovskikh is used to determine the reflection coefficient under the elastic seabed, so as to obtain Equation (6).

[0175] It should be further noted that the solid seabed is the elastic seabed.

[0176] It should be noted that the relationship characterized by Equation (7) is obtained from Snell's theorem.

[0177] It should be further noted that when the structural size is more than 10 times smaller than the distance between the structure and the seabed, the change range of the incident angle θ is very small when solving the intensity of the equivalent sound source, and the corresponding reflection coefficient γ hardly changes. Therefore, in this specific embodiment, the reflection coefficient corresponding to the incident angle θ = 0° is used to calculate the corresponding Green's function G k (r S , r 0 ) and its partial derivatives, so as to achieve the purpose of simplifying the calculation under the condition of ensuring the accuracy.

[0178] It should be noted that the reason why γ 1 = -1 is that here the sea surface is regarded as a pressure release boundary.

[0179] It should be noted that in the waveguide space, sound waves will be reflected multiple times between the seabed and the sea surface. By applying the mirror theory, the Green's function in the waveguide space can be regarded as a superposition form contributed by an infinite number of virtual sources, and thus Equation (5) can be obtained.

[0180] S130. Use Euler's equation to construct the velocity-equivalent sound source relationship equation; then use the velocity-equivalent sound source relationship equation to calculate the equivalent sound source intensity.

[0181] In this specific embodiment, the velocity-equivalent sound source relationship equation is expressed as Equation (8):

[0182]

[0183] where: u(r S ) represents the normal displacement of the surface of the structure; is the Laplace operator; n S is the outer normal direction of the structure, is the normal gradient operator.

[0184] It should be noted that so far, substituting u(r S ) into Equation (8), the equivalent sound source intensity can be solved.

[0185] S200. Construct an axisymmetric finite element coupling calculation model for the equivalent source intensity; specifically, it includes the following steps:

[0186] S210. Construct an axisymmetric finite element coordinate system in cylindrical coordinates.

[0187] In the axisymmetric finite element coordinate system in cylindrical coordinates, the cylindrical coordinate system is expressed as (r, θ, z), which is used to describe the stress and strain of the element; the overall seawater and seabed coupling solution domain takes r = 0 as the axis of symmetry; the equivalent sound source is located on r = 0.

[0188] It should be noted that as Figure 4 shows, in this embodiment, r = 0 is artificially preset, and since r = 0, it can be known that the stress and strain are both independent of the circumferential angle θ.

[0189] S220. Establish the axisymmetric finite element corresponding to the structure.

[0190] In the axisymmetric finite element corresponding to the structure in this specific embodiment, the partial differential equation of the structure in the cylindrical coordinate system with the equivalent sound source located on the axis of symmetry is expressed as Equation (9):

[0191]

[0192] where: σ rr is the normal stress of the axisymmetric element in the r direction; σ zz is the normal stress of the axisymmetric element in the z direction σθθ is the normal stress of the axisymmetric element in the θ direction; τ rz is the shear stress; u is the displacement along the r direction; w is the displacement along the z direction; the axisymmetric element includes an acoustic axisymmetric element and a structural axisymmetric element.

[0193] The strain matrix in the cylindrical coordinate system is expressed by Equation (10):

[0194]

[0195] where: ε rr is the normal strain of the axisymmetric element in the r direction; ε θθ is the normal strain of the axisymmetric element in the z direction; ε zz is the normal strain of the axisymmetric element in the θ direction; γ rz is the shear strain.

[0196] According to Hooke's law, the corresponding linear elastic stress matrix is expressed by Equation (11):

[0197] σ = [σ rr σ θθ σ zz τ rz T = Dε (11)

[0198] where: D is the matrix coefficient, expressed by Equation (12):

[0199]

[0200] where: E is the elastic modulus; v is the Poisson's ratio.

[0201] The radial displacement u and the axial displacement w of the element are interpolated using the same shape function. Thus, the displacement on each axisymmetric element is expressed by Equation (13):

[0202]

[0203] where: q is the nodal displacement vector of the axisymmetric element, and q = [u 1 w 1 … u m w m T ; N s is the matrix composed of the shape function N(ξ, η), expressed by Equation (14):

[0204]

[0205] where: N(ξ, η) is the shape function; ξ is the interpolation function in the ξ direction of the local coordinate system; η is the interpolation function in the η direction of the local coordinate system.​​

[0206] As shown Figure 5 in this specific embodiment, a second-order shape function is used to interpolate and represent the axisymmetric element, and each axisymmetric element corresponds to 9 calculation nodes; the matrix N s formed by the shape function N(ξ,η) has m = 9, and each shape function corresponds to a calculation node; the shape function N(ξ,η) corresponding to each computer point is expressed as follows:

[0207]

[0208]

[0209]

[0210] The strain matrix of the axisymmetric element is expressed by Equation (15):

[0211] ε = Bq (15)

[0212] It should be noted that Equation (15) is obtained by substituting Equation (13) into (10).

[0213] where: B is the stress-displacement matrix, and B = [B 1 , B 2 , …, B m ; B m is expressed by Equation (16):

[0214]

[0215] The stress matrix corresponding to the axisymmetric element is expressed by Equation (17):

[0216] σ = DBq (17)

[0217] It should be noted that Equation (17) is obtained by substituting Equation (15) into (11).

[0218] According to Hamilton's principle, the stiffness matrix corresponding to the axisymmetric element is expressed by Equation (18):

[0219]

[0220] According to Hamilton's principle, the mass matrix corresponding to the axisymmetric element is expressed by Equation (19):

[0221]

[0222] where: J is the Jacobian matrix for converting from the global coordinates (r, z) to the local coordinates (ξ, η), and is expressed by Equation (20):

[0223]

[0224] det(J) is the determinant corresponding to matrix J.

[0225] The elastic seabed response matrix is obtained by calculating the matrices and respectively using Gauss integration, and then combining each axisymmetric element according to each discrete node.

[0226] S230. Establish the axisymmetric finite element corresponding to seawater.

[0227] In the axisymmetric finite element corresponding to seawater in this specific embodiment:

[0228] In the two-dimensional axisymmetric coordinate system, when the equivalent sound source is located in seawater, the control equation of the acoustic axisymmetric element can be expressed by Equation (21) according to the Helmholtz equation:

[0229]

[0230] where: z s is the distance from the equivalent sound source to the seabed; p is the sound pressure on each acoustic axisymmetric element, which is obtained by interpolating the sound pressure using the same shape function as the structure, and is expressed by Equation (22):

[0231] p = N f p (22)

[0232] where: N f is the matrix composed of the shape function N(ξ, η), and N f = [N 1 (ξ, η) N 2 (ξ, η) … N m (ξ, η)]; p is the sound pressure vector of the acoustic axisymmetric element node, and p = [p 1 p 2 … p m T .

[0233] According to Hamilton's principle, the stiffness matrix corresponding to the acoustic axisymmetric element is expressed by Equation (23):

[0234]

[0235] According to Hamilton's principle, the mass matrix corresponding to the acoustic axisymmetric element is expressed by Equation (24):

[0236]

[0237] where: ​

[0238] The seawater corresponding matrix is ​​obtained by converting the matrix and The Guass integral is used for calculation, and then the axisymmetric units of each structure are combined according to each discrete node.

[0239] In the axisymmetric finite element corresponding to the acoustic-solid coupling boundary surface in S240:

[0240] From the compatibility boundary conditions between the structure and the fluid, it can be obtained that at the acoustic-solid interface, the stress continuity condition is expressed as follows:

[0241] p| Γ =(σ·n)| Γ (25)

[0242] At the acoustic-solid interface, the displacement continuity condition is expressed as follows:

[0243]

[0244] Where: Γ is the sound-solid interface; p| Γ is the sound pressure at the interface of the acoustic axisymmetric unit; (σ·n)| Γ is the normal stress at the interface of the axisymmetric unit of the structure; is the projection of the acoustic pressure gradient of the acoustic axisymmetric unit in the normal direction of the interface; (u T ·n)| Γ is the normal displacement of the axisymmetric unit at the interface.

[0245] According to equation (25), equation (26) and the prior art, the coupling matrix is ​​expressed as equation (27):

[0246]

[0247] Where: n is the external normal direction of the acoustic-solid coupling interface structure.

[0248] It should be noted that the construction of the acoustic-solid coupling unit is now complete.

[0249] In the matching infinite cell of the truncated infinite outer domain in S250:

[0250] The mapping relationship of the coordinate position of the matching infinite element from the global coordinate system to the local coordinate system is expressed by formula (28):

[0251] r=M 1 (ξ)r 1 +M 2 (ξ)r 2 (28)

[0252] It should be noted that the matching infinite element, as a type of infinite element, can essentially be regarded as a special form of the finite element. Its function is to truncate the spatial domain of the infinite calculation domain to well solve the problem of solving complex infinite domains. The mapping relationship between the infinite element in the global coordinate system and the local coordinate system is shown in Figure 6 as shown.

[0253] It should be further noted that Figure 6 in which r 0 is the position of the pole.

[0254] Where: M 1 (ξ) and M 2 (ξ) are both shape functions and are expressed by Equation (29):

[0255]

[0256] The shape function M 1 (ξ) and the shape function M 2 (ξ) for matching the infinite element in the global coordinate system are expressed by Equation (30):

[0257]

[0258] Where: the position of the pole is the position of the axis of symmetry; r 1 is the connection point between the matching infinite element and the finite element, and r 3 is the place where the matching infinite element extends to infinity.

[0259] It should be noted that as Figure 7 shown, based on this, 2D matching infinite elements can be constructed. For the 9-node axisymmetric element, 9-node infinite elements as Figure 7 shown can be used for matching. In this way, we can solve and obtain the stiffness matrix and mass matrix corresponding to the corresponding elements.

[0260] S240. Establish the axisymmetric finite element corresponding to the acoustic-solid coupling boundary surface.

[0261] S250. Establish the matching infinite elements for truncating the infinite outer domain.

[0262] S260. Establish the acoustic-vibration coupling dynamic equation under axisymmetric conditions; then establish the axisymmetric finite element coupling calculation model for the equivalent source strength based on the acoustic-vibration coupling dynamic equation under axisymmetric conditions.

[0263] In this specific embodiment, S260 specifically includes the following steps:

[0264] S261. Assemble all the axisymmetric structural elements into an axisymmetric structural element stiffness matrix and an axisymmetric structural element mass matrix; assemble all the axisymmetric acoustic elements into an axisymmetric acoustic element stiffness matrix and an axisymmetric acoustic element mass matrix; assemble all the fluid-structure coupling elements into a fluid-structure coupling element stiffness matrix and a fluid-structure coupling element mass matrix; assemble all the matched infinite elements into a matched infinite element stiffness matrix and a matched infinite element mass matrix.

[0265] S262. According to the corresponding relationships between the nodes, assemble the axisymmetric structural element stiffness matrix, the axisymmetric acoustic element stiffness matrix, the fluid-structure coupling element stiffness matrix, and the matched infinite element stiffness matrix into a coupled stiffness matrix.

[0266] S263. According to the corresponding relationships between the nodes, assemble the axisymmetric structural element mass matrix, the axisymmetric acoustic element mass matrix, the fluid-structure coupling element mass matrix, and the matched infinite element mass matrix into a coupled mass matrix.

[0267] S264. According to the pressure-release boundary condition on the sea surface and the displacement boundary condition related to the axis of symmetry, respectively trim the coupled stiffness matrix and the coupled mass matrix to obtain the acoustic-vibration coupling dynamic equation under axisymmetric conditions.

[0268] It should be noted that, as Figure 8 shown, after obtaining the acoustic-vibration coupling dynamic equation under axisymmetric conditions, the submarine seismic waves caused by the equivalent point source under the elastic seabed condition can be modeled and calculated. Both the sea water and the seabed are truncated in the infinite space by using infinite elements, so as to establish an axisymmetric finite element coupling calculation model for the equivalent source strength; Figure 8 where IE in

[0269] It should be noted that the pressure-release boundary condition on the sea surface and the displacement boundary condition related to the axis of symmetry are expressed by Equation (31):

[0270]

[0271] S265. According to the acoustic-vibration coupling dynamic equation under axisymmetric conditions, respectively truncate the sea water and the seabed in the infinite space by using infinite elements, so as to establish an axisymmetric finite element coupling calculation model for the equivalent source strength.

[0272] S300. Calculate the submarine seismic waves induced by the axisymmetric body target according to the equivalent sound source strength and the axisymmetric finite element coupling calculation model for the equivalent source strength.

[0273] In this specific embodiment, S300 specifically includes the following steps:

[0274] S310. Calculate the intensity corresponding to the equivalent sound source by using the mirror image method and the equivalent source method simultaneously.

[0275] S320. Substitute the position and intensity of the equivalent sound source into the axisymmetric finite element coupling calculation model for the equivalent source intensity, calculate the sound pressure corresponding to the seabed surface in the near field and the sound pressure corresponding to the seabed surface in the far field respectively, and then superimpose the sound pressure corresponding to the seabed surface in the near field and the sound pressure corresponding to the seabed surface in the far field to finally obtain the submarine seismic wave induced by the axisymmetric body target.

[0276] To verify the effectiveness of the method of the present invention, two experiments are also carried out in this specific embodiment to verify the technical effects of the present invention, as follows:

[0277] As follows Figure 9 In the calculation model and coordinate system shown below, first use a pulsating sphere to verify the correctness of the method of the present invention. Assume that the vibration velocity on the surface of the pulsating sphere is v = 1 m / s, the position where the sound source center is located is at a seawater depth h = 50 m, and a corresponding coordinate system is established with the position where the sphere center is located.

[0278] The seawater depth is set to d = 100 m, the sound speed profile corresponding to the shallow sea is a constant value, the seabed is set to conglomerate and basalt respectively, and the specific parameters of the seabed and seawater are shown in Table 1. Calculate the seismic waves at the interface between the seawater and the seabed induced by the pulsating sphere under the conditions of a hard seabed (where the shear wave velocity of the seabed medium is greater than the sound speed in seawater) and a soft seabed (where the shear wave velocity of the seabed medium is less than the sound speed in seawater) respectively. The reference sound pressure level is taken as p ref = 10 -6 Pa, and the analysis frequencies are selected as 30 Hz and 50 Hz.

[0279] Table 1. List of medium types and parameters

[0280] Medium type Longitudinal wave velocity (m / s) Shear wave velocity (m / s) <![CDATA[Density (g / cm 3 )]]> Seawater 1500 0 1.0 Basalt 3500 1800 2.3 Conglomerate 2000 800 2.1

[0281] As Figure 10 shown, distribute the equivalent point sources on the axis coinciding with the z-axis, and distribute the field points on the generatrix of the spherical shell surface. Assume that the radius of the sphere is R = 2.5 m; the distribution diagrams of the corresponding equivalent sources and field points are as Figure 10 shown.

[0282] The number of equivalent point sources and the corresponding nodes on the structure surface are both taken as n = N j = 20, the r of the source point coordinates distributed on the z-axis is 0, and the z coordinate of the n j (n j = 1, 2,..., 20)th equivalent point source is expressed by Equation (32):

[0283]

[0284] Correspondingly, the r and z coordinates of the n j (n j = 1, 2, …, 20) field points are expressed by equations (33) and (34) respectively.

[0285]

[0286] In the axisymmetric finite element model, the equivalent point source Q is located on the axis of symmetry of the model, that is, at r = 0. Its horizontal distance is 2000 m, the seawater depth is 100 m, and the seabed depth is 100 m. Its calculation grid is discretized using quadratic quadrilateral elements, and the specific grid size is 5 m × 5 m.

[0287] First, the mirror method and the equivalent source method are used to calculate the corresponding intensity of the sound source. Then, the information related to the equivalent source position and intensity is substituted into the axisymmetric finite element to calculate the sound pressure corresponding to the near field and the far field on the seabed surface and perform superposition, so as to calculate the seismic waves radiated by the spherical elastic seabed on the seabed surface.

[0288] The obtained results are shown in Figures 11a - 11d below; it should be noted that the solid line "-" in the figure represents the result of simulation calculation by Comsol software, and "*" represents the result of calculation by the combined method.

[0289] Thus, it can be seen from Figures 11a - 11d that under the conditions of hard seabed and soft seabed, the sound pressure propagation loss on the seawater - seabed interface calculated by the combined method of equivalent source and axisymmetric finite element is compared with the result of Comsol finite element software. The overall change trends of the two are consistent, and the results are in good agreement. Therefore, the method of the present invention has good calculation accuracy.

[0290] In the above detailed description, various features are combined in a single embodiment to simplify the present disclosure. This method of disclosure should not be interpreted as reflecting the intention that the embodiments of the claimed subject matter require more features than those clearly stated in each claim. On the contrary, as reflected in the appended claims, the present invention is in a state with fewer features than all the features of the disclosed single embodiment. Therefore, the appended claims are hereby clearly incorporated into the detailed description, where each claim stands alone as a separate preferred embodiment of the present invention.

[0291] In order to enable any person skilled in the art to implement or use the present invention, the above - described disclosed embodiments are described. For those skilled in the art; various modification methods of these embodiments are obvious, and the general principles defined herein can also be applied to other embodiments without departing from the spirit and protection scope of the present disclosure. Therefore, the present disclosure is not limited to the embodiments given herein, but is consistent with the broadest scope of the principles and novel features disclosed in this application.

[0292] The foregoing description includes examples of one or more embodiments. Of course, it is not possible to describe all possible combinations of components or methods for the purpose of describing the above embodiments, but those of ordinary skill in the art should recognize that the various embodiments can be further combined and arranged. Accordingly, the embodiments described herein are intended to cover all such changes, modifications, and variations that fall within the scope of the appended claims. In addition, with respect to the term "comprising" as used in the specification or claims, the term is inclusive in a manner similar to the term "including," as is explained when "including" is used as a transitional word in a claim. Further, any use of the term "or" in a claim of the specification is to mean "non-exclusive or."

[0293] The specific embodiments described above further elaborate on the object, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only the specific embodiments of the present invention and is not used to limit the protection scope of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A fast calculation method for the seabed seismic wave field induced by an axisymmetric volume sound source in water, characterized in that: It includes the following steps: S100. Calculate the equivalent sound source intensity; specifically, it includes the following steps: S110. Construct an equation for the relationship between equivalent sound source intensity and sound pressure; S120. Construct a series form of the Green's function; then calculate the reflection coefficient under the elastic seabed; S130. Construct an equation for the relationship between vibration velocity and equivalent sound source; then use the equation for the relationship between vibration velocity and equivalent sound source to calculate the equivalent sound source intensity; S200. Construct an axisymmetric finite element coupling calculation model for the equivalent source intensity; specifically, it includes the following steps: S210. Construct an axisymmetric finite element coordinate system in cylindrical coordinates; S220. Establish an axisymmetric finite element corresponding to the structure; S230. Establish an axisymmetric finite element corresponding to the seawater; S240. Establish an axisymmetric finite element corresponding to the acoustic-solid coupling boundary surface; S250. Establish a matching infinite element for truncating the infinite outer domain; S260. Establish an acoustic-vibration coupling dynamic equation under axisymmetric conditions; then establish the axisymmetric finite element coupling calculation model for the equivalent source intensity according to the acoustic-vibration coupling dynamic equation under axisymmetric conditions; S300. Calculate the seabed seismic waves induced by the axisymmetric body target according to the equivalent sound source intensity and the axisymmetric finite element coupling calculation model for the equivalent source intensity.

2. The fast calculation method for the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 1, characterized in that: Constructing the equation for the relationship between equivalent sound source intensity and sound pressure in S110 includes the following steps: S111. Establish a sound radiation calculation model for the structure in the shallow sea channel; the sound radiation calculation model for the structure in the shallow sea channel includes the following parts: The area where the structure is located in the shallow sea channel, the area where the shallow sea fluid medium is located, the area where the solid seabed is located, the coupling boundary between the elastic structure and the fluid, the coupling boundary between the seabed and the seawater, the outer normal direction of the structure, the outer normal direction of the seabed, the virtual surface where the equivalent sound source is located, the coordinate origin, the coordinate of the position vector of the equivalent sound source, the relative position vector coordinate between the equivalent sound source and the point where the sound field is to be obtained, and the seawater depth; S112. Use the sound radiation calculation model for the structure in the shallow sea channel to construct the equation for the relationship between equivalent sound source intensity and sound pressure; The equation for the relationship between equivalent sound source intensity and sound pressure is expressed as follows: Where: p(r S ) is the sound pressure on the surface of the structure in the shallow sea channel; r S is the surface of the structure in the shallow sea channel; r 0 = is the position vector coordinate of the equivalent sound source, expressed as r 0 = (x 0 , y 0 , z 0 ), located above the virtual plane where the equivalent sound source is located; k is the wave number of the sound wave, expressed by the following formula: k = ω / c f where: ω is the angular frequency; c f is the sound wave propagation speed; ρ f is the density of the fluid medium; i is the unit imaginary number; c is the acoustic wave speed; s(r 0 ) is the equivalent sound source strength; G k (r S , r 0 ) is the Green's function of the waveguide space, expressed by the following formula: where: δ(r - r 0 ) is the Delta function, is the Laplace operator; M is the number of equivalent sources; m is the number marking symbol of the equivalent source, and m ∈ [0, M].

3. The fast calculation method for the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 2, characterized in that: The series form of the Green's function is expressed as follows: where: l is the order of the virtual source, and l ∈ [0, ∞); R l1 、R l2 、R l3 and R l4 The calculation formulas are as follows; e is the natural logarithm; Where: where: x, y, z are the relative position vector coordinates between the equivalent sound source and the point to be determined in the sound field, and are uniformly expressed as r = (x, y, z); h is the depth of the sea water; h 1 is the distance from the equivalent sound source to the sea surface; h 2 is the distance from the equivalent sound source to the seabed; γ 1 = -1; γ is the reflection coefficient under the elastic seabed, and is expressed by the following formula: Where: Z is an intermediate variable; ρ s2 is the seabed density; θ is both the incident wave angle and the reflected wave angle; β is the refracted longitudinal wave angle; α is the refracted shear wave angle; θ, β, and α satisfy the following relationship: where: c p2 is the shear wave velocity of the seabed medium; c s2 is the compressional wave velocity of the seabed medium.

4. The fast calculation method for the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 3, characterized in that: The equation for the relationship between vibration velocity and equivalent sound source in S130 is expressed as follows: where: u(r S ) represents the normal displacement of the surface of the structure; is the Laplace operator; n S is the outer normal direction of the structure, is the normal gradient operator.

5. The fast calculation method for the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 1, characterized in that: In the axisymmetric finite element coordinate system in cylindrical coordinates in S210, the cylindrical coordinate system is expressed as (r, θ, z); the overall seawater and seabed coupled solution domain has r = 0 as the axis of symmetry; the equivalent sound source is located on r = 0.

6. The rapid calculation method of the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 1, characterized in that: In the axisymmetric finite element corresponding to the structure in S220, when the equivalent sound source is located on the axis of symmetry, the partial differential equation of the structure in the cylindrical coordinate system is expressed as follows: Where: σ rr is the normal stress in the r direction of the axisymmetric element; σ zz is the normal stress in the z direction of the axisymmetric element σ θθ is the normal stress in the θ direction of the axisymmetric element; τ rz is the shear stress; u is the displacement along the r direction; w is the displacement along the z direction; the axisymmetric element includes an acoustic axisymmetric element and a structural axisymmetric element; The strain matrix in the cylindrical coordinate system is expressed as follows: Where: ε rr is the normal strain of the axisymmetric element in the r direction; ε θθ is the normal strain of the axisymmetric element in the z direction; ε zz is the normal strain of the axisymmetric element in the θ direction; γ rz is the shear strain; The corresponding linear elastic stress matrix is expressed as follows: σ = [σ rr σ θθ σ zz τ rz T = Dε​ where: D is the matrix coefficient, expressed as follows: where: E is the elastic modulus; v is the Poisson's ratio; The displacement on each axisymmetric unit is expressed as follows: where: q is the nodal displacement vector of the axisymmetric element, and q = [u 1 w 1 …u m w m T ; N s is the matrix composed of the shape functions N(ξ, η), expressed by the following formula:​ where: N(ξ, η) is the shape function; ξ is the interpolation function in the ξ direction in the local coordinate system; η is the interpolation function in the η direction in the local coordinate system; The strain matrix of the axisymmetric unit is expressed as follows: ε = Bq Where: B is the stress-displacement matrix, and B = [B 1 , B 2 , …, B m ; B m is expressed by the following formula: The stress matrix corresponding to the axisymmetric unit is expressed as follows: σ = DBq The stiffness matrix corresponding to the axisymmetric unit is expressed as follows: The mass matrix corresponding to the axisymmetric unit is expressed as follows: where: J is the Jacobian matrix for converting from the global coordinates (r, z) to the local coordinates (ξ, η), expressed as follows: det(J) is the determinant corresponding to the matrix J; The elastic seabed response matrix is obtained by calculating the matrix and the matrix using Gauss integration respectively, and then combining the axisymmetric elements according to each discrete node.

7. The rapid calculation method of the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 6, characterized in that: Each axisymmetric unit corresponds to 9 computational nodes; in the matrix N composed of the shape functions N(ξ, η), m = 9, and each shape function corresponds to one of the computational nodes; the shape functions N(ξ, η) are expressed as follows: s where m = 9, and each shape function corresponds to one of the computational nodes; the shape functions N(ξ, η) are expressed as follows:

8. The rapid calculation method of the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 7, characterized in that: In the axisymmetric finite element corresponding to the seawater in S230: When the equivalent sound source is located in the seawater, the control equation of the acoustic axisymmetric unit is expressed as follows: where: p is the sound pressure on each acoustic axisymmetric unit, expressed as follows: p = N f p Where: N f is a matrix composed of the shape function N(ξ,η), and N f = [N 1 (ξ,η) N 2 (ξ,η)…N m (ξ,η)]; p is the acoustic pressure vector at the nodes of the axisymmetric acoustic element, and p = [p 1 p 2 …p m T ;​ The stiffness matrix corresponding to the acoustic axisymmetric unit is expressed as follows: The mass matrix corresponding to the acoustic axisymmetric unit is expressed as follows: Wherein: The seawater response matrix is obtained by combining the matrix with which is calculated using Guass integration, and then combining the axisymmetric elements of each structure according to each discrete node; In the axisymmetric finite element corresponding to the acoustic-solid coupling boundary surface in S240: On the acoustic-solid interface, the stress continuity condition is expressed as follows: p| Γ =(σ·n)| Γ On the acoustic-solid interface, the displacement continuity condition is expressed as follows: Where: Γ is the acoustic-solid interface; p| Γ is the sound pressure at the interface of the acoustic axisymmetric element; (σ·n)| Γ is the normal stress at the interface of the structural axisymmetric element; is the projection of the sound pressure gradient of the acoustic axisymmetric element in the normal direction of the interface; (u T ·n)| Γ is the normal displacement of the structural axisymmetric element at the interface; The coupling matrix is expressed as follows: where: n is the outer normal direction of the acoustic-solid coupling interface structure; In the matching infinite element of the truncated infinite outer domain in S250: The mapping relationship of the coordinate position of the matching infinite element from the global coordinate system to the local coordinate system is expressed as follows: r = M 1 (ξ)r 1 +M 2 (ξ)r 2 where: M 1 (ξ) and M 2 (ξ) are both shape functions and are expressed by the following formula: Shape function M 1 (ξ) and shape function M 2 (ξ) for matching the infinite element in the global coordinate system is expressed by the following formula: Wherein: the position of the pole is the position of the axis of symmetry; r 1 is for matching the connection between the infinite element and the finite element, r 3 is for matching the infinite element extending to infinity.

9. The rapid calculation method of the seabed seismic wave field induced by an axisymmetric volume sound source in water according to claim 8, characterized in that: S260 specifically includes the following steps: S261. Assemble all the structural axisymmetric units into a structural axisymmetric unit stiffness matrix and a structural axisymmetric unit mass matrix; assemble all the acoustic axisymmetric units into an acoustic axisymmetric unit stiffness matrix and an acoustic axisymmetric unit mass matrix; Assemble all the acoustic-solid coupling units into an acoustic-solid coupling unit stiffness matrix and an acoustic-solid coupling unit mass matrix; Assemble all the matching infinite units into a matching infinite unit stiffness matrix and a matching infinite unit mass matrix; S262. Assemble the structural axisymmetric element stiffness matrix, the acoustic axisymmetric element stiffness matrix, the fluid-structure coupling element stiffness matrix, and the matched infinite element stiffness matrix into a coupling stiffness matrix according to the corresponding relationships between nodes; S263. Assemble the structural axisymmetric element mass matrix, the acoustic axisymmetric element mass matrix, the fluid-structure coupling element mass matrix, and the matched infinite element mass matrix into a coupling mass matrix according to the corresponding relationships between nodes; S264. Cut the coupling stiffness matrix and the coupling mass matrix respectively according to the pressure release boundary condition on the sea surface and the displacement boundary condition related to the axis of symmetry, so as to obtain the acoustic-vibration coupling dynamic equation under the axisymmetric condition; S265. According to the acoustic-vibration coupling dynamic equation under the axisymmetric condition, truncate the seawater and the seabed using infinite elements for the infinite space, so as to establish the axisymmetric finite element coupling calculation model for the equivalent source strength.

10. The fast calculation method of the submarine seismic wave field induced by an axisymmetric volume sound source in water according to claim 9, characterized in that: S300 specifically includes the following steps: S310. Calculate the intensity corresponding to the equivalent sound source by using the mirror image method and the equivalent source method simultaneously; S320. Substitute the position and intensity of the equivalent sound source into the axisymmetric finite element coupling calculation model for the equivalent source strength, calculate the sound pressure on the seabed surface corresponding to the near field and the sound pressure on the seabed surface corresponding to the far field respectively, and then superimpose the sound pressure on the seabed surface corresponding to the near field and the sound pressure on the seabed surface corresponding to the far field, and finally obtain the submarine seismic wave induced by the axisymmetric body target.

Images