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

Method, system and medium for forecasting ocean vector sound field Download PDF

Info

Publication number
CN112254798B
CN112254798B CN202011083604.9A CN202011083604A CN112254798B CN 112254798 B CN112254798 B CN 112254798B CN 202011083604 A CN202011083604 A CN 202011083604A CN 112254798 B CN112254798 B CN 112254798B
Authority
CN
China
Prior art keywords
sound
vibration velocity
vector
sound pressure
horizontal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011083604.9A
Other languages
Chinese (zh)
Other versions
CN112254798A (en
Inventor
刘巍
王勇献
张理论
程兴华
肖汶斌
马树青
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202011083604.9A priority Critical patent/CN112254798B/en
Publication of CN112254798A publication Critical patent/CN112254798A/en
Application granted granted Critical
Publication of CN112254798B publication Critical patent/CN112254798B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H3/00—Measuring characteristics of vibrations by using a detector in a fluid
    • G01H3/04—Frequency
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H3/00—Measuring characteristics of vibrations by using a detector in a fluid
    • G01H3/10—Amplitude; Power
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

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

Method, system and medium for forecasting 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 forecasting an ocean vector sound field.
Background
The sound wave is the main means of underwater communication, marine environment and target detection at present, and has important application value in the military and economic fields. The receiving sensor of the underwater sound wave is generally called a hydrophone, and the traditional hydrophone is a scalar hydrophone which can only measure scalar parameters (such as sound pressure) in a sound field. The current advanced hydrophone is a vector hydrophone, can measure scalar parameters and vector parameters (such as particle vibration speed, vibration speed for short) in a sound field, and has important significance for improving the performance of an underwater detection system.
Because factors such as submarine topography, marine environment (mainly sound velocity and density and having time-varying characteristics), sound source frequency and position and the like all have important influence on underwater sound propagation, in the processes of hydrophone layout and site selection, array shape optimization, received signal processing and analysis and the like, a marine sound field needs to be forecasted, namely, sound field numerical simulation is carried out by combining the dynamically-changing marine environment, and the forecast result is fused into the hydrophone signal processing process, so that the marine sound field forecasting has important effects on improving the application level of the hydrophone and enhancing the performance of an underwater detection system.
However, the traditional underwater acoustic model (such as wave number integration method, normal wave method, etc.) can only predict a scalar sound field, and if a vector sound field is to be obtained, dense grid points need to be arranged in the sound field, after sound pressures of various points are solved, a finite difference method is adopted to calculate a sound pressure derivative of each point, and finally the sound pressure derivative is converted into a vibration velocity. The above conventional method for calculating the sound field vector (vibration velocity) by using the finite difference method has the following disadvantages:
(1) the sound field grid spacing is limited by finite difference methods. 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 greatly in the interval of adjacent grid points), the grid spacing in each direction generally needs to be a fraction (for example, a tenth) of the reference wavelength, and in the case that only a few spatial positions around the hydrophone need to be calculated, the limitation condition will have an adverse effect on the calculation amount of the sound field vector and the flexibility of the calculation method.
(2) The finite difference method of calculating the sound pressure derivative introduces additional numerical errors. In practical applications, a finite difference method with a finite order (such as second-order precision) is generally adopted, and when the grid spacing is finite (not approaching zero), the process of calculating derivatives by difference will introduce the dissipation and dispersion errors of the finite difference method.
Disclosure of Invention
The technical problem to be solved by the invention is as follows: aiming at the technical problems in the prior art, the invention provides a method, a system and a medium for forecasting an ocean vector sound field, which can realize the forecasting of a vibration velocity vector based on ocean environment measurement data, and can improve the precision of the vibration velocity vector forecasting based on a vibration velocity calculation method for directly performing horizontal wave number integration.
In order to solve the technical problems, the technical scheme provided by the invention is as follows:
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 with 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, calculating sound pressure by using a Hankel inverse transformation integral formula according to the sound pressure kernel function, calculating a vertical vibration velocity by using a vertical derivative based on the Hankel inverse transformation integral formula according to the vertical vibration velocity kernel function, and converting a derivative of a zeroth-order Bessel function into a first-order Bessel function to obtain a horizontal vibration velocity integral formula so as to calculate a horizontal vibration velocity and obtain a vibration velocity vector;
and S4, calculating the underwater sound propagation loss and the sound intensity vector of the ocean to be forecasted according to the sound pressure and the vibration velocity vector solved in the step S3.
Further, the step of step S1 includes:
s11, establishing an underwater sound Helmholtz equation of the cylindrical coordinate system according to the following formula:
Figure BDA0002719530090000021
wherein P (r, z) is relative sound pressure in frequency domain, ρ is the density of the sound propagation medium, k is wave number and k is 2 π f/c,f is the sound source frequency, c is the medium sound velocity, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, z issIs 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.:
Figure BDA0002719530090000022
wherein phi (k)rZ) is the sound pressure kernel, krIs a horizontal wavenumber, J0Is a Bessel function;
integrating the two sides of the cylindrical coordinate underwater sound Helmholtz equation simultaneously
Figure BDA0002719530090000023
The depth equation can be found as:
Figure BDA0002719530090000024
further, when the depth equation is solved in step S2, the method includes a step of forming a joint vector by using a sound pressure kernel function and a vertical vibration velocity kernel function, and the specific steps include:
s211, in any medium layer without a sound source, the general solution form of the sound pressure kernel function is as follows:
Figure BDA0002719530090000031
wherein phi (k)rZ) is the sound pressure kernel, kzIs a vertical wavenumber and
Figure BDA0002719530090000032
A+(kr) Representing a downwardly propagating item, A-(kr) 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:
Figure BDA0002719530090000033
wherein Γ ═ kz(ρ ω), ρ is the sound propagation medium density, ω ═ 2 π f is the sound source vibration angular frequency, and f is the sound source vibration frequency;
then the sound pressure kernel function and the vertical vibration velocity kernel function form a joint vector as follows:
Figure BDA0002719530090000034
s212, obtaining an upper interface z of the mth layer according to the combined vector formed in the step S211m-1The joint vector of (a) is:
Figure BDA0002719530090000035
and the lower interface z of the mth layermThe joint vector of (a) is:
Figure BDA0002719530090000036
s213, combining and eliminating the upper and lower interface joint vector expressions of the mth layer obtained in the step S212
Figure BDA0002719530090000037
And then, obtaining a transfer formula of the joint vector from bottom to top as follows:
vm(kr,zm-1)＝Mm(kr)vm(kr,zm)
wherein, Mm(kr) Is a transfer matrix of the mth layer medium from bottom to top, if the thickness of the mth layer is hm＝zm-zm-1Then M ism(kr) The expression of (c) is:
Figure BDA0002719530090000038
and the transfer formula and the transfer matrix from top to bottom are:
Figure BDA0002719530090000039
Figure BDA00027195300900000310
further, the step S2 of obtaining the sound pressure kernel function and the vertical vibration velocity kernel function of the interface of each dielectric layer specifically includes:
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
Figure BDA0002719530090000041
The upper boundary is denoted by subscript "0", and the acoustic pressure kernel function φ at the upper boundary is obtained0(kr,z0) Vertical vibration velocity kernel function w0(kr,z0) Respectively as follows:
Figure BDA0002719530090000042
Figure BDA0002719530090000043
and the joint vector at the upper boundary is:
Figure BDA0002719530090000044
wherein w0(kr,z0) Vertical velocity kernel function as the upper boundary, vector (1, B)0)TIs an upper bound acoustic vector, wherein
Figure BDA0002719530090000045
The term which is only possible to propagate downwards and upwards according to the interface acoustic energy under the calculation domain is zero, namely
Figure BDA0002719530090000046
Using the subscript "N" to denote the lower boundary, the joint vector at the lower boundary is found to be:
Figure BDA0002719530090000047
wherein, wN(kr,zN) Is the vertical vibration velocity kernel function of the lower boundary, phiN(kr,zN) Is the sound pressure kernel function of the lower boundary (1, B)N)TIs a lower boundary acoustic vector, wherein
Figure BDA0002719530090000048
S222, acoustic vector is divided from lower boundary zNTo the depth z of the sound sourcesPassing, i.e. from the lower boundary zNStarting, transmitting and calculating acoustic vector layer by layer until the depth z of the acoustic sourcesWherein the joint vector calculation formula immediately below the sound source depth is:
Figure BDA0002719530090000049
in the formula, wN(kr,zN) For the undetermined lower boundary vertical vibration velocity kernel function, the sound vector immediately below the sound source depth
Figure BDA00027195300900000410
S223, enabling the sound vector to be in the upper boundary z0To a sound source depth zsAnd (3) transferring, wherein a joint vector calculation formula obtained by transferring immediately above the sound source depth is as follows:
Figure BDA0002719530090000051
in the formula, w0(kr,z0) 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
Figure BDA0002719530090000052
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:
Figure BDA0002719530090000053
i.e. vs+1(kr,zs)-vs(kr,zs)＝Δv(kr,zs) And unfolding to obtain:
Figure BDA0002719530090000054
then solving to obtain the kernel functions of the vertical vibration velocity of the lower and upper boundaries as follows:
Figure BDA0002719530090000055
Figure BDA0002719530090000056
s225, according to the calculated lower and upper boundary vertical vibration velocity kernel function wN+1(kr,zN) And w0(kr,z0) Calculating sound pressure kernel function phi (k) of each dielectric layer interface positionrZ) vertical vibration velocity kernel function w (k)rZ) value.
Further, when the sound pressure is solved in step S3, a specific Hankel inverse transformation integral equation is as follows:
Figure BDA0002719530090000057
the horizontal wave number k in the Hankel inverse transformation integral expression is usedrThe sound pressure dispersion calculation formula obtained after dispersion is as follows:
Figure BDA0002719530090000058
wherein, Δ krIs a horizontal wave number step size and Δ kr＝2π/(rmaxnw)，rmaxMaximum horizontal distance of sound field, nwFor the minimum number of sampling points, k, of the Bessel function in a 2 pi oscillation periodr,nIs a discrete horizontal wavenumber and kr,n＝nΔkr-iεkI is an imaginary unit, εkIs a complex offset of ∈k＝3Δkr/(2πlog10e) M is the maximum index number of the discrete horizontal wavenumber and M ═ kmax/Δkr，kmaxThe cut-off wavenumber. And 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).
Further, 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 depth z direction to obtain a first transformation expression as follows:
Figure BDA0002719530090000061
wherein P (r, z) is relative sound pressure in frequency domain, phi (k)rZ) is the sound pressure kernel, krIs 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 verticalKernel function of vibration velocity, J0The 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 particle is (V)r(r,z),Vz(r, z)) of the vertical vibration velocity component Vz(r, z) performing an integral calculation according to the first transformation equation to obtain:
Figure BDA0002719530090000062
and for horizontal wave number krAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the vertical vibration velocity as follows:
Figure BDA0002719530090000063
and 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 speed to obtain the vertical vibration speed corresponding to each position (r, z).
Further, 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 depth r direction, and converting the derivative of a zero-order Bessel function into a first-order Bessel function, namely
Figure BDA0002719530090000064
The second transformation is obtained as:
Figure BDA0002719530090000065
wherein P (r, z) is relative sound pressure in frequency domain, phi (k)rZ) is the sound pressure kernel, krIs the horizontal wave number, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, J0Is a Bessel function of the zero order,
Figure BDA0002719530090000066
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 vectorr(r, z) is:
Figure BDA0002719530090000067
wherein ρ is the density of the sound propagation medium, ω ═ 2 π f is the sound source vibration angular frequency, and f is the sound source vibration frequency;
and for horizontal wave number krAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the horizontal vibration velocity as follows:
Figure BDA0002719530090000068
and S323, performing horizontal wave number integral calculation on each position (r, z) of the sound field by using the discrete integral calculation formula of the horizontal vibration speed to obtain the horizontal vibration speed corresponding to each position (r, z).
Further, the step S4 includes calculating a propagation loss according to the sound pressure obtained in step S3, forming a sound field propagation loss scalar cloud, and obtaining a (V, z) vector V (r, z) from the obtained sound pressure P (r, z) and the vibration velocity vector V (r, z)r(r,z),Vz(r, z)) calculating a time-averaged sound intensity vector of each position (r, z) of the sound field according to the following formula;
Figure BDA0002719530090000071
wherein I (r, z) is a time-averaged sound intensity vector, i.e. the time average value of the sound energy passing through the unit area on a plane perpendicular to the vibration velocity direction,
Figure BDA0002719530090000072
respectively, a complex sound pressure P (r, z) and a complex conjugate of a vibration velocity V (r, z), where r is a coordinate in the horizontal direction and z is a coordinate in the vertical or depth direction.
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, the processor being adapted to execute the computer program and the processor being adapted to execute the computer program to perform the method described above.
A computer readable storage medium having stored thereon marine environmental data, sound source parameters, and a computer program programmed or configured to perform the above-described method of predicting a marine vector sound field.
Compared with the prior art, the invention has the advantages that: the invention discloses a method, a system and a medium for forecasting an ocean vector sound field, which are characterized in that after field measurement data and sound source parameter information of an ocean area to be forecasted are acquired, on the basis of calculating sound pressure by using a wave number integration method underwater acoustic model, a wave number integration formula of vertical vibration velocity and horizontal vibration velocity is obtained by differentiating a Hankel inverse transformation integration formula, the horizontal wave number can be directly dispersed and integrated to calculate the vibration velocity vector of any point of the sound field, the limitation of the traditional method on the distance between grid points of the sound field caused by calculating the vibration velocity by differentially calculating the sound pressure derivative is broken through, and meanwhile, the numerical error of a finite difference method can be avoided, so that the vibration velocity vector forecasting precision is improved.
Drawings
Fig. 1 is a detailed implementation flowchart of the method for forecasting an ocean vector sound field according to the embodiment.
Fig. 2 is a schematic structural diagram of a system for implementing a marine vector sound field forecasting method in a specific application embodiment of the present invention.
Fig. 3 is a sound field propagation loss scalar cloud diagram obtained in a specific application embodiment of the invention.
Fig. 4 is a sound intensity vector flow chart obtained in a specific application embodiment of the present invention.
Detailed Description
The invention is further described below with reference to the drawings and specific preferred embodiments of the description, without thereby limiting the scope of protection of the invention.
The embodiment is particularly applied to the ocean environment with absolute hard seabed constant-speed waveguideFor example, the sound field prediction includes the following parameters: uniform sea water density rhow＝1.0g/cm3The sound velocity of the water body is uniform cw1500m/s, seabed level and sea depth zN100m, 1m in the step length dz in the depth z direction, and r is the maximum solving distance in the r directionmax1000m, and a step Δ r of 1 m. Sound source frequency f 100Hz, sound source depth zs25m, upper boundary (sea surface z)00m) pressure release boundary condition (zero sound pressure), lower boundary (seabed z)N100m) takes the absolute hard boundary condition (the z-derivative of the sound pressure is zero).
As shown in fig. 1, the detailed steps of the method for forecasting the marine vector sound field of the present embodiment include: 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 with a sound pressure kernel function as a variable after integral transformation;
s2, solving a 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;
and S4, calculating the underwater sound propagation loss and the sound intensity vector of the ocean to be forecasted according to the sound pressure and the vibration velocity vector obtained in the step S3.
The acoustic pressure can be calculated by using a wave number integration method underwater acoustic model, meanwhile, on the basis of using a wave number integration method, a Hankel inverse transformation integral formula is used for derivation in the vertical direction, a vertical vibration velocity integral formula can be obtained, then the vertical vibration velocity can be solved, and the Bessel function is considered to have the characteristics:
Figure BDA0002719530090000081
wherein J0Is a Bessel function of the zero order,J′0is J0Derivative of (A), J1The first-order Bessel function is a zero-order Bessel function, namely the zero-order Bessel function can be converted into the first-order Bessel function, and on the basis of the Hankel inverse transformation integral expression, the horizontal vibration velocity integral expression can be obtained by utilizing the characteristics of the Bessel function after derivation in the horizontal direction, and further the horizontal vibration velocity can be solved. In the embodiment, after the field measurement data of the marine area to be forecasted and the sound source parameter information are acquired, the sound pressure is calculated by using a wave number integration method underwater sound model, meanwhile, the wave number integration formulas of the vertical vibration speed and the horizontal vibration speed are respectively obtained by deriving the Hankel inverse transformation integral formula, wherein the vertical vibration speed is calculated by using the vertical vibration speed kernel function integration, the horizontal vibration speed is solved by using the characteristics of the Bessel function and the Hankel inverse transformation integral formula horizontal direction derivative, and the vibration speed vector is obtained, so that the forecasting of underwater sound propagation loss and sound intensity vector can be realized. According to the method, the vibration velocity vector of any point of the sound field can be directly calculated and obtained based on the vibration velocity calculation mode of directly performing horizontal wave number integration, the limitation on the distance between grid points of the sound field caused by the fact that the vibration velocity is calculated by differentially calculating the sound pressure derivative in the traditional method is broken through, and meanwhile, the numerical error of the finite difference method can be avoided, so that the vibration velocity vector prediction precision is improved, and the application technical level of the vector hydrophone is further improved.
In step S1 of this embodiment, the field measurement data of the ocean area to be forecasted and the sound source parameter information are first obtained, where the field measurement data of the ocean area includes data such as ocean depth, sound velocity, and density (in this embodiment, the sea depth is 1000 meters, the sound velocity in water is 1500m/S, and the total field density is 1.0g/cm3) The sound source parameter information includes data such as sound source frequency and position (in this embodiment, the sound source frequency is 100Hz, and the depth is 25 meters), a cylindrical coordinate system underwater sound Helmholtz equation under the horizontal layered marine environment is established according to the acquired data, and then Hankel transformation is performed on the cylindrical coordinate system underwater sound Helmholtz equation to obtain a depth equation with the sound pressure wave number kernel function as a variable.
In this embodiment, the step S1 includes the following specific steps:
s11, establishing an underwater sound Helmholtz equation of a cylindrical coordinate system according to the following formula:
Figure BDA0002719530090000091
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 densitysIs 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 the 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.:
Figure BDA0002719530090000092
wherein phi (k)rZ) is the sound pressure kernel, krIs a horizontal wavenumber, J0Is a Bessel function;
integrating the two sides of cylindrical coordinate underwater sound Helmholtz equation simultaneously
Figure BDA0002719530090000093
The depth equation can be found as:
Figure BDA0002719530090000094
in this embodiment, an underwater sound propagation medium is divided into N100 layers in the vertical direction, each layer has a thickness dz of 1m and is a uniform medium in each layer, Hankel transformation is performed on the cylindrical coordinate system underwater sound Helmholtz equation according to the above formula (2), and integration is performed on both sides of the cylindrical coordinate underwater sound Helmholtz equation simultaneously
Figure BDA0002719530090000095
Obtaining depth with sound pressure wave number kernel function as variableDegree equation (3).
In this embodiment, 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 and a vertical vibration velocity kernel, and the specific steps include:
s211, in any medium layer without a sound source, the general solution form of the sound pressure kernel function is as follows:
Figure BDA0002719530090000096
wherein phi (k)rZ) is the sound pressure kernel, kzIs a vertical wavenumber and
Figure BDA0002719530090000097
index term
Figure BDA0002719530090000098
Respectively have the physical significance of propagating upwards and downwards, A+(kr) Representing a downwardly propagating item, A-(kr) Representing an upwardly propagating term, w (k)rAnd z) is a vertical vibration velocity kernel function in the z direction (vertical direction) and satisfies:
Figure BDA0002719530090000101
wherein Γ ═ kz(ρ ω), ρ is the sound propagation medium density, ω ═ 2 π f is the sound source vibration angular frequency, and f is the sound source vibration frequency;
the sound pressure kernel function and the vertical vibration velocity kernel function form a combined vector as follows:
Figure BDA0002719530090000102
s212, obtaining an upper interface z of the mth layer according to the combined vector formed in the step S211m-1The joint vector of (a) is:
Figure BDA0002719530090000103
and the lower interface z of the mth layermThe joint vector of (a) is:
Figure BDA0002719530090000104
s213, combining and eliminating the upper and lower interface joint vector expressions of the mth layer obtained in the step S212
Figure BDA0002719530090000105
Then, the transfer formula of the obtained joint vector from bottom to top is as follows:
vm(kr,zm-1)＝Mm(kr)vm(kr,zm) (9)
wherein M ism(kr) Is a transfer matrix of the mth layer medium from bottom to top, if the thickness of the mth layer is hm＝zm-zm-1Then M ism(kr) The expression of (a) is:
Figure BDA0002719530090000106
and the transfer formula and the transfer matrix from top to bottom are respectively:
Figure BDA0002719530090000107
according to the steps, a combined vector consisting of the sound pressure kernel function and the vertical vibration velocity kernel function can be established, and meanwhile, transfer formulas of the combined vector from bottom to top and from top to bottom are obtained, so that the depth equation can be conveniently solved in the follow-up process.
In this embodiment, the step S2 of solving the depth equation to obtain the sound pressure kernel function and the vertical vibration velocity kernel function of the interface of each dielectric layer includes:
s221, only the boundary acoustic energy on the calculation domain can be usedThe terms propagating upwards and downwards are zero, i.e.
Figure BDA0002719530090000108
The upper boundary is denoted by subscript "0", and the acoustic pressure kernel function φ at the upper boundary is obtained0(kr,z0) Vertical vibration velocity kernel function w0(kr,z0) Respectively as follows:
Figure BDA0002719530090000111
Figure BDA0002719530090000112
and the joint vector at the upper boundary is:
Figure BDA0002719530090000113
wherein, w0(kr,z0) Vertical vibration velocity kernel function of upper boundary, vector (1, B)0)TIs an upper bound acoustic vector, wherein
Figure BDA0002719530090000114
This embodiment is specific
Figure BDA0002719530090000115
ρ0Representing the density of the medium above the sea surface, which can be seen as a vacuum, i.e. ρ, under absolutely soft (pressure release) boundary conditions0＝0g/cm3；
Similarly, the term of the boundary acoustic energy in the calculation domain is zero according to the fact that the term is only possible to propagate downwards and upwards, namely
Figure BDA0002719530090000116
Using the subscript "N" to denote the lower boundary, the joint vector at the lower boundary is found to be:
Figure BDA0002719530090000117
wherein, wN(kr,zN) Is the vertical vibration velocity kernel function of the lower boundary, phiN(kr,zN) Is the sound pressure kernel of the lower boundary, (1, B)N)TIn order to be the lower boundary acoustic vector,
Figure BDA0002719530090000118
this embodiment is specific
Figure BDA0002719530090000119
ρN+1Respectively representing the density of the medium under the sea bottom, which can be regarded as a rock with very high density under the condition of an absolute hard boundary, namely taking rhoN+1＝1099g/cm3；
S222, acoustic vector is divided from lower boundary zNTo the depth z of the sound sourcesPassing, i.e. passing the formula from bottom to top using vectors, from the lower boundary zNStarting, transmitting and calculating acoustic vector layer by layer until the depth z of the acoustic sourcesZ is specificallyN＝100m、zs25m, where the joint vector calculation immediately below the sound source depth is:
Figure BDA00027195300900001110
in the formula, wN(kr,zN) For the undetermined lower boundary vertical vibration velocity kernel function, the sound vector immediately below the sound source depth
Figure BDA00027195300900001111
S223, enabling the sound vector to be in the upper boundary z0To the depth z of the sound sourcesDelivery of, in particular z00m, where the joint vector calculation passed immediately above the source depth is:
Figure BDA0002719530090000121
in the formula, w0(kr,z0) For the kernel function of the vertical vibration velocity of the upper boundary to be determined, the sound vector immediately above the sound source depth
Figure BDA0002719530090000122
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:
Figure BDA0002719530090000123
i.e. vs+1(kr,zs)-vs(kr,zs)＝Δv(kr,zs) And unfolding to obtain:
Figure BDA0002719530090000124
then solving to obtain the kernel functions of the vertical vibration velocity of the lower and upper boundaries as follows:
Figure BDA0002719530090000125
Figure BDA0002719530090000126
s225, according to the calculated lower and upper boundary vertical vibration velocity kernel function wN+1(kr,zN) And w0(kr,z0) Calculating sound pressure kernel function phi (k) of each dielectric layer interface positionrZ) vertical vibration velocity kernel function w (k)rZ) value.
Through the steps, firstly, the acoustic vectors of all layers are transferred and solved layer by layer from the acoustic vector of the lower boundary upwards by using a transfer function matrix method until the acoustic vector of the position immediately below the sound source interface is solved; then, transmitting and solving each layer of sound vector layer by layer from the upper boundary sound vector downwards until solving the sound vector immediately above the sound source interface; at a sound source interface, establishing a joint vector equation of a lower part and an upper part by utilizing sound source conditions and solving the vertical vibration speed of a lower boundary and an upper boundary; and finally, multiplying the vertical vibration velocity of the lower boundary and the upper boundary by the sound vector stored in each layer of interface to obtain a sound pressure kernel function and a vertical vibration velocity kernel function.
In this embodiment, the sound pressure kernel function solved by each layer interface is firstly utilized to solve the sound pressure according to the Hankel inverse transformation integral formula; then, solving a vertical vibration velocity kernel function solved by each layer interface according to a vertical derivative of a Hankel inverse transformation integral formula; and then, the horizontal vibration velocity is solved according to the horizontal derivative of the Hankel inverse transformation integral expression by converting the derivative of the zero-order Bessel function into a first-order Bessel function, and the details are shown as follows.
In this embodiment, when the sound pressure is solved in step S3, the specific Hankel inverse transformation integral equation is as follows:
Figure BDA0002719530090000131
inverse Hankel transform horizontal wave number k in integral expressionrThe sound pressure dispersion calculation formula obtained after dispersion is as follows:
Figure BDA0002719530090000132
wherein, Δ krStep size of horizontal wave number and Δ kr＝2π/(rmaxnw)，rmaxFor maximum horizontal distance of sound field (r is taken in this embodiment)max＝3000m)，nwMinimum number of sampling points in a 2 pi oscillation period for Bessel function (n is taken in this embodiment)w＝10)，kr,nIs a discrete horizontal wavenumber and kr,n＝nΔkr-iεkI is an imaginary unit, εkFor complex offsets for preventing singularities, epsilon, in the process of solving the depth equationk＝3Δkr/(2πlog10e) M is the maximum index number of the discrete horizontal wavenumber and M ═ kmax/Δkr，kmaxTo cut off the wavenumber, this example takes k specificallymax＝20k0Wherein the wave number of seawater k02 pi f/1500. And 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 in this embodiment, the specific 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:
Figure BDA0002719530090000133
wherein P (r, z) is relative sound pressure in frequency domain, phi (k)rZ) is the sound pressure kernel, krIs the horizontal wavenumber, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, w (k)rZ) is the kernel function of the vertical vibration velocity, J0The 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),Vz(r, z)) of the vertical vibration velocity component Vz(r, z) is obtained by performing integral calculation according to a first transformation expression:
Figure BDA0002719530090000134
and the horizontal wave number k is calculated by the same method as the above-mentioned sound pressure calculationrAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the vertical vibration velocity as follows:
Figure BDA0002719530090000135
and S313, performing horizontal wave number integral calculation on each position (r, z) of the sound field by adopting a discrete integral calculation formula of the vertical vibration speed to obtain the vertical vibration speed corresponding to each position (r, z).
In the present embodiment, when calculating the vibration velocity vector in step S3, the step of solving the horizontal vibration velocity includes:
s321, obtaining a derivative of the 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
Figure BDA0002719530090000141
The second transformation is obtained as:
Figure BDA0002719530090000142
wherein P (r, z) is relative sound pressure in frequency domain, phi (k)rZ) is the sound pressure kernel, krIs the horizontal wavenumber, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, J0Is a zero order Bessel function, J1Is a first order Bessel function;
s322, carrying out integral calculation according to a second transformation expression to obtain a horizontal vibration velocity component V in the vibration velocity vectorr(r, z) is:
Figure BDA0002719530090000143
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 the horizontal wave number k is calculated by the same method as the above-mentioned sound pressure calculationrAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the horizontal vibration velocity as follows:
Figure BDA0002719530090000144
and S323, performing horizontal wave number integral calculation on each position (r, z) of the sound field by adopting a discrete integral calculation formula of the horizontal vibration speed to obtain the horizontal vibration speed corresponding to each position (r, z).
Through the steps, the integral expression of the vertical vibration speed and the horizontal vibration speed is obtained by carrying out derivation on the Hankel inverse transformation integral expression in the vertical direction and the horizontal direction, the horizontal wave number is dispersed by adopting the same method as the sound pressure integral expression, the vertical vibration speed and the horizontal vibration speed can be solved, and the vibration speed vector V (r, z) ═ V (V, z) is realizedr(r,z),Vz(r, z)).
In step S4 of this embodiment, the method further includes the step S3 of calculating propagation loss from sound pressure to form a sound field propagation loss scalar cloud, where the calculation formula of calculating propagation loss (scalar) from sound pressure is specifically:
TL(r,z)＝-20log10|P(r,z)| (30)
and according to the solved sound pressure P (r, z) and vibration velocity vector V (r, z) ═ Vr(r,z),Vz(r, z)) calculating a time-averaged sound intensity vector of each position (r, z) of the sound field according to equation (31);
Figure BDA0002719530090000145
wherein I (r, z) is a time-averaged sound intensity vector, i.e. the time average value of the sound energy passing through the unit area on a plane perpendicular to the vibration velocity direction,
Figure BDA0002719530090000146
respectively, the complex sound pressure P (r, z) and the conjugate complex of the vibration velocity V (r, z).
The scalar cloud chart of the propagation loss of the sound field obtained in the specific application embodiment is shown in fig. 3, and the streamline of the sound field sound intensity vector obtained by calculation (background color is the propagation loss) is shown in fig. 4, as can be seen from fig. 4, the regional sound intensity vector lines with small propagation loss (high energy) are relatively dense, and the regional sound intensity vector lines with large propagation loss (low energy) are relatively sparse.
As shown in fig. 2, when the method of this embodiment is applied in a specific application embodiment, a high performance computer workstation is combined, a program module capable of implementing the function of the method for predicting an ocean vector sound field of this embodiment is loaded in a data storage medium of the high performance computer workstation, the high performance computer workstation receives field measurement data including ocean depth, sound velocity, density, and the like, and sound source parameter information including sound source frequency, position, and the like, and after the method for predicting an ocean vector sound field, an ocean vector sound field prediction graphic image is generated, and further, the marine vector sound field prediction graphic image can be provided for subsequent vector hydrophone sound signal processing and analysis.
The embodiment also provides a system for forecasting a marine vector sound field, which comprises a processor and a memory, wherein the memory is used for storing marine environment data, sound source parameters and computer programs, the processor is used for executing the computer programs, and the processor is used for executing the computer programs so as to execute the method for forecasting the marine 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 processor and the memory.
In addition, the present embodiment also provides a computer readable storage medium having stored thereon marine environmental data, sound source parameters, and a computer program programmed or configured to perform the above-described method of enhancing a forecasted marine vector sound field.
As will be appreciated by one skilled in the art, embodiments of the present application may be provided as a method, system, or 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, and the like) having computer-usable program code embodied therein. The present application is directed to methods, apparatus (systems), and computer program products according to embodiments of the application wherein instructions, which execute via a flowchart and/or a processor of the computer program product, create means for implementing functions specified in the flowchart and/or block diagram block or blocks. These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks. These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
The foregoing is illustrative of the preferred embodiments of the present invention and is not to be construed as limiting the invention in any way. Although the present invention has been described with reference to the preferred embodiments, it is not intended to be limited thereto. Therefore, any simple modification, equivalent change and modification made to the above embodiments according to the technical spirit of the present invention should fall within the protection scope of the technical scheme of the present invention, unless the technical spirit of the present invention departs from the content of the technical scheme 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:
Figure FDA0003616838800000011
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 densitysIs 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.:
Figure FDA0003616838800000012
wherein phi (k)rZ) is the sound pressure kernel, krIs a horizontal wavenumber, J0Is a Bessel function;
integrating the two sides of the cylindrical coordinate underwater sound Helmholtz equation simultaneously
Figure FDA0003616838800000013
The depth equation can be found as:
Figure FDA0003616838800000014
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:
Figure FDA0003616838800000021
wherein phi (k)rZ) is the sound pressure kernel, kzIs a vertical wavenumber and
Figure FDA0003616838800000022
A+(kr) Representing a downwardly propagating item, A-(kr) 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:
Figure FDA0003616838800000023
wherein Γ ═ kz(ρ ω), ρ 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:
Figure FDA0003616838800000024
s212, obtaining an upper interface z of the mth layer according to the combined vector formed in the step S211m-1The joint vector of (a) is:
Figure FDA0003616838800000025
and the lower interface z of the mth layermThe joint vector of (a) is:
Figure FDA0003616838800000026
s213, combining and eliminating the upper and lower interface joint vector expressions of the mth layer obtained in the step S212
Figure FDA0003616838800000027
And then, obtaining a transfer formula of the joint vector from bottom to top as follows:
vm(kr,zm-1)＝Mm(kr)vm(kr,zm)
wherein M ism(kr) Is a transfer matrix of the mth layer medium from bottom to top, if the thickness of the mth layer is hm＝zm-zm-1Then M ism(kr) The expression of (a) is:
Figure FDA0003616838800000028
and the transfer formula and the transfer matrix from top to bottom are:
Figure FDA0003616838800000029
Figure FDA00036168388000000210
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
Figure FDA00036168388000000211
Using subscript "0" to denote the upper boundary, the sound pressure kernel at the upper boundary is obtained0(kr,z0) Hang downKernel function w of direct vibration velocity0(kr,z0) Respectively as follows:
Figure FDA0003616838800000031
Figure FDA0003616838800000032
and the joint vector at the upper boundary is:
Figure FDA0003616838800000033
wherein, w0(kr,z0) Vertical vibration velocity kernel function of upper boundary, vector (1, B)0)TIs an upper bound acoustic vector, wherein
Figure FDA0003616838800000034
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
Figure FDA0003616838800000035
Using the subscript "N" to denote the lower boundary, the joint vector at the lower boundary is found to be:
Figure FDA0003616838800000036
wherein, wN(kr,zN) Is the vertical vibration velocity kernel function of the lower boundary, phiN(kr,zN) Is the sound pressure kernel function of the lower boundary (1, B)N)TIs a lower boundary acoustic vector, wherein
Figure FDA0003616838800000037
S222, acoustic vector is divided from lower boundary zNTo the depth z of the sound sourcesPassing, i.e. from the lower boundary zNStarting, transmitting and calculating acoustic vector layer by layer until the depth z of the acoustic sourcesWherein the joint vector calculation formula immediately below the sound source depth is:
Figure FDA0003616838800000038
in the formula, wN(kr,zN) For the undetermined lower boundary vertical vibration velocity kernel function, the sound vector immediately below the sound source depth
Figure FDA0003616838800000039
S223, enabling the sound vector to be in the upper boundary z0To the depth z of the sound sourcesAnd (3) transferring, wherein a joint vector calculation formula obtained by transferring immediately above the sound source depth is as follows:
Figure FDA00036168388000000310
in the formula, w0(kr,z0) 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
Figure FDA0003616838800000041
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:
Figure FDA0003616838800000042
i.e. vs+1(kr,zs)-vs(kr,zs)＝Δv(kr,zs) And unfolding to obtain:
Figure FDA0003616838800000043
then solving to obtain the kernel functions of the vertical vibration velocity of the lower and upper boundaries as follows:
Figure FDA0003616838800000044
Figure FDA0003616838800000045
s225, according to the calculated lower and upper boundary vertical vibration velocity kernel function wN+1(kr,zN)、w0(kr,z0) Calculating sound pressure kernel function phi (k) of each dielectric layer interface positionrZ) 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:
Figure FDA0003616838800000046
converting the horizontal wave number k in the integral formula of the Hankel inverse transformationrThe sound pressure dispersion calculation formula obtained after dispersion is as follows:
Figure FDA0003616838800000047
wherein, Δ krIs a horizontal wave number step size and Δ kr＝2π/(rmaxnw)，rmaxMaximum horizontal distance of sound field, nwFor the minimum number of sampling points, k, of the Bessel function in a 2 pi oscillation periodr,nIs a discrete horizontal wavenumber and kr,n＝nΔkr-iεkI is an imaginary unit, εkIs a complex offset ofk＝3Δkr/(2πlog10e) M is the maximum index number of the discrete horizontal wavenumber and M ═ kmax/Δkr，kmaxIs 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:
Figure FDA0003616838800000048
wherein P (r, z) is relative sound pressure in frequency domain, phi (k)rZ) is the sound pressure kernel, krIs 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, J0The 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)，Vz(r, z)) of the vertical vibration velocity component Vz(r, z) performing integral calculation according to the first transformation formula to obtain:
Figure FDA0003616838800000051
and for horizontal wave number krAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the vertical vibration velocity as follows:
Figure FDA0003616838800000052
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
Figure FDA0003616838800000053
The second transformation is obtained as:
Figure FDA0003616838800000054
wherein P (r, z) is relative sound pressure in frequency domain, phi (k)rZ) is the sound pressure kernel, krIs the horizontal wave number, r is the coordinate in the horizontal direction, z is the coordinate in the vertical or depth direction, J0Is a Bessel function of the zero order,
Figure FDA0003616838800000055
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 vectorr(r, z) is:
Figure FDA0003616838800000056
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 krAnd (3) carrying out dispersion to obtain a dispersion integral calculation formula of the horizontal vibration velocity as follows:
Figure FDA0003616838800000057
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)＝-20log10|P(r,z)|
and according to the solved sound pressure P (r, z) and vibration velocity vector V (r, z) ═ Vr(r,z),Vz(r, z)) calculating a time-averaged intensity vector of each position (r, z) of the sound field according to the following formula;
Figure FDA0003616838800000058
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,
Figure FDA0003616838800000059
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.
CN202011083604.9A 2020-10-12 2020-10-12 Method, system and medium for forecasting ocean vector sound field Active CN112254798B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011083604.9A CN112254798B (en) 2020-10-12 2020-10-12 Method, system and medium for forecasting ocean vector sound field

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011083604.9A CN112254798B (en) 2020-10-12 2020-10-12 Method, system and medium for forecasting ocean vector sound field

Publications (2)

Publication Number Publication Date
CN112254798A CN112254798A (en) 2021-01-22
CN112254798B true CN112254798B (en) 2022-07-12

Family

ID=74241902

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011083604.9A Active CN112254798B (en) 2020-10-12 2020-10-12 Method, system and medium for forecasting ocean vector sound field

Country Status (1)

Country Link
CN (1) CN112254798B (en)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113108897B (en) * 2021-04-23 2022-09-06 自然资源部第三海洋研究所 Ocean environment noise field forecasting method based on non-uniform air seal sound source
CN113641954B (en) * 2021-07-20 2022-04-05 中国科学院声学研究所 Method and system for rapidly forecasting three-dimensional sound field in complex marine environment
CN114510848B (en) * 2022-04-20 2022-08-02 自然资源部第一海洋研究所 Offshore wind farm underwater noise calculation method, software and measurement device
CN116341408A (en) * 2023-03-13 2023-06-27 中国人民解放军国防科技大学 Ocean mode data coordinate conversion method for distance-related underwater acoustic application
CN116068903B (en) * 2023-04-06 2023-06-20 中国人民解放军国防科技大学 Real-time optimization method, device and equipment for robustness performance of closed-loop system

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004212121A (en) * 2002-12-27 2004-07-29 Kobayashi Rigaku Kenkyusho Object sound detection method and device therefor
CN102141431A (en) * 2010-02-01 2011-08-03 鸿远亚太科技（北京）有限公司 Method for measuring and transforming sound field in double-layer medium space
CN102997988A (en) * 2012-11-16 2013-03-27 哈尔滨工程大学 Pool testing method of low-frequency acoustic directivity of large submerged buoy vector hydrophone
JP2017227489A (en) * 2016-06-21 2017-12-28 Ｎｅｃネットワーク・センサ株式会社 Test system, waveform simulator device, test method and program
CN107576388A (en) * 2017-08-22 2018-01-12 哈尔滨工程大学 Three-dimensional structure sound source radiation sound field forecasting procedure under a kind of shallow sea channel
CN109489796A (en) * 2018-09-01 2019-03-19 哈尔滨工程大学 A kind of underwater complex structural radiation noise source fixation and recognition based on unit radiation method and acoustic radiation forecasting procedure
CN109556701A (en) * 2018-11-01 2019-04-02 浙江海洋大学 A kind of shallow sea geoacoustic inversion method based on broadband vertical wave impedance
CN110750934A (en) * 2019-11-01 2020-02-04 哈尔滨工程大学 Deep sea elastic structure and environment coupling sound radiation forecasting method
CN111639429A (en) * 2020-05-29 2020-09-08 中国人民解放军国防科技大学 Underwater sound field numerical simulation method, system and medium based on Chebyshev polynomial spectrum

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7814936B2 (en) * 2005-11-16 2010-10-19 Fisher Controls International Llc Sound pressure level feedback control

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004212121A (en) * 2002-12-27 2004-07-29 Kobayashi Rigaku Kenkyusho Object sound detection method and device therefor
CN102141431A (en) * 2010-02-01 2011-08-03 鸿远亚太科技（北京）有限公司 Method for measuring and transforming sound field in double-layer medium space
CN102997988A (en) * 2012-11-16 2013-03-27 哈尔滨工程大学 Pool testing method of low-frequency acoustic directivity of large submerged buoy vector hydrophone
JP2017227489A (en) * 2016-06-21 2017-12-28 Ｎｅｃネットワーク・センサ株式会社 Test system, waveform simulator device, test method and program
CN107576388A (en) * 2017-08-22 2018-01-12 哈尔滨工程大学 Three-dimensional structure sound source radiation sound field forecasting procedure under a kind of shallow sea channel
CN109489796A (en) * 2018-09-01 2019-03-19 哈尔滨工程大学 A kind of underwater complex structural radiation noise source fixation and recognition based on unit radiation method and acoustic radiation forecasting procedure
CN109556701A (en) * 2018-11-01 2019-04-02 浙江海洋大学 A kind of shallow sea geoacoustic inversion method based on broadband vertical wave impedance
CN110750934A (en) * 2019-11-01 2020-02-04 哈尔滨工程大学 Deep sea elastic structure and environment coupling sound radiation forecasting method
CN111639429A (en) * 2020-05-29 2020-09-08 中国人民解放军国防科技大学 Underwater sound field numerical simulation method, system and medium based on Chebyshev polynomial spectrum

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
声场波数积分最大截止波数自动选取算法;刘巍 等;《国防科技大学学报》;20190831;第41卷(第4期);第177-180页 *

Also Published As

Publication number Publication date
CN112254798A (en) 2021-01-22

Similar Documents

Publication Publication Date Title
CN112254798B (en) Method, system and medium for forecasting ocean vector sound field
CN112254797B (en) Method, system and medium for improving ocean sound field forecasting precision
CN108089155B (en) Passive positioning method for single hydrophone sound source in deep sea environment
CN107340540B (en) Direction wave decomposition method, device and the computer storage medium of elastic wave field
CN112577592B (en) Finite space plane near-field acoustic holography measurement method based on space Fourier transform
EP3353577B1 (en) Determining node depth and water column transit velocity
Robinson et al. Prediction systems with data assimilation for coupled ocean science and ocean acoustics
Randeni P et al. Determining the horizontal and vertical water velocity components of a turbulent water column using the motion response of an autonomous underwater vehicle
Golubev et al. Taking into account fluid saturation of bottom sediments in marine seismic survey
CN110068309B (en) Offshore water depth measurement method and device based on novel high-order dispersion relation
Van Komen et al. A CNN for range and seabed estimation on normalized and extracted time-series impulses
CN103575927B (en) The estimation method of the water speed of sound node
Ji et al. Performance evaluation and analysis for dipole source localization with lateral line sensor arrays
Czapiewska et al. Analysis of Impulse Responses Measured in Motion in a Towing Tank
Bacopoulos et al. Physics-based sizing functions for ocean-to-estuary shallow-water models
Vedenev et al. Airborne and Underwater Noise Produced by a Hovercraft in the North Caspian Region: Pressure and Particle Motion Measurements
Vidar Grindheim et al. Comparison of two models for prediction of seismic streamer state using the ensemble Kalman filter
Margolina et al. BRS Sound Exposure Modeling Tool: a system for planning, visualization and analysis
KR101938724B1 (en) Seismic imaging of water column structure apparatus and method using frequency-domain reverse-time migration with analytic Green's function
Meyer et al. Remote sensing of Tyrrhenian shallow waters using the adjoint of a full-field acoustic propagation model
JP2005024354A (en) Parallel sound field calculator
Chin-Bing et al. Analysis of coupled oceanographic and acoustic soliton simulations in the Yellow Sea: a search for soliton-induced resonances
Liu et al. Analysis of the influence of double reflection on the bathymetric function restoration
Bender et al. Estimates of the water column sound-speed field in three dimensions using a linearized perturbative technique
Martins et al. Bayesian acoustic prediction assimilating oceanographic and acoustically inverted data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant