CN112254798B - Method, system and medium for forecasting ocean vector sound field - Google Patents

Abstract

The invention discloses a method, a system and a medium for forecasting an ocean vector sound field, wherein the method comprises the following steps: s1, acquiring field measurement data and sound source parameter information of a to-be-forecasted ocean area, establishing a cylindrical coordinate system underwater sound Helmholtz equation, and obtaining a depth equation after transformation; s2, solving a depth equation to obtain a sound pressure kernel function and a vertical vibration velocity kernel function of each medium layer interface; s3, calculating sound pressure according to a sound pressure kernel function, calculating a vertical vibration velocity by using a vertical direction derivative based on a Hankel inverse transformation integral formula according to a vertical vibration velocity kernel function, and calculating a horizontal vibration velocity based on a Hankel inverse transformation integral formula horizontal direction derivative; and S4, calculating the underwater sound propagation loss and the sound intensity vector of the ocean area to be forecasted according to the solved sound pressure and vibration velocity vector. The method can realize the prediction of the vibration velocity vector based on the marine environment measurement data, and can improve the precision of the vibration velocity vector prediction.

Description

A method, system and medium for predicting ocean vector sound field

technical field

The invention relates to the technical field of underwater sound field detection, in particular to a method, a system and a medium for predicting an ocean vector sound field.

Background technique

Acoustic waves are the main means of underwater communication, marine environment and target detection, and have important application value in military and economic fields. The receiving sensors of underwater sound waves are generally called hydrophones, and traditional hydrophones are scalar hydrophones, which can only measure scalar parameters (such as sound pressure) in the sound field. At present, the more advanced hydrophone is the vector hydrophone, which can measure both scalar parameters and vector parameters in the sound field (such as particle vibration velocity, referred to as vibration velocity), which is of great significance for improving the performance of underwater detection systems.

Due to the seabed topography, marine environment (mainly sound speed, density, with time-varying characteristics), sound source frequency and location and other factors have an important impact on underwater sound propagation, in the placement of hydrophones, array shape optimization, signal reception In the process of processing and analysis, it is necessary to forecast the ocean sound field, that is, to carry out numerical simulation of the sound field in combination with the dynamic changing ocean environment, and integrate the forecast results into the hydrophone signal processing process. It plays an important role in the application level and enhancing the performance of the underwater detection system.

However, traditional underwater acoustic models (such as wavenumber integration method, normal wave method, etc.) can only predict the scalar sound field, and to obtain the vector sound field, it is necessary to arrange denser grid points in the sound field. Then the finite difference method is used to calculate the sound pressure derivative of each point, and finally the sound pressure derivative is converted into vibration velocity. The above-mentioned traditional method for calculating the sound field vector (vibration velocity) using the finite difference method has the following shortcomings:

(1) The grid spacing of the sound field is limited by the finite difference method. In order to ensure the correctness of the derivative difference calculation, the grid spacing needs to be small enough (the sound pressure value does not change much in the interval between adjacent grid points), and the grid spacing in each direction generally takes several points of the reference wavelength. One (for example, one tenth), in the case that only a few spatial positions around the hydrophone need to be calculated, this limitation will adversely affect the calculation amount of the sound field vector and the flexibility of the calculation method.

(2) The calculation of the sound pressure derivative by the finite difference method will introduce additional numerical errors. In practical applications, the finite difference method of finite order (such as second-order accuracy) is generally used. When the grid spacing is a finite value (not tending to zero), the process of differential calculation of the derivative will introduce the dissipation and dispersion errors of the finite difference method. .

SUMMARY OF THE INVENTION

The technical problem to be solved by the present invention is: aiming at the technical problems existing in the prior art, the present invention provides a method, system and medium for predicting the ocean vector sound field, which can realize the prediction of the vibration velocity vector based on the marine environment measurement data, and based on the The vibration velocity calculation method that directly integrates the horizontal wavenumber can improve the accuracy of the vibration velocity vector prediction.

In order to solve the above-mentioned technical problems, the technical scheme proposed by the present invention is:

A method for predicting an ocean vector sound field, the steps comprising:

S1. Obtain on-site measurement data and sound source parameter information of the ocean area to be predicted, establish a cylindrical coordinate system hydroacoustic Helmholtz equation in a horizontally layered ocean environment, and obtain a depth equation with a sound pressure kernel function as a variable after integral transformation;

S2. Use the transfer function matrix method to solve the depth equation, and obtain the sound pressure kernel function and the vertical vibration velocity kernel function of the interface of each medium layer;

S3. Calculate the sound pressure using the Hankel inverse transform integral formula according to the sound pressure kernel function, calculate the vertical vibration velocity using the vertical direction derivative based on the Hankel inverse transform integral formula according to the vertical vibration velocity kernel function, and calculate the vertical vibration velocity based on the Hankel inverse transform integral formula Take the directional derivative and transform the derivative of the zero-order Bessel function into a first-order Bessel function to obtain the horizontal vibration velocity integral formula to calculate the horizontal vibration velocity and obtain the vibration velocity vector;

S4. Calculate the underwater sound propagation loss and sound intensity vector of the ocean to be predicted according to the sound pressure and vibration velocity vector obtained in step S3.

Further, the steps of step S1 include:

S11. Establish the hydroacoustic Helmholtz equation of the cylindrical coordinate system according to the following formula:

Among them, P(r,z) is the relative sound pressure in the frequency domain, ρ is the density of the sound propagation medium, k is the wave number and k=2πf/c, f is the frequency of the sound source, c is the sound velocity of the medium, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, z _s is the depth of the sound source, and δ is the Dirac function;

S12. Divide the underwater acoustic propagation medium into N layers in the depth direction, and approximate each layer as a homogeneous medium; in each divided layer, perform Hankel transformation on the hydroacoustic Helmholtz equation of the cylindrical coordinate system to convert The sound pressure P(r,z) in the (r,z) space is converted to the (k _r ,z) space, that is:

Among them, φ(k _r , z) is the sound pressure kernel function, k _r is the horizontal wave number, and J ₀ is the Bessel function;

可得到深度方程为：Integrate both sides of the cylindrical hydroacoustic Helmholtz equation at the same time

The depth equation can be obtained as:

Further, when solving the depth equation in the step S2, it includes the step of forming a joint vector between the sound pressure kernel function and the vertical vibration velocity kernel function, and the specific steps include:

S211. In any medium layer without sound source, the general solution form of the sound pressure kernel function is:

A⁺(k_r)表示向下传播的项，A^-(k_r)表示向上传播的项，w(k_r,z)为z方向的垂直振速核函数，且满足：where φ(k _r ,z) is the sound pressure kernel function, k _z is the vertical wave number and

A ⁺ (k _r ) represents the downward propagating term, A ^- (k _r ) represents the upward propagating term, and w(k _r ,z) is the vertical vibration velocity kernel function in the z direction, and satisfies:

Among them, Γ=k _z /(ρω), ρ is the density of the sound propagation medium, ω=2πf is the angular frequency of the sound source vibration, and f is the sound source vibration frequency;

Then the joint vector formed by the sound pressure kernel function and the vertical vibration velocity kernel function is:

S212. According to the joint vector formed in step S211, the joint vector obtained at the upper interface z _m-1 of the mth layer is:

And the joint vector at the lower interface z _m of the mth layer is:

后，得到所述联合矢量由下至上传递公式为：S213. Combine and eliminate the joint vector expressions of the upper and lower interfaces of the mth layer obtained in step S212

After that, the bottom-up transfer formula of the joint vector is obtained as:

v _m (k _r ,z _m-1 )=M _m (k _r )v _m (k _r ,z _m )

Among them, M _m (k _r ) is the transfer matrix of the m-th layer medium from bottom to top. If the thickness of the m-th layer is h _m =z _m -z _m-1 , the expression of M _m (k _r ) is:

And the top-to-bottom transfer formula and transfer matrix are:

Further, the specific steps of obtaining the sound pressure kernel function and the vertical vibration velocity kernel function of the interface of each medium layer in the step S2 include:

使用下标“0”表示上边界，得到上边界处的声压核函数φ₀(k_r,z₀)、垂直振速核函数w₀(k_r,z₀)分别为：S221. According to the boundary surface acoustic energy on the computational domain, only the items that can propagate upward and propagate downward are zero, that is,

Using the subscript "0" to represent the upper boundary, the sound pressure kernel function φ ₀ (k _r , z ₀ ) and the vertical vibration velocity kernel function w ₀ (k _r , z ₀ ) at the upper boundary are obtained as:

and the joint vector at the upper boundary is:

Among them, w ₀ (k _r , z ₀ ) is the vertical vibration velocity kernel function of the upper boundary, and the vector (1, B ₀ ) ^T is the sound vector of the upper boundary, where

使用下标“N”表示下边界，得到下边界处的联合矢量为：According to the boundary of the computational domain, the boundary acoustic energy can only propagate downward and the upward propagation term is zero, that is,

Using the subscript "N" to denote the lower boundary, the joint vector at the lower boundary is:

Among them, w _N (k _r , z _N ) is the vertical vibration velocity kernel function of the lower boundary, φ _N (k _r , z _N ) is the sound pressure kernel function of the lower boundary, (1, B _N ) ^T is the lower boundary sound vector, where

S222. Transfer the sound vector from the lower boundary _zN to the sound source depth _zs , that is, starting from the lower boundary _zN , transfer and calculate the sound vector layer by layer until the sound source depth _zs , where the joint vector immediately below the sound source depth is calculated The formula is:

In the formula, w _N (k _r , z _N ) is the undetermined lower boundary vertical vibration velocity kernel function, and the sound vector is immediately below the sound source depth

S223. Transfer the sound vector from the upper boundary z ₀ to the sound source depth z _s , wherein the joint vector calculation formula obtained by the transfer immediately above the sound source depth is:

In the formula, w ₀ (k _r , z ₀ ) is the undetermined upper boundary vertical vibration velocity kernel function, and the sound source depth is just above the sound vector

S224. Calculate the vertical vibration velocity kernel functions of the lower and upper boundaries according to the sound source interface conditions, which satisfy at the sound source interface:

That is, v _s+1 (k _r ,z _s )-v _s (k _r ,z _s )=Δv(k _r ,z _s ), expand to get:

Then, the lower and upper boundary vertical vibration velocity kernel functions are obtained respectively as:

_S225 . _Calculate the sound pressure _kernel function _{φ ( k r} , z), the value of the vertical vibration velocity kernel function w(k _r , z).

Further, when solving the sound pressure in the step S3, the specific Hankel inverse transform integral formula is:

After discretizing the horizontal wavenumber k _r in the inverse Hankel transform integral formula, the sound pressure discrete calculation formula is:

Among them, Δk _r is the horizontal wave number step size and Δk _r =2π/(r _max n _w ), r _max is the maximum horizontal distance of the sound field, n _w is the minimum number of sampling points of the Bessel function in a 2π oscillation period, k _r,n is Discrete horizontal wave number and k _r,n =nΔk _r -iε _k , i is an imaginary unit, ε _k is a complex offset and ε _k =3Δk _r /(2πlog ₁₀ e), M is the maximum index of the discrete horizontal wave number and M=km _max /Δk _r , where _{km max} is the cutoff wave number. Using the discrete calculation formula of the sound pressure at each position (r, z) of the sound field to perform horizontal wavenumber integration calculation, the sound pressure corresponding to each position (r, z) is obtained.

Further, when calculating the vibration velocity vector in the step S3, the step of solving the vertical vibration velocity includes:

S311. Calculate the derivative of the integral formula of the Hankel inverse transform in the depth z direction, and obtain the first transform formula:

Among them, P(r, z) is the relative sound pressure in the frequency domain, φ(k _r , z) is the sound pressure kernel function, k _r is the horizontal wave number, r is the coordinate in the horizontal direction, and z is the coordinate in the vertical or depth direction. , w(k _r , z) is the vertical vibration velocity kernel function, J ₀ is the zero-order Bessel function, ρ is the sound propagation medium density, ω=2πf is the sound source vibration angular frequency, and f is the sound source vibration frequency;

S312. The vertical vibration velocity component V _z (r, z) in the vibration velocity vector V(r, z)=(V _r (r, z), V _z (r, z)) of the particle is transformed according to the first Integrate the formula to get:

And the horizontal wave number k _r is discretized, the discrete integral calculation formula of vertical vibration velocity is obtained as:

S313. Perform horizontal wavenumber integral calculation at each position (r, z) of the sound field using the discrete integral calculation formula of the vertical vibration velocity to obtain the vertical vibration velocity corresponding to each position (r, z).

Further, when calculating the vibration velocity vector in the step S3, the specific steps of solving the horizontal vibration velocity include:

得到第二变换式为：S321. Calculate the derivative of the Hankel inverse transform integral in the depth r direction, and convert the derivative of the zero-order Bessel function into a first-order Bessel function, namely

The second transformation is obtained as:

为一阶Bessel函数；Among them, P(r, z) is the relative sound pressure in the frequency domain, φ(k _r , z) is the sound pressure kernel function, k _r is the horizontal wave number, r is the coordinate in the horizontal direction, and z is the coordinate in the vertical or depth direction. , J ₀ is a zero-order Bessel function,

is a first-order Bessel function;

S322. The horizontal vibration velocity component V _r (r, z) in the vibration velocity vector obtained by integral calculation according to the second transformation formula is:

Among them, ρ is the density of the sound propagation medium, ω=2πf is the angular frequency of the sound source vibration, and f is the sound source vibration frequency;

And the horizontal wave number k _r is discretized, and the discrete integral calculation formula of the horizontal vibration velocity is obtained as:

S323. Perform horizontal wavenumber integral calculation at each position (r, z) of the sound field using the discrete integral calculation formula of the horizontal vibration velocity to obtain the horizontal vibration velocity corresponding to each position (r, z).

Further, in the step S4, it also includes calculating the propagation loss according to the sound pressure obtained in the step S3, forming a scalar cloud map of the sound field propagation loss, and calculating the sound pressure P(r, z) and the vibration velocity vector V(r according to the obtained sound pressure P(r, z). , z)=(V _r (r, z), V _z (r, z)) according to the following formula to calculate the time-averaged sound intensity vector of each position (r, z) of the sound field;

分别表示复数声压P(r,z)、振速V(r,z)的共轭复数，r为水平方向的坐标，z为竖直或深度方向的坐标。Among them, I(r,z) is the time-averaged sound intensity vector, that is, on the plane perpendicular to the direction of the vibration velocity, the time-averaged sound energy passing through a unit area,

Represents the complex conjugate of the complex sound pressure P(r,z) and the vibration velocity V(r,z) respectively, r is the coordinate in the horizontal direction, and z is the coordinate in the vertical or depth direction.

A system for predicting an ocean vector sound field, comprising a processor and a memory, the memory is used for storing marine environment data, sound source parameters and a computer program, the processor is used for executing the computer program, the processor is used for executing all the The computer program described above is used to perform the above method.

A computer-readable storage medium storing marine environment data, sound source parameters, and a computer program programmed or configured to execute the above-described method for predicting an oceanic vector sound field.

Compared with the prior art, the advantages of the present invention are: the method, system and medium for predicting the ocean vector sound field of the present invention, after acquiring the on-site measurement data and sound source parameter information of the ocean area to be predicted, use the wave number integration method for underwater sound. On the basis of the sound pressure calculated by the model, the integral formula of the vertical vibration velocity and the horizontal vibration velocity can be obtained by derivation of the integral formula of the Hankel inverse transformation. The limitation on the distance between the grid points of the sound field caused by the differential calculation of the sound pressure derivative to calculate the vibration velocity caused by the traditional method can be avoided, and the numerical error of the finite difference method can be avoided, thereby improving the prediction accuracy of the vibration velocity vector.

Description of drawings

FIG. 1 is a schematic diagram of a detailed implementation flow of the method for predicting an ocean vector sound field in this embodiment.

FIG. 2 is a schematic diagram of a system structure for implementing a method for predicting an ocean vector sound field in a specific application embodiment of the present invention.

FIG. 3 is a scalar cloud diagram of sound field propagation loss obtained in a specific application embodiment of the present invention.

FIG. 4 is a streamline diagram of a sound intensity vector obtained in a specific application example of the present invention.

Detailed ways

The present invention will be further described below with reference to the accompanying drawings and specific preferred embodiments, but the protection scope of the present invention is not limited thereby.

This embodiment specifically takes the sound field prediction in the marine environment of absolute hard seabed isokinetic waveguides as an example, and the parameters are: uniform seawater density ρ _w =1.0g/cm ³ , uniform sound speed of water body c _w =1500m/s, seabed level And the sea depth z _N =100m, the depth z direction step size dz=1m, the maximum solution distance in the r direction is r _max =1000m, and the step size Δr=1m. The sound source frequency f=100Hz, the sound source depth z _s =25m, the upper boundary (sea surface z ₀ =0m) takes the pressure release boundary condition (sound pressure is zero), and the lower boundary (the seabed z _N =100m) takes the absolute hard boundary condition (the sound pressure z derivative is zero).

As shown in FIG. 1 , the detailed steps of the method for predicting an ocean vector sound field in this embodiment include: S1. Obtain on-site measurement data and sound source parameter information of the ocean area to be predicted, and establish a cylindrical coordinate system underwater acoustic in a horizontally layered ocean environment Helmholtz equation, and after integral transformation, the depth equation with sound pressure kernel function as variable is obtained;

S2. Use the transfer function matrix method to solve the depth equation, and obtain the sound pressure kernel function and the vertical vibration velocity kernel function of the interface of each medium layer;

S3. Calculate the sound pressure using the Hankel inverse transform integral formula according to the sound pressure kernel function, calculate the vertical vibration velocity using the vertical direction derivative based on the Hankel inverse transform integral formula according to the vertical vibration velocity kernel function, and calculate the vertical vibration velocity based on the Hankel inverse transform integral formula Take the directional derivative and transform the derivative of the zero-order Bessel function into a first-order Bessel function to obtain the horizontal vibration velocity integral formula to calculate the horizontal vibration velocity and obtain the vibration velocity vector;

S4. Calculate the underwater sound propagation loss and sound intensity vector of the ocean to be predicted according to the sound pressure and vibration velocity vector obtained in step S3.

其中J₀为零阶Bessel函数，J′₀为J₀的导数，J₁为一阶Bessel函数，即可以将零阶Bessel函数转换为一阶Bessel函数，则在Hankel反变换积分式的基础上，在水平方向求导后利用Bessel函数的上述特性，又可以得到水平振速积分式，进而可以求解出水平振速。本实施例利用上述特性，在获取到待预报海洋区域的现场测量数据以及声源参数信息后，使用波数积分法水声模型计算声压，同时通过对Hankel反变换积分式求导，分别获得垂直振速与水平振速的波数积分式，其中利用垂直振速核函数积分计算出垂直振速，利用Bessel函数的特性以及Hankel反变换积分式水平方向导数求解出水平振速，得到振速矢量，从而可实现水声传播损失与声强矢量的预报。本实施例上述方法，基于直接进行水平波数积分的振速计算方式，可以直接计算获得声场任意点的振速矢量，突破了传统方法因差分计算声压导数求振速引起的对声场网格点间距的限制，同时可避免有限差分法的数值误差从而提高振速矢量预报精度，进而提升矢量水听器的应用技术水平。The sound pressure can be calculated using the wavenumber integration method for the underwater acoustic model. At the same time, on the basis of the wavenumber integration method, the vertical vibration velocity integration formula can be obtained by using the Hankel inverse transform integral formula to obtain the vertical vibration velocity integral formula. , while considering that the Bessel function has the property:

Among them, J ₀ is the zero-order Bessel function, J′ ₀ is the derivative of J ₀ , and J ₁ is the first-order Bessel function, that is, the zero-order Bessel function can be converted into a first-order Bessel function, then on the basis of the Hankel inverse transform integral formula , after the derivation in the horizontal direction, the above-mentioned characteristics of the Bessel function can be used, and the integral formula of the horizontal vibration velocity can be obtained, and then the horizontal vibration velocity can be solved. In this embodiment, using the above characteristics, after acquiring the on-site measurement data and sound source parameter information of the ocean area to be predicted, the underwater acoustic model of the wave number integration method is used to calculate the sound pressure. The wavenumber integral formula of the vibration velocity and the horizontal vibration velocity, in which the vertical vibration velocity is calculated by using the vertical vibration velocity kernel function integral, and the horizontal vibration velocity is solved by using the characteristics of the Bessel function and the horizontal direction derivative of the Hankel inverse transformation integral formula, and the vibration velocity vector is obtained, Thus, the prediction of underwater acoustic propagation loss and sound intensity vector can be realized. The above method of this embodiment, based on the vibration velocity calculation method of directly performing horizontal wavenumber integration, can directly calculate and obtain the vibration velocity vector of any point in the sound field, breaking through the traditional method of calculating the vibration velocity by differentially calculating the sound pressure derivative. At the same time, it can avoid the numerical error of the finite difference method and improve the prediction accuracy of the vibration velocity vector, thereby improving the application technology level of the vector hydrophone.

In step S1 of this embodiment, the on-site measurement data and sound source parameter information of the marine area to be predicted are first obtained, wherein the on-site measurement data of the marine area includes data such as ocean depth, sound speed, density (in this embodiment, the sea depth is 1000 meters, The speed of sound in water is uniformly 1500m/s, and the density of the entire field is 1.0g/cm ³ ), and the sound source parameter information includes data such as sound source frequency and location (in this embodiment, the sound source frequency is 100 Hz and the depth is 25 meters). The hydroacoustic Helmholtz equation of the cylindrical coordinate system in the horizontal layered ocean environment is established, and then Hankel transform is performed on the hydroacoustic Helmholtz equation of the cylindrical coordinate system to obtain the depth equation with the sound pressure wavenumber kernel function as the variable.

In this embodiment, the specific steps of step S1 include:

S11. Establish the hydroacoustic Helmholtz equation of the cylindrical coordinate system according to the following formula:

Among them, P(r,z) is the relative sound pressure in the frequency domain, ρ is the density of the sound propagation medium, k is the wave number and k=2πf/c, f is the frequency of the sound source, c is the sound velocity of the medium, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, z _s is the depth of the sound source, and δ is the Dirac function;

S12. Divide the underwater acoustic propagation medium into N layers in the depth direction, and approximate each layer as a homogeneous medium; in each divided layer, perform Hankel transformation on the hydroacoustic Helmholtz equation in the cylindrical coordinate system to convert (r ,z) The sound pressure P(r,z) in the space is converted to the (k _r ,z) space, that is:

Among them, φ(k _r , z) is the sound pressure kernel function, k _r is the horizontal wave number, and J ₀ is the Bessel function;

可得到深度方程为：Simultaneous Integration of Both Sides of Cylindrical Hydroacoustic Helmholtz Equation

The depth equation can be obtained as:

得到以声压波数核函数为变量的深度方程(3)。In this embodiment, the underwater acoustic transmission medium is divided into N=100 layers in the vertical direction, the thickness of each layer is dz=1m, and each layer is a homogeneous medium. In each layer, the cylindrical coordinate system underwater acoustic Helmholtz The equation is Hankel transformed according to the above formula (2), and both sides of the cylindrical coordinate hydroacoustic Helmholtz equation are integrated at the same time.

The depth equation (3) is obtained with the sound pressure wavenumber kernel function as the variable.

When solving the depth equation in step S2 of this embodiment, it includes the step of forming a joint vector between the sound pressure kernel function and the vertical vibration velocity kernel function, and the specific steps include:

S211. In any medium layer without sound source, the general solution form of the sound pressure kernel function is:

指数项

分别具有向上、下传播的物理意义，A⁺(k_r)表示向下传播的项，A^-(k_r)表示向上传播的项，w(k_r,z)为z方向(垂直方向)的垂直振速核函数，且满足：where φ(k _r ,z) is the sound pressure kernel function, k _z is the vertical wave number and

index term

They have the physical meaning of upward and downward propagation, respectively, A ⁺ (k _r ) represents the item of downward propagation, A ^- (k _r ) represents the item of upward propagation, and w(k _r , z) is the z direction (vertical direction) vertical vibration velocity kernel function, and satisfy:

Among them, Γ=k _z /(ρω), ρ is the density of the sound propagation medium, ω=2πf is the angular frequency of the sound source vibration, and f is the sound source vibration frequency;

The sound pressure kernel function and the vertical vibration velocity kernel function form a joint vector as:

S212. According to the joint vector formed in step S211, the joint vector at the upper interface z _m-1 of the mth layer is obtained as:

And the joint vector at the lower interface z _m of the mth layer is:

后，得到联合矢量由下至上传递公式为：S213. Combine and eliminate the joint vector expressions of the upper and lower interfaces of the mth layer obtained in step S212

After that, the bottom-up transfer formula of the joint vector is obtained as:

v _m (k _r ,z _m-1 )=M _m (k _r )v _m (k _r ,z _m ) (9)

Among them, M _m (k _r ) is the transfer matrix of the m-th layer medium from bottom to top. If the thickness of the m-th layer is h _m =z _m -z _m-1 , the expression of M _m (k _r ) is:

And the transfer formula and transfer matrix from top to bottom are:

According to the above steps, a joint vector composed of the sound pressure kernel function and the vertical vibration velocity kernel function can be established, and the bottom-up and top-down transfer formulas of the joint vector can be obtained at the same time, which is convenient for the subsequent solution of the depth equation.

In this embodiment, the depth equation is solved in step S2, and the specific steps of obtaining the sound pressure kernel function and the vertical vibration velocity kernel function of the interface of each medium layer include:

使用下标“0”表示上边界，得到上边界处的声压核函数φ₀(k_r,z₀)、垂直振速核函数w₀(k_r,z₀)分别为：S221. According to the boundary surface acoustic energy on the computational domain, only the items that can propagate upward and propagate downward are zero, that is,

Using the subscript "0" to represent the upper boundary, the sound pressure kernel function φ ₀ (k _r , z ₀ ) and the vertical vibration velocity kernel function w ₀ (k _r , z ₀ ) at the upper boundary are obtained as:

and the joint vector at the upper boundary is:

本实施例具体

ρ₀表示海面上方介质密度，在绝对软(压力释放)边界条件下，海面上方可视为真空，即ρ₀＝0g/cm³；Among them, w ₀ (k _r , z ₀ ) is the vertical vibration velocity kernel function of the upper boundary, and the vector (1, B ₀ ) ^T is the sound vector of the upper boundary, where

This embodiment is specific

ρ ₀ represents the density of the medium above the sea surface, and under the absolute soft (pressure release) boundary condition, the above sea surface can be regarded as a vacuum, that is, ρ ₀ =0g/cm ³ ;

使用下标“N”表示下边界，得到下边界处的联合矢量为：In the same way, according to the boundary surface acoustic energy under the computational domain, only the items that can propagate downward and propagate upward are zero, that is,

Using the subscript "N" to denote the lower boundary, the joint vector at the lower boundary is:

本实施例具体

ρ_N+1分别表示海底下方介质密度，在绝对硬边界条件下，海底下方可视为密度非常大的岩石，即取ρ_N+1＝10⁹⁹g/cm³；Among them, w _N (k _r , z _N ) is the vertical vibration velocity kernel function of the lower boundary, φ _N (k _r , z _N ) is the sound pressure kernel function of the lower boundary, (1, B _N ) ^T is the lower boundary sound vector,

This embodiment is specific

ρ _N+1 represents the density of the medium under the seabed respectively. Under the absolute hard boundary condition, the bottom of the seabed can be regarded as a very dense rock, that is, ρ _N+1 = 10 ⁹⁹ g/cm ³ ;

S222. Transfer the sound vector from the lower boundary _zN to the sound source depth _zs , that is, use the vector bottom-up transfer formula, start from the lower boundary _zN , and transfer the calculated sound vector layer by layer until the sound source depth _zs , the specific z _N =100m, z _s =25m, where the joint vector calculation formula immediately below the sound source depth is:

In the formula, w _N (k _r , z _N ) is the undetermined lower boundary vertical vibration velocity kernel function, and the sound vector is immediately below the sound source depth

S223. Transfer the sound vector from the upper boundary z ₀ to the sound source depth z _s , specifically z ₀ =0m, where the joint vector calculation formula obtained by the transfer immediately above the sound source depth is:

In the formula, w ₀ (k _r , z ₀ ) is the undetermined upper boundary vertical vibration velocity kernel function, and the sound source depth is just above the sound vector

S224. Calculate the vertical vibration velocity kernel functions of the lower and upper boundaries according to the sound source interface conditions, which satisfy at the sound source interface:

That is, v _s+1 (k _r ,z _s )-v _s (k _r ,z _s )=Δv(k _r ,z _s ), expand to get:

Then, the lower and upper boundary vertical vibration velocity kernel functions are obtained respectively as:

_S225 . _Calculate the sound pressure _kernel function _{φ ( k r} , z), the value of the vertical vibration velocity kernel function w(k _r , z).

Through the above steps, the transfer function matrix method is used to first transfer the sound vector from the lower boundary upward layer by layer to solve the sound vector of each layer until the sound vector immediately below the sound source interface is solved; The sound vector is obtained until the sound vector immediately above the sound source interface is obtained; at the sound source interface, the joint vector equations immediately below and immediately above are established using the sound source conditions, and the vertical vibration velocities of the lower and upper boundaries are solved; The vertical vibration velocity of the boundary is multiplied by the sound vector stored at the interface of each layer to obtain the sound pressure kernel function and the vertical vibration velocity kernel function.

In this embodiment, the sound pressure kernel function obtained by the interface of each layer is used first, and the sound pressure is solved according to the integral formula of Hankel inverse transformation; The vertical vibration velocity is solved by the directional derivative; then the horizontal vibration velocity is solved according to the horizontal direction derivative of the Hankel inverse transformation integral formula by converting the derivative of the zero-order Bessel function into a first-order Bessel function, as shown in the following details.

When solving the sound pressure in step S3 of this embodiment, the specific integral formula of Hankel inverse transformation is:

The discrete calculation formula of sound pressure obtained by discretizing the horizontal wavenumber k _r in the inverse Hankel transform integral formula is:

Among them, Δk _r is the horizontal wave number step size and Δk _r =2π/(r _max n _w ), r _max is the maximum horizontal distance of the sound field (r _max =3000m in this embodiment), and n _w is the Bessel function in a 2π oscillation period The minimum number of sampling points (n _w =10 in this embodiment), k _r,n is a discrete horizontal wave number and k _r,n =nΔk _r -iε _k , i is an imaginary unit, and ε _k is a complex offset to use In order to prevent singularity in the solution process of the depth equation, ε _k =3Δk _r /(2πlog ₁₀ e), M is the maximum index number of the discrete horizontal wave number and M=km _max /Δk _r , km _max is the cutoff wave number, this embodiment is specific Take km _max =20k ₀ , where the seawater wave number k ₀ =2πf/1500. Using the discrete calculation formula of the sound pressure at each position (r, z) of the sound field to perform horizontal wavenumber integration calculation, the sound pressure corresponding to each position (r, z) is obtained.

When calculating the vibration velocity vector in step S3 of the present embodiment, the specific steps of solving the vertical vibration velocity include:

S311. Calculate the derivative of the integral formula of Hankel's inverse transformation in the vertical z direction, and obtain the first transformation formula as:

Among them, P(r, z) is the relative sound pressure in the frequency domain, φ(k _r , z) is the sound pressure kernel function, k _r is the horizontal wave number, r is the coordinate in the horizontal direction, and z is the coordinate in the vertical or depth direction. , w(k _r , z) is the vertical vibration velocity kernel function, J ₀ is the zero-order Bessel function, ρ is the sound propagation medium density, ω=2πf is the sound source vibration angular frequency, and f is the sound source vibration frequency;

S312. The vertical vibration velocity component V _z (r, z) in the vibration velocity vector V(r, z)=(V _r (r, z), V _z (r, z)) is performed according to the first transformation formula The integral is calculated to get:

And the horizontal wavenumber k _r is discretized by the same method as the above sound pressure calculation, and the discrete integral calculation formula of the vertical vibration velocity is obtained as:

S313. Use the discrete integral calculation formula of vertical vibration velocity at each position (r, z) of the sound field to perform horizontal wavenumber integral calculation to obtain the vertical vibration velocity at each position (r, z).

When calculating the vibration velocity vector in step S3 of the present embodiment, the step of solving the horizontal vibration velocity includes:

得到第二变换式为：S321. Calculate the derivative of the Hankel inverse transform integral in the horizontal r direction, and convert the derivative of the zero-order Bessel function into a first-order Bessel function, namely

The second transformation is obtained as:

Among them, P(r, z) is the relative sound pressure in the frequency domain, φ(k _r , z) is the sound pressure kernel function, k _r is the horizontal wave number, r is the coordinate in the horizontal direction, and z is the coordinate in the vertical or depth direction. , J ₀ is a zero-order Bessel function, and J ₁ is a first-order Bessel function;

S322. Perform integral calculation according to the second transformation formula to obtain the horizontal vibration velocity component V _r (r, z) in the vibration velocity vector as:

Among them, ρ is the density of the sound propagation medium, ω=2πf is the angular frequency of the sound source vibration, and f is the sound source vibration frequency;

And the horizontal wavenumber k _r is discretized by the same method as the above sound pressure calculation, and the discrete integral calculation formula of the horizontal vibration velocity is obtained as:

S323. Use the discrete integral calculation formula of the horizontal vibration velocity at each position (r, z) of the sound field to perform the horizontal wavenumber integral calculation to obtain the horizontal vibration velocity corresponding to each position (r, z).

Through the above steps, the integral formulas of the vertical vibration velocity and the horizontal vibration velocity are obtained by derivating the integral formula of the Hankel inverse transform in the vertical and horizontal directions, and the horizontal wavenumber can be discretized by the same method as the integral formula of sound pressure, and the solution can be obtained. Vertical vibration speed and horizontal vibration speed, so as to realize the prediction of vibration speed vector V(r,z)=( _Vr (r,z), _Vz (r,z)).

In step S4 of this embodiment, the sound pressure calculation propagation loss obtained in step S3 is further included to form a scalar cloud map of sound field propagation loss. The calculation formula for calculating the propagation loss (scalar) according to the sound pressure is specifically:

TL(r,z)=-20log ₁₀ |P(r,z)| (30)

And according to the calculated sound pressure P(r, z), the vibration velocity vector V(r, z) = (V _r (r, z), V _z (r, z)) Calculate each position of the sound field according to formula (31) (r,z) time-averaged sound intensity vector;

分别表示复数声压P(r,z)、振速V(r,z)的共轭复数。Among them, I(r,z) is the time-averaged sound intensity vector, that is, on the plane perpendicular to the direction of the vibration velocity, the time-averaged sound energy passing through a unit area,

They represent the complex conjugates of the complex sound pressure P(r,z) and the vibration velocity V(r,z), respectively.

The scalar cloud diagram of the sound field propagation loss obtained in the specific application example is shown in Figure 3, and the sound field sound intensity vector streamline (the background color is the propagation loss) obtained by the calculation is shown in Figure 4. It can be seen from Figure 4 that the propagation loss The sound intensity vector lines of small (high energy) regions are relatively dense, and the sound intensity vector lines of regions with large propagation loss (low energy) are relatively sparse.

As shown in FIG. 2 , when the above method of this embodiment is applied in a specific application embodiment, in combination with a high-performance computer workstation, the data storage medium of the high-performance computer workstation is loaded with functions capable of realizing the above-mentioned method for predicting the ocean vector sound field of this embodiment. The high-performance computer workstation receives on-site measurement data including ocean depth, sound speed, density, etc., as well as sound source parameter information including sound source frequency, location, etc. After the above method and steps for predicting ocean vector sound field, the ocean vector sound field is generated. The forecast graphic image can be further provided for the subsequent processing and analysis of the vector hydrophone acoustic signal.

This embodiment also provides a system for predicting an ocean vector sound field, including a processor and a memory, where the memory is used to store marine environment data, sound source parameters and a computer program, the processor is used to execute the computer program, and the processor is used to execute the computer program to Execute the above method for predicting the ocean vector sound field. The system of this embodiment may specifically adopt the structure shown in FIG. 2 , and the high-performance computer workstation is configured with the above-mentioned processor and memory.

In addition, the present embodiment also provides a computer-readable storage medium storing marine environment data, sound source parameters, and a computer program programmed or configured to execute the above-mentioned method for improving and forecasting an oceanic vector sound field.

As will be appreciated by those skilled in the art, the embodiments of the present application may be provided as a method, a system, or a computer program product. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) having computer-usable program code embodied therein. The present application refers to flowcharts of methods, apparatus (systems), and computer program products according to embodiments of the present application and/or processor-executed instructions generated for implementing a process or processes and/or block diagrams in a flowchart. A means for the function specified in a block or blocks. These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory result in an article of manufacture comprising instruction means, the instructions The apparatus implements the functions specified in the flow or flow of the flowcharts and/or the block or blocks of the block diagrams. These computer program instructions can also be loaded on a computer or other programmable data processing device to cause a series of operational steps to be performed on the computer or other programmable device to produce a computer-implemented process such that The instructions provide steps for implementing the functions specified in the flow or blocks of the flowcharts and/or the block or blocks of the block diagrams.

The above are only preferred embodiments of the present invention, and do not limit the present invention in any form. Although the present invention has been disclosed above with preferred embodiments, it is not intended to limit the present invention. Therefore, any simple modifications, equivalent changes and modifications made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solutions of the present invention should fall within the protection scope of the technical solutions of the present invention.

Claims (4)

1. A method of predicting a marine vector sound field, the steps comprising:

s1, acquiring field measurement data and sound source parameter information of a to-be-forecasted marine area, establishing a cylindrical coordinate system underwater sound Helmholtz equation under a horizontal layered marine environment, and obtaining a depth equation taking a sound pressure kernel function as a variable after integral transformation;

s2, solving the depth equation by adopting a transfer function matrix method to obtain a sound pressure kernel function and a vertical vibration velocity kernel function of each medium layer interface;

s3, according to the sound pressure kernel function, using a Hankel inverse transformation integral type to calculate sound pressure, according to the vertical vibration velocity kernel function, using a vertical direction derivative based on the Hankel inverse transformation integral type to calculate vertical vibration velocity, and based on the Hankel inverse transformation integral type horizontal direction derivative, converting a derivative of a zero-order Bessel function into a first-order Bessel function to obtain a horizontal vibration velocity integral type so as to calculate horizontal vibration velocity, and obtaining a vibration velocity vector;

s4, calculating underwater sound propagation loss and sound intensity vectors of the ocean to be forecasted according to the sound pressure and vibration velocity vectors obtained in the step S3;

the step of step S1 includes:

s11, establishing an underwater sound Helmholtz equation of the cylindrical coordinate system according to the following formula:

wherein P (r, z) is relative sound pressure in frequency domain, ρ is acoustic propagation medium density, k is wave number and k is 2 π f/c, f is sound source frequency, c is medium sound velocity, r is coordinate in horizontal direction, z is coordinate in vertical or depth direction, P (r, z) is relative sound pressure in frequency domain, ρ is acoustic propagation medium density, k is wave number and k is 2 π f/c, f is sound source frequency, c is medium sound velocity, r is coordinate in horizontal direction, z is coordinate in vertical or depth direction, z is acoustic propagation medium density, and z is acoustic medium density_sIs the sound source depth, delta is the dirac function;

s12, dividing the underwater sound transmission medium into N layers in the depth direction, and enabling each layer to be similar to a uniform medium; within each layer of division, Hankel transformation is carried out on the cylindrical coordinate system underwater sound Helmholtz equation to convert the sound pressure P (r, z) of the (r, z) space into (k)_rZ) space, i.e.:

wherein phi (k)_rZ) is the sound pressure kernel, k_rIs a horizontal wavenumber, J₀Is a Bessel function;

integrating the two sides of the cylindrical coordinate underwater sound Helmholtz equation simultaneously

The depth equation can be found as:

when the depth equation is solved in step S2, the method includes a step of forming a combined vector by using a sound pressure kernel function and a vertical vibration velocity kernel function, and includes the specific steps of:

s211, in any medium layer without a sound source, the general solution form of the sound pressure kernel function is as follows:

wherein phi (k)_rZ) is the sound pressure kernel, k_zIs a vertical wavenumber and

A⁺(k_r) Representing a downwardly propagating item, A^-(k_r) Representing an upwardly propagating term, w (k)_rAnd z) is a vertical vibration velocity kernel function in the z direction and satisfies the following conditions:

wherein Γ ═ k_z(ρ ω), ρ is the sound propagation medium density, ω ═ 2 π f is the sound source vibration angular frequency, and f is the sound source vibration frequency;

forming a joint vector by the sound pressure kernel function and the vertical vibration velocity kernel function as follows:

s212, obtaining an upper interface z of the mth layer according to the combined vector formed in the step S211_m-1The joint vector of (a) is:

and the lower interface z of the mth layer_mThe joint vector of (a) is:

s213, combining and eliminating the upper and lower interface joint vector expressions of the mth layer obtained in the step S212

And then, obtaining a transfer formula of the joint vector from bottom to top as follows:

v_m(k_r,z_m-1)＝M_m(k_r)v_m(k_r,z_m)

wherein M is_m(k_r) Is a transfer matrix of the mth layer medium from bottom to top, if the thickness of the mth layer is h_m＝z_m-z_m-1Then M is_m(k_r) The expression of (a) is:

and the transfer formula and the transfer matrix from top to bottom are:

the specific steps of solving the depth equation in step S2 to obtain the sound pressure kernel and the vertical vibration velocity kernel of the interface of each dielectric layer include:

s221, according to the fact that the sound energy of the boundary on the calculation domain can only be upwards propagated and downwards propagated, the term is zero, namely

Using subscript "0" to denote the upper boundary, the sound pressure kernel at the upper boundary is obtained₀(k_r,z₀) Hang downKernel function w of direct vibration velocity₀(k_r,z₀) Respectively as follows:

and the joint vector at the upper boundary is:

wherein, w₀(k_r,z₀) Vertical vibration velocity kernel function of upper boundary, vector (1, B)₀)^TIs an upper bound acoustic vector, wherein

The term of downward propagation and upward propagation is zero according to the fact that the boundary acoustic energy under the calculation domain is only possible to be downward propagation, namely

Using the subscript "N" to denote the lower boundary, the joint vector at the lower boundary is found to be:

wherein, w_N(k_r,z_N) Is the vertical vibration velocity kernel function of the lower boundary, phi_N(k_r,z_N) Is the sound pressure kernel function of the lower boundary (1, B)_N)^TIs a lower boundary acoustic vector, wherein

S222, acoustic vector is divided from lower boundary z_NTo the depth z of the sound source_sPassing, i.e. from the lower boundary z_NStarting, transmitting and calculating acoustic vector layer by layer until the depth z of the acoustic source_sWherein the joint vector calculation formula immediately below the sound source depth is:

in the formula, w_N(k_r,z_N) For the undetermined lower boundary vertical vibration velocity kernel function, the sound vector immediately below the sound source depth

S223, enabling the sound vector to be in the upper boundary z₀To the depth z of the sound source_sAnd (3) transferring, wherein a joint vector calculation formula obtained by transferring immediately above the sound source depth is as follows:

in the formula, w₀(k_r,z₀) For the kernel function of the vertical vibration velocity of the upper boundary to be determined, the sound vector just above the sound source depth

S224, calculating a lower boundary and an upper boundary vertical vibration velocity kernel function according to the sound source interface conditions, wherein the sound source interface conditions are as follows:

i.e. v_s+1(k_r,z_s)-v_s(k_r,z_s)＝Δv(k_r,z_s) And unfolding to obtain:

then solving to obtain the kernel functions of the vertical vibration velocity of the lower and upper boundaries as follows:

s225, according to the calculated lower and upper boundary vertical vibration velocity kernel function w_N+1(k_r,z_N)、w₀(k_r,z₀) Calculating sound pressure kernel function phi (k) of each dielectric layer interface position_rZ) vertical vibration velocity kernel function w (k)_rZ) value;

when sound pressure is solved in the step S3, the specific Hankel inverse transformation integral formula is as follows:

converting the horizontal wave number k in the integral formula of the Hankel inverse transformation_rThe sound pressure dispersion calculation formula obtained after dispersion is as follows:

wherein, Δ k_rIs a horizontal wave number step size and Δ k_r＝2π/(r_maxn_w)，r_maxMaximum horizontal distance of sound field, n_wFor the minimum number of sampling points, k, of the Bessel function in a 2 pi oscillation period_r,nIs a discrete horizontal wavenumber and k_r,n＝nΔk_r-iε_kI is an imaginary unit, ε_kIs a complex offset of_k＝3Δk_r/(2πlog₁₀e) M is the maximum index number of the discrete horizontal wavenumber and M ═ k_max/Δk_r，k_maxIs the cut-off wave number;

performing horizontal wave number integral calculation on each position (r, z) of the sound field by adopting the discrete calculation formula of the sound pressure to obtain the sound pressure corresponding to each position (r, z);

when the vibration velocity vector is calculated in step S3, the step of solving the vertical vibration velocity includes:

s311, solving a derivative of the Hankel inverse transformation integral expression in the vertical z direction to obtain a first transformation expression as follows:

wherein P (r, z) is relative sound pressure in frequency domain, phi (k)_rZ) is the sound pressure kernel, k_rIs the horizontal wave number, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, w (k)_rZ) is the vertical vibration velocity kernel function, J₀The acoustic source is a zero-order Bessel function, rho is the density of an acoustic propagation medium, omega-2 pi f is the vibration angular frequency of the acoustic source, and f is the vibration frequency of the acoustic source;

s312. the vibration velocity vector V (r, z) of the mass point is (V)_r(r,z)，V_z(r, z)) of the vertical vibration velocity component V_z(r, z) performing integral calculation according to the first transformation formula to obtain:

and for horizontal wave number k_rAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the vertical vibration velocity as follows:

s313, performing horizontal wave number integral calculation on each position (r, z) of the sound field by using the discrete integral calculation formula of the vertical vibration velocity to obtain the vertical vibration velocity corresponding to each position (r, z);

when the vibration velocity vector is calculated in step S3, the specific step of solving the horizontal vibration velocity includes:

s321, solving a derivative of a Hankel inverse transformation integral expression in the horizontal r direction, and converting the derivative of a zero-order Bessel function into a first-order Bessel function, namely the first-order Bessel function

The second transformation is obtained as:

wherein P (r, z) is relative sound pressure in frequency domain, phi (k)_rZ) is the sound pressure kernel, k_rIs the horizontal wave number, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, J₀Is a Bessel function of the zero order,

is a first order Bessel function;

s322, carrying out integral calculation according to the second transformation expression to obtain a horizontal vibration velocity component V in the vibration velocity vector_r(r, z) is:

wherein ρ is the density of the sound propagation medium, ω ═ 2 π f is the vibration angular frequency of the sound source, and f is the vibration frequency of the sound source;

and for horizontal wave number k_rAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the horizontal vibration velocity as follows:

s323, performing horizontal wave number integral calculation on each position (r, z) of the sound field by adopting the discrete integral calculation formula of the horizontal vibration speed to obtain the horizontal vibration speed corresponding to each position (r, z);

in step S3, the calculation formula for calculating the propagation loss from the sound pressure is specifically:

TL(r,z)＝-20log₁₀|P(r,z)|

and according to the solved sound pressure P (r, z) and vibration velocity vector V (r, z) ═ V_r(r,z),V_z(r, z)) calculating a time-averaged intensity vector of each position (r, z) of the sound field according to the following formula;

wherein I (r, z) is a time-averaged sound intensity vector, i.e. the time-averaged value of the sound energy passing through the unit area on a plane perpendicular to the vibration velocity direction,

respectively, the complex sound pressure P (r, z) and the complex conjugate of the vibration velocity V (r, z).

2. A method as claimed in claim 1, wherein the step S4 further includes calculating propagation loss according to the sound pressure solved in step S3, and forming a sound field propagation loss scalar cloud.

3. A system for forecasting a marine vector sound field, comprising a processor and a memory for storing marine environmental data, sound source parameters and a computer program, wherein the processor is adapted to execute the computer program to perform the method according to any of claims 1-2.

4. A computer readable storage medium having stored thereon marine environmental data, acoustic source parameters, and a computer program programmed or configured to perform the method of predicting a marine vector sound field according to any one of claims 1-2.

Images