CN110196407B - Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation - Google Patents

Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation Download PDF

Info

Publication number
CN110196407B
CN110196407B CN201910370347.8A CN201910370347A CN110196407B CN 110196407 B CN110196407 B CN 110196407B CN 201910370347 A CN201910370347 A CN 201910370347A CN 110196407 B CN110196407 B CN 110196407B
Authority
CN
China
Prior art keywords
signal
frequency
incoming wave
wave direction
vector hydrophone
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
CN201910370347.8A
Other languages
Chinese (zh)
Other versions
CN110196407A (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.)
Southeast University
Original Assignee
Southeast University
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 Southeast University filed Critical Southeast University
Priority to CN201910370347.8A priority Critical patent/CN110196407B/en
Publication of CN110196407A publication Critical patent/CN110196407A/en
Application granted granted Critical
Publication of CN110196407B publication Critical patent/CN110196407B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G01S3/802—Systems for determining direction or deviation from predetermined direction
    • G01S3/803—Systems for determining direction or deviation from predetermined direction using amplitude comparison of signals derived from receiving transducers or transducer systems having differently-oriented directivity characteristics

Abstract

The invention discloses a single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation, which comprises the following steps: acquiring sound pressure of a single-vector hydrophone signal to be processed and a vibration velocity signal acquisition sequence in the X, Y direction, and calculating discrete Fourier transform; carrying out initial estimation on the incoming wave direction of the signal to obtain the incoming wave direction of the signal; selecting different frequency deviations to estimate the signal frequency according to the initially estimated signal incoming wave direction; respectively calculating the single-point Fourier transform of the sound pressure signal, the vibration velocity signal acquisition sequences in the X direction and the Y direction at the estimated signal frequency; and respectively carrying out conjugate multiplication on the single-point Fourier transform result of the sound pressure signal acquisition sequence and the single-point Fourier transform of the X, Y direction vibration velocity signal acquisition sequence, calculating to obtain a real part of the conjugate multiplication, and substituting into an inverse trigonometric function to calculate to obtain the incoming wave direction of the single-vector hydrophone signal. The method can more accurately extract the frequency spectrums of the three signals of the vector hydrophone, and has higher estimation accuracy under the condition of the same signal-to-noise ratio.

Description

Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation
Technical Field
The invention relates to a single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation, and belongs to the technical field of signal processing.
Background
Under the non-cooperative condition, the accurate estimation of the incoming wave direction of a single-vector hydrophone signal polluted by noise is one of research hotspots in signal processing, and the technology has wide application in the fields of speech signal processing, radar, sonar, electronic warfare and the like.
The key of the incoming wave direction estimation of the single vector hydrophone signal is to accurately extract angle information carried by the vector hydrophone receiving signal. The current commonly used single-vector hydrophone direction finding methods mainly comprise two methods: the mean-power method and the cross-spectrum power method.
The average sound energy flow method is to multiply the X, Y directional vibration velocity signal received in a period of time with the sound pressure signal in time domain, and then directly substitute the division result into the inverse trigonometric function for calculation. The method is simple to implement and clear in principle, but the estimation precision of the method is sharply reduced under the condition of low signal-to-noise ratio, so that the direction estimation error in a practical environment is large, and the engineering practicability is poor.
The cross-spectrum acoustic energy flow method is characterized in that a sound pressure signal and vibration velocity signals in the X and Y directions are subjected to DFT conversion, then cross-spectrum calculation is respectively carried out on the DFT conversion of the vibration velocity signals in the X and Y directions and the sound pressure signals, and a real part of a calculation result is taken and then substituted into an inverse trigonometric function for calculation. The processing method has higher pointing accuracy than the average acoustic energy flow method and is widely applied to practical engineering. However, when the conventional cross-spectrum method is used for obtaining a DFT conversion result, there may be a problem that the output signal-to-noise ratio is reduced due to spectrum leakage.
Disclosure of Invention
The invention aims to solve the technical problem that the conventional cross-spectrum acoustic energy flow method has the defect of low output signal-to-noise ratio caused by frequency spectrum leakage, and provides a single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation.
The invention specifically adopts the following technical scheme to solve the technical problems:
a single vector hydrophone signal incoming wave direction estimation method based on frequency estimation comprises the following steps:
the method comprises the following steps: sound pressure signal acquisition sequence p (n) for acquiring single vector hydrophone signal to be processed1) And a vibration velocity signal acquisition sequence v in the direction X, Yx(n2)、vy(n3) Wherein p (n)1)，n1＝0，1，...N-1，vx(n2)，n2＝0，1，...N-1，vy(n3)，n3N-1, said N1、n2、n3Respectively representing sampling points, wherein N represents the number of the sampling points, is an integer power of 2, and is more than or equal to 4;
step two: respectively calculating sound pressure signal acquisition sequences p (n)1) And a vibration velocity signal acquisition sequence v in the direction X, Yx(n2)、vy(n3) Discrete Fourier transforms X (k), V ofx(k) And Vy(k) Wherein k represents X (k), Vx(k) And Vy(k) The discrete frequency index of (a);
step three: according to calculated discrete Fourier transforms X (k), Vx(k) And Vy(k) The incoming wave direction of the signal is initially estimated to obtain an initially estimated incoming wave direction theta of the signal1；
Step four: according to the initially estimated signal incoming wave direction theta1Estimating signal frequency
Figure BDA0002049725920000021
The method comprises the following specific steps:
step (4.1) by discrete Fourier transform X (k), Vx(k) And Vy(k) Respectively calculating to obtain frequency deviation deltap、δy、ξx；
Step (4.2) according to the initially estimated signal incoming wave direction theta1Using different frequency deviations deltap、δy、δxCalculating the relative frequency deviation δ:
if sin θ1|＞|cosθ1Selecting frequency deviation of sound pressure and vibration speed in Y direction to calculate relative frequency deviation delta:
Figure BDA0002049725920000022
if | cos θ1|＞|sinθ1Selecting frequency deviation of sound pressure and vibration speed in the X direction to calculate relative frequency deviation delta:
Figure BDA0002049725920000023
if | cos θ1|＝|sinθ1Selecting sound pressure, and calculating a relative frequency deviation delta from the frequency deviation of the vibration speed in the X direction and the vibration speed in the Y direction:
Figure BDA0002049725920000024
step (4.3) estimating the signal frequency by using the relative frequency deviation delta
Figure BDA0002049725920000025
Figure BDA0002049725920000026
Where Δ f is the frequency resolution of a discrete fourier transform of length N, where Δ f ═ fs/N，fsThe signal is a sampling frequency;
step five: respectively calculating sound pressure signals, X-direction vibration velocity signal acquisition sequences p (n) and Y-direction vibration velocity signal acquisition sequences p (n)1)、vx(n2) And vy(n3) At the estimated signal frequency
Figure BDA0002049725920000035
Result of single point fourier transform of (Z)p，Zx，Zy；
Step six: the single-point Fourier transform result Z of the sound pressure signal acquisition sequencepRespectively with the single-point Fourier transform result Z of the X, Y direction vibration velocity signal acquisition sequencex，ZyCarrying out conjugate multiplication and calculating to obtain a real part xi of the conjugate multiplicationxAnd xiy；
Step seven: according to the real part xi of conjugate multiplicationxAnd xiySubstituting the inverse trigonometric function to calculate the incoming wave direction of the single-vector hydrophone signal
Figure BDA0002049725920000036
Further, as a preferred technical solution of the present invention, in the first step, a sound pressure signal acquisition sequence p (n) of a single vector hydrophone signal to be processed is obtained1) And a vibration velocity signal acquisition sequence v in the direction X, Yx(n2)、vy(n3) Comprising receiving real-time acquisition data of N sampling points from a single-vector hydrophone as a signal acquisition sequence p (N) to be processed1)、vx(n2)、vy(n3) Or extracting data of N sampling points from the moment of detecting the signal from a memory as a signal acquisition sequence p (N) to be processed1)、vx(n2)、vy(n3)。
Further, as a preferred technical solution of the present invention, in the second step, discrete fourier transforms x (k), V are calculatedx(k) And Vy(k) The formula is adopted:
Figure BDA0002049725920000031
Figure BDA0002049725920000032
Figure BDA0002049725920000033
wherein k represents X (k), Vx(k) And Vy(k) J denotes an imaginary unit, i.e.
Figure BDA0002049725920000034
Further, as a preferred technical solution of the present invention, in the third step, the initial estimation is performed to obtain an initial estimated signal incoming wave direction θ1The method specifically comprises the following steps:
step (3.1) calculating discrete Fourier transform X (k) and V respectivelyx(k) And Vy(k) Square of cross-power spectral modulus of (2) Px(k)，Py(k)：
Px(k)＝(|X(k)||Vx(k)|)2，k＝0，1，…N-1
Py(k)＝(|X(k)||Vy(k)|)2，k＝0，1，...N-1
Wherein | | | represents a modulo value operation;
step (3.2) of subjecting said P tox(k) And Py(k) Adding to obtain a signal power function P (k):
P(k)＝Px(k)+Py(k)，k＝0，1，…N-1
step (ii) of(3.3) searching a discrete frequency index k corresponding to the maximum value of the signal power function P (k)0：
Figure BDA0002049725920000041
Wherein arg max1≤k≤N/2-1[P(k)]Representing that the discrete frequency index corresponding to the maximum value of | P (k) | is searched in the range of 1 ≦ k ≦ N/2-1;
step (3.4) taking discrete Fourier transform X (k), Vx(k) And Vy(k) At k0Value X (k) of0)、Vx(k0)、Vy(k0) Respectively conjugate multiplied to obtain the real part xix1And xiy1：
ξx1＝Re[X(k0)]Re[Vx(k0)]+Im[X(k0)]Im[Vx(k0)]
ξy1＝Re[X(k0)]Re[Vy(k0)]+Im[X(k0)]Im[Vy(k0)]
Step (3.5) is to calculate the real part xix1And xiy1Substituting the arctangent function into the signal to perform initial estimation to obtain the initial estimated signal incoming wave direction theta1：
Figure BDA0002049725920000042
Wherein, tan-1() Is an arctan operation.
Further, as a preferred technical solution of the present invention, in the step (4.1), discrete fourier transform x (k), V is performedx(k) And Vy(k) Respectively calculating to obtain frequency deviation deltap、δy、δxThe formula is adopted:
Figure BDA0002049725920000043
Figure BDA0002049725920000044
Figure BDA0002049725920000045
wherein, X (k)0)、Vx(k0)、Vy(k0) Respectively, discrete Fourier transforms X (k), Vx(k) And Vy(k) At k is0A value of (b), and X (k)0-1)、Vx(k0-1)、Vy(k0-1) the discrete Fourier transforms X (k), V, respectivelyx(k) And Vy(k) At k0-a value at 1, and X (k)0+1)、Vx(k0+1)、Vy(k0+1 represents the discrete Fourier transforms X (k), V, respectivelyx(k) And Vy(k) At k0A value at + 1.
Further, as a preferred technical solution of the present invention, in the fifth step, a sound pressure signal, an X-direction and Y-direction vibration velocity signal acquisition sequence p (n) is calculated1)、vx(n2) And vy(n3) In estimating the signal frequency
Figure BDA0002049725920000051
Result of single point fourier transform of (Z)p，Zx，ZyAnd the formula is adopted:
Figure BDA0002049725920000052
Figure BDA0002049725920000053
Figure BDA0002049725920000054
further, the air conditioner is provided with a fan,as a preferred embodiment of the present invention, in the sixth step, a real part ξ of conjugate multiplication is calculatedxAnd xiyThe formula is adopted:
ξx＝Re[Zp]Re[Zx]+Im[Zp]Im[Zx]
ξy＝Re[Zp]Re[Zy]+Im[Zp]Im[Zy]
wherein Re [ ] represents a real part; im [ ] denotes the imaginary part.
Further, as a preferred technical solution of the present invention, in the seventh step, the incoming wave direction of the single-vector hydrophone signal is calculated
Figure BDA0002049725920000055
The formula is adopted:
Figure BDA0002049725920000056
by adopting the technical scheme, the invention can produce the following technical effects:
the method provided by the invention is an improvement on the basis of a conventional cross-spectrum method, and comprises the steps of firstly carrying out initial estimation on the incoming wave direction of a signal by utilizing sound pressure, X-direction and Y-direction signal sampling data sequences received by a single-vector hydrophone, and then accurately estimating the signal frequency according to the initially estimated incoming wave direction of the signal, namely selecting different sampling sequence combinations according to the initially estimated angle to correct the signal frequency to obtain the finally estimated signal frequency; and then, performing single-point Fourier transform on the estimated frequency point, solving a real part of the conjugate multiplication of the X-direction vibration velocity signal and the sound pressure signal single-point Fourier transform, and substituting the real part of the conjugate multiplication of the Y-direction vibration velocity signal and the sound pressure signal single-point Fourier transform into an inverse trigonometric function for calculation to obtain the incoming wave direction of the signal.
Therefore, compared with the prior art, the invention has the following advantages:
1. the estimation method of the invention utilizes the discrete Fourier transform X (k) of the sound pressure signal sampling data sequence and the vibration velocity signals in the X direction and the Y directionDiscrete Fourier transform V of a sequence of sampled datax(k) And Vy(k) The cross-power spectrum module value sum of squares is used for searching the discrete frequency index corresponding to the maximum value, and then the signal frequency is estimated.
2. The estimation method can more accurately extract the frequency spectrums of the three signals of the vector hydrophone through frequency estimation, avoids the problem of frequency spectrum leakage in the conventional method, and has higher estimation accuracy under the condition of the same signal-to-noise ratio; the lower limit of the signal-to-noise ratio which can be achieved by the method is lower under the condition of different signal-to-noise ratios. The method can effectively overcome the defect of low output signal-to-noise ratio caused by frequency spectrum leakage in the conventional cross-spectrum acoustic energy flow method, is simple to realize and has better engineering practicability.
Drawings
Fig. 1 is a flow chart of a method for estimating the incoming wave direction of a single-vector hydrophone signal based on frequency estimation.
Fig. 2 shows a normalized cross-power spectral module sum of squares of signals of the simulated single-vector hydrophone in embodiment 1 of the present invention.
Fig. 3 shows a normalized cross-power spectral module sum of squares of signals of the simulated single-vector hydrophone in embodiment 2 of the present invention.
Detailed Description
The following describes embodiments of the present invention with reference to the drawings.
As shown in fig. 1, the present invention provides a method for estimating an incoming wave direction of a single-vector hydrophone signal based on frequency estimation, which comprises the following steps:
the method comprises the following steps: acquiring a single-vector hydrophone signal sound pressure signal acquisition sequence to be processed and a two-dimensional vibration velocity signal acquisition sequence p (n)1)，n1＝0，1，...N-1，vx(n2)，n2＝0，1，...N-1，vy(n3)，n3N-1, wherein p (N)1) Sound pressure signal acquisition sequence, v, for a single vector hydrophonex(n2) For the acquisition of sequences of the vibration velocity signals of a single-vector hydrophone in the X direction, vy(n3) Acquiring a sequence of vibration velocity signals of the single-vector hydrophone in the Y direction; preferably, real-time acquisition data of N sampling points can be received from the single-vector hydrophone as a data sequence p (N) to be processed1)，n1＝0，1，...N-1，vx(n2)，n2＝0，1，...N-1，vy(n3)，n3N-1 or extracting data of N sampling points from the time of detecting a signal from a memory as a data sequence p (N) to be processed1)，n1＝0，1，...N-1，vx(n2)，n2＝0，1，...N-1，vy(n3)，n3N-1, said N1、n2、n3Respectively representing sampling points; n is the number of sampling points corresponding to the detected pulse width length of the signal, the value is an integral power of 2, and N is more than or equal to 4.
Step two: respectively calculating sound pressure, X-direction vibration velocity and Y-direction vibration velocity signal acquisition data sequence p (n)1)、vx(n2) And vy(n3) Discrete Fourier transform X (k), V of (2)x(k) And Vy(k) The calculation process is as follows:
Figure BDA0002049725920000061
Figure BDA0002049725920000062
Figure BDA0002049725920000071
wherein k represents X (k), Vx(k) And Vy(k) J denotes an imaginary unit, i.e.
Figure BDA0002049725920000072
In step two, a data sequence p (n) is acquired1)，vx(n2) And vy (n)3) The discrete fourier transform of (1), (2) and (3) is realized by fast fourier transform, and the computation amount of the algorithm can be reduced by using the fast fourier transform, thereby improving the computation efficiency of the algorithm.
Step three: according to the discrete Fourier transform X (k), Vx(k) And Vy(k) The incoming wave direction of the signal is initially estimated to obtain an initially estimated incoming wave direction theta of the signal1The estimation process is as follows:
step (3.1) calculating discrete Fourier transform X (k) and V respectivelyx(k) And Vy(k) Square P of cross-power spectral modulus ofx(k)，Py(k)：
Px(k)＝(|X(k)||Vx(k)|)2K is 0, 1, … N-1 formula (4)
Py(k)＝(|X(k)||Vy(k)|)2N-1 formula (5) with k being 0, 1
Wherein | | | represents a modulo value operation;
step (3.2) of adding Px(k) And Py(k) Adding to obtain a signal power function P (k):
P(k)＝Px(k)+Py(k) k is 0, 1, … N-1 formula (6)
Step (3.3) of searching discrete frequency index k corresponding to maximum value of signal power function P (k)0：
k0＝arg max1≤k≤N/2-1[P(k)]Formula (7)
Wherein argmax1≤k≤N/2-1[P(k)]Representing that the discrete frequency index corresponding to the maximum value of | P (k) | is searched in the range of 1 ≦ k ≦ N/2-1;
step (3.4) taking the discrete Fourier transform X (k), V of the three-way signalx(k) And Vy(k) At k is0Value X (k) of0)、Vx(k0)、Vy(k0) Respectively conjugate multiplying to obtain real part xix1And xiy1：
ξx1＝Re[X(k0)]Re[Vx(k0)]+Im[X(k0)]Im[Vx(k0)]Formula (8)
ξy1＝Re[X(k0)]Re[Vy(k0)]+Im[X(k0)]Im[Vy(k0)]Formula (9)
Step (3.5) is to calculate the real part xix1And xiy1Substituting the arctangent function, and carrying out initial estimation to obtain an initial estimated signal incoming wave direction:
Figure BDA0002049725920000081
wherein, tan-1() Is an arctan operation.
In the third step, the method is realized in five steps: firstly, calculating the discrete Fourier transform X (k) of the sound pressure signal acquisition sequence and the discrete Fourier transform V of the vibration velocity signal acquisition sequence in the X direction and the Y directionx(k) And Vy(k) Cross power spectrum module value square Px(k) And Py(k) (ii) a The second step: calculating Px(k) And Py(k) To obtain a power function p (k); thirdly, searching a discrete frequency index k corresponding to the maximum value of P (k)0(ii) a The fourth step: taking discrete Fourier transform of three-way signal at k0The values of (A) are respectively subjected to conjugate multiplication to obtain a real part xix1And xiy1(ii) a The fifth step, according to xix1And xiy1And obtaining the initial estimation of the incoming wave direction of the signal.
Step four: according to the initially estimated incoming wave direction theta of the signal1Estimating signal frequency
Figure BDA0002049725920000082
The calculation steps are as follows:
step (4.1) discrete Fourier transform X (k), V of three signalsx(k) And Vy(k) Separately calculating the frequency deviation deltap、δy、δxNamely:
Figure BDA0002049725920000083
Figure BDA0002049725920000084
Figure BDA0002049725920000085
wherein, X (k)0)、Vx(k0)、Vy(k0) Respectively, discrete Fourier transforms X (k), Vx(k) And Vy(k) At k0A value of (b), and X (k)0-1)、Vx(k0-1)、Vy(k0-1) the discrete Fourier transforms X (k), V, respectivelyx(k) And Vy(k) At k is0-a value at 1, and X (k)0+1)、Vx(k0+1)、Vy(k0+1 represents the discrete Fourier transforms X (k), V, respectivelyx(k) And Vy(k) At k0A value at + 1.
Step (4.2) according to the initially estimated signal incoming wave direction theta1Using different frequency deviations deltap、δy、δxCalculating the relative frequency deviation δ:
if sin θ1|＞|cosθ1Selecting frequency deviation of sound pressure and vibration speed in Y direction to calculate relative frequency deviation delta:
Figure BDA0002049725920000086
if | cos θ1|＞|sinθ1Selecting frequency deviation of sound pressure and vibration speed in the X direction to calculate relative frequency deviation delta:
Figure BDA0002049725920000087
if | cos θ1|＝|sinθ1Selecting sound pressure, and calculating the relative frequency deviation delta according to the frequency deviation of the vibration speed in the X direction and the vibration speed in the Y direction:
Figure BDA0002049725920000091
step (4.3) estimating the signal frequency by using the relative frequency deviation delta
Figure BDA0002049725920000092
Figure BDA0002049725920000093
Where Δ f is the frequency resolution of a discrete Fourier transform of length N, where Δ f ═ fs/N，fsThe signal is a sampling frequency;
in the fourth step, calculating the relative frequency deviation delta, and taking the relative frequency deviation delta as an estimated value of the relative frequency deviation of the signal of the single-vector hydrophone; the method is realized by three steps: firstly, respectively calculating frequency deviation through discrete Fourier transform of three signals; secondly, according to the estimated signal incoming wave direction theta1Selecting different signals to calculate relative frequency deviation delta; thirdly, estimating the signal frequency by using the relative frequency deviation delta
Figure BDA0002049725920000094
Step five: respectively calculating sound pressure, X direction and Y direction vibration velocity signal acquisition sequences p (n)1)、vx(n2) And vy(n3) At the estimated frequency point
Figure BDA0002049725920000095
Result of single point fourier transform of (Z)p，Zx，ZyThe calculation process is as follows:
Figure BDA0002049725920000096
Figure BDA0002049725920000097
Figure BDA0002049725920000098
step six: the single-point Fourier transform result Z of the sound pressure signal acquisition sequencepThe result Z of single-point Fourier transform of the vibration velocity signal acquisition sequence in the X direction and the Y direction respectivelyx，ZyCarrying out conjugate multiplication and calculating to obtain a real part xi of the conjugate multiplicationxAnd xiyThe calculation process is as follows:
ξx＝Re[Zp]Re[Zx]+Im[Zp]Im[Zx]formula (21)
ξy＝Re[Zp]Re[Zy]+Im[Zp]Im[Zy]Formula (22)
Step seven: according to the real part xi of conjugate multiplicationxAnd xiySubstituting the obtained data into an inverse trigonometric function to calculate and estimate the incoming wave direction of the single-vector hydrophone signal
Figure BDA0002049725920000099
The process is as follows:
Figure BDA00020497259200000910
wherein, tan-1() Is an arctan operation.
The method provided by the invention utilizes the sound pressure, X-direction and Y-direction signal acquisition sequences received by the single-vector hydrophone to carry out primary estimation on the incoming wave direction of the signal; and selecting different sampling sequence combinations according to the initial estimation angle to correct the signal frequency to obtain the final estimated signal frequency, performing single-point Fourier transform of sound pressure, X-direction and Y-direction signal acquisition sequences at the frequency point, then solving a conjugate multiplication real part of the X-direction vibration velocity signal and the sound pressure signal and a conjugate multiplication real part of the Y-direction vibration velocity signal and the sound pressure signal through the single-point Fourier transform, substituting the conjugate multiplication real parts into an inverse trigonometric function to calculate, and solving the final estimated angle of the incoming wave direction of the signal.
In the invention, a signal receiving model of the single-vector hydrophone is as follows:
Figure BDA0002049725920000101
Figure BDA0002049725920000102
where a is the amplitude of the signal received by the single vector hydrophone,
Figure BDA0002049725920000103
initial phase of signal received for single vector hydrophone, N is number of signal sampling points, f0Is the signal frequency, fsFor the sampling frequency, θ is the incoming wave direction of the signal, i.e. the value to be estimated, np(n1)，nx(n2)，ny(n3) The three paths of signals respectively receive mutually independent white Gaussian noise, the mean value of the three paths of signals is 0, and the variance is sigma2And the magnitude of the variance is determined by the signal-to-noise ratio SNR:
SNR＝10log10[A2/(2σ2)]
in order to verify that the method of the present invention can more accurately extract the frequency spectrum of the three-way signal of the vector hydrophone, and avoid the frequency spectrum leakage existing in the conventional method, two embodiments are listed for verification and explanation.
Examples 1,
In this embodiment, the simulation signal parameters are respectively set as: signal amplitude a 1, initial phase
Figure BDA0002049725920000105
Number of sampling points N is 1024, sampling frequency fs4000Hz, the frequency resolution of the discrete fourier transform, Δ f, is fs3.9063Hz, signal frequency f0277Hz, 50 ° of the incoming signal direction θ, and-3 dB of the SNR.
Fig. 2 is a normalized cross-power spectral module sum of squares of the simulated single-vector hydrophone signal with the frequency of 277Hz shown in this embodiment, and as can be seen from fig. 2, the index corresponding to the maximum value of the cross-power spectral module sum of squares is 71.
Searching an index corresponding to the square sum maximum of cross-power spectrum modulus to obtain k 071. The initial estimated angle θ obtained by using the equations (8) to (10)1＝49.6824°。
Due to | sin θ1|＞|cosθ1Estimating the frequency deviation delta by using the frequency deviation of the sound pressure and the vibration speed in the Y direction, namely the formula (14):
Figure BDA0002049725920000104
further calculating an estimated frequency according to equation (17)
Figure BDA0002049725920000111
Figure BDA0002049725920000112
Then, the following equation (18), equation (19), and equation (20) are used to obtain:
Zp＝-0.0150-0.5035i
Zx＝-0.0195-0.3100i
Zy＝-0.0120-0.3685i
further, the following equation (21) and equation (22) are used to obtain:
ξx＝0.1564，ξy＝0.1853
then, the following is obtained in place of formula (23):
Figure BDA0002049725920000113
estimating relative error of incoming wave direction
Figure BDA0002049725920000114
Example 2
In this embodiment, the simulation signal parameters are respectively set as: signal amplitude a 1, initial phase
Figure BDA0002049725920000115
Number of sampling points N is 1024, sampling frequency fs4000Hz, the frequency resolution of the discrete fourier transform, Δ f, is fs3.9063Hz, signal frequency f0314Hz, the incoming wave direction theta of the signal is 35 DEG, and the SNR is 3 dB.
Fig. 3 is a normalized cross-power spectral module sum of squares of the simulated single vector hydrophone signal with the frequency of 277Hz shown in this embodiment, and as can be seen from fig. 3, the index corresponding to the maximum value of the cross-power spectral module sum of squares is 80.
Searching an index corresponding to the square sum of the cross-power spectrum modulus and the maximum value to obtain k 080. The initial estimated angle θ obtained by using the equations (8) to (10)1＝34.2210°。
Due to | cos θ1|＞|sinθ1Estimating the frequency deviation delta by using the frequency deviation of the sound pressure and the vibration speed in the X direction, namely the formula (15):
Figure BDA0002049725920000116
further calculating an estimated frequency according to equation (17)
Figure BDA0002049725920000117
Figure BDA0002049725920000118
Then, the following equation (18), equation (19), and equation (20) are used to obtain:
Zp＝-0.0108-0.5114i
Zx＝-0.0028-0.4097i
Zy＝-0.0092-0.2876i
further, the following equation (21) and equation (22) are used to obtain:
ξx＝0.2095，ξy＝0.1472
then, the following is obtained in place of formula (23):
Figure BDA0002049725920000121
estimating relative error of incoming wave direction
Figure BDA0002049725920000122
It can be seen from the above embodiments 1 and 2 that the estimation method of the present invention can obtain good estimation accuracy, and is simple in calculation, small in calculation amount, and suitable for the occasion of fast estimating the incoming wave direction of the signal of the single-vector hydrophone with high accuracy.
Therefore, the method can more accurately extract the frequency spectrums of the three signals of the vector hydrophone through frequency estimation, avoids the problem of frequency spectrum leakage in the conventional method, and has higher estimation accuracy under the condition of the same signal-to-noise ratio; the lower limit of the signal-to-noise ratio which can be achieved by the method is lower under the condition of different signal-to-noise ratios.
The embodiments of the present invention have been described in detail with reference to the drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge of those skilled in the art without departing from the gist of the present invention.

Claims (8)

1. A single vector hydrophone signal incoming wave direction estimation method based on frequency estimation is characterized by comprising the following steps:
the method comprises the following steps: sound pressure signal acquisition sequence p (n) for acquiring single vector hydrophone signal to be processed1) And a sequence v of velocity signal acquisitions in the direction X, Yx(n2)、vy(n3) Wherein p (n)1)，n1＝0，1，...N-1，vx(n2)，n2＝0，1，...N-1，vy(n3)，n3N-1, said N1、n2、n3Respectively representing sampling points, wherein N represents the number of the sampling points, the value of N is an integer power of 2, and N is more than or equal to 4;
step two: respectively calculating sound pressure signal acquisition sequences p (n)1) And a vibration velocity signal acquisition sequence v in the direction X, Yx(n2)、vy(n3) Discrete Fourier transforms X (k), V ofx(k) And Vy(k) Wherein k represents X (k), Vx(k) And Vy(k) The discrete frequency index of (a);
step three: according to calculated discrete Fourier transforms X (k), Vx(k) And Vy(k) The incoming wave direction of the signal is initially estimated to obtain an initially estimated incoming wave direction theta of the signal1；
Step four: according to the initially estimated signal incoming wave direction theta1Estimating signal frequency
Figure FDA0003554233320000011
The method specifically comprises the following steps:
step (4.1) by discrete Fourier transform X (k), Vx(k) And Vy(k) Respectively calculating to obtain frequency deviation deltap、δy、δx；
Step (4.2) according to the initially estimated signal incoming wave direction theta1Using different frequency deviations deltap、δy、δxCalculating the relative frequency deviation δ:
if sin θ1|＞|cosθ1Selecting frequency deviation of sound pressure and vibration speed in Y direction to calculate relative frequency deviation delta:
Figure FDA0003554233320000012
if | cos θ1|＞|sinθ1Selecting sound pressure and X-direction vibrationFrequency deviation of velocity calculation relative frequency deviation δ:
Figure FDA0003554233320000013
if | cos θ1|＝|sinθ1Selecting sound pressure, and calculating the frequency deviation delta of the frequency deviation of the vibration speed in the X direction and the vibration speed in the Y direction:
Figure FDA0003554233320000014
step (4.3) estimating the signal frequency by using the relative frequency deviation delta
Figure FDA0003554233320000015
Figure FDA0003554233320000016
Where Δ f is the frequency resolution of a discrete fourier transform of length N, where Δ f ═ fs/N，fsThe signal is a sampling frequency; k is a radical of0A discrete frequency index corresponding to the maximum value of the signal power function P (k);
step five: respectively calculating sound pressure signals, X-direction vibration velocity signal acquisition sequences and Y-direction vibration velocity signal acquisition sequences p (n)1)、vx(n2) And vy(n3) In estimating the signal frequency
Figure FDA0003554233320000021
Result of single point fourier transform of (Z)p，Zx，Zy；
Step six: the single-point Fourier transform result Z of the sound pressure signal acquisition sequencepRespectively with the single-point Fourier transform result Z of the X, Y direction vibration velocity signal acquisition sequencex，ZyCarrying out conjugate multiplication and calculating to obtain a real part xi of the conjugate multiplicationxAnd xiy；
Step seven: according to the real part xi of conjugate multiplicationxAnd xiySubstituting the inverse trigonometric function to calculate the incoming wave direction of the single-vector hydrophone signal
Figure FDA0003554233320000022
2. The method for estimating the incoming wave direction of the single-vector hydrophone signal based on the frequency estimation as claimed in claim 1, wherein in the first step, a sound pressure signal acquisition sequence p (n) of the single-vector hydrophone signal to be processed is obtained1) And a sequence v of velocity signal acquisitions in the direction X, Yx(n2)、vy(n3) By receiving real-time acquisition data of N sampling points from a single-vector hydrophone as a signal acquisition sequence p (N) to be processed1)、vx(n2)、vy(n3) Or extracting data of N sampling points from the moment of detecting the signal from a memory as a signal acquisition sequence p (N) to be processed1)、vx(n2)、vy(n3)。
3. The method for estimating the incoming wave direction of the single-vector hydrophone signal based on the frequency estimation as claimed in claim 1, wherein in the second step, discrete Fourier transforms X (k), V (V) are calculatedx(k) And Vy(k) The formula is adopted:
Figure FDA0003554233320000023
Figure FDA0003554233320000024
Figure FDA0003554233320000025
wherein k represents X (k), Vx(k) And Vy(k) J denotes an imaginary unit, i.e.
Figure FDA0003554233320000026
4. The method for estimating the incoming wave direction of a single-vector hydrophone signal based on frequency estimation as claimed in claim 1, wherein the initial estimation performed in the third step is performed to obtain the initial estimated incoming wave direction θ1The method specifically comprises the following steps:
step (3.1) calculating discrete Fourier transform X (k) and V respectivelyx(k) And Vy(k) Square P of cross-power spectral modulus ofx(k)，Py(k)：
Px(k)＝(|X(k)||Vx(k)|)2，k＝0，1，...N-1
Py(k)＝(|X(k)||Vy(k)|)2，k＝0，1，...N-1
Wherein | | | represents a modulo value operation;
step (3.2) of subjecting said P tox(k) And Py(k) Adding to obtain a signal power function P (k):
P(k)＝Px(k)+Py(k)，k＝0，1，...N-1
step (3.3) searching discrete frequency index k corresponding to maximum value of signal power function P (k)0：
Figure FDA0003554233320000031
Wherein argmax1≤k≤N/2-1[P(k)]Representing that the discrete frequency index corresponding to the maximum value of | P (k) | is searched in the range of 1 ≦ k ≦ N/2-1;
step (3.4) taking discrete Fourier transform X (k), Vx(k) And Vy(k) At k is0Value X (k) of0)、Vx(k0)、Vy(k0) Respectively obtained after conjugate multiplicationReal part xix1And xiy1：
ξx1＝Re[X(k0)]Re[Vx(k0)]+Im[X(k0)]Im[Vx(k0)]
ξy1＝Re[X(k0)]Re[Vy(k0)]+Im[X(k0)]Im[Vy(k0)]
Step (3.5) is to calculate the real part xix1And xiy1Substituting the arctangent function into the signal to perform initial estimation to obtain the initial estimated signal incoming wave direction theta1：
Figure FDA0003554233320000032
Wherein, tan-1() Is an arctan operation.
5. The method for estimating the incoming wave direction of a single-vector hydrophone signal based on frequency estimation in claim 1, wherein the step (4.1) is performed by discrete Fourier transform (X (k), V)x(k) And Vy(k) Respectively calculating to obtain frequency deviation deltap、δy、δxThe formula is adopted:
Figure FDA0003554233320000033
Figure FDA0003554233320000034
Figure FDA0003554233320000035
wherein, X (k)0)、Vx(k0)、Vy(k0) Respectively, discrete Fourier transforms X (k), Vx(k) And Vy(k) At k0A value of (b), and X (k)0-1)、Vx(k0-1)、Vy(k0-1) the discrete Fourier transforms X (k), V, respectivelyx(k) And Vy(k) At k0-a value at 1, and X (k)0+1)、Vx(k0+1)、Vy(k0+1 represents the discrete Fourier transforms X (k), V, respectivelyx(k) And Vy(k) At k0A value at + 1.
6. The method for estimating the incoming wave direction of the single-vector hydrophone signal based on frequency estimation as claimed in claim 1, wherein the sound pressure signal, the vibration velocity signal acquisition sequences p (n) in the X direction and the vibration velocity signal acquisition sequences p (n) in the Y direction are calculated in the fifth step1)、vx(n2) And vy(n3) In estimating the signal frequency
Figure FDA0003554233320000041
Result of single point fourier transform of (Z)p，Zx，ZyThe formula is adopted:
Figure FDA0003554233320000042
Figure FDA0003554233320000043
Figure FDA0003554233320000044
7. the method for estimating the incoming wave direction of a single-vector hydrophone signal based on frequency estimation as recited in claim 1, wherein a real part xi of conjugate multiplication is calculated in the sixth stepxAnd xiyAnd the formula is adopted:
ξx＝Re[Zp]Re[Zx]+Im[Zp]Im[Zx]
ξy＝Re[Zp]Re[Zy]+Im[Zp]Im[Zy]
wherein Re [ ] represents a real part; im [ ] represents the imaginary part.
8. The method according to claim 1, wherein the incoming wave direction of the single-vector hydrophone signal is calculated in step seven
Figure FDA0003554233320000045
The formula is adopted:
Figure FDA0003554233320000046
CN201910370347.8A 2019-05-06 2019-05-06 Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation Active CN110196407B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910370347.8A CN110196407B (en) 2019-05-06 2019-05-06 Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910370347.8A CN110196407B (en) 2019-05-06 2019-05-06 Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation

Publications (2)

Publication Number Publication Date
CN110196407A CN110196407A (en) 2019-09-03
CN110196407B true CN110196407B (en) 2022-06-24

Family

ID=67752418

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910370347.8A Active CN110196407B (en) 2019-05-06 2019-05-06 Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation

Country Status (1)

Country Link
CN (1) CN110196407B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112415467B (en) * 2020-11-06 2022-10-25 中国海洋大学 Single-vector subsurface buoy target positioning implementation method based on neural network
CN112816937B (en) * 2020-12-23 2022-10-28 中国船舶重工集团有限公司第七一0研究所 Helicopter scale estimation method and device based on fundamental frequency spectrum azimuth angle change rate
CN112881970A (en) * 2021-01-18 2021-06-01 湖南国天电子科技有限公司 Target direction finding method and system suitable for sonar buoy platform

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102981146A (en) * 2012-11-19 2013-03-20 哈尔滨工程大学 Passive positioning method of single-vector hydrophone
CN105182345A (en) * 2015-09-26 2015-12-23 哈尔滨工程大学 Continuous spectrum signal target automatic detection method based on single vector subsurface buoy
CN106066468A (en) * 2016-05-25 2016-11-02 哈尔滨工程大学 A kind of based on acoustic pressure, the vector array port/starboard discrimination method of vibration velocity Mutual spectrum
CN106546949A (en) * 2016-10-28 2017-03-29 东南大学 A kind of double array element sinusoidal signal arrival bearing's methods of estimation based on frequency estimation meter
CN107728109A (en) * 2017-09-18 2018-02-23 哈尔滨工程大学 A kind of noncooperative target radiated noise measurement and positioning technology

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102981146A (en) * 2012-11-19 2013-03-20 哈尔滨工程大学 Passive positioning method of single-vector hydrophone
CN105182345A (en) * 2015-09-26 2015-12-23 哈尔滨工程大学 Continuous spectrum signal target automatic detection method based on single vector subsurface buoy
CN106066468A (en) * 2016-05-25 2016-11-02 哈尔滨工程大学 A kind of based on acoustic pressure, the vector array port/starboard discrimination method of vibration velocity Mutual spectrum
CN106546949A (en) * 2016-10-28 2017-03-29 东南大学 A kind of double array element sinusoidal signal arrival bearing's methods of estimation based on frequency estimation meter
CN107728109A (en) * 2017-09-18 2018-02-23 哈尔滨工程大学 A kind of noncooperative target radiated noise measurement and positioning technology

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于频率预估的单矢量水听器测向方法;方子彦 等;《声学技术》;20190430;第38卷(第2期);全文 *

Also Published As

Publication number Publication date
CN110196407A (en) 2019-09-03

Similar Documents

Publication Publication Date Title
CN110196407B (en) Single-vector hydrophone signal incoming wave direction estimation method based on frequency estimation
CN103941089B (en) Sinusoidal signal frequency method of estimation based on DFT
CN110007148B (en) Single-frequency signal frequency estimation method based on comprehensive interpolation of discrete spectrum phase and amplitude
CN101893698B (en) Noise source test and analysis method and device
CN105589056A (en) Multi-objective near-and-far field mixed source positioning method
CN106646350B (en) A kind of modification method when each channel amplitude gain of single vector hydrophone is inconsistent
CN101150345A (en) Direction measurement method applicable to phase interference signal source under non stabilized noise background
CN111474521A (en) Sound source positioning method based on microphone array in multipath environment
CN111610503B (en) Linear frequency modulation signal parameter estimation method based on improved LVD
CN104714235A (en) Ranging method and system for double low-frequency vector hydrophone arrays
CN109239680B (en) Parameter estimation method for low interception probability radar LFM signal
CN109061599B (en) STAP method based on cyclostationarity and symmetric prior knowledge
CN109541573A (en) A kind of element position calibration method being bent hydrophone array
CN111781590B (en) Efficient FRFT-based radar target parameter estimation method
CN108957389A (en) A kind of real number field multi channel signals method for estimating target azimuth
CN109541572B (en) Subspace orientation estimation method based on linear environment noise model
CN109117698B (en) Noise background estimation method based on minimum mean square error criterion
CN111722178B (en) Far-field narrow-band signal incoming wave direction estimation method based on numerical solution of directivity model
CN103185879A (en) Method for detecting single-channel synthetic aperture radar moving target with high radial velocity target
CN112162235B (en) Smooth segmented stochastic resonance enhanced acoustic vector signal orientation method
CN107238813B (en) Method and device for determining direction of arrival and time of arrival of near-field signal source
CN111007473A (en) High-speed weak target detection method based on distance frequency domain autocorrelation function
Wang et al. An improved phase correction algorithm in extended towed array method for passive synthetic aperture
CN100520442C (en) Arrival time difference positioning method by total least square equalization algorithms
Deng et al. A method of extracting underwater acoustic beaconing signal

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