Direction finding method for multiple targets in shallow water using single vector hydrophone based on particle swarm optimization compensation

Abstract

The invention discloses a shallow sea multi-target direction finding method of a single vector hydrophone based on particle swarm optimization compensation, comprising: 1) extracting channel information of a single vector hydrophone; 2) searching for a point source signal as a cooperative sound source; 3) initialization; 4) updating the fitness of each particle, and updating the change law of the inverse incomplete gamma function; 5) calculating the individual optimal value of the particle and the global optimal value of the group; 6) updating the particle speed and position; 7) repeating 4)-6) until the conditions are met, and obtaining the optimal compensation angles in the horizontal and vertical dimensions, corresponding to the maximum compensation of the target source; 8) extracting the instantaneous phase of each frequency unit of the four channels, and performing instantaneous phase difference weighted processing; 9) obtaining the acoustic energy flow beam domain information; 10) distinguishing and extracting the target in the acoustic energy flow beam domain, and obtaining the target characteristic spectrum information; 11) counting the number of target appearances; 12) performing phase compensation on the signal; and 13) performing azimuth estimation on multiple targets.

Description

Technical Field

[0001] The invention belongs to the technical field of signal processing, and relates to a target direction finding method in a shallow sea channel, in particular to a single-vector hydroacoustic multi-target direction finding method based on particle swarm optimization compensation. Background Art

[0002] In recent years, azimuth estimation has been a hot topic in the field of underwater acoustic signal processing. In this field, vector hydrophone technology, as an important means, has been widely recognized and received more and more attention. Compared with traditional scalar hydrophones, vector hydrophones can accurately estimate the azimuth of underwater targets by picking up three orthogonal components of sound pressure and vibration velocity in the sound field, which provides significant spatial advantages for underwater small platforms equipped with acoustic sensors to detect targets. Existing direction-finding algorithms are mainly based on deep-sea free fields. In the research of related algorithms, scholars such as Hui Junying discussed the basic principle of joint processing of sound pressure and vibration velocity, that is, under far-field conditions, the sound pressure and vibration velocity responses of underwater point source targets are completely coherent, but for isotropic background noise, they are incoherent; Yao Zhixiang et al. proposed the complex sound intensity method and the improved bar graph method based on this method. These methods are widely used in deep-sea free-field models, but the above methods have the problems of low accuracy and large variance under shallow sea conditions, and the comprehensive performance cannot fully meet the requirements of engineering applications. According to the analysis of the shallow sea model, in the shallow sea channel, each velocity component in the sound field and the sound pressure show different degrees of lack of correlation. When the existing traversal search method is used for phase compensation, the time complexity is positively correlated with the required accuracy. For every increase in accuracy, the number of searches will increase to ten times the original. This leads to excessive time complexity and reduced efficiency of the algorithm under the demand for high-precision direction finding. In addition, in the shallow sea environment, due to the limitations of low signal-to-noise ratio and multi-target conditions, the estimation accuracy and stability of the pitch angle are further reduced. In view of the problems existing in the above-mentioned prior art, the present invention has developed a single-vector hydrophone particle swarm optimization compensation high-precision direction finding technology under shallow sea channel conditions. Summary of the invention

[0003] The purpose of the present invention is to provide a method for shallow sea multi-target direction finding using a single vector hydrophone based on particle swarm optimization compensation. The present invention needs to first search for a cooperative sound source, and then based on the cooperative sound source, use the vector hydrophone sound pressure and vibration velocity correlation coefficient as the objective function, and use an improved particle swarm algorithm (IPSO) to iteratively search for a three-channel vibration velocity phase compensation angle, in the hope of achieving the optimal correlation coefficient. The present invention can achieve adaptive optimization of the vibration velocity channel phase compensation angle, thereby maximizing the correlation between sound pressure and vibration velocity. Combined with the instantaneous phase difference weighted processing method to improve the signal-to-noise ratio of the line spectrum and narrowband channel, and the screening model of high-precision feature spectrum extraction and multi-feature joint tracking is used to eliminate wild values, which significantly improves the direction finding accuracy under multi-target and low signal-to-noise ratio conditions, improves the algorithm robustness, and ultimately achieves high-precision azimuth estimation of the horizontal and vertical dimensions of multiple targets in shallow waters.

[0004] To achieve the above object, the present invention adopts the following technical solutions:

[0005] A single vector hydrophone shallow sea multi-target direction finding method based on particle swarm optimization compensation specifically includes the following steps:

[0006] Step 1: Extract the sound pressure channel and three-dimensional vibration velocity channel information of the single vector hydrophone, and then perform discretization processing.

[0007] Step 2: Search for a point source signal as the collaborative sound source and obtain its sound pressure p and horizontal vibration velocity v x And the vertical vibration velocity v z .

[0008] Step 3: Initialize the particle swarm algorithm parameters, initialize the speed and position of each particle, where the position follows the Beta distribution.

[0009] Step 4: Calculate the fitness of each particle and obtain the current iteration step number s. Dynamically update the inertia factor according to the iteration step number, and update it in accordance with the change law of the inverse incomplete gamma function.

[0010] Step 5: Update the individual optimal value of each particle and the global optimal value of the entire group.

[0011] Step 6: Optimize and update the particle speed and position, judge the boundary conditions, and adjust the data.

[0012] Step 7: Repeat steps 4 to 6 until the end condition is met and the optimal compensation angle in the horizontal dimension is obtained. and the optimal compensation angle in the vertical dimension Maximum compensation of the sound pressure-velocity correlation of the target source.

[0013] Step 8: Perform FFT transformation on the discrete signals of the target four-channel sound pressure velocity and convert them into corresponding frequency domain signals, then extract the instantaneous phase of each frequency unit of the four channels and perform instantaneous phase difference weighting processing.

[0014] Step 9: Perform horizontal conjugate cross-spectrum on the target characteristic spectrum, calculate the azimuth information of each frequency point, perform clustering processing on the acoustic energy in each direction, and obtain the acoustic energy flow beam domain information.

[0015] Step 10: According to the multi-feature joint model of left and right width boundaries, slope limitation and peak height threshold, the acoustic energy flow beam domain is discriminated and the target is extracted to obtain the target characteristic spectrum information, and the information is recorded.

[0016] Step 11: Count the number of times the target appears in each orientation window and eliminate suspected targets.

[0017] Step 12: Calculate the horizontal and vertical dimensions of the processed signal. Phase compensation.

[0018] Step 13: Perform high-precision position estimation on the horizontal and vertical dimensions of multiple targets.

[0019] Preferably, the sound pressure velocity signal received by the single vector hydrophone in the shallow sea environment described in step 2 and step 8 is the superposition of various sound fields. The point source signal conforms to the normalized vector output model of the shallow sea channel, and the expression is:

[0020]

[0021] Among them, p, v x 、v y and v z They are sound pressure and three-channel vibration velocity signals, the horizontal azimuth angle between the sound source and the receiving point is θ, and the propagation path length of the i-th sound line is r i , the corresponding pitch angle is φ i , the corresponding attenuation factor is a i , i=1,2,···, j is the imaginary unit, e is the base of the natural logarithm, the wave number k=ω / c, ω is the angular velocity of the signal, the acoustic impedance ratio Z=ρc is the constant of the ratio of sound pressure to the vibration velocity mode, and ρ is the density of the propagation medium.

[0022] Preferably, in step 3, the initial distribution of IPSO, i.e., the distribution of Beta distribution is as follows:

[0023]

[0024] Where x is the normalized particle distribution position and the denominator is the β function, which is defined as follows:

[0025]

[0026] Wherein, t is the time scale of the function, a and b are two parameters for controlling the probability density, and a=b<1 is adjusted. In this case, the shape of the density function is a symmetrical U-shape. In the present invention, a=b=0.8, which achieves a better effect.

[0027] Preferably, in step 4, the dynamic change mode of the inertia factor follows the inverse incomplete gamma function, and the specific adjustment mode is as follows:

[0028]

[0029] In the formula, ω max and ω min are the maximum and minimum values ​​of the inertia factor, respectively, where the incomplete gamma function γ is expressed as:

[0030]

[0031] The inverse incomplete gamma function is the inverse function of this function, where α>0 is the resolution parameter, α=s / Maxs, Maxs is the limited number of iteration steps, and λ is one of the constants. When λ is a constant, the inverse incomplete gamma function conforms to the rule of close to linear decline in the early stage and close to exponential decline in the later stage, which conforms to the best overall optimization strategy.

[0032] Preferably, in step 6, the particle velocity is updated as follows:

[0033] v ij (s+1)=v ij (s)+c1r1(s)[p ij (s)-x ij (s)]+c2r2(s)[g j (s)-x ij (s)]

[0034] Among them, c1 and c2 are learning factors, also known as acceleration constants, r1 and r2 are uniform random numbers in the range of [0,1], which are used to increase the randomness of particle flight, and v ij is the velocity of the particle, i=1,2,···,N, N is the number of particles, j=1,2,···,D, D is the unknown vector dimension, v ij ∈[-v max ,v max ],v max A constant used to limit the particle velocity, p ij is the optimal position (individual extreme value) searched by the current particle, x ij ∈[x min ,x max ] is the current particle position, x min and x maxUsed to limit the position of the particle, g j is the global optimal value currently searched by the entire group.

[0035] The way to update the particle position is:

[0036] x ij (s+1)=x ij (s)+v ij (s+1)

[0037] The speed and position are judged and adjusted by the following method:

[0038]

[0039] Preferably, in step 8, the instantaneous phase extraction method taking the sound pressure channel as an example is as follows:

[0040]

[0041] Where n'=1,2,···,N', N' is the data length, p is the discrete time domain sound pressure signal, P(f m ) is the frequency domain signal of the sound pressure channel, m=1,2,···,M, where M is the number of all frequency units, To extract the instantaneous phase, the spectrum expression obtained after weighting the frequency domain signal is:

[0042]

[0043] in, is the weighted value obtained after the variance average fitting and mean filtering of the instantaneous phase difference values ​​of the four channels, and is the spectrum weighted by the instantaneous phase difference.

[0044] Preferably, in step 10, the left and right boundary widths are less than 10; the target azimuth beam slope exceeds 1; the amplitude exceeds 3 times the average acoustic energy flow size, and the acoustic energy flow beam domain information is automatically extracted. The instantaneous frequency corresponding to the target azimuth is obtained as f g , g=1,2,...,G, the corresponding frequency unit is g, and the corresponding instantaneous frequency under interference noise is f y , y=1,2,...,Y, the corresponding frequency unit is y.

[0045] At horizontal speed For example, set the initial value and subsequent value of each frequency window to:

[0046]

[0047] Where W is the window weight, n = n'Ts , T s is the sampling period;

[0048] To change the window weight, use the following formula:

[0049]

[0050] Among them, W t is the size of the window weight at time t, a1 and a2 are the weight changes of the target instantaneous frequency and the noise instantaneous frequency in each iteration, and ε1 is the target instantaneous frequency window width.

[0051] Preferably, in step 12, the method for compensating the horizontal and vertical channel phases is:

[0052]

[0053] in, and is the phase compensation angle in the horizontal and vertical dimensions, v' x (n), v' y (n) and v' z (n) is the out-of-phase wave sequence after compensation.

[0054] Preferably, the precise direction finding sound pressure velocity conjugate cross-spectrum method in step 13 is as follows:

[0055]

[0056] in, and V x (f) V y (f) and V z (f) Conjugate spectrum.

[0057] After averaging the cross-spectrum results over time and taking the real part, the acoustic energy flow I of each frequency energy is obtained. x (f) I y (f) and I z (f) Use the acoustic energy flow to estimate the horizontal azimuth and elevation angle at the frequency point:

[0058]

[0059] What you get and are the estimated values ​​of the horizontal azimuth and elevation angle respectively.

[0060] Divide [0°, 360°] into several intervals Δθ, and divide [-90°, 90°] into several intervals Δφ, then the total number of intervals allowed for the spatial angle is:

[0061]

[0062] Scan the entire frequency band to obtain the corresponding angle. Assume that the number of estimated azimuths of each frequency point in the horizontal dimension that fall into each statistical angle interval is m1, m2, m3, ..., The number of vertical dimensions are m1, m2, m3, ..., Then calculate the following statistic R(θ):

[0063]

[0064] The distribution of acoustic energy flow in each statistical interval is obtained; where R θ (f) and R φ (f) is the modulus of the acoustic energy flow that conforms to the angle range of the statistical interval, and the specific calculation formula is as follows:

[0065]

[0066] The modulus value is used as the weight of the corresponding azimuth count, and the estimated values ​​of each frequency point are counted. The final result is the accurate horizontal azimuth and pitch angle.

[0067] The invention can be used for direction finding in shallow sea environment, and is generally suitable for estimating the direction of targets with a sound source frequency of more than 200 Hz.

[0068] The present invention has the following beneficial effects:

[0069] 1. The present invention aims at the problem of missing correlation between sound pressure and velocity when finding direction under shallow sea channel conditions. By searching for cooperative sound sources and using particle swarm algorithm to perform adaptive search calculation, the phase compensation angle is determined to achieve the maximum correlation between sound pressure and velocity. The present invention can achieve the optimal compensation of the correlation between the sound pressure and velocity of the target sound source, and can greatly improve the azimuth estimation accuracy of the single vector hydrophone in shallow sea environment.

[0070] 2. The present invention adopts an improved particle swarm algorithm to realize adaptive compensation of the sound pressure-velocity correlation. Compared with the existing traversal search method, the present invention has lower time complexity and is suitable for the needs of high-precision direction finding. Compared with the standard particle swarm algorithm, the algorithm adopted by the present invention makes it easier for particles to surround the target in the setting of initial values, and the addition of dynamically changing inertia weights enables the algorithm to obtain excellent global search capabilities in the early stage, which is conducive to jumping out of the local optimum, and obtaining excellent local search capabilities in the later stage, which is conducive to improving accuracy.

[0071] 3. The present invention is combined with the instantaneous phase difference weighted processing method to improve the signal-to-noise ratio of the line spectrum and narrowband signal, and uses a high-precision feature spectrum extraction and multi-feature combined target screening model to eliminate outliers, which significantly improves the direction finding accuracy under multi-target and low signal-to-noise ratio conditions. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 It is a schematic diagram of the sound propagation model in shallow water channels.

[0073] Figure 2 This is the flow chart of particle swarm algorithm to obtain the optimal phase compensation angle.

[0074] Figure 3 It is a flow chart of the frequency locking and tracking method.

[0075] Figure 4 It is a flow chart of a shallow sea multi-target direction finding method of a single vector hydrophone based on particle swarm optimization compensation in a preferred embodiment of the present invention.

[0076] Figure 5 It is a schematic diagram of the experimental system deployment.

[0077] Figure 6 This is a convergence comparison chart of particle swarm algorithm, improved particle swarm algorithm and traversal search method.

[0078] Figure 7 It is a scatter plot of three-dimensional orientation estimation under different signal-to-noise ratios of the present invention.

[0079] Figure 8 It is a comparison diagram of the signal-to-noise ratio-pitch angle estimation error of the present invention and the conventional complex sound intensity method under shallow sea multi-target conditions.

[0080] FIG9(a) and FIG9(b) are horizontal dimension time-azimuth history diagrams of the conventional complex sound intensity method and the method of the present invention, respectively.

[0081] FIG. 10( a ) and FIG. 10( b ) are vertical dimension time-azimuth history diagrams of the conventional complex sound intensity method and the method of the present invention, respectively. DETAILED DESCRIPTION

[0082] The preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0083] The present invention is a shallow-water multi-target direction finding method using a single vector hydrophone based on particle swarm optimization compensation. The existing complex sound intensity azimuth estimation method is suitable for a deep-sea free-field model with basically consistent sound pressure velocity phase differences, but the existing method often has the problem of low direction finding accuracy in shallow-water channels, and the multi-target resolution is low in the vertical dimension pitch angle estimation. In view of the above problems, the present invention improves the correlation of the sound pressure velocity channel efficiently and stably by improving the particle swarm optimization algorithm; improves the target characteristics of multiple targets by referencing multi-feature line spectrum extraction, and proposes a shallow-water multi-target direction finding method using a single vector hydrophone based on particle swarm optimization compensation.

[0084] Firstly, the principle of single vector hydrophone direction finding in shallow sea channels is introduced.

[0085] 1.1 Shallow water channel acoustic vector model

[0086] The most significant difference between the shallow sea multipath model and the deep sea free field model is the consideration of the sea surface and seabed boundary conditions. This boundary condition leads to the multipath effect, that is, the diversity of sound wave propagation paths caused by the boundary. The sound field at the receiving point is no longer a single direct sound line, but a collection of multiple sound lines. The sound propagation model established according to the ray acoustics theory, such as Figure 1 As shown in the figure, each reflected sound line can be regarded as a direct sound line emitted by a virtual source. Based on this theory, a ray model is established. Suppose the sound source in a shallow sea multipath environment is located at point O, the receiving point is located at point R, the horizontal azimuth angle from the sound source to the receiving point is θ, and the corresponding pitch angle is φ i , the propagation path of the sound is r i In three-dimensional space, the description of the vibration velocity signal follows the Cartesian coordinate definition. According to acoustic theory, the normalized vector output model at the receiving point can be expressed as:

[0087]

[0088] Where p is the sound pressure signal picked up by the vector hydrophone, v x and v y are the horizontal X-axis and Y-axis vibration velocity signals picked up by the vector hydrophone, v z is the vertical Z-axis vibration velocity signal picked up by the vector hydrophone, where wave number k = ω / c, ω is the angular velocity of the signal, c is the underwater sound propagation speed, and a i is the normalized attenuation factor, which is caused by the combined effect of the sea surface and the seabed. The acoustic impedance ratio Z = ρc is a proportional constant of the ratio of the sound pressure to the vibration velocity mode, and this constant is only related to the properties of the medium itself.

[0089] Will Figure 1 The xΟy plane shown is defined as the sea surface, and it is assumed to be a flat free interface. The coordinates of the sound source O are set to (0, h), and the depth of the seabed is set to H. For the convenience of calculation, the sound line is classified as BS i ,BB i , SS i , S.B. i There are four types (i=1,2...), the first B and S respectively indicate that the first reflection is the sea surface reflection or the seabed reflection, the second B and S respectively indicate that the last reflection is the sea surface reflection or the seabed reflection, and each time a sea surface and seabed reflection is added, the corresponding number of sound line groups n increases by one.

[0090] Each sound ray corresponds to a pitch angle. According to the law of mirror reflection and the equivalent virtual source method, the pitch angle of each type of sound ray can be obtained as shown in formula (2):

[0091]

[0092] Where z is the depth of the hydrophone and x is the horizontal distance between the sound source and the hydrophone.

[0093] The path length of each type of sound line is shown in formula (3), where r D is the path length of the direct sound ray.

[0094]

[0095] The number of reflections of each type of sound line on the sea surface and the seabed is different, so the corresponding reflection coefficient of the sound line is also different. It is a multiple superposition of the sea surface and the seabed. The normalized attenuation coefficient can be calculated and summarized as:

[0096]

[0097] Where a is the seawater absorption coefficient, r s 、r b are the sea surface reflection coefficient and the seabed reflection coefficient respectively.

[0098] 1.2. Single vector hydrophone cross-spectrum direction finding method

[0099] In the prior art, the average sound intensity method is widely used to estimate the direction of a single sound source. This method can achieve the maximum likelihood estimation of the azimuth angle. However, if there are multiple sound sources, the average sound intensity method will fail because it can only provide the direction of the synthetic sound energy flow of multiple sound sources, but cannot distinguish them. In contrast, the complex sound intensity method based on cross-spectral direction finding is a method for estimating the direction of sound energy flow in the frequency domain. As long as the power spectra of multiple sound sources are not completely consistent, the complex sound intensity method can distinguish the directions of multiple targets. The specific method is as follows:

[0100] First, the sound pressure and three-channel vibration velocity signals are cross-spectral, and the results are:

[0101]

[0102] Among them, P(f), V x (f) V y (f) and V z (f) are p and v respectively x 、v y and v z The spectrum of the signal, and V x (f) V y (f) and V z (f) Conjugation spectrum.

[0103] Use a sliding window to perform time averaging on the cross-spectrum results, and the average periodogram is:

[0104]

[0105] According to the basic characteristics of Fourier transform, the energy of two signals with the same phase is mainly concentrated in the real part. In the shallow sea channel, the horizontal vibration velocity v x and v y The correlation loss with the sound pressure is small and is approximately in phase, and the vertical vibration velocity and sound pressure also basically conform to this characteristic after correlation compensation. Therefore, the sound energy flow I of each frequency energy is obtained by taking the real part of each cross-spectrum result. y (f) I y (f) and I z (f), the horizontal azimuth and elevation angles at this frequency are estimated using the acoustic energy flow as shown in formula (7):

[0106]

[0107] In this embodiment, a method for shallow-water multi-target direction finding using a single vector hydrophone based on particle swarm optimization compensation is specifically described as follows.

[0108] High-precision direction finding in shallow waters must first solve the problem of matching the sound pressure velocity channel with the shallow water channel. Through the analysis of sound pressure velocity under the shallow water multi-path model, it is found that there is a different degree of correlation loss between the velocity and the sound pressure in this environment, while the traditional direction finding method is generally based on the deep sea free field and the joint processing method of the sound pressure velocity is completely correlated. Therefore, in order to achieve high-precision direction finding, it is first necessary to adaptively compensate the phase of the velocity channel. The present invention introduces the concept of collaborative sound source and uses particle swarm algorithm optimization. Through this method, the maximum phase compensation angle is iteratively searched to maximize the correlation coefficient, and finally the optimal compensation of the correlation is achieved; then, the instantaneous phase difference weighted processing is combined to improve the signal-to-noise ratio of the narrowband signal, and the target line spectrum is locked and tracked by introducing a high-precision feature spectrum extraction and a multi-feature combined target positioning model, thereby improving the robustness of the algorithm under low signal-to-noise ratio and multi-target conditions, and improving the pitch angle estimation accuracy of the algorithm in this environment. The overall process of the method is as follows: Figure 4 As shown, the specific steps are as follows:

[0109] Step 1: The three-dimensional vector hydrophone receives the sound pressure and vibration velocity signals of the target source, first filters and amplifies the received signal, the frequency band of the bandpass filter is 200Hz~400Hz, and then samples it through the A / D circuit with a sampling rate of f s , and get the discrete signals p(n), v x (n), v y (n), v z (n), where n = n'T s , T s =1 / f s , n'=1,2,···,N', N' is the data length.

[0110] Step 2: Search for a target as a cooperative sound source. All non-primary targets with a certain correlation between sound pressure and vibration velocity can be used as cooperative sound sources. In the marine environment, isotropic noise is uncorrelated with sound pressure and vibration velocity, while anisotropic noise is correlated. Therefore, in the absence of an obvious target, anisotropic noise can also be used as a cooperative target. Then obtain the sound pressure p and horizontal vibration velocity v of the cooperative sound source. x and vertical velocity channel signal v z , where the horizontal two-dimensional vibration velocity v x 、v y The phase compensation angles are approximately equal. x Perform horizontal dimension compensation.

[0111] Step 3: Use the Pearson correlation coefficient of the sound pressure-vibration velocity signal as the objective function of the particle swarm algorithm search. The correlation coefficient calculation method is as follows:

[0112]

[0113] Where p and v' conform to the shallow sea vector output model, v' is the vibration velocity signal after phase compensation, and is the mean of the corresponding sequence. First, the Hilbert transform is used to resolve the original signal into a complex signal for subsequent processing:

[0114]

[0115] Perform phase offset compensation on the horizontal and vertical vibration velocities to obtain out-of-phase waves. The specific method is as follows:

[0116]

[0117] in, and is the phase offset angle of the horizontal dimension and the vertical dimension, and is also the search target of the particle swarm algorithm, v' x (n), v' y (n) and v' z (n) is the out-of-phase wave sequence after compensation.

[0118] First, the initial position and speed of each particle are set. In order to form an "encirclement" of the optimal solution, this embodiment proposes a population initialization method based on Beta distribution. The Beta distribution function is defined as:

[0119]

[0120] Where x is the normalized particle distribution position and the denominator is the β function, which is defined as follows:

[0121]

[0122] Among them, a and b are two parameters that control the probability density. When a=b<1, the shape of the density function is a symmetrical U-shape. At this time, the candidate solution is most likely to be located near the boundary of the search space, so the global optimal solution is well "surrounded" by the initial particle population. Experiments show that when a=b=0.8, the initialization effect is better.

[0123] Let the number of particles be N, the optimization space be D-dimensional, and the candidate solution after initialization be Xi. According to the Beta distribution, each dimension of Xi is assigned an initial value as follows.

[0124] X ij =X min +(X max -X min )×betarnd (13)

[0125] Among them, i=1,2,···,N, j=1,2,···,D, and betarnd is a random number that conforms to the Beta distribution.

[0126] Step 4: Calculate the fitness of the current particle, obtain the current particle swarm algorithm iteration step s, and iterate the inertia factor ω according to the inverse incomplete gamma function, where the incomplete gamma function is:

[0127]

[0128] The inverse incomplete gamma function is the inverse function of this function, where α>0 is the resolution parameter, α=s / Maxs, Maxs is the limited number of iteration steps. When λ takes a constant value, the inverse incomplete gamma function conforms to the dynamic descent law, making the inertia factor larger in the early stage, and the algorithm's global search ability stronger, and smaller in the later stage, and the algorithm's local optimization ability stronger.

[0129] The inertia factor is adjusted as follows:

[0130]

[0131] In the formula, ω max and ω min are the maximum and minimum values ​​of the inertia factor, respectively.

[0132] Step 5: Update the individual optimal value p of each particle best and the global optimal value g of the entire group best , their expressions are as follows:

[0133]

[0134] Among them, p best and gbest Variables under the conditions of D-dimensional optimization space and N particles.

[0135] Step 6: Optimize and update the particle's speed and position. The update formula is as follows:

[0136] v ij (s+1)=v ij (s)+c1r1(s)[p ij (s)-x ij (s)]+c2r2(s)[g j (s)-x ij (s)] (17)

[0137] x ij (s+1)=x ij (s)+v ij (s+1) (18)

[0138] Among them, c1 and c2 are learning factors, also known as acceleration constants, r1 and r2 are uniform random numbers in the range of [0,1], which are used to increase the randomness of particle flight, and v ij is the velocity of the particle, v ij ∈[-v max ,v max ],v max A constant used to limit the particle velocity, p ij is the optimal position (individual extreme value) searched by the current particle, x ij ∈[x min ,x max ] is the current particle position, x min and x max Used to limit the position of the particle, g j is the global optimal value currently searched by the entire group. Then judge the boundary conditions and adjust the data that crosses the boundary conditions:

[0139]

[0140] Step 7: Repeat steps 4 to 6 until the end condition is met and the optimal compensation angles in the two horizontal dimensions are obtained. and the optimal compensation angle in the vertical dimension The maximum compensation of the sound pressure-velocity correlation of the target source, where the maximum number of iterations Maxs, the calculation accuracy or the maximum number of stagnation steps Δs of the optimal solution are all end conditions. The process involved in step 3 to step 7 is as follows Figure 2 shown.

[0141] Step 8: Perform FFT transformation on the four channels, convert the time domain to the frequency domain, and obtain the frequency domain signal of each channel as P(f m )、Vx(fm )、Vy(f m ) and Vz(f m ), m = 1, 2, ···, M, where M is the number of all frequency units, and the instantaneous phase of all frequency units is Extraction, taking the sound pressure channel as an example, the extraction method is as follows:

[0142]

[0143] After extracting the instantaneous phase difference, the variance of the difference is calculated, the variance of the four channels is averaged and fitted, and the weighted value is obtained after mean filtering. The signal spectrum is weighted by phase difference, and the spectrum expression obtained is:

[0144]

[0145] Step 9: Perform horizontal conjugate cross-spectrum on the target characteristic spectrum according to equations (5) to (7), calculate the azimuth information of each frequency point, perform clustering on the acoustic energy in each azimuth, and obtain the acoustic energy flow beam domain information.

[0146] Step 10: According to the multi-feature joint model of left and right width boundaries, slope limitation and peak height threshold, the acoustic energy flow beam domain is discriminated and the target is extracted to obtain the target characteristic spectrum information, and the information is recorded;

[0147] Step 11: Count the number of times the target appears in each orientation window and remove suspected targets;

[0148] The overall process of step 10 and step 11 is as follows Figure 3 As shown in the figure, the left and right boundary widths should be less than 10; the target azimuth beam slope should be greater than 1; the amplitude should be greater than 3 times the average acoustic energy flow size, and the acoustic energy flow beam domain information is automatically extracted. The instantaneous frequency corresponding to the target azimuth is f g , g=1,2,...,G, the corresponding frequency unit is g, and the corresponding instantaneous frequency under interference noise is f y , y=1,2,...,Y, the corresponding frequency unit is y.

[0149] by For example, set the initial value and change value of each frequency window as:

[0150]

[0151]

[0152] Update the sliding window weight according to formula (24), where a1 and a2 are the weight changes of the target instantaneous frequency and the noise instantaneous frequency in each iteration, and ε1 is the target instantaneous frequency window width:

[0153]

[0154] Set the bandwidth threshold to 10, and treat the next frame and the previous frame as the same target if the azimuth difference is within the bandwidth threshold. If not, continue to record the previous frame for comparison with subsequent data. Delete suspected targets whose target beam domain extraction times are less than half of the frames.

[0155] Step 12: Calculate the horizontal and vertical dimensions of the processed signal. The phase compensation is shown in formula (10).

[0156] Step 13. The azimuth and elevation angles corresponding to each frequency point can be counted. When the frequency range of the target radiation noise is unknown, the target azimuth cannot be directly obtained. Therefore, it is necessary to use the histogram statistics method to count the azimuths corresponding to each frequency point. When there is a broadband target signal, the azimuths corresponding to each frequency point of the signal are within a certain range, while the azimuths corresponding to each frequency point in the noise band are randomly distributed. Divide [0°, 360°] into several intervals Δθ, and divide [-90°, 90°] into several intervals Δφ, then the total number of intervals allowed for the spatial angle is:

[0157]

[0158] Scan the entire frequency band to obtain the corresponding angle. Assume that the number of estimated azimuths of each frequency point in the horizontal dimension that fall into each statistical angle interval is m1, m2, m3, ..., The number of vertical dimensions are m1, m2, m3, ..., Then calculate the following statistic R(θ):

[0159]

[0160] The distribution of acoustic energy flow in each statistical interval is obtained by formula (27). θ (f) and R φ (f) is the modulus of the acoustic energy flow that conforms to the angle range of the statistical interval, and the specific calculation formula is as follows:

[0161]

[0162] These modulus values ​​are used as weights of the corresponding azimuth counts, and the estimated values ​​of each frequency point are counted. The final result is the accurate horizontal azimuth and pitch angle.

[0163] Simulation and performance analysis:

[0164] In order to test the performance of the improved particle swarm algorithm of the present invention, a simulation experiment of optimal phase compensation angle search was first carried out. The parameters in the experiment were set as follows: the population size N of IPSO and PSO was 50, the maximum number of iterations Maxs was set to 100, and the particle value boundary x min and x max -90 and 90 respectively, the maximum particle speed v max is 10, where IPSO uses Beta distribution to set the initial value and follows the descending law of the inverse incomplete gamma function to realize the dynamic change of the inertia factor. The traversal search method (search) searches for targets within the range of [-90°, 90°] with a search accuracy of 1°. The search target is the vertical vibration velocity v in the shallow sea environment. z The optimal phase compensation angle, the loss value is the absolute value of the difference between the corresponding target correlation coefficient and the optimal correlation coefficient that can be achieved, given by Figure 6 As shown in the figure, compared with PSO, IPSO converges faster in the early stage, and in the later stage, due to the smaller inertia factor, its local optimization ability is better than PSO, which makes the results obtained by IPSO have smaller errors. When the error value is 0.1, the number of iterations of the traversal search method is 1800 times, while the IPSO obtained from multiple tests only needs 41.84 times on average, which is 912.4% faster than the search method in terms of time.

[0165] In order to test the robustness of the method at various azimuth angles, a total of 9 points were arranged at 30, 60 and 90 degrees horizontal azimuth, -20, 0 and 20 degrees elevation, and 100 separate direction finding were performed at different signal-to-noise ratios. Their azimuth estimates were recorded and scatter plots were drawn. The results are shown in Figure 2. Figure 7 As shown, the present invention has a large position estimation variance under low signal-to-noise ratio, but can perform high-precision position estimation on targets with different three-dimensional spatial distributions under high signal-to-noise ratio.

[0166] Monte Carlo comparison experiments were conducted between the method of the present invention and the traditional complex sound intensity method in the signal-to-noise ratio range of -20 to 30 dB, and 40,000 independent direction finding experiments were conducted for each spectral level. The sampling length was T = 10 s, and the sampling frequency was f s is 1024Hz. In each experiment, two target sound sources are randomly set to simulate multi-target direction finding. The frequency of each target line spectrum is uniformly distributed within 200Hz~400Hz, the horizontal azimuth angle is uniformly distributed within 0°~360°, the sound source depth h, the hydrophone depth and the horizontal distance x between the two are randomly distributed. Under this condition, the pitch angle is uniformly distributed within -90°~90°, so that the target signal is in the ocean environment noise simulated by Gaussian white noise. The simulation results are shown in Figure 2. Figure 8As shown, a relationship diagram between the root mean square error (RMSE) of azimuth estimation and the signal-to-noise ratio is obtained. Under different signal-to-noise ratios, the method of the present invention has a lower root mean square error, which is 34.05% lower than the root mean square error of the traditional complex sound intensity method on average. Experiments show that the method of the present invention has higher direction finding accuracy in shallow sea environments.

[0167] The horizontal dimension time-azimuth history diagram of the traditional complex sound intensity method and the method of the present invention under shallow sea multi-target conditions is shown in Figure 9, where Figure 9 (a) is the result of the traditional complex sound intensity method, and Figure 9 (b) is the result of the method of the present invention. In the simulation, the sea surface is an ideal mirror, that is, r s =1, the seafloor reflection coefficient is r b =0.7, the speed of sound is c=1500m / s, the seabed depth H is 60m, the depth of sound source 1 is h1=45m, the depth of sound source 2 is h2=10m, the depth of the receiving point hydrophone is z=46m, the horizontal distance between the target sound source and the hydrophone is x=100m, the target signal is a broadband signal superimposed on the main frequency signal of 200-400Hz and the ocean background noise simulated by Gaussian white noise, the horizontal azimuth angles of target 1 and target 2 are 60° and 170° respectively, the pitch angles of target 1 and target 2 are 0° and 20° respectively, the sampling frequency is 1024Hz, the signal-to-noise ratio is 20dB, and the number n of each type of sound line is 10. As shown in the figure, the time-azimuth history diagram of the method of the present invention is clearer than that of the complex sound intensity method; the vertical dimension time-azimuth history diagrams of the two methods are shown in Figure 10, wherein Figure 10(a) is the result of the traditional complex sound intensity method, and Figure 10(b) is the result of the method of the present invention. The complex sound intensity method can hardly distinguish the pitch angles of multiple targets in a shallow sea environment. In comparison, the present method has good resolution.

[0168] In summary, the present invention discloses a shallow-water multi-target direction finding method using a single vector hydrophone based on particle swarm optimization compensation. The traditional vector hydrophone direction finding method is based on a deep-sea free-field model, and has the problems of low direction finding accuracy and large variance in shallow-water environment. The present invention first proposes an efficient search method based on IPSO to obtain the optimal phase compensation angle. After the compensation phase angle is obtained with the assistance of a collaborative sound source, the maximum correlation compensation of the three-channel vibration velocity of the target sound source can be achieved. In addition, the present invention improves the signal-to-noise ratio of narrowband signals through instantaneous phase weighted processing, and improves the robustness of the algorithm under multi-target conditions through high-precision line spectrum extraction and multi-feature tracking models. The simulation analysis results prove the effectiveness of the method of the present invention and its good engineering practicality.

[0169] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any technician familiar with the technical field can easily think of changes or substitutions within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the scope of the present application should be based on the protection scope of the claims.

Claims

1. A single vector hydrophone shallow sea multi-target direction finding method based on particle swarm optimization compensation, characterized by The following steps are involved: Step 1, extracting the sound pressure channel and three-dimensional vibration velocity channel information of the single vector hydrophone and performing discretization processing; Step 2: Search for a point source signal as a collaborative sound source and obtain the sound pressure p and horizontal vibration velocity v x And the vertical vibration velocity v z ; Step 3: Initialize the particle swarm algorithm parameters, initialize the speed and position of each particle, where the position follows the Beta distribution; Step 4: Update the fitness of each particle and obtain the current iteration step number s. Dynamically update the inertia factor according to the iteration step number, and update in accordance with the change law of the inverse incomplete gamma function. Step 5: Calculate the individual optimal value of each particle and the global optimal value of the entire group; Step 6: Optimize and update the speed and position of the particles, judge the boundary conditions, and trim the data; Step 7: Repeat steps 4 to 6 until the end condition is met to obtain the optimal compensation angle in the horizontal dimension. and the optimal compensation angle in the vertical dimension Maximum compensation of the sound pressure-velocity correlation of the target source; Step 8: Perform FFT transformation on the discrete signals of the target four-channel sound pressure velocity to convert them into corresponding frequency domain signals, then extract the instantaneous phase of each frequency unit of the four channels, and perform instantaneous phase difference weighting processing; Step 9, perform horizontal conjugate cross-spectrum on the target characteristic spectrum, calculate the azimuth information of each frequency point, perform clustering processing on the acoustic energy in each direction, and obtain the acoustic energy flow beam domain information; Step 10: According to the multi-feature joint model of left and right width boundaries, slope limitation and peak height threshold, the acoustic energy flow beam domain is discriminated and the target is extracted to obtain the target characteristic spectrum information, and the information is recorded; Step 11: Count the number of times the target appears in each orientation window and remove suspected targets; Step 12: Calculate the horizontal and vertical dimensions of the processed signal. Phase compensation; Step 13: Estimate the horizontal and vertical directions of multiple targets.

2. The single vector hydrophone shallow sea multi-target direction finding method based on particle swarm optimization compensation as claimed in claim 1, characterized in that: In step 2 and step 8, the sound pressure velocity signal received by the single vector hydrophone in the shallow sea environment is the superposition of various sound fields, and the point source signal conforms to the normalized vector output model of the shallow sea channel, and the expression is as follows: Among them, p, v x 、v y and v z They are sound pressure and three-channel vibration velocity signals, the horizontal azimuth angle between the sound source and the receiving point is θ, and the propagation path length of the i-th sound line is r i , the corresponding pitch angle is ϕ i , the corresponding attenuation factor is a i , i=1,2,···, j is the imaginary unit, e is the base of the natural logarithm, the wave number k=ω / c, ω is the angular velocity of the signal, c is the underwater sound propagation speed, the acoustic impedance Z=ρc is the constant of the ratio of sound pressure to the vibration velocity mode, and ρ is the density of the propagation medium.

3. The single vector hydrophone shallow sea multi-target direction finding method based on particle swarm optimization compensation as claimed in claim 1, characterized in that: In step 3, the distribution of Beta distribution is as follows: Where x is the normalized particle distribution position and the denominator is the β function, which is defined as follows: Among them, t is the time scale of the function, a and b are two parameters that control the probability density, and a=b<1 is adjusted.

4. The method for shallow sea multi-target direction finding using single vector hydrophone based on particle swarm optimization compensation according to any one of claims 1 to 3, characterized in that: In step 4, the dynamic change of the inertia factor follows the inverse incomplete gamma function, and the specific adjustment method is as follows: In the formula, ω max and ω min are the maximum and minimum values ​​of the inertia factor, respectively, where the incomplete gamma function γ is expressed as: The inverse incomplete gamma function is the inverse function of the above function, where α>0 is the resolution parameter, α=s / Maxs, Maxs is the limited number of iteration steps, and λ is a constant. When λ is a constant, the inverse incomplete gamma function conforms to the rule of close to linear decline in the early stage and close to exponential decline in the later stage.

5. The method for shallow sea multi-target direction finding using single vector hydrophone based on particle swarm optimization compensation as claimed in claim 4, characterized in that: In step 6, the particle velocity is updated as follows: v ij (s+1)=v ij (s)+c1r1(s)[p ij (s)-x ij (s)]+c2r2(s)[g j (s)-x ij (s)] Among them, c1 and c2 are learning factors, r1 and r2 are uniform random numbers in the range of [0,1], which are used to increase the randomness of particle flight, and v ij is the velocity of the particle, i=1,2,···,N, N is the number of particles, j=1,2,···,D, D is the unknown vector dimension, v ij ∈[-v max ,v max ],v max A constant used to limit the particle velocity, p ij is the optimal position searched by the current particle, x ij ∈[x min ,x max ] is the current particle position, x min and x max Used to define the position of particles, g j The global optimal value currently searched by the entire group; The particle position is updated as follows: x ij (s+1)=x ij (s)+v ij (s+1) The speed and position are judged and adjusted by the following method:

6. The method for shallow sea multi-target direction finding using single vector hydrophone based on particle swarm optimization compensation as claimed in claim 5, characterized in that: In step 8, the instantaneous phase of the sound pressure channel is extracted as follows: Where n'=1,2,···,N', N' is the data length, p is the discrete time domain sound pressure signal, P(f m ) is the frequency domain signal of the sound pressure channel, m = 1, 2, ···, M, where M is the number of all frequency units, To extract the instantaneous phase, the spectrum expression obtained after weighting the frequency domain signal is: Among them, Vx(f m )、Vy(f m ) and Vz(f m ) is the channel frequency domain signal, is the weighted value obtained after the variance average fitting and mean filtering of the instantaneous phase difference values ​​of the four channels, and is the spectrum weighted by the instantaneous phase difference.

7. The method for shallow sea multi-target direction finding using single vector hydrophone based on particle swarm optimization compensation as claimed in claim 6, characterized in that: In step 10, the left and right boundary widths are less than 10; the target azimuth beam slope exceeds 1; the amplitude exceeds 3 times the average acoustic energy flow size, and the acoustic energy flow beam domain information is automatically extracted; the instantaneous frequency corresponding to the target azimuth is obtained as f g , g=1,2,...,G, the corresponding frequency unit is g, and the corresponding instantaneous frequency under interference noise is f y , y=1,2,...,Y, the corresponding frequency unit is y; At horizontal speed For example, the initial value and subsequent value of each frequency window are set as follows: Where W is the window weight, n = n'T s , T s is the sampling period; To change the window weight, use the following formula: Among them, W t is the size of the window weight at time t, a1 and a2 are the weight changes of the target instantaneous frequency and the noise instantaneous frequency in each iteration, and ε1 is the target instantaneous frequency window width.

8. The method for shallow sea multi-target direction finding using single vector hydrophone based on particle swarm optimization compensation as claimed in claim 7, characterized in that: In step 12, the phase compensation method for the horizontal dimension and the vertical dimension is as follows: in, and is the phase compensation angle in the horizontal and vertical dimensions, v' x (n), v' y (n) and v' z (n) is the out-of-phase wave sequence after compensation.

9. The method for shallow sea multi-target direction finding using single vector hydrophone based on particle swarm optimization compensation as claimed in claim 8, characterized in that: In step 13, in the high-precision azimuth estimation of the multi-target horizontal and vertical dimensions, the direction-finding sound pressure velocity conjugate cross-spectrum method is as follows: in, and V x (f) V y (f) and V z (f) The conjugate spectrum; After averaging the cross-spectrum results, the real part is taken to obtain the acoustic energy flow I of each frequency energy. x (f) I y (f) and I z (f) Use the acoustic energy flow to estimate the horizontal azimuth and elevation angle at the frequency point: What you get and are the estimated values ​​of horizontal azimuth and pitch angle respectively; Divide [0°, 360°] into several intervals Δθ, and divide [-90°, 90°] into several intervals Δφ, then the total number of intervals allowed for the spatial angle is: Scan the entire frequency band to obtain the corresponding angle. Assume that the number of estimated azimuths of each frequency point in the horizontal dimension that fall into each statistical angle interval is The number of vertical dimensions is Then calculate the following statistic R(θ): Obtain the distribution of acoustic energy flow in each statistical interval; In the above formula, R θ (f) and R φ (f) is the modulus of the acoustic energy flow that conforms to the angle range of the statistical interval, and the specific calculation formula is as follows: The modulus value is used as the weight of the corresponding azimuth count, and the estimated values ​​of each frequency point are counted. The final result is the horizontal azimuth and pitch angle.

Images